Why do spin correlation functions in Ising Models decay exponentially below the critical temperature? I'm trying to form a better understanding of the 2D Ising Model, in particular the behaviour of the correlation functions between spins of distance $r$.
I've found a number of explanatory texts that seem to indicate that at both above and below the critical temperature $T_c$, the correlation function $C(r)$ decays exponentially at some correlation length $\xi$, and this helps us determine typical domain sizes. At $T_c$, the correlation length diverges to infinity. Here's an image from (http://math.arizona.edu/~tgk/541/chap1.pdf) to illustrate what I mean ($\beta$ of course is inverse temperature).



*

*$T > T_c$: This makes sense to me - adjacent spins are practically independent so domains are tiny, and the correlation length tends to zero as temperature increases.

*$T < T_c$: This doesn't make intuitive sense to me - I was under the impression that below $T_c$ the domains were large enough to ensure spontaneous magnetisation was observed. However the above is indicating that the typical correlation length ($\xi$) tends to zero as the temperature rises - and therefore domains shrink? I would have thought:


*

*The correlation function $C_r$ should not decay exponentially at all, but remain at a constant (as is indicated on Chapter 1, page 6 of Paul Fendley's "Modern Statistical Mechanics" draft)

*If they do decay exponentially, they would decay to a value > 0 (as is indicated on page 216 of Sethna's "Entropy, Order Parameters and Complexity")

*If they do decay exponentially and to 0, then surely as $T \rightarrow 0$, the correlation length $\xi$ should tend to infinity, as at $T = 0$ we know that all spins must perfectly correlate with arbitrarily far spins?


*$T = T_c$: If the correlation length is infinity, surely we would see perfect correlation and total magnetisation in the region close to the critical temperature? Instead at the critical temperature references indicate the correlation function takes the form $C_r \sim r^ {- \lambda} $ - the intuition of how this leap is made isn't clear to me.
Clearly, I have fundamentally misunderstood either how the correlation function and critical length work, how it relates to domain sizes, or both. I would really appreciate if someone could indicate where I have misunderstood.
For reference, I have an economics, not physics background, but need to understand the intuition of these models for my doctoral research in opinion dynamics. My primary reference text is "Dynamical Processes on Complex Systems" by Barrat et al.
 A: First note that, as you say, the 2-point function $\langle \sigma_i \sigma_j\rangle$ does not tend to zero as $\|j-i\|\to\infty$ when $T<T_{\rm c}$; namely,
$$
\lim_{\|j-i\|\to\infty} \langle \sigma_i \sigma_j\rangle^+ = (m^*(T))^2,
$$
where $m^*(T)=\langle\sigma_0\rangle^+$ is the spontaneous magnetization (the $+$ superscript indicates that the expectation is taken in the $+$ state).
So, when one says that correlations decay exponentially when $T<T_{\rm c}$, one is actually speaking of the truncated correlation functions. The truncated 2-point function is defined as
$$
\langle \sigma_i; \sigma_j\rangle^+ =
\bigl\langle (\sigma_i - \langle \sigma_i\rangle^+) (\sigma_j - \langle \sigma_j\rangle^+) \bigr\rangle^+ =
\langle \sigma_i \sigma_j\rangle^+ - \langle \sigma_i\rangle^+ \langle \sigma_j\rangle^+.
$$
It measures how correlated the fluctuations at $i$ and $j$ are. (In probabilistic terms, this is simply the covariance between the random variables $\sigma_i$ and $\sigma_j$.)
The truncated 2-point function converges to $0$ at all temperatures $T$. This is a general fact, true in the pure phases of any models, at any temperature.
Now, turning to your question, the truncated correlations indeed decay exponentially fast in the 2d Ising model for all $T\neq T_{\rm c}$. 
When $T<T_{\rm c}$, you should think about it as follows: at low temperatures, spins typically take the same values (say $+1$ in the $+$ phase), with only rare fluctuations. One useful way to look at these fluctuations is by taking a geometrical viewpoint on the configurations: draw a unit length segment separating each pair of nearest-neighbor vertices at which the spins take opposite values. The union of these segments form the Peierls contours of the configuration.

It is easy to check that the energetic cost associated to each such contour is proportional to its length. Note also, that the contours provide a complete description of the configuration (if you know that you are in the $+$ phase).
Now, what has this to do with the exponential decay of correlations? As mentioned above, the truncated 2-point function $\langle \sigma_i; \sigma_j\rangle^+$ measures how correlated the fluctuations at $i$ and $j$ are. What is the event that will lead to a simultaneous flip of both the spin at $i$ and $j$? It is not difficult to convince oneself that this should occur when a large contour surrounds simultaneously the vertices $i$ and $j$. Namely, using correlation inequalities, one can check that
$$
0\leq \langle \sigma_i; \sigma_j\rangle^+ \leq {\rm Prob} (\text{there exists a contour surrounding both $i$ and $j$}),
$$
where the $+$ superscript indicates that the computation is done in the $+$ phase. But the event in the right-hand side has a probability that is exponentially small in the distance $\|j-i\|$.
All this can be made rigorous. At sufficiently low temperatures, you can use cluster expansion techniques (this works in any dimension $d\geq 2$). You can find a detailed argument, for example, in Theorem 5.27 in this book. In dimension $2$, this can also be proved non-perturbatively, either through explicit computations, or by relating the cost of a large contour to the surface tension in the model; see, for example this paper.
As a final remark : the above argument suggests, and you can make this rigorous (see for example the cluster expansion proof I referred to above) that the correlation length tends to zero when $T\downarrow 0$. This is because it becomes extremely unlikely to have a contour surrounding two distant vertices and fluctuations at very low temperatures become essentially purely local events, thus occurring (approximately) independently at different vertices.
