Solution (source code)

= Solution
{c}

Let $s=\langle S_i^z\rangle$ be uniform, with dimensionless spin operators $S=1/2$. The factor $1/2$ in the source Hamiltonian removes double counting of pairs. Decoupling $\mathbf S_i\cdot\mathbf S_j$ about a uniform polarization gives the one-site <mean-field theory> Hamiltonian
$$
H_i^{\mathrm{MF}}=-sJ_0S_i^z+\mathrm{constant},\qquad J_0=\sum_jJ_{ij}.
$$
A body-centred-cubic site has eight nearest neighbors and six next-nearest neighbors, so
$$
\frac{J_0}{k_B}=8(1\,\mathrm{mK})+6(-0.5\,\mathrm{mK})=5\,\mathrm{mK}.
$$
For a spin one half in this effective field, the <mean-field self-consistency equation> is
$$
s=\frac12\tanh\left(\frac{J_0s}{2k_BT}\right).
$$
Linearizing at the continuous onset gives $s=J_0s/(4k_BT)+O(s^3)$, hence
$$
\boxed{T_C^{\mathrm{MF}}=\frac{J_0}{4k_B}=1.25\,\mathrm{mK}.}
$$
Equivalently, the general spin-$S$ result is $k_BT_C^{\mathrm{MF}}=S(S+1)J_0/3$. The antiferromagnetic next-nearest-neighbor exchange reduces the uniform ordering scale rather than contributing its absolute magnitude. The negative nuclear <g-factor> does not change this exchange instability; aligning the spins also aligns their moments with one another.

For these particular competing couplings, the uniform mode is indeed the largest exchange eigenvalue. Writing $x=\cos(q_xa/2)$, $y=\cos(q_ya/2)$ and $z=\cos(q_za/2)$ gives
$$
\frac{J(\mathbf q)}{k_B\,\mathrm{mK}}=8xyz-2(x^2+y^2+z^2)+3\le5,
$$
with equality at $\mathbf q=0$ modulo <reciprocal lattice vectors>. To see the bound, replace a negative $xyz$ by its absolute value for an upper bound, let $t=(|xyz|)^{1/3}$, and use $x^2+y^2+z^2\ge3t^2$. The expression is at most $8t^3-6t^2+3$, whose maximum on $0\le t\le1$ is five. Thus another modulation does not preempt the stated ferromagnetic onset at this level of <mean-field theory>. Fluctuations can change the actual transition temperature.