Validity of mean-field approximation

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.

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).

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$.