Correlations in Ising mean-field theory I am reading the book "Critical Dynamics - A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior" (sections 1.1.2 and 1.1.3) and have been somewhat confused about the meaning of a mean-field theory (related questions: 1, 2).
To summarise, in deriving the Ising mean-field theory, the spin variable at site $i$ is written as
$$\sigma_i = m + (\sigma_i - \langle \sigma_i \rangle),$$
where $m$ represents the mean magnetisation, and then the Ising Hamiltonian,
$$ H = -\frac{1}{2} \sum_{i,j}J \sigma_i \sigma_j - h \sum_i \sigma_i$$
($h$ is a magnetic field) is expanded by neglecting the term $(\sigma_i - \langle \sigma_i\rangle)(\sigma_j - \langle \sigma_j \rangle)$. Then the partition function is computed from which a consistency equation for $m$ is derived, namely
$$ m = \tanh \big(\beta (h+J m) \big).$$
So far so good.
Next, the author proceeds to compute the correlation function via the response function by considering a site-dependent magnetic field by writing
$$\chi_{ij} = \frac{\partial m_i}{\partial h_j} = \beta \langle (\sigma_i - \langle \sigma_i\rangle)(\sigma_j - \langle \sigma_j \rangle )\rangle.$$
The response function $\chi$ on the lhs can be computed by taking appropriate derivatives from the consistency relation for $m$, from which the so-called Ornstein-Zernicke form for the correlation funciton on the rhs is obtained. This form is also later obtained from the coarse-grained Landau free energy.
It appears to me that this is an inconsistent approach, since we first ignore the correlations between spatial fluctuations (by omitting the $(\sigma_i - \langle \sigma_i \rangle)(\sigma_j - \langle \sigma_j \rangle)$ ) in order to obtain the partition function, from which free energy and $m$ are evaluated; and then we go back to compute the correlations that we had assumed to be zero and get a non-zero value for it!
What am I missing here?
PS: I also wonder whether "mean-field" theories, in general, should fully ignore any correlations between the variables, or do they still account for the correlations in some sense?
 A: You are not missing anything. There is an inconsistency, but it is not in the approach.
A point that is not frequently appreciated is that the whole machinery of Statistical Mechanics formalism is based on the assumption that the Hamiltonian is independent of the thermodynamic state. Many statistical mechanics formulae establish relations between averages and response functions, i.e., derivatives, assuming that the Hamiltonian depends only on the microscopic degrees of freedom. For example, a simple formula in the canonical ensemble like
$$
\langle H\rangle = -\frac{\partial}{\partial \beta} \log Z,
$$
where $Z$ is the canonical partition function and $\beta=\frac{1}{k_BT}$, ceases to be valid if the Hamiltonian would depend on $\beta$.
However, some macrostate dependence of the Hamiltonian is the unavoidable result of the Mean Field approximation, whatever is its formal introduction (self-consistent relations, minimum of approximate free energy, variational approach).
Consequently, there is a built-in inconsistency in many formulae of Statistical Mechanics we are used to considering equivalent. Such inconsistency is only due to the Mean Field approximation, and would not be there in an exact treatment of a honest, state-independent Hamiltonian.
In practice, once we realize that, due to the MF approximation, different routes to the same physical quantity are not equivalent anymore, we may investigate which among them perform better to approximate the exact results.
That is the case, for example, of the spin-spin correlation function. On the one hand, Statistical Mechanics shows that it may be evaluated either as averages ( $ \chi_{ij}= \langle S_i S_j \rangle  -  \langle S_i \rangle  \langle S_j \rangle$ ), or from the response function of the local magnetization to a change of a local magnetic field ( $ \chi_{ij}=\frac{\partial{m_i}}{h_j} $ ). On the other hand, the two formulae are not equivalent in MF. The former is trivially zero, while the latter is different from zero. The implicit assumption, which is done by taking the second formula as representative of the spin-spin correlation, is that, among two approximate formulae, one uses the non-zero one as the best approximation of the exact result.
