Solution (source code)

= Solution

Differentiate the <homogeneous spin-field effective action>:
$$
\frac{d\mathcal S_{\rm eff}}{dM}=\beta\left[\frac{NM}{U}-\sum_{\mathbf k,\sigma}\sigma n_F(\epsilon_{\mathbf k}-\sigma M)\right].
$$
At a stationary field, $M=(U/N)\sum_{\mathbf k,\sigma}\sigma n_F(\epsilon_{\mathbf k}-\sigma M)$. Define the <spin-summed density of states> per lattice site by
$$
\nu(\epsilon)=\frac2N\sum_{\mathbf k}\delta(\epsilon-\epsilon_{\mathbf k}).
$$
Then the requested form is
$$
\boxed{M=\frac U2\sum_{\sigma=\pm1}\sigma\int d\epsilon\,\nu(\epsilon)n_F(\epsilon-\sigma M).}
$$
The factor one half fixes the spin convention. With a <single-spin density of states> $\nu_s=\nu/2$, the prefactor would instead be $U$. In terms of site occupations, the same equation is $M=U(n_\uparrow-n_\downarrow)=2U\langle S^z\rangle$. The <Fermi distribution> includes the <chemical potential> in its argument definition; it must not be subtracted a second time.