\frac{i}{2}(q_x+G_x)(q_y+G_y)}\,\,,
\label{Udis}$$ $${\rm U_H}({{\bf q}})=\frac{1}{2\pi}\sum_{{{\bf G}}} e^{-\frac{1}{2}|{{\bf q}}+{{\bf G}}|^2}
\frac{2\pi e^2}{|{{\bf q}}+{{\bf G}}|}(1-\delta_{{{\bf q}}+{{\bf G}},0})\,\,,$$ and $${\rm U_F}({{\bf q}})=-\frac{1}{A }\sum_{{{\bf q}}'}
e^{-\frac{1}{2}|{{\bf q}}'|^2}e^{i{q'}_x q_y-iq_x q'_y}
\frac{2\pi e^2}{|{{\bf q}}'|}(1-\delta_{{{\bf q}}',0})\,,
\label{U_F}$$ with ${{\bf G}}=(\frac{2\pi N_\phi}{L_x} n_x,\frac{2\pi N_\phi}{L_y} n_y)$.

For $v_E({{\bf r}})$, we choose a white noise potential without spatial correlation, $\langle\langle v_E({{\bf r}})v_E({{\bf r}}')\rangle\rangle=\sigma^2 \delta({{\bf r}}-{{\bf r}}')$. The density of states in the non-interacting limit has been calculated exactly by Wegner for this distribution of the disorder potential,[@Wegner] yielding a full width at half maximum of approximately $1.06 \sigma$. In our calculations the parameter $\gamma=e^2/\epsilon\sigma$ specifies the relative strength of interactions and disorder broadening. This type of disorder potential distribution is characterized by a single parameter $\sigma$, which we use as our unit of energy. As we discuss later, our results are insensitive to correlations in the disorder potential on length scales smaller than $l$, but would change in some respects for disorder potentials which are smooth on the magnetic length scale. The HF equations are solved by an iterative approach which can create difficulties which must be addressed. The HF equations generally have many solutions that correspond to different, usually metastable, extrema of the HF energy functional. The challenge is 