Sum of all two point correlation functions in the Ising model The two point correlation function for the Ising model is defined as $\left[\langle\sigma_i\sigma_j\rangle  -\langle\sigma_i\rangle\langle\sigma_j\rangle\right]$. Then the sum over $i$ $j$ of that function gives:
\begin{equation}
\sum_{ij}\left[\langle\sigma_i\sigma_j\rangle  -\langle\sigma_i\rangle\langle\sigma_j\rangle\right] = \langle M^2\rangle - \langle M\rangle^2
\end{equation}
Where $\sigma_i = \pm 1$ and $M  =\sum_i \sigma_i$. The expression can be normalized by dividing it by $M^2 = N^2$ where $N$ is the number of spins of the system. My question is, I am assuming I can justify the long range correlations in the lattice by computing the sum of the correlation functions, is that correct? I should expect it to vanish when most of the two point correlation vanish, and to have  peak when most two point correlation reach its maximum value.
 A: Let us be a bit more precise, so that the question can be answered accurately.
Let us thus assume that your system is composed of $N^d$ spins attached to the vertices of the box $\Lambda_N=\{1,\dots,N\}^d$. I am assuming that there is a positive magnetic field $h>0$ acting on the spins (this is technically useful to ensure that we are considering the proper state below, but the field will soon be set equal to $0$). Let me denote by $\langle\cdot\rangle_{N,\beta,h}$ the corresponding expectation.
Then, the susceptibility is given by
$$
\chi_N(\beta,h) = 
\frac1{N^d}\sum_{i,j\in\Lambda_N} \Bigl[ \langle \sigma_i\sigma_j \rangle_{N,\beta,h} - \langle \sigma_i \rangle_{N,\beta,h} \langle \sigma_j \rangle_{N,\beta,h} \Bigr].
$$
We are really interested in the thermodynamic limit of this quantity,
$$
\chi(\beta,h) = \lim_{N\to\infty} \chi_N(\beta,h) =
\sum_{i\in\mathbb{Z}^d} \Bigl[ \langle \sigma_0\sigma_i \rangle_{\beta,h} - \langle \sigma_0 \rangle_{\beta,h} \langle \sigma_i \rangle_{\beta,h} \Bigr],
$$
where I used the fact that the resulting infinite-volume state is translation invariant. We can now get rid of the magnetic field and define
$$
\chi(\beta) = \lim_{h\downarrow 0} \chi(\beta,h)
= \sum_{i\in\mathbb{Z}^d} \Bigl[ \langle \sigma_0\sigma_i \rangle_{\beta}^+ - \langle \sigma_0 \rangle_{\beta}^+ \langle \sigma_i \rangle_{\beta}^+ \Bigr],
$$
where $\langle \cdot \rangle_\beta^+$ denotes expectation with respect to the $+$ state.
It is known, in any dimension $d$, that the truncated 2-point function decays exponentially with the distance when $\beta\neq\beta_{\rm c}$. Namely, for all  $\beta\neq\beta_{\rm c}$, there exists $c=c(\beta,d)>0$ such that
$$
0\leq \langle \sigma_0\sigma_i \rangle_{\beta}^+ - \langle \sigma_0 \rangle_{\beta}^+ \langle \sigma_i \rangle_{\beta}^+ \leq e^{-c\|i\|}.
$$
(In other words, the correlation length is finite at all non-critical temperatures.) This immediately implies that the susceptibility is finite away from the critical point:
$$
\chi(\beta) < \infty\qquad\forall\beta\neq\beta_{\rm c}.
$$
Moreover, it is also known that, in all dimensions $d\geq 2$, the 2-point function does not decay exponentially fast when $\beta=\beta_{\rm c}$ (much more precise information is available when $d=2$ and when $d$ is large enough).
Finally, it is known that the susceptibility diverges as $\beta\uparrow\beta_{\rm c}$:
$$
\lim_{\beta\uparrow\beta_{\rm c}} \chi(\beta) = +\infty.
$$
