quantum mechanics - Mean field equations in the BCS theory of superconductivity

In BCS theory, one takes the model Hamiltonian

$$\sum_{k\sigma} (E_k-\mu)c_{k\sigma}^\dagger c_{k\sigma} +\sum_{kk'}V_{kk'}c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger c_{-k'\downarrow} c_{k'\uparrow}$$

This Hamiltonian clearly conserves particle number. Thus, we expect the ground state to have a definite particle number. It's possible the ground state is degenerate, but that could be lifted by perturbing $\mu$.

Then, one makes a mean field approximation. One replaces

$$c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger c_{-k'\downarrow} c_{k'\uparrow}$$

with

$$\langle c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger\rangle c_{-k'\downarrow} c_{k'\uparrow}+c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger \langle c_{-k'\downarrow} c_{k'\uparrow}\rangle-\langle c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger\rangle \langle c_{-k'\downarrow} c_{k'\uparrow}\rangle$$

This doesn't make any sense to me. This seems to be saying that we know the terms $c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger$ and $c_{-k'\downarrow} c_{k'\uparrow}$ don't fluctuate much around their mean values. But we also know that in the actual ground state, the mean values are given by $\langle c_{-k'\downarrow} c_{k'\uparrow}\rangle=\langle c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger\rangle=0$ since the ground state will have definite particle number. Thus the fluctuations about the mean value aren't small compared to the mean value.

How can this mean field treatment be justified?