small) local magnetic field at the center of our simulation cells which points in the x-direction. It is this space-dependent field, required on purely technical grounds, which has motivated developing the formalism in a manner which permits non-constant Zeeman coupling strengths. To ensure that this artificially field and the finite temperature do not affect the final solution, we lower the magnitude of these terms until no change is seen in local charge and spin densities, or in HF quasiparticle energies.

The phase diagrams discussed in Secs. III and IV, were obtained by starting from the non-interacting case and incrementing the interaction strength $\gamma$, taking as the starting densities the self-consistent densities from the previous $\gamma$ value. There is, of course, some hysteresis involved in this process so we do a backwards sweep on $\gamma$ once we have reached the maximum interaction strength for a given run. If they differ, we use the smaller of the values obtained in upward and downward sweeps for the energy per particle. The energy per particle is obtained using the expression $$\begin{aligned}
&&\frac{E}{N}=\frac{1}{\nu}\sum_{{{\bf q}}\in {\rm BZ}}{\rm U}_{\rm dis}({{\bf q}}) {\Delta_0}^*({{\bf q}})
+\frac{\gamma}{2\nu}\sum_{{{\bf q}}\in {\rm BZ}}{\rm U}_{\rm H}({{\bf q}}) |\Delta _0({{\bf q}})|^2\\
&&+{\rm U}_{\rm F}({{\bf q}})
