Exact diagonalization as an impurity solver in dynamical mean field theory

The dynamical mean-field theory (DMFT) maps a correlated lattice problem onto an impurity problem of a single correlated site coupled to an uncorrelated bath. Most implementations solve the DMFT equations using quantum Monte-Carlo sampling on the imaginary time and frequency (Matsubara) axis. We will here review alternative methods using exact diagonalization, i.e., representing the many-body ground state of the impurity as a sum over Slater determinants and calculating Green’s functions using iterative Lanczos procedures. The advantage being that these methods have no sign problem, can handle involved multi-orbital Hamiltonians (low crystal symmetry, spin-orbit coupling) and – when working completely on the real axis – do not need a mathematically ill-posed analytical continuation. The disadvantage of traditional implementations of exact diagonalization has been the exponential scaling of the calculation problem as a function of number of bath discretization points. In the last part we will review how recent advances in exact diagonalization can evade the exponential barrier thereby increasing the number of bath discretization points to reach the thermodynamic limit.


Introduction
The DMFT equations map a many-body problem on an infinite lattice to a local impurity problem. The potentials of the uncorrelated bath are chosen such that a self-consistency condition for the local one particle Green's function is fulfilled. In spite of being a great simplification, the resulting many-body impurity theory is still not trivial to solve. Finding solutions to the correlated Anderson impurity model has a long history including calculations of the Kondo problem and X-ray edge singularities [1][2][3][4][5][6][7][8]. With the development of DMFT, the numerical challenges encounter in solving the many-body impurity problem gained more attention and great progress has been made. One of the most successful methods is the continuous-time diagrammatic quantum Monte-Carlo impurity solver [9]. While this method is among the most efficient methods available for the single-band Hubbard model, difficulties can arise when facing multi-band Hamiltonians with low-symmetry interactions, rotational invariant Coulomb repulsion or spin-orbit coupling that lead to inherently non-diagonal elements in the Green's function. An active field of research continuously pushes further the limits of what is possible to calculate [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] but more work needs to be done before one can reach the chemical accuracy needed for many applications. A further complication of current quantum Monte-Carlo methods is the ill-posed analytical continuation needed to transform the calculations on the Matsubara axis to the real-frequency axis. This transformation sets a limit to the amount of detail one can derive about the spectral shape. Exact diagonalization has long been employed as a method to solve the Anderson impurity model [1][2][3][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]. One of its major disadvantages is the need to discretize the bath Green's function. The exponential scaling of storage and computation time in the number of discretization points for most exact-diagonalization impurity solvers only allows for the inclusion of a few bath sites. The great advantages of the exact diagonalization, on the other hand, are that it can handle any form of Hamiltonian (there is no sign problem) and calculate the Green's functions directly on the real frequency axis -thereby avoiding the analytical continuation. In addition, it can straightforwardly calculate different types of spectroscopy and response functions.
Recent advances in numerical methods now allow for the inclusion of many more bath orbitals, thereby overcoming many of the disadvantages the early implementations of this method had. In this paper we will first review literatures available on different exact-diagonalization schemes, starting with the formulation of DMFT in a tight-binding language. We further discuss the problems that arise when representing involved Green's functions on the Matsubara axis. We then continue with the method introduced by Lu et al. [30] in order to find a stable self-consistency cycle that functions fully on the real-frequency axis which also allows one to increase the number of bath sites such that (quasi-)continuous spectral functions can be computed. In the last section we show examples of different spectral functions one can calculate using this method.

DMFT in a tight-binding language
We start by defining a Hamiltonian on an infinite lattice. Each lattice site i can have several spin-orbitals which we label with τ ranging from 1 to N τ . There are only one-particle interactions between different lattice sites. On each site we include oneand two-particle interactions. The resulting lattice Hamiltonian is: The parameters include the onsite energy, the crystal-field splitting between the orbitals, the hopping between sites as well as spin-orbit coupling and possible magnetic fields. The parameters U include the orbital and spin dependent onsite Coulomb repulsion. Although the Coulomb repulsion in its bare form decays as 1/r as a function of distance r and is thus long ranged, the locality of the Coulomb interaction can be justified by the strong screening one generally expects to take place in solids. The DMFT approximates the full lattice Hamiltonian by an impurity model where only one site i = 1 remains fully correlated. At all other sites the Coulomb interaction is replaced by a site independent local self-energy Σ τ1,τ2 (ω). This approximation allows one to take the self-energy calculated from the impurity model to be the selfenergy of the lattice model. The self-energy is ω dependent and thus cannot be entered as a mere potential change in the tight-binding Hamiltonian. Nonetheless it is possible to add energy dependent interactions into a tight-binding description by adding Fig. 1. Graphical representation of the DMFT approximation. A full locally interacting lattice model (left pannel) is mapped to an impurity model with one correlated site and all other sites approximated by a local self-energy (middle pannel). The local self-energy is chosen such that a full lattice model with all sites approximated by a local-self energy (right pannel) give the same local Green's function as the impurity site. In blue we depict the one particle tight binding Hamiltonian defined by the parameters i 1 ,i 2 ,τ 1 ,τ 2 in equations (1) and (4). In red we depict the local Coulomb interaction defined by the parameters U τ 1 ,τ 2 ,τ 3 ,τ 4 in equations (1) and (4). In Green we depict the self energy defined by the parameters α j and β j in equation (2) and equation (4).
additional auxiliary lattice sites. One way to determine the auxiliary lattice sites and their Hamiltonian is by writing the self-energy as a continued fraction: with α j and β j matrices of the dimension N τ by N τ whose numerical values determine the actual shape of the self energy. The index j runs from 0 to N j where N j is the number of discretization points used to represent the energy dependence of Σ τ1,τ2 (ω).
In this tridiagonal form one can create a tight-binding Hamiltonian acting on auxiliary sites whose one-particle Green's function reproduces equation (2). This auxiliary Hamiltonian has 1 site that represents the site where the self energy is added and N j auxiliary sites with j running from 1 to N j . α 0 defines an additional onsite potential that has to be added to each real site where a self energy is added. The matrices α h define the onsite energy of the auxiliary site j. The matrices β j define the hopping between site j and j + 1. For each site where a self-energy needs to be included one can now add the auxiliary Hamiltonian and sites to the single particle part of the original tight-binding Hamiltonian. The impurity Hamiltonian to solve thus becomes: with the operators b † i,j,τ (b i,j,τ ) creating (annihilating) a particle in the j-th auxiliary site representing the self energy acting on site i with local spin-orbital τ . The parameters and U used in the first line of equation (4) are equivalent to those from the original lattice Hamiltonian as given in equation (1). The parameters α and β are those of the tridiagonal representation of the self-energy as given in equation (2).
In Figure 1 we show the DMFT approximation graphically following the review of Held [43]. On the left we depict the original lattice model, where each site, represented by a square, has full local Coulomb interactions. The middle panel shows the DMFT approximation whereby only on one site the full interactions remain. On all other sites the interactions are replaced by a self-energy. The self-consistent condition on the self-energy is such that the impurity Green's function must be the same as the Green's function obtained from the full lattice using the impurity self-energy.
We now can formulate the DMFT self consistency loop (see also Koch et al. [44]).
(1) Starting from an approximate local impurity self energy Σ i one can calculate the full lattice Green's function G l using the local self-energy and the noninteracting part of the lattice Hamiltonian (right panel of Fig. 1) (2) The lattice Green's function together with the local self-energy define the bath Green's G b function as: wich combines the non-interacting part of the Hamiltonian as well as the selfenergy of the Anderson impurity problem as depicted in the middle panel of Figure 1. (3) Representing the bath Green's function as: with A i and B i Hermitian matrices of dimension N τ by N τ and N b + 1 the number of discretisation points used to represent the energy dependence of G b (ω), one can write the Anderson impurity Hamiltonian as: How to write any Green's function in the form of equation (6) can be found in appendix B of the work by Lu et al. [30]. The matrices A j and B j are not equivalent to, but related, by a near unitary transformation, to the matrices j , α j and β j of equation (4). (4) Exact diagonalization methods without basis set optimisations used to equation (7) scale exponentially in terms of memory storage and computation time needed in the number of bath sites used (N b ). Traditionally only a few bath sites are included and the number of bath sites is reduced by fitting a small number of parameters representing the bath Green's function on the imaginary axis. See for example Koch et al. [44] or Georges et al. [26]. Alternative reduction schemes have been proposed, such as truncated continues fractions by Si et al. [45,46]. Both methods work well for simple models, but have problems with numerical stability and general applicability once the Hamiltonian and Green's functions become more complex. In Section 4 we will review progress in this respect made in the last years. (5) Once the Anderson impurity model has been defined, one needs to solve it and calculate the local impurity Green's function G i which is in general an N τ by N τ matrix. This can be done using the Lanczos routines as for example explained in Appendix C of Lu et al. [30]. These calculations consist of two steps. In the first step the lowest energy eigen-state of the impurity model is calculated, in the second step the Green's function is calculated starting from the ground-state. (6) After the local impurity Green's function has been obtained, a new self-energy needs to be calculated, (7) Restart at step (1) with the new impurity self-energy until convergence is reached.

Numerical challenges in DMFT
There are two parts in the algorithm presented in the previous section that pose numerical challenges. Mostly discussed in the literature is that the computation time as well as memory cost of the solution to the Anderson impurity problem scales exponentially in the number of bath orbitals included. Traditional exact-diagonalization methods can handle about 16 total basis states (sites times orbitals), which translates to about 15 bath sites for a single band model, 7 bath sites per orbital for a two band model and 4 bath sites per orbital for a three band model. Recalling that the number of bath sites is equal to the number of discretisation points in energy used to represent the Green's functions it is clear that this number is far from enough to create continuous spectral functions. The other numerical challenge for the DMFT algorithm lies in the calculations during the self-consistency loop. The calculations of the self-energy and bath Green's function as done in step (2) and (6) are numerically unstable, subject to severe loss of significance at large energies and can lead to unphysical (non-causal) representations of the Green's function or self energy.

Loss of significance in calculating the self energy
In step (6) in the algorithm outlined above, the self energy is calculated as For large values of |ω| both G i and G b go to zero and their inversions thus become infinite while the self-energy, i.e. their difference remains finite. This poses a problem once implemented in a computer where only finite precision math is available (about 16 digits in most calculations). The direct calculation of the self-energy is numerical unstable and leads to loss of significance. Furthermore due to the discretisation of both G b and G i at possible different energy meshes, one might find that for some energy points the imaginary part of Σ i becomes positive thereby breaking causality. Traditionally, this problem is solved by calculating the self energy on the imaginary axis, where the later instability is absent and the first is solved by analytical relations between lim ω→∞ G(iω) and moments of the Green's function [47]. The use of Green's functions on the Matsubara axis requires one to transfer back from the imaginary axis to the real axis at some point in the calculation which leads to an ill-posed inversion problem. This problem will be discussed in more detail in Section 3. In Section 5, we will review the recent progress made to do numerically stabel DMFT self-consistency loops on the real frequency axis.

Anderson impurity problem
In step (5) of the algorithm, an Anderson impurity problem needs to be solved. If no further approximations are made, the Hilbert space scales exponentially with the number of bath sites and one is restricted by a total of bath plus impurity spinorbitals of about 16. This puts a severe limitation on the accuracy one can reach using traditional exact-diagonalization methods. It is however to be expected that not all elements of the Hilbert space (determinants) will be important. Based on this notion one can device further approximations that allow one to include many more bath orbitals. These methods applied to impurity problems go back to the work of Gunnarsson and Schönhammer in the mid 80's [1][2][3]. In Section 4 we will review recent progress made that allows one to include several hundreds of bath orbitals in exact-diagonalization routines.
3 Matsubara frequencies and the ill posed fitting to functions on the real axis Exact-diagonalization methods require that the bath Green's function is represented on a finite mesh. Traditionally this is done by optimizing the parameters of the discretized Green's function by comparing the discretized Green's function to the original one in the imaginary-frequency domain. See for example Section VI.A.2 of the review of Georges et al. [26]. As noted by Koch et al. [44], details of the fitting procedure can be important and are accordingly discussed in the literature [48][49][50][51]. Before we review how stable self-consistency DMFT algorithms can be implemented on the real-frequency axis, including solving the loss-of-significance problem inherent to the calculation of a self-energy from two Green's functions, we briefly review the mathematics behind Green's functions on the Matsubara axis and the corresponding transformations.
The relation between the fermionic Matsubara Green's function and the Green's function on the real frequency axis is given as: with ω the energy, β = 1/(k B T ) the inverse temperature, τ ∈ [0, β] the imaginary time, and A(ω) = −(1/π)Im[G(ω)] the spectral density. In order to work with exactdiagonalization methods on the real frequency axis we represent our Green's functions as a sum over delta functions: with β i and α i the parameters that determine the Green's function. Note that η → 0 + is an infinitesimal positive number and should not be confused with a finite broadening, i.e. the Green's function is given as a true sum of delta functions. For a single delta function with α i = ω 0 and β i = 1 the transformation to Matsubara frequencies yields: The transformation from G(ω) to G(τ ) is bijective, i.e., for each delta function at energy α i there is a unique function G αi (τ ). However, the transformation does not conserve norm or orthogonality of the transformed functions. The set of delta functions on the real frequency axis is complete and orthonormal ( δ(ω , α i )δ(ω, α j )dω = δ(i, j)), while the transformed functions on the Matsubara axis are not. In order to study the implications of this loss of orthonormality in more detail we plot in Figure 2 the function G(τ ) as a function of ω 0 , the energy position of the delta function. We separate the total Green's function in a symmetric and antisymmetric part: The symmetric Green's function (G s (τ, ω 0 )) is the transform of the sum of two delta functions at ±ω 0 , the antisymmetric Green's function is the transform of the difference between two delta functions at ±ω 0 : As shown in Figure 2, the Green's function of a delta function when transformed to the imaginary axis varies only slowly as a function of the position of the delta function. The transformation is ill posed in the sense that the detail and information content in G(ω), represented by different values of G(ω) at each value of ω, is only represented in G(τ ) if this function is known to a high precision. One of the implications of representing the Green's function on the imaginary axis becomes clear when one considers the derivative of G(τ, ω 0 ) with respect to ω 0 (Fig. 3). This derivative illustrates how sensitive the error in the position of the pole depend on the numerical accuracy with which the Green's function on the imaginary axis is known. As can be seen in Figure 3, for large values of ω 0 , dG(τ, ω 0 )/dω 0 goes to zero. The Matsubara representation of the Green's function is thus largely insensitive to the exact position of the pole at large energies. More important, for many calculations is the observation that the derivative for the symmetric Green's function goes to zero for ω 0 = 0. The Matsubara representation of the Green's function is not sensitive to the exact spectral weight at the chemical potential. One can always open a small gap in the Green's function without numerically changing the Matsubara representation of the Green's function by a significant amount.
Most self-consistency loops used in exact-diagonalization schemes need to fit a set of poles on the real-frequency axis to a known Green's function on the Matsubara axis. Knowing that this transformation is ill conditioned raises the question how many poles on the real-frequency axis one can fit to a given function known on the Matsubara axis within certain numerical accuracy. Naturally one can always fit more poles, but at some point the weight and position of additional poles become ill defined or arbitrary. Standard methods to regularize an ill-posed inversion problem use singular value decomposition techniques and only keep those eigen values larger than the known error in the original data. The simple form of G(τ, ω 0 ) for the transformation of a delta function allows for an analytical expresion. One can expand G(τ, ω 0 ) on a complete set of Legendre polynomials, which yields analytical prefactors for the expansion coefficients as a function of ω 0 , i.e. the position of the delta function on the realfrequency axis. In Figure 4 we show the first 22 expansion coefficients as a function of the position of the pole on the real frequency axis. For the symmetric part of the Green's function (even expansion coefficients) we notice that all expansion coefficients are positive due to causality and they decay exponentially with increasing order of the Legendre polynomial.
For a given numerical accuracy with which G(τ ) is known, one can now define the number of poles n on the real frequency axis that represents G(τ ) with the given accuracy. All expansion coefficients smaller than the accuracy of the Green's function for all values of ω 0 are not important in determining the Green's function. For an accuracy of 10 −4 and a maximum energy of the pole at 5/β the even expansion coefficients α 0 to α 6 are large enough to be determined. This leaves 4 numbers, which can be uniquely created by linearly combining 4 poles at different energy. In fact there might be solutions with less poles, but there will at least be one solution with 4 poles. Naturally there will be infinitely many different solutions with more than 4 poles.
In Figure 5 we show the number of poles that completely determine the Green's function on the Matsubara axis for different maximal pole energies. We here assume that the Green's function on the Matsubara axis is known upto 10 digits accuracy. As one can see the Green's functions can be represented by a relatively small number of poles. From ligand-field theory or finite-size cluster calculations of valence band photo emission spectroscopy for transition metal compounds, we know that several hundred poles are needed to converge the multiplets and charge-transfer excitations, which can further increase to a thousand poles for spectral functions of rare earth compounds. This kind of detail is not easily represented using Matsubara frequencies unless one is able to use arbitrary precision math on the imaginary axis. As the later is not suited for heavy numerical work we will review methods that discard completely the transformation to the imaginary axis. This then does pose the question on how to calculate the self-energy in a numerically stable fashion avoiding the problem of loss of significance as outlined above.

DMFT on the real frequency axis
To prevent the numerical problems inherent to the Matsubara representation of involved spectral functions, Lu et al. [30] devised a scheme fully on the real-frequency axis, which can be seen as an extension of the work on a Bethe lattice by Si et al. [45]. All Green's functions and self-energies are represented as a sum over poles that reside exactly on the real axis: with a j and b j matrices of dimension N τ and a j proportional to the identity matrix. Summing Green's functions in this notation is trivial. Inverting the Green's function, which is needed in order to calculate a self-energy can be done by rewriting the Green's function as a resolvent of a Matrix. Following appendix B of the work of Lu et al. [30] one can rewrite the Green's function as: with α j and β j matrices of dimension N τ determined by the matrices a j and b j . The algorithm involves once a block Lanczos tri-diagonalization routine and once a diagonalization of a block tri-diagonal matrix, both nummerical stable algorithms, it's inverse is then easily calculated as: In order to reduce the number of poles one can combine two closely-located poles into one, such that the first moments of the spectral function are retained. Appendix F of the work of Lu et al. [30] shows the numerical detail on how this can be implemented. Several different choices can be made, whereby one can choose for a fixed spaced mesh, a mesh that has more dense points close to the Fermi energy or optimize both the mesh spacing as well as the pole intensity at the same time. The optimal choice will depend on the quantity of interest. The resulting set of equations allow one to do the full DMFT self-consistency loop directly on the real frequency axis.

Enlarging the number of bath sites
The remaining challenge to tackle is how to solve the Anderson impurity model with several hundred of bath sites without running into the problem of an exponentially large Hilbert space. Here one can use ideas from quantum chemistry and renormalization group theory that allow one to optimize the basis set such that at each step of the calculation only a finite number of determinants is important. Using standard quantum chemistry methods whereby the bath orbitals are transformed to their natural impurity basis, one can converge the ground-state with only a few (about 10) bath orbitals per spin-orbital [28,29,52]. The convergence of the Green's function, which plays a central part in DMFT, requires a much larger basis, i.e., the optimized basis for the ground-state does not represent the excited states well.
Within a Lanczos routine the spectral function is represented on a Krylov basis defined by the states that span the Hilbert space given by ψ n = Hψ n−1 with ψ 1 = T ψ 0 . Here ψ 0 is the ground-state, T is an operator that either removes or adds an electron to an orbital τ , and H is the interacting Hamiltonian. For an optimal representation of the excited states, the basis set should be optimized at each step during this iterative procedure. This allows one not only to converge the ground-state, but also to give a good representation of the spectral function.
The result is an algorithm that can handle several hundred of bath orbitals with all the advantages of exact-diagonalization routines. In Figure 6 we show the DMFT spectral function of the Hubbard model for a semi-circular G 0 (ω) as a function of the value of U and the number of included bath orbitals. With small broadening (applied after the calculation is finished), one can observe discrete peaks in the Green's function for a small number of bath sites, which become quasi-continuous when the number of bath sites is increased. Once a larger broadening is applied, one finds that the Green's functions with more than 30 bath orbitals are reasonably converged. One should be very careful though with calculations that do not have enough bath orbitals. The critical value of U separating metallic from insulating solutions changes with the number of bath orbitals included. This is understandable, as the bath Green's function is discretized and the impurity actually always hybridizes with a gaped and thus insulating state. Once there is a gap at the Fermi energy it becomes easier for correlations to widen the gap and thus create an insulator.

Spectroscopy in DMFT
One of the great advantages of real-frequency DMFT using a Lanczos-based solver is that once a self-consistent solution is obtained, several types of spectroscopy can be directly calculated. The spectral function (I) of any linear-response spectroscopy is given as: 2560 The European Physical Journal Special Topics  and with ψ 0 the Anderson impurity ground state, H the Anderson impurity Hamiltonian and T an operator describing the excitation made during the experiment. This spectrum can be calculated in the same fashion as the Green's function. The only restriction (inherent to the single-site approximation we discuss here) is that T can only act on the impurity site, therefor one can calculate, for example, the total magnon density of states but not its dispersion without going beyond the single-site approximation.
For neutron scattering, one would replace the operator T by the spin operator acting on the impurity site, which allows one to calculate both the total spin density of states as well as orbital excitations in e.g., rare-earth compounds. Several forms of core level spectroscopy can be calculated, including core level photo emmission (T annihilates a core electron and H includes the core-valence interaction explicitly), X-ray absorption spectroscopy (T annihilates a core electron and creates a valence electrons and again H includes the core-valence interaction explicitly), or inelastic X-ray scattering where T has the same form as in X-ray absorption spectroscopy with different orbital-dependent matrix elements. Higher order response functions such as resonant raman spectroscopy or the emerging technique of resonant inelastic X-ray scattering can also be calculated by adding an additional interaction: where T 1 and T 2 describe the excitation and de-excitation process in a resonant Raman experiment. ω and Ω are the incident photon energy and Raman energy loss, respectively.
The great advantage of DMFT is that it can capture band excitations as well as excitons and the hard-to-describe intermediate regime of resonances. In spectroscopy one can divide excitations into three different categories depending on the interaction strength between the electron-hole pair created by promoting an electron in the  occupied bands into the empty ones. If this interaction is weak, one observes basically a convolution of the band structure of the occupied and unoccupied bands. If the interaction is very strong a bound state is formed and one observes a single excitonic peak outside the continuum. Often this peak has an internal structure due to multiplet effects. In the intermediate regime, one observes resonances with typical asymmetric line-shapes that are notoriously hard to describe theoretically. All these types of excitations are found in different materials with important properties resulting from the interplay between local and itinerant states.
In Figure 7 we show an example of the transition from band excitations to excitions via resonances with the variation of U and Q. With U being the Coulomb repulsion within the Hubbard band and Q the Coulomb repulsion between a deeper dispersionless occupied band and the valence band. We show both the absorption spectra (the electron is transfered from the core level into the valence band) as well as the photo emission spectra (the electron is removed from the core level). More details on the description of excitons, resonances and band transitions using DMFT can be found in the original paper by Haverkort et al. [53].
An example of a real material and experiment where both excitons and resonances are observed is the case of 3p to 3d excitations as observed in non-resonant inelastic X-ray scattering of MnO. At low energy transfer one finds a few strongly bound excitons, whereas at higher energy transfer one finds a broad continuum. The excitonic bound states are well modeled using a ligand-field model, in which the solid is approximated by a single atom interacting with it nearest neighbor ligand atoms. In the case of MnO this equivalents to a single Mn d shell interacting with six O 2p shells of the neighbors. Slight deviations may occur as, for example, one needs to include the Mn 4s orbital into the tails of O 2p orbitals. While such approximations generally work very well for excitonic bound states in transition metal compounds [56], they can not describe band transitions nor resonances. As shown in Figure 8, the resonance at hihger transfer energy in MnO is poorly described in ligand field theory. In DMFT, however, such resonances can be captured and the transition into the resonance between 50 to 60 eV is well described in addition to the excitons.
To summarize, we have presented how DMFT self-consistency loops can be fully implemented on the real-frequency axis using an exact-diagonalization impurity solver, and demonstrated that it is a method well suited to describe several types of spectroscopy on correlated transition metal compounds with a reasonable computational effort.
We gratefully acknowledge support of the Deutsche Forschungsgemeinschaft through FOR 1346. Open access funding provided by Max Planck Society.
Open Access This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.