Solution (source code)

= Solution

Let the coordination number be $z_c=4$, and choose the spin-space $z$ axis along $\mathbf S=S\hat{\mathbf z}$. The mean-field decoupling of a bond is
$$
\mathbf S_i\cdot\mathbf S_j\simeq
\mathbf S_i\cdot\mathbf S+\mathbf S\cdot\mathbf S_j-S^2.
$$
The original directed neighbor sum counts each bond twice, so its factor $1/2$ cancels that double counting. The resulting energy is
$$
\boxed{E_{\rm MF}=-z_cJS\sum_i S_{iz}+\frac{Nz_cJ}{2}S^2.}
$$
The constant prevents double counting of the mean interaction energy and is necessary if <free energy> is minimized with respect to $S$. With the source's unnormalized solid-angle measure, the single-site integral is
$$
Z_1(x)=2\pi\int_{-1}^{1}e^{x u}\,du
=4\pi\frac{\sinh x}{x},\qquad x=\beta z_cJS.
$$
Therefore
$$
\boxed{Z_{\rm MF}(S)=e^{-\beta Nz_cJS^2/2}
\left(4\pi\frac{\sinh(\beta z_cJS)}{\beta z_cJS}\right)^N,
\qquad z_c=4.}
$$
At $S=0$ take $\sinh x/x\to1$. A normalized solid-angle measure would remove the factor $(4\pi)^N$ without changing the <self-consistency equation>. These are three-component classical <spins> on a two-dimensional positional lattice, not <spins> constrained to a plane.

The PDF's hint is a correct angular-integration identity for a term with coefficient $\beta J/2$. That coefficient is not the full single-site molecular field after summing all four neighbors and decoupling both entries of each bond. The net coefficient is $\beta z_cJ$, as follows from the explicitly stated <Hamiltonian>; inserting the hint's coefficient as the total field would lose the coordination and double-counting factors.