Consider the following three nearest-neighbor stat-mech Hamiltonians on the 2$d$ square lattice:
$$H_{1} = -J\sum_{\langle ij \rangle} \left|\sin\left(\frac{\theta_i - \theta_j}{2}\right)\right| $$
$$H_{2} = -J\sum_{\langle ij \rangle} \sin\left(\frac{\theta_i - \theta_j}{2}\right)^2$$
$$H_{4} = -J\sum_{\langle ij \rangle} \sin\left(\frac{\theta_i - \theta_j}{2}\right)^4$$
At low temperatures $T$ and large distances $r$ in these three models, how does the spin-spin correlation function $\cos(\theta_i - \theta_{j})$ for $i$ and $j$ a distance $r$ apart depend on $r$?

***
The above models are all ferromagnetic, and it's tempting that at low temperatures, $\theta_i$ and $\theta_j$ will be close together for $i$ and $j$ nearest neighbors. My motivation for this question comes from trying to understand whether Taylor expansions of the Hamiltonian, and resulting classical field theories, are indeed valid.

Note that my $H_2$ matches the regular classical $XY$ model of $H_{XY} = -J \sum_{\langle ij \rangle} \cos(\theta_i - \theta_j)$ up to a multiplicative factor of $2$ and a constant offset. $H_2$ is well-known for its [BKT transition,][1] where, at low temperatures, the model has algebraically decaying correlation functions. To see this behvaior, in lecture notes and textbooks (e.g. Altland and Simons chapter 8.6), one begins by considering the low temperature behavior of the model in the continuum limit. Here, one integrates over the functional integrand of $e^{-S[\theta]}$ with the action $S$ corresponding to 
$$S = K\int\mathrm{d}^2x\ |\vec{\nabla}\theta|^2$$
for an appropriate constant $K$ proportional to $\beta$ and $J$. At low temperatures, where one can ignore the periodicity of $\theta$, it is straightforward to see that spin-spin correlation functions decay algebraically with distance.

I feel a naive variant of the above action, motivated by Taylor-expanding the Hamiltonians above, would be 

$$S_1 = K_1 \left( \left|\frac{d \theta}{dx}\right| + \left|\frac{d \theta}{dy}\right| \right)$$
$$S_2 = K_2 \left( \left|\frac{d \theta}{dx}\right|^2 + \left|\frac{d \theta}{dy}\right|^2 \right) = K_2 |\vec{\nabla} \theta |^2$$
$$S_4 = K_4 \left( \left|\frac{d \theta}{dx}\right|^4 + \left|\frac{d \theta}{dy}\right|^4 \right)$$

where $K_i$ are constants proportional to $\beta$ and $J$.

I believe these three field theories predict rather different behavior for the correlation functions. $S_2$ corresponds to the $XY$ action, with an algebraic decay. $S_4$ has very weak penalties for spin-waves and vortices, so it's tempting the correlation function would in fact feature an exponential decay. $S_1$ appears to me to strongly penalize spin-waves and vortices, and so it's tempting that this might correspond to a symmetry-broken state with the leading asymptotic behavior being constant.

However, I am relatively unfamiliar with going from classical lattice Hamiltonians to field theories, and I'm also unfamiliar with extracting the relevant information from (non-quadratic!) field theories. It's possible I have flaws in my argument at multiple steps. Thus, I hope that in learning the answer to the behavior of the spin-spin correlation function $\cos(\theta_i - \theta_j)$ at far-apart $i$ and $j$, I'll better understand some of these steps.


  [1]: https://physics.stackexchange.com/questions/255909/what-is-the-kosterlitz-thouless-transition