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?
 A: I will try to answer the question : Why do we have to enlarge the number
of possibilities for the ground state ? I guess this is at the heart of
your problem. Before doing so, let us quickly discuss the notion of ground state itself, which are called Fermi sea in the case of fermion.
At the beginning: the Fermi surface or quasiparticles
So at the beginning is a collection of fermions, and the notion of Fermi
surface. An important element of the forthcoming discussion is the notion of quasiparticle. In condensed matter systems, the Fermi surface is not built from free electrons, or bare electrons, or genuine electrons. To understand how this notion comes in, let's take a free electron, i.e. a particle following the Dirac equation with charge $e$ and spin $1/2$, and put it into any material (say a semi-conductor or a metal for instance). By electrons any condensed matter physicist means that in a complex system (s)he studies, it is impossible to take into account all the possible interactions acting on the bare electron. Most of the interaction will be of bosonic nature (think especially about phonons). So one hopes that taking all the complicated bosonic interactions on the fermionic bare electron results in a composite particle which still has a fermionic statistics, and hopefully behaves as a Schrödinger equation, possibly with kinetic energy $E_{c}=\dfrac{p^{2}}{2m_{2}^{*}}+\dfrac{p^{4}}{4m_{4}^{*}}+\cdots$ with effective masses $m_{2}^{*}$ that one can suppose $p$-independent and possibly $m_{2}^{*}\ll m_{4}^{*}$, such that one has the usual Schrödinger equation.
Note that 


*

*the Schrödinger equation is anyways invariant with respect to the statistics of the particle it describes

*the above construction can be made a bit more rigorous using the tools of effective field theory and renormalisation  


The important thing is that these electrons are still fermions, and so they pile up to form a Fermi sea. From now on we write electrons or quasiparticle without making distinction, since in condensed matter there is only quasiparticle. A quasiparticle of charge $e$ and spin $1/2$ will be called an electron. A Fermi surface is a stable object, as has been discussed in an other
question. Stable ? Well, not with respect to the Cooper mechanism,
which allows bound states of electrons to be generated on top of a Fermi
sea. These bound states are of charge $2e$ and their total spin makes
them some kind of bosons. They are as well quasiparticles, but we will call them Cooper pairs, instead of quasiparticle of charge $2e$ and spin $1$ or $0$.
Now we identified the ground state of a metal as a Fermi sea of fermionic quasiparticles called electrons, we can try to understand how this ground state becomes unstable and why we should then take into account several ground states, of possibly different statistical nature, as the bosonic versus fermionic ground state in a superconductor. The reason why we need to enlarge the available ground states is clearly
due to the symmetry breaking, as we review now.
Symmetry breaking and space of ground states
First, think about the para-ferromagnetic transition. Before the
transition (paramagnetic phase) you can choose the orientation of the
spin the way you want: they are random in the electron gas. A nice
picture about that is to say: the Fermi surface is just the same for all
the electrons.
Now comes the ferromagnetic phase: the system chooses either to align
all the spins in the up or down direction (of course the direction is
not fixed and the system is still rotationally invariant unless a
magnetic field is applied, but clearly the electron spins are
polarised). What about the Fermi surface ? Well, it becomes two-fold ...
There is now one Fermi surface for the spin-up electrons and one Fermi
surface for the spin-down electron. So the number of available ground
states increases.
The link to the symmetry breaking is clear: the more you want of
symmetries, the less possible states you allow. Say in the other way:
breaking the symmetry allows for more states to exist. This is also
quite straightforward from the following argument: once you allow an
interaction responsible for the parra-to-ferro transition, you must
first answer the question: is the unpolarised Fermi surface or one of
the polarised Fermi surfaces the true ground state ? So you need a way to
compare the spin-unpolarised and the spin-polarised ground states. So
clearly the number of accessible ground states must be greater once a
phase transition and a symmetry breaking is under the scope.
Now, about superconductivity and the relation to symmetry breaking, I
refer to this (about particle number conservation) and this
(about the $\text{U }\left(1\right)\rightarrow\mathbb{Z}_{2}$ symmetry
breaking in superconductors) answers. The important thing is that the
Cooper mechanism makes the Fermi surface instable. What results ? A kind
of Bose-Einstein condensate of charged particles (the Cooper pairs of
charge $2e$) and some electrons still forming a Fermi-Dirac condensate
and so a Fermi sea, with less fermions than before the transition (hence
the number of electrons is not conserved). So now the available ground
states are i) the genuine Fermi sea made of electrons, ii) the charged
Bose-Einstein condensate made with all the electrons paired up via
Cooper mechanism and iii) a mixture of the two Fermi-Dirac and
Bose-Einstein condensates (be careful, the terminology is misleading, a
Bose-Einstein condensate and a Cooper pair condensate are not really the
same thing).
Unfortunately, the real ground state is a kind of mixture, but at
zero-temperature, one might suppose that all the conduction electrons
have been transformed in Cooper pairs (in particular, this can not be
true if one has an odd number of electrons to start with, but let forget
about that). Let us call this complicated mixture the Cooper
condensate, for simplicity.
In any case, we have to compare the Fermi sea with the Cooper
pairs condensate. That's precisely what we do by supposing a term like $\left\langle
cc\right\rangle \neq0$.
A bit of mathematics
We define a correlation as $\left\langle c_{1}c_{2}\right\rangle $ with
creation or annihilation operators. In a paramagnetic phase, we have
$$\left\langle N\right|c^{\dagger}c\left|N\right\rangle =\left\langle
n_{\uparrow}\right|c_{\uparrow}^{\dagger}c_{\uparrow}\left|n_{\uparrow}\right\rangle +\left\langle
n_{\downarrow}\right|c_{\downarrow}^{\dagger}c_{\downarrow}\left|n_{\downarrow}\right\rangle $$ as the only non-vanishing correlation, with
$\left|N\right\rangle$ a Fermi sea filled with $N$ electrons. Note in
that case the polarised Fermi seas
$\left|n_{\uparrow,\downarrow}\right\rangle$ make no sense, since there
is no need for an internal degree of freedom associated to the
electrons ; or these two seas are the two shores of the unpolarised Fermi ocean... Now, in the ferromagnetic phase, the correlations
$\left\langle
n_{\uparrow}\right|c_{\uparrow}^{\dagger}c_{\uparrow}\left|n_{\uparrow}\right\rangle$ and $\left\langle
n_{\downarrow}\right|c_{\downarrow}^{\dagger}c_{\downarrow}\left|n_{ \downarrow}\right\rangle$ starts to make sense individually, and the ground states
with polarised electrons as well. In addition, we must compare all of
these inequivalent ground states. One way to compare all the possible
ground states is to construct the matrix $$\left\langle
\begin{array}{cc} c_{\uparrow}c_{\downarrow}^{\dagger} &
c_{\uparrow}c_{\uparrow}^{\dagger}\\
c_{\downarrow}c_{\downarrow}^{\dagger} &
c_{\downarrow}c_{\uparrow}^{\dagger} \end{array}\right\rangle
=\left\langle \left(\begin{array}{c} c_{\uparrow}\\ c_{\downarrow}
\end{array}\right)\otimes\left(\begin{array}{cc}
c_{\downarrow}^{\dagger} &
c_{\uparrow}^{\dagger}\end{array}\right)\right\rangle $$ where the
ground state is a bit sloppily defined (i.e. I did not refer to
$\left|n_{\uparrow,\downarrow}\right\rangle$, and I simply put the
global $\left\langle \cdots\right\rangle $ for simplicity). The fact
that the construction is a tensor product (the symbol $\otimes$ in the
right-hand-side just does what appears in the left-hand-side) clearly
shows that you can restaure the different ground states as you wish. In
a sense, the problem of defining the different ground states is now put
under the carpet, and you just have to deal with the above matrix.
Clearly the diagonal elements exist only in the paramagnetic phase and
the off-diagonal elements exists only in the ferromagnetic case, but
this is no more a trouble, since we defined a tensorial product of
several ground states and we are asking: which one is the good one?
Now, for the superconductor, one does not polarise the spins of the
electrons, one creates some bosonic correlations on top of two fermionic
excitations. So the natural choice for the matrix is $$\left\langle
\begin{array}{cc} cc^{\dagger} & cc\\ c^{\dagger}c^{\dagger} &
c^{\dagger}c \end{array}\right\rangle =\left\langle
\left(\begin{array}{c} c\\ c^{\dagger}
\end{array}\right)\otimes\left(\begin{array}{cc} c^{\dagger} &
c\end{array}\right)\right\rangle $$ where still, the diagonal part
exists already in a simple metal, and the off-diagonal shows up once the
system transits to the superconducting phase. Clearly, if you call
$\left|N\right\rangle $ the Fermi-Dirac condensate with $N$ electron of
charge $e$, and $c$ the operator destroying an electron in this Fermi
sea, then you must define $\left\langle cc\right\rangle \equiv\left\langle N-2\right|cc\left|N\right\rangle \propto\left\langle N\right|cc\left|N+2\right\rangle \propto\left\langle N-1\right|cc\left|N+1\right\rangle $ (note you can do the way you want, and this has profound implications for Josephson physics, but this is not the story today) for the correlation to exist. What are the states $\left|N-2\right\rangle $ then ? Well, it is clear from the context: a Fermi sea with two electrons removed. Without the Cooper mechanism, we would have no idea what is this beast, but thanks to him, we know this is just the instable Fermi sea with one Cooper pair removed due to the action of a virtual phonon.
Whereas the parra-ferromagnetic transition was seen in the competition
between unpolarised versus polarised spins, the normal-superconducting
transition can be seen in the competition between particle and hole
versus particle-hole mixtures.
About the mean-field construction
Now, how do we construct the mean-field Hamiltonian ? We simply use the
two body interaction term, and we apply the Wick's theorem. That is,
one does $$\left\langle
c_{1}c_{2}c_{3}^{\dagger}c_{4}^{\dagger}\right\rangle =\left\langle
c_{1}c_{2}\right\rangle \left\langle
c_{3}^{\dagger}c_{4}^{\dagger}\right\rangle -\left\langle
c_{1}c_{3}^{\dagger}\right\rangle \left\langle
c_{2}c_{4}^{\dagger}\right\rangle +\left\langle
c_{1}c_{4}^{\dagger}\right\rangle \left\langle
c_{2}c_{3}^{\dagger}\right\rangle $$ valid for any average taken over
Gaussian states. Clearly, one has (replace the numbers by spin vectors
eventually) : the Cooper pairing terms, the Heisenberg-ferromagnetic
coupling and some anti-ferromagnetic coupling (not discussed here).
Usually, since a system realises one ground state, we do not need to try
all the different channels. For superconductivity we keep the first term
on the right-hand-side.
As explained in many other pages, the mean-field treatment can be
justified carefully in the case of conventional superconductivity. See
e.g. this, or this questions, or this answer (on this
page).
A: 
First, why do we choose those terms to put the brakets around? Do we
  somehow know that these are the terms that won't fluctuate much?

The heuristic motivation is that the superconductivity phenomenon can be understood as the non-relativistic analog of the Higgs mechanism for the EM $U_{\text{EM}}(1)$ (this idea explains the main properties of the superconductor, such as absence of resistivity, magnetic field expulsion and other). For the Higgs mechanism one need a condensate of some particles. In the solid body, the only candidates are electrons. The simplest choice of the condensate is the scalar condensate (very-very heuristically, we require the Galilei invariance into the superconductor). Therefore the only allowed simplest scalar bounded state $|\psi\rangle$ is constructed from the electron pair with zero total spin, which can be imagine as two fermions moving in opposite directions:
$$
|\psi\rangle = c^{\dagger}_{\mathbf k, \uparrow}c^{\dagger}_{-\mathbf k, \downarrow}|0\rangle
$$
This is indeed realized when we calculate the four-fermion vertex in the microscopic theory of electron-phonon interactions. Indeed, we start from the interaction Hamiltonian
$$
H_{\text{int}} = g\psi^{\dagger}\psi \varphi, 
$$ 
where $\psi$ is electron field, while $\varphi$ is phonon field, and assume the second order in $g$ diagrams of electron-electron scattering; let's denote their momenta as $\mathbf p_{i}$. After tedious calculations of the corresponding diagram one finds that the vertex $\Gamma^{(4)}(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4})$ (with $\mathbf{q} = \mathbf{p}_{1} + \mathbf{p}_{2}$) has the pole at zero momentum $\mathbf q$ and zero total spin of fermions $1,2$ and $3,4$ or ($1,4$ and $2,3$). The pole quickly disappears once one enlarges $|\mathbf q|$. 
Finally, since 4-fermion vertex is expressed through two-particles Green $G^{(2)}$ function, then the presence of the pole in $\Gamma^{(4)}$ means the presence of the pole in $G^{(2)}$. This means that bounded states of electrons appear. 

Second, why can't we just immediately say $\langle c_{-k'\downarrow}c_{k'\uparrow}\rangle=\langle c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger\rangle=0$ since we know the ground state will
  have definite particle number, and $c_{-k'\downarrow} c_{k'\uparrow}$
  and $c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger$ don't conserve
  particle number?

(Added) Let's again use the idea of Higgs mechanism. It is clearly that the number of particles isn't strictly conserved in superconducting state. Really, because of the SSB of $U_{\text{EM}}(1)$ group the system loses the definition of the conserved particle number (being associated to the global phase invariance current) in the ground state.
Microscopically the formation of the condensate is argued in a following way: the presence of the two-fermions bounded states leads to instability of the electron gas inside the superconductor, and fermions are converted in Cooper pairs. Really, for bosons the lowest energy state can be filled by infinitely many number of particles, while for fermions this is not true due to Pauli principle. This (under some assumptions, see below) leads to formation of new vacuum state on which the operator $\hat{c}^{\dagger}_{\mathbf k,\uparrow}\hat{c}^{\dagger}_{-\mathbf k, \downarrow}$ has non-zero VEV:
$$
\langle \text{vac} | \hat{c}^{\dagger}_{\mathbf p,\uparrow}\hat{c}^{\dagger}_{-\mathbf p, \downarrow}|\text{vac}\rangle = \Psi \neq 0
$$
It's not hard to construct this state:
$$
|\text{vac}\rangle = \prod_{\mathbf k = \mathbf k_{1}...k_{\frac{N}{2}}}(u_{\mathbf k} + v_{\mathbf k}\hat{c}^{\dagger}_{\mathbf k,\uparrow}\hat{c}^{\dagger}_{-\mathbf k, \downarrow})|0\rangle, \quad |v_{\mathbf k}|^{2} + |u_{\mathbf k}|^{2} = 1
$$
This state called the coherent state really doesn't have definite number of particles. 
For the particles number operator $\hat{N} = \sum_{\mathbf k}\hat{c}^{\dagger}_{\mathbf k,\uparrow}\hat{c}^{\dagger}_{-\mathbf k, \downarrow}$ the VEV $\langle \text{vac}|\hat{N}|\text{vac}\rangle$ is
$$
\tag 1 N = \langle \text{vac}|\hat{N}|\text{vac}\rangle \sim \sum_{\mathbf k}|v_{\mathbf k}|^{2},
$$
Why the number of particles is not so indefinite
However, the root $\sqrt{\Delta N^{2}}$ of quadratic deviation from
while the quadratic deviation $\langle \text{vac}|(\hat{N} - N)^{2}|\text{vac}\rangle$ is
$$
\tag 2 \Delta N^{2} = \langle \text{vac}|(\hat{N} - N)^{2}|\text{vac}\rangle = \sum_{\mathbf k}|u_{\mathbf k}|^{2}|v_{\mathbf k}|^{2}
$$
is negligible in the limit of large particle number (and volume). Really, both of $(1)$ and $(2)$ are proportional to the volume $V$ (note only one summation over wave-numbers $\mathbf k$ in $(2)$), and hence $\sqrt{\Delta N^{2}} \sim \sqrt{V}$ which is negligible in the large volume limit. This argument is seemed to be first published in the BCS article.
(Added) This means that the state of superconductor with good accuracy can be projected on the state with definite particle number $N$. Precisely, define the phase dependent ground state
$$
|\text{vac}(\theta)\rangle = |\text{vac}[c^{\dagger}_{\mathbf k, \uparrow/\downarrow} \to c^{\dagger}_{\mathbf k, \uparrow/\downarrow}e^{-i\theta}]\rangle
$$ 
The state with definite particle number $N$ is defined as
$$
|N\rangle = \int \limits_{0}^{2\pi} d\theta e^{i\hat{N}\theta}|\text{vac}(\theta)\rangle
$$
A: After some thoughts, I think an other (and hopefully better) answer might go along the following lines, perhaps more rigorous. Nevertheless, I will not give the full mathematical details, since they require too much writing.
At the heart of perturbative treatment of many-body problems lies the Wick theorem. For two body interaction of fermionic nature, it states that 
$$\left\langle c_{1}^{\dagger}c_{2}^{\dagger}c_{3}c_{4}\right\rangle =\left\langle c_{1}^{\dagger}c_{2}^{\dagger}\right\rangle \left\langle c_{3}c_{4}\right\rangle -\left\langle c_{1}^{\dagger}c_{3}\right\rangle \left\langle c_{2}^{\dagger}c_{4}\right\rangle +\left\langle c_{1}^{\dagger}c_{4}\right\rangle \left\langle c_{2}^{\dagger}c_{3}\right\rangle $$
As any theorem it can be rigorously proven. In fact, the proof is not that cumbersome, and I refer to 


*

*Danielewicz, P. (1984). Quantum theory of nonequilibrium processes, I. Annals of Physics, 152(2), 239–304. See Appendix A.

*Gaudin, M. (1960). Une démonstration simplifiée du théorème de Wick en mécanique statistique. Nuclear Physics, 15, 89–91. In french.
where the demonstration is the same in both papers. The key point is that one needs a Gaussian state to demonstrate the theorem, namely one needs a statistical operator in the form $\rho=e^{-\sum\epsilon_{n}c_{n}^{\dagger}c_{n}}$. In any other case the Wick's decomposition does not work (or at least I'm not aware of such a decomposition). Importantly, the demonstration in the above references uses only the anti-commutation properties between the $c$'s operators. So any time one has a statistical average made over Gaussian states of some anti-commuting (i.e. fermionic) operators, the Wick's theorem applies.
Now, the strategy of the mean-field treatment is the following. Write the Hamiltonian (I use the symbol $\sim$ to withdraw integrals and/or sums over all degrees of freedom)
$$H\sim H_{0}+c_{1}^{\dagger}c_{2}^{\dagger}c_{3}c_{4} \sim H_{0}+H_{\text{m.f.}}+\delta H$$
with $H_{0}$ the free particle Hamiltonian (which can be written in the form $H_{0}\sim\epsilon_{n}c_{n}^{\dagger}c_{n}$ in a convenient basis), $H_{\text{m.f.}}=\Delta c_{1}^{\dagger}c_{2}^{\dagger}+\Delta^{\ast}c_{3}c_{4}$ the mean-field Hamiltonian, and $\delta H=c_{1}^{\dagger}c_{2}^{\dagger}c_{3}c_{4}-H_{\text{m.f.}}$ the correction to the mean-field Hamiltonian.
The mean-field Hamiltonian can be diagonalised using a Bogliubov transformation, namely, one can write
$$H_{0}+H_{\text{m.f.}}=\sum_{n}\epsilon_{n}\gamma_{n}^{\dagger}\gamma_{n}$$
where the $\gamma$'s are some fermionic operators verifying 
$$c_{i}=u_{ik}\gamma_{k}+v_{ik}\gamma_{k}^{\dagger}\\c_{i}^{\dagger}=u_{ik}^{\ast}\gamma_{k}^{\dagger}+v_{ik}^{\ast}\gamma_{k}$$
note: there are some constraints imposed on the $u$'s and $v$'s in order for the above Bogoliubov transformation to preserve the (anti-)commutation relation ; in that case one talks about canonical transformation, see 


*

*Fetter, A. L., & Walecka, J. D. (1971). Quantum theory of many-particle systems. MacGraw-Hill.


for more details about canonical transformations. Once one knows a few properties of the mean-field ground state, one can show that $\lim_{N\rightarrow\infty}\left\langle \delta H\right\rangle \rightarrow0$ when $N$ represents the number of fermionic degrees of freedom, and the statistical average $\left\langle \cdots\right\rangle $ is performed over the mean-field ground state, namely $\left\langle \cdots\right\rangle =\text{Tr}\left\{ e^{-\sum\epsilon_{n}\gamma_{n}^{\dagger}\gamma_{n}}\cdots\right\} $. See 


*

*De Gennes, P.-G. (1999). Superconductivity of metals and alloys. Advanced Book Classics, Westview Press


for a clear derivation of this last result. In fact deGennes shows how to choose the $u$ and $v$ in order for the Bogoliubov transformation to diagonalise the mean-field Hamiltonian, and then he shows that this choice leads to the best approximation of the zero-temperature ground state of the interacting (i.e. full) Hamiltonian. The idea to keep in mind is that the mean-field treatment works. Up to now, we nevertheless restrict ourself to operator formalism. When dealing with quantum field theory and its methods, one would prefer to use the Wick's theorem.
So, now we come back to the Wick's decomposition. One can manipulate the statistical average and show that 
$$\left\langle c_{1}^{\dagger}c_{2}^{\dagger}\right\rangle =\text{Tr}\left\{ \rho c_{1}^{\dagger}c_{2}^{\dagger}\right\} \\=\text{Tr}\left\{ e^{-\varepsilon_{n}c_{n}^{\dagger}c_{n}}c_{1}^{\dagger}c_{2}^{\dagger}\right\} =\text{Tr}\left\{ c_{1}^{\dagger}c_{2}^{\dagger}e^{-\varepsilon_{n}c_{n}^{\dagger}c_{n}}\right\} =e^{-2\varepsilon_{n}}\text{Tr}\left\{ e^{-\varepsilon_{n}c_{n}^{\dagger}c_{n}}c_{1}^{\dagger}c_{2}^{\dagger}\right\} \\=e^{-2\varepsilon_{n}}\left\langle c_{1}^{\dagger}c_{2}^{\dagger}\right\rangle =0$$
using the cyclic property of the trace and the Bakker-Campbell-Hausdorf formula to pass the exponential from right to left of the two creation operators. The quantity should be zero because it is equal to itself multiplied by a positive quantity $e^{-2\varepsilon_{n}}$. The manipulation is the same as one does for the Wick's theorem, so I do not write it fully, check the above references.
The conclusion is that clearly, the quantity $\left\langle c_{1}^{\dagger}c_{2}^{\dagger}\right\rangle $ is zero in the Wick's decomposition. But due to the Cooper problem (or what we learned from the mean-field treatment), one knows that we should not take the statistical average over the free fermions $e^{-\sum\varepsilon_{n}c_{n}^{\dagger}c_{n}}$, but over the Bogoliubov quasi-particles $e^{-\sum\epsilon_{n}\gamma_{n}^{\dagger}\gamma_{n}}$, since they represent the approximately true (in the limit of large $N$ it is exact) ground state. Say differently, the Fermi surface is not the ground state of the problem, and the statistical average should not be taken over this irrelevant ground state. We should choose instead the statistical average over the Bogoliubov's quasiparticles $\gamma$. Namely, one has 
$$\left\langle c_{1}^{\dagger}c_{2}^{\dagger}\right\rangle =\text{Tr}\left\{ e^{-\varepsilon_{n}\gamma_{n}^{\dagger}\gamma_{n}}c_{1}^{\dagger}c_{2}^{\dagger}\right\} \neq0$$
and this quantity is not zero, as you can check when writing the $c$'s in term of the $\gamma$'s. Importantly, note that the Wick's theorem can be prove rigorously for this statistical operator (see above).
This complete the proof that the mean-field approximation can be made rigorous, and that one can take the quantity $\left\langle c_{1}^{\dagger}c_{2}^{\dagger}\right\rangle$ as the order parameter of a phase whose ground state is not filled with free electrons, but with Bogoliubov quasi-particles generated by the $\gamma$ operators. One usually refers to this ground state as the Cooper sea. If you expand the exponential with the $\gamma$'s over the zero-electron state noted $\left|0\right\rangle $, you will see you end up with something like $e^{-\sum\epsilon_{n}\gamma_{n}^{\dagger}\gamma_{n}}\left|0\right\rangle \sim\prod\left(u+vc_{1}^{\dagger}c_{2}^{\dagger}\right)\left|0\right\rangle $ which is the BCS Ansatz. Using rigorous notations you can make a rigorous proof of this last statement.
Now if you want to make the derivation completely rigorous, it requires to demonstrate the Wick's theorem, to manipulate the mean-field Hamiltonian in order to show the Bogoliubov mean-field transformation works quite well, and to calculate the statistical average over the Bogoliubov ground state. This goes straightforwardly ... over many many pages which I'm too lazy to write here.
