Quantum-embedding description of the Anderson lattice model with the ghost Gutzwiller Approximation

We present benchmark calculations of the Anderson lattice model based on the recently-developed"ghost Gutzwiller approximation". Our analysis shows that, in some parameters regimes, the predictions of the standard Gutzwiller approximation can be incorrect by orders of magnitude for this model. We show that this is caused by the inability of this method to describe simultaneously the Mott physics and the hybridization between correlated and itinerant degrees of freedom (whose interplay often governs the metal-insulator transition in real materials). Finally, we show that the ghost Gutzwiller approximation solves this problem, providing us with results in remarkable agreement with dynamical mean field theory throughout the entire phase diagram, while being much less computationally demanding. We provide an analytical explanation of these findings and discuss their implications within the context of ab-initio computation of strongly-correlated matter.

We present benchmark calculations of the Anderson lattice model based on the recently-developed "ghost Gutzwiller approximation". Our analysis shows that, in some parameters regimes, the predictions of the standard Gutzwiller approximation can be incorrect by orders of magnitude for this model. We show that this is caused by the inability of this method to describe simultaneously the Mott physics and the hybridization between correlated and itinerant degrees of freedom -whose interplay often governs the metal-insulator transition in real materials. Finally, we show that the ghost Gutzwiller approximation solves this problem, providing us with results in remarkable agreement with dynamical mean field theory throughout the entire phase diagram, while being much less computationally demanding. We provide an analytical explanation of these findings and discuss their implications within the context of ab-initio computation of strongly-correlated matter.
Understanding and simulating quantitatively the electronic behavior of strongly correlated matter is one of the most fundamental problems in condensed-matter science. The substantial progress achieved today in this direction largely owes to quantum embedding methods [1,2]. In particular, the development of dynamical mean field theory (DMFT) [3][4][5][6][7][8][9][10][11][12][13] constituted a great leap in our understanding of strong-correlation phenomena, which advanced dramatically our ability of describing the properties of real materials. In the past decade, the perspective of expanding the predictive power of simulations within the blooming field of theory-assisted materials-bydesign [14,15] contributed to stimulate the development of alternative computational frameworks, capable of taking into account strong correlations at a lower computational cost. Within this context, particularly promising approaches are the Gutzwiller approximation (GA) [16][17][18][19][20][21] -or, equivalently [22,23], the rotationally-invariant slave-boson mean-field theory [24][25][26] -and density matrix embedding theory [27,28]. These frameworks have similar algorithmic structures. In fact, as in density matrix embedding theory, the GA equations can be cast in terms of ground-state calculations of auxiliary impurity models called embedding Hamiltonians (EH), where the bath has the same number of degrees of freedom as the impurity [20]. Furthermore, density matrix embedding theory can be formally derived from the GA equations, setting to unity the parameters encoding the quasiparticle mass-renormalization weights [29,30]. More recently, a more accurate extension of the GA, called "ghost Gutzwiller approximation" (g-GA) has been de-veloped [31], based on the idea of extending the GA variational space introducing auxiliary (ghost) fermionic degrees of freedom.
Here we present benchmark calculations of the Anderson lattice model (ALM) and demonstrate that, by construction, the GA cannot capture the interplay between Mott physics and the hybridization between correlated and itinerant degrees of freedom -which generally coexist and whose interplay often governs the metal-insulator transition in real materials. We also show that, in some parameters regimes, this limitation of the GA can result in overestimating the Mott critical point by orders of magnitude. Finally, we demonstrate, both numerically and analytically, that the g-GA method resolves these problems, while remaining much less computationally demanding than DMFT. Furthermore, we show that this method allows us to describe semi-analytically the spectral properties (both at low and high energies) throughout the entire phase diagram of the ALM, facilitating the physical interpretation of the numerical results.
Model.-We consider the ALM on a Bethe lattice, in the limit of infinite coordination number [32]: where: p iσ and d iσ are Fermionic annihilation operators, p † iσ and d † iσ are Fermionic creation operators, i and j are site labels, σ is the spin, < i, j > indicates that the corresponding summation is restricted to first nearest neighbours, the hopping matrix t ij is uniform, µ is the chemical potential, the columns of U are the eigenvectors of t and k are its eigenvalues. From now on we fix the hopping matrix t by using the half-bandwidth of k (corresponding to a semi-circular density of states) as the energy unit.
Method.-Here we summarize the algorithmic structure of the g-GA and the GA, pointing out the key differences between these two methods, from a quantumembedding perspective [20,26,31]. For simplicity, below we focus on the ALM introduced above, while the general theory for arbitrary multi-orbital systems is summarized in the supplemental material [33]. For both the g-GA and the GA, the solution is obtained by calculating recursively the ground state of 2 auxiliary systems: (1) the so-called "quasiparticle Hamiltonian" (QPH) and (2) the "embedding Hamiltonian" (EH).
The EH can be expressed in the following form: where thed operators correspond to the EH impurity degrees of freedom,n d = σd † σdσ , thef operators correspond to the EH bath degrees of freedom and the parameters D and λ c are determined self-consistently [33]. Note that, while in standard GA the bath of the EH has the same size of the impurity, within the g-GA it contains a larger number of sites (B > 1). As in Refs. [31,34], here we will set B = 3 (the effect of increasing B, which would enlarge further the variational space, will be subject of future work). After convergence, the expectation value of any local operatorÔ[d † iα , d iα ] can be calculated from the ground state |Φ of the EH as: The QPH can be expressed as: where the f iaσ operators are "quasiparticle modes" residing in an auxiliary (enlarged) Hilbert space [33]. Once the parameters l and r are determined self-consistently in the form above [33], the resulting Green's function for the modes φ kσ is: where the only non-zero entry of Σ(ω) is the dd component, which is given by the following equations: As demonstrated in the supplemental material [33] (and shown in the calculations below): (i) the poles of G(k, ω) are located on top of the eigenvalues of the corresponding QPHs [35][36][37][38], and (ii) the resulting total spec-tral weight of the d degrees of freedom is given by: where A g-GA dd and A GA dd are the g-GA and GA d-electron spectral functions, respectively, and r 2 is the GA delectron quasiparticle weight [39]. From the parameters R, λ and the ground state |Ψ 0 of the corresponding QPH, it is also possible to calculate the expectation values of all non-local quadratic operators. In particular: For completeness, the generalization of all analytical results listed above to arbitrary multi-orbital systems is given in the supplemental material [33].
Results.-In Fig. 1 we show the g-GA phase diagram of the ALM (in the paramagnetic phase) for total occupancy N i = 3. The g-GA results are compared with DMFT -with the continuous time quantum Monte Carlo (CTQMC) impurity solver [40][41][42], at temperature T = 0.01-and with the bare GA. Our benchmark calculations show that the g-GA phase diagram is consistent with previous work [43][44][45] and in remarkable agreement with DMFT. As expected, both the g-GA method and the bare GA capture the fact that the Mott metalinsulator transition point U c2 vanishes for ε p → −∞corresponding to the limit where the p degrees of freedom are gapped out. However, the interaction U c of the metal-insulator transition is largely overestimated within the GA, especially for p −1. Note that the phase di-agram for N i = 1 can be automatically inferred from our calculations above, as they are related to each other by a particle-hole transformation.
In Fig. 2 we show the behavior of the g-GA total energy E, the occupancies n p = n pi and n d = n di , the d-electron double occupancy D = n di↑ndi↓ , the p-d "hybridization energy" H = σ d † iσ p iσ + c.c., and the delectron quasiparticle weight: .
The g-GA results are shown in comparison with the bare GA and DMFT. While the GA solution is considered accurate only for weak interactions, we find that the agreement between g-GA and DMFT is remarkable for all observables, in all parameters regimes. Let us now analyze the single-particle Green's function. In Fig. 3 we consider p = −1 and 3 values of U , showing the total g-GA energy-resolved spectral function: and the p and d local density of states. The DMFT spectra were obtained by performing analytical continuation with the maximum entropy method [46]. Interestingly, the g-GA captures systematically the main features of the DMFT spectra (including the Hubbard bands, and the hybridization between the p and d degrees of freedom). Note that, to interpret the behavior of the g-GA spectra, it is possible to exploit its relation with the bands of the QPH [Eq. (6)], previously discussed in the methods section. In particular, the relative position of the p band with respect to the Fermi level is approximately encoded in the corresponding on-site energy * p = p − µ, while the positions of the d-electron low-energy and high-energy excitations are approximately encoded in the variational parameters l a (a = 1, 2, 3).
A key fact emerging from our benchmark calculations is that the GA can overestimate U c dramatically (especially for p −1), see Fig. 1. To explain this result we note that, by construction, the GA Mott transition occurs when the quasiparticle weight Z = r 2 vanishes, see Eq. (10). Therefore, the corresponding approximation to the Mott phase is such that d † iσ p iσ GA = 0, see Eq. (14). In other words, this method cannot describe simultaneously the Mott phase and the p-d charge fluctuations. But this is unrealistic for the ALM, where the p-d hybridization effects are generally very large, not only in the weakly-interacting regime, but also for U U c2 and U > U c2 , see Figs. 2,3. Because of the variational principle, this results in a systematic overestimation of the total energy as we approach the Mott phase (see Fig. 2), causing an overestimation of the metal-insulator transition point. This point is documented in further detail in the supplemental material, where the GA overestimation of U c at p −1 is explained in relation to a qualitative pathological behavior of the GA variational parameter r in the narrow-bandwidth limit (t → 0).
Remarkably, since the g-GA captures the existence of the Hubbard bands, the right side of Eq. (11) never vanishes [31]. Therefore, similar to DMFT, d † iσ p iσ g-GA remains finite even in the Mott phase (see Eq. (13)). This shows that the ability of the g-GA of describing simultaneously the Mott physics and the d-p hybridization is directly connected with its ability of describing the transfer of d-electron spectral weight to the Hubbard bands (which the bare GA lacks).
Conclusions.-We performed benchmark calculations of the ALM, showing that the g-GA provides us with results with accuracy comparable to DMFT, both for the ground-state and the spectral properties. In particular, we showed that the g-GA is capable of describing accurately the interplay between the Mott physics and the hybridization between correlated and itinerant degrees of freedom, while the GA cannot describe simultaneously these effects. This is particularly relevant for real-material calculations in combination with density functional theory [17,20,47,48], where the correlated orbitals are generally very localized around their atomic positions [49,50] -so that the interactions with their environment are mainly mediated by the itinerant modes (as in the ALM studied here). In fact, it is well possible that the limitation of the GA here uncovered explains why, in some cases, simulating the properties of real materials with the GA requires to use unphysically-large Hubbard U [51], and suggests that multi-orbital implementations of the g-GA will resolve these problems.
From the computational standpoint, the g-GA is more expensive than the bare GA (as the bath of the EH contains additional degrees of freedom). On the other hand, its computational complexity remains much lower than DMFT. In fact, the g-GA requires to calculate only the ground state of a finite-size impurity model (while in DMFT it is necessary to calculate the spectra of an impurity model with an infinite bath). For example, in our DMFT calculations each CTQMC iteration (performed using 5 × 10 8 Monte Carlo steps, in parallel, on 72 cores) required about 2 minutes of computational time. In-stead, within our g-GA calculations each EH iteration (performed on a single core) required about 0.2-0.3 seconds. Note also that the difference in computational complexity between DMFT and g-GA grows exponentially as a function of the impurity size.
A particularly promising perspective is the possibility of solving the g-GA equations with hybrid quantumclassical frameworks [52], employing impurity solvers based on quantum algorithms such as variational quantum eigensolvers [53][54][55][56]. In fact, within the g-GA, realizing such program for real-material applications may require devices consisting of only tens of qubits, while it has been estimated that quantum computers with at least 100 logical qubits will be necessary for applications within DMFT [57]. Furthermore, since the number of parameters characterizing the EH is finite, the recentlydeveloped approach based on machine learning for the GA [58] will be applicable also to the g-GA, as we hope to show in future work. Note also that the g-GA can be equivalently formulated in terms of the rotationallyinvariant slave-boson mean-field theory [31], which is based on an exact reformulation of the many-body problem. This line of interpretation may open the possibility of developing beyond-mean-field schemes, providing us with new routes for high-precision calculations.
For completeness, here we provide a comprehensive derivation of the g-GA method, presenting the theory from a slightly different perspective with respect to Ref. [1] (more direct and concise, but equivalent).

The Hamiltonian
We consider a generic multi-orbital Fermionic Hamiltonian represented as follows: where R indicates the unit-cell label, k is the corresponding crystal-momentum, while i and α classify the Fermionic degrees of freedom (both spin and orbital) according to the following convention: (1) all modes with i = 0 are "uncorrelated", i.e., they appear only in the quadratic part ofĤ (as the p modes of the ALM, introduced in the main text); (2) all modes with i ≥ 1 are "correlated" (as the d modes of the ALM), i.e.,Ĥ loc Ri are generic operators (including interaction terms) constructed from c † Riα c Rjβ . As in the multi-orbital GA, we assume that: and thatĤ loc Ri includes both the two-body local terms of the Hamiltonian and the one-body local term represented as:Ĥ Note that in first-principle calculations of real materials based on density functional theory in combination with many-body techniques [2][3][4][5], such as the GA or DMFT, the solution is obtained by solving recursively multi-orbital Hamiltonian with the structure of Eq. (1) [4]. It is important to note that, within this context, the correlated degrees of freedom correspond to localized orbitals [6,7]. Therefore, as pointed out in the conclusions of the main text, the hopping parameters between correlated and uncorrelated modes is much larger with respect to the hopping between correlated modes (as in the ALM studied in the main text). Therefore, the pathology of the GA method uncovered in this work and our observation that the g-GA extension resolves this problem are very relevant in ab-initio calculations of correlated materials. arXiv:2106.05985v2 [cond-mat.str-el] 24 Jul 2021 Figure 1. Schematic representation of the g-GA variational ansatz. The wavefunction |ΨG is constructed by mapping a generic single-particle wavefunction |Ψ0 constructed in an auxiliary Hilbert space, using the operatorPG. Both |Ψ0 andPG are determined variationally The g-GA variational ansatz The g-GA consists in minimizing the expectation value ofĤ with respect to a wave function represented as follows: where |Ψ 0 is a single-particle wavefunction constructed in an auxiliary Hilbert space, generated byν i > ν i degrees of freedom f † Ria (a = 1, ..,ν i ) for each (R, i), while: is an operator mapping the local auxiliary-space states into the physical space, see Fig. 1, where q i (j) represents the i-th digit of j in binary representation, and the matrix Λ i controls howP Ri modifies the weight of the local electronic configurations. The key reason why the g-GA variational ansatz generalizes the standard multi-orbital GA theory is thatν i > ν i . In the present work we focus on the normal phase. Note that, to enforce the condition that |Ψ G is an eigenstate of the number operator, we do not need to assume thatP Ri commutes with the number operator, but only that: where m i is integer. Our choice of m i (that, in principle, could be considered as an arbitrary variational parameter) will be specified below. As in the standard multi-orbital GA, the variational wave function is restricted by the following conditions: which are commonly called "Gutzwiller constraints". Furthermore, the so-called "Gutzwiller Approximation", which becomes exact in the limit of infinite coordination number [8,9] -where Dynamical Mean Field Theory (DMFT) is exact [10]-is assumed.

Matrix of slave-boson amplitudes
Following Refs. [4,11,12], we introduce the so-called matrix of slave-boson (SB) amplitudes [13][14][15]: By substituting Eq. (25) in Eqs. (18)- (21), it can be readily verified that the Gutzwiller constraints can be rewritten as follows: and that Eq. (15) can be rewritten as follows: Embedding mapping Following Ref. [4], we introduce the so called "embedding states," which are related to the SB amplitudes as follows: where U PH is a particle-hole transformation acting over the |n; i states and Note that the set of all embedding states represented as in Eq. (29) constitute a Fock space, corresponding to an "impurity" (generated by the Fermionic degrees of freedomĉ iα , α ∈ {1, .., ν i }) and a "bath" (generated by the Fermionic degrees of freedomf ia , a ∈ {1, ..,ν i }).
The case B = 1 corresponds to the standard GA theory, where the bath has the same size of the impurity. In the present work we assumed thatν i = Bν i and set B = 3, see Fig. 2. Furthermore, we set m i = (ν i − ν i )/2 = 2 (see Eq. (9)). Note that, from the definitions [Eqs. (25), (29)], it follows that this is equivalent to assume that |Φ i has a total of (ν i + ν i )/2 =4 electrons (i.e., that the embedding states are half-filled).
It can be readily verified by inspection that the Gutzwiller constraints can be rewritten as follows: the expectation value of the local terms ofĤ can be calculated as: and that Eq. (15) can be rewritten as: In summary, by substituting Eqs. (13) and (14) in Eq. (12) and using the equations above, we deduce that the goal is to minimize the following expression for the variational energy: while fulfilling the Gutzwiller constraints [Eqs. (33), (34)].
Lagrange formulation of the g-GA Following Refs. 1, 4, and 15, the constrained energy-minimization problem described above can be cast in terms of the following Lagrange function: where N is the total number of unit cells, E and E c are real numbers, ∆ i , λ c i and λ i areν i ×ν i Hermitian matrices, D i and R i are generic (rectangular)ν i × ν i matrices. The auxiliary HamiltoniansĤ qp andĤ emb , which are called "quasiparticle Hamiltonian" and "Embedding Hamiltonian," respectively, are defined as follows: The saddle-point of the Lagrangian L defined in Eq. (38) is given by the following equations: where f is the zero-temperature Fermi function and in Eqs. (41), (42) we introduced the following block-matrices: where [1] n×n is the n × n identity matrix and [0] n×n is the n × n zero matrix. We also introduced the following projectors: Finally, we introduced the following expansions of the matrices ∆ i , λ i , λ c i , in terms of an orthonormal basis of Hermitian matrices {[h i ] s } (with respect to the canonical scalar product (A, B) = Tr A † B ): where d 0 i s , [l i ] s and [l c i ] s are real-valued coefficients.
Algorithmic structure of the g-GA The equations (41)- (46) can be solved numerically using the following numerical procedure (represented schematically in Fig. 2 where, as expected, Σ(ω) is local (k-independent) and Σ i (ω) are the ν i × ν i matrices ∀ i ≥ 1, given by the following equation: which reduces to Eq. 10 of the main text for single-orbital impurities. Note that, as pointed out in the main text, from Eq. (64) it follows that the poles of the Green's function lie on top of the eigenvalues of the quasi-particle Hamiltonian, see Eq. (39). However, because of the matrices R † i and R j at the left and the right of Eq. (64) -which are rectangular, -only a portion of the quasi-particle spectral weight is physical (see Fig. 3 of the main text).

PATHOLOGY OF THE BARE GA IN THE NARROW-BANDWIDTH LIMIT OF THE ALM
As pointed out in the main text, within the GA the Mott transition of the ALM occurs (by construction) only when d † iσ p iσ GA = r Ψ 0 | f † iσ p iσ |Ψ 0 = 0, i.e., when the variational parameter r vanishes. However, as shown in Fig. 1 of the main text, the GA critical value U c is overestimated dramatically for the ALM, especially for p −1. To explain this behavior, it is insightful to inspect the GA solution also in the narrow-bandwidth limit, i.e., setting the half-bandwidth W of the Bethe-lattice hopping matrix to 0 -corresponding to a series of decoupled p-d dimers.
The resulting evolution of U c (i.e., the interaction strength such that r vanishes) is shown here in Fig. 3, where we set V = 1 and consider p −1. We observe that, not surprisingly (since V = 1), d † iσ p iσ GA does not vanish at U ∼ 0 (neither for W = 0 nor for W = 1). In fact, at U = 0 the GA wavefunction is exact. On the other hand, in this limit any infinitesimal U should induce a Mott transition for W → 0. This simple observation showcases very clearly how, as pointed out in the main text, the p-d charge fluctuations cannot be neglected in the Mott phase -as they must coexist in the narrowbandwidth limit of the ALM. In other words, the argument above shows that -for small W -the GA behavior of