Ignoring $(\sigma_i-M)(\sigma_j-M)$ in mean field theory? A way to do mean field theory for the Ising model is as follows.



*

*First take the Ising Hamiltonian: $$H=-J \sum_{\left<i,j\right>} \sigma_i\sigma_j$$

*Let $\sigma_i=\sigma_i-M+M$ and likewise for $\sigma_j$ to get: $$H=-J \sum_{\left<i,j\right>} (M^2 +(\sigma_i-M)M+(\sigma_j-M)M+\underbrace{(\sigma_i-M)(\sigma_j-M)}_{\bigstar})$$

*Ignore the stared ($\star$) term.

*Write down the partition function, apply a self-consistency condition etc.

Given that in the Ising model $\sigma_i=\pm1$ thus for any given $i$ and $j$, the ($\star$) term is not going to be small. What is the standard justificiation for then ignoring it?
 A: Even though $\sigma_i-M$ is not small, expectation value of its square is small as that corresponds to the variance, hence fluctuations, which are assumed to be next order terms in the mean-field-theory. That's why summation $\sum (\sigma_i-M)(\sigma_j-M)$, which is basically autocovariance function along lattice sites, should be small. By the way, here we should have $M=<\sigma_i>$.
A: I could not find any references that explains this in detail so if you think it is wrong please let me know.
This issue is a lot more subtle then it first appears - especially in the case where the external magnetic field, $h=0$, like I have put in the question.
Why is it a subtle issue? Well for $h=0$ our Hamiltonian has symmetry under $\sigma_i \leftrightarrow -\sigma_i$, meaning for any configuration $\{\sigma_i\}$ the configuration $\{-\sigma_i\}$ has exactly the same energy and therefore contributes the same waiting to the partition function.  If the $\star$ term is small for $\{\sigma_i\}$ then it clearly won't be for $\{-\sigma_i\}$ (in general it will be large and positive).
To see this think about the case where $M$ is close to $+1$, for the configuration $\{\sigma_i=1\;\forall i\}$ this last term will indeed be small but for $\{\sigma_i=-1\;\forall i\}$ it won't be. 
Thus for any case except $M=0$ ignoring the $\star$ term dramatically decreases the probability of the flipped spins for the case of $h=0$. 
Therefore to make the mean field theory valid in this case and for the argument as presented in Soner's answer to work we must a spontaneous symmetry breaking field present (or at least have it's present in the back of our minds). This means that the inverted configurations already have much smaller probability then those with mean value near $M$ and as such ignoring $\star$ does not have a big effect on our partition function. 
