Density of states in $k$ space for the BCS Hamiltonian I would like to understand how we can find the density of states in $k$ space for the BCS Hamiltonian.
First, let's talk about free electrons. When we deal with free electrons, the Hamiltonian is $\frac{p^2}{2m}$ and it is more simple to write it in $k$ space: $$H=\frac{h^2k^2}{2m}.$$
An eigenvector is then $\psi_k=\delta(k-k_0)$ and we know how many modes we have in an interval $dk$ because of Fourier series. If we have a function confined in a box of size $L$, we can periodise the function with period $L$ on all the space, and two Fourier modes will be separated by $\frac{2 \pi}{L}$.
So, we have $\rho(k)dk= \frac{Ldk}{2 \pi}$ modes in $dk$.
To summarize: the $k$ is really a Fourier mode here. And because of Fourier theory I know how many modes I have in $dk$.
Now, lets talk about the BCS Hamiltonian,
$$ H=\sum_{k,\sigma} E_k \hat{a}^{\dagger}_{k, \sigma}\hat{a}_{k, \sigma} -|g_{eff}|^2 \sum_{k_1, k_2, \sigma_1, \sigma_2}\hat{a}^{\dagger}_{k1+q, \sigma_1}\hat{a}^{\dagger}_{k2-q, \sigma_2}\hat{a}_{k2, \sigma_2}\hat{a}_{k1, \sigma_1}$$
and after diagonalisation, we have
$$ H= \sum_k E_k \gamma^{\dagger}_{k+}\gamma_{k+} -E_k \gamma^{\dagger}_{k-}\gamma_{k-}.$$
I note $$|\phi^+_k\rangle=\gamma^{\dagger}_{k+}|0\rangle_{BCS}.$$
$k$ is now just a quantum number used to "note" an eigenvector, it is not a Fourier mode. How can one determine the density of states in $k$ space now? 
 A: You almost get the point, although the "Fourier mode" is not the right terminology here. The problem with BCS quasiparticles is that they are not electron modes, but electron-hole superposition modes. This means the quasiparticle is partially electron and partially hole (their probabilities add up to one). On the other hand, the density of state (as measured in STM or ARPES experiments) only counts the density of electron component, therefore each quasiparticle mode in the spectrum must be weighted by the portion (the probability) of electron in that mode. So knowing the BCS quasiparticle dispersion is not sufficient to calculate the density of state, we must also know the amount of electronic component at each point. One can go on with this idea and carry out the calculation, but this undergraduate-level approach is not systematic and is sometimes tedious.
So let us introduce the quantum field theory approach to calculate the density of state systematically, which applies to both normal metals and BCS superconductors. The key observation is that the electron band (dispersion relation) is not well defined in BCS superconductors anymore, it must be replaced by the concept of electron spectral weight in the momentum-frequency space. The spectral weight (or the spectral function) $A(\omega,\boldsymbol{k})$ describes the probability density of observing an electron mode at frequency $\omega$ and momentum $\boldsymbol{k}$. Then the density of state (DOS) is just the integral of the spectral weight over the momentum space or over the Brillouin zone (BZ)
$$\text{DOS}(\omega)=\int_\text{BZ} A(\omega,\boldsymbol{k})\mathrm{d}^d\boldsymbol{k}.$$
The spectral function can be calculated as the imaginary part of the electron Green's function $A(\omega,\boldsymbol{k})=-2\text{Im}G(\omega+\text{i}0_+,\boldsymbol{k})$, where the Green's function is defined as $G(t-t',\boldsymbol{x}-\boldsymbol{x}')=-\langle c(t,\boldsymbol{x})c^\dagger(t',\boldsymbol{x}')\rangle$ in the space-time, and given by $G(\omega,\boldsymbol{k})=(\omega-h(\boldsymbol{k}))^{-1}$ after Fourier transforming to the momentum-frequency space.
Now we apply this approach to the BCS superconductor described by the following Hamiltonian
$$H_\text{BCS}=\sum_{\boldsymbol{k}}c_{\boldsymbol{k}\sigma}^\dagger \epsilon_{\boldsymbol{k}} c_{\boldsymbol{k}\sigma}+\Delta(c_{\boldsymbol{k}\uparrow}c_{-\boldsymbol{k}\downarrow}+\text{h.c.}),$$
with $\epsilon_{\boldsymbol{k}}=\boldsymbol{k}^2/(2m)-\mu$ and $\Delta\in\mathbb{R}$ for example. In the Nambu basis $\psi_{\boldsymbol{k}}=(c_{\boldsymbol{k}\uparrow},c_{-\boldsymbol{k}\downarrow}^\dagger)$, we have $$\begin{split}H_\text{BCS}&=\sum_{\boldsymbol{k}}\psi_{\boldsymbol{k}}^\dagger h(\boldsymbol{k})\psi_{\boldsymbol{k}},\\
h(\boldsymbol{k})&=\epsilon_{\boldsymbol{k}}\sigma^3+\Delta\sigma^1,\end{split}$$
where $\sigma^3$ and $\sigma^1$ denote the Pauli matrices acting on the Nambu spinor $\psi_{\boldsymbol{k}}$. The Green's function of the Nambu spinor is given by
$$G(\omega,\boldsymbol{k})=(\omega-h(\boldsymbol{k}))^{-1}=\frac{\omega\sigma^0+\epsilon_{\boldsymbol{k}}\sigma^3+\Delta\sigma^1}{\omega^2-\epsilon_{\boldsymbol{k}}^2-\Delta^2}.$$ 
Note that this $G(\omega,\boldsymbol{k})$ is a $2\times 2$ matrix at each momentum-frequency point. The (1,1)-component $G(\omega,\boldsymbol{k})_{11}$ corrsponds to the correlation $-\langle c_{\boldsymbol{k}\uparrow}c_{\boldsymbol{k}\uparrow}^\dagger\rangle$ and is therefore the Green's function of the up-spin electron. We expect the Green's function of the down-spin electron to be the same as the up-spin electron due to the time-reversal symmetry. With this, we can compute the electron spectral function
$$A(\omega,\boldsymbol{k})_{11}=-2\text{Im}G(\omega+\text{i}0_+,\boldsymbol{k})_{11}.$$
The result is shown in the following figure. Darker color means larger spectral weight and higher probability of observing the electron at this momentum-frequency point.

One can see the electron dispersion (opening upwards) is kind of "broken" and gradually transition into the hole dispersion (opening downwards) across the Fermi momentum with the spectral weight fading away as we go deeper into the hole band. This is typically what can be seen in the ARPES experiment. Now the density of state can be obtained by simply collapsing the spectral weight to the frequency space (by integrating out the momentum), which gives something like the following (let us say in 2D),

We can observe the BCS gap of $2\Delta$ across the Fermi energy ($\omega=0$). Two superconducting coherence peaks appear on the gap edges, as electrons are depleted in the gap and piled up to the coherence peaks. In this example, the gap has a U-shape because the Fermi surface is fully gapped. In $d$-wave superconductors like cuprates where the spectrum is nodal, there will be a V-shape gap instead in the density of state. This is what can be seen in the STM $\text{d}I/\text{d}V$ (differential conductance) spectrum.
