Solution (source code)

= Solution

For the uniform <mean-field approximation> field $\mathbf M=M\hat{\mathbf e}_z$, the inverse fermion propagator is diagonal in <momentum>, frequency and spin: $-i\omega_l+\xi_k-\sigma M$. The <Grassmann Gaussian integral> therefore gives
$$
Z_{\rm uniform}=\int d\mu(M)\,e^{-\beta NM^2/(2U)}\prod_{k,\sigma=\pm1}\left[1+e^{-\beta(\xi_k-\sigma M)}\right].
$$
Here $d\mu(M)$ denotes the residual uniform-field measure; integration over its orientation only affects a nonextensive prefactor. Hence
$$
\boxed{S_{\rm eff}(M)=\frac{\beta NM^2}{2U}-\sum_{k,\sigma=\pm1}\ln\left[1+e^{-\beta(\epsilon_k-\mu-\sigma M)}\right].}
$$
The exchange field separates the two spin bands by $2M$. Integrating out the fermions is exact for a specified auxiliary field; restricting that field to be constant is the <mean-field approximation>.