Three versions of the classical $XY$ model and their low-temperature behaviors 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$$
Consider $i$ and $j$ a large distance $r$ apart. How does the spin-spin correlation function $\langle \cos(\theta_i - \theta_{j})\rangle$ for $i$ and $j$ a distance $r$ apart depend on $r$ at small temperature $T$?
As I explain below, I naively expect that in the three respective situations above, the correlator tends to a constant, decays algebraically as $\sim \frac{1}{r^{\eta}}$, and decays exponentially $\sim e^{-\frac{r}{\xi}}$. Here, I only write the leading functional behavior. However, I'm prepared to be surprised!

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, 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, and I'm open to learning new methods of deducing how correlation functions decay (or don't decay) at low temperature. I am also happy to accept references discussing related models, if the methods in those papers or textbooks are valid and can be applied to the models above, particularly $H_1$ and $H_4$.
 A: Edit: I've added some references (modified slightly the answer accordingly), some rigorous thanks to Y. Velenik, and some reporting numerical simulations.
I will not comment on $H_1$. It is non-analytic, but it has been shown that this does not change the fact that long-range order is impossible also in this case, see arxiv:math/0110127.
As said by the OP, $H_2$ is equivalent to the standard XY model up to trivial relabelling.
The physics of $H_4$ can be very interesting. It can be rewritten as
$$
H_4\propto-J\sum_{\langle ij\rangle} (\cos(\theta_i-\theta_j)-4\cos(2\theta_i-2\theta_j)).
$$
There is a nice competition between the two terms. We see that the second term does not care if neighboring spins differ by $\pi$. It corresponds to the so-called $\phi/2\phi$ model, which reads
$$
H_{\phi/2\phi}=-\sum_{\langle ij\rangle} (J_1\cos(\phi_i-\phi_j)+J_2\cos(2\phi_i-2\phi_j)).
$$
I will comment on that model in the case $J_1,J_2>0$, since I know it better, but the physics for J_2<0 will probably be similar up to anti-ferromagnetic order.
The most interesting regime is for $J_2\gg J_1$. At high temperatures, $T\gg J_2$, the system is disordered. When $T\simeq J_2$, the spins quasi-order up to a $\pi$-phase, i.e. there is some kind of nematic ordering (the spins are aligned, up to a free $\pi$ difference). This transition is BKT-like, but with a jump of the stiffness which is $4$ times the standard one.
Then, at lower temperatures, there is a 2D Ising transition, when the spins finally completely align.
Here are a few references: there is a rigorous result in agreement with discrete symmetry breaking in a similar model. See also this paper for numerical simulations, and this one for a thorough discussion of the phase diagram of this model and a slight generalization of it.
