Why are single-electron states occupied at $T=0$ in superconductors? Some textbooks show a density of states where the single-electron states below the Fermi energy are occupied at $T=0$ (see picture). However, I thought that at $T=0$ all electrons are paired. Hence this would mean that they are bosons and should not show up in the density of states for single electrons. How is that possible?

 A: In few words, I believe that this picture is misleading, since there is a totally equivalent picture which makes your point more intuitive. However, I need to go through the details.
Mathematical analysis
You can see this very clearly from the BCS theory. The BCS Hamiltonian reads
$$
H = \sum_k \Psi_k^{\dagger} \left( 
\begin{array}{cc}
\varepsilon_k & \Delta \\
\Delta & -\varepsilon_k
\end{array}
\right) \Psi_k,
$$
where I have introduced the Nambu spinor containing usual electron operators
$$
\Psi_k = \left( \begin{array}{c} c_{k\uparrow} \\ c^{\dagger}_{-k\downarrow} \end{array} \right).
$$
The BCS Hamiltonian does not commute with the particle number operator, meaning that the number of fermions in the system is not a good quantum number, and is not fixed.
Now performing a suitable unitary transformation, you can diagonalize the matrix and get a diagonal Hamiltonian in terms of new fermionic operators $\alpha_{k\sigma}$, $\sigma=\uparrow, \downarrow$:
$$
H = \sum_{k\sigma} \lambda_k \alpha^{\dagger}_{k\sigma} \alpha_{k\sigma} - \sum_k \lambda_k
$$
where $\lambda_k = \sqrt{\varepsilon_k^2 + \Delta^2}$.
Notice that $\lambda_k >0$ for every $k$, as long as the system is gapped ($\Delta\neq 0$).
So what is the ground state of such system, provided that you don't have any constraint on the particle number?
Well, it's the state with no $\alpha$-fermions at all, since the presence of any $\alpha$-fermion would increase the energy by an amount $\lambda_k$!
Moreover, the ground state energy is $E_0 = - \sum_k \lambda_k$.
This makes totally sense: all the electrons are in the condensate, $\alpha$-fermions represent quasi-particle excitations on top of it, and in the ground state there are no such excitations.
Now the misleading - but equivalent - picture.
There is a trivial trick to make things more symmetric: let's write the identity $\lambda_k \alpha^{\dagger}_{k\sigma} \alpha_{k\sigma} = (\lambda_k/2) (\alpha^{\dagger}_{k\sigma} \alpha_{k\sigma} + 1 - \alpha_{k\sigma} \alpha^{\dagger}_{k\sigma})$, so that the Hamiltonian becomes
$$
H = \sum_k \left( \frac{\lambda_k}{2} \alpha^{\dagger}_{k\sigma} \alpha_{k\sigma} - \frac{\lambda_k}{2} \alpha_{k\sigma} \alpha^{\dagger}_{k\sigma} \right),
$$
Finally let's rewrite the second term introducing the operator $\beta^{\dagger}_{k\sigma} = \alpha_{k\sigma}$ (notice that a $\beta$-fermion can be interpreted as a hole of $\alpha$-fermions):
$$
H = \sum_k \left( \frac{\lambda_k}{2} \alpha^{\dagger}_{k\sigma} \alpha_{k\sigma} - \frac{\lambda_k}{2} \beta^{\dagger}_{k\sigma} \beta_{k\sigma} \right).
$$
The plot you are showing clearly comes out from this Hamiltonian (notice the symmetryc energy bands $\pm \lambda_k/2$). What is the ground state now? Again the state with no $\alpha$-particles, corresponding to the state with the maximum amount of $\beta$-particles! The ground state energy is $E_0 = -\sum_k \lambda_k$, which is again consistent.
The latter approach is more nice and symmetric, but also misleading, because it suggests the wrong idea that the ground state is made by electrons filling up the lowest energy band. The correct interpretation is that the system is filled up with $\beta$-fermions, or equivalently it is completely empty of $\alpha$-fermions: in either case, if you write down the state in terms of the original electron operators $c_{k\sigma}$, you will get the BCS ground state, i.e. a coherent state of condensed electronic pairs :)
