Solution (source code)

= Solution

The spin-length sum rule and Fourier inversion give
$$
1=\langle|\sigma_0|^2\rangle
=\frac1{|\Lambda_L|}\sum_k\widehat G_L(k).
$$
The zero mode contains the square of the field-directed magnetization, while the <infrared bound> controls all nonzero modes. After $L\to\infty$ and then $h\downarrow0$,
$$
m(\beta)^2
\geq1-\frac{n}{2\beta}
\int_{[-\pi,\pi]^d}\frac{dk}{(2\pi)^d\,\varepsilon(k)}.
$$
Near zero, $\varepsilon(k)\asymp|k|^2$, so the integral is finite exactly when $d\geq3$. Choose $\beta_0$ larger than the resulting finite constant. Then for $\beta>\beta_0$ the right side is positive, and $m(\beta)\geq c>0$.

If $h\downarrow0$ before $L\to\infty$, every finite torus retains global $O(n)$ symmetry and its magnetization is zero. The reversed iterated limit is therefore zero. The order of limits is what permits <spontaneous magnetization>.