Variational exact diagonalization method for Anderson impurity models
aa r X i v : . [ c ond - m a t . s t r- e l ] J un Variational exact diagonalization method for Anderson impurity models
M. Sch¨uler,
1, 2, ∗ C. Renk,
1, 2 and T. O. Wehling
1, 2 Institut f¨ur Theoretische Physik, Universit¨at Bremen, Otto-Hahn-Allee 1, 28359 Bremen, Germany Bremen Center for Computational Materials Science,Universit¨at Bremen, Am Fallturm 1a, 28359 Bremen, Germany (Dated: June 3, 2015)We describe a variational approach to solving Anderson impurity models by means of exact di-agonalization. Optimized parameters of a discretized auxiliary model are obtained on the basisof the Peierls-Feynman-Bogoliubov principle. Thereby, the variational approach resolves ambigui-ties related with the bath discretization, which is generally necessary to make Anderson impuritymodels tractable by exact diagonalization. The choice of variational degrees of freedom made hereallows systematic improvements of total energies over mean field decouplings like Hartree-Fock.Furthermore, our approach allows us to embed arbitrary bath discretization schemes in total energycalculations and to systematically optimize and improve on traditional routes to the discretizationproblem such as fitting of hybridization functions on Matsubara frequencies. Benchmarks in termsof a single orbital Anderson model demonstrate that the variational exact diagonalization methodaccurately reproduces free energies as well as several single- and two-particle observables obtainedfrom an exact solution. Finally, we demonstrate the applicability of the variational exact diagonal-ization approach to realistic five orbital problems with the example system of Co impurities in bulkCu and compare to continuous-time Monte Carlo calculations. The accuracy of established bathdiscretization schemes is assessed in the framework of the variational approach introduced here.
PACS numbers: 72.80.Rj; 73.20.Hb; 73.61.Wp
I. INTRODUCTION
The Anderson impurity model (AIM) is a generalmodel for the description of interacting impurities inmetallic host systems. Originally, it was developed to de-scribe single atoms with open d- or f-shells embedded inbulk materials and to understand the formation of theirmagnetic moments . Furthermore the model includes thewidely discussed Kondo physics . Multi orbital variantsof the AIM gained considerable attention in the contextof rare-earth impurity systems as well as more recentlymagnetic adatoms or molecules on surfaces . Finally,dynamical mean field theory (DMFT) links correlatedbulk systems as well as nanostructures to Anderson im-purity models.To address the electronic structure of realistic corre-lated electron materials one often resorts to LDA++approaches , where quantum lattice or impurity mod-els are derived from first principles calculations. The re-sulting models are typically multi-orbital models includ-ing complex hybridization between the impurity and acontinuous bath of states from the surrounding material,which brings along two challenges: First, the numericalsolution of the impurity models and second the interpre-tation of the physics contained in these generally complexmodels in more simple terms. Experiments are for in-stance often interpreted in terms of atomic spins, crystalfield, ligand field or cluster approaches , which typi-cally involve a small discrete set of bath states or no bathstates at all. The link of the complex, ab initio derivedmodels and simpler phenomenological models is a prioriunclear and relates to the so-called bath discretizationproblem of exact diagonalization solvers of the AIM. The solution of the Anderson impurity model for gen-eral parameters has to be done numerically by means ofe.g. quantum Monte Carlo (QMC), numerical renor-malization group (NRG), or exact diagonalization (ED)methods. While NRG and QMC are in principle nu-merically exact methods, they become computationallyvery demanding, when dealing with many orbitals, hy-bridization functions with low symmetry, spin orbit cou-pling and general fermionic four operator Coulomb ver-tices. ED methods deal with low symmetries and generalCoulomb vertices at no additional computational costbut suffer from the so-called bath discretization prob-lem: Due to the exponential growth of the many particleFock space with the system size, it can handle only afew bath levels per orbital. A mapping of the continu-ous bath to a discrete version has to be found. Severalapproaches to this task have been introduced. One is tofit the hybridization function of the continuous bath onMatsubara frequencies , another is to represent the hy-bridization function by a continued fraction and to linkits coefficients to the parameters of the bath . Theseschemes are systematic in the sense that they convergeto the full model when including more and more bathsites. However, in the multi-orbital case, the number ofbath sites is limited (typically on the order of three orless for a five orbital impurity problem), so the quality ofthe mapping can hardly be checked by an analysis of theconvergence.Basically two different strategies have been laid out tocircumvent this problem. First, the many body Hilbertspace can be truncated in the sense of configuration in-teraction (CI) expansions, which have a long traditionin the context of quantum impurity problems and aresubject of recent developments. CI expansions arevariational, i.e. they deliver upper bounds for total ener-gies, but they do not provide simplified auxiliary Hamil-tonians. On the other hand, there are several approachestowards optimized cluster approximations to Andersonimpurity problems. In this context, self-energy functionaltheory is based on an extremal principle but it is notvariational regarding total energies and does not allowfor variations of interaction parameters or the interactingorbitals. More general optimizations are possible in theframework of the so-called self-energy embedding theory(SEET) , which is however not variational.In this paper, we combine ideas of variational ap-proaches and optimized cluster approximations to theAIM. We introduce a strictly variational method of ap-proximating an AIM with continuous bath by an AIMwith finite strongly reduced number of bath sites, whichwe call variational ED method. It guaranties an optimalapproximation to the AIM for a given number of bathsites in the sense of thermodynamic ground state proper-ties. The method is based on the well-known Peierls-Feynman-Bogoliubov variational principle , whichfinds optimal effective models on the basis of an optimaldensity matrix by minimizing a free energy functional.We will introduce the AIM and the variational prin-ciple in Sec. II, where we also explain the details ofcalculating the Peierls-Feynmann-Bogoluibov free energyfunctional and how to minimize it efficiently. By treatinga single orbital model with the variational ED methodin Sec. III we analyze its performance in comparisonto an exact treatment, established bath discretizationmethods as well as Hartree-Fock theory. In Sec. IV,we demonstrate the applicability of the method to realis-tic five orbital system with the example of Co impuritiesin bulk Cu and compare to QMC simulations. We showthat the variational ED method leads to systematicallylower, i.e. more accurate, free energy estimates than un-restricted Hartree-Fock and traditional bath discretiza-tion schemes also in the multiorbtial case. Conclusionsand outlook are given in Sec. VI. II. MODEL AND METHOD
After introducing the Anderson impurity model we willrecapitulate the Peierls-Feynman-Bogoliubov variationalprinciple and show how to apply it to discretize Andersonimpurity models in an optimal manner.
A. The Anderson impurity model
The Hamiltonian of the initial AIM (termed “originalmodel” hereafter) reads H = H bath + H hyb + H imp . (1) The bath is described by H bath = X αk,σ ε αk n cαkσ , (2)where ε αk is the energy of the bath state withband/orbital index α and some additional quantum num-ber k . n cαkσ = c † αkσ c αkσ is the corresponding particlenumber operator. The hybridization part H hyb = X αk,σ V αk (cid:16) c † αkσ d ασ + d † ασ c αkσ (cid:17) (3)couples the bath sites of one band to an orbital of theimpurity with a coupling strength V αk . The bath elec-trons with spin σ are created and annihilated by c † αkσ and c αkσ , respectively, while d † ασ ( d ασ ) denote the cre-ation (annihilation) operators of the impurity electrons.The impurity site is described by H imp = X α,σ ε dα n dασ + X α,β,γ,δ,σσ ′ U αβγδ d † ασ d † βσ ′ d γσ ′ d δσ , (4)which contains the on-site Coulomb interaction U αβγδ and the on-site energies ε dα . By integrating out all bathdegrees of freedom we arrive at the hybridization function∆ α ( ω ) = X k V ∗ αk V αk ω + i + − ε αk , (5)which describes the energy dependent coupling of the im-purity to the bath. B. Peierls-Feynman-Bogoliubov variationalprinciple
Given a Hamiltonian H , which is “difficult” to solve,we search for an optimal approximation to H within aset of simpler effective Hamiltonians ˜ H . The Peierls-Feynman-Bogoliubov variational principle providesus with a prescription on how to fix the parameters of ˜ H in a thermodynamically optimal way, i.e. such that thecanonical density matrix resulting form ˜ H approximatesthe density matrix corresponding to H as close as possi-ble. More strictly speaking: the canonical density opera-tor ρ ˜ H = 1 /Z ˜ H exp( − β ˜ H ) of the auxiliary system, where Z ˜ H = Tr exp( − β ˜ H ) is the partition function, approxi-mates the exact density operator ρ derived from H asclose as possible, when the Peierls-Bogoliubov-Feynmanfunctional ˜Φ[ ρ ˜ H ] = Φ ˜ H + h H − ˜ H i ˜ H , (6)becomes minimal. Here Φ ˜ H = − β ln Z ˜ H is the free en-ergy of the effective system. h H − ˜ H i ˜ H = Tr ρ ˜ H ( H − ˜ H )denotes a thermodynamic expectation value with respect Figure 1. (Color online) Illustration of the original and effec-tive model for the case of one orbital and six bath sites. Bluerepresents bath character and red impurity character: In theeffective model bath and impurity states can be mixed. ε unc n are eigenvalues of h R kk ′ . to the effective system. In the case of ρ ˜ H = ρ the func-tional ˜Φ[ ρ ˜ H ] becomes minimal and coincides with the freeenergy Φ H of the original system. In our case H repre-sents the full AIM, Eqs. (1)-(4), and ˜ H is the model withdiscretized bath, which is now introduced. C. Effective Hamiltonian
The structure of the effective Hamiltonian for the caseof a single impurity orbital is depicted in the right panelof Fig. 1. In contrast to the original model (left panelof Fig. 1), the effective model consists of two decoupledparts: First, the effective impurity coupled to one bathsite only and second the remaining bath sites. I.e. wepartition the full Hilbert space H into a correlated sub-space C (first part) and an uncorrelated rest R (secondpart). In this work, we consider for concreteness a clusterconsisting of a multi-orbital impurity and one bath siteper impurity orbital for the correlated space but otherchoices are similarly possible. The single particle statesof the effective model are related to those of the origi-nal model by a unitary transformation, which allows formixing of original “bath” and “impurity” character inthe effective model.The optimal matrix elements of the effective model, aswell as the optimal unitary transformation are found byminimizing the functional (6). The states spanning C aredefined by | ˜ d α i = u d α d α | d α i + X k u d α c αk | c αk i , (7) | ˜ c α i = u c α d α | d α i + X k u c α c αk | c αk i , (8)where the coefficients u are chosen such that | ˜ d α i and | ˜ c α i form an orthonormal basis of C . An orthonormal basis spanning R is defined by | ˜ c αk i = u c αk d α | d α i + X k ′ u c αk c αk ′ | c αk ′ i , k > . (9)As a whole, the coefficients u form a unitary matrix. Inpractice, we obtain the elements of this matrix from theQR decomposition of a matrix, in which the first two rowsare defined by the coefficients of | ˜ d α i and | ˜ c α i and allother elements are zero. This leads to a new orthonormalbasis for the full space H , which provides the partitioningaccording to H = C ⊕ R . The ansatz for the effectiveHamiltonian in this new basis explicitly reads˜ H = ˜ H C + ˜ H R , (10)with˜ H C = X α ˜ V α (cid:16) ˜ c † α ˜ d α + ˜ d † α ˜ c α (cid:17) + X α ˜ ε α ˜ c † α ˜ c α + X α ˜ ε dα ˜ d † α ˜ d α + X αβγδ,σσ ′ ˜ U αβγδ ˜ d † ασ ˜ d † βσ ′ ˜ d γσ ′ ˜ d δσ (11)and ˜ H R = X α, ( k,k ′ ) > h R αkk ′ ˜ c † αk ˜ c αk ′ . (12)It is stressed, that the new states are linear combi-nations of the original impurity and bath states, lead-ing to mixed basis states. The new impurity states canhave some amount of bath character and vice versa. TheHamiltonian in Eq. (11) states a many-body problemwhich can be solved by exact diagonalization, as long asits Hilbert space is sufficiently small. In contrast, theHamiltonian (12) states a one-particle problem and canbe solved by diagonalizing the matrix h R αkk ′ . In sum-mary, the Hamiltonian ˜ H = ˜ H C + ˜ H R defines an effectiveHamiltonian, which can be solved exactly and thus thefunctional (6) can be calculated.This ansatz implies several approximations. First, allcouplings between C and R are neglected. Second inter-action terms are restricted to new effective impurity or-bitals ˜ d α within C , which is motivated by the fact that theoriginal model includes only on-site interactions too. Thelatter approximation can be relaxed to include arbitraryinteractions within C , but we keep it here for simplicity.Finally, we note that the amount of variational degreesof freedom in the variational ED approach is such thatit includes Hartree-Fock as the limiting case ˜ U αβγδ → D. Implementation
In order to perform the minimization in practice, thenumber of free parameters has to be kept sufficiently low.First, for the rest of this work it is assumed that theCoulomb tensor ˜ U αβγδ is not varied. We choose it to bethe same as in the original model. Test calculations haveshown, that the variation of the Coulomb tensor is notcrucial, as this can mostly be absorbed into the variationof the impurity level. The single particle matrix elementsof ˜ H C are assumed to be free parameters. In principle,the parameters of the uncorrelated Hamiltonian are freeparameters, too. However, to further reduce the numberof free parameters, we define h R αkk ′ by a projection ofa Hartree-Fock solution of the original Hamiltonian ontothe states | ˜ c αk i . The Hartree-Fock solution of the originalHamiltonian (1) can be written as H HF = X n ε HF n c † n c n , (13)where the eigenstates | n i and energies ε HF n are found byapplying the Hartree-Fock decoupling d † ασ d † βσ ′ d γσ ′ d δσ →h d † ασ d δσ i d † βσ ′ d γσ ′ + h d † βσ ′ d γσ ′ i d † ασ d δσ −h d † ασ d γσ ′ i d † βσ ′ d δσ − h d † βσ ′ d δσ i d † ασ d γσ ′ (14)to (4) and solving the resulting non-interacting prob-lem self-consistently. The single particle matrix elementswithin the uncorrelated space R explicitly read h R αkk ′ = X n ε HF n h ˜ c αk | n ih n | ˜ c αk ′ i . (15)In order to not break any spin rotation symmetries, re-stricted Hartree-Fock is used.The functional ˜Φ[ ρ ˜ H ] now depends on the unitarytransformation and on the matrix elements of ˜ H C . Theminimum of the functional is searched by iterative meth-ods. Thus, the functional ˜Φ[ ρ ˜ H ] has to be calculated forvarious points of the variational space with the computa-tionally most expensive part being here the diagonaliza-tions of ˜ H C . Therefore, we first search for fixed parame-ters in ˜ H C a corresponding optimal unitary transforma-tion matrix defining the optimal partitioning H = C ⊕ R using an SLSQP algorithm . The search of the mini-mum w.r.t. the parameters of ˜ H C is then done by theNelder-Mead algorithm . The number of independentparameters can be further reduced when the original sys-tem shows symmetries like orbital degeneracies which areassumed not to be broken in the effective model. III. BENCHMARK FOR A SINGLE ORBITALAIM
In this section the variational ED method is tested forits performance in reproducing the density operator aswell as observables such as the occupation number, dou-ble occupancy and crystal orbital overlap populations ofa simple original model. The original model we considerhere is a single orbital model with only 6 bath sites, which itself can be solved by exact diagonalization. The de-tailed setup of the model is as follows: The impurity levelis ε d = − . U = 4 . ε b in an interval of 2 eV (i.e. the bandwidth of thebath). The coupling is V k = 0 . ε b is swept from − . . ε F = 0. The systemis solved for T = 0. The model is first solved exactly, sec-ond by the variational ED method, third by unrestrictedHartree-Fock and finally by ED using reduced bath sitesobtained by fitting of the hybridization functions on theimaginary Matsubara frequency axis . The latter typeof approaches require generally the introduction of a so-called weight function W n for the fitting procedure, asexplained in the appendix (A).The central object for the assessment of the qualityof the methods is the difference between ˜Φ[ ρ ˜ H ] and theexact free energy Φ H , as shown in Fig. 2a). For bathsites energetically far away from the Fermi level and fromsingle particle excitation energies of the impurity ( | ε b | > | ε b | < | ε b | > | ε b | < | ε b | < W n = 1 /ω n leads to the lowest free energies.Otherwise, the constant weight function W n = 1 showssmallest deviations of the free energy functional from theexact solution. The variational ED method is generallyvery close to the exact solution. Only for the specialcase of a strictly symmetric distribution of the bath sitesaround the Fermi energy ( ε b = 0 eV) a deviation on theorder of meV occurs.Fig. 2b)-e) shows a comparison of several observables(chemical bond strength b), bath occupation c), impurityoccupation d) and double occupation e)) calculated withthe different methods. For the outer most regions ( | ε b | > ε b ≈ | ε b | < −6 −4 −2 0 2 4 6 ε b (eV)0.000.050.100.150.200.250.300.35 ˜ Φ − Φ H ( e V ) a) varfit0fit1HF −2 210 -5 -3 -1 −3 −2 −1 0 1 ε k (eV) −0.20.00.20.40.60.81.0 u i j ε b =−0.3eV f) u dc k u dd u c c k u c d −6 −4 −2 0 2 4 6 ε b (eV) (cid:0) i (cid:2) n i (cid:1) c) origvarHFfit0fit1 −6 −4 −2 0 2 4 6 ε b (eV) (cid:2) n d (cid:1) d) origvarHFfit0fit1−6 −4 −2 0 2 4 6 ε b (eV) (cid:2) n d ↑ n d ↓ (cid:1) e) origvarHFfit0fit1 −6 −4 −2 0 2 4 6 ε b (eV) (cid:0) i | (cid:2) c † i d (cid:1) | b) origvarHFfit0fit1 Figure 2. (Color online) Benchmark of different ED ap-proaches and spin-polarized Hartree-Fock theory against anexact solution for single orbital Anderson impurity modelswith a mean bath energy ε b . (a) Difference between the freeenergy functional obtained by different approximate meth-ods according to Eq. (6) and the free energy of the originalmodel. The inset shows a close up view for the vicinity ofthe Fermi energy on a logarithmic scale. (b)-(e): Compari-son of local and non-local observables obtained from an exactsolution (“orig”, bold cyan) and calculated by the four dif-ferent approximate methods, i.e. the variational ED method(“var”, solid black), Hartree-Fock (“HF”, dashed red) and fitsof hybridization functions on the imaginary axis with differentweight functions ( W n = 1: “fit0”, dotted blue; W n = 1 /ω n :“fit1”, dashed blue.) Panel b) shows the chemical bondstrength, c) shows the total bath occupation, d) the impu-rity occupation and e) the double occupation. Panel e) showsthe coefficients of the unitary transformation linking the orig-inal model to the optimized effective model with one bath-siteper spin-orbital in C for the example of ε b = − . states spanning the correlated space | ˜ d i and | ˜ c i (see Eqs.(7) and (8)) for an original model with the bath centeredaround the energy ε b = − . | d i character with small bath admixtureand can approximately be interpreted as the old impuritystate. The coupled effective bath state is nearly a purelinear combination of old bath states, where states closerto the Fermi energy contribute stronger than those fur-ther away. This behavior is very reminiscent of effectivebath wave functions obtained in variational approacheslike the Varma-Yafet or the Gunnarsson-Sch¨onhammerexpansion. For the treatment of original models with far more bathsites, it is important to note, that the coefficients definingthe unitary transformation from the original bath states to the effective impurity and bath orbitals, i.e. u c c k and u dc k , vary smoothly as function of the bath energies oneither side of the Fermi energy. IV. CO IMPURITIES IN CU: APPLICATIONTO A REALISTIC FIVE ORBITAL SYSTEMA. The AIM derived from LDA
In this section the variational ED method is applied toa realistic model of Co impurities in bulk Cu, which hasbeen obtained from super-cell DFT calculations and hasbeen analyzed using a QMC impurity solver in Ref. 7.The cubic symmetry of the Cu crystal leads to a split-ting of the Co 3 d -orbitals into blocks of t and e g sym-metry. From the DFT hybridization function, which isa continuous function, we obtain our initial model as-suming some large number of bath sites, here 100 perorbital. (This number does not present a limiting fac-tor and could be chosen arbitrarily larger). The bathsites are assumed to be equidistantly distributed between −
10 eV and 10 eV, and the hybridization terms V αk arethen found by fitting the imaginary part of a discretizedhybridization function∆ disc ( ω ) = X k V ∗ αk V αk ω − ε k + iδ , (16)with some broadening δ = 0 . ω ) on the real axis. The V ik areplotted in Figure 3a). The crystal field obtained from theDFT calculation is ε d e g − ε d t = 0 .
136 eV. As in Ref. 7,we consider a rotationally invariant Coulomb interactiondefined by U αβγδ = l X k =0 a k ( α m β m , γ m δ m ) F k , (17)where a k ( α m β m , γ m δ m ) are the Gaunt coefficients and where F = U , F = 14 / (1 + 0 . J and F =0 . F are Slater parameters with the average Coulombinteraction U = 4 . J = 0 . µ = 27 eV as in Ref. 7.All data is obtained at a inverse temperature of β = 40,like in the case of the QMC simulations. Finally, we as-sume that the cubic symmetry of the system prevails,which means that only two independent sets of matrixelements (for the t and e g states) have to be variedduring the minimization of ˜Φ[ ρ ˜ H ]. B. Implementation of the variational ED methodfor the 5 orbital AIM
We compare two different sets of variational degreesof freedom for the optimization of the one particle basis,which we refer to as “bath” and “all”. In the “bath” case,only bath sites are optimized, i.e. we fix the expansioncoefficients u c α d α = 0, u d α c αk = 0 and u d α d α = 1. This leadsto considerably less variational parameters and a muchsmaller amount of expectation values to be calculated ineach step of the iteration. In the second approach, “all”,which is computationally more demanding because thefull two-particle density matrix of the effective systemhas to be calculated, we optimize the full one particlebasis of the bath and that of the impurity.Because a full optimization of the parameters of theeffective model is computationally challenging, it is cru-cial to start the optimization from a good initial guess.We obtain such initial guesses for the parameters of thebath by fitting of hybridization functions on Matsubarafrequencies as introduced in the appendix A. We choose˜ ε dα = ε dα as the initial guess for the parameters of the im-purity. The resulting first guesses using different weightfunctions are summarized in the Table I. While all weightfunctions lead to setups with the t bath sites above theFermi energy and the e g bath sites below, the details oftheir energetic positions and the hybridization strengthsdepend strongly on the form of W n . Adding more weighton features on small Matsubara frequencies shifts the ef-fective bath parameters to smaller values. The qualityof these starting guesses in the context of the variationalprinciple is discussed in the next section. Table I. Parameters of the effective model (see Eq. (11))obtained by the fit of hybridization functions on imagi-nary frequencies using different weight functions ( W n =1 , /ω n , /ω n , see appendix A) and the iterative optimization(“var”).weight function 1 1 /ω n /ω n var˜ ε t (eV) 3 .
203 0 .
775 0 .
068 2 . V t (eV) 1 .
563 0 .
606 0 .
223 1 . ε e g (eV) − . − . − . − . V e g (eV) 1 .
049 0 .
170 0 .
156 1 . ε d t (eV) − . − . − . − . ε d e g (eV) − . − . − . − . The large number of bath sites (100 per orbital) in theoriginal model leads to 202 variational parameters defin-ing the unitary transformation in the “all” case or 100parameters in the “bath” case for each orbital. The ob-servation, that u d α c αk , u c α c αk are smooth functions of energy(c.f. Fig. 2e) below and above E F , leads to the possibil-ity of expanding them in a set of smooth functions andthereby reducing the number of variational parametersconsiderably. Here, we chose five Chebyshev polynomials T n ( k ) per orbital for bath sites above and five for those Table II. The free energy functional ˜Φ, the total impurity oc-cupancy n d and the local spin S as obtained from simulationsof the AIM for Co impurities in Cu. The values of ˜Φ areshown as differences to the results from unrestricted Hartree-Fock (UHF): ∆ ˜Φ = ˜Φ − ˜Φ UHF . Total impurity occupation andspin calculated with the variational ED method are comparedto QMC solutions of the AIM from Ref. 7. Different flavors ofthe variational ED method are considered: first “bath” andsecond “all” with the model parameters obtained from thefits of the hybridization function on the imaginary frequenciesusing different weight functions ( W n = 1 “fit0”, W n = 1 /ω n “fit1” and W n = 1 /ω n “fit2”) and finally full optimizationof transformation and model parameters labeled “var”. Re-stricted Hartree-Fock (RHF) and unrestricted HF (UHF) re-sults are also shown.∆ ˜Φ(eV) h n d i S QMC 7.78 ± ± below the Fermi energy. Therefore only 22 (or 10 in thecase of “bath”) parameters per orbital have to be variedto find the optimal unitary transformation to embed theeffective model into the full Hilbert space. C. Results
We will first compare free energy estimates as well asdifferent local observables obtained from variational EDtreatments to unrestricted Hartree-Fock as well as QMCcalculations. Afterwards, we investigate the nature of theoptimized effective bath and impurity states as obtainedfrom the variational ED treatment.
1. Free energy functional and local observables
Table II shows the free energy functional ˜Φ (relative tounrestricted Hartree-Fock (UHF)), as obtained with dif-ferent starting points and different amounts of variationaldegrees of freedom in the variational ED approach. Inthe case of “bath”, the constant weight function (“fit0”, W n = 1) leads to the lowest values of the ˜Φ[ ρ ]. The mod-els derived using weight functions W n = 1 /ω n (“fit1”)and W n = 1 /ω n (“fit2”) lead to free energy estimateswhich are about 1 to 3 eV higher in energy. The situ-ation for the case of “all” is similar. On this basis, wehave chosen the starting guess obtained with the constantweight function for the full optimization of the effectivemodel parameters. The resulting parameters are shownin the last column of Tab. I and are close to the start-ing guess. The full optimization schemes (“var,bath” and“var,all”) find parameters which lower the functional ˜Φconsiderably for “all” and slightly for “bath”.Regarding the impurity occupation ( h n d i , see Tab. II),we see that the description by unrestricted Hartree-Fockis rather close to QMC, whereas restricted Hartree-Fockoverestimates the occupation. All versions of exact di-agonalization lead to occupations close to the QMC re-sults, and many cases within the QMC error bars. Thespin S (defined as h ˆ S i = S ( S + 1)), which is a two-particle observable, reveals the problems of the Hartree-Fock description. S is vastly overestimated by unre-stricted Hartree-Fock. The variational ED methods, es-pecially the “all” case for the constant weight function(“fit0”) and the full optimized ED model, lead to resultsclose to QMC.To compare the results of the variational ED methodwith those ED methods based on fitting of the hybridiza-tion function on the imaginary axis, we should comparethe “fit0,bath”,“fit1,bath”, and “fit2,bath” cases to thecorresponding “all” and “var,all” cases. We see that hav-ing more variational degrees of freedom leads to an im-proved description of the free energies, as it should be.In general, we learn that only in the case of optimizingboth effective bath and effective impurity states (termed“all”) we reach lower values of the free energy functional˜Φ than with unrestricted Hartree-Fock: The freedom toform mixtures of bath and impurity states in the effec-tive model is important to describe the free energy andlocal observables of the system adequately. As the vari-ational ED method provides more accurate (free) energyestimates than unrestricted Hartree-Fock, the approachintroduced here could be a way to improve LDA+U totalenergy schemes.
2. The effective basis states
Now we analyze the unitary transformation relatingthe optimized basis states of the effective model and theoriginal basis states. The transformation obtained for themodels from the starting guesses with weight functions W n = 1, 1 /ω n , and 1 /ω n is shown in Fig. 3 b), c),and d), respectively. We observe a clear trend that theadmixture of the original bath states into the effectivebath state ( u c αk c α , blue lines) is strongest in the vicinity ofthe effective bath site energy ˜ ε α . For bath states closeto the Fermi energy we get a sharp cut off for states onthe opposite side of the Fermi energy. This is very similarto first order configuration interaction treatments of theAIM . The original impurity admixture in the effectiveimpurity ( u dd α ) rises with the distance of the effectivebath site from the Fermi energy. − − ε k (eV) . . . . V α k a) t2geg − − ε k (eV) − . . . . . . . . u c α c α k / u d α c α k b) W n = 1 u c d t2g = 0 . u c d eg = − . u dd t2g = 0 . u dd eg = 0 . ε e g ˜ ε t c t2gc egd t2gd eg − − ε k (eV) − . . . . . . . . u c α c α k / u d α c α k c) W n = 1 /ω n ˜ ε e g ˜ ε t u c d t2g = 0 . u c d eg = − . u dd t2g = 0 . u dd eg = 0 . c t2gc egd t2gd eg − − ε k (eV) − . . . . . . . . u c α c α k / u d α c α k d) W n = 1 /ω n ˜ ε e g ˜ ε t u c d t2g = 0 . u c d eg = − . u dd t2g = 0 . u dd eg = 0 . c t2gc egd t2gd eg Figure 3. (Color online) a) Hopping matrix elements betweenthe impurity orbitals and the bath from the original AIM forCo impurities in Cu (solid t , dashed e g ). b)-d) Coefficientsdefining the optimal transformation from original bath statesto effective bath states (blue/dark gray) and effective impu-rity states (red/light gray), c.f. Eqs. (7) and (8). Optimizedtransformations for different effective models defined throughfits of the hybridization with weight functions W n = 1 b), W n = 1 /ω n c), and W n = 1 /ω n d) are shown. The ener-gies of the effective coupled bath sites ˜ ε α are depicted as thinvertical lines. The numerical values of the transformation co-efficients defining the admixture of original impurity states tothe effective bath and impurity states are given as insets. V. SPECTRAL FUNCTIONS
The variational principle results in an effective modelwhich represents thermodynamic ground state proper-ties in an optimal manner. This is a necessary but nota sufficient condition to give a good approximation alsofor excitation spectra. In the following, we study theone particle spectral function for single orbital impu-rity benchmark systems from Sec. III with the bathstates centered around ε b = 0 . V k = 0 . V k = 0 . G α ( ω ) = 1 Z X µν (cid:12)(cid:12) h µ | d † α | ν i (cid:12)(cid:12) ω + E ν − E µ − i + (cid:0) e − βE ν − e − βE µ (cid:1) , (18)where in our calculations 0 + is replaced by a broadeningof δ = 0 . β = 3200,which is very close to the T = 0 calculations of expecta-tion values in Sec. III.We assess the quality of the spectra obtained from thevariational ED method, from ED with the hybridizationfunction fitted on the Matsubara axis (with two differ-ent weight functions W n = 1 and W n = 1 /ω n ) and froman unrestricted Hartree-Fock treatment by comparing tothe exact spectrum of the original model. For the caseof W n = 1 /ω n , we additionally compare the spectra fordifferent amounts of variational freedom in choosing thebasis states of the effective model, where we either op-timized the bath states only (termed “bath”, c.f. Sec.IV) or allowed impurity and bath states to mix (termed“all”).The dominant features of the original spectrum in thecase of strong hybridization ( V k = 0 . − . W n = 1)shows a similar picture but with major peaks and satel-lites shifted considerably towards the Fermi energy. Theresult for the weight function emphasizing small Matsub-ara frequencies ( W n = 1 /ω n ) leads to a good represen-tation of the peaks around − . V k = 0 . − . . W n = 1 /ω n , the variational methodgives satisfactory results in both cases. The spectra fromvariational ED are in both cases at least as close tothe spectra of the original model as the best spectrumobtained with any of the two bath fitting procedures W n = 1 or 1 /ω n ) under investigation. − − − − ω ( eV )0 . . . . I n t en s i t y ( a . u . ) a) V k = 0 . origHFvarfit0fit1,“all”fit1,“bath” − − − − ω ( eV )0 . . . . I n t en s i t y ( a . u . ) b) V k = 0 . origHFvarfit0fit1,“all”fit1,“bath” Figure 4. (Color online) One particle impurity spectral func-tions for the single-orbital benchmark model introduced inSec. III with ε b = − . V k = 0 . V k = 0 . W n = 1),dotted blue (“fit1”, W n = 1 /ω n optimizing all states) anddashed blue (“fit1”, W n = 1 /ω n optimizing only bath states). VI. CONCLUSION
In conclusion, we present a variational exact diago-nalization method which provides self-consistently opti-mized parameters of discretized Anderson impurity mod-els. The method is based on an optimal partitioning ofthe system into a correlated part, where electronic inter-actions are explicitly taken into account, and an uncor-related rest and on finding optimal effective Hamiltoni-ans for both parts of the system. To this end, a varia-tion of a free energy functional w.r.t. one-particle basisstates spanning the correlated subspace and the matrixelements of the effective Hamiltonians is performed.A benchmark of the variational ED method against anexact solution of a one orbital Anderson model demon-strates its excellent performance in reproducing groundstate observables of the impurity and bath and addi-tionally a sound performance in reproducing the im-purity spectral function. A comparison with Hartree-Fock and established bath discretization schemes for EDshows that the variational approach introduced here evenworks for difficult cases, i.e. when the bath is symmet-ric around the Fermi energy. Furthermore, applicabilityof the variational ED method to realistic multi-orbitalcases is demonstrated with the example of Co impuritiesin bulk Cu. Also here, the variational method leads to anaccurate description of local one and two particle observ-ables like the impurity occupation and the spin. Energet-ically the method outperforms unrestricted Hartree-Fock,which suggests that the variational ED approach couldbe useful to improve total-energy approaches to corre-lated systems beyond LDA+U. Finally, the method intro-duced, here, can be used to embed established bath dis-cretization schemes such as the fit of hybridization func-tions on Matsubara frequencies into a variational frame-work and to reach an unbiased decision of e.g. whichweight function to choose in the fits of the hybridiza-tion functions. In the example studied, here, the con-stant weight function leads to best results in terms of thefree energy, whereas W n = 1 /ω n leads to results quali-tatively similar to a first order configuration interactionexpansion . For the description of systems closer to theatomic limit, like Co on Cu or 4f systems we expect CIexpansions to converge faster and weight functions like W n = 1 /ω n or W n = 1 /ω n might be a better choice. Thepresented variational approach will allow for an unbiaseddecision in any case.Additionally, the variational method is universal anddifferent choices of the effective correlated space are pos-sible. To gain yet higher accuracy, more bath levels could be included. On the other hand, to make contact withligand field theory or crystal field theory descriptions ofmagnetic impurity systems and nanostructures one couldconsider correlated subspaces with very few or withoutany effective bath orbitals at all. VII. ACKNOWLEDGMENTS
The authors thank M. Katsnelson, A. Lichtenstein, M.Potthoff, P. Bl¨ochl, R. Schade ans S. Barthel for usefuldiscussions as well as the Central Research DevelopmentFund of the University of Bremen and the DFG via FOR1346 for financial support.
Appendix A: Fit of hybridization functions on theMatsubara axis
The method of fitting hybridization functions is shortlyintroduced for the sake of completeness. The basic ideais to minimize a cost function for the inverse impurityGreen’s function (or equivalently the hybridization func-tion) of the discretized model and that of the originalmodel, both defined on the imaginary frequency axis .In the case of one effective bath site, the discrete impurityGreen’s function is defined as g ( iω n ) = iω n − ǫ d − µ − ˜ V iω n − ˜ ε ! − (A1)and the Green’s function of the original model as G ( iω n ) = iω n − ǫ d − µ − X k V k iω n − ε k ! − . (A2)The cost function then reads χ = 1 n max + 1 n max X n =0 W n (cid:12)(cid:12) G − ( iω n ) − g − ( iω n ) (cid:12)(cid:12) , (A3)where W n is a weight function. Popular choices forthe weight function are W n = 1 , W n = 1 /ω n and W n = 1 /ω n . Different weight functions put different em-phasis of low/higher Matsubara frequencies . Through-out this work, we have chosen β = 40 and n max = 1000.This method only provides the effective parameters ˜ ε and ˜ V . However, in order to calculate the functional˜Φ[ ρ ˜ H ] an optimal unitary transformation in above senseis calculated and the h R ikk ′ are found by a projection ofa Hartree-Fock solution onto the basis states of R . Weassume that the effective energy of the impurity site isthe same as in the original model (˜ ε d = ε d ). ∗ [email protected] P. W. Anderson, Phys. Rev. , 41 (1961). H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson,Phys. Rev. B , 1003 (1980). C. M. Varma and Y. Yafet, Phys. Rev. B , 2950 (1976). O. Gunnarsson and K. Sch¨onhammer,Phys. Rev. B , 4315 (1983). D. Jacob, K. Haule, and G. Kotliar,Phys. Rev. Lett. , 016803 (2009). T. O. Wehling, A. V. Balatsky, M. I. Kat-snelson, A. I. Lichtenstein, and A. Rosch,Phys. Rev. B , 115427 (2010). B. Surer, M. Troyer, P. Werner, T. O. Wehling,A. M. L¨auchli, A. Wilhelm, and A. I. Lichtenstein,Phys. Rev. B , 085114 (2012). J. K¨ugel, M. Karolak, J. Senkpiel, P.-J. Hsu, G. Sangio-vanni, and M. Bode, Nano Letters , 3895 (2014). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Reviews of Modern Physics , 13 (1996). A. I. Lichtenstein and M. I. Katsnelson,Phys. Rev. B , 6884 (1998). C. J. Ballhausen, MacGraw-Hill, New York (1962). F. M. F. Groot, J Elec. Spectrosc. Relat. Phenom. , 529 (1994). E. Gull, A. J. Millis, A. I. Lichtenstein,A. N. Rubtsov, M. Troyer, and P. Werner,Rev. Mod. Phys. , 349 (2011). R. Bulla, T. Costi, and T. Pruschke,Rev. Mod. Phys. , 395 (2008). M. Caffarel and W. Krauth,Phys. Rev. Let. , 1545 (1994). Q. Si, M. J. Rozenberg, G. Kotliar, and A. E. Ruckenstein,Phys. Rev. Lett. , 2761 (1994). D. Zgid, E. Gull, and G. K.-L. Chan,Phys. Rev. B , 165128 (2012). Y. Lu, M. H¨oppner, O. Gunnarsson, and M. W. Haverkort,Phys. Rev. B , 085102 (2014). C. Lin and A. A. Demkov,Phys. Rev. B , 235122 (2014). M. Potthoff, M. Aichhorn, and C. Dahnken,Phys. Rev. Lett. , 206402 (2003). A. A. Kananenka, E. Gull, and D. Zgid,Phys. Rev. B , 121111 (2015). R. E. Peierls, Phys. Rev. , 918 (1938). N. N. Bogoliubov., Dokl. Akad. Nauk SSSR , 244(1958). R. P. Feynman,
Statistical Mechanics (Benjamin, ReadingMass., 1972). A sequential least squares programming algorithm as im-plemented in the package scipy.optimize.minimize . A simplex algorithm as implemented in the python modulescipy.optimize.minimize . J. C. Slater,
Quantum theory of atomic structure , Vol. 1(McGraw-Hill New York, 1960). R. Eder, in
Correlated electrons: from models to materials ,Schriften des Forschungszentrums J¨ulich. Reihe Modelingand simulation, edited by E. Pavarini, E. Koch, F. Anders,and M. Jarrell (Forschungszentrum J¨ulich GmbH, Institutefor Advanced Simulation, 2012). D. S´en´echal, Phys. Rev. B , 235125 (2010). D. Kraft,
A software package for sequential quadratic pro-gramming (DFVLR K¨oln, Germany, 1988). J. A. Nelder and R. Mead,The Computer Journal7