Validity of mean-field approximation In mean-field approximation we replace the interaction term of the Hamiltonian by a term, which is quadratic in creation and annihilation operators. For example, in the case of the BCS theory, where
$$
\sum_{kk^{\prime}}V_{kk^{\prime}}c_{k\uparrow}^{\dagger}c_{-k\downarrow}^{\dagger}c_{-k^{\prime}\downarrow}c_{k^{\prime}\uparrow}\to\sum_{k}\Delta_{k}c_{k\uparrow}^{\dagger}c_{-k\downarrow}^{\dagger} + \Delta_{k}^{\star}c_{-k\downarrow}c_{k\uparrow}\text{,}
$$ 
with $\Delta_{k}=\sum_{k^{\prime}}V_{kk^{\prime}}\langle c_{-k^{\prime}\downarrow}c_{k^{\prime}\uparrow}\rangle\in\mathbb{C}$. Then, in books, like this by Bruss & Flensberg, there is always a sentence like "the fluactuations around $\Delta_{k}$ are very small", such that the mean-field approximation is a good approximation. But we known for example in the case of the 1D Ising model the mean-field approximation is very bad.
My question: Is there a inequality or some mathematical conditions which says something about the validity of the mean-field approach? Further, is there a mathematical rigoros derivation of the mean-field approximation and the validity of it?
 A: Mean field theory is only good when fluctuations are small, which means that the free energy of a fluctuation must be much smaller than the total free energy.
The free energy of the typical fluctuation is of order $kT$ and its size is determined by the correlation length $\xi$, and is of order $\xi^d$, with $d=$ dimension:
$$F_{fluct}\sim \frac{kT}{\xi^d}\sim \mid t\mid^{\nu d}$$
Having used $\xi \sim \mid t \mid^{-\nu}$ where $t=(T-T_c)/T_c$ and $T_c$ is the critical temperature.
To obtain the total free energy we must integrate twice the specific heat, $c\sim \mid t \mid^{-\alpha}$. Imposing that $F_{fluct}/F$ goes to $0$ for $t\to 0$ we obtain
$$d \nu > 2-\alpha$$
For example in the Ising model we have $\alpha=0$ and $\nu=1/2$, so that the condition is $d>4$. This is why the mean field approximation is bad for the Ising Model in less than four dimensions.
A: You can introduce the ``would-be'' bosonic mean field exactly, using the Hubbard-Stratonich (a.k.a partial bosonization) method, see wikipedia and Interacting fermions on a lattice and Hubbard-Stratonovich transformation and mean-field approximation .
The mean field approximation correspond to performing the integral over the bosonic field using the stationary phase approximation. The fermion action is bilinear and can be performed exactly. Corrections to the stationary phase approximation correspond to fluctuations around the mean field. These are small if the four-fermion coupling in the BCS theory is small. If the coupling is strong it may be possible to justify the mean field approximation by a large N (N component vector field) or large d (number of dimensions)
approximation. 
Schmatically, the effective action of the bosonized theory is 
$$
 S = {\rm Tr}\{\log[G_0^{-1}G(\phi)]\} -\frac{\phi^2}{g}
$$
where $g$ is the coupling, $\phi$ is the bosonic field, and 
$$
G(\phi) = \left(
\begin{matrix}
p_0-\epsilon_p  & \phi \\ 
\phi    & p_0+\epsilon_p
\end{matrix}
\right)
$$
is the propagator. I also define $G_0=G(0)$. Now 
$$
\left. \frac{\delta S}{\delta\phi}\right|_{\phi_0} = 0
$$
is the MFA gap equation 
$$
\phi_0 = g\int \frac{d^3p}{(2\pi)^3}
\frac{\phi_0}{\sqrt{\epsilon_p^2+\phi_0^2}} .
$$
Corrections can be found by expanding $S$ around
$\phi=\phi_0+\delta \phi$. This will give higher loops containing $G(\phi_0)$. The expansion parameter is g. In physical units $1/g$ is the logarithm of the Fermi energy over the gap, $g\sim [\log(E_F/\phi_0)]^{-1}$.
Near $T_c$ corrections to mean field (=Landau-Ginsburg) are controlled by the Ginsburg criterion, as explained in Valrio92's answer. In weak coupling BCS the Ginsburg window is small, and mean field is accurate, except very close to $T_c$. 
A: I have recently noticed the reference given by @akhmeteli, which greatly helps me understand MFA. I shall do some math here (with respect to the original problem):
$$
\newcommand{\expect}[1]{\langle #1 \rangle}
\newcommand{\Tr}{\text{Tr}}
H = \sum \epsilon_k c^\dagger_k c_k + \sum_{k k'} V_{kk^{\prime}}c_{k\uparrow}^{\dagger}c_{-k\downarrow}^{\dagger}c_{-k^{\prime}\downarrow}c_{k^{\prime}\uparrow}
$$
MFA, as suggested by Bogoliubov, can be understand completely within variational principle, therefore one may get rid of the confusing "order parameter" and "fluctuation", with the following procedure:

*

*step 1: choose one non-interacting Hamiltonian with as many adjustable parameters as you wish, say
$$
 H_0 = H_0(E, \Delta) = \sum E_k c^\dagger_{k\sigma} c_{k\sigma} + \sum_{k, \sigma\sigma'} (\Delta_{k, \sigma\sigma'}c_{k\sigma}^{\dagger}c_{- k\sigma'}^{\dagger} + \text{H.C.})
 $$
note that the form here implies the spatial translation symmetry. You can of course using a more general form, which may give you more observations (say the density wave state which breaks translational symmetry)


*step 2: Bogoliubov inequality helps one to find the optimal parameters set, defining:
$$
 F(H) = -\frac{1}{\beta} \log \left(\Tr(e^{-\beta H})\right)
 $$
$$
 \expect{\hat{O}}_0 = \frac{\Tr( \hat{O} e^{-\beta H_0})}{\Tr( e^{-\beta H_0})}
 $$
the inequality (variational principle):
$$
 F(H) \leq F(H_0) + \expect{H-H_0}_0, \quad \forall H_0
 $$
optimal $H_0$ is thus given by:
$$
 \min_{\text{all parameters}} F(H_0) + \expect{H-H_0}_0
 $$
with some additional math, especially the Feynman-Hellmann equation:
$$
\begin{aligned}
\frac{\partial F(H_0)}{\partial E_k} &=  \sum_{\sigma} \expect{c^\dagger_{k\sigma}c_{k\sigma}}_0 \\
\frac{\partial F(H_0)}{\partial \Delta_{k,\sigma\sigma'}} &=  \expect{c^\dagger_{k\sigma} c^\dagger_{-k,\sigma'}}_0 \\
& \cdots
\end{aligned}
$$
and the Wick's theorem for expectation calculated on a non-interacting ground state:
$$
\expect{c_1^\dagger c_2^\dagger c_3 c_4}_0 = \expect{\dagger\dagger}_0\expect{..}_0 - \expect{}_0\expect{}_0 + \expect{}_0\expect{}_0
$$
It shall be easy to restore your desired result.
A: Mean-field theory is exact (in the thermodynamic limit) in the case of long-range interaction (which is not the case for the nearest-neighbor Ising model). Therefore, mean-field theory is exact for BCS, where you have an effective long-range interaction. As for rigorous results, Bogoliubov rigorously proved that in the ground state (zero temperature) the energies per particle coincide for BCS and the mean-field theory in the thermodynamic limit. Later, Bogoliubov Jr. rigorously proved the same for arbitrary temperatures (and free energies per particle). See the bibliography, e.g., at http://arxiv.org/abs/1507.00563 (the discussion is around page 43).
