# Superperturbation solver for quantum impurity models

###### Abstract

We present a very efficient solver for the general Anderson impurity problem. It is based on the perturbation around a solution obtained from exact diagonalization using a small number of bath sites. We formulate a perturbation theory which is valid for both weak and strong coupling and interpolates between these limits. Good agreement with numerically exact quantum Monte-Carlo results is found for a single bath site over a wide range of parameters. In particular, the Kondo resonance in the intermediate coupling regime is well reproduced for a single bath site and the lowest order correction. The method is particularly suited for low temperatures and alleviates analytical continuation of imaginary time data due to the absence of statistical noise compared to quantum Monte-Carlo impurity solvers.

###### pacs:

71.27.+a## I Introduction

Quantum impurity models have been widely used in condensed matter physics. Examples comprise the study of the Kondo-effectHewson (1993), or of adatoms on surfacesSavkin et al. (2005). In particular, the success of the dynamical mean-field theory (DMFT) to describe strongly correlated systems has triggered efforts to develop efficient solvers for the impurity problem. In DMFT, the lattice problem is mapped onto an effective quantum impurity problem which needs to be solved repeatedly to satisfy a self-consistency condition. Generalizations of the DMFT, such as cluster DMFT or dynamical cluster approximations require highly efficient methods for the solution of the multiorbital quantum impurity problem.

The combination of the DMFT with electronic structure methods, such as the local density approximation (LDA+DMFT approach)Kotliar et al. (2006) requires the solution of impurity models with general types of the interactions, such as the full Coulomb vertex. The problem has been solved using approximate methods such as a multiorbital version of the fluctuation-exchange approximation (FLEX) for weak couplingKatsnelson et al. (2008). A strong-coupling solver based on an expansion in the impurity-bath hybridization has been proposed in order to address systems with open d or f-shells or Mott insulatorsDai et al. (2005). Recently, next-generation Monte-Carlo methods provide a numerically exact solutionRubtsov et al. (2005); Werner et al. (2006); Haule (2007). While being suitable for general interactions, these methods however require a sizeable numerical effort and analytical continuation is complicated through statistical noise in the imaginary time data.

In this article we propose an efficient solver for the impurity problem which is suitable for both weak and strong coupling and general interactions. It is essentially based on a “superperturbation”, i.e. a perturbation on top of a solution obtained by exact diagonalization (ED) for a small number of bath levels. The method alleviates analytical continuation due to the absence of noise and is suitable for the study of, e.g., multiplet effects in solids.

## Ii Formalism

For notational convenience, we introduce the formalism for the single-orbital case. It was introduced earlier in the context of lattice fermion models to include spatial correlations beyond DMFT in Ref. Rubtsov et al., 2008. A generalization of the underlying concepts to the multiorbital case can be found in Refs. Hafermann et al., 2007; Brener et al., 2008. In a complementary approach, named dynamical vertex approximation, the single- and two-particle Green functions of the impurity were computed using EDToschi et al. (2007). The Hamiltonian of the Anderson impurity model (AIM) reads

(1) |

where the local impurity degrees of freedom represented by , couple to a bath of free conduction electrons with dispersion via a hybridization . stands for any local interaction. Since is quadratic in the bath operators ,, they can be integrated out giving rise to the action (henceforth we switch to the path integral representation):

(2) |

Here is the chemical potential and the sum is over Matsubara frequencies , where is the inverse temperature. In ED the continuous dispersion of the bath is approximated by a finite number of bath levels, which corresponds to replacing the hybridization function by its discrete counterpart

(3) |

The problem now is to determine the parameters , such that is the best approximation to the original hybridization . Mathematically, this corresponds to projecting the hybridization function onto the discrete subspace spanned by the parameters ,. Various methods have been proposed for this purposeGeorges et al. (1996); Koch et al. (2008). One way is to perform a conjugate gradient minimization of, e.g., the distance function

(4) |

where the sum is over Matsubara frequencies. The parameter , if chosen large (e.g. ), enhances the importance of the lowest Matsubara frequencies in the minimization procedure. This parameter is particularly important when a small number of bath sites is used to approximate the original hybridization.

The quality of the approximation can be measured by and is determined by the number of bath sites used to represent the original hybridization. On the other hand, the number of bath sites is limited due to the exponential growth of the Hilbert space. Instead of approximating by a large number of bath sites, we shall follow a different route and rewrite the action Eqn. 2 in terms of an exactly solvable part and an additional bilinear term:

(5) | |||||

The exactly solvable part (the first line of Eqn. 5) we will henceforth refer to as . Clearly, the difference can be made arbitrarily small by including more and more bath sites. The main point is that the number of bath sites which we will employ to solve (5) and hence the Hilbert space is much smaller than for the case of conventional ED. Since in general is non-Gaussian, Wick’s theorem is not directly applicable. We formulate a perturbation theory in by introducing new fermionic degrees of freedom in the path integral for the partition function by means of the following identity:

(6) |

A similar approach has been used to derive a strong coupling expansion for the Hubbard modelPairault et al. (2000). In the above equation denotes the impurity Green function w.r.t. the discrete hybridization , which can be calculated exactly. With this substitution, the action becomes

(7) | |||||

From the expression for the partition function , we now integrate out and by expanding the exponential and using the fact that the functional integral over produces correlation functions that can be obtained exactly from the ED, such as the two-particle Green function:

(8) |

where we have used the shorthand notation . A compact expression for the Lehmann representation of the two-particle Green function is given in appendix A. Here we carry out this expansion up to fourth order (note that odd terms drop out of the expansion), with the result

(9) |

Here is the two-particle vertex constructed from the two-particle Green function as (henceforth we omit the superscript “” on and )

(10) |

with .

The auxiliary Green function is given by

(11) |

This function has some remarkable properties. Let us discuss these for the case and Hubbard interaction . Then we have and we can approximate for strong hybridization by . In the weak coupling limit, i.e. , approaches the bare Green function and , so that the dual perturbation theory becomes equivalent to conventional perturbation theory. On the other hand, for small, itself is small and can obviously be approximated as . This generates a fast converging strong coupling perturbation expansion around the atomic limit. In fact, one can show that in this case the expression for Green’s function obtained by a hybridization expansion of the Green function given in Ref. Dai et al., 2005 is contained in the lowest order diagram of our expansion. To see this, we write the auxiliary Green function as , where the self-energy correction of diagram a) in Fig. 1 is given by . In order to compare this to known results, we have to relate the auxiliary Green function to the Green function of the impurity, . The fact that we have used an exact identity to introduce the auxiliary fermions allows us to establish an exact relation between and , i.e.

(12) |

(see Ref. Rubtsov et al., 2008).

Inserting the approximation to into this equation yields the following expression for after some straightforward algebra:

where we have used that fact that the vertex is expressed in terms of the two-particle Green function by Eqn. 10. Considering the first term in the expression for yields a contribution . In order to compare with Ref. Dai et al., 2005, we take , as above. In this case, corresponds to the Green function for the atomic limit. For small we have and . Gathering the results we can approximate in the limit of small as

(14) |

which is essentially Eqn. (13) of Ref. Dai et al., 2005.

Hence the dual perturbation theory has the correct limiting behavior in both and weak and strong coupling limits. It therefore interpolates between these limits, so that sensible results can be anticipated even in the intermediate coupling regime. In the following we will demonstrate that this indeed is the case. For intermediate coupling, we further exploit the possibility to improve the starting point of the perturbation theory by expanding around the ED solution for a finite number of bath sites. In this case, low energy Kondo physics and the high energy incoherent features are correctly reproduced.

## Iii Results

The results shown in the following were performed using the diagrams a) to c) depicted in Fig. 1. If not otherwise stated, the results were obtained by making use of the Dyson iterations illustrated in the same figure: The self-energy was calculated using the bare auxiliary Green function, Eqn. 11, on the first iteration. Inserting the self-energy into the Dyson equation yields a new Green function which is subsequently used in the diagrams on the next iteration. This procedure is carried out until self-consistency and converges in typically less than ten iterations.

In the following we will discuss results obtained for the case of Hubbard interaction . In order to test our approach, we performed calculations for up to bath sites. The calculations were done for a semielliptical density of states of bandwidth , with . We take the half-bandwidth as the energy unit: .

In Fig. 2 we present results for , obtained for different number of bath sites up to . The quality of the representation of by is shown in the inset for . In order to determine the parameters and , we have minimized the distance function, Eqn. 4 for . This choice enhances the importance of the lowest frequencies in the minimization procedure. As can be seen in the inset, this condition results in and being equal on the first Matsubara frequencies. This turns out to be a good starting point for the perturbation theory. Using a smaller leads to a which better represents the tail of the hybridization function and generally leads to worse results.

The superperturbation results are compared to those of a numerically exact continuous-time quantum Monte-Carlo (CTQMC) calculation. One can see that while the expansion around the atomic limit lacks accuracy, a drastic improvement occurs for a single bath site, although the difference between and is still significant. For , the results are essentially converged. We find qualitatively the same behavior for a wide-range of values. The fast convergence w.r.t the number of bath sites significantly reduces the computational effort compared to CTQMC calculations.

In the figure we have plotted the results obtained from summing up skeleton diagrams. The use of skeleton diagrams is theoretically relevant since in this case the result is conserving in the Baym-Kadanoff senseBaym and Kadanoff (1961). The results obtained from the first Dyson iteration or from the lowest-order approximation, i.e. , however, achieve the same quality of approximation (not shown here).

Let us now investigate the role of the different diagrams in the perturbation expansion. To this end, we have plotted the results for the same parameters as in the previous figure, for different combinations of the diagrams shown in Fig. 3. We have used only a single bath site, for which the exact solution (open circles) and the initial solution obtained by ED (filled circles) differ substantially. One can clearly see that diagram a) yields by far the largest correction to the initial solution. The value on the first Matsubara frequency is almost exactly reproduced. This might be expected, since is identical to on the first Matsubara frequency. The largest deviations occur for the second and third frequency. We have zoomed into this region in the inset. One can see that all diagrams give a correction in the right direction, whereby the correction by diagram c) is negligible.

From Fig. 4 it is obvious that also spectral properties are correctly reproduced. We show the maximum entropy density of statesJarrell and Gubernatis (1996) (DOS) obtained from the imaginary time data. The analytical continuation of the quantum Monte-Carlo data (dashed-dotted line) exhibits the two Hubbard bands at and shows the Kondo resonance at the Fermi level. We cannot reproduce the Kondo physics by perturbing around the atomic limit (, solid line), in accordance with the findings in Ref. Dai et al., 2005. However, perturbation around the ED solution for a single bath-site already captures the Kondo resonance and yields good agreement compared to the exact solution.

In order to demonstrate that the approach also works for low temperatures, we present results for in Fig. 5. Although the expansion around the solution for a single bath site () shows small deviations in the imaginary time Green function , the approximation appears insufficient as seen in the density of states. The superperturbation around the two bath-site solution however is almost exact.

As already mentioned, the formalism of introducing auxiliary degrees of freedom in the path integral representation presented here has originally been introduced in the context of lattice models, termed “dual fermion approach”Rubtsov et al. (2008). In that case the action is decomposed into an impurity part with hybridization and a term which contains the bare dispersion of the lattice fermions. Since this decomposition is independent of the form of the hybridization, it can be represented by the discrete hybridization . Hence the present framework can be used to combine dual fermion calculations with ED for the efficient solution of lattice models.

## Iv Conclusions

In conclusion, we have presented an efficient approximate solver for the quantum impurity problem. The formalism straightforwardly generalizes to the multiorbital case and is suitable for a general form of the interaction. The perturbation expansion has been shown to be convergent in both weak and strong coupling limits. It was further demonstrated that the solver also works well in the intermediate coupling regime and down to low temperatures. The approximation is controlled in the sense that the final result can be judged on the basis of how well the input hybridization is approximated by a finite number of bath sites. The method can be used for the study of multiplet effects in solids in the context of realistic LDA+DMFT calculations.

###### Acknowledgements.

We would like to thank Erik Koch for valuable discussions. This work was supported by DFG Grant No. SFB 668-A3 (Germany), the Leading scientific schools program and the “Dynasty” foundation (Russia).## Appendix A Lehmann representation of the two-particle Green function

In this appendix, we derive a compact expression for the Lehmann representation of the two-particle Green function (2PGF). A similar expression was given in Ref. Toschi et al. (2007) without explicit consideration of the singular contributions. By definition, the 2PGF in Matsubara space is given by

(15) | |||||

Here time translation invariance of the imaginary time 2PGF has been used. Note that here the frequencies in the exponential corresponding to annihilation and creation operators have the same sign in contrast to the usual definition for the Fourier transform. Correspondingly, energy conservation requires . By restricting the range of integration such that time ordering is explicit, one obtains different terms. These can be brought into the same form by permuting the operators and corresponding frequencies. By the anticommutation relations, each term picks up the sign of the permutation. After introducing the sum over eigenstates, the 2PGF can be written as

where the first sum is over the eigenstates and the second over all permutations of the indices . We further have defined , and and e.g. denotes the permutation of the first index. Here the different choice of convention for the Fourier transform simplifies the notation since otherwise the sign of the frequency associated with the creation operator would have to be permuted. The function is given by the integral

(17) | |||||

The latter expression can be evaluated by taking care of the delta functions that arise from equal energies, with the final result

(18) |

## Appendix B Two-particle vertex in the atomic limit

For the calculations without bath-sites, we have used an explicit expression for the two-particle vertex in the atomic limit. This can be obtained from result of Appendix A together with the definition of the vertex, Eqn. 10. Here we reintroduce explicitly the dependence on , which yields a more symmetric form of the result. We also switch back to the usual convention for the Fourier transform, for which the energy conservation reads .

Using that the eigenenergies of for the impurity states ,,, are given by , respectively, after some algebra the vertex is obtained as

(19) |

## References

- Hewson (1993) A. C. Hewson, The Kondo Problem to Heavy Fermions (Cambridge Univ. Press, Cambridge, 1993).
- Savkin et al. (2005) V. V. Savkin, A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. Lett. 94, 026402 (pages 4) (2005).
- Kotliar et al. (2006) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (pages 87) (2006).
- Katsnelson et al. (2008) M. I. Katsnelson, V. Y. Irkhin, L. Chioncel, A. I. Lichtenstein, and R. A. de Groot, Rev. Mod. Phys. 80, 315 (pages 64) (2008).
- Dai et al. (2005) X. Dai, K. Haule, and G. Kotliar, Phys. Rev. B 72, 045111 (pages 6) (2005).
- Rubtsov et al. (2005) A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (pages 9) (2005).
- Werner et al. (2006) P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (pages 4) (2006).
- Haule (2007) K. Haule, Phys. Rev. B 75, 155113 (pages 12) (2007).
- Rubtsov et al. (2008) A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 77, 033101 (pages 4) (2008).
- Hafermann et al. (2007) H. Hafermann, S. Brener, A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, JETP Lett. 86, 677 (2007).
- Brener et al. (2008) S. Brener, H. Hafermann, A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 77, 195105 (pages 12) (2008).
- Toschi et al. (2007) A. Toschi, A. A. Katanin, and K. Held, Phys. Rev. B 75, 045118 (pages 8) (2007).
- Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- Koch et al. (2008) E. Koch, G. Sangiovanni, and O. Gunnarsson (2008), arXiv:0804.3320.
- Pairault et al. (2000) S. Pairault, D. Sénéchal, and A.-M. S. Tremblay, Eur. Phys. J. B 16, 85 (2000).
- Baym and Kadanoff (1961) G. Baym and L. P. Kadanoff, Phys. Rev. 124, 287 (1961).
- Jarrell and Gubernatis (1996) M. Jarrell and J. E. Gubernatis, Phys. Rep. 269, 133 (1996).