Multifractality and quantum-to-classical crossover in the Coulomb anomaly at the Mott–Anderson metal–insulator transition

We study the interaction-driven localization transition, which a recent experiment (Richardella et al 2010 Science 327 665) in Ga1−xMnxAs has shown to come along with the multifractal behavior of the local density of states (LDoS) and the intriguing persistence of critical correlations close to the Fermi level. We show that the bulk of these phenomena can be understood within a Hartree–Fock (HF) treatment of disordered, Coulomb-interacting spinless fermions. A scaling analysis of the LDoS correlation demonstrates multifractality with the correlation dimension d2 ≈ 1.57, which is significantly larger than at a non-interacting Anderson transition and is compatible with the experimental value dexp2 = 1.8 ± 0.3. At the interaction-driven transition, the states at the Fermi level become critical, while the bulk of the spectrum remains delocalized up to substantially stronger interactions. The mobility edge stays close to the Fermi energy in a wide range of disorder strength, as the interaction strength is further increased. The localization transition is concomitant with the quantum-to-classical crossover in the shape of the pseudo-gap in the tunneling density of states, and with the proliferation of metastable HF solutions that suggest the onset of a glassy regime with poor screening properties.


Introduction
The recent tunneling experiments on Ga 1−x Mn x As by Richardella et al [1] have demonstrated critical, multifractal correlations in the local tunneling density of states (DoS) at the metal-insulator (MI) transition in this semiconductor. While this would be expected from the theory of Anderson localization in non-interacting fermions, the experiment bears clear signs of the relevance of electron-electron interactions. Most strikingly, the critical correlations typical of the mobility edge ε ≈ ε m were observed very close to the Fermi level ε ≈ ε F , even upon doping further into the metallic regime. This phenomenon clearly originates in interactions, which single out the Fermi level as a special energy. In a non-interacting system, ε m , which is controlled by disorder, and ε F , which is controlled by the electron density, are generally very different.
The effect of interaction has been studied on either side of the MI transition [2], but rather little is known about its role close to criticality. In a seminal work, Efros and Shklovskii [3,4] showed that deep in the insulator the Coulomb interactions between localized electrons create a pseudogap in the DoS near the Fermi level ε F , where the DoS vanishes as ρ(ε) ∼ |ε − ε F | 2 in three-dimensional (3D). This is reflected in the Efros-Shklovskii law of variable-range hopping conductivity at low temperatures, ln σ (T ) ∝ − √ T 0 /T . In the opposite limit of weakly disordered metals, Altshuler and Aronov [5] discovered interaction corrections to both the DoS near ε F and the low T conductance. In particular, for spinless fermions, disorder and repulsive interactions both enhance the tendency to localize. In the weakly localized regime, the tunneling DoS has a dip at ε F with [5] √ ω. However, in contrast to the classical Efros-Shklovskii pseudo-gap, the Altshuler-Aronov corrections are of purely quantum (exchange) origin: δρ ∝h 3/2 , if the diffusion coefficient and the Fermi-velocity are held fixed. The quantum corrections of Altshuer et al [5] can be effectively summed to obtain a nonperturbative result near d = 2 by using the formalism of the nonlinear sigma-model due to Finkel'stein [7] and Belitz and Kirkpatrick [8]. An effective action approach was suggested by Levitov and Shytov [9] and Kamenev and Andreev [10] to derive the non-perturbative expression for the tunneling DoS in a weakly disordered two-dimensional system near the Fermi energy. Remarkably, in the lower critical dimension d = 2, the DoS at ε F vanishes exactly in the thermodynamic limit. In higher dimensions with d > 2, the above results suggest the following qualitative picture [8,12,14]: the pseudo-gap in the one-particle DoS gradually grows with increasing disorder or repulsion strength. ρ(ε F ) eventually vanishes at the localization transition, and remains zero in the insulator. The shape of the pseudo-gap evolves from the quantum behavior δρ ∝ √ ω in the metal to some non-trivial power ρ(ε F + ω) ∝ ω µ , µ > 0 at the Anderson transition point, to the classical ω 2 behavior in the deep insulator. A power law suppressed DoS at criticality, ρ(ω) ∝ ω µ is also predicted within an -expansion in d = 2 + dimensions [8]. However, the actual scenario might be more complex if the localization transition is accompanied, or even preceded, by a transition to glassy or other density-modulated phases, aspects which so far have not been taken into account by existing theories. A fascinating property of electronic eigenfunctions near the localization transition of noninteracting particles is their multifractality [15]. This is an exact property of critical states at the mobility edge ε = ε m , but, off-critical states also exhibit a multifractal character inside their localization or correlation radius ξ [16]. This fractality has important effects on interactions (both repulsive and attractive), as it enhances their local matrix elements. In the case of predominantly attractive interactions, it may induce local pairing gaps in weak Anderson insulators. In more-conducting systems, it may lead to enhanced superconducting transition temperatures [17][18][19][20].
It remains an unresolved theoretical question whether such subtle wavefunction correlations survive in the presence of Coulomb interactions. The reason to doubt their survival is most easily seen on the level of a Hartree-Fock (HF) approach, where the best one-particle HF orbitals in the presence of interaction are a linear combination of non-interacting wavefunctions of different energies. If one naively assumes that the fractal patterns of such wavefunctions are only weakly correlated, one may expect partial or even complete degradation of the fractal structure in the HF wavefunctions, due to a superposition of a large number of random uncorrelated fractal patterns. In reality the fractal patterns of non-interacting wavefunctions are strongly correlated even at a large energy separation [17,18,20]. Nevertheless, the question remains whether such interaction-induced superpositions give rise to a change of the fractal dimension of HF wavefunctions or may destroy fractality completely, despite the correlations. There is indeed a subtle interplay between the strength of correlations and the effective number of non-interacting wavefunctions which superpose in the HF wavefunctions. Another naïve argument, advocating the opposite conclusion, puts forward that the HF Hamiltonian is essentially a one-particle Hamiltonian of the same basic symmetry as in the non-interacting case. Hence, invoking universality, one would expect the same statistics for both non-interacting and HF wavefunctions. The flaw in this argument is that the matrix elements of the HF Hamiltonian, which are self-consistently determined, possess correlations which could be long-range in the presence of badly screened long-range interactions between particles. Thus the two different naïve arguments lead to two opposite conclusions. According to one of them, the fractal pattern of the HF wavefunction should be smeared out while the other one advocates unchanged fractal patterns. One of the main results of this study is to show that the first argument is in fact closer to reality.
On the other hand, the direct observation of multifractality in the tunneling spectra of Richardella [1] strongly suggests that the multifractality survives interaction. Indeed, the measured auto-correlation function of the local DoS (LDoS) showed a well-established powerlaw decay with distance on the sample surface. Surprisingly, this critical behavior appeared to be nearly pinned to the Fermi energy without any fine-tuning of the Mn impurity concentration, implying that the mobility edge ε m remains close to ε F in a broad range of disorder strengths. As mentioned above, this indicates the importance of interactions, since they single out ε F as the center of the pseudogap. In this letter, we address the problem of multifractality at the interacting localization transition theoretically, and study the mechanism by which interactions pin the mobility edge ε m nearly to ε F in a broad parameter range.
As a naive rationale for this pinning of the mobility edge one may consider that localization occurs earlier where the DoS is lower. Thus localization is naturally prone to occur first within the interaction-induced pseudogap, making ε m track ε F rather closely. However, it is only the single-particle (tunneling) DoS that has a pronounced dip near ε F , while the global thermodynamic DoS, dn/dµ, that enters the conductivity via the Einstein relation, usually shows a different behavior. Hence, it is not obvious which notion of DoS is relevant for localization and transport purposes (cf the discussions in [11,12,21] about global versus local (dn/dµ) −1 ). A more-detailed analysis is thus required in order to show that indeed the LDoS close to ε F can be critical, while in the bulk of the spectrum the correlations are still metallic.

Model
To address these questions we consider a model of spinless fermions on a 3D cubic lattice of size N = 10 3 with the tight-binding Hamiltonian interacting via long-range Coulomb repulsion We employ periodic boundary conditions and choose t = 1 as the unit of energy. The onsite energies i are random, independently and uniformly distributed in i ∈ [− W 2 , W 2 ]. The chemical potential µ depends on interaction and is chosen so as to keep the average density 1/2. For non-interacting particles (U = 0), the localization transition is known to occur at the disorder strength W c = 16.5 in this model [22]. In the present work we choose W = 14 < W c (not particularly close to W c ), so as to mimic conditions of Richardella [1] where the impurity concentration was not specially tuned.
We attack this problem numerically by considering the interactions in the HF approximation. This amounts to studying an effective single-particle model with self-consistent on-site energies and hopping amplitudes. In order to clarify the role of long-range interactions, we truncated the Coulomb interaction at a finite range, and then progressively increased its range up to the size L = 10 of the 3D system, defining r i j as the shortest distance on the torus. We first took into account the Hartree terms (occupation numbers n j in the sum U j r −1 i j n j ) and the Fock terms (expectation values of c † i c j ) up to the fifth nearest neighbors. Then we considered the Fock terms up to the fifth nearest neighbors while the Hartree terms were considered up to the 20th nearest neighbors. Finally we tackled the full self-consistent problem for all neighbors.

Overview of results
The main result of our paper is to establish the persistence of multifractality in the presence of a full-range Coulomb interaction. Notably, the fractal dimension we find, d 2 ≈ 1.57 ± 0.05, appears to be significantly larger than in the non-interacting case. With a decreasing range of interaction, the effective d 2 in a finite sample decreases (d 2 = 1.38 ± 0.05 for interaction up to the fifth nearest neighbor) until it reaches its value d 2 = 1.35 ± 0.05 for the non-interacting case. This marks essential progress in comparison to earlier works based on the HF approach [14]. As we will describe below, the critical behavior exhibits various further interesting features that are specific to the interacting case. Most importantly, we establish that within the insulating phase, even considerably far from the MI transition, the mobility edge remains very close to the Fermi level.
Furthermore, we study the evolution of the pseudo-gap in the HF DoS ρ(ω) as the increasing interaction drives the system towards the localization transition. In particular we confirm (within our accuracy) the scaling relationships suggested by McMillan [23], Shklovskii and co-workers [12] which relate the critical power law of the pseudo-gap ρ ∝ ω µ with the dynamical scaling exponent η and the exponent that describes the dependence of the static dielectric constant κ 0 ∝ ξ η−1 on the localization radius ξ in the insulator phase.
Finally, for the first time we address the question of multiplicity of HF solutions, and the competition of the related glassy features and localization. We show that within our accuracy the onset of a multiplicity of solutions with increasing interaction strength (an indication of an emerging glassy energy landscape) coincides with the localization transition. In contrast, the charge ordering (typical for a Mott transition) occurs at much stronger interaction. Thus we argue that the localization transition should better be called an Anderson-glass transition, rather than an Anderson-Mott transition.
A preliminary version of this paper was published as a preprint [13].

Hartree-Fock calculations
The effective HF Hamiltonian, which corresponds to the model given in equations (1) and (2), is HereṼ where t i j = t is the bare nearest-neighbor hopping and . . . 0 denotes the quantum-mechanical ground-state expectation value evaluated on the Slater determinant formed by the lowest N /2 HF levels. The effective on-site energyṼ j contains the interaction-induced Hartree term, which leads to correlated on-site energies (potentially at long range), while the effective hoppingt i j contains the Fock term, which may be long-range as well. We carried out calculations on a cubic 3D lattice of size L = 10, using three ranges of interactions. The first one took into account up to the 20th nearest neighbors (460 sites j nearest to i) in the Hartree term, and Fock terms corresponding to the fifth nearest neighbors (the 56 nearest sites up to distance √ 5). A second calculation restricted the Hartree and the Fock terms equally to the fifth nearest neighbors in order to check the importance of the Hartree terms and to ensure a correct implementation of the Pauli principle. A third calculation performed the self-consistent HF calculations with the full range of Coulomb interactions.
Note that the role of the HF terms is not the same in the metal and the insulator. Deep in the insulator side, it is important to keep the longer range Hartree terms to obtain the full classical Efros-Shklovskii gap, while long range Fock terms are negligible due to a strong localization of the wavefunctions. In the metal, however, the role of the Fock terms is expected to be more significant, while the Hartree terms incorporate an effective screening at long distances.

Numerical implementation
Even though completely standard, we briefly review the main steps involved in finding solutions of the HF equations. The set X of parameters to be found self-consistently comprises all the c † i c j 0 , (∼ L 6 /2 parameters for the full-scale Coulomb interaction) plus the L 3 diagonal parameters n i 0 = c † i c i 0 (i.e. ∼ 500 000 parameters for L = 10). The chemical potential µ is always adjusted to assure half filling, as described below.
In order to find a self-consistent solution we begin with a random initial guess for all the parameters X (0) in satisfying the condition where the number of particles N e = N /2 is fixed in our calculation (half-filling). Diagonalizing the effective Hamiltonian equations (3)-(5) using the initial X (0) in one obtains the eigenfunctions ψ m (r) and eigenvalues ε m , from which we compute the output parameters X (0) where f (ε) is the Fermi distribution function with the chemical potential µ. It is to be found from the condition At T = 0, considered in this paper, there is an uncertainty of the position of µ between the two energy levels ε N e +1 and ε N e . For most of the calculations we have chosen µ = 1 2 (ε N e +1 + ε N e ). However, to avoid an artificial hard minigap at the bottom of the DoS dip, to study the latter we used a parametric mixing µ = (1 − a) ε N e +1 + a ε N e with 0 < a < 1; a was fixed for a given disorder realization, but taken at random for different disorder realizations.
An updated set of parameters to be used as initial parameters X (n+1) in for the next, i.e. (n + 1)th, iteration is chosen as follows (n = 0, 1, . . .): The parameter α ∈ [0, 1] is chosen such that the iteration process is stable and leads to a convergent solution. The iteration procedure is terminated and the output set of parameters is taken as the converged solution if the absolute value of the difference between the values of all the parameters X of the previous and the final iteration is less than 10 −4 for the truncated Coulomb interaction, and 10 −5 for the full-range Coulomb interaction. Once a converged solution is obtained, and the final set of HF eigenfunctions ψ m (r) and eigenvalues ε m has been calculated, one can compute any quantity expressible in terms of ψ m (r) and ε m . The procedure is then repeated for different realizations of disorder to obtain disorder averaged quantities, such as the DoS and the LDoS correlation functions. We point out that the solution of the HF equations is unique only at small enough U within the metallic regime. At large U , one expects a number of solutions that grows exponentially with the volume. In this regime, we analyze typical solutions of the HF equations, without optimizing the HF energy among different solutions. This choice will be discussed and justified in section 8.
The typical number of iterations needed to obtain an HF solution for one realization of disorder was ∼2000 for the full-scale Coulomb interaction. The total computational time to obtain one HF solution was mostly limited by the time ∼ L 9 of diagonalization of a matrix Hamiltonian of the size L 3 × L 3 needed in each iteration. With a typical number of disorder realizations ∼2000 the total time at L = 10 was of the order of 1000/(#cores) hours for each parameter set of interaction and disorder strengths. For all values (∼20) of interaction strengths U necessary for our scaling analysis and the average number of cores ∼50 used for parallel computing, the total time was of the order of 400 h for L = 10.

McMillan-Shklovskii scaling
The MI transition in a disordered system is expected to occur as a second-order phase transition at some interaction strength U c . Close to criticality, where τ ≡ |1 − U/U c | 1, one expects a scaling form for the DoS as In the critical regime, whereas in the metal and the insulator which capture the shape of the Altshuler-Aronov and Efros-Shklovskii pseudogaps as limiting cases. The exponents in equation (12) are chosen such that the dependence on the critical parameter τ disappears at criticality. The scaling [12,23] is based on an assumption about the potential of a point charge within the critical regime, i.e. at a distance a r ξ . Here ξ is the correlation length, which diverges at the transition, as and a is a certain microscopic length (e.g. the distance between donors in doped semiconductors [12]). In our simulations it can be taken to be equal to the lattice spacing. The assumption is that the potential behaves as a modified power law, As equation (17) is essentially the relationship between the length scale r and the energy (or inverse-time) scale V , the exponent η should coincide with the dynamical scaling exponent z.
For the non-interacting case the dynamic exponent takes its maximum value η = d = 3, while the minimal theoretically admissible value cannot be smaller than the exponent of the Coulomb potential [23], η 1.
The exponent η also governs the scaling of the static dielectric constant in the insulator [23] To show this, it is enough to assume that in the insulator at distances r ξ the potential V (r ) takes the usual form of dielectric screening and matches with the potential in equation (17) at distances r ∼ ξ . The characteristic energy scale δ in equation (11) is set by the potential V (r ) at r ∼ ξ : In a system of finite size L, in the critical region where ξ L, one should replace δ by the mean level spacing δ L ∼ V (L).
Finally, the characteristic scale of the DoS can be expressed through δ and ξ by a relationship following from dimensional arguments: Equations (20) and (21), as well as the scaling equation (11), are valid for From equations (20), (21) and (12) one immediately obtains the following scaling relations [12,23]: in terms of η and the correlation length exponent ν. Note that the above scaling assumes only one critical scale ξ separating different regimes of V (r ). Should an additional scale (e.g. one related to a 'screening transition') appear, the exponent µ will be independent of the dynamical scaling exponent η.
In the next section we analyze the evolution of the DoS across the transition in the light of the above scaling assumptions.

Pseudo-gap in the density of states (DoS)
In figure 1 we present the DoS of the HF levels. One can see that deep in the metallic and in the insulating regimes, HF correctly captures the Altshuler-Aronov and Efros-Shklovskii pseudogap features discussed above, while it provides a non-trivial mean field approach to describing various interesting phenomena happening at and close to the MI transition. The curvature of ρ(ε) at small ω = ε − ε F is seen to change sign as U increases.
In figure 2 we verified the scalings equations (11), (13) and (14) by collapsing the 'lowenergy' data (in regime (22)) for ρ(ε) close to ε F onto the universal scaling functions f ρ,M (x) and f ρ,I (x). From the power-law behavior of f ρ,M at x 1, we found for the exponent µ which is very close to the value experimentally observed in the tunneling DoS close to criticality [24]. We cross-checked this result in figure 3 by plotting ln versus ln (1/δ). We obtained an almost linear curve in accordance with equation (12), with the same slope µ M = 0.53 as in figure 2(a). In the insulator, the best collapse corresponds to but the reliability of this exponent is not as high as the one on the metallic side (for instance a test similar to that in figure 3 yields a slope of 0.54, consistent rather with µ M , but smaller than µ I = 0.68). However, to be conservative we may conclude that the critical exponent µ lies in the range From the obtained exponent µ we can estimate the dynamical scaling exponent η This yields the exponent ζ = 0.9 ± 0.2, characterizing the divergence of the dielectric constant in equation (18), in a reasonably good agreement with the experimental value [12] ζ ≈ 0.71. 5 Thus we conclude that our results are compatible with the McMillan-Shklovskii scaling as well as with the available experimentally obtained values of the exponents µ and ζ .

Auto-correlation of the local DoS and the fractal dimension d 2
To study the multifractality of the LDoS, we have computed the spatial correlations of the HF wavefunctions |ψ n (r)| 2 : K (R; ε) = n,ε n ∈ (ε) r |ψ n (r)| 2 |ψ n (r + R)| 2 n,ε n ∈ (ε) where ε n are the associated eigenvalues, . . . denotes the ensemble average over random realizations of on-site energies n , and (ε) is a narrow interval of energies of the order of the mean level spacing δ, centered at ε. Multifractal correlations [15,18] imply that in the range of distances 0 < R < ξ , the correlation function is the same as at the localization transition, Here ξ is the localization or correlation length which diverges at the transition, and 0 is of the order of the lattice constant. d = 3 is the dimensionality of space and d 2 < d is the correlation fractal dimension. For R > ξ , the correlation function K distinguishes delocalized and localized regimes, saturating to a constant in a metal and decreasing exponentially in an insulator. Close to criticality one expects scaling behavior-the correlations should collapse to a single curve upon rescaling R by α (a finite size corrected version of ξ ), and amplitudes by β, and expressing K (R; ε) = β −1 f M,I (R/α), where f M,I (x) are universal scaling functions on the metallic and insulating sides of the transition, respectively. In order to optimize the choice of α, β (which both depend on U and ε, while W = 14 is fixed) and to determine the correlation dimension d 2 , we use a simple analytical ansatz for f M,I , which captures the multifractal characteristics of the eigenfunction correlations The fits were optimized by the value B ≈ 2. We first discuss the DoS correlations K (R; ε F ) at the Fermi level, upon varying the strength of the interaction U , cf figure 4. The good quality of the collapse demonstrates that the behavior of K (R; ε) is consistent (using full-range Coulomb interactions) with multifractal correlations with dimension d 2 ≈ 1.57 ± 0.05, full-range Coulomb.
This is significantly larger than the fractal dimension found for the non-interacting case in the limit of large sample sizes d 2 = 1.29 ± 0.05 [15] from the multifractal analysis of the moments  This result is in full agreement with the qualitative picture outlined in the introduction. It is also in line with recent results obtained via the = d − 2 expansion in the unitary ensemble [6]. According to that study Although the expansion fails to give an accurate prediction for d = 3, the tendency for 1 is clearly that d inter An increase of d 2 is also expected from studies of systems with frustrating interactions on the Bethe lattice, where the Efros-Shklovskii-(or Hartree-) type suppression of the DoS around the chemical potential is found to reduce the abundance of resonances (i.e. small denominators in the locator expansion). Therefore, the wavefunctions have less tendency to follow rare paths to increase the number of resonant sites visited, and thus form less sparse fractals [31].
The same analysis was repeated at the higher interaction strengths U = 1.5 and 3.0, where the critical HF states appear away from the Fermi energy. The result is that away from the Fermi energy, d 2 is practically indistinguishable from the non-interacting case.
To assess the effect of the range of the interactions, we also computed d 2 at the Fermi level when the Coulomb interaction was truncated as described in section 2 (the critical interaction strength in this case was U c ≈ 0.75, a bit smaller than for the full-scale Coulomb interaction). With the Coulomb interaction restricted to five nearest neighbors in the Fock terms, the correlation dimension was d 2 ≈ 1.39 ± 0.05, both when the Coulomb interactions in the Hartree terms were restricted to 5 or to 20 nearest neighbors. This result shows that Coulomb interactions of full range are essential to change the fractal dimension d 2 significantly. With a truncated Coulomb interaction, the effective d 2 in a finite sample gradually decreases and approaches its non-interacting value. Figure 5 shows the evolution of the finite size-corrected correlation length α(U ), as obtained from the scaling collapse of figure 4. It exhibits strong non-monotonicity, indicating a localization transition: for U < U < ≈ 0.79, α(U ) increases with increasing U while for U > U > ≈ 0.89, it decreases. The best fits to critical power laws yield ξ(U ) = a |U − U <(>) | −ν M(I ) , with ν M ≈ 0.50 ± 0.05, ν I ≈ 0.96 ± 0.05, full-range Coulomb.

Metal-insulator transition
(33) The difference in the fit exponents is too big to be a mere result of statistical errors, or of systematic errors related to the fitting procedure. Indeed, as an independent check we computed the critical exponents for a non-interacting system of equal size, using the same method. We obtained 6 much closer exponents ν M ≈ 1.20 ± 0.05 and ν I ≈ 1.08 ± 0.08. We thus believe that the difference in the fit exponents (33) is a genuine interaction effect, which persists to fairly large scales. Possible interpretations of these findings are discussed further below. In this context it is interesting to note that the exponent ν M ≈ 0.5 has been reported in earlier experiments on Si:P, which remained a puzzle for theorists for a long time [8]. Equations (33) might reflect the degradation of the multifractal pattern due to the interaction-induced mixing of non-interacting wavefunctions, which we expect to be much stronger in the delocalized than in the localized phase (where fewer non-interacting wavefunctions involved have a significant overlap). Another phenomenon that undoubtedly influences the MI transition is the gradual breakdown of screening in the metallic phase. The interactions, which are well screened deep in the metal, must become long range somewhere on the way to the insulator [12]. This entails a crossover, or even a phase transition, to a glassy phase [25,26]. We observe a trace of the latter via the onset of non-uniqueness of the HF solutions roughly at the same point as the MI transition, but our resolution is not sufficient to determine whether the two phenomena coincide. It also remains an interesting open question whether screening breaks down at the MI transition only, or already within the metal, as it happens in mean field models with similar ingredients [27].
It is interesting to compare our results for the exponent ν I with the -expansion [6] obtained from Finkel'stein's theory [7] in the unitary ensemble. According to this theory, These expressions are meaningless at = 1, where they evaluate to negative ν. One may think, however, that for small they give a correct relationship between ν inter and ν non−inter , as was the case for the fractal dimension d 2 in equation (32). However, the relationship between ν inter and ν non−inter is ambiguous in the region of < 1 where both of them are still positive. Indeed, one can see that for very small , one has ν inter > ν non−inter . However, as increases, ν non−inter catches up with ν inter , and at > 0.55 we have ν inter < ν non−inter , as in our results for d = 3. This may indicate that two competing mechanisms are at play, whose relative importance depends on the dimensionality.
We also note that the exponent ν I increases when the Coulomb interaction is truncated or when the mobility edge moves away from the Fermi level In contrast, the exponent ν M is almost insensitive to truncation, but decreases with increasing interaction strength The large error bars in the interval U ∈ [0.8, 0.9] do not allow us to determine ξ(U ) by an accurate treatment of the finite-size scaling α(U ) = ξ(U ) f α (ξ/L). 7 Two different scenarios may be envisioned to reconcile the fit exponents (33) with standard theoretical considerations: (a) there is a single localization transition close to U c = 0.9 with a shoulder in the dependence of ξ(U ) on the metallic side, due to an additional phase transition or a crossover in a different sector, such as the breakdown of screening or the onset of glassiness. In that case, ν M would be expected to approach ν I sufficiently close to U c and on large scales. (b) There are two separate transitions at U c1 ≈ 0.79 and U c2 ≈ 0.89, with critical wavefunctions in an entire finite interval U ∈ [U c1 , U c2 ]. In this more exotic (and less probable) scenario, there would be no a priori reason for the two exponents to coincide.

Finite-energy mobility edge in the insulator
A similar scaling analysis as above may be performed for ε away from ε F , which determines a critical line-the mobility edge ε m (U ). Of course, such a mobility edge is defined sharply only at the mean field level of the HF equations, which neglect the finite lifetime of higher energy excitations due to inelastic processes involving either phonons or delocalized excitations of purely electronic origin [29]. Nevertheless, the phase space for decay processes at low energies is strongly suppressed, and due to the pseudogap even more severely so than in an ordinary Fermi liquid. This gives us confidence that the features of single-particle HF levels are representative of the fully interacting system. In particular, the statement that ε m remains close to ε F is a result which we believe to be robust beyond the HF approximation. Hereby ε m should be interpreted as the (approximate) location where the LDoS correlations become critical up to the relevant length scale set by inelastic decay processes. Our result is shown in figure 6, together with the band edge, defined as the energy where ρ(ε) drops to half of its maximal value. Figure 6 demonstrates that the mobility edge is indeed trapped in a narrow range around ε F . This holds for values of U nearly all the way up to U * ≈ 4, where the last states around the maximum of the DoS localize (at W = 14). This confirms the expectations of Richardella et al [1] that in a relatively broad region of the parameters W and U , states near ε F are almost critical. Note also that U * is almost five times larger than U c ≈ 0.85, where the MI transition occurs. In fact U > U * brings the system already very close to a Motttype transition where a charge-density wave order sets in. It remains an interesting question to study how such charge ordering effects and glassiness (i.e. the multiplicity of HF solutions) affect the localization as the interaction strength increases. Figure 6. Phase diagram. In a wide range of the interaction strength U , the mobility edge (solid blue line) stays close to ε F as compared to the band edges (green dashed line). This behavior is in qualitative agreement with a conjecture [28] that the mobility

Multiplicity of Hartree-Fock solutions: glassiness and charge ordering
Finally we briefly address the issue of multiple solutions of the HF equations and its relation to the onset of glassy behavior as the interaction U increases.
Our iterative procedure to solve the HF equations begins with an input configuration of occupation numbers n in (r) on each of the N lattice sites. At sufficiently large U , the converged output HF solution n out (r) is generally different for runs with different inputs. To quantify this difference statistically we studied the quantity where the superscript m labels the set of N sol different solutions which were obtained from initial density patterns n (m) in (r), while 0 denotes a reference solution. In the simplest test we have chosen N sol = 10 solutions, out of which eight were obtained from random inputs and two had a checkerboard order as the input. The results for D(U ) are presented in figure 7. One can see that for U < 4, the average deviation of the solutions from the reference solution is small. It is thus reasonable to assume that physical properties evaluated on the various solutions are statistically very similar. However, starting from U ≈ 4, the function D (U ) sharply increases.
In order to check whether this increase is due to a significant variation between random HF solutions or whether this increase in D (U ) is due to the stabilization of a checkerboard density pattern in the HF solution, we consider the solution obtained from initializing with a checkerboard input and plot the difference between the solution and the corresponding input pattern.
The result for U = 4 shows (see figure 8(a)) that the difference has a clear checkerboard structure which implies that the output was random. However, as U increases to U = 5, the difference reduces significantly (see figure 8(b)), which signals the tendency to retain the checkerboard order in the solution. Comparing also with free energies of random solutions, we concluded that the transition to a checkerboard structure (charge density wave) occurs somewhere in the range 4 < U < 5.
A more-precise identification of the onset of the multiplicity of solutions shows that it starts at much smaller values U ≈ 0.7, which roughly coincides with the U c at which localization at the Fermi energy occurs 8 . In order to show this we generated ten different HF solutions at U = 1.0, and characterized them globally by the total energy E per site. The fact that the iterative HF procedure at U = 1.0 converges to different values of E is a manifestation of the existence of multiple local minima. Physically, one may expect that this will reflect in the onset of glassy behavior associated with the slow dynamics or relaxation between the minima that correspond to the various HF solutions. We then used the above solutions as inputs at the slightly decreased interaction strength U = 0.9, the resulting solutions still being of different total energy (see figure 9(a)). Upon decreasing the interaction in steps of 0.1, and using the solutions of the previous step as an initial condition, we found that at U ≈ 0.7 the total energies, after a large number of iterations, coincided ( figure 9(b)). That value of U can thus be interpreted, at this HF mean field level, as the border of a glassy regime. Upon further decrease of U the solutions did not diverge anymore, implying the existence of a unique HF solution (figure 9(c)).
The fact that in the insulating regime the HF equations develop a number of solutions that grows exponentially with the volume is to be expected, as this is well-known to be the case in the classical limit of vanishing hopping, t = 0. One may wonder whether and how key features like the Efros-Shklovskii Coulomb gap are present in typical solutions to the HF equations, or whether they occur only in the lowest energy solutions, which are very difficult to find. We argue here that all typical solutions are expected to exhibit a parabolic Coulomb gap, as we indeed observed numerically.
To understand how such a Coulomb gap comes about, consider the HF equations in the limit of vanishing hopping, t = 0, where all HF orbitals are completely localized. An HF solution consists in an assignment of occupation numbers n i ∈ {0, 1} to the sites i, according to whether the local potential is above (→ n i = 0) or below (→ n i = 1) the chemical potential µ (µ ≈ 0 is always adjusted to assure half filling). In the classical limit, the HF procedure consists of updating the occupation numbers until convergence to a stable point is reached. The final HF solution is a minimum of the HF energy with respect to the change of any of the n i , if the HF energy is written in grandcanonical form, including the term −µ i n i . That is, there is stability with respect to single particle addition or removal.  Figure 9. Convergence of the HF procedure for the total energy per site: (a) U = 0.9, the total energy converges to different output values, depending on the initial configuration of on-site occupation numbers; (b) U = 0.7, different initial configurations of occupation numbers obtained in the previous step converge to the same value of the total energy, and thus to the same HF solution; and (c) U = 0.5, no further divergence in the total energy is observed.
As long as there is no suppression in the low energy distribution of the local potentials E i , there are lots of rearrangements in each update. Those die out only once at least a parabolic pseudo-gap develops in the distribution of the E i 's, such that the probability of a change of occupation triggering other rearrangements becomes small. Note that the E i s are just the classical limit of the HF eigen-energies. Thus the convergence of the HF procedure guarantees essentially the presence of a Coulomb gap in the LDoS. This happens, even though we do not impose explicitly the stability of HF solutions with respect to single particle moves, i.e. swaps between configurations (n i , n j ) = (0, 1) and (1, 0). The latter are the elementary moves considered in standard arguments for the Efros-Shklovskii Coulomb gap in classical Coulomb glasses, but the above reasoning shows that one does not need to impose that extra stability constraint to obtain a well-developed Coulomb gap 9 .
When, on top of the HF equations, stability with respect to particle-hole excitations and more complex rearrangements is imposed, the Coulomb gap hardens a bit, but no essential new features appear [30]. For this reason, we contented ourselves with an analysis of typical solutions of the HF equations, without further minimizing the HF energy among the exponentially many solutions. We expect that the localization properties, multifractality, etc evolve only very weakly as one biases the considered HF solutions toward lower-lying and more-stable solutions. Indeed, our scalings work well when evaluating them in typical solutions on the insulating side, and key physical observables behave as we expected, even in the limit t = 0.

Conclusion
We have studied numerically the localization transition in a 3D Anderson model of spinless fermions, with Coulomb interactions treated within the HF approximation. The MI transition was identified via the localization at the Fermi level, determined from a detailed study of the auto-correlation function of the HF eigenfunctions. Our main results are: (i) multifractal power law scalings in the LDoS survive the presence of interactions, and extend up to a (large) correlation length ξ(U, ε). (ii) A critical Coulomb gap in the weakly insulating phase pins the mobility edge close to ε F for a wide range of parameters, while most higher energy excitations are still delocalized. At disorder strength W = 14 (moderately close, but not finetuned, to the non-interacting critical disorder W = 16.5) the critical U c (ε F ) for the MI transition is ∼ 5 times smaller than the U * required for localization of the entire HF spectrum. This is in qualitative agreement with the experimental observations of Richardella et al [1]. (iii) A scaling analysis of the DoS reveals a critical Coulomb anomaly ρ(ω) ∼ ω 0.6±0. 15 , and scaling laws as anticipated in references [12,23]. (iv) The apparent correlation length exponents display a significant asymmetry between the metallic and insulating sides, similar to tendencies reported in experiments. We conjecture that they arise from crossover phenomena in the metallic phase related to the breakdown of screening or the onset of glassy metastability seen in HF. Those deserve further study.