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.
![enter image description here](https://i.stack.imgur.com/Etqbu.png)
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.
Best Answer
There are many thorough references on this topic.
A "spin wave" / "magnon" shows up in the spin-spin correlator of the system, which is seen by neutron scattering. Specifically, it measures the correlator $$ \langle S^+(q, \omega)S^-(-q,0)\rangle $$ which in linear spin wave theory [in a nutshell: one equates spin raising (lowering) with the creation (annihilation) of a boson] is proportional to
$$ \langle a^\dagger(q, \omega)a(-q,0)\rangle $$
i.e. the propagator of a free Boson field (that can be calculated analytically). One may then dress these magnons by interactions in the usual perturbative way.
Triple-axis single crystal measurements can measure sharp, dispersive excitations in most conventional magnets. We are inclined to call these 'magnons' because they