# Hierarchical equations of motion for impurity solver in dynamical mean-field theory

###### Abstract

A nonperturbative quantum impurity solver is proposed based on a formally exact hierarchical equations of motion (HEOM) formalism for open quantum systems. It leads to quantitatively accurate evaluation of physical properties of strongly correlated electronic systems, in the framework of dynamical mean-field theory (DMFT). The HEOM method is also numerically convenient to achieve the same level of accuracy as that using the state-of-the-art numerical renormalization group impurity solver at finite temperatures. The practicality of the novel HEOM+DMFT method is demonstrated by its applications to the Hubbard models with Bethe and hypercubic lattice structures. We investigate the metal-insulator transition phenomena, and address the effects of temperature on the properties of strongly correlated lattice systems.

###### pacs:

71.10.Fd, 71.27.+a, 71.30.+h## I Introduction

Strongly correlated electronic systems (SCS) generally refer to materials with d- or f-electrons. These localized electrons strongly interact with each other and with the surrounding itinerant electrons, leading to various intriguing many-body phenomena, such as Kondo phenomena, metal-insulator transition in transition metal oxides, and heavy fermion in rare-earth compounds. SCS quickly becomes an active research area in condensed matter physics.Ful95 ()

Exact treatment of SCS for classic band theory and model methods is difficult because of the nonperturbative many-body nature and the complexity of SCS materials. Since the proposal of investigating SCS in the limit of infinite dimension in 1989,Met89324 () the dynamical mean-field theory (DMFT)Jar92168 (); Geo926479 (); Geo9613 () has been rapidly developed and applied, leading to a dramatic progress in understanding the properties of SCS. In this limit, the self-energy is local, so that the spatial fluctuations in strongly correlated systems can be neglected while only the on-site Coulomb interaction is taken into consideration. As a result, the many-body lattice problem is mapped to an effective impurity problem, in which the correlated electrons on the impurity site interacts with a frequency dependent mean field. This mean field is then updated in a self-consistent way via the solution of the impurity model in order to obtain the self-energy of the lattice Green’s function.

The DMFT method was firstly used with Hubbard model, and revealed the three-peak spectral function induced by the strong electron-electron correlation.Jar92168 (); Geo926479 () These three peaks are comprised of one central quasi-particle peak around the Fermi energy and two single-particle Hubbard peaks. In some transition metal oxides, the central quasi-particle peak gradually transfers its weight to the Hubbard peaks with increasing Coulomb interaction,Zha931666 () finally leading to metal-insulator transition, which is usually referred to as Mott transition.Geb97 () It is one of the most fascinating phenomena of SCS that has been observed in various transition metal oxides,Ima981039 () and extensively investigated by DMFT methods.Geo937167 (); Roz9410181 (); Roz993498 (); Bul99136 (); Bul01045103 (); Dai05045111 (); Ono03035119 (); Flo02165111 () In recent years, the DMFT methods have been extended to the investigations of some complicated SCS with multi orbital, nonlocal Coulomb and exchange interaction, or taking into account spatial correlations.Si963391 (); Chi003678 (); Kot2001 (); Drc0561 (); Kus06054713 (); Tos07045118 () Especially, the combination of DMFT method with density functional theory makes it possible to simulate complicated real compounds with strong electron correlations.Ani977359 (); Kot06865 () The DMFT method has already been a powerful tool for the investigation of strongly correlated physics.

The accuracy and efficiency of DMFT calculation is determined by the key component, the impurity solver, that solves the effective impurity problem. A vast amount of theoretical efforts have been devoted to achieving this goal. The perturbative methods such as the iterated perturbation theoryKaj9616214 () and non-crossing approximationBic87845 () sum the perturbative series of interaction diagrams to different orders, resulting in different levels of accuracy. Extensions of these methods to away from half filling (iterated perturbation theory) or to very low temperature regime (non-crossing approximation) are technically complicated, which largely limit the applications of these methods. The currently used nonperturbative methods include the exact diagonalization approach,Dag94763 (); Caf941545 (); Si942761 () quantum Monte Carlo approach,Hir862521 (); Gul11349 () and numerical renormalization group (NRG) method.Wil75773 (); Bul08395 (); Hof001508 (); Wei07076402 (); Pet06245114 ()

Although the existing methods have been successfully used in characterizing the fundamental characteristics of SCS, they have their own limitations.And6141 (); Hub63238 (); Lee8699 () Both the exact diagonalization approach and the configuration interaction methods diagonalize the Hamiltonian in a finite dimensional Hilbert space. In DMFT applications, the continuous bath degrees of freedom are discretized and a small number of bath sites (orbitals) are chosen to represent the full bath. Consequently, the local density of states is composed of discrete peaks which are then broadened artificially. The energy resolution in the spectrum depends on the number of bath sites. Therefore, although the quasiparticle weight can be extracted from the Matsubara Green’s function, it is difficult to obtain the real frequency quantities,Ani2010 () including the low energy Kondo peak and the high energy Hubbard peaks; see for instance, Fig. S11 in the Supplemental Material of Ref. Li12266403, . The computational cost of some quantum Monte Carlo (QMC) approaches, including the Hirsch-Fye algorithm and the continuous time QMC, increases dramatically as the temperature decreases. Moreover, numerical methods such as the maximum entropy method are often required to extract real frequency spectral function from the imaginary time domain via some analytic continuation operation. This may introduce additional errors.Jar96133 () In recent years the NRG method has achieved significant progress, thanks to the improvement of its algorithm and advancement of computer hardware. However, some of its basic features have remained unchanged: It utilizes logarithmic discretization and truncation of energy spectrum during iterative diagonalization. Therefore, quantities produced by the NRG method are less accurate at high energy or high temperatures than at low energy or low temperatures.Bul08395 () Therefore, an accurate and efficient DMFT impurity solver for the investigation of the strong correlation effects at finite temperatures is highly desirable.

In this work we propose to use a hierarchical equations of motion (HEOM) approachJin08234703 () as the impurity solver of DMFT. The HEOM method treats quantum impurity systems from the perspective of open dissipative dynamics. As will be discussed in Sec. II.1, in principle the HEOM formalism is formally exact, as long as the bath environment satisfies Gaussian statistics, which is true for noninteracting electron reservoirs. The HEOM theory can be established based on the Feynman–Vernon path-integral formalism,Jin08234703 (); Fey63118 (); Kle09 (); Xu05041103 (); Xu07031107 () in which all the system-bath correlations are taken into consideration.

In practice, the HEOM method is implemented without invoking any approximation or any tricky controlling parameter. The numerical results are quantitatively accurate for a wide range of impurity systems, as long as they converge with respect to the increasing truncation level .

Usually, the HEOM results converge uniformly and rapidly with the increasing . A higher is needed to achieve numerical convergence at a lower temperature. Consequently, the computational cost increases substantially as the temperature lowers. Currently for a symmetric single impurity Anderson model, the HEOM approach reaches quantitatively accuracy for temperature (Kondo temperature) with the computational resources at our disposal.Li12266403 () In particular, the HEOM approach remains accurate at high temperature and large frequency (energy) range. For all the results presented in this work, is chosen to be sufficiently large so that numerical convergence is always guaranteed.

The HEOM approach is applicable to a general quantum impurity system. It is capable of treating strongly correlated impurities with more than one orbitals. The applicability of HEOM approach to various equilibrium and nonequilibrium properties of quantum dot systems has been demonstrated.Zhe08184112 (); Wan13035129 () These include the studies on dynamical Coulomb blockadeZhe08093016 () and dynamical Kondo transitions in quantum dots.Zhe09164708 (); Zhe13086601 () Previous studies have shown that the HEOM approach is capable of achieving the accuracy of the latest high-level NRG approach, as demonstrated by the impurity spectral functionsLi12266403 () and the local magnetic susceptibilityWan13035129 () of Anderson impurity model systems at finite temperatures. All these successful applications suggest that the HEOM approach is very suitable to be used as an impurity solver in the framework of DMFT method.

The remainder of this paper is organized as follows. The HEOM based impurity solver is introduced in Sec. II, together with its main features and current limitations. In Sec. III, the results for the Mott transition from DMFT using HEOM impurity solver are systematically compared with previous NRG results to show its accuracy and efficiency. Concluding remarks are finally given in Sec. IV.

## Ii Hierarchical equations of motion based DMFT impurity solver

### ii.1 A formally exact HEOM formalism for quantum impurity systems

The HEOM approach is based on quantum dissipation theory, which can be used to characterize the quantum impurity systems as open systems embedded in surrounding environments composed of itinerant electron reservoirs. The derivation of the HEOM formalism has been detailed in Refs. Jin08234703, , Zhe09164708, , and Zhe121129, . Here, we introduce some of its basic features.

The Hamiltonian of the quantum impurity system can be expressed as

(1) |

where represents the quantum impurity systems containing strong correlation interactions, is the noninteracting electron reservoirs, and accounts for the impurity-reservoir couplings. In these terms, and are the creation and annihilation operators for the impurity state ; and are those for the reservoir state of energy . The hybridization functions of the electron reservoir are defined as .

In the quantum dissipation theory, the quantity of primary interest is the reduced density matrix of the impurity system, . Here, is the density matrix of the total system (impurity plus electron reservoir); and denotes the trace over all reservoir’s degrees of freedom. The dynamics of the reduced density matrix can be characterized by a Liouville-space quantum propagator as follows,

(2) |

The specific form of can be expressed using the Feynmann–Vernon path integral representation ofFey63118 ()

(3) |

Here, is the classical action functional associated with . For the convenience of evaluating expectation values, usually an initial factorization ansatz is adopted for the total system at :

(4) |

with being the equilibrium reduced density matrix of reservoir environment. We emphasize that the uncorrelated total system at is chosen only as a reference for the construction of . For the calculation of any physical quantity or process, we first solve for the fully correlated stationary state at a finite time , and then use as the initial condition for subsequent dynamics of the impurity system.

It is that accounts for the influence of the dissipative reservoir environment to the properties of the impurity, and it is thus termed as the influence functional, which has the following form:Jin08234703 ()

(5) | ||||

(6) | ||||

(7) | ||||

(8) |

where and . It should be emphasized that the above Eq. (6) is an exact formula, as long as the reservoir environment satisfies Gaussian statistics, which is true for reservoirs consisting of noninteracting electrons.Jin08234703 ()

Although has a rather complicated form, it is apparent that the impurity-reservoir coupling enters (exclusively) through the reservoir correlation function (or memory kernel) . Therefore, with the exact at hand, we can obtain the exact Liouville propagator and hence the exact via the path integral of Eq. (3).

For a noninteracting reservoir, the correlation function is fully determined by the reservoir hybridization function via the following fluctuation-dissipation theorem:

(9) |

Here, is the time-dependent bias voltage applied to the reservoir, . is the Fermi function for electron () or hole (), , and is the equilibrium chemical potential of reservoir. In fact, are exactly equivalent to the “embedding” self energies: and .Jin08234703 ()

Clearly, can be obtained via the formally exact path integral formalism by combining Eqs. (3)–(9). The only requirement is the reservoir should consist of noninteracting electrons, and hence the reservoir environment satisfies Gaussian statistics. With the resulted , the expectation value of any system operator can then be calculated as .

In practice, the path integral formalism will lead to an integro-differential quantum master equation for , and hence very difficult to solve. To circumvent this problem, an HEOM formalism is instead proposed. In the HEOM formalism, the integro-differential equation is replaced by a hierarchical set of linear differential equations. As a consequence, the numerical calculations become much more convenient.

The central step towards the establishment of HEOM is the decomposition of reservoir memory kernel into exponential functions, i.e.,

(10) |

where is the characteristic memory time of th dissipation mode. A number of memory decomposition schemes have been developed, including the conventional Matsubara decomposition scheme,Jin08234703 () a hybrid spectrum decomposition and frequency dispersion scheme,Zhe09164708 () the partial fraction decomposition scheme,Cro09073102 () and the Padé spectrum decomposition scheme.Oza07035123 (); Hu10101106 (); Hu11244106 ()

Most of the existing memory decomposition schemes are based on a contour integral algorithm for Eq. (9) with the use of a residual theorem,Zhe09164708 () and each exponential term in Eq. (10) corresponds to a pole of complex function or in Eq. (9), via the following expansions:

(11) | ||||

(12) |

Therefore, the total number of memory components is , with and being the number of poles for and in the upper or lower -plane, respectively.

Note that Eq. (11) is a formal expansion of the Fermi function, and the choice of is not unique. At any finite temperature, can always be chosen so that the right-hand side of Eq. (11) and their contributions to converge to a preset precision. To the best of our knowledge, the Padé expansion of requires a smallest , and thus is by far the most efficient scheme for Eq. (11).Hu11244106 () In contrast, Eq. (12) is realized by a least square fit of to a linear combination of Lorentzian functions, with being the fitted parameters. The accuracy of the resulted memory decomposition is determined by the quality of the least square fit. For a continuous the fitting of Eq. (12) is usually satisfactory; see for instance Fig. 2 in Sec. II.4. Therefore, the combined Padé-Lorentzian schemeWan13205126 () achieves an accurate and efficient exponential decomposition of via Eq. (10).

With the exponential decomposition of reservoir memory, the HEOM can be cast into a compact form asJin08234703 ()

(13) |

Here, is the reduced density matrix, and are the auxiliary density matrices, with being the truncation level. The multi-component index . The Grassmann superoperators and are defined via their actions on a fermionic/bosonic operator as and , respectively, with denoting the opposite sign of . The electron correlation interaction is contained in the Liouvillian of impurities, .

Moreover, the subscript index set () in an -level auxiliary density matrix belongs to an ordered set of distinct -indices. The number of distinct -indices, , is equal to not only the number of the first-level auxiliary density matrices in total, but also the maximum hierarchical level (). Therefore, the number of the -level auxiliary density matrices is . Then, the total number of unknowns to solve , which dominates the computational cost of the HEOM approach, is the summation over all the truncation level, , as . In this expression, can be evaluated by , in which is determined by the resolution of the bath memory, and is the number of the system states (orbitals) that couple directly to electron reservoirs. This scheme provides an optimal basis to exploit the different characteristic time scales associated with system-reservoir dissipation processes. The hierarchy terminates automatically at for noninteracting ;Jin08234703 () while for involving electron correlation interactions, the solution of Eq. (II.1) must go through systematic tests to confirm its convergence versus . In practice, recent investigations on a general quantum impurity system have indicated that a relatively low () is usually sufficient to yield quantitatively converged results.Li12266403 ()

In the framework of HEOM, the density matrix of the total system is effectively “folded” into the impurity subsystem – represented by the reduced density matrix along with all the auxiliary density matrices . The influence of reservoir environment is fully accounted for by the reservoir memory kernel . As long as the reservoir environment () is noninteracting (and hence satisfies Gaussian statistics), the HEOM formalism is formally exact, without any approximation.

In practice, the numerical results of HEOM approach achieve quantitative accuracy as long as they converge with respect to the number of memory components (of Eq. (10)) and the truncation level . There is no other controlling parameter. For all the results presented in this work, the numerical convergence has been affirmed.

### ii.2 Linear response theory via the HEOM dynamics

The HEOM (II.1) formally define a quantum Liouville propagator , which associates the reduced and the auxiliary density matrices at time to those at time :

(14) |

Here, a bold symbol denotes a quantity in the HEOM Liouville space defined by Eq. (II.1). For instance, is a vector in the HEOM space and is comprised of all the density matrices involved in Eq. (II.1), i.e., .

Denote as an equilibrium-state solution of Eq. (II.1) at a given temperature . The auxiliary density matrices are nonzero, , which reflect the presence of correlations between the impurity and reservoirs. Consider an impurity system initially at thermal equilibrium, its dynamics starts with . In the absence of any probe field, the system equilibrium propagator satisfies the translational invariance in time, , and hence .

If a probe field is applied to the impurity (the perturbation Hamiltonian is and the corresponding HEOM Liouvillian is ), the response of the reduced and auxiliary density matrices can be obtained formally as

(15) |

Equation (15) is formally analogous to the celebrated Hilbert-space time-dependent perturbation theory in the interaction picture. Such formal analogy highlights a straightforward equivalence mapping between the conventional Hilbert-space and the HEOM Liouville-space formulations for the response properties of quantum impurity systems.

Consider, for example, the equilibrium-state two-time correlation function between two arbitrary dynamical variables and of the impurity system. We have

(16) |

Here, denotes inner product between vectors. The application of Eq. (15) leads to the evaluation of correlation and response functions of impurity systems via the linear response theory in the HEOM Liouville space.

Regarding the linear response theory, we are primarily interested in how system properties (such as the expectation value of a system observable ) respond to an external perturbation in the physical subspace of the system. Suppose the perturbation Hamiltonian assumes the form of

(17) |

where is a Hermitian system operator, and is the time-dependent perturbative field which is turned on from . The induced variation of the expectation value of a system observable due to the presence of is

(18) |

Here, has the same vector form as and . Note that all components of the vector are zero.

Taking the first-order term of Eq. (15), we arrive at

(19) |

We then define the time-independent HEOM-space superoperator as

(20) |

whose action can be determined as

(21) |

Inserting Eq.(21) into (19), we have

(22) |

from which we obtain the response function in the HEOM space as

(23) |

We also obtain the correlation function in the HEOM space as

(24) |

Here, the same propagator as that used in Eq. (14) is applied to the HEOM-space vector . This amounts to taking the as the initial condition at time and then solving Eq. (II.1) for the final state vector at time . Then, the expectation value of the system observable is evaluated with such a final state. By setting while noting , Eq. (24) recovers Eq. (II.2) immediately.

In this work, the dynamical observables of our primary interest are Green’s function and spectral function of the impurity. The spectral function is associated with the system correlation function through the fluctuation-dissipation theorem at thermal equilibrium. Consider the retarded single-electron Green’s functions in terms of the correlation functions,

(25) |

where the Green’s function is defined for two arbitrary operators and . We now focus on cases in which (e.g., and for conventional fermionic Green’s functions). In such circumstances, follows by definition. We can further define the spectral function as

(26) |

in which is obtained by Fourier transform of

To summarize, within the framework of HEOM the system Green’s functions are obtained through the following procedures: (1) Solve HEOM (II.1) for the equilibrium state vector consisting of all the reduced and auxiliary density matrices. (2) Construct the HEOM-space vector . (3) Propagate the HEOM (II.1) from to time by using the vector as the initial condition. (4) Evaluate the expectation value of the system observable to obtain the correlation function . (5) The system Green’s functions are obtained by combining the correlation functions, such as via Eq. (25). Clearly, the evaluation of system Green’s functions is exactly based on the linear response theory, and no approximation is made throughout the above procedures.

### ii.3 HEOM evaluation of dynamical quantities

The system dynamical observable for strongly correlated electronic systems, such as the spectral function , can also be evaluated separately for each frequency point. This is achieved by the half-Fourier transform of the correlation function via,

(27) |

Here, the HEOM propagator is formally cast into the form of . At a fixed frequency , the HEOM-space vector is determined by solving the following linear problem with the quasi-minimal residual algorithm:Fre91315 (); Fre93470 ()

(28) |

The system spectral function is then

(29) |

Usually, the band width of the total system is finite.Ani2010 () Consequently, a finite number of frequency points within the bandwidth are sufficient to accurately characterize the system observables. In this work, the number of frequency points used ranges from 64 to 128. Moreover, calculations on each individual frequency points are parallelized to improve the efficiency of HEOM solver.

The retarded Green’s function of the impurity system can be evaluated as follows.

(30) |

Apparently, one can recast Eq. (29) as

(31) |

### ii.4 HEOM for a DMFT impurity solver

The simplest model to investigate the correlation effect in lattice systems is the one-band Hubbard model

(32) |

where represents creation (annihilation) of an electron on site with spin , and is the hopping parameter between site and . Under the approximation of infinite dimension, the hopping parameter is scaled as with being the lattice dimension or the number of nearest neighbors and a constant. Meanwhile, the electron-electron self-energy becomes local, which is represented by . Then the local Green’s function of the lattice model is

(33) |

where , with being a positive minimal, is the noninteracting lattice Green’s function component of the specified energy . Considering the corresponding density of states , the above local Green’s function can be also calculated by the integration over energy variable

(34) |

The locality of the lattice Green’s function makes it possible to map the lattice problem to effective impurity problem with exactly the same on-site electron-electron self-energy. Kot2001 () The Hamiltonian of the single impurity model is expressed as

(35) |

In the HEOM formalism the effects of coupling reservoir dynamics, as from the last term of Eq. (II.4), are accounted for by the hybridization function . This is just the imaginary part of the reservoir self-energy, , that can also be evaluated from the noninteracting Green’s function . The impurity Green’s function is given by , where denotes electron-electron self-energy. By mapping the lattice model to an effective impurity model, the following relations should be satisfied:Ani2010 ()

(36) |

In this way the self-consistent equations can be derived as

(37) |

The overall framework of HEOM-based DMFT is presented in the flow chart of Fig. 1(a), while the numerical procedures involved in the HEOM impurity solver is sketched in Fig. 1(b). The HEOM impurity solver starts with the pre-determined system Hamiltonian and the input hybridization function . First, the HEOM of Eq. (II.1) is solved to obtain the equilibrium-state reduced density matrix and auxiliary density matrices . Then, the system correlation functions can be obtained by solving the HEOM-space linear problem of Eq. (28) for each frequency point. Finally, the impurity spectral function and Green’s function are evaluated through Eqs. (29) and (30), which are used to construct a new hybridization function via the lattice Green’s function.

For achieving the optimal efficiency of HEOM, a multi-Lorentzian decomposition schemeZhe10114101 (); Xie12044113 () for the hybridization function is adopted. In this way, the hybridization function is spanned by a set of Lorentzian functions . The parameters and for each Lorentzian function are obtained by a least square fit. As the hybridization function is updated in each DMFT iteration, the number of the Lorentzian functions is tuned case by case to find the minimal with sufficient fitting accuracy.

As shown in Fig. 2, two typical hybridization functions are fitted by Lorentzian functions. For a hybridization function with a well-defined peak structure, is sufficient to reach a reasonable accuracy; see Fig. 2(a). However, for a hybridization function possessing a complex peak structure, a larger is necessary to reach the same level of fitting accuracy. As shown in Fig. 2(b), the presence of a small and narrow peak at the Fermi energy requires for an accurate fitting. At present, the major bottleneck of the HEOM approach is the physical memory required to store all the auxiliary density matrices (rather than the CPU time), with the total number of , where , as described earlier (see Sec. II.1) for single impurity systems, with being the number of the Lorentzian function plus the number of Padé decomposition term for the Fermi function. For a typical calculation in this work, the maximum memory required using Lorentzian number is about 2361 MB, while the memory required for is 11366 MB (see Table 1).

The nonperturbative feature of HEOM ensures it is applicable to SCS of a wide range of system parameters. It goes with a well-defined convergence scheme determined by the truncation level , without any tricky parameters to deal with. The convergence of the dynamical quantities for quantum impurity systems is guaranteed once convergence with respect to truncation level is reached. The minimal truncation level required to achieve convergence is closely dependent on the system configuration and the surrounding environment, such as the strength of electron correlation, the system-reservoir coupling, and the temperature. It is difficult to have an a priori estimation for the required minimal . In practice, the convergence of HEOM with respect to is tested case by case. In Fig. 3, we demonstrate such a convergence test for different temperatures. The spectral function is monitored for various combinations of and . For the lowest temperature involved in this work, is necessary (data for are not shown). As the temperature increases to and , and are sufficient for convergence, respectively. In this work, all the results presented are affirmed as already converged at .

In fact, the HEOM solver retains its accuracy and becomes even more efficient at high temperatures. In practice, the HEOM solver is applicable to an arbitrary finite temperature, provided with sufficient computational resources. Therefore, it is possible to investigate temperature sensitive electron correlation effects by continuously varying the temperature. Furthermore, the HEOM approach works in the real frequency domain, and hence avoids invoking the numerically ill-defined analytical continuation problems for approaches working in the imaginary frequency domain.

CPU time (s) | Memory (MB) | ||
---|---|---|---|

1 | 9 | 2 | 0.2 |

2 | 9 | 34 | 9 |

3 | 9 | 959 | 175 |

4 | 9 | 48934 | 2361 |

5 | 9 | 302174 | 21378 |

4 | 14 | 99356 | 5584 |

4 | 18 | 1637256 | 11366 |

The major limitation of the current HEOM implementation is that the computational cost increases dramatically as the system temperature decreases. This is because for a lower temperature, a larger number of memory components ( in Eq. (10)) and a higher truncation level () are necessary to ensure numerical convergence, leading to a rapid growth of the total number of auxiliary density matrices involved in the HEOM formalism. In Table 1, we collect the CPU time and memory usage information for some typical HEOM impurity solver calculations. In practical calculations, parallel programming techniques have been employed to reduce the computational time. The lowest temperature we are able to access with the computational resources at our disposal is , in which is the effective bandwidth of the noninteracting lattice density of states. Actually, the current HEOM implementation is mainly limited by the insufficient physical memory to store all the auxiliary density matrices (rather than CPU time), as mentioned earlier. Fortunately, it is possible to improve substantially the efficiency of HEOM by removing part of the auxiliary density matrices those are exactly zero based on physical consideration. This forthcoming improvement may realize the investigations of SCS using the HEOM based DMFT at further lower temperatures.

## Iii results and discussions

The Mott metal-insulator transition is one of the best test beds for the accuracy and efficiency of the DMFT implementation. In this work we focus on the Mott transition in the half-filling paramagnetic phase at finite temperature. Employing our developed DMFT approach using the HEOM impurity solver, we calculate the spectral functions of SCS, and compare our results with those obtained by Bulla et al.Bul99136 (); Bul01045103 () who used an NRG impurity solver along with a “patching method” for the evaluation of spectral functions.

The calculations are performed on both the hypercubic lattice with infinite dimension and the Bethe lattice with infinite coordination number. The noninteracting density of states for these two lattices are

(38) |

and

(39) |

The effective bandwidth, defined as , is calculated to be for both the hypercubic and the Bethe lattices.Bul99136 () In this work, is set to unity, and is set as the energy scale. To ensure the half-filling of the impurity system, a symmetric Anderson impurity model is employed, and the chemical potential is set to zero.

The convergence of the DMFT method using the HEOM impurity solver versus the truncation level is examined, and shown in Fig. 4. At the lowest temperature performed in our study, , the calculated spectral function is converged at . At a truncation level lower than (), the resulted are not accurate enough, even qualitatively. Sometimes even unphysical results may appear, such as the negative spectral function at . The typical time for a single DMFT iteration at truncation level is about 6 hours on a PC with a 3.0GHz frequency CPU. When the temperature increases to , is sufficient for the convergence of the resulting spectral function. At this truncation level, the typical time for a single DMFT iteration dramatically decreases to less than 30 minutes. In this work, the convergence of all the data represented has been carefully checked.

The DMFT iteration, calculated by using the proposed method, is usually converged within 10 iterations for systems away from the metal-insulator transition region. As an example, we show in Fig. 5 the convergence of the spectral function and self-energy in hypercubic lattice by iteration, with the parameter and . The initial hybridization function can be obtained by setting the initial electron-electron self-energy to zero. Alternatively, in this work we choose to use the initial hybridization in the form of a single Lorentzian function centered at the Fermi level. We have confirmed that they converge to the same final self-energy and spectral function. By adopting the single Lorentzian type initial hybridization function, numerical convergence is reached within 10 DMFT iterations. However, more iterations are necessary to achieve convergence when the system is near to the phase transition area.

We report in Fig. 6 the spectral functions , with different values of at finite temperature , for both the hypercubic lattice and the Bethe lattice. These two lattices have similar behaviors. For both show metallic three-peak characteristics, with a quasi-particle peak at the Fermi level. Further increasing the values of , the height of quasi-particle peak decreases abruptly to a nearly zero value at a critical (see the line in the inset of Fig. 9), leading to a transition from metal to insulator. The transition behavior is very similar to the zero temperature spectral function,Bul99136 () except for that the quasi-particle peaks are slightly broadened and lowered at finite temperature.

The corresponding converged self-energies are shown in Fig. 7 for both the hypercubic and the Bethe lattice systems at the same temperature. For , the two-peak structure of the imaginary part of the self-energy gives rise to the three-peak structure in the spectral function. This is because the peaks in the imaginary part of the self-energy induce dips in the spectral function. With the increase of the value of , the peaks in the imaginary self-energy around the Fermi energy show up, leading to the vanishing of the quasi-particle peaks in the spectral function. The peaks around the Fermi energy in the imaginary self-energy are also broadened by the finite temperature. Thus the critical value of for metal-insulator transition would not depend on temperature.Bul01045103 ()

For the numerical aspect, the computational cost (time and memory) of present HEOM approach grows rapidly with the lowered as higher truncation level is needed for convergence. Currently the lowest temperature that is available for the present HEOM algorithm as DMFT impurity solver is .

As will be shown in the following part of this section, this currently HEOM achievable temperature is already lower than the critical temperature for first-order Mott-Hubbard metal-insulator transition of the Bethe lattice.Bul01045103 () One of the main feature for is the coexistence of both the metallic and insulating solutions within a certain range of .Bul01045103 (); Fel04136405 () This phenomenon is also found in our calculation, as shown in Fig. 8. When at , a metallic three-peak spectral function is obtained if the initial spectral function is also metallic. In this case, the converged spectral function of a lattice system with a smaller is used as the initial guess. Whereas if a converged insulating spectral function corresponding to a lattice system with a larger is feeded as the initial guess, the resulted spectral function has a dramatically depressed component around the Fermi energy. It is noted that for the “insulating” solution in Fig. 8, the spectral function around the Fermi energy is nonzero although the value is very small. This is possibly because that the finite temperature of is very close to .

One of the advantages of the HEOM approach is that it is even more convenient to achieve quantitative accuracy at higher temperatures. This feature ensures the continuous variation of the temperature to investigate possible temperature sensitive phase transitions. It should be noted that when the temperature increases, the efficiency of the HEOM approach is dramatically enhanced. For example, the calculation speed for single DMFT iteration at is more than ten times faster than that at . Figure 9 depicts the spectral functions and both the real and imaginary parts of the self-energies of Bethe lattice. The value of is fixed to , while the temperature is varied from to . At , the spectral function has the typical three-peak feature. When the temperature increases, the quasi-particle peak shrinks continuously. This change is also reflected in the corresponding self-energy. For example, by increasing the temperature, the two peaks in the imaginary part of the self-energy at low temperature become smaller and move toward the Fermi energy. Finally, these peaks disappear at high temperature , and a peak around the Fermi energy shows up.

Although the quasi-particle peak shrinks in the spectral function as shown in Fig. 9, the system of remains in the metallic state. However, by tuning the values of , metal-insulator transition may occur as the temperature varies.

In the following, we focus on the dependence of the spectral function at the Fermi energy at different temperatures, and compare to the results obtained by the NRG method.Bul01045103 () The comparison is shown in Fig. 10. At some low temperatures, for instance, or , the HEOM results agree quantitatively with those obtained by the NRG method. As increases from zero, decreases to zero at a certain value of . Such a value depends on the temperature. Moreover, as shown in the inset, while the HEOM calculated gradually decreases to zero at , at the curve exhibits an abrupt decrease around . This indicates a first-order phase transition at .

As the temperature increases further, the HEOM results start to deviate from the NRG data; see Fig. 10. The deviation is one-sided, meaning that the HEOM calculated is always smaller than that obtained by the NRG method, and the magnitude of deviation increases consistently versus temperature.

It has been elaborated in the Sec. II that the HEOM approach is in principle numerically exact. Meanwhile, the NRG method is also considered as a numerically exact impurity solver. Therefore, it is intriguing why the two in-principle exact approaches would give different results at the somewhat high temperature . To clarify this issue, we emphasize here that the numerical exactness of a certain approach is achieved only when the computation results converge quantitatively with respect to all the involving parameters.

To understand the discrepancy between the HEOM and NRG results obtained at the high temperature, we now check the numerical convergence of both impurity solvers, for the same system and with the same hybridization function.

For the HEOM approach, there are only two controlling parameters, the memory component in Eq. (10) and the truncation level . Its convergence with respect to these two parameters have already been affirmed; see Figs. 2, 3 and 4. It can be seen that the convergence of the HEOM approach is much easier at high temperature than that at low temperature, as by construction it is intrinsically more favorable at higher temperature cases. The convergence of all the HEOM results shown in Fig. 10 and Fig. 11 have also been carefully checked.

In contrast, the NRG approach used in Ref. Bul01045103, involves several controlling parameters. The most important two are the discretization parameter and the number of kept states . The NRG results are considered to be numerically exact if they converge as and . However, it is usually very difficult to reach at high temperature, due to the exponential increase of the computational cost. The NRG dataBul01045103 () we cited in Fig. 10 were obtained under , and the convergence was not reported. Therefore, we carry out a check of its convergence for both low and high temperatures; see Fig. 11.

Figure 11 shows that the NRG results converge rapidly with respect to and at low temperature, but rather slowly at high temperature, which is very different (or opposite to) the HEOM approach. At the low temperature , the overall shape of the best NRG curve with and is very close to the HEOM converged spectral function, which is consistent with the fact that at low temperature the data in Fig. 10 are almost identical. Here means that the actual value of is chosen between 2048 and 2100 with the largest energy gap. At high temperature , the best NRG spectral function with and in our calculation is still not guaranteed to be converged. Moreover, its overall shape dramatically differs from the converged HEOM one. It indicates that large uncertainties reside in the NRG data shown in Fig. 10.

It should be pointed out that the NRG data of Ref. Bul01045103, used for the comparison in Fig. 10 do not represent the best possible NRG results. In the past decade, the quality of the spectral function produced by NRG has been improved dramatically by the new algorithms such as the density matrix NRGHof001508 () and the full density matrix NRG.Wei07076402 (); Pet06245114 () We have also checked the convergence of the corresponding spectral function at using the full density matrix NRG approach; see Supplemental Material.Sup () The resulted spectral function is overall much closer to the HEOM counterpart, including the heights of the Kondo and the Hubbard peaks, and the overall line shape of the dip between them.Sup ()

At this stage, the origin of the remaining difference between the HEOM and NRG spectral functions displayed in Fig. 11(a) is unclear. Clarifying this issue requires a more careful and comprehensive assessment of the NRG methods (since the HEOM data have converged). However, this could be technically very difficult, and is beyond the scope of this paper. Therefore, this problem is left open for further investigations.

## Iv Concluding Remarks

To summarize, the HEOM approach is used as the impurity solver for DMFT method. This method is employed to investigate the Mott-Hubbard metal-insulator transition in both the hypercubic and the Bethe lattices at finite temperature. At low temperatures the results obtained by the HEOM based DMFT method agree well with those using the NRG method as the impurity solver. At higher temperatures, the HEOM approach becomes much more efficient. Thus the HEOM approach provides a complement to the NRG impurity solver, in the sense that it is highly efficient and converges rapidly at finite temperatures.

Currently the developed approach is limited to a temperature no lower than . The reason is that at lower temperature the numeric convergence requires larger truncation level. It is however possible to design more efficient reservoir memory decomposition schemes to dramatically reduce the computational resources requirements. The Lorentzian fit scheme for evaluating the hybridization function may also introduce minor residual error under certain circumstances. For example, the quasi-particle peak at the Fermi level becomes quite narrow when the system is very near to the metal-insulator transition point. To distinguish this type of peak in the process of Lorentzian fit, much finer frequency meshes near the Fermi level is necessary. It will sometimes lead to minor fitting error for the narrow quasi-particle peak.

The HEOM approach adopts a general form of the system Hamiltonian. It is quite convenient to extend the current HEOM+DMFT approach to systems other than the half-filled single-band situation. For example, the system occupation number can be tuned away from half-filling by the adjustment of the chemical potential, to simulate the influence of charge doping to SCS. In this situation, the asymmetric Anderson impurity model is to be employed in the HEOM impurity solver, and an additional search is required for the correct chemical potential that reproduces the true occupation number on a lattice site.

The HEOM approach is also applicable to complicated multi-orbital impurity models. For instance, the HEOM approach has been applied to a three-impurity Anderson model; see Fig. S8 in the Supplemental Material of Ref. Li12266403, . This confirms the potential applicability of the HEOM method to more complex systems. Moreover, it is possible to improve further the efficiency of HEOM approach based on physical considerations and making use of the sparse nature of HEOM. This may extend further the applicability of HEOM+DMFT method to further complex strongly correlated systems.

###### Acknowledgements.

The support from the National Science Foundation of China (Grants No. 21303175, No. 21103157, No. 21033008, No. 21233007, No. 11074302, No. 11374363, and No. 21322305), the Strategic Priority Research Program (B) of the CAS (Grant No. XDB01020000), the National Program on Key Basic Research Project of China (Grant No. 2012CB921704), and the Hong Kong UGC (Grant No. AoE/P-04/08-2) and RGC (Grant No. 605012) is gratefully appreciated.## References

- (1) P. Fulde, Electron Correlations in Molecules and Solids, Springer-Verlag, Berlin Heidelberg, 1991.
- (2) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
- (3) M. Jarrell, Phys. Rev. Lett. 69, 168 (1992).
- (4) A. Georges and G. Kotliar, Phys. Rev. B 45, 6479 (1992).
- (5) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- (6) X. Y. Zhang, M. J. Rozenberg, and G. Kotliar, Phys. Rev. Lett. 70, 1666 (1993).
- (7) F. Gebhard, The Mott Metal-Insulator Transition: Models and Methods. (Springer Tracts in Modern Physics), Springer-Verlag, Berlin Heidelberg, 1997.
- (8) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
- (9) A. Georges and W. Krauth, Phys. Rev. B 48, 7167 (1993).
- (10) M. J. Rozenberg, G. Kotliar, and X. Y. Zhang, Phys. Rev. B 49, 10181 (1994).
- (11) M. J. Rozenberg, R. Chitra, and G. Kotliar, Phys. Rev. Lett. 83, 3498 (1999).
- (12) R. Bulla, Phys. Rev. Lett. 83, 136 (1999).
- (13) R. Bulla, T. A. Costi, and D. Vollhardt, Phys. Rev. B 64, 045103 (2001).
- (14) X. Dai, K. Haule, and G. Kotliar, Phys. Rev. B 72, 045111 (2005).
- (15) Y. Ōno, M. Potthoff, and R. Bulla, Phys. Rev. B 67, 035119 (2003).
- (16) S. Florens and A. Georges, Phys. Rev. B 66, 165111.
- (17) Q. Si and J. L. Smith, Phys. Rev. Lett. 77, 3391 (1996).
- (18) R. Chitra and G. Kotliar, Phys. Rev. Lett. 84, 3678 (2000).
- (19) G. Kotliar and S. Savrasov, in New Theoretical Approaches to Strongly Correlated Systems, edited by A. M. Tsvelik, pages 259–301, Kluwer Academic, Dordrecht, 2001.
- (20) V. Drchal, V. Jani, J. Kudrnovsk, V. S. Oudovenko, X. Dai, K. Haule, and G. Kotliar, J. Phys.: Condens. Matter 17, 61 (2005).
- (21) H. Kusunose, Journal of the Physical Society of Japan 75, 054713 (2006).
- (22) A. Toschi, A. A. Katanin, and K. Held, Phys. Rev. B 75, 045118 (2007).
- (23) V. Anisimov, A. Poteryaev, M. Korotin, A. Anokhin, and G. Kotliar, Journal of Physics: Condensed Matter 9, 7359 (1997).
- (24) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
- (25) H. Kajueter, G. Kotliar, and G. Moeller, Phys. Rev. B 53, 16214 (1996).
- (26) N. E. Bickers, Rev. Mod. Phys. 59, 845 (1987).
- (27) E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
- (28) M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 (1994).
- (29) Q. Si, M. J. Rozenberg, G. Kotliar, and A. E. Ruckenstein, Phys. Rev. Lett. 72, 2761 (1994).
- (30) J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 (1986).
- (31) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
- (32) K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
- (33) R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
- (34) W. Hofstetter, Phys. Rev. Lett. 85, 1508 (2000).
- (35) A. Weichselbaum and J. von Delft, Phys. Rev. Lett. 99, 076402 (2007).
- (36) R. Peters, T. Pruschke, and F. B. Anders, Phys. Rev. B 74, 245114 (2006).
- (37) P. W. Anderson, Phys. Rev. 124, 41 (1961).
- (38) J. Hubbard, Proc. Roy. Soc. A 276, 238 (1963).
- (39) P. A. Lee, T. M. Rice, J. W. Serene, L. J. Sham, and J. W. Wilkins, Comments on Cond. Matter Phys. 12, 99 (1986).
- (40) V. Anisimov and Y. Izyumov, Electronic Structure of Strongly Correlated Materials, Springer-Verlag, Berlin Heidelberg, 2010.
- (41) Z. H. Li, N. H. Tong, X. Zheng, D. Hou, J. H. Wei, J. Hu, and Y. J. Yan, Phys. Rev. Lett. 109, 266403 (2012).
- (42) M. Jarrell and J. Gubernatis, Phys. Rep. 269, 133 (1996).
- (43) J. S. Jin, X. Zheng, and Y. J. Yan, J. Chem. Phys. 128, 234703 (2008).
- (44) R. P. Feynman and F. L. Vernon, Jr., Ann. Phys. 24, 118 (1963).
- (45) H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets, World Scientific, Singapore, 2009, 5th ed.
- (46) R. X. Xu, P. Cui, X. Q. Li, Y. Mo, and Y. J. Yan, J. Chem. Phys. 122, 041103 (2005).
- (47) R. X. Xu and Y. J. Yan, Phys. Rev. E 75, 031107 (2007).
- (48) X. Zheng, J. S. Jin, and Y. J. Yan, J. Chem. Phys. 129, 184112 (2008).
- (49) S. K. Wang, X. Zheng, J. S. Jin, and Y. J. Yan, Phys. Rev. B 88, 035129 (2013).
- (50) X. Zheng, J. S. Jin, and Y. J. Yan, New J. Phys. 10, 093016 (2008).
- (51) X. Zheng, J. S. Jin, S. Welack, M. Luo, and Y. J. Yan, J. Chem. Phys. 130, 164708 (2009).
- (52) X. Zheng, Y. J. Yan, and M. Di Ventra, Phys. Rev. Lett. 111, 086601 (2013).
- (53) X. Zheng, R. X. Xu, J. Xu, J. S. Jin, J. Hu, and Y. J. Yan, Prog. Chem. 24, 1129 (2012).
- (54) A. Croy and U. Saalmann, Phys. Rev. B 80, 073102 (2009).
- (55) T. Ozaki, Phys. Rev. B 75, 035123 (2007).
- (56) J. Hu, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 133, 101106 (2010).
- (57) J. Hu, M. Luo, F. Jiang, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 134, 244106 (2011).
- (58) R. Wang, D. Hou, and X. Zheng, Physical Review B 88, 205126 (2013).
- (59) R. W. Freund and N. M. Nachtigal, SIAM J. Numer. Math. 60, 315 (1991).
- (60) R. W. Freund, SIAM J. Sci. Comput. 14, 470 (1993).
- (61) X. Zheng, G. Chen, Y. Mo, S. K. Koo, H. Tian, C. Y. Yam, and Y. J. Yan, J. Chem. Phys. 133, 114101 (2010).
- (62) H. Xie, F. Jiang, H. Tian, X. Zheng, Y. Kwok, S. Chen, C. Y. Yam, Y. J. Yan, and G. Chen, J. Chem. Phys. 137, 044113 (2012).
- (63) M. Feldbacher, K. Held, and F. F. Assaad, Phys. Rev. Lett. 93, 136405 (2004).
- (64) See Supplemental Material xxx for a comparison between the HEOM and full density matrix NRG impurity solvers.