Solution (source code)

= Solution

Let $x=g\mu_BB/(k_BT)$. The single-ion <partition function> and mean projected <magnetic moment> are
$$
Z_J(x)=\sum_{m_J=-J}^{J}e^{xm_J},\qquad
\langle\mu_z\rangle=g\mu_B\frac{\sum m_J e^{xm_J}}{\sum e^{xm_J}}
=g\mu_B\frac{\partial\ln Z_J}{\partial x}.
$$
For $J>0$, the geometric sum also gives
$$
\langle\mu_z\rangle=g\mu_B\left[
\left(J+\tfrac12\right)\coth\!\left(\left(J+\tfrac12\right)x\right)
-\tfrac12\coth(x/2)\right].
$$
In the high-temperature, weak-polarization regime $|x|J\ll1$, expand the <Boltzmann factors>. The numerator's zero-order term vanishes because $\sum m_J=0$, and the first-order term uses $\sum m_J^2=J(J+1)(2J+1)/3$. The denominator is $2J+1+O(x^2)$. Therefore
$$
\boxed{\langle\mu_z\rangle\simeq\frac{(g\mu_B)^2J(J+1)B}{3k_BT}
=\frac CT,\qquad
C=\frac{(g\mu_B)^2J(J+1)B}{3k_B}.}
$$
This <fixed-field per-ion Curie coefficient> describes \b[mean moment per ion at fixed field], so it includes $B$. It is not the field-independent volume-susceptibility <Curie constant>. For <number density> $n$ and $B\simeq\mu_0H$, the latter is $C_\chi=n\mu_0(g\mu_B)^2J(J+1)/(3k_B)$, giving $\chi_H=C_\chi/T$. At $J=0$, the specified moment is zero and both expressions vanish.