Quantum Coulomb glass on the Bethe lattice
Izabella Lovas, Annamária Kiss, Cătălin Paşcu Moca, Gergely Zaránd
QQuantum Coulomb glass on the Bethe lattice
Izabella Lovas,
1, 2
Annam´aria Kiss,
3, 4
C˘at˘alin Pa¸scu Moca,
5, 6 and Gergely Zar´and
5, 7 Department of Physics and Institute for Advanced Study,Technical University of Munich, 85748 Garching, Germany Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 M¨unchen Institute for Solid State Physics and Optics, Wigner Research Centre for Physics, POBox 49, H-1525 Budapest, Hungary Department of Physics, Budapest University of Technology and Economics and MTA-BMELend¨ulet Spintronics Research Group (PROSPIN), POBox 91, H-1521 Budapest, Hungary MTA-BME Quantum Dynamics and Correlations Research Group, Institute of Physics,Budapest University of Technology and Economics, Budafoki ´ut 8., H-1111 Budapest, Hungary Department of Physics, University of Oradea, 410087, Oradea, Romani BME-MTA Exotic Quantum Phases ’Lend¨ulet’ Research Group, Institute of Physics,Budapest University of Technology and Economics, Budafoki ´ut 8., H-1111 Budapest, Hungary
We study the Coulomb glass emerging from the interplay of strong interactions and disorderin a model of spinless fermions on the Bethe lattice. In the infinite coordination number limit,strong interactions induce a metallic Coulomb glass phase with a pseudogap structure at the Fermienergy. Quantum and thermal fluctuations both melt this glass and induce a disordered quantumliquid phase. We combine self-consistent diagrammatic perturbation theory with continuous timequantum Monte-Carlo simulations to obtain the complete phase diagram of the electron glass, andto characterize its dynamical properties in the quantum liquid, as well as in the replica symmetrybroken glassy phase. Tunneling spectra display an Efros-Shklovskii pseudogap upon decreasingtemperatures, but the density of states remains finite at the Fermi energy due to residual quantumfluctuations. Our results bear relevance to the metallic glass phase observed in Si inversion layers.
I. INTRODUCTION
Describing the localization of disordered electrons inthe presence of long-ranged Coulomb interactions and themelting of the Coulomb glass due to quantum-fluctuationsrepresent some of the most challenging unsolved and elu-sive problems in modern condensed matter physics [1].In the absence of interactions, disorder tends to suppressquantum fluctuations, and leads to Anderson localiza-tion [2, 3].The presence of interactions, however, changes thestructure of localization transition entirely: unscreenedCoulomb interactions lead to stronger and stronger anoma-lies on the metallic side as one approaches the phasetransition [4], amount in the formation of curious spinfluctuations [5], and lead to the emergence of the Coulombgap [6, 7] on the insulating side, accompanied by glassydynamics and memory effects [8–11]. A major step to-wards understanding this quantum phase transition hasbeen made by A. M. Finkel’stein, who developed a scalingtheory in the presence of Coulomb interactions and weakdisorder [12]. Certain implications of this scaling theoryregarding the critical behavior have been verified exper-imentally [13], but a perturbative scaling theory leavesthe structure of the localized phase unrevealed, and haslittle to say about properties of the localized phase suchas the formation of the pseudogap or the glassy structureof the localized phase, not to mention the connection withmany-body localization [14, 15].The influence of quantum quantum tunneling on theCoulomb gap has been addressed initially by means ofnumerical approaches. A configuration interaction ap-proximation based computation predicted a considerable reduction of the width of the Coulomb gap [16], whileHartree-Fock calculations predicted a modification of thestructure of the classical Efros-Shklovskii pseudogap [6]close to the Fermi surface [17, 18].A major step towards understanding the quantum melt-ing of the Coulomb glass is the construction of a solvablemean field theory, similar to the Sherrington-Kirkpatrick(SK) model, the standard mean field model of classicaland quantum spin glass transitions [19–22]. Such a meanfield model has been proposed by Pastor and Dobrosavl-jevi´c in their seminal work [23], possibly inspired by theextended dynamical mean field approach applied to cleancorrelated systems [24, 25]. In the spinless version of theirmodel, electrons move on a Bethe lattice of coordination z → ∞ , experience some onsite disorder, ε i , and interactwith each other through a repulsive nearest neighbor inter-action, V ij = V / √ z , mimicking the long-ranged Coulomb FIG. 1. Sketch of the mean field electron glass model, Eq. (1).Electrons move on a disordered Bethe lattice of coordinationnumber z , can hop between neighboring sites, and interactthrough a Coulomb interaction with nearest neighbors. a r X i v : . [ c ond - m a t . d i s - nn ] S e p interaction (see Fig. 1),ˆ H = − t √ z (cid:88) (cid:104) i,j (cid:105) (cid:16) ˆ c † i ˆ c j + h . c (cid:17) + V √ z δ ˆ n i δ ˆ n j + (cid:88) i ε i δ ˆ n i . (1)Here δ ˆ n i = ˆ c † i ˆ c i − / ε i are drawn from a Gaussian distribution, P ( ε ) ∼ e − ε / (2 W ) . In the rest of the paper, we shall referto this model as the disordered t–V model .As shown in Ref. [23], even though the interaction isuniform, spontaneous density fluctuations lead to theemergence of a glass phase, and the model (1) mapsonto the Sherrington-Kirkpatrick model in the absenceof quantum tunneling, t = 0. We emphasize that thistransition is structural in the sense that it takes placeeven in the absence of disorder, W →
0, and is driven byinteractions rather than disorder. Here, in contrast to theSK model, frustrations do not originate from a frustratedinteraction: rather, a fluctuation in the occupation ofsome levels creates a ”frustration by choice”, leading tothe glass transition.Later works revealed a number of key properties ofthe disordered t–V model, Eq. (1). The quantum criticalbehavior has been analyzed for small disorder in termsof a Landau theory [26, 27], following a line similar tothe work of Read, Sachdev, and Ye [28], and it has beenargued that for finite coordination numbers, z , a glassymetallic phase should separate the glassy insulating phasefrom the disordered Fermi liquid [29], as observed on lowmobility Si inversion layers [30].Nevertheless, in spite of all these achievements andefforts, a complete solution of (1) is still missing, even inthe mean field limit, z → ∞ . Here we attempt to givesuch an accurate and extensive numerical solution of thedisordered t–V model in the z → ∞ , mean field limit. Inaddition to determining the complete phase diagram, thedistribution of local levels, the order parameter, the freeenergy, and the entropy, we also determine the spectralproperties and the tunneling spectra of the electrons, aswell as their scaling properties away from the criticalpoint.The solution of (1) represents a quite challenging task:to enter the glassy phase and capture the formation ofthe pseudogap, we must allow for complete replica sym-metry breaking, –– accounting for the distribution oflocal (renormalized) energy levels, — and, at the sametime, we must solve an ensemble of quantum impurityproblems coupled self-consistently back to the spin glassorder parameter [31–33]. We derive the appropriate equa-tions using a path integral formalism, and solve themnumerically.We apply two different numerical methods: In theFermi liquid phase, we use an extended continuous timeQuantum Monte-Carlo (CTQMC) dynamical mean fieldapproach. This method provides us the numerically exact,self-consistent solution, however, is numerically demand-ing, and is only appropriate to give us a solution at arelatively small number of points in parameter space. We FIG. 2. (a) Boundary separating the Fermi liquid and theelectron glass phases. Cuts along the solid and dashed linesare presented in Fig. 6. (b) Spectral function computed at thequantum phase transition point, indicated in panel (a) by thearrow. A correlation hole preceding the pseudogap structurestarts to form already at the phase boundary. therefore combine this approach with an Iterative Pertur-bation Theory (IPT), similar in spirit to the one used todescribe the Mott transition in a pioneering work by A.Georges et al. [34]. A combination of these two approachesallows us to obtain a coherent picture, summarized inFig. 2.For convenience, we measure all energy scales in Fig. 2in the disorder strength, W . The electron glass formsat large interactions, and is destroyed both by thermal( ∼ T ) and by quantum ( ∼ t ) fluctuations. Typical spec-tral densities are presented in Fig. 2(b) at a transitionpoint, where T (cid:28) t , and therefore quantum fluctuationsdrive the quantm glass to quantum liquid phase transi-tion. A marked correlation hole structure starts to formalready at the critical point, This anomaly gradually de-velops into a pseudogap that gets deeper and deeper aswe enter the glass phase, but the density of states remainsfinite at the Fermi energy in the glass phase for any finitequantum tunneling, even in the T = 0 temperature limit.This is a peculiarity of the z = ∞ limit, where no Ander-son localization takes place. The glass state we find istherefore identified as a metallic (spinless) electron glass ,observed in several experiments [30, 35, 36]. A similarglassy metallic phase has been predicted to emerge in itin-erant fermionic systems with cavity mediated long-rangeinteractions, based on a replica symmetric effective fieldtheory approach [37].The rest of the paper is organized as follows. We discussthe mean field equations for model (1) in Sec. II. First wepresent the intuitive cavity approach in Sec. II A, then weturn to the more technical replica method in Sec. II B. Wepresent the replica symmetric solution of the mean fieldequations in Sec. III. We first show the self-consistencyequations in Sec. III A, then we provide details on thenumerical solution in Sec. III B. Here we discuss two dif-ferent approaches; the faster but approximate iterativeperturbation theory, and the exact continuous time quan-tum Monte Carlo method. In Sec. III C we study thespectral function in the replica symmetric Fermi liquidphase, and we also determine the phase boundary separat-ing this region from the electron glass phase. We turn tothe properties of the electron glass phase, displaying fullreplica symmetry breaking, in Sec. IV. First we discussthe distribution of the overlaps between different replicasin Sec. IV A, then we turn to the distribution of Hartreeenergies and to the properties of the spectral function inSec. IV B. In Sec. V we discuss the thermodynamic prop-erties of the model, and we show that we find a weaklyfirst order phase transition with latent heat. Finally, wesummarise our main findings in Sec. VI. Technical detailson the continuous time quantum Monte Carlo simulationsand on the numerical solution of the model in the replicasymmetry broken phase, as well as additional numericalresults on the universal scaling of the Coulomb gap arerelegated to appendices. II. MEAN FIELD EQUATIONS
The mean field equations of disordered t–V model (1)have been derived for z → ∞ using the replica methodin Ref. [23]. To our knowledge, these equations have,however, never been solved before in their full power. Infact, to obtain a full solution and to capture the formationof the Coulomb gap, one must consider the structure of full replica symmetry breaking, just as for the Sherrington-Kirckpatrick model [33, 38–40].Before we discuss the quite technical replica method,let us start with a cavity consideration. This allows oneto understand the ultimate structure of the equations tobe solved. A. Cavity approach and effective local action
Let us focus on site i = 0, and write the action corre-sponding to the Hamiltonian (1) as follows: S = (cid:90) τ c τ ( ∂ τ + (cid:15) ) c τ + V √ z (cid:88) i (cid:48) (cid:90) τ δn ,τ δn i,τ − t √ z (cid:88) i (cid:48) (cid:90) τ ( c τ c i τ + h.c. ) + S i (cid:54) =0 . (2)Here we used the shorthand notation, (cid:82) τ = (cid:82) β d τ , with β the inverse temperature, and have separated those pieces,which involve site 0. Primes indicate summations over nearest neighbors only. We can then formally expand e − S , andintegrate out all i (cid:54) = 0 sites to obtain an effective action for site 0, S eff0 = (cid:90) τ c τ ( ∂ τ + (cid:15) ) c τ − t (cid:90) τ (cid:90) τ (cid:48) c τ c τ (cid:48) z (cid:88) i (cid:48) (cid:104) c i τ c i τ (cid:48) (cid:105) cav + V √ z (cid:88) i (cid:48) (cid:90) τ δn ,τ (cid:104) δn i,τ (cid:105) cav − V (cid:90) τ (cid:90) τ (cid:48) δn ,τ δn ,τ (cid:48) z (cid:88) i (cid:48) (cid:104) δn i,τ δn i,τ (cid:48) (cid:105) cav + . . . . (3)Here the (cid:104) . . . (cid:105) cav denote cavity averages, i.e. averages computed in the absence of site i = 0. Higher order contributions,that are not displayed, vanish in the z → ∞ limit on the Bethe lattice. The third term in this expansion represents arandom chemical potential, arising from charge fluctuations at neighboring sites. We can rewrite the above action as S eff0 = (cid:90) τ (cid:90) τ (cid:48) (cid:110) c τ (cid:0) δ ( τ − τ (cid:48) )( ∂ τ + ˜ (cid:15) ) − t G ( τ − τ (cid:48) ) (cid:1) c τ (cid:48) − V χ ( τ − τ (cid:48) ) δn ,τ δn ,τ (cid:48) (cid:111) , (4)with G and χ denoting a nearest neighbour average overcavity Green’s functions and dynamical susceptibilities,and the random field ˜ (cid:15) incorporating charge fluctuationson neighbouring sites into the bare local field, (cid:15) . Since thepresence of site 0 induces a perturbation of order ∼ / √ z on its nearest neighbors, G and χ can be replaced by the lattice average of the local Green’s function and dynamicalcharge susceptibility, respectively. Would we know thedistribution of ˜ (cid:15) , (cid:101) P (˜ (cid:15) ), we could replace these spacialaverages by an average over ˜ (cid:15) . We could thus solve theaction (4) for G ˜ (cid:15) ( τ ) and χ ˜ (cid:15) ( τ ), and obtain G ( τ ) and χ ( τ )by averaging over P (˜ (cid:15) ), thereby closing a dynamical mean field theory (DMFT) cycle.Unfortunately, it is not so simple to obtain (cid:101) P (˜ (cid:15) ). Thedifficulty is related to spontaneous symmetry breaking.Even for a given set of the on site energies, { (cid:15) i } , each’leg’ attached to the cavity has namely many symme-try broken states. We should pick a symmetry brokencharge distribution on each of these legs. However, we cannot choose these independently of each other, sincethe central site i = 0 creates correlations between the legs.Adding/removing the site 0 induces namely a correlatedcharge shift of order ∼ / √ z on neighboring charge dis-tributions. This, in turn, amounts in a change of O (1) inthe value of ˜ (cid:15) and, more importantly, correlates the oc-cupation of neighboring sites through charge fluctuationsat site i = 0. The situation is quite similar to that ofthe ferromagnetic phase of an Ising magnet on the Bethelattice, where the direction of magnetization on each leggets correlated through the central site.The appropriate distribution (cid:101) P (˜ (cid:15) ) follows from stabilitycriteria , usually formulated in terms of the replica method,discussed in the next subsection. We shall also follow this– somewhat formal – route to determine (cid:101) P (˜ (cid:15) ). B. Replica approach
The action S eff0 in Eq. (4) can also be obtained via thereplica trick, whereby we express the logarithm of thepartition function aslog Z = lim n → Z n − n . (5)We therefore take n → z → ∞ limit, where a systematic1 /z cumulant expansion leads to a simple (extended) dy-namical mean field theory structure [41] with the effectiveaction [42] S rep = (cid:90) β dτ (cid:90) β dτ (cid:48) (cid:40) n (cid:88) a =0 (cid:18) c aτ (cid:2) δ ( τ − τ (cid:48) ) ∂ τ (cid:48) − t G ( τ − τ (cid:48) ) (cid:3) c aτ (cid:48) − V χ ( τ − τ (cid:48) ) δn aτ δn aτ (cid:48) (cid:19) − n (cid:88) a (cid:54) = b V Q ab δn aτ δn bτ (cid:48) − n (cid:88) a,b =0 W δn aτ δn bτ (cid:48) (cid:41) . (6)The action (6) is supplemented by the self-consistency conditions, G ( τ − τ (cid:48) ) = (cid:104) c aτ c aτ (cid:48) (cid:105) S rep , χ ( τ − τ (cid:48) ) = (cid:104) δn aτ δn aτ (cid:48) (cid:105) S rep , Q a (cid:54) = b = (cid:104) δn aτ δn bτ (cid:48) (cid:105) S rep . (7)Disorder appears at this point only through the term ∼ W , coupling (correlating) different replicas, and theoff-diagonal structure of the glass order parameter, Q a (cid:54) = b ,capturing density fluctuation correlations between differ-ent replicas subject to the same disorder. This replica-replica coupling in S rep may lead to spontaneous replicasymmetry breaking, characteristic to the glassy phase,and signaling that replicas break ergodicity individuallyand differently. III. THE REPLICA SYMMETRICAL SOLUTION
In general, the non-trivial replica structure of Q ab leadsto difficulties when taking the limit, n →
0. The equationssimplify, however, considerably in the (non-glassy) replicasymmetrical phase, where all replicas behave in the sameway, and Q a (cid:54) = b = Q RS for all a (cid:54) = b . This phase isidentified as a disordered Fermi liquid phase [29]. A. Selfconsistency equations
In this case, we can decouple the off-diagonal part ofthe last term of the effective action (6) with a Hubbard- Stratonovitch field, (cid:101) ε , leading to the local effective action, S (cid:101) ε = (cid:82) τ (cid:82) τ (cid:48) (cid:110) c τ (cid:104) δ τ,τ (cid:48) [ ∂ τ (cid:48) + (cid:101) ε ] − t G ( τ − τ (cid:48) ) (cid:105) c τ (cid:48) − V χ ( τ − τ (cid:48) ) − Q RS ) δn τ δn τ (cid:48) (cid:111) − β (cid:101) ε , (8)with the Hubbard-Stratonovic field ˜ (cid:15) a Gaussian variableof distribution (cid:101) P RS ( (cid:101) ε ) ∼ exp (cid:0) − (cid:101) ε / ( W + V Q RS ) / (cid:1) .The self-consistency equations (7) are now replaced bythe conditions, (cid:26) G ( τ ) χ ( τ ) (cid:27) = (cid:90) d (cid:101) ε (cid:101) P RS ( (cid:101) ε , Q RS ) (cid:26) G (cid:101) ε ( τ ) χ (cid:101) ε ( τ ) (cid:27) , (9)and Q RS is also determined selfconsistently by Q RS = (cid:104) δn (cid:105) = (cid:90) d (cid:101) ε (cid:101) P RS ( (cid:101) ε , Q RS ) (cid:104) δn (cid:105) (cid:101) ε . (10)with G (cid:101) ε ( τ ), χ (cid:101) ε ( τ ), and (cid:104) δn (cid:105) (cid:101) ε computed by the effective(local) action, Eq. (8).In the replica symmetrical case, we thus convertedthe problem into an ensemble of local, self-interactingfermion levels. The width of the distribution of the level ˜ (cid:15) as well as the fermion’s self-energy ( ∼ t G ( τ )) and its self-interaction ( ∼ V [ χ ( τ ) − Q RS ]) must all be determinedself-consistently.Remarkably, the local action has exactly the same struc-ture as (4). However, the replica approach also providesus the self-consistent distribution function ˜ P ( (cid:101) ε ): in thereplica symmetrical Fermi liquid phase, the ”Hartree field”distribution retains the Gaussian structure of the baredisorder ε i , and interactions only renormalize the varianceof the effective field.Importantly, in the classical limit, t = 0, we can set G → δn τ δn τ (cid:48) =1 /
4. Then we simply obtain (cid:104) δn (cid:105) (cid:101) ε = − tanh( (cid:101) ε / (2 T )) / Q RS (cid:54) = 0 is intrinsically unstable, herereplica symmetry is stabilized by finite disorder as well asfinite quantum fluctuations, and a valid replica symmetricphase exists. B. Numerical solution
To solve the action, Eq. (8), i.e., to compute the quan-tities G (cid:101) ε ( τ ), χ (cid:101) ε ( τ ), and (cid:104) δn (cid:105) (cid:101) ε , and then iterate theself-consistency equations Eqs. (9) and (10), we employedtwo different methods: for a fast but approximate solu-tion, we used Iterative Perturbation Theory (IPT), whichallowed us to get a complete solution in the replica sym-metric Fermi liquid phase as well as to access the glassphase. To complement this approach, we have also ob-tained a ”numerically exact” solution by the ContinuousTime Quantum Monte Carlo method (CTQMC), whichwe used to obtain reference solutions at many points in theparameter space, and also to verify the phase boundaries.
1. Iterative Perturbation Theory
The effective action Eq. (8) describes a fermion prop-agating with the unperturbed propagator G (cid:101) ε associatedwith the first term of Eq. (8), G − (cid:101) ε ( τ ) = δ τ,τ (cid:48) [ ∂ τ (cid:48) + (cid:101) ε ] − t G ( τ − τ (cid:48) ) , (11)and interacting through the retarded interaction V ˜ χ ( τ − τ (cid:48) ) ≡ V ( χ ( τ − τ (cid:48) ) − (cid:104) δn (cid:105) ) . To compute all needed Green’s functions and suscepti-bilities in a systematic way, it is worth introducing thelocal (negative) free energy Φ (cid:101) ε as e β Φ (cid:101) ε ( T ) ≡ (cid:90) D c D c e − S loc [ c,c ] , (12)and express it as Φ (cid:101) ε = Φ (0) (cid:101) ε + ∆Φ (cid:101) ε (13) with the second term accounting for the interaction-induced part of Φ (cid:101) ε , and Φ (0) (cid:101) ε being the non-interactingfree energy, Φ (0) (cid:101) ε = (cid:101) ε β Tr log G − (cid:101) ε . (14)The interacting part ∆Φ (cid:101) ε can be considered as a func-tional of the dressed propagator. Then its functionaldifferential with respect to the dressed propagator is justthe self-energy.Within IPT, we simply replace the local free energy(12) by the second order perturbative expression,∆Φ HF (cid:101) ε ( (cid:101) ε ) = V (cid:0) G (cid:101) ε (0 − ) + 1 / (cid:1) (cid:90) β dτ ˜ χ ( τ ) (15) − V (cid:90) β dτ ˜ χ ( τ ) G (cid:101) ε ( τ ) G (cid:101) ε ( − τ ) , represented by the free energy diagrams in Fig. 3(a) [43].Although not constructed in terms of the full Green’sfunction, we shall also refer to this approximation as”Hartree-Fock” approximation, as also inferred by thelabels, ”HF”. For the self-energy, we use a similar approx-imation, represented in Fig. 3(b)Σ HF (cid:101) ε ( τ ) = δ ( τ ) V (cid:0) G (cid:101) ε (0 − ) + 1 / (cid:1)(cid:90) β dτ (cid:48) (cid:101) χ (cid:101) ε ( τ (cid:48) ) − V (cid:101) χ (cid:101) ε ( τ ) G (cid:101) ε ( τ ) . (16)These expressions can also be obtained by functionaldifferentiation of (15) with respect to the unperturbedpropagators. The term 1/2 originates from normal order-ing, and is just the average occupation.Formally, the occupation (cid:104) δn (cid:105) (cid:101) ε and for the local com-pressibility χ (cid:101) ε ( τ ) can be computed by inserting a timedependent energy in the action, (cid:101) ε → (cid:101) (cid:15) τ , and taking thefunctional derivatives of Φ loc with respect to (cid:101) ε → (cid:101) (cid:15) τ .We use this procedure to obtain the IPT expressions for (a)(b) FIG. 3. (a) ”Hartree-Fock” local free energy contributionsand (b) corresponding self-energy diagrams. (Counter-termdiagrams are not shown.) Wavy lines represent the effec-tive interaction, V (cid:101) χ ( τ − τ (cid:48) ), while continuous lines standfor the unperturbed local propagator, G (cid:101) ε ( τ ) = (cid:104) c ( τ )¯ c (0) (cid:105) (0)loc ,computed from the non-interacting part of Eq. (8) . (a)(b) FIG. 4. Hartree-Fock diagrams determining (a) the Hartree-Fock occupation, (cid:104) δn (cid:105) HFloc , and (b) the Hartree-Fock response, χ HF (cid:101) ε ( τ ). Cuts indicate functional derivatives, counter-termdiagrams are omitted. the occupation (cid:104) δn (cid:105) (cid:101) ε and for the local compressibility χ (cid:101) ε ( τ ), consistent with the approximations above, by justdifferentiating Φ HFloc = Φ (0)loc + ∆Φ
HFloc as (cid:104) δn (cid:105) HF (cid:101) ε = β δ ∆Φ HFloc δ (cid:101) (cid:15) τ (cid:12)(cid:12)(cid:12)(cid:12) (cid:101) (cid:15) τ → (cid:101) ε (17)and χ HF (cid:101) ε ( τ − τ (cid:48) ) = β δ Φ HFloc δ (cid:101) (cid:15) τ δ (cid:101) (cid:15) τ (cid:48) (cid:12)(cid:12)(cid:12)(cid:12) (cid:101) (cid:15) τ → (cid:101) ε . (18)The resulting expressions are quite lengthy, we thereforedo not display them here, but the corresponding diagrams,shown in Fig. 4, have a quite transparent structure, andit is easy to construct the explicit formulas from them byfollowing standard diagramatic rules.The IPT iteration is then straightforward in the RSphase. Assuming some ansatz for Q RS , χ ( τ ), and G ( τ )we use Eq. (16) to compute G (cid:101) ε ( τ ), and the diagrams inFig. 4 to determine (cid:104) δn (cid:105) (cid:101) ε and χ (cid:101) ε ( τ ) for a dense set ofenergies, (cid:101) ε . We then determine Q RS , χ ( τ ), and G ( τ )iteratively by means of Eqs. (9) and (10).
2. Continuous time quantum Monte Carlo
An alternative route to compute G ˜ ε ( τ ), (cid:104) δn (cid:105) ˜ ε , and χ ˜ ε ( τ ) within the dynamical mean-field theory is to per-form a continuous-time quantum Monte Carlo (CTQMC)computation with the effective local action S loc givenin Eq. (8). We used an extension of the hybridization-expansion CTQMC algorithm that can treat retardedinteractions in action formalism [44, 45]. In this approachwe expand the partition function Z = Tr e − S loc in the hybridization function F ( τ − τ (cid:48) ) = t G ( τ − τ (cid:48) ) while wetreat the level energies ˜ ε and interaction V exactly.In general, the hybridization-expansion CTQMCmethod[46, 47] relies on the expansion of the partitionfunction Z in the hybridization into a series of diagramsand sampling these diagrams stochastically, where Z canbe written as a sum of configurations z k with weight w ( z k ) as Z = (cid:80) z k w ( z k ). In the segment picture a MonteCarlo configuration z k with expansion order k is rep-resented by k segments with imaginary time intervals { τ , τ (cid:48) } , ..., { τ k , τ (cid:48) k } where the particle number is 1, andit is 0 where there is no segment.In our case the creation operators c τ i at times τ i areconnected to annihilation operators c τ (cid:48) j at times τ (cid:48) j bythe hybridization function F ( τ i − τ (cid:48) j ), and the collectionof these k ! diagrams corresponding to the hybridizationlines F is summed up into a determinant of a matrixˆ F ( k ) composed of the hybridization functions. The weight w ( z k ) is expressed as w ( z k ) = det ˆ F ( k ) w ˜ ε w ˜ χ where thecontributions w ˜ ε and w ˜ χ corresponding to the level energy˜ ε and the interaction term V ˜ χ are given in Eqs. (A7)and (A10) in the appendix, respectively. Further detailsare given in Appendix A.We used this extended CTQMC impurity solver withthe combined weight w ( z k ) by means of the Metropolisalgorithm to solve the effective local action given in Eq. (8)for G ˜ ε ( τ ), (cid:104) δn (cid:105) ˜ ε , and χ ˜ ε ( τ ) self-consistently in the Fermi-liquid (replica symmetric) phase. We proceeded to obtainthe self-consistent replica symmetric solution throughthe following iteration steps: We start by an arbitraryansatz for G [0] ( τ ), ˜ χ [0] ( τ ), and Q [0]RS at the zeroth iterationstep (see Appendix A for our choice), and compute thequantities G [1]˜ ε ( τ ), ˜ χ [1]˜ ε ( τ ), and (cid:104) δn (cid:105) ˜ ε at the first iterationstep with the effective local action, Eq. (8), using theCTQMC impurity solver for a wide range of level energies˜ ε . The averaged quantities G [1] ( τ ), ˜ χ [1] ( τ ), and Q [1]RS are obtained by (numerical) integration over ˜ ε with thedistribution ˜ P (˜ ε ) as it is given in Eqs. (9) and (10). Theyare used for the next, second iteration step, and we repeatthis procedure until we reach convergence.We calculated several points of the phase boundary byCTQMC using the stability condition given in Eq. (20)below for various parameter values for t , V and T , andfound excellent agreement between the IPT and CTQMCcalculations. The spectral functions are also comparedand found to show similar energy dependence betweenIPT and CTQMC. However, difference arises around zeroenergy ω ∼ V . C. Spectral functions and phase boundary
We used both approaches described in the previoussubsection to compute the Green’s function G ( τ ) and thesusceptibility ˜ χ ( τ ). The average local tunneling densityof states (DOS) can then be computed from the Fourier t / W = 0.1, IPT t / W = 0.1, QMC W () t / W = 0.18, IPT t / W = 0.18, QMC / W V / W = 0.455 T / W = 0.01 t = t c critical point t / W = 0.25, IPT t / W = 0.25, QMC FIG. 5. Average low temperature density of states in theFermi liquid phase, as computed by CTQMC and by IPT.The density of states develops a remarkable zero bias anomalyalready in the Fermi liquid, though the distribution of Hartreelevels is still featureless. Increasing quantum fluctuations washaway this correlation hole. transform G ( iω ) as ρ ( ω ) = 1 π Im G ( iω → ω + i + ) . (19)In the very last step, we have used a Pad´e constructionto carry out the analytical continuation.Fig. 5 shows the spectral functions in the Fermi liquidphase for a moderate interaction, V /W = 0 . T /W = 0 .
01, as we drive thesystem closer and closer to the Fermi liquid - electronglass quantum phase transition. For large t , quantumfluctuations destroy the electron glass, and a dirty Fermiliquid is formed. There the density of states is almostfeatureless. As we decrease t/W , quantum fluctuationsget suppressed, and a plasma dip structure starts to formin the middle of the band. IPT and CTQMC are in verygood agreement, and yield both very similar structures.Differences can be attributed to the approximations madewithin IPT, to the limited CTQMC accuracy, and touncertainties related to the analytical continuation.The presence of a plasma dip (or correlation hole) in theFermi liquid phase reflects short range charge correlationsdue to the repulsive interactions between neighboringsites. This correlation hole is a manifestation of Onsager’sback reaction, and it is not directly related to the Efros-Shklovskii gap of the glassy phase [32, 33]. Indeed, in theliquid phase, replica symmetry is maintained, implyingthat the distribution of the renormalized Hartree-Focklevels, ˜ P (˜ (cid:15) ) is still a featureless Gaussian, in contrast tothe tunneling density of states.The boundary of the electron glass phase is determined by a stability (marginality) condition against replica sym-metry breaking. This is essentially identical to the stabil-ity condition appearing in the SK model [23],1 = V (cid:90) d (cid:101) ε (cid:101) P RS ( (cid:101) ε ) χ ( (cid:101) ε ) , (20)with the static local susceptibility defined as χ stat ( (cid:101) ε ) ≡ ∂ (cid:101) ε (cid:104) δn (cid:105) (cid:101) ε .The resulting phase diagram has been presented inFig. 2 for a finite disorder, W . At any temperature and forany hopping t , replica symmetry is broken at interactionslarger than some critical value, V ≥ V C ( T, t, W ). In theclassical limit, t = 0, in particular, an interaction-drivenglass phase emerges at low temperatures for small disorder.Contrary to naive expectations, strong disorder destroys the glassy phase, and leads to a trivial strongly disor-dered phase without replica symmetry breaking (RSB):fluctuations of the bare levels ε i are so large that eachlevel becomes occupied or unoccupied essentially indepen-dently, leaving no room to interaction-induced frustration.For sufficiently strong interaction, however, a Coulombglass phase emerges.The glass phase can be destroyed not only by extremedisorder, but also by thermal and quantum fluctuations,induced by the temperature, T , or the tunneling, t . Thisis demonstrated in the cuts shown in Fig. 6 (indicatedas dashed lines in Fig. 2), where we also compare theCTQMC results with those of IPT. The excellent agree-ment of these two approaches validates the latter, approx-imate method.The role of thermal and quantum fluctuations is notquite identical. In the classical ( t →
0) limit, V class C ∼√ T W , while in the quantum case ( T → ∼ t /W takes over the role of the temperature,and V quant C ∼ t .At finite temperatures, quantum fluctuations and ther-mal fluctuations compete with each other. As demon-strated in Fig. 6(b), at a finite temperatures, small quan-tum fluctuations with Γ (cid:46) T do not change the criticalinteraction strength, V C , and the transition is mostlyinduced by just thermal fluctuations. For Γ (cid:38) T , i.e., t/W (cid:38) (cid:112) T /W , however, quantum fluctuations play thedominant role, as evidenced by the almost linear shift of V C with increasing t/W . IV. ELECTRON GLASS PHASE: FULLREPLICA SYMMETRY BREAKING
In the electron glass phase, replica symmetry is fullybroken. Fortunately, the construction of the previoussection can be generalized to incorporate full replica sym-metry breaking, thereby yielding a complete descriptionof the glassy phase as well. Although derivations mayseem cumbersome, the interpretation of the final resultsis relatively straightforward.The local effective action Eq. (8), supplemented bythe self-consistency conditions Eq. (9), remains unaltered, T / W V c / W IPTQMC 0.0 0.2 0.4 0.6 0.8 1.0 t / W V c / W T / W = 0.05(a) Electron GlassFermi Liquid t / W = 0.1 (b) Electron GlassFermi Liquid
IPTQMC
FIG. 6. Phase boundary between the electron glass and the Fermi liquid, as computed by continuous time Monte Carlo andIPT along the lines indicated in Fig. 2. Interactions lead to the formation of the glass. Quantum fluctuations as well as thermalfluctuations melt the electron gas. expressing that electrons at each site experience a different”Hartree field”, ε i → (cid:101) ε i , due to the conspiracy of randomonsite energies and nearest neighbor Coulomb interactions.Only the ”Hartree field’s” distribution (cid:101) P ( (cid:101) ε ) acquires amore complicated, non-Gaussian structure, that must bedetermined self-consistently together with the averagepropagators and susceptibilities (see Appendix B).The solution of the effective action Eq. (8) is carriedout exactly the same way as in the replica symmetricalphase. Only the last, least intuitive step of this derivation,the determination of the distribution of the renormalized”Hartree” energies (cid:101) P ( (cid:101) ε ) is much more difficult. Thederivation of this distribution follows similar lines as thesolution of the classical spin-glass problem [19], apart fromthe fact that here we need to work with the quantumaction.We parametrize Q ab using Parisi’s variables [21] asa function Q ( x ), with the replica variable x ∈ [0 , family of effective actions, parametrized by x . Physical quantities at different layers are related byso-called ”flow equations”. The final structure of theselatter is outlined in Appendix B).Fortunately, the flow equations are decoupled from thequantum solution in the sense that the quantum problemonly provides boundary conditions for them. In fact, asinput one only needs the (negative) free energy of theembedded level, Φ loc ( (cid:101) ε ) ≡ k B T ln Z loc ( (cid:101) ε ), determined byEq. (12), and approximated within IPT by Eqs. (13), (14),and (15). The solution of the flow equations provides thenbetter and better approximations for the renormalizeddistribution, ˜ P ( (cid:101) ε ), and the order parameter Q ( x ).The exact solution of these self-consistency equations isa demanding task. One first needs to solve the non-localquantum impurity problem, Eq. (8) for a relatively largeset of (cid:101) ε values, extract the expectation values (cid:104) δn (cid:105) (cid:101) ε aswell as the dressed local Green’s functions and suscep- tibilities. Then one needs to solve the above-mentionedflow equations in replica space to update the distribution (cid:101) P ( (cid:101) ε ), compute the average susceptibilities and Green’sfunctions using (9), and then close the cycle by Eq. (9).Although this is, in principle, possible at a given pointin parameter space using, e.g., continuous time quantumMonte-Carlo methods [48], it appears to be unavoidableto use an approximate scheme such as IPT if one aims atdetermining the complete phase diagram. Below, we sum-marize the results of IPT computations. Further CTQMCresults shall be published elsewhere [48]. T / W = 0.1 P ( Q ) T / W = 0.05 Q V / W = 2.0 t / W = 0.5 T / W = 0.008 FIG. 7. Overlap distribution P ( Q ) in the electron glass phase,as a function of temperature. A. Overlap distribution
The differential of the inverse function x ( Q ) turns outto be just the distribution of the overlaps between differentreplicas, Q ab = lim N →∞ N (cid:88) i (cid:104) δn i (cid:105) a (cid:104) δn i (cid:105) b ,P ( Q ab = Q ) = d x d Q . (21)The numerically computed function P ( Q ) is presentedin Fig. 7. In the Fermi liquid phase (not shown), P ( Q )is trivial, and consists of a delta function, P RS ( Q ) = δ ( Q − Q RS ). In the electron glass phase, this distributionbecomes non-trivial, and possible overlaps have a range, Q ∈ [ Q min , Q max ]. This overlap window becomes broaderand broader as the temperature is lowered, and at thesame time, the distribution gets depleted, and has ahight ∼ T . This is in line with the observation, that at T = 0 temperature, replica symmetry is restored. (It is,however, not so clear if a valid expansion around this limitexists [28].) Notice that the maximal value, Q max ( T → /
4; this is a consequence of quantumfluctuations, which tend to reduce the overlaps.
B. Distribution of Hartree-Fock levels andtunneling spectra
As in classical spin glasses [19, 32, 33], a clear signatureof the glass transition is the emergence of a Coulombgap structure in the distribution of Hartree-Fock energies, (cid:101) P ( (cid:101) ε ), shown in Fig. 8. The Coulomb gap starts to openup gradually after crossing the phase transition, and afully developed Coulomb gap appears only deep in theglassy phase. Although the pseudogap gets deeper anddeeper as the temperature decreases, (cid:101) P ( (cid:101) ε = 0) remainsfinite even at T = 0 temperature for any finite t . Thisis a property of the infinite coordination limit, z → ∞ , / W W P () V / W = 2.0 t / W = 0.5 T T C T T C T < T C T / W = 0.200 T / W = 0.156 T / W = 0.133 T / W = 0.067 T / W = 0.050 T / W = 0.025 T / W = 0.008 T / W = 0.005 FIG. 8. Formation of the pseudogap in the distribution (cid:101) P ( (cid:101) ε )as a function of temperature. At T > T C the distributionremains Gaussian, but as we decrease the temperature below T C a pseudogap develops gradually. The thick dashed linerepresents the distribution at the critical temperature T = T C . / W W () T / W = 0.005 t / W = 0.5 V / W = 0.25 V / W = 1.60 V / W = 2.00 V / W = 2.20 V / W = 2.40 V / W = 2.60 V / W = 2.80 / V () FIG. 9. (Main panel) Evolution of the density of states forvarious interaction strengths V . The critical coupling is V c ≈ . W . When V > V C a pseudogap appears in (cid:37) ( ω ), thatgrows with increasing V, but the density of states at the Fermilevel (cid:37) (0) always remains finite. (Inset) Universal collapseof (cid:37) ( ω ). The dashed curve indicates the scaling curve in theclassical limit, t → where Anderson localization is absent, and a disorderedmetal state emerges in the absence of interactions at T = 0 temperature. Nevertheless, interactions largerthan a critical value drive a phase transition to a replicasymmetry broken phase, where the density of states isstrongly suppressed but finite even at T = 0 temperature.We can interpret this phase as a metallic Coulomb glass.While the distribution of the renormalized energies (cid:101) ε is conceptually interesting, excepting the classical limit,their distribution is not directly measurable. What is,however, measurable is the tunneling density of states ata given site of renormalized energy, (cid:101) ε , ρ (cid:101) ε ( ω ) = 1 π Im G (cid:101) ε ( ω + ) , (22)and the average density of states, ρ ( ω ) = Im G ( ω + ) /π , ρ ( ω ) = (cid:90) d (cid:101) ε ρ (cid:101) ε ( ω ) ˜ P ( (cid:101) ε ) . (23)Fig. 9 shows the formation of the pseudogap in ρ ( ω )at very small temperatures, as interactions are increased.The density of states at the Fermi energy is finite, anddefines a natural energy scale ∆ ≡ (cid:37) − (0). This scalebecomes smaller and smaller upon increasing interactions,while ρ ( ω ) develops universal scaling as a function of ω · ∆ /V at low energies, where it crosses over froma constant to a linear regime, ρ ( ω ) ∝ ω/V (see insetin Fig. 9). Notice that the presence of disorder doesnot influence this slope, also indicating that the phasetransition we observe is driven by interactions and notby disorder. The classical scaling function correspondingto t = 0, also displayed in the inset of Fig. 9, yieldsthe same slope as the quantum version, but the twoscaling functions clearly differ, thereby demonstratingthe difference between the role of thermal and quantumfluctuations. As shown in Appendix C, the distribution (cid:101) P ( (cid:101) ε ) exhibits similar universal scaling structure.0 / W W ( , ) T / W = 0.005 t / W = 0.5 V / W = 2.2 / W = 0.00/ W = 0.15/ W = -0.15/ W = 0.50/ W = -0.50/ W = 1.00/ W = -1.00 FIG. 10. Unaveraged density of states deep in the glass phasefor various values of (cid:101) ε . It is instructive to investigate the structure of individualtunneling spectra, ρ (cid:101) ε ( ω ), shown for a set of levels deep inthe quantum glass regime in Fig. 10. The local density ofstates displays peaks at around the renormalized level, (cid:101) ε ,which is broadened by quantum fluctuations. Levels closeto the Fermi level become sharp since surrounding siteshave a suppressed density of states at the pseudogap. V. THERMODYNAMICS
To determine the free energy of the glass we first needto determine the (negative) free energy density Φ loc ( T )of the effective replica action S rep , Eq. (6),Φ loc ( T ) ≡ lim n → n k B T log (cid:26) (cid:90) D c D c e − S rep [ c,c ] (cid:27) . (24)This is slightly different from the physical free energydensity of the lattice model, Eq. (1), which we denote byΦ phys ( T ), since we must restore some terms that we threwaway in course of the Hubbard-Stratonovic transformation.Restoring these terms, which depend on the local Green’sfunction and susceptibility, we obtainΦ latt ( T ) = Φ loc ( T ) + t (cid:90) β d τ G ( τ ) G ( − τ ) − V (cid:90) β d τ χ ( τ ) − βV n (cid:88) b (cid:54) = a Q ab Q ba . (25)In the replica symmetrical (Fermi liquid) case, Eq.(24)simplifies, and we obtainΦ RSloc ( T ) = (cid:90) d (cid:101) ε (cid:101) P RS ( (cid:101) ε ) Φ loc ( (cid:101) ε , T ) , (26)with the free energy Φ loc ( (cid:101) ε ) computed from the localeffective action, Eq. (8), and (cid:101) P RS the Gaussian Hartreelevel distribution, displayed below Eq. (8). In this case,the last term of (25) also simplifies to − n βV (cid:88) b (cid:54) = a Q ab Q ba → β V Q T / W ( T ) t / W = 0.75V / W = 2.1 T T C RSRSB
FIG. 11. Temperature dependence of the entropy for thereplica symmetric (RS) and replica symmetry beaking (RSB)solution. The replica symmetry breaking solution seems todisplay a first order entropy jump at the phase transition, andan entropy scaling to zero as T → in the n → loc ( T ) becomes more complex, since one cannot decouplereplicas with a single Hubbard-Stratonovich transforma-tion. One still has to solve local effective action S loc ( (cid:101) ε )at the start, (with ˜ χ = χ − Q RS replaced by χ − Q (1)),and solve the so-called flow equations in replica space(see Appendix B for details) to obtain an expression forΦ loc ( T ) analogous to Eq. (26).Once Φ loc ( T ) and the converged Green’s functions andthe order parameter Q ( x ) at hand, the thermodynamicquantities of the lattice model can then be computed fromΦ latt ( T ). In particular, we determined the temperaturedependent entropy density S ( T ), given by S = ∂ Φ latt ∂T . (27)Our results for S ( T ) are displayed the temperature de-pendent entropy S ( T ) in Fig. 11.Though the distribution (cid:101) P is apparently continuousthrough the Coulomb glass phase transition, the entropyshows a clear jump, and signals a first order transitionwith a latent heat. We note that this first order transitioncontradicts previous studies applying Landau theory toexamine the glass transition, assuming a continuous phasetransition [26, 27]. However, distinguishing a weakly firstorder transition from a continuous one is an extremelydifficult problem, and given the approximations involvedin our calculation, more evidence would be needed toconfirm our findings.Although we cannot decrease the temperature verydeep down into the RSB phase, but the numerical dataare consistent with the entropy remaining positive andgoing quadratically to zero, and a corresponding quadraticspecific heat.1 VI. DISCUSSION
We presented here a detailed study of the mean fieldCoulomb glass (disordered t–V) model of Ref. [23] in thequantum regime, in the Fermi liquid (replica symmetrical)as well as deep in the glassy (replica symmetry) brokenphase. The combination of continuous time quantumMonte-Carlo approach with Iterative Perturbation Theory(IPT) allowed us to map accurately the phase boundariesseparating the interaction induced glassy phase from theFermi liquid phase in the classical as well as in the quan-tum regime, and to determine spectral functions as wellas thermodynamic properties. Having validated our IPTscheme in the metallic regime, we used it to enter theelectron glass phase, where complete replica symmetrybreaking must be incorporated in the theory.In the spectral function, we observe the formation of aplasmonic correlation hole in the average tunneling density(local density of states) already in the Fermi liquid phase.This correlation hole smoothly develops into an Efros-Shklovskii pseudogap when we enter the electron glassphase. The Efros-Shklovskii pseudogap gap is, however,not fully developed in this mean field model even at T = 0 temperature: similar to thermal fluctuations, smallquantum fluctuations induce a finite density of states evenat the Fermi energy. This is a peculiarity of the infinitecoordination limit, where localization is absent, and aglassy Fermi liquid state emerges rather than a glassylocalized phase. For small tunneling, the average densityof states exhibits universal scaling at at low temperaturesand low energies.We have also computed the local density of states. Inthe electron glass phase, this typically consists of sharpresonances, located around some renormalized Hartreeenergies, whose distribution also exhibits a pseudogap.These resonances become sharper and sharper as oneapproaches the Fermi energy, but retain their finite width,even at the Fermi surface, indicating again that thesestates remain extended even in the glass phase.We have also constructed the full thermodynamic meanfield description of the disordered t–V model, and ana-lyzed its thermodynamic properties in both phases. Ourpresent results are not yet accurate enough to access thecritical behavior and to confirm the results of a Landaufunctional approach [26, 27]. Nevertheless, while we doobtain a continuous free energy, our numerical resultsare more consistent with a jump in the entropy, and aresuggestive of a first order transition into the quantumglass phase, clearly conflicting with a Landau approach.The observed behavior of the entropy, however, may wellbe a numerical artifact of IPT. Further, very accurate con-tinuous time quantum Monte-Carlo computations wouldbe necessary within the replica symmetry broken phaseto resolve (or confirm) this controversy.As mentioned above, the absence of the insulating elec-tron glass phase is an artifact of the infinite coordination limit. However, the unavoidably metallic glass phaseemerging in this model is relevant for many metallic disor-dered systems, which exhibit glassy behavior. Amorphouspolycrystaline solids [11] or granular metals [49] are suchexamples, but metallic electron glass phases can be ob-served in certain two dimensional systems [30, 50]. Thor-ough studies of Na + doped silicon MOSFETs reveal ametal insulator transition at a carrier density, n ≈ n c , andan intermediate metallic glass phase emerges on the metal-lic side of the transition at concentrations n g > n > n c ,as evidenced by low frequency resistance noise [50] andageing [51] experiments on low mobility samples. A metal-lic glass phase could also be experimentally realized inultracold atomic settings, by placing fermionic atoms intoa multi-mode cavity [37].The understanding and solution of the mean fieldCoulomb glass model, Eq. (1), is just a first step in con-structing a mean field theory of the real Coulomb glass.In fact, it is quite unclear, how one could incorporate lo-calization and long-ranged interactions at the same timein a mean field model. To have an Anderson-localizedphase, one should impose a finite coordination number, z , and thereby exclude the presence of infinite numberof nearest neighbors. Of course, similarly to dynamicalmean field theory (DMFT) and the coherent potentialapproximation (CPA), one could use the present schemeas a local approximation to describe the glassy phasetransition in a finite dimensional system [33].Another open question is that of glassy dynamics.Global charge response should reflect the emergence of aglassy phase through anomalously slow response and abroad distribution of scales [9, 10]. It remains an openquestion, how the present approach is able to explain thisbehavior. Finally, spin degrees of freedom have been com-pletely neglected in this work. The role of Mott-Andersonphysics and spontaneous spin formation should be furtherexplored and elucidated.Although in this work we only focused on the descrip-tion of the Coulomb glass phase, the method and formal-ism presented pave the way to study quantum correlationsin the glassy phase of many mean field quantum glass mod-els, such as the transverse field Sherrington-Kirkpatrickmodel, and the disordered Dicke model. All these ques-tions are and should be subject of future research. ACKNOWLEDGMENTS
We thank Vladimir Dobrosavljevi´c and Pascal Simon forinsightful discussions. This work has been supported bythe National Research, Development and Innovation Of-fice (NKFIH), through Grant No. K124176, and throughthe Hungarian Quantum Technology National ExcellenceProgram, project no. 2017-1.2.1-NKP-2017-00001, by theBME-Nanotechnology FIKP grant (BME FIKP-NAT)and by the European Research Council (ERC) under theEuropean Union’s Horizon 2020 research and innovationprogramme (grant agreement No. 771537).2 [1] M. Pollak, M. Ortu˜na, and A. Frydman,
The ElectronGlass (Cambridge University Press, 2013).[2] F. Evers and A. D. Mirlin, Rev. Mod. Phys. , 1355(2008).[3] B. Kramer and A. MacKinnon, Reports on Progress inPhysics , 1469 (1993).[4] B. L. Altshuler, A. G. Aronov, and P. A. Lee, Phys. Rev.Lett. , 1288 (1980).[5] M. A. Paalanen, S. Sachdev, R. N. Bhatt, and A. E.Ruckenstein, Phys. Rev. Lett. , 2061 (1986).[6] A. L. Efros and B. I. Shklovskii, Journal of Physics C:Solid State Physics , L49 (1975).[7] M. Amini, V. E. Kravtsov, and M. M¨uller, New Journalof Physics , 015022 (2014).[8] J. Jaroszy´nski and D. Popovi´c, Phys. Rev. Lett. ,046405 (2007).[9] A. Amir, Y. Oreg, and Y. Imry, Annual Review of Con-densed Matter Physics , 235 (2011).[10] A. Vaknin, Z. Ovadyahu, and M. Pollak, Phys. Rev. B , 134208 (2002).[11] Z. Ovadyahu, Phys. Rev. Lett. , 226603 (2007).[12] A. M. Finkel’stein, Zeitschrift f¨ur Physik B CondensedMatter , 189 (1984).[13] R. J. Zeches, M. D. Rossell, J. X. Zhang, A. J. Hatt,Q. He, C.-H. Yang, A. Kumar, C. H. Wang, A. Melville,C. Adamo, G. Sheng, Y.-H. Chu, J. F. Ihlefeld, R. Erni,C. Ederer, V. Gopalan, L. Q. Chen, D. G. Schlom, N. A.Spaldin, L. W. Martin, and R. Ramesh, Science , 977(2009).[14] R. Nandkishore and D. A. Huse, Annual Review of Con-densed Matter Physics , 15 (2015).[15] K. S. Tikhonov and A. D. Mirlin, Phys. Rev. B , 214205(2018).[16] G. Vignale, Phys. Rev. B , 8192 (1987).[17] F. Epperlein, M. Schreiber, and T. Vojta, Phys. Rev. B , 5890 (1997).[18] T. Vojta, T. Vojta, F. Epperlein, and M. Schreiber,physica status solidi (b) , 53 (1998).[19] D. Panchenko, The Sherrington-Kirkpatrick Model (Springer, 2000).[20] K. Binder and A. P. Young, Rev. Mod. Phys. , 801(1986).[21] M. Mezard, G. Parisi, and M. A. Virasoro, Spin glasstheory and beyond (World Scientific, 1987).[22] K. H. Fisher and J. A. Hertz,
Spin glasses (CambridgeUniversity Press, 1991).[23] A. A. Pastor and V. Dobrosavljevi´c, Phys. Rev. Lett. ,4642 (1999).[24] Q. Si and J. L. Smith, Phys. Rev. Lett. , 3391 (1996).[25] R. Chitra and G. Kotliar, Phys. Rev. Lett. , 3678(2000).[26] D. Dalidovich and V. Dobrosavljevi´c, Phys. Rev. B ,081107 (2002).[27] L. Arrachea, D. Dalidovich, V. Dobrosavljevi´c, and M. J.Rozenberg, Phys. Rev. B , 064419 (2004).[28] N. Read, S. Sachdev, and J. Ye, Phys. Rev. B , 384(1995).[29] V. Dobrosavljevi´c, D. Tanaskovi´c, and A. A. Pastor,Phys. Rev. Lett. , 016402 (2003).[30] S. c. v. Bogdanovich and D. Popovi´c, Phys. Rev. Lett. , 236401 (2002). [31] M. M¨uller and L. B. Ioffe, Phys. Rev. Lett. , 256403(2004).[32] S. Pankov and V. Dobrosavljevi´c, Phys. Rev. Lett. ,046402 (2005).[33] M. M¨uller and S. Pankov, Phys. Rev. B , 144201 (2007).[34] A. Georges and G. Kotliar, Phys. Rev. B , 6479 (1992).[35] P. V. Lin, X. Shi, J. Jaroszynski, and D. Popovi´c, Phys.Rev. B , 155135 (2012).[36] B. H. Moon, J. J. Bae, M.-K. Joo, H. Choi, G. H. Han,H. Lim, and Y. H. Lee, Nature Communications , 2052(2018).[37] M. M¨uller, P. Strack, and S. Sachdev, Phys. Rev. A ,023604 (2012).[38] A. Crisanti and T. Rizzo, Phys. Rev. E , 046137 (2002).[39] R. Oppermann and D. Sherrington, Phys. Rev. Lett. ,197203 (2005).[40] S. Pankov, Phys. Rev. Lett. , 197204 (2006).[41] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996).[42] Note that Eq. (6) neglects certain (temperature de-pendent) constant terms resulting from integrating outfermions at all lattice sites except one. These terms, notinvolving the Grassmann variables c aτ and c aτ , do not enterthe correlators of the local effective model, however, theywill become important when we study thermodynamicproperties of the full lattice Hamiltonian.[43] For clarity, we omitted diagrams incorporating thecounter-terms arising from the chemical potential of thehalf-filled model in Figs. 3 and 4. These diagrams can beconstructed by replacing propoagator loops G (cid:101) ε (0 − ) witha constant factor -1/2.[44] T. Ayral, S. Biermann, and P. Werner, Phys. Rev. B ,125149 (2013).[45] P. Werner and A. J. Millis, Phys. Rev. Lett. , 146401(2010).[46] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, andA. J. Millis, Phys. Rev. Lett. , 076405 (2006).[47] P. Werner and A. J. Millis, Phys. Rev. B , 155107(2006).[48] A. Kiss, I. Lovas, C. P. Moca, and G. Zarand, unpub-lished.[49] J. Delahaye, T. Grenet, and F. Gay, The EuropeanPhysical Journal B , 5 (2008).[50] J. Jaroszy´nski, D. Popovi´c, and T. M. Klapwijk, Phys.Rev. Lett. , 276401 (2002).[51] J. Jaroszy´nski and D. Popovi´c, Phys. Rev. Lett. ,046405 (2007).[52] H. J. Sommers and W. Dupont, Journal of Physics C:Solid State Physics , 5785 (1984). Appendix A: Details of continuous time quantumMonte Carlo method
In this appendix we present the derivation of the com-bined Monte Carlo weight w ( z k ) = det ˆ F ( k ) w ˜ ε w ˜ χ basedon the implementation presented in Ref. [44].We separate the effective local action given in Eq. (8)3as S loc ≡ S F + S , (A1) S F = − (cid:90) τ (cid:90) τ (cid:48) c τ t G ( τ − τ (cid:48) ) c τ (cid:48) , (A2) S = (cid:90) τ c τ ( ∂ τ + ˜ ε ) c τ − V (cid:90) τ (cid:90) τ (cid:48) ( χ ( τ − τ (cid:48) ) − Q RS ) δn τ δn τ (cid:48) , (A3)and expand the partition function Z = Tr e − S loc in termsof the hybridiztion part S F , which gives Z = Tr e − ( S F + S ) = (A4)= (cid:88) k (cid:90) β dtτ ...dτ k (cid:90) β dtτ (cid:48) ...dτ (cid:48) k det ˆ F ( k ) × Tr (cid:104) e − S c τ c τ (cid:48) ...c τ k c τ (cid:48) k (cid:105) = (cid:90) D ( k ) w ( z k ) , (A5)where we introduced the notation (cid:82) D = (cid:80) k (cid:82) β dtτ ...dτ k (cid:82) β dtτ (cid:48) ...dτ (cid:48) k , and therefore the weight w ( z k ) is expressed as w ( z k ) = det ˆ F ( k ) Tr (cid:104) e − S c τ c τ (cid:48) ...c τ k c τ (cid:48) k (cid:105) = det ˆ F ( k ) (cid:104) c τ c τ (cid:48) ...c τ k c τ (cid:48) k (cid:105) . (A6)By evaluating the first term of S in Eq. (A3) in thesegment picture, we obtain the weight w ˜ ε as w ˜ ε = e − ˜ ε l , (A7)where l = (cid:80) ki =1 l i is the total length of the segments l i with l i = τ (cid:48) i − τ i .The weight contribution from the second term inEq. (A3) is expressed as w ˜ χ = e V (cid:82) β dτ (cid:82) β dτ ˜ χ ( τ − τ ) δn τ δn τ (cid:48) . (A8)Defining a function K ( τ ) as K ( τ ) (cid:48)(cid:48) = V ˜ χ ( τ ) with theconditions K (0) = 0 and K (cid:48) (0) = 1 / V (cid:82) β dτ ˜ χ ( τ ), theintegral in Eq. (A8) is evaluated as w ˜ χ = exp (cid:18) (cid:88) k ,k (cid:2) − K ( τ (cid:48) k − τ (cid:48) k ) + K ( τ k − τ (cid:48) k )+ K ( τ (cid:48) k − τ k ) − K ( τ k − τ k ) (cid:3) + K (cid:48) (0) l (1 − (cid:104) n (cid:105) ) (cid:19) , (A9)which can be rewritten as w ˜ χ = exp (cid:18) − (cid:88) i>j s i s j [ K (˜ τ i − ˜ τ j ) − K (0)]+ K (cid:48) (0) l (1 − (cid:104) n (cid:105) ) (cid:19) = exp − (cid:88) i>j s i s j K (˜ τ i − ˜ τ j ) , (A10) where the times are ordered as 0 < ˜ τ < ˜ τ < ... < β ,and s is +1 for a creation operator and − w ( z k ) can finally be expressed in thecompact form w ( z k ) = det ˆ F ( k ) w ˜ ε w ˜ χ . We note that quan-tities can be measured without additional computationcost with this modified weight.The Green’s function G ( τ ) = (cid:104) T τ c (0) c ( τ ) (cid:105) is evaluatedin the same way in our Monte Carlo procedure as in theabsence of the retarded interaction V ˜ χ ( τ − τ (cid:48) ). Namely,for measuring the Green’s function we need a configura-tion where operators c τ (cid:48) and c τ are unconnected. In thehybridization method it is achieved by removing one ofthe hybridization lines which gives G ( τ ) = (cid:42) β k (cid:88) i,j (cid:16) ˆ F ( k ) (cid:17) − ji δ ( τ (cid:48) i − τ j + τ ) (cid:43) MC (A11)with using the hybridization matrix F ( k ) i,j = F ( τ i − τ (cid:48) j ).Both averaged susceptibilities χ ( τ ) = (cid:104) δn τ δn (cid:105) and˜ χ ( τ ) = (cid:104) n τ n (cid:105) − (cid:104) n (cid:105) can be sampled in the Monte Carlosimulation, and therefore the properties χ ( τ = 0) = 1 / χ ( τ = 0) = 1 / − q RS , and ˜ χ ( τ ) = χ ( τ ) − q RS can becheck-points for the correctness of the CTQMC code.Our choice for the zeroth order ansatz for G [0] ( τ ),˜ χ [0] ( τ ), and Q [0]RS in obtaining the self-consistent replicasymmetric solution are the non-interacting ones as G [0] ( iω n ) = i t (cid:16) ω n − (cid:112) t + ω n (cid:17) FFT −−−→ G [0] ( τ ) , (A12)˜ χ [0] ( τ ) = G [0] ( τ ) G [0] ( − τ ) , (A13) Q [0]RS = 0 . (A14)We note that the choice of the ansatz does not affect theconverged result. Appendix B: Replica symmetry breaking
In the glassy phase, replica symmetry is broken, and Q ab acquires a non-trivial structure in replica space. Inthe limit n →