\epsilon_{\rm mf}=
J+ n_{\rm mf}\int_{L \times L}  \frac{A}{r^\alpha} d{\bf r} =
J+ n_{\rm mf} T L_\alpha,$$ where we use the notation $$L_{\alpha}=c_\alpha \beta A L^{2-\alpha}.$$ The constant $c_\alpha$ is a geometrical factor of order 1, given by the integration of $1/r^\alpha$ on a unit square centered at the origin. In particular, $c_0 = 1$. On the other hand, we have $n_{\rm mf}=1/(e^{\beta\epsilon_{\rm mf}}+1)$ since the occupation numbers $n_p$ can only assume the values $0$ or $1$. By using Eq.  to calculate $n_{\rm mf}$, we find the self-consistent equation $$\label{meanfield}
n_{\rm mf} = \frac{1}{e^{\beta J+ n_{\rm mf} L_\alpha }+1},$$ with the following expansion at large $L_\alpha$ $$\label{leadingNmf}
n_{\rm mf}=\frac{1}{L_\alpha}\left[\ln L_\alpha - \ln\ln L_\alpha
- \beta J+\dots \right].$$ Higher order terms in the square brackets are small if $\ln L_{\alpha}\gg \beta J, |\ln \ln L_{\alpha}|$. For fixed temperature $T$ and interaction strength $A$, these conditions are always satisfied at sufficiently large $L$ since $L_\alpha \propto L^{2-\alpha}$.

We have confirmed the validity of the mean-field approximation by Monte Carlo simulations. By using the Metropolis algorithm [@Metropolis1953] to sample the probability distribution $\propto e^{-\beta/2 \sum{p,p'}U_{pp'}n_p n_{p'}}$, see Eq. (\[H0\]), the equilibrium number of excited plaquettes can be approximated with arbitrary accuracy. This can be 