# Hierarchical Liouville-space approach for accurate and universal characterization of quantum impurity systems

## 1General remarks on the HEOM formalism

The hierarchical equations of motion (HEOM) formalism is established based on the Feynman-Vernon influence functional path-integral theory,[1] and implemented with the Grassmann algebra for fermionic dissipatons.[2] Thus, it is formally exact for a general open system coupled to reservoirs bath which satisfies Grassmann Gaussian statistics. The mathematical construction of HEOM has been discussed comprehensively in Refs. . We shall not repeat the tedious derivations, but just address some key features of the HEOM formalism, particularly for those pertinent to the quantum impurity systems studied in the main text of this work. By elucidating these basic features we will discuss the reasons why the HEOM formalism is not just a formally exact theory, but also leads to an accurate, efficient and universal approach for characterization of quantum impurity systems at finite temperatures.

### 1.1Accuracy and efficiency of the HEOM approach

The efficiency of the HEOM approach is rooted mainly at its nonperturbative nature, as discussed later. As a result, it converges often at a relatively low tier level (. This appealing feature is related to the fact that the HEOM formalism is of finite support, as seen below.

() The HEOM formalism is of finite support by a maximum tier level , although numerically it is usually too high to reach. Let us start with the HEOM, Eq. (1) of main text, in which the basic variables are a set of ADOs, , with denoting the converged or truncated tier level in practice. The reduced system density operator is set to be the zeroth-order ADO; *i.e.*, , with representing the trace over all reservoir environment degrees of freedom. In constructing HEOM, the influence of reservoirs enters via the reservoir correlation functions , which are related to the inverse Laplace-Fourier transform of the self-energy functions via and . To achieve formally close HEOM, are expanded in exponential series. Without loss of generality, consider , with the total number of distinct exponential terms being , where the factor comes from the choice of or , is determined by resolution of reservoir memory, and denotes the number of system states through which the itinerary electrons transferring into () and out of () the central system. In writing the HEOM (1) and the involving ADOs, we adopt the abbreviation , so that . Thus, the damping term in Eq. (1) collects all involving exponents. Denote also , with being the opposite sign of . Involved are also the superoperators, and , defined via their actions on an arbitrary operator of fermionic or bosonic nature, or , respectively, by

In particular, while the reduced system density operator and its associated even-tier ADOs are bosonic, the odd-tier are fermionic.

Moreover, the subscription index set () in a generic -tier ADO belongs to an ordered set of *distinct* -indices. Swapping any two of them leads to a minus sign to the ADO, *i.e.*, . As a result, the number of distinct -indices is equal to not only the number of the first-tier ADOs in total, but also the maximum hierarchical level () as by the full HEOM theory. Apparently, the number of the -tier ADOs is , while the total number of unknowns to solve, up to the truncated tier level , is , as . The overall computational cost increases dramatically with both and .

To minimize the computational expenditure while maintaining the quantitative accuracy, various reservoir spectrum decomposition schemes (in relation to the size) have been developed, including the Matsubara spectrum decomposition,[5] the hybrid Matsubara decomposition-frequency dispersion scheme,[6] the partial fractional decomposition,[7] and the recently proposed Padé spectrum decomposition (PSD) scheme.[8] Among all these schemes, the PSD scheme provides an optimal basis to exploit the different characteristic time scales associated with system-reservoir dissipation processes, and hence has the best performance. It leads to an optimal construction of HEOM, with the minimal size of .

() The HEOM formalism is nonperturbative. An important feature of the formalism is that for noninteracting electronic systems (*e.g.* the single-impurity Anderson model (SIAM) with ) that are characterized completely by single-particle properties, the resultant HEOM terminate automatically at second tier level () without approximation.[5] This feature highlights the nonperturbative nature of the HEOM approach. The underlying hierarchy construction resolves nonperturbatively the combined effects of *e-e* interaction, system-reservoir dissipative couplings, and non-Markovian memory of reservoir. For quantum impurity systems with nonzero *e-e* interactions (*e.g.* the SIAM with ), in principle the hierarchy extends to the -tier level, as discussed earlier. In practice, the calculation results usually converge quantitatively at a low truncation tier (), for weak to medium system-reservoir coupling strength. It is important to verify the numerical convergence with the increasing , and such testing procedure has been carried out for all calculations presented in the main text and also in the Supplemental Materials.

The above features imply the following analogies of the HEOM formalism to conventional many-particle quantum mechanics theories. The full HEOM theory (with ) having a total of unknowns resembles the -particle full-space configuration interaction theory, but now for dissipatons resulted from decomposition of reservoir correlation functions. On the other hand, practically the HEOM formalism resembles a certain coupled-cluster method, as it converges rapidly at a relatively low truncation tier level.

### 1.2A universal quantum mechanical method for open systems

The versatility of the HEOM formalism is rooted at the fact that reduced system dynamics is completely representable in a linear space, referred as the HEOM space hereafter. Such a linear space can be deemed as an extension of conventional Liouville space to open systems. In this subsection, we will establish the HEOM space and elaborate its associated algebra. We will demonstrate that the HEOM space is naturally compatible with the Schrödinger picture, the Heisenberg picture, and the interaction picture of the HEOM formalism. Linear response theory can be established in the HEOM space for open quantum systems. Its application to the evaluation of system correlation functions and response functions will be presented in Section 1.3. In particular, the calculation of various dynamical properties of quantum impurity systems, such as the spectral function, dynamic admittance, and local magnetic susceptibility, will be discussed in Section 1.4.

#### The Schrödinger picture versus the Heisenberg picture

Let us start with the Schrödinger/Heisenberg picture in the conventional Liouville (or Hilbert) space of isolated systems. Consider the expectation value of a system dynamical variable,

The last two quantities are the evaluations in the Schrödinger picture and the Heisenberg picture, respectively. The former is described by the Liouville-von Neumann equation, , for isolated systems with . The Heisenberg picture is usually defined only for time-independent system Hamiltonian cases. The corresponding Liouville-space propagator of the isolated system is of the translational invariance in time and given by . The Heisenberg picture is then or equivalently in isolated systems.

To identify the HEOM space and the corresponding Heisenberg picture, we recast Eq. with its HEOM-space evaluation,

The above identities go by

and

with the (by-definition) initial conditions of and . The HEOM-space inner product is defined by

The on the right-hand-side (rhs) has been given by the second identity of Eq. . The bra has the involving operators in a row vector, while the ket a column vector. Consequently, the -indexes set () in the former is of the reversed order from the () in the latter.

Let us write the Schrödinger picture of the HEOM formalism as

It is just the matrix form of HEOM (1); i.e.,

The HEOM-space generator is a matrix of superoperators defined by Eq. . It dictates the HEOM-space propagator , for both and . The HEOM in the Heisenberg picture assumes therefore

with the row vector defined in Eq. .

To have the explicit Heisenberg HEOM expression, consider the matrix form of Eq. , with the part involving explicitly being highlighted as (denoting )

In writing the second row of the above equation for , we exploit the identity of . The Heisenberg counterpart is therefore of

It reads explicitly as

This is the Heisenberg counterpart of Eq. . The Grassmann-parity associated superoperators and act now from the right or backward side on the specified operators that are of either fermionic or bosonic in nature, dictated by that of and then alternates in every subsequent tier. We have the following equivalence to Eq. :

#### The interaction picture and the Dyson equation

The HEOM space interaction picture can also be readily established straightforwardly from its Hilbert-space counterpart. Let us revisit Eq. , with and thus the HEOM generator . Here, is the impurities system Liouvillian that may be time-dependent, for example, in the presence of external pump fields. On the other hand, denotes the additional contribution from the external “probe” field interrogation on the system. Let be the total HEOM propagator in the presence of the probe field, while be the probe-free counterpart; i.e., and , respectively. The standard interaction picture technique leads then readily to the Dyson’s equation with the HEOM dynamics the form of

where denotes *any time before* the specified probe field interrogation. In terms of the interrogated , or , the above Dyson equation can be recast as (setting )

In writing the second identity, we emphasize the fact that appears in only but every diagonal matrix element of the HEOM space generator, so that , where is the unit matrix. Hence, , as inferred from Eq. .

### 1.3Linear response theory via the HEOM dynamics

Consider now the equilibrium HEOM dynamics where the impurities system is initially at thermal equilibrium, at a given temperature , and there is no time-dependent pump field. Therefore, the probe-free propagator is of the translational invariance in time, , and further that in the integrand of Eq. . Together with the Dyson equation, Eq. for the equilibrium dynamics reads

This resembles the celebrated Hilbert-space time-dependent perturbation theory via interaction picture. Note that the thermal equilibrium is a steady-state solution of the HEOM (Equation 7), in the absence of time-dependent external fields; see its evaluation in Section 2. The resulting nonzero ADOs, , reflect the initial system-reservoirs correlations.

The above observations highlight a straightforward equivalence mapping between the conventional Hilbert space and the HEOM space formulations for the evaluations of various correlation and response functions of impurity systems. Consider, for example, the steady-state two-time correlation function between two arbitrary dynamical variables and of the impurity system of interest. We have [*cf*. Eqs. and (Equation 2)]

Similarly, the response function will be

The proof will be exemplified with the response function of Eq. as follows.

Application of Eq. is to be illustrated with the HEOM evaluation of correlation and response functions of impurity systems via the linear response theory (LRT) which applies to the HEOM linear space defined in Section 1.2.

Regarding the LRT, 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 following form

where is a Hermitian system operator, and is the time-dependent perturbative field which takes effect from []. The induced change in the expectation value of a system observable due to the presence of is

By using the first order perturbation in Eq. , we arrive at

To proceed, we define the time-independent HEOM-space superoperator as

whose action can be determined as

Inserting Eq. into , we finally get

from which we are ready to define the response function in HEOM space as

The correlation function in HEOM space can be similarly defined as

where . By setting while noting , Eqs. and (Equation 22) recover Eqs. and (Equation 13) immediately.

The dynamical observables of our primary interest in the present work are the system correlation function and spectral function, which associate with each other through the fluctuation-dissipation theorem (FDT) at thermal equilibrium. Consider, for example, the retarded and advanced single-particle Green’s functions in terms of the correlation functions,

and

Here, the Green’s functions in our theory is defined for two arbitrary operators and . We now focus on cases where (*e.g.*, and for conventional fermionic Green’s functions). In such circumstances, follows by definition. Using this relation in the Fourier transform of Eqs. (Equation 24) leads to . We can further define the spectral function as

Let , which satisfies , the detailed-balance relation at thermal equilibrium. By combining Eqs. (Equation 25), we arrive at the FDT of

### 1.4Evaluation of dynamical observables of quantum impurity systems

In principle, the dynamical correlation of any two system operators ( and ) can be calculated by using the HEOM-space LRT via Eq. . This allows us to evaluate at quantitative accuracy a variety of dynamical observables of experimental significance, such as the Green’s functions and spectral function of the system, local electric and magnetic susceptibilities, and dynamic admittance (inverse of complex impedance), *etc.* The calculation of system Green’s functions and spectral function will be presented in Section 2. Here, we take the local magnetic susceptibility as an example to illustrate how the response properties are evaluated by using the HEOM approach.

The local magnetic susceptibility describes the magnetization of an impurity in linear-order response to an external magnetic field applied locally at the impurity site. For simplicity, we consider the case that the magnetic field is applied only along the -direction, so is the induced spin polarization for electrons on the impurity.

Here, , , and are the local magnetization, local magnetic field, and local magnetic susceptibility of frequency along -direction, respectively. Note that , with the local magnetization operator being

Here, is the gyromagnetic ratio, the Bohr magnetization, and is the electron spin polarization operator along -direction for the impurity.

Upon application of an *ac* local magnetic field, , the system perturbation Hamiltonian is

Apparently, Eq. assumes the same form as Eq. , with playing the role of . The physical observable of primary interest is . Based on the HEOM-space LRT, the local magnetic susceptibility is expressed as

Equation is formally analogous to Eq. , and hence can be readily calculated by solving the quantum Liouville propagator in HEOM space. The *dynamic* (frequency-dependent) local magnetic susceptibility is then obtained via Fourier transform, *i.e.*, . In particular, the *static* local magnetic susceptibility is just its zero-frequency component, . The above formulas can be easily extended to scenario of noncollinear spin.

The above procedures are applicable to a variety of dynamical observables, such as the local electric susceptibility. However, the evaluation of dynamic admittance is somewhat different. It involves the calculation of response function, , where and are the current operator for –reservoir and electron number operator for –reservoir, respectively. Although and are not “system” operators, can still be evaluated by making use of the quantum Liouville propagator in HEOM space. The detailed formulation and numerical demonstrations will be published elsewhere.

It is important to emphasize that the HEOM approach is capable of addressing the nonequilibrium electronic dynamics under arbitrary time-dependent external fields. This thus provides an alternative way for the evaluation of dynamical observables, *i.e.*, through the real-time evolution of system reduced density matrix and the ADOs. The expectation value of physical observable is then extracted from . For instance, the dynamic local magnetic susceptibility can be obtained via , where is the spin polarization for electrons on the impurity under the applied magnetic field . Another example is the dynamic admittance, which can be evaluated via , with and being the time-dependent current and voltage for – and –reservoir, respectively. Note that is extracted from first-tier ADOs; see Ref. []. In our previous work, we have carried out comprehensive studies on the dynamic admittance of open impurity systems; see Refs. []. Moreover, such an approach allows us to go beyond the linear regime and explore system properties in far-from-equilibrium situations. Furthermore, the static response () to external perturbation can be obtained straightforwardly by solving stationary states corresponding to fixed external fields. For instance, the differential conductance data () displayed in Fig. 3 of main text are extracted from steady-state current under different bias voltages via a finite-field treatment.

To summarize, there are two ways to evaluate the dynamical observables of quantum impurity systems within the framework of HEOM theory. One is to calculate the relevant correlation/response function via the HEOM-space LRT; while the other is to calculate the nonequilibrium static/dynamic electronic response to applied external fields. These two ways are completely equivalent in the linear response regime.

## 2Numerical aspects of the HEOM formalism

The numerical implementation of the HEOM method for calculations of dynamical properties of open electronic system, such as transient transport current through quantum dot systems under time-dependent applied voltages, has been addressed in our previous publications; see Refs. []. Our recent improvement on numerical algorithms and codes concerns mainly the follows: (i) improve the numerical efficiency by making use of various advanced numerical techniques, so that the quantum impurity systems in strong correlation regime can be treated accurately with the presently available computational resources; and (ii) extend the HEOM method to evaluation of dynamical observables at equilibrium state.

To reduce the computational cost, we have employed the Padé spectrum decomposition scheme[8] to construct the HEOM, the quasi-minimal residual (QMR) method[13] to solve the large linear problem for stationary states, and the fourth-order Runge-Kutta[15] or Chebyshev[16] method for the time evolution of quantum Liouville propagator in HEOM space. We have also used sparse matrix algorithms and parallel computing techniques to accelerate the calculations. The detailed descriptions of our numerical code will be published elsewhere. Here, we demonstrate the numerical capability of our HEOM code through the calculations on the same single-impurity Anderson model (SIAM) as studied in Fig. 1 of the main text. Tables Table 1 and Table 2 show the computational cost and numerical convergence for the SIAM at relatively high () and low () temperatures, respectively. In particular, the low temperature brings the SIAM system into the strong correlation (Kondo) regime, and the fourth-tier truncation () predicts the number of occupied electrons with a minor uncertainty of less than electron.

A variety of system dynamical observables, such as the spectral density function associated with th single-electron level of system, , can be evaluated in two ways: either with a time-domain scheme or by calculations in the frequency domain.

The time-domain scheme starts with the evaluation of system correlation functions and through the time evolution of the quantum Liouville propagator at . The spectral function is then obtained straightforwardly by a half Fourier transform as follows,

The averaged system spectral function is . However, in practice the time-domain scheme can be very time-consuming, especially for strongly correlated systems at low temperatures, where exhibits prominent Kondo signatures at around zero frequency. In such cases, both and decay rather slowly at large , due to the significant long-time reservoir memory component. Consequently, to calculate at low frequency accurately, it is necessary to have the correlation functions propagate to a long enough time.

An alternative scheme is based on the half-Fourier transformed correlation function via

At a fixed frequency , the HEOM-space vector is determined by solving the following linear problem with the QMR algorithm:

The system spectral function is then evaluated via

Within the frequency-domain scheme, is calculated at one frequency point by solving Eqs. and (Equation 34) once. Therefore, to obtain the whole spectrum of , a full scan for the whole frequency domain is needed.

The equivalence between the two schemes is exemplified with Figure 1, where the spectral function of a SIAM is calculated with both time- and frequency-domain schemes. Apparently, the results out of both schemes agree perfectly with each other, as they should. In practice, a hybrid time- and frequency-domain scheme can be employed. The overall profile of can be obtained with the time-domain scheme, while the detailed structures of the resonance and Kondo peaks which are difficult to resolve by the time-domain algorithm can be exploited with the frequency-domain scheme. This hybrid scheme provides an accurate and efficient method to calculate the full spectrum of for a strongly correlated quantum impurity system in the Kondo regime.

Similarly, the retarded Green’s function is evaluated based on system correlation functions with the time- or frequency-domain scheme as follows.

For an SIAM system, the single-electron retarded Green’s function is formally expressed as

Here, is the self-energy due to the dissipative couplings between the impurity and reservoirs. Its imaginary part gives the hybridization function , which is calculated from the surface Green’s function of the isolated electron reservoirs, or provided as a known property of the reservoirs (such as the Lorentzian lineshape model adopted in this work). is the self-energy due to the presence of electron-electron interaction (the -term) which reflects the information of electron correlation. With evaluated by the HEOM approach via Eq. , can be extracted from Eq. as the only unknown variable.

## 3Remarks on application of HEOM approach to strongly correlated quantum impurity systems

### 3.1Application to single-impurity Anderson model

#### Spectral function of SIAM

The numerical convergence of HEOM is exemplified via the calculated of an asymmetric SIAM system at different truncation tier ; see Fig. 1 in the main text, where the resulting converges uniformly to the numerically exact solution with the increasing . Figure 2 depicts the derivation of calculated at various from the result obtained at . It is observed that the magnitude of difference indeed decreases as approaches . In particular, there is only rather minor deviation between the calculated at and that at , which is barely recognizable at around zero frequency.

We also investigate the temperature dependence of . Figure 3 plots the calculated of the same asymmetric SIAM system at various temperatures. Apparently, as the temperature increases over an order of magnitude, the two resonance peaks at and , which are due to transition between different occupancy states, remain intact. In contrast, the Kondo peak at almost vanishes at the higher temperature. This clearly highlights the significant role of strong electron correlation in quantum impurity systems.

#### Differential conductance and Kondo temperature

We proceed to examine the zero-bias differential conductance of various two-terminal SIAM systems. This is to verify that the low-temperature physical properties obtained by the HEOM approach indeed follows the correct Kondo scaling relation.

The analytical expression for Kondo temperature of a symmetric SIAM is[17]

Equation is derived based on Bethe Ansatz,[17] and is expected to be valid provided that the relation is satisfied. Here, is the impurity-reservoir coupling strength, and is the bandwidth of electron reservoirs.

Figure 5 and Figure 3(a) of main text depict the temperature dependence of zero-bias differential conductance () of symmetric SIAM systems with different values of . Apparently, the versus curves for different values of display rather scattered distribution; see Figure 5(a). In contrast, with the temperature scaled by of Eq. , all curves merge into one at , showing the expected universal Kondo scaling behavior; see Figure 5(b). Whereas the curves become separated at , indicating the systems are outside the Kondo regime. It is important to point out that the second term in the exponent of Eq. is unnegligible to achieve the universal scaling behavior quantitatively. Moreover, note that as the value of approaches , of Eq. is no longer the true Kondo energy scale. This is clearly manifested by the curve in Figure 5(b), where the versus data deviate significantly from the universal scaling curve. Instead, it is found that by taking as the Kondo energy scale, the universal scaling is retrieved; see the pink squares in Figure 5(b). It is thus inferred that the Kondo energy scale depends sensitively on the relative magnitudes of relevant parameters (, , and ). For lower temperatures, a higher truncation level is necessary to achieve quantitative convergence, which is computationally too demanding for present HEOM approach.

By varying the value of or while fixing the others, the SIAM is tuned away from the particle-hole symmetry. The degree of asymmetry can be characterized by a parameter . It is well known that for an asymmetric SIAM (corresponding to ), its Kondo temperature deviates from the form of Eq. .[18] In particular, it had been derived by Krishna-murthy, Wilkins, and Wilson[19] that the Kondo temperature of an asymmetric SIAM in the local moment regime ( and ) can be evaluated via

Here, is the renormalized on-site energy. We examine a series of SIAMs in this regime, with fixed but different values of (and hence different ). The corresponding curves versus and are depicted in Figure 6(a) and (b), respectively. The scattered curves in Figure 6(a) merge into one in Figure 6(b), indicating that of Eq. indeed provides a universal Kondo energy scale for asymmetric SIAM systems. Note that in Figure 6(b) the meV data (represented by pink triangles) deviate slightly from the universal scaling relation at low temperatures, which is possibly due to the fact that such a value of approaches to the reservoir bandwidth of meV.

#### Spectral tail for Kondo resonance

It has been predicted that for symmetric SIAM systems, the Kondo peak of scaled spectral function has a slowly varying logarithmic tail as follows,[20]

Here, with being the high-energy cut-off. It has been verified that Eq. fits well with the NRG result at .[20]

Figure 3(b) of main text plots the spectral function around the Kondo peak calculated by the HEOM approach, which is compared with the analytic curve of Eq. . Clearly, at relatively large scaled frequencies, (say, at ) the numerical results agree remarkably with the analytic curve under all temperatures studied; while in the low-frequency region, the HEOM results converge consistently to the analytic curve (corresponding to ) as the temperature decreases. At even larger frequencies (), the Hubbard peak starts to rise, which overwhelms the tail of Kondo resonance peak, and hence this part of does not show up in the comparison.

Figure 7 depicts the comparison between the HEOM calculated for symmetric SIAMs with various values of and the analytic expression of Eq. . Apparently, the numerical data agree better with the analytic formula as increases. This is due to the fact that the exponent involved in Eq. is asymptotically exact in the strong coupling limit.[20]

The above comparisons clearly affirm that the HEOM approach correctly reproduces the asymptotic logarithmic scaling of Kondo resonance tail.

### 3.2Application to two-impurity Anderson model

As described in the main text, we have applied the HEOM approach to two-impurity Anderson model (TIAM) system. For clarity, a symmetric parallel-coupled two-impurity system is considered. For the TIAM, we focus on how the inter-impurity coupling affects the system dynamical properties (*e.g.* the spectral function of the impurities) in the presence of strong *e-e* interaction. Fig. 3 of the main text reveals that with the increasing , the system exhibits a continuous crossover transition from the Kondo singlet state of individual impurity to the singlet spin-state between the two impurities. This is affirmed by the impressive agreement between the essential features of calculated and differential conductance ().

Actually, with our universal HEOM approach, the full spectrum of can be obtained for any coupling strength . As an example, Figure 4 depicts the full spectral functions of a symmetric TIAM at and , calculated by the HEOM approach at . The same set of parameters as that to plot Fig. 3 in the main text are adopted (in unit of ): , , , and . The sum rule is satisfied to numerical precision. One can see that the full spectrum resolves more resonance peaks at high frequency range, which provides rich information about the evolution of the energetic structure of the impurities system with the variation in .

In Figure 4, the three-peak-structure (the black line) at is nothing but the duplication of the result for individual impurity, with the itinerant electrons screening the local spin moment at each impurity separately. At (the red line), the anti-ferromagnetic coupling () between the local spin moments splits the Kondo peak at . Meanwhile, the inter-impurity coupling splits the high-frequency peaks at and . Then, the three-peak structure of changes to a six-peak one, as shown in Figure 4. The details will be systematically studied in our future work.

### 3.3Application to multi-impurity Anderson models

The HEOM approach is applicable to a general quantum impurity system, and the numerical procedures are easily extended to more complex systems than the single- and two-impurity Anderson models as discussed in Section 3.1 and Section 3.2. To exemplify the applicability of the HEOM approach, we also perform calculations on a three-impurity Anderson model (3IAM) with the system Hamiltonian as follows,

where and Figure 8 exhibits the calculated system spectral functions of 3IAM for two cases: (i) , and (ii) and