I'll only answer the first part of your question.
Away from a critical point, the typical decay of the (truncated) 2-point function in lattice systems follows the Ornstein-Zernike asymptotics: as $\|y-x\|\to\infty$,
$$
\langle \phi(x);\phi(y) \rangle = \frac{C(n_{y-x})}{\|y-x\|^{(d-1)/2}} e^{-\nu(n_{y-x})\|y-x\|}\, (1+o(1)),
$$
where $\|\cdot\|$ denotes the Euclidean norm and $n_{y-x}=(y-x)/\|y-x\|$. The function $C$ and the mass $\nu$ are anisotropic, but analytic, strictly positive functions of the direction.
The above has been proved rigorously in many models in perturbative regimes, as well as non-perturbatively in simple systems (for instance, the Ising or Potts models in any dimension above the critical temperature, or the Ising model in any dimension in the presence of a nonzero magnetic field).
The behavior stated above is for "ferromagnetic" systems. For more general systems, there is usually an oscillatory term in the prefactor.
Note that Ornstein-Zernike behavior occurs very often, but not always. For instance, in the planar Ising model below the critical temperature, it is known that the power in the prefactor is $2$ rather than the usual $1/2$. This is a pathology related to planarity.
More importantly, the range of the interaction plays a crucial role. Namely, if you consider, for instance, a ferromagnetic Ising model with a formal Hamiltonian
$$
\mathcal{H} = -\sum_{x,y\in\mathbb{Z^d}} J_{y-x} \sigma_x\sigma_y,
$$
then:
- Ornstein-Zernike behavior applies when the coupling constants $J_z$ decay faster than any exponential of $\|z\|$ (that is, $\lim_{\|z\|\to\infty}\frac{1}{\|z\|}\log J_z = -\infty$).
- When $J_z$ decays slower than any exponential of $\|z\|$ (that is, $\lim_{\|z\|\to\infty}\frac{1}{\|z\|}\log J_z = 0$), the behavior is different, namely one can prove that, for all $\beta<\beta_c$, as $\|y-x\|\to\infty$,
$$
\langle \sigma_x\sigma_y \rangle = \beta\chi^2 J_{y-x}\, (1+o(1)),
$$
where $\chi$ is the magnetic susceptibility.
- When $J_z$ decay exponentially with $\|z\|$ (that is, $-\infty < \lim_{\|z\|\to\infty}\frac{1}{\|z\|}\log J_z < 0$), it is the behavior of the sub-exponential prefactor to the exponential decay of $J_z$ that will determine which of the asymptotic behaviors is relevant. There is a known "if and only if" criterion, but it is a bit too complicated to state here.
Finally, the above is for the spin-spin correlations. The behavior of other correlations can be more complicated. For instance, for the finite-range Ising model above the critical temperature, the "energy-energy" correlations behave like
$$
\langle \sigma_a\sigma_b;\sigma_{x+a}\sigma_{x+b} \rangle \sim \frac{1}{\psi(\|x\|)} e^{-2\nu(n_x)\|x\|},
$$
(for fixed $a,b\in\mathbb{Z^d}$) as $\|x\|\to\infty$, where
$$
\Psi(r) =
\begin{cases}
r^2 & \text{if }d=2,\\
(r\log r)^2 & \text{if }d=3\\
r^{d-1} & \text{if }d\geq 4.
\end{cases}
$$
This rather recent survey presents a panorama of some of the known (mathematically rigorous) results for the Ising model. You can also have a look at this older paper, which includes a brief review of the phenomenology for a larger class of models.
Some additions, answering your further questions in the comments.
[T]he authors define the mass (or inverse correlation length) by $m=\xi−1=\lim_{|x|\to\infty}\frac{1}{|x|}\log\langle \phi(x)\phi(y)\rangle$, which basically tells us that $\langle \phi(x)\phi(y)\rangle \sim e^{−|x−y|/\xi}$.
Notice that if $\langle \phi(x)\phi(y)\rangle \sim f(|x-y|)e^{-|x-y|/\xi}$ with $f$ sub-exponential (say, $f(r)=r^{-(d-1)/2}$, for instance), then
$$
\lim_{|x-y|\to\infty} \frac{1}{|x-y|}\log \bigl(f(|x-y|)e^{-|x-y|/\xi}\bigr) = - \frac{1}{\xi} + \lim_{|x-y|\to\infty} \frac{\log f(|x-y|)}{|x-y|} = -\frac{1}{\xi},
$$
so this is perfectly compatible with what the authors state.
In particular, taking the latter definition, if $\xi\to\infty$ at criticality, then $\langle \phi(x)\phi(y)\rangle \sim 1$, not $|x−y|^{−a}$. It seems inconsistent to the critical behavior.
Apart from the confusion already discussed just above, there is another important problem. The asymptotic behavior described in the first part of my answer (in particular the OZ asymptotics) are long-distance phenomena. In particular, they only apply when $|x-y| \gg \xi$. When approaching criticality, the correlation length $\xi$ diverges and (often) the relevant behavior is between spins located at vertices such that $|x-y| \ll \xi$. The exponents governing the latter are different from the exponents in the prefactor to the exponential decay when $|x-y|\gg\xi$. This is in particular true at the critical point, since $\xi=\infty$ in that case. (For instance, the planar Ising model for any $\beta<\beta_c$ has $\langle\sigma_0\sigma_x\rangle_\beta \sim |x|^{-1/2} e^{-\nu_\beta(n_x)|x|}$ (when $|x|$ is much larger than the correlation length), while $\langle\sigma_0\sigma_x\rangle_{\beta_c} \sim |x|^{-1/4}$.)
In order to understand better how these various regimes follow from a common asymptotic behavior, I strongly recommend to have a look at this paper. The latter provides detailed computations in the simpler case of the lattice Green function. In particular, they obtain a formula for the correlation valid in all regimes. One can then see how this general behaviors crosses to OZ behavior when the distance is large compared to the correlation length, and to a power law decay far below the correlation length (in particular, at the critical point).
