# Variational exact diagonalization method for Anderson impurity models

###### Abstract

We describe a variational approach to solving Anderson impurity models by means of exact diagonalization. Optimized parameters of a discretized auxiliary model are obtained on the basis of the Peierls-Feynman-Bogoliubov principle. Thereby, the variational approach resolves ambiguities related with the bath discretization, which is generally necessary to make Anderson impurity models tractable by exact diagonalization. The choice of variational degrees of freedom made here allows 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 energy calculations and to systematically optimize and improve on traditional routes to the discretization problem such as fitting of hybridization functions on Matsubara frequencies. Benchmarks in terms of a single orbital Anderson model demonstrate that the variational exact diagonalization method accurately reproduces free energies as well as several single- and two-particle observables obtained from an exact solution. Finally, we demonstrate the applicability of the variational exact diagonalization approach to realistic five orbital problems with the example system of Co impurities in bulk Cu and compare to continuous-time Monte Carlo calculations. The accuracy of established bath discretization schemes is assessed in the framework of the variational approach introduced here.

###### pacs:

72.80.Rj; 73.20.Hb; 73.61.Wp## I Introduction

The Anderson impurity modelAnderson (1961) (AIM) is a general model for the description of interacting impurities in metallic host systems. Originally, it was developed to describe single atoms with open d- or f-shells embedded in bulk materials and to understand the formation of their magnetic momentsAnderson (1961). Furthermore the model includes the widely discussed Kondo physics Krishna-murthy et al. (1980). Multi orbital variants of the AIM gained considerable attention in the context of rare-earth impurity systemsVarma and Yafet (1976); Gunnarsson and Schönhammer (1983) as well as more recently magnetic adatoms or molecules on surfacesJacob et al. (2009); Wehling et al. (2010); Surer et al. (2012); Kügel et al. (2014). Finally, dynamical mean field theoryGeorges et al. (1996) (DMFT) links correlated bulk systems as well as nanostructures to Anderson impurity models.

To address the electronic structure of realistic correlated electron materials one often resorts to LDA++ approachesLichtenstein and Katsnelson (1998), where quantum lattice or impurity models are derived from first principles calculations. The resulting models are typically multi-orbital models including complex hybridization between the impurity and a continuous bath of states from the surrounding material, which brings along two challenges: First, the numerical solution of the impurity models and second the interpretation of the physics contained in these generally complex models in more simple terms. Experiments are for instance often interpreted in terms of atomic spins, crystal field, ligand fieldBallhausen (1962) or cluster approachesGroot (1994), which typically involve a small discrete set of bath states or no bath states at all. The link of the complex, ab initio derived models and simpler phenomenological models is a priori unclear and relates to the so-called bath discretization problem of exact diagonalization solvers of the AIM.

The solution of the Anderson impurity model for general parameters has to be done numerically by means of e.g. quantum Monte CarloGull et al. (2011) (QMC), numerical renormalization groupBulla et al. (2008) (NRG), or exact diagonalization (ED) methods.Georges et al. (1996) While NRG and QMC are in principle numerically exact methods, they become computationally very demanding, when dealing with many orbitals, hybridization functions with low symmetry, spin orbit coupling and general fermionic four operator Coulomb vertices. ED methods deal with low symmetries and general Coulomb vertices at no additional computational cost but suffer from the so-called bath discretization problem: Due to the exponential growth of the many particle Fock space with the system size, it can handle only a few bath levels per orbital. A mapping of the continuous bath to a discrete version has to be found. Several approaches to this task have been introduced. One is to fit the hybridization function of the continuous bath on Matsubara frequencies Caffarel and Krauth (1994), another is to represent the hybridization function by a continued fraction and to link its coefficients to the parameters of the bath Si et al. (1994). These schemes are systematic in the sense that they converge to the full model when including more and more bath sites. However, in the multi-orbital case, the number of bath sites is limited (typically on the order of three or less for a five orbital impurity problem), so the quality of the mapping can hardly be checked by an analysis of the convergence.

Basically two different strategies have been laid out to circumvent this problem. First, the many body Hilbert space can be truncated in the sense of configuration interaction (CI) expansions, which have a long tradition in the context of quantum impurity problemsVarma and Yafet (1976); Gunnarsson and Schönhammer (1983) and are subject of recent developments.Zgid et al. (2012); Lu et al. (2014); Lin and Demkov (2014) CI expansions are variational, i.e. they deliver upper bounds for total energies, but they do not provide simplified auxiliary Hamiltonians. On the other hand, there are several approaches towards optimized cluster approximations to Anderson impurity problems. In this context, self-energy functional theoryPotthoff et al. (2003) is based on an extremal principle but it is not variational regarding total energies and does not allow for variations of interaction parameters or the interacting orbitals. More general optimizations are possible in the framework of the so-called self-energy embedding theory (SEET)Kananenka et al. (2015), which is however not variational.

In this paper, we combine ideas of variational approaches and optimized cluster approximations to the AIM. We introduce a strictly variational method of approximating an AIM with continuous bath by an AIM with finite strongly reduced number of bath sites, which we call variational ED method. It guaranties an optimal approximation to the AIM for a given number of bath sites in the sense of thermodynamic ground state properties. The method is based on the well-known Peierls-Feynman-Bogoliubov variational principle Peierls (1938); Bogoliubov. (1958); Feynman (1972), which finds optimal effective models on the basis of an optimal density matrix by minimizing a free energy functional.

We will introduce the AIM and the variational principle in Sec. II, where we also explain the details of calculating the Peierls-Feynmann-Bogoluibov free energy functional and how to minimize it efficiently. By treating a single orbital model with the variational ED method in Sec. III we analyze its performance in comparison to an exact treatment, established bath discretization methodsCaffarel and Krauth (1994) as well as Hartree-Fock theory. In Sec. IV, we demonstrate the applicability of the method to realistic five orbital system with the example of Co impurities in bulk Cu and compare to QMC simulations. We show that the variational ED method leads to systematically lower, i.e. more accurate, free energy estimates than unrestricted Hartree-Fock and traditional bath discretization schemes also in the multiorbtial case. Conclusions and outlook are given in Sec. VI.

## Ii Model and method

After introducing the Anderson impurity model we will recapitulate the Peierls-Feynman-Bogoliubov variational principle and show how to apply it to discretize Anderson impurity models in an optimal manner.

### ii.1 The Anderson impurity model

The Hamiltonian of the initial AIM (termed “original model” hereafter) reads

(1) |

The bath is described by

(2) |

where is the energy of the bath state with band/orbital index and some additional quantum number . is the corresponding particle number operator. The hybridization part

(3) |

couples the bath sites of one band to an orbital of the impurity with a coupling strength . The bath electrons with spin are created and annihilated by and , respectively, while () denote the creation (annihilation) operators of the impurity electrons. The impurity site is described by

(4) |

which contains the on-site Coulomb interaction and the on-site energies . By integrating out all bath degrees of freedom we arrive at the hybridization function

(5) |

which describes the energy dependent coupling of the impurity to the bath.

### ii.2 Peierls-Feynman-Bogoliubov variational principle

Given a Hamiltonian , which is “difficult” to solve, we search for an optimal approximation to within a set of simpler effective Hamiltonians . The Peierls-Feynman-Bogoliubov variational principlePeierls (1938); Bogoliubov. (1958); Feynman (1972) provides us with a prescription on how to fix the parameters of in a thermodynamically optimal way, i.e. such that the canonical density matrix resulting form approximates the density matrix corresponding to as close as possible. More strictly speaking: the canonical density operator of the auxiliary system, where is the partition function, approximates the exact density operator derived from as close as possible, when the Peierls-Bogoliubov-Feynman functional

(6) |

becomes minimal. Here is the free energy of the effective system. denotes a thermodynamic expectation value with respect to the effective system. In the case of the functional becomes minimal and coincides with the free energy of the original system. In our case represents the full AIM, Eqs. (1)-(4), and is the model with discretized bath, which is now introduced.

### ii.3 Effective Hamiltonian

The structure of the effective Hamiltonian for the case of a single impurity orbital is depicted in the right panel of Fig. 1. In contrast to the original model (left panel of Fig. 1), the effective model consists of two decoupled parts: First, the effective impurity coupled to one bath site only and second the remaining bath sites. I.e. we partition the full Hilbert space into a correlated subspace (first part) and an uncorrelated rest (second part). In this work, we consider for concreteness a cluster consisting of a multi-orbital impurity and one bath site per impurity orbital for the correlated space but other choices are similarly possible. The single particle states of the effective model are related to those of the original model by a unitary transformation, which allows for mixing of original “bath” and “impurity” character in the effective model.

The optimal matrix elements of the effective model, as well as the optimal unitary transformation are found by minimizing the functional (6).

The states spanning are defined by

(7) | ||||

(8) |

where the coefficients are chosen such that and form an orthonormal basis of . An orthonormal basis spanning is defined by

(9) |

As a whole, the coefficients form a unitary matrix. In practice, we obtain the elements of this matrix from the QR decomposition of a matrix, in which the first two rows are defined by the coefficients of and and all other elements are zero. This leads to a new orthonormal basis for the full space , which provides the partitioning according to . The ansatz for the effective Hamiltonian in this new basis explicitly reads

(10) |

with

(11) |

and

(12) |

It is stressed, that the new states are linear combinations of the original impurity and bath states, leading to mixed basis states. The new impurity states can have some amount of bath character and vice versa. The Hamiltonian in Eq. (11) states a many-body problem which can be solved by exact diagonalization, as long as its Hilbert space is sufficiently small. In contrast, the Hamiltonian (12) states a one-particle problem and can be solved by diagonalizing the matrix . In summary, the Hamiltonian defines an effective Hamiltonian, which can be solved exactly and thus the functional (6) can be calculated.

This ansatz implies several approximations. First, all couplings between and are neglected. Second interaction terms are restricted to new effective impurity orbitals within , which is motivated by the fact that the original model includes only on-site interactions too. The latter approximation can be relaxed to include arbitrary interactions within , but we keep it here for simplicity.

Finally, we note that the amount of variational degrees of freedom in the variational ED approach is such that it includes Hartree-Fock as the limiting case . Thus, we expect that variational ED will generally give more accurate energy estimates than Hartree-Fock.

### ii.4 Implementation

In order to perform the minimization in practice, the number of free parameters has to be kept sufficiently low. First, for the rest of this work it is assumed that the Coulomb tensor is not varied. We choose it to be the same as in the original model. Test calculations have shown, that the variation of the Coulomb tensor is not crucial, as this can mostly be absorbed into the variation of the impurity level. The single particle matrix elements of are assumed to be free parameters. In principle, the parameters of the uncorrelated Hamiltonian are free parameters, too. However, to further reduce the number of free parameters, we define by a projection of a Hartree-Fock solution of the original Hamiltonian onto the states . The Hartree-Fock solution of the original Hamiltonian (1) can be written as

(13) |

where the eigenstates and energies are found by applying the Hartree-Fock decoupling

(14) |

to (4) and solving the resulting non-interacting problem self-consistently. The single particle matrix elements within the uncorrelated space explicitly read

(15) |

In order to not break any spin rotation symmetries, restricted Hartree-Fock is used.

The functional now depends on the unitary transformation and on the matrix elements of . The minimum of the functional is searched by iterative methods. Thus, the functional has to be calculated for various points of the variational space with the computationally most expensive part being here the diagonalizations of . Therefore, we first search for fixed parameters in a corresponding optimal unitary transformation matrix defining the optimal partitioning using an SLSQP algorithm ^{1}^{1}1A sequential least squares programming algorithm as implemented in the package scipy.optimize.minimize Kraft (1988). The search of the minimum w.r.t. the parameters of is then done by the Nelder-Mead algorithm ^{2}^{2}2A simplex algorithm as implemented in the python module scipy.optimize.minimizeNelder and Mead (1965). The number of independent parameters can be further reduced when the original system shows symmetries like orbital degeneracies which are assumed not to be broken in the effective model.

## Iii Benchmark for a single orbital AIM

In this section the variational ED method is tested for its performance in reproducing the density operator as well as observables such as the occupation number, double occupancy and crystal orbital overlap populations of a simple original model. The original model we consider here is a single orbital model with only 6 bath sites, which itself can be solved by exact diagonalization. The detailed setup of the model is as follows: The impurity level is , the interaction strength is . The 6 bath levels are equally aligned around a mean bath energy in an interval of (i.e. the bandwidth of the bath). The coupling is . The mean bath energy is swept from to . All energies are measured w.r.t. to the Fermi energy . The system is solved for . The model is first solved exactly, second by the variational ED method, third by unrestricted Hartree-Fock and finally by ED using reduced bath sites obtained by fitting of the hybridization functions on the imaginary Matsubara frequency axisCaffarel and Krauth (1994). The latter type of approaches require generally the introduction of a so-called weight function for the fitting procedure, as explained in the appendix (A).

The central object for the assessment of the quality of the methods is the difference between and the exact free energy , as shown in Fig. 2a). For bath sites energetically far away from the Fermi level and from single particle excitation energies of the impurity (), all methods lead essentially to the correct free energy. Deviations occur, however, for bath levels closer to the Fermi energy. The Hartree-Fock free energy differs from the exact thermodynamical potential on the order of basically in the whole range of . The fitting of the hybridization function function on the imaginary axis leads to rather accurate free energies as long as all bath sites are above or below the Fermi energy (), while for deviations from the exact thermodynamical potential on the order of to occur. The choice of an optimal weight function (see appendix A) depends on details of the bath: For the case of bath sites on both sides of the Fermi level () leads to the lowest free energies. Otherwise, the constant weight function shows smallest deviations of the free energy functional from the exact solution. The variational ED method is generally very close to the exact solution. Only for the special case of a strictly symmetric distribution of the bath sites around the Fermi energy () a deviation on the order of meV occurs.

Fig. 2b)-e) shows a comparison of several observables (chemical bond strength b), bath occupation c), impurity occupation d) and double occupation e)) calculated with the different methods. For the outer most regions () all methods describe the observables accurately. The fitting of hybridization functions on the Matsubara axis leads to deviations depending on the weight function, especially for the double occupation and chemical bond strength if the bath is centered around the Fermi energy (). Hartree-Fock systematically overestimates the double occupancy for . The variational ED method shows nearly no deviations from the exact solution at all.

It is instructive to examine the unitary transformation linking the basis of the original and effective model. Fig. 2f) shows the coefficients of the linear combination of the states spanning the correlated space and (see Eqs. (7) and (8)) for an original model with the bath centered around the energy . The effective impurity has mainly character with small bath admixture and can approximately be interpreted as the old impurity state. The coupled effective bath state is nearly a pure linear combination of old bath states, where states closer to the Fermi energy contribute stronger than those further away. This behavior is very reminiscent of effective bath wave functions obtained in variational approaches like the Varma-YafetVarma and Yafet (1976) or the Gunnarsson-Schönhammer expansion.Gunnarsson and Schönhammer (1983)

For the treatment of original models with far more bath sites, it is important to note, that the coefficients defining the unitary transformation from the original bath states to the effective impurity and bath orbitals, i.e. and , vary smoothly as function of the bath energies on either side of the Fermi energy.

## Iv Co impurities in Cu: Application to a realistic five orbital system

### iv.1 The AIM derived from LDA

In this section the variational ED method is applied to a realistic model of Co impurities in bulk Cu, which has been obtained from super-cell DFT calculations and has been analyzed using a QMC impurity solver in Ref. Surer et al., 2012.

The cubic symmetry of the Cu crystal leads to a splitting of the Co -orbitals into blocks of t and e symmetry. From the DFT hybridization function, which is a continuous function, we obtain our initial model assuming some large number of bath sites, here 100 per orbital. (This number does not present a limiting factor and could be chosen arbitrarily larger). The bath sites are assumed to be equidistantly distributed between and , and the hybridization terms are then found by fitting the imaginary part of a discretized hybridization function

(16) |

with some broadening to the ab initio hybridization function on the real axis. The are plotted in Figure 3a). The crystal field obtained from the DFT calculation is . As in Ref. Surer et al., 2012, we consider a rotationally invariant Coulomb interaction defined by

(17) |

where are the Gaunt coefficients Slater (1960); Eder (2012) and where , and are Slater parameters with the average Coulomb interaction and Hunds exchange interaction . Due to the so-called double counting problem inherent to LDA++ approaches, the filling of the impurity d-levels is not exactly known. Here, we consider the double counting potential as in Ref. Surer et al., 2012. All data is obtained at a inverse temperature of , like in the case of the QMC simulations. Finally, we assume that the cubic symmetry of the system prevails, which means that only two independent sets of matrix elements (for the t and e states) have to be varied during the minimization of .

### iv.2 Implementation of the variational ED method for the 5 orbital AIM

We compare two different sets of variational degrees of 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 expansion coefficients , and . This leads to considerably less variational parameters and a much smaller amount of expectation values to be calculated in each step of the iteration. In the second approach, “all”, which is computationally more demanding because the full two-particle density matrix of the effective system has to be calculated, we optimize the full one particle basis of the bath and that of the impurity.

Because a full optimization of the parameters of the effective model is computationally challenging, it is crucial to start the optimization from a good initial guess. We obtain such initial guesses for the parameters of the bath by fitting of hybridization functions on Matsubara frequencies as introduced in the appendix A. We choose as the initial guess for the parameters of the impurity. The resulting first guesses using different weight functions are summarized in the Table 1. While all weight functions lead to setups with the t bath sites above the Fermi energy and the e bath sites below, the details of their energetic positions and the hybridization strengths depend strongly on the form of . Adding more weight on features on small Matsubara frequencies shifts the effective bath parameters to smaller values. The quality of these starting guesses in the context of the variational principle is discussed in the next section.

weight function | var | |||
---|---|---|---|---|

(eV) | ||||

(eV) | ||||

(eV) | ||||

(eV) | ||||

(eV) | ||||

(eV) |

The large number of bath sites (100 per orbital) in the original model leads to 202 variational parameters defining the unitary transformation in the “all” case or 100 parameters in the “bath” case for each orbital. The observation, that are smooth functions of energy (c.f. Fig. 2e) below and above , leads to the possibility of expanding them in a set of smooth functions and thereby reducing the number of variational parameters considerably. Here, we chose five Chebyshev polynomials per orbital for bath sites above and five for those below the Fermi energy. Therefore only 22 (or 10 in the case of “bath”) parameters per orbital have to be varied to find the optimal unitary transformation to embed the effective model into the full Hilbert space.

### iv.3 Results

We will first compare free energy estimates as well as different local observables obtained from variational ED treatments to unrestricted Hartree-Fock as well as QMC calculations. Afterwards, we investigate the nature of the optimized effective bath and impurity states as obtained from the variational ED treatment.

#### iv.3.1 Free energy functional and local observables

(eV) | |||
---|---|---|---|

QMC | 7.78 0.05 | 0.92 0.02 | |

UHF | 0 | 7.78 | 1.78 |

RHF | 0.52 | 8.20 | 1.06 |

fit0,bath | 0.05 | 7.75 | 1.06 |

fit1,bath | 1.15 | 7.84 | 1.04 |

fit2,bath | 2.77 | 7.92 | 1.03 |

fit0,all | -0.22 | 7.71 | 1.03 |

fit1,all | 0.32 | 7.65 | 1.05 |

fit2,all | 1.04 | 7.60 | 1.00 |

var,bath | 0.00 | 7.76 | 1.05 |

var,all | -0.30 | 7.75 | 1.02 |

Table 2 shows the free energy functional (relative to unrestricted Hartree-Fock (UHF)), as obtained with different starting points and different amounts of variational degrees of freedom in the variational ED approach. In the case of “bath”, the constant weight function (“fit0”, ) leads to the lowest values of the . The models derived using weight functions (“fit1”) and (“fit2”) lead to free energy estimates which are about to higher in energy. The situation for the case of “all” is similar. On this basis, we have chosen the starting guess obtained with the constant weight function for the full optimization of the effective model parameters. The resulting parameters are shown in the last column of Tab. 1 and are close to the starting 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 (, see Tab. 2), we see that the description by unrestricted Hartree-Fock is rather close to QMC, whereas restricted Hartree-Fock overestimates the occupation. All versions of exact diagonalization lead to occupations close to the QMC results, and many cases within the QMC error bars. The spin (defined as ), which is a two-particle observable, reveals the problems of the Hartree-Fock description. is vastly overestimated by unrestricted Hartree-Fock. The variational ED methods, especially the “all” case for the constant weight function (“fit0”) and the full optimized ED model, lead to results close to QMC.

To compare the results of the variational ED method with those ED methods based on fitting of the hybridization function on the imaginary axis, we should compare the “fit0,bath”,“fit1,bath”, and “fit2,bath” cases to the corresponding “all” and “var,all” cases. We see that having more variational degrees of freedom leads to an improved description of the free energies, as it should be.

In general, we learn that only in the case of optimizing both effective bath and effective impurity states (termed “all”) we reach lower values of the free energy functional than with unrestricted Hartree-Fock: The freedom to form mixtures of bath and impurity states in the effective model is important to describe the free energy and local observables of the system adequately. As the variational ED method provides more accurate (free) energy estimates than unrestricted Hartree-Fock, the approach introduced here could be a way to improve LDA+U total energy schemes.

#### iv.3.2 The effective basis states

Now we analyze the unitary transformation relating the optimized basis states of the effective model and the original basis states. The transformation obtained for the models from the starting guesses with weight functions , , and is shown in Fig. 3 b), c), and d), respectively. We observe a clear trend that the admixture of the original bath states into the effective bath state (, blue lines) is strongest in the vicinity of the effective bath site energy . For bath states close to the Fermi energy we get a sharp cut off for states on the opposite side of the Fermi energy. This is very similar to first order configuration interaction treatments of the AIMGunnarsson and Schönhammer (1983); Varma and Yafet (1976). The original impurity admixture in the effective impurity () rises with the distance of the effective bath site from the Fermi energy.

## V Spectral functions

The variational principle results in an effective model which represents thermodynamic ground state properties in an optimal manner. This is a necessary but not a sufficient condition to give a good approximation also for excitation spectra. In the following, we study the one particle spectral function for single orbital impurity benchmark systems from Sec. III with the bath states centered around and two different hybridization strengths, and , respectively. The impurity spectral function is obtained from the Lehmann representation of the impurity Green’s function

(18) |

where in our calculations is replaced by a broadening of and the inverse temperature is , which is very close to the calculations of expectation values in Sec. III.

We assess the quality of the spectra obtained from the variational ED method, from ED with the hybridization function fitted on the Matsubara axis (with two different weight functions and ) and from an unrestricted Hartree-Fock treatment by comparing to the exact spectrum of the original model. For the case of , we additionally compare the spectra for different amounts of variational freedom in choosing the basis states of the effective model, where we either optimized 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 the case of strong hybridization (, Fig. 4 a)) are two major peaks at about and (stemming from bonding and anti-bonding combinations of impurity and bath orbitals), two satellite peaks far from the Fermi energy and additional smaller peaks around the Fermi energy. The spectral function from Hartree-Fock reproduces the bonding/anti bonding peaks and those close to the Fermi energy very well, while the satellites are missing. The variational ED method describes the positions of the main peaks well and also reproduces the satellite peaks, whereas the minor peaks around the Fermi energy are not present. The spectral function from a fit on the imaginary axis with a constant weight function () shows a similar picture but with major peaks and satellites shifted considerably towards the Fermi energy. The result for the weight function emphasizing small Matsubara frequencies () leads to a good representation of the peaks around and and some minor peaks around the Fermi energy but the satellites are missing completely. The resulting spectrum obtained by not mixing bath and impurity states (“bath”) shows only little resemblance to the original spectrum. I.e. in this case the mixing of bath and impurity basis states can not only improve total energies / thermodynamic potentials but also spectra quite significantly.

In the case of weaker hybridization (, Fig. 4 b)) the original impurity spectral function shows two Hubbard peaks at about and and in comparison to the former case more spectral weight and additional features close to the Fermi energy. Here, the Hartree-Fock description results in a spin-polarized ground state and describes the positions of the Hubbard peaks in the spectrum correctly. The enhanced spectral weight at the Fermi energy is however not reproduced. Exact diagonalization of the effective models obtained from the variational method and from the fits of the hybridization functions on imaginary axis give similar results. In addition to the upper and lower Hubbard peaks the ED methods also reproduce enhanced spectral weight at the Fermi level.

While the performance of the fit methods in reproducing the spectra of the original model differs between the case with strong and weak hybridization particularly for the weight function , the variational method gives satisfactory results in both cases. The spectra from variational ED are in both cases at least as close to the spectra of the original model as the best spectrum obtained with any of the two bath fitting procedures or ) under investigation.

## Vi Conclusion

In conclusion, we present a variational exact diagonalization method which provides self-consistently optimized parameters of discretized Anderson impurity models. The method is based on an optimal partitioning of the system into a correlated part, where electronic interactions are explicitly taken into account, and an uncorrelated rest and on finding optimal effective Hamiltonians for both parts of the system. To this end, a variation of a free energy functional w.r.t. one-particle basis states spanning the correlated subspace and the matrix elements of the effective Hamiltonians is performed.

A benchmark of the variational ED method against an exact solution of a one orbital Anderson model demonstrates its excellent performance in reproducing ground state observables of the impurity and bath and additionally a sound performance in reproducing the impurity spectral function. A comparison with Hartree-Fock and established bath discretization schemes for ED shows that the variational approach introduced here even works for difficult cases, i.e. when the bath is symmetric around the Fermi energy. Furthermore, applicability of the variational ED method to realistic multi-orbital cases is demonstrated with the example of Co impurities in bulk Cu. Also here, the variational method leads to an accurate description of local one and two particle observables like the impurity occupation and the spin. Energetically the method outperforms unrestricted Hartree-Fock, which suggests that the variational ED approach could be useful to improve total-energy approaches to correlated systems beyond LDA+U. Finally, the method introduced, here, can be used to embed established bath discretization schemes such as the fit of hybridization functions on Matsubara frequencies into a variational framework and to reach an unbiased decision of e.g. which weight function to choose in the fits of the hybridization functions. In the example studied, here, the constant weight function leads to best results in terms of the free energy, whereas leads to results qualitatively similar to a first order configuration interaction expansionGunnarsson and Schönhammer (1983); Varma and Yafet (1976). For the description of systems closer to the atomic limit, like Co on Cu or 4f systems we expect CI expansions to converge faster and weight functions like or might be a better choice. The presented variational approach will allow for an unbiased decision in any case.

Additionally, the variational method is universal and different choices of the effective correlated space are possible. To gain yet higher accuracy, more bath levels could be included. On the other hand, to make contact with ligand field theory or crystal field theory descriptions of magnetic impurity systems and nanostructures one could consider correlated subspaces with very few or without any effective bath orbitals at all.

## Vii Acknowledgments

The authors thank M. Katsnelson, A. Lichtenstein, M. Potthoff, P. Blöchl, R. Schade ans S. Barthel for useful discussions as well as the Central Research Development Fund of the University of Bremen and the DFG via FOR 1346 for financial support.

## Appendix A Fit of hybridization functions on the Matsubara axis

The method of fitting hybridization functions is shortly introduced for the sake of completeness. The basic idea is to minimize a cost function for the inverse impurity Green’s function (or equivalently the hybridization function) of the discretized model and that of the original model, both defined on the imaginary frequency axis Caffarel and Krauth (1994). In the case of one effective bath site, the discrete impurity Green’s function is defined as

(19) |

and the Green’s function of the original model as

(20) |

The cost function then reads

(21) |

where is a weight function. Popular choices for the weight function are , and . Different weight functions put different emphasis of low/higher Matsubara frequencies Sénéchal (2010). Throughout this work, we have chosen and . This method only provides the effective parameters and . However, in order to calculate the functional an optimal unitary transformation in above sense is calculated and the are found by a projection of a Hartree-Fock solution onto the basis states of . We assume that the effective energy of the impurity site is the same as in the original model ().

## References

- Anderson (1961) P. W. Anderson, Phys. Rev. 124, 41 (1961).
- Krishna-murthy et al. (1980) H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson, Phys. Rev. B 21, 1003 (1980).
- Varma and Yafet (1976) C. M. Varma and Y. Yafet, Phys. Rev. B 13, 2950 (1976).
- Gunnarsson and Schönhammer (1983) O. Gunnarsson and K. Schönhammer, Phys. Rev. B 28, 4315 (1983).
- Jacob et al. (2009) D. Jacob, K. Haule, and G. Kotliar, Phys. Rev. Lett. 103, 016803 (2009).
- Wehling et al. (2010) T. O. Wehling, A. V. Balatsky, M. I. Katsnelson, A. I. Lichtenstein, and A. Rosch, Phys. Rev. B 81, 115427 (2010).
- Surer et al. (2012) B. Surer, M. Troyer, P. Werner, T. O. Wehling, A. M. Läuchli, A. Wilhelm, and A. I. Lichtenstein, Phys. Rev. B 85, 085114 (2012).
- Kügel et al. (2014) J. Kügel, M. Karolak, J. Senkpiel, P.-J. Hsu, G. Sangiovanni, and M. Bode, Nano Letters 14, 3895 (2014).
- Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Reviews of Modern Physics 68, 13 (1996).
- Lichtenstein and Katsnelson (1998) A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B 57, 6884 (1998).
- Ballhausen (1962) C. J. Ballhausen, MacGraw-Hill, New York (1962).
- Groot (1994) F. M. F. Groot, J Elec. Spectrosc. Relat. Phenom. 67, 529 (1994).
- Gull et al. (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
- Bulla et al. (2008) R. Bulla, T. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
- Caffarel and Krauth (1994) M. Caffarel and W. Krauth, Phys. Rev. Let. 72, 1545 (1994).
- Si et al. (1994) Q. Si, M. J. Rozenberg, G. Kotliar, and A. E. Ruckenstein, Phys. Rev. Lett. 72, 2761 (1994).
- Zgid et al. (2012) D. Zgid, E. Gull, and G. K.-L. Chan, Phys. Rev. B 86, 165128 (2012).
- Lu et al. (2014) Y. Lu, M. Höppner, O. Gunnarsson, and M. W. Haverkort, Phys. Rev. B 90, 085102 (2014).
- Lin and Demkov (2014) C. Lin and A. A. Demkov, Phys. Rev. B 90, 235122 (2014).
- Potthoff et al. (2003) M. Potthoff, M. Aichhorn, and C. Dahnken, Phys. Rev. Lett. 91, 206402 (2003).
- Kananenka et al. (2015) A. A. Kananenka, E. Gull, and D. Zgid, Phys. Rev. B 91, 121111 (2015).
- Peierls (1938) R. E. Peierls, Phys. Rev. 54, 918 (1938).
- Bogoliubov. (1958) N. N. Bogoliubov., Dokl. Akad. Nauk SSSR 119, 244 (1958).
- Feynman (1972) R. P. Feynman, Statistical Mechanics (Benjamin, Reading Mass., 1972).
- (25) A sequential least squares programming algorithm as implemented in the package scipy.optimize.minimize Kraft (1988).
- (26) A simplex algorithm as implemented in the python module scipy.optimize.minimizeNelder and Mead (1965).
- Slater (1960) J. C. Slater, Quantum theory of atomic structure, Vol. 1 (McGraw-Hill New York, 1960).
- Eder (2012) R. Eder, in Correlated electrons: from models to materials, Schriften des Forschungszentrums Jülich. Reihe Modeling and simulation, edited by E. Pavarini, E. Koch, F. Anders, and M. Jarrell (Forschungszentrum Jülich GmbH, Institute for Advanced Simulation, 2012).
- Sénéchal (2010) D. Sénéchal, Phys. Rev. B 81, 235125 (2010).
- Kraft (1988) D. Kraft, A software package for sequential quadratic programming (DFVLR Köln, Germany, 1988).
- Nelder and Mead (1965) J. A. Nelder and R. Mead, The Computer Journal 7, 308 (1965).