Improved strong-coupling perturbation theory of the symmetric Anderson impurity model

TL;DR

This paper enhances strong-coupling perturbation theory for the symmetric Anderson impurity model by improving self-consistent calculations and resummation techniques, resulting in better predictions of impurity spin responses and antiferromagnetism in the Hubbard model.

Contribution

It introduces a more accurate self-consistent scheme and updated continued fraction coefficients, advancing the theoretical description of the Anderson impurity model.

Findings

More accurate impurity spin response to magnetic bias.

Improved description of antiferromagnetism in the Hubbard model.

Enhanced agreement with expected physical behavior.

Abstract

In a previous work (N. H. Tong, Phys. Rev. B 92, 165126 (2015)), an equation-of-motion based series expansion formalism was used to do the second-order strong-coupling expansion for the single-particle Green function of the Anderson impurity model. In this paper, we improve this theory in two aspects. We first use a more accurate scheme to self-consistently calculate the averages that appear in G1. In the resummation process, we use updated coefficients for the continued fraction, guided by the formally exact continued fraction from the Mori-Zwanzig theory. These changes lead to more accurate impurity spin response to the magnetic bias of the bath. Combined with the dynamical mean-field theory, our theory gives improved description for the antiferromagnetism of Hubbard model at half filling.

Equations

$$G^r\left[A(t)|B(t^{\prime })\right]\equiv -i\theta (t-t^{\prime })⟨[A(t),B(t^{\prime }){]}_{\pm }⟩.$$

$$G(A|B{)}_{\omega +i\eta }\equiv {\int }_{-\infty }^{+\infty }{G}^r\left[A(t)|B(t^{\prime })\right]e^{i(\omega +i\eta )(t-t^{\prime })}d(t-t^{\prime })$$

$$\omega G(A|B{)}_{\omega }=⟨[A,B{]}_{\pm }⟩+G{\left([A,H]|B\right)}_{\omega },$$

$$\omega G(A|B{)}_{\omega }=⟨[A,B{]}_{\pm }⟩-G{\left(A|[B,H]\right)}_{\omega }.$$

$$G(A|B{)}_{\omega }$$

$$\omega G_0(A|B{)}_{\omega }=⟨\{A,B\}⟩+G_0([A,H_0]|B{)}_{\omega }.$$

$$\omega G_i(A|B{)}_{\omega }=G_{i-1}([A,H_1]|B{)}_{\omega }+G_i([A,H_0]|B{)}_{\omega }.$$

$$\omega {\Gamma }_n(A|B{)}_{\omega }=G_n([A,H_1]|B{)}_{\omega }+{\Gamma }_n([A,H]|B{)}_{\omega },$$

$$G(A|B{)}_{\omega }\approx G_0(A|B{)}_{\omega }+G_1(A|B{)}_{\omega }+\cdots +G_m(A|B{)}_{\omega }.$$

$$G_{CF}(A|A^{†}{)}_{\omega }=\frac{a_0}{\omega +b_1-\frac{a_1}{\omega +b_2-...}}.$$

$$⟨BA⟩=-\frac{1}{\pi }{\int }_{-\infty }^{\infty }\text{Im}G(A|B{)}_{\varpi +i\eta }\frac{1}{e^{\beta \varpi }+1}d\varpi .$$

$$H_{AIM}=\sum _{k\sigma }{ϵ}_{k\sigma }{c}_{k\sigma }^{†}{c}_{k\sigma }+U{n}_{\uparrow }{n}_{\downarrow }+{ϵ}_d\sum _{\sigma }{n}_{\sigma }+\sum _{k\sigma }{V}_{k\sigma }(c_{k\sigma }^{†}{d}_{\sigma }+h.c.).$$

$$H_0=\sum _{k\sigma }{ϵ}_{k\sigma }{c}_{k\sigma }^{†}{c}_{k\sigma }+U{n}_{\uparrow }{n}_{\downarrow }+{ϵ}_d\sum _{\sigma }{n}_{\sigma },$$

$$H_1=\sum _{k\sigma }{V}_{k\sigma }(c_{k\sigma }^{†}{d}_{\sigma }+h.c.).$$

$${\Delta }_{\sigma }(\omega )=\frac{\Gamma {\omega }_c^2}{(\omega +\sigma \delta \omega {)}^2+{\omega }_c^2}.$$

$${\Gamma }_{\sigma }(\omega )$$

$${\Lambda }_{\sigma }(\omega )$$

$$\{A_{\alpha \beta }^{\sigma },d_{{\sigma }^{\prime }}^{†}\}=\sum _{\mu \nu }{M}_{\alpha \beta ,\mu \nu }^{\sigma {\sigma }^{\prime }}{A}_{\mu \nu },$$

$$\{A_{\alpha \beta }^{\sigma },d_{{\sigma }^{\prime }}\}=\sum _{\mu \nu }{N}_{\alpha \beta ,\mu \nu }^{\sigma {\sigma }^{\prime }}{A}_{\mu \nu },$$

$$M_{\alpha \beta ,\mu \nu }^{\sigma {\sigma }^{\prime }}={\delta }_{\mu \alpha }{f}_{\nu \beta }^{{\sigma }^{\prime }\ast }+{\delta }_{\nu \beta }{f}_{\alpha \mu }^{{\sigma }^{\prime }\ast },$$

$$N_{\alpha \beta ,\mu \nu }^{\sigma {\sigma }^{\prime }}={\delta }_{\mu \alpha }{f}_{\beta \nu }^{{\sigma }^{\prime }}+{\delta }_{\nu \beta }{f}_{\mu \alpha }^{{\sigma }^{\prime }}.$$

$$\left[A_{\alpha \beta }^{\sigma },H_0\right]=(E_{\beta }-E_{\alpha })A_{\alpha \beta }^{\sigma },$$

$$\left[A_{\alpha \beta }^{\sigma },H_1\right]=\sum _{k{\sigma }^{\prime }}\sum _{\mu \nu }{V}_{k{\sigma }^{\prime }}\left(M_{\alpha \beta ,\mu \nu }^{\sigma {\sigma }^{\prime }}{A}_{\mu \nu }{c}_{k{\sigma }^{\prime }}-N_{\alpha \beta ,\mu \nu }^{\sigma {\sigma }^{\prime }}{c}_{k{\sigma }^{\prime }}^{†}{A}_{\mu \nu }\right).$$

$$G(A_{\alpha \beta }^{\sigma }|A_{\gamma \delta }^{\sigma †}{)}_{\omega }\approx G_0(A_{\alpha \beta }^{\sigma }|A_{\gamma \delta }^{\sigma †}{)}_{\omega }+G_1(A_{\alpha \beta }^{\sigma }|A_{\gamma \delta }^{\sigma †}{)}_{\omega }+G_2(A_{\alpha \beta }^{\sigma }|A_{\gamma \delta }^{\sigma †}{)}_{\omega }$$

$$G(c_{k{\sigma }^{\prime }}|A_{\gamma \mu }^{\sigma †}{)}_{\omega }$$

$$⟨A_{\lambda \gamma }{c}_{p{\sigma }^{\prime \prime }}{c}_{k{\sigma }^{\prime }}⟩\approx ⟨A_{\lambda \gamma }⟩⟨c_{p{\sigma }^{\prime \prime }}{c}_{k{\sigma }^{\prime }}{⟩}_0=0,$$

$$⟨c_{p{\sigma }^{\prime \prime }}^{†}{c}_{k{\sigma }^{\prime }}{A}_{\lambda \gamma }⟩\approx ⟨c_{p{\sigma }^{\prime \prime }}^{†}{c}_{k{\sigma }^{\prime }}{⟩}_0⟨A_{\lambda \gamma }⟩.$$

$$G(c_{k{\sigma }^{\prime }}|A_{\gamma \mu }^{\sigma †}{)}_{\omega }=\frac{V_{k{\sigma }^{\prime }}}{\omega -{ϵ}_{k{\sigma }^{\prime }}}G(d_{{\sigma }^{\prime }}|A_{\gamma \mu }^{\sigma †}{)}_{\omega },$$

$$G_1(A_{\alpha \beta }^{\sigma }|A_{\gamma \delta }^{\sigma †}{)}_{\omega }=\frac{K_{\alpha \beta ,\gamma \delta }^{\sigma }(\omega )}{\omega +E_{\alpha }-E_{\beta }},$$

$$K_{\alpha \beta ,\gamma \delta }^{\sigma }(\omega )$$

$${\Phi }_{\mu \gamma }^{\sigma }(\omega )$$

$$G(d_{\sigma }|d_{\sigma }^{†}{)}_{\omega }$$

$$+$$

Full text

1 Introduction

The Anderson impurity model (AIM) is one of the best studied quantum many-body models in condensed matter physics. It describes local electron orbitals with Coulomb interaction (impurities) hybridizing with the itinerant non-interacting electron states (bath) RefPW:1 . AIM has been used to describe metals doped with dilute magnetic impurities RefAC:2 , mesoscopic quantum dot systems, molecular conductors RefTK:3 ; RefLI:4 ; RefYM:4 ; RefRZ:5 ; RefRH:5 , and absorption of atoms onto the surface of materials RefRB:6 ; RefDC:7 . In the framework of the dynamical mean-field theory (DMFT), AIM is the effective model for describing the temporal fluctuations of Hubbard model in large spatial dimensions RefAG:8 . Many theoretical and numerical methods are used to solve the AIM, including Bethe-ansatz RefPB:6 , perturbation expansion RefKY:9 ; RefSE:10 , Green’s function equation of motion RefCL:11 ; RefHG:12 , conserving slave boson theory RefJK:13 , noncrossing approximation RefNE:13 ; RefKH:14 , quantum Monte Carlo (QMC) RefJE:15 ; RefMF:16 ; RefEG:17 , renormalization group RefHR:18 ; RefTA:19 ; RefRB:20 ; RefRB:21 , hierarchical equations of motion (EOM) RefZH:22 , $etc$. Through decades of efforts, the physics of single impurity AIM has been well understood, though the treatment of multi-orbital AIM is still a challenge RefEG:17 ; RefRB:21 . The single-impurity Anderson model is therefore a priority model for testing the validity of new theoretical methods.

In a previous work RefNH:22 (hereafter referred to as I), an EOM based series expansion method for Green’s functions (GF) was developed and applied to AIM with a spin bias in the bath. The local GF $G(d_{\sigma}|d_{\sigma}^{\dagger})_{\omega}$ was expanded to the second order of hybridization $V_{k\sigma}$. Resummation was then carried out in the continued fraction (CF) formalism to approximately sum up the series to infinite order, recovering the correct analytical structure of GF. The obtained static averages and local density of states are in quantitative agreement with those from numerical renormalization group (NRG) calculations at intermediate to high temperatures and for small spin bias of the bath.

However, it is notable that in this theory the impurity spin polarization by the magnetic bias of the bath is significantly underestimated, especially at large $U$, large bias, and low temperature regime. This undermines the value of this theory in the study of quantum dots with ferromagnetic leads and of the magnetic phase of lattice models via DMFT. This deficiency is due to the over simplified approximation, i.e., the atomic approximation, used in I to calculate the higher order static correlation functions. We also found that the CF used for resummation, though recovers the exact form in the atomic limit, does not apply at finite hybridization and arbitrary polarization. Proper treatment of correlations beyond the $V^{2}_{k\sigma}$ order is therefore vital for quantitative accuracy, even in the strong-coupling regime where a second-order expansion in $V_{k\sigma}$ is supposed to be adequate.

In this paper, we aim to improve the second-order strong-coupling series expansion in the two aspects stated above, by making the following changes to the original theory: (i) We do the exact self-consistent calculation for the averages appearing in the first order GF $G_{1}(d_{\sigma}|d_{\sigma}^{{\dagger}})_{\omega}$ (instead of using the atomic approximation for them); (ii) We use improved CF resummation formalism, containing the exact coefficient $a_{1}$ from the Mori-Zwanzig theory RefMZ:19 and revised the expression for the coefficient $b_{2}$; (iii) All averages are calculated self-consistently from $G(d_{\sigma}|d_{\sigma}^{\dagger})_{\omega}$ alone (instead of being calculated from separate resumed GFs). With these changes in the theory, we achieve quantitative improvement in the impurity spin polarization under the magnetic bias of bath.

The outline of the paper is as follows. In Sec.2, a brief summary is given for the self-consistent EOM series expansion and the CF resummation. In Sec.3, the improved strong-coupling expansion is carried out for the single-impurity AIM to $V^{2}_{k\sigma}$ order. In Sec.4, we compare the numerical results from the improved theory, the original theory, and NRG. We also use the improved theory as an impurity solver for DMFT and study the local density of states of the half-filled Hubbard model in the antiferromagnetic phase. The discussion and summary are given in Sec.5.

2 General formalism

The general formalism of the EOM-based self-consistent series expansion for GFs was developed in I. For the sake of completeness, here we give a brief overview. We consider the following retarded GF of operators $A(t)$ and $B(t^{\prime})$,

[TABLE]

Here, $A(t)=e^{iHt}Ae^{-iHt}$. $[A,B]_{\pm}\equiv AB\pm BA$. $\langle O\rangle$ represents the thermodynamic average of the operator $O$ in the equilibrium state of $H$. Here and below, we use the natural unit $\hbar=1$. The Fourier transform of $G^{r}[A(t)|B(t^{\prime})]$ is denoted as

[TABLE]

where $\eta$ is an infinitesimal positive number. On the real frequency axis, the EOM of this double-time GF reads

[TABLE]

Below, we only use the Fermion-type GF with the anti-commutator $[A,B]_{+}=\{A,B\}$ in the EOM.

2.1 Self-consistent series expansion

We decompose the given Hamiltonian $H$ as $H=H_{0}+\lambda H_{1}$. $H_{0}$ is the exactly solvable part and $H_{1}$ is regarded as the perturbation. $\lambda$ is the formal expansion parameter and will be set as unity after the calculation. To develop an EOM-based self-consistent series expansion for GF, we formally write

[TABLE]

Here, $G_{i}(A|B)_{\omega}\propto\lambda^{i}$ ($i=0,1,...,n$) is the $i$-th order GF. $\Gamma_{n}(A|B)_{\omega}$ is the residue of this expansion at order $n$ .

From the left-side EOM Eq.(3), the zeroth-order GF is defined by the EOM

[TABLE]

Here, $\langle\{A,B\}\rangle$ is the full average to be calculated self-consistently after the GF is formally obtained. By putting Eq.(5) into Eq.(3) and comparing the two sides of equation order by order, we obtain the EOM for the $i$-th order GF $G_{i}(A|B)_{\omega}$ ( $i{\geqslant}1$) as

[TABLE]

Each order of GF can be solved exactly if the commutator series $[A,H_{0}]$, $[[A,H_{0}],H_{0}]$, $...$ closes automatically. The $n$th-order residue $\Gamma_{n}(A|B)_{\omega}$ satisfies

[TABLE]

which cannot be solved exactly in general. The series expansion could be deduced also from the right-side EOM Eq.(4), which gives different $G_{i}(A|B)_{\omega}$’s. Hereafter, we will use the series from the left-side EOM, Eqs.(6)-(8). This theory is named self-consistent series expansion.

2.2 CF resummation

A physically acceptable GF should be Lehmann representable, $i.e.$, it consists of real simple poles. Suppose we have obtained the self-consistent series expansion of $G(A|B)_{\omega}$ up to the $m$-th order,

[TABLE]

Using the CF resummation RefSP:23 , we can obtain the GF $G_{CF}(A|A^{\dagger})_{\omega}$ from Eq.(9) (taking $B=A^{\dagger}$)

[TABLE]

It was proven that for $a_{l}$, $b_{l+1}$ real and $a_{l}\geqslant 0$ ($l=0,1,...$), Eq.(10) is Lehmann representable RefJG:24 . In order to determine $a_{l}$ and $b_{l}$, we can expand $G(A|A^{\dagger})_{\omega}$ and $G_{CF}(A|A^{\dagger})_{\omega}$ into Taylor series of $1/\omega$ and require that for every $n\in[1,+\infty]$, the $\omega^{-n}$ terms of $G(A|B)_{\omega}$ and $G_{CF}(A|B)_{\omega}$ equal on the level of $\lambda^{m}$. Here, $m$ is the order of perturbation. In this paper, we take $\lambda=V_{k}$ and $m=2$. To calculate the averages (all real in this work) appearing in the expansion, we use the fluctuation-dissipation theorem,

[TABLE]

3 Strong-coupling expansion for Anderson impurity model

We consider the single-impurity Anderson model

[TABLE]

Here, $c_{k\sigma}$/$d_{\sigma}$ is the bath/impurity electron annihilation operator. We split $H_{AIM}=H_{0}+H_{1}$ as

[TABLE]

and

[TABLE]

$H_{0}$ is exactly solvable due to the decoupling of impurity from bath. $H_{1}$ is treated as a perturbation. For the hybridization function $\Delta_{\sigma}(\omega)\equiv\sum_{k}V_{k\sigma}^{2}\delta(\omega-\epsilon_{k\sigma})$, we use a Lorentzian form

[TABLE]

Here, $\Gamma$ is the strength of the hybridization. $\delta\omega$ is the spin bias on the bath electrons, which is used to describe the magnetic electrode in quantum dot systems or the magnetic phase of a lattice Hamiltonian in DMFT. $\sigma=1$ for spin up and $-1$ for spin down. $\omega_{c}$ is set as the energy unit. As in paper I, here we only focus on the particle-hole symmetric case. Generalization to the particle-hole asymmetric case is straightforward. For the above hybridization function that fulfils $\Delta_{\sigma}(-\omega)=\Delta_{\bar{\sigma}}(\omega)$, the particle-hole symmetric point is located at $\epsilon_{d}=-U/2$. We also define two related intermediate functions

[TABLE]

They fulfil the relations $\Gamma_{\sigma}(-\omega)=-\Gamma_{\bar{\sigma}}(\omega)$ and $\Lambda_{\sigma}(-\omega)=\Lambda_{\bar{\sigma}}(\omega)-\Gamma_{\bar{\sigma}}(\omega)$.

3.1 Second order $V_{k\sigma}$ expansion of GF

In I, a self-consistent expansion of the impurity GF was developed up to $V_{k\sigma}^{2}$ order in the standard basis operator (SBO) formalism. In this paper, these expressions are kept intact and we only modify the schemes for the subsequent resummation and self-consistent calculation. For completeness, in this part we summarize the notations and results of I.

The standard basis operators (SBOs) RefSB:25 ; RefHL:26 used in I are the excitation operators of the local Hamiltonian $h_{0}=Un_{\uparrow}n_{\downarrow}+\epsilon_{d}\sum_{\sigma}n_{\sigma}$. Denoting $|\mu\rangle$ ($\mu=1\sim 4$) as the eigenstates of $h_{0}$, i.e., $h_{0}|\mu\rangle=E_{\mu}|\mu\rangle$, we have $|1\rangle=d_{\uparrow}^{{\dagger}}|0\rangle$, $|2\rangle=d_{\downarrow}^{{\dagger}}|0\rangle$, $|3\rangle=|0\rangle$, and $|4\rangle=d_{\uparrow}^{{\dagger}}d_{\downarrow}^{{\dagger}}|0\rangle$. The corresponding eigen energies are $E_{1}=E_{2}=\epsilon_{d}$, $E_{3}=0$, and $E_{4}=U+2\epsilon_{d}$. The SBOs $\{A_{\alpha\beta}\}$ are defined as $A_{\alpha\beta}\equiv|\alpha\rangle\langle\beta|$. Obviously, we have $\sum_{\alpha}A_{\alpha\alpha}=1$ and $A_{\alpha\beta}A_{\mu\nu}=\delta_{\beta\mu}A_{\alpha\nu}$. The impurity annihilation operator is expressed as $d_{\sigma}=\sum_{\mu\nu}f_{\mu\nu}^{\sigma}A_{\mu\nu}^{\sigma}$. $f^{\sigma}_{\mu\nu}=0$ except for $f_{31}^{\uparrow}=f_{24}^{\uparrow}=1$ and $f_{32}^{\downarrow}=-f_{14}^{\downarrow}=1$. Below we use $A^{\sigma}_{\mu\nu}$ to denote the SBO which net annihilates or creates an electron with spin $\sigma$. We use $A_{\mu\nu}$ without superscript to denote the SBO that changes the number of electrons by an even number. Therefore, $A^{\sigma}_{\mu\nu}$ is Grassman odd and $A_{\mu\nu}$ is Grassman even. We have $\{A_{\mu\nu}^{\sigma},c_{k\sigma}\}=0$ and $\left[A_{\mu\nu},c_{k\sigma}\right]=0$.

With these definitions, we have

[TABLE]

and

[TABLE]

The commutators between SBOs and $H_{0}$ and $H_{1}$ are respectively

[TABLE]

and

[TABLE]

The impurity GF $G(d_{\sigma}|d_{\sigma}^{{\dagger}})_{\omega}$ can be decomposed into $G(d_{\sigma}|d_{\sigma}^{{\dagger}})_{\omega}=\sum_{\alpha\beta\gamma\delta}f_{\alpha\beta}^{\sigma}f_{\gamma\delta}^{\sigma\ast}G(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$. The self-consistent strong-coupling expansion

[TABLE]

was obtained in I. We will not repeat the expressions but refer the readers to I for details: Eq.(84) for $G_{0}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$, Eq.(88) for $G_{1}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$, and Eqs.(89)-(91) for $G_{2}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$.

3.2 Improved self-consistent calculation

These GF components in Eq.(21) were expressed in terms of the unknown averages such as $\langle A_{\delta\beta}\rangle$ (in $G_{0}$), $\langle A_{\gamma\mu}^{\sigma{\dagger}}c_{k\sigma^{\prime}}\rangle$ (in $G_{1}$), and $\langle A_{\lambda\gamma}c_{p\sigma^{\prime\prime}}c_{k\sigma^{\prime}}\rangle$, $\langle A_{\lambda\gamma}c_{p\sigma^{\prime\prime}}^{{\dagger}}c_{k\sigma^{\prime}}\rangle$ (in $G_{2}$). In I, they were calculated from the following minimum approximations that ensures Eq.(3.1) being exact at the order $V_{k\sigma}^{2}$.

$\langle A_{\delta\beta}\rangle$ in $G_{0}$ was calculated from the resummed GF $G_{CF}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$. The averages of the type $\langle Ac\rangle$ in $G_{1}$ were calculated from the corresponding GFs at $V_{k\sigma}$ order. For instance, $\langle A_{\gamma\mu}^{\sigma{\dagger}}c_{k\sigma^{\prime}}\rangle$ was calculated via the fluctuation-dissipation theorem from the approximate GF

[TABLE]

Other averages in $G_{1}$ can be obtained by symmetry transformations from it.

The averages of the type $\langle Acc\rangle$ or $\langle Ac^{{\dagger}}c\rangle$ in $G_{2}$ describe the spin exchange and pair hopping between impurity and bath, being important for the description of Kondo effect. In I, they were calculated by a truncation valid at the order $V_{k}^{0}$,

[TABLE]

$\langle O\rangle_{0}$ means the zero-th order average. Eqs.(22) and (3.2) guarantees that the truncated series of $G(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$ is exact at $V_{k\sigma}^{2}$ order.

As analyzed in the introduction, merely keeping the exact $V_{k\sigma}^{2}$ order is not sufficient for certain purpose even in the strong-coupling limit, for which the accuracy of higher order terms are vital. In this section, we go one step forward than the previous minimum approximation. Instead of using the approximation Eq.(22), we calculate $\langle A_{\gamma\mu}^{\sigma{\dagger}}c_{k\sigma^{\prime}}\rangle$ in $G_{1}$ from the following exact relation

[TABLE]

with $G(d_{\sigma^{\prime}}|A_{\gamma\mu}^{\sigma{\dagger}})_{\omega}$ to be calculated self-consistently from the CF-resummed impurity GFs. Compared with Eq.(22), the present scheme treats higher order hybridization effect more accurately.

With the new self-consistent scheme Eq.(24), we obtain an updated expression for $G_{1}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$ as

[TABLE]

with

[TABLE]

The function $\Phi_{\mu\gamma}^{\sigma}(\omega)$ is given by

[TABLE]

Here, for a given function $g(\varpi)$, we define $\langle g(\varpi)\rangle\equiv-(1/\pi)\int_{-\infty}^{\infty}{\text{I}m}g(\varpi+i\eta)/(e^{\beta\varpi}+1)d\varpi$. If we replace $G(d_{\sigma}|A_{\gamma\mu}^{\sigma{\dagger}})_{\varpi}$ in Eq.(27) with the zero-th order quantity $G_{0}(d_{\sigma}|A_{\gamma\mu}^{\sigma{\dagger}})_{\varpi}$, Eq.(25) reduces to the old result, Eq.(97) of I. Therefore, the new scheme Eq.(24) modifies only the $V_{k\sigma}^{n\geqslant 3}$ order contributions in $G_{1}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$ but keeps the $V_{k\sigma}^{n\leqslant 2}$ order intact. Note that Eq.(25) has an explicit pole at $E_{\beta}-E_{\alpha}$, while the corresponding old one, Eq.(97) of I, has two explicit poles. This suggests that although $K_{\alpha\beta,\gamma\delta}^{\sigma}(\omega)$ in Eq.(26) contains real simple poles apparently, an additional pole emerges in the expansion to $V_{k\sigma}^{2}$ order. This could make a second-order pole in $K_{\alpha\beta,\gamma\delta}^{\sigma}(\omega)$ in the strong coupling limit. Due to this observation, when we do CF resummation, $K_{\alpha\beta,\gamma\delta}^{\sigma}(\omega)$ should not be regarded as a constant (see below).

For the averages in $G_{2}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$, in this paper we still use the old scheme Eq.(3.2). The resulting expression for $G_{2}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$ is not changed and the expression was given by Eq.(98) of I. We will discuss how to improve Eq.(3.2) in Sec.5.

3.3 Improved CF resummation

Besides the self-consistent calculation, there is room for improvement also in the CF resummation process. At the particle-hole symmetric point, the second-order truncated series of the single-particle local GF is simplified into (from Eq.(3.1))

[TABLE]

The expressions for $W_{1}^{\sigma}\sim W_{5}^{\sigma}(\omega)$ are summarized in Appendix A. Due to the improved self-consistent calculation scheme Eq.(24), $W_{3\sim 5}^{\sigma}(\omega)$ are different from those in I. The following exact relations still hold, $W_{1}^{\sigma}+W_{2}^{\sigma}=1$, $W_{1}^{\sigma}-W_{2}^{\sigma}=1-2\langle n_{\bar{\sigma}}\rangle$, and $W_{3}^{\sigma}(\omega)+W_{4}^{\sigma}(\omega)+W_{5}^{\sigma}(\omega)=\Gamma_{\sigma}(\omega)$. They give out certain exact properties of GF such as sum rule and the exactness in the non-interacting limit.

We used the following two level extended CF to do the resummation for Eq.(29),

[TABLE]

Here, $a_{0}$ and $a_{1}$ are real and positive $\omega$-independent constants. $G_{CF}(d_{\sigma}|d_{\sigma}^{{\dagger}})_{\omega}$ is Lehmann representable if $b_{1}(\omega)$ and $b_{2}(\omega)$ contain only real simple poles. Eq.(30) is an extension of Eq.(10) to $\omega$-dependent $b_{i}$’s, which can produce a continuous spectral function from a finite level CF. In I, we expanded Eqs.(29) and (30) into Taylor series of $1/\omega$ and equated the corresponding coefficients to produce

[TABLE]

In the expansion, $W_{3\sim 5}^{\sigma}(\omega)$ and $b_{1\sim 2}(\omega)$ are regarded as constants and the exact relations among $W_{i}^{\sigma}$’s were used. Since Eq.(29) contains second-order pole at most and hence a two level extended CF is sufficient. In I, the coefficients $a_{i}$’s and $b_{i}(\omega)$’s were solved from the above equations at the $V_{k\sigma}^{2}$ order. In particular, in solving the third equation of Eq.(31) for $a_{1}$, $b_{1}(\omega)^{2}$ and $W_{3}^{\sigma}(\omega)-W_{4}^{\sigma}(\omega)$ were discarded since they are at the order $V_{k\sigma}^{4}$ (Note $\langle n_{\bar{\sigma}}\rangle=1/2+\mathcal{O}(V_{k\sigma}^{2})$). This produced

[TABLE]

The $b_{1}(\omega)$ and $b_{2}(\omega)$ obtained above contain real simple poles only and satisfy the Lehmann representability requirement.

The above resummation scheme used in I has some degrees of freedom to modify. First, the forms of $W_{1}^{\sigma}\sim W_{5}^{\sigma}(\omega)$ cannot be uniquely determined by the expression of Eq.(29). Second, the terms of order $V_{k\sigma}^{n\geqslant 3}$ could be kept in the final $a_{i}$ ($i=0,1$) or $b_{i}(\omega)$ ($i=1,2$) to improve the accuracy. In this work, we first find a suitable expression for $W_{i}^{\sigma}$’s which are now different from those of I. Considering that in Eq.(25), $K_{\alpha\beta,\gamma\delta}^{\sigma}(\omega)$ contains the second-order pole at $\omega=E_{\delta}-E_{\gamma}$ in the strong coupling limit, we multiply $\left(\omega+E_{\gamma}-E_{\delta}\right)$ simultaneously on the nominator and the denominator of Eq.(25). This leads to the modified expression for $W_{i}^{\sigma}$’s in Appendix A, which in any cases only have real simple poles. In mapping GF to CF, these $W_{i}^{\sigma}$’s are regarded as constants.

Second, in solving the third equation of Eq.(3.3) for $a_{1}$, instead of neglecting $b_{1}(\omega)^{2}$ completely, we now keep the static part of it $U/2(W_{1}^{\sigma}-W_{2}^{\sigma})=U/2(1-2\langle n_{\bar{\sigma}}\rangle)$. This is motivated by the observation that though each coefficient in Eq.(3.3) is exact at the $V_{k\sigma}^{2}$ order, the response of $\langle n_{\sigma}\rangle$ to $\delta\omega$ is unsatisfactory, especially at large $U$ RefNH:22 . The exact local self-energy in the atomic limit $\Sigma_{\sigma}(\omega)=\epsilon_{d}+U\langle n_{\bar{\sigma}}\rangle+U^{2}\langle n_{\bar{\sigma}}\rangle(1-\langle n_{\bar{\sigma}}\rangle)/\left[\omega-\epsilon_{d}-U(1-\langle n_{\bar{\sigma}}\rangle)\right]$ reminds us that to accurately describe the polarization of impurity spin by the bias of bath, one needs to use the full term $a_{1}=U^{2}\langle n_{\bar{\sigma}}\rangle(1-\langle n_{\bar{\sigma}}\rangle)$ instead of the $V_{k\sigma}^{2}$ order one $a_{1}=(U/2)^{2}$. Preserving the full term $U/2(W_{1}^{\sigma}-W_{2}^{\sigma})$ in the static part of $b_{1}(\omega)$ recovers the atomic self-energy for arbitrary impurity polarization.

Using the above modifications, we obtain the improved CF coefficients as

[TABLE]

Eq.(3.3) satisfies the Lehmann representability and recovers Eq.(3.3) in the unbiased limit ($\delta\omega=0$). In the biased case, Eq.(3.3) is expected to describe the impurity polarization more accurately.

The improved CF resummation scheme can be justified by the exact extended CF formalism of the local GF of AIM. Using the Mori-Zwanzig operator projection method RefMZ:19 , we obtain the formally exact expression

[TABLE]

The frequencies in the above equation are given by $\Omega_{0}=\epsilon_{d}+U\langle n_{\bar{\sigma}}\rangle$ and $\Omega_{1}=\epsilon_{d}+U\left(1-\langle n_{\bar{\sigma}}\rangle\right)+\tilde{\beta}_{\sigma}$, with

[TABLE]

The memory function $iK_{2c}(-i\omega)$ is a GF of higher order operators, whose EOM cannot be solved exactly. The derivation of Eq.(34) is summarized in Appendix B. Comparing Eq.(33) to Eq.(34), we see that our improved CF resummation gives the exact coefficients $a_{0}$, $a_{1}$, and $b_{1}(\omega)$. For $b_{2}(\omega)$, the constant term is also exact, given $W_{5}^{\sigma}(\omega=\infty)=2\sum_{k}V_{k\bar{\sigma}}\langle c_{k\bar{\sigma}}^{{\dagger}}d_{\bar{\sigma}}\left(1-2n_{\sigma}\right)\rangle$. The $\omega$-dependence of $b_{2}(\omega)$ is approximate. Eq.(34) is also recovered at various levels by the recently developed projective truncation approximations with different operator bases Fan1 ; Fan2 .

The improved second-order strong-coupling expansion formulas involve the unknown averages $\langle A_{\alpha\alpha}\rangle$ and the functions $\Phi_{\alpha\beta}^{\sigma}(\omega)$. To calculate them self-consistently, we first write $\langle A_{11}\rangle=\langle n_{\uparrow}\rangle-\langle n_{\uparrow}n_{\downarrow}\rangle$, $\langle A_{22}\rangle=\langle n_{\downarrow}\rangle-\langle n_{\uparrow}n_{\downarrow}\rangle$, $\langle A_{33}\rangle=\langle A_{44}\rangle=(1-\langle A_{11}\rangle-\langle A_{22})/2$. To calculate $\langle n_{\uparrow}n_{\downarrow}\rangle$ and $\Phi_{\alpha\beta}^{\sigma}(\omega)$ , we use

[TABLE]

Consequently, all averages are calculated from the CF-resummed single-particle GF $G_{CF}(d_{\sigma}|d_{\sigma}^{{\dagger}})_{\omega}$. This is different from I where the CF-resummed GFs $G_{CF}(A_{31}^{\uparrow}|A_{31}^{\uparrow{\dagger}})_{\omega}$, $G_{CF}(A_{24}^{\uparrow}|A_{24}^{\uparrow{\dagger}})_{\omega}$, $etc$ are required.

4 Numerical results

In this section, we present numerical results for the improved strong-coupling expansion (denoted as ISC below) and compare them with those from the original strong-coupling expansion (denoted as SC below) and NRG. We use the full density matrix NRG algorithm RefAW:27 supplemented with the self-energy trick for the local density of states (LDOS) RefRB:20 to provide a reference for comparison. We use logarithmic discretization parameter $\Lambda=2.0$ and keep $M_{s}=[356,380]$ states in the NRG calculation RefZH .

Fig.1 presents impurity electron occupancies $\langle n_{\uparrow}\rangle$ and the double occupancy $\langle n_{\uparrow}n_{\downarrow}\rangle$ as functions of $U$ for a spin-polarized bath at $T=0.1$. In Fig.1(a), it is seen that $\langle n_{\uparrow}\rangle$ calculated by the two methods are in good agreement with NRG in the small $U$ regime. For SC, the quantitative agreement is maintained only upto $U=0.5$. The quantity is less accurate for ISC in the small $U$ regime but the better qualitative agreement sustains to larger interaction $U\sim 2.0$. In larger $U$ regime, both results deviate from NRG significantly, with ISC overestimating and SC underestimating the polarization, respectively. In Fig.1(b), the double occupancy obtained by the two methods are all in good agreement with NRG, except for slight deviations at medium $U$ value.

Fig.2 presents $\langle n_{\uparrow}\rangle$ and $\langle n_{\uparrow}n_{\downarrow}\rangle$ as functions of temperature for a negative bias $\delta\omega=-0.2$ and intermediate interaction $U=2.0$. In Fig.2(a), the value of $\langle n_{\uparrow}\rangle$ from ISC is almost same as that of NRG except for $T<0.1$, showing dramatical improvement over SC. In Fig.2(b), $\langle n_{\uparrow}n_{\downarrow}\rangle$ from ISC has no improvement but decreases slightly from the SC result at low temperatures. Neither ISC nor SC produces the upturn of $\langle n_{\uparrow}n_{\downarrow}\rangle$ in the NRG data for $T<0.2$ which is associated with the establishment of Kondo screening and Fermi liquid state. Such deviation is expected because neither SC nor ISC can accurately describe the Kondo effect at low temperatures due to the crude truncation approximation Eq.(3.2) for the averages of the type $\langle Acc\rangle$ and $\langle Ac^{{\dagger}}c\rangle$ in $G_{2}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma\dagger})_{\omega}$.

Fig.3 shows $\langle n_{\uparrow}\rangle$ and $\langle n_{\uparrow}n_{\downarrow}\rangle$ as functions of $\delta\omega$ for an intermediate $U=2.0$. In Fig.3(a), ISC produces a curve of $\langle n_{\uparrow}\rangle(\delta\omega)$ almost identical to NRG upto the large bias regime, significantly surpassing SC. Similar to Fig.2(b), the accuracy of the double occupancy data from ISC deteriorates slightly with respect to SC, both being much smaller than that of NRG. This shows that ISC, with improved evaluation of the averages of type $\langle Ac\rangle$ and with the exact CF coefficients $a_{1}$, only improves the response of the impurity spin to the bias of bath but gains little in the hybridization-induced correlation effect, i.e., the Kondo physics. As stated above, such effect is encoded in the averages of the type $\langle Acc\rangle$ and $\langle Ac^{{\dagger}}c\rangle$ in $G_{2}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma\dagger})_{\omega}$.

To further demonstrate the improvement of ISC over SC, we use them as impurity solvers to study the antiferromagnetism of the half-filled Hubbard model within DMFT. The description of antiferromagnetism in DMFT relies crucially on the accurate description of the impurity density of states under the bias of bath. As is well known, many strong-coupling based approximations such as the Hubbard-I approximation and the alloy analogy approximation cannot give out the antiferromagnetic phase of the Hubbard model at half filling RefGebhard . ISC is expected to do better than SC in this respect due to the improved description of the impurity spin response to the bias of bath. The Hubbard model Hamiltonian reads

[TABLE]

DMFT is exact for this model on the Bethe lattice with coordination number $z=\infty$, which has a semicircular free density of states RefAG:8

[TABLE]

$W$ is the half-bandwidth which is set as the energy unit. We used the standard DMFT self-consistent equations for the Bethe lattice with $z=\infty$ in the antiferromagnetic phase RefAG:8 .

Fig.4 shows the LDOS of the half-filled Hubbard model at a low temperature $T=0.01$, obtained from SC (Fig.4(a)) and ISC (Fig.4(b)), respectively. At $U=0$, the LDOS’s from SC and ISC are both exact, corresponding to a paramagnetic metal phase. With gradual increase of $U$, distinctions appear in two aspects. First, the LDOS from SC opens an energy gap only at intermediate $U$ while that from ISC opens a gap almost immediately as $U>0$. This contrast is seen clearly in the $U=1.0$ curves: the SC curve is clearly a metal while the ISC curve has a dip at $\omega=0$, being more close to a true insulator. The ISC results agree better with the fact that the ground state of Hubbard model on a bipartite lattice is a Slater insulator in the weak coupling regime, due to the Fermi surface nesting RefAG:8 ; RefTP:27 . We checked that the small finite $\rho_{\uparrow}(\omega=0)$ in the ISC curve of $U=1.0$ is not due to finite temperature, but due to the ISC approximation.

In the second aspect, for intermediate to large $U$ regimes, SC only gives out very weak antiferromagnetism as shown by the weakly asymmetric shape of the LDOS. The antiferromagnetism disappears for $U\geqslant 4.0$ in Fig.4(a). In contrast, in Fig.4(b), ISC gives much more robust antiferromagnetism with a strongly asymmetric LDOS and a sharp peak on the shoulder of the upper sub-band. The antiferromagnetism persists to the large $U$ limit. This shows that ISC correctly produces the antiferromagnetic ground state of this system in the large $U$ regime. The Néel temperature $T_{N}$ obtained from ISC varies with $U$ in the same qualitative trend as the QMC data. However, both the ground state order parameter and the magnitude of Néel temperature are overestimated by the ISC solver, compared to the QMC results RefQMC:28 . This observation is consistent with the poor description of the Kondo screening by ISC since the unscreened local moments order easier and generate stronger magnetism.

5 Discussions and Summary

Besides the improved self-consistent calculation of the averages of the type $\langle Ac\rangle$ in $G_{1}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma\dagger})_{\omega}$, we also tried to improve the truncation approximation Eq.(3.2) for the averages of the types $\langle Acc\rangle$ and $\langle Ac^{\dagger}c\rangle$ in $G_{2}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma\dagger})_{\omega}$. For an example, for $\langle A_{\lambda\gamma}c_{p\sigma^{\prime\prime}}c_{k\sigma^{\prime}}\rangle$, we use \langle A_{\lambda\gamma}c_{p\sigma^{\prime\prime}}c_{k\sigma^{\prime}}\rangle=\Big{\langle}G(c_{k\sigma^{\prime}}|A_{\lambda\gamma}c_{p\sigma^{\prime\prime}})_{\omega}\Big{\rangle} and expand the latter to first-order of $V_{k\sigma}$, generating

[TABLE]

The averages of the type $\langle Ac\rangle$ in the second equation are calculated self-consistently by scheme Eq.(24). Compared to the truncation approximation Eq.(3.2) which is exact at $V_{k\sigma}^{0}$, the above approximation is exact at $V_{k\sigma}^{2}$ and partially takes into account the spin exchange between impurity and the bath electrons. Replacing Eq.(3.2) with this approximation, we obtain the improved $G_{2}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma\dagger})_{\omega}$ which modifies the expression of functions $F_{\alpha\beta,\gamma\delta}^{\sigma}(\omega)$ in $W_{3\sim 5}^{\sigma}(\omega)$.

However, we find that this improvement, together with the improvement in calculating $\langle Ac\rangle$ (Eq.(24)) and in the CF resummation (Eq.(3.3)), does not produce uniformly improved results over ISC. This may be attributed to the fact that the Kondo singlet are formed by degenerate many-body states and is hence singular at $V_{k\sigma}=0$. Any expansion from the atomic limit will break the spin exchange process. To handle the Kondo screening accurately, a truly non-perturbative calculation of the averages of the types $\langle Acc\rangle$ and $\langle Ac^{\dagger}c\rangle$ are necessary. We note that the decoupling approximation of the Lacroix type RefCL:11 , such as $\langle AcAc\rangle\approx\langle Ac\rangle\langle Ac\rangle$ was proved to be able to describe the Kondo effect. Exploration in this direction will be carried out in the forthcoming work.

One of the motivations of the present study is to develop an impurity solver for the multi-orbital AIM, which is of more practical interest in the community of DMFT. An accurate and fast multi-orbital impurity solver is valuable for DMFT study of strongly correlated $d$ or $f$ electron materials. Most of our formula used in this paper can be generalized to multi-orbital AIM, due to the generality of the SBO formalism. However, the CF resummation will become more complicated for multi-orbital AIM because the zero-th order GF contains larger number of poles instead of two in the single-orbital case. Consequently, a multi-level extended CF must be used to do the resummation. This makes the analytical derivation more difficult, if not impossible at all. Therefore, a computer-aided CF resummation or the matrix generalization of the present strong-coupling expansion theory are required for this purpose. Note that our formalism, though demonstrated in this paper only for the particle-hole symmetric case, can be generalized to particle-hole asymmetric case without problem.

Besides the EOM-based series expansion of GF used in this paper, the Mori-Zwanzig CF formalism RefMZ:19 ; RefLee is also a useful framework to carry out systematic strong-coupling expansion, as was done for Hubbard model on the Bethe lattice in Ref.RefHong . In the present work, the GF is first expanded into a series of a small parameter and an extended CF is then constructed from the truncated series. Ref.RefHong used an inverse procedure, i.e., one first write down the exact formal expression for the CF via the recursive relation, and then the coefficients of this CF are expanded into a series and truncated to a given order. It thus avoids the difficulty of mapping GF to CF in our approach. The exact formal expression for $G(d_{\sigma}|d_{\sigma}^{\dagger})_{\omega}$ Eq.(3.3) obtained from the Mori-Zwanzig formalism suggests that this approach has certain advantages and may provide an alternative to the present approach. In both approaches, the self-consistent calculation of averages that encode the important ground state correlation is unavoidable. Quantitative comparison of the two approaches are under study.

In summary, we improve the second-order EOM-based strong-coupling perturbation theory for the symmetric AIM, which was first developed in I. The improvement consists of two aspects. First, the averages of the type $\langle Ac\rangle$ in $G_{1}(A_{\alpha\beta}^{\sigma}|A_{\gamma\delta}^{\sigma{\dagger}})_{\omega}$ is calculated from an exact relation instead of a $V_{k\sigma}^{1}$ order approximation. Second, contributions from $V_{k\sigma}^{n\geqslant 3}$ are retained in the CF coefficients $a_{1}$ which now assumes the exact form. The averages are calculated self-consistently through the CF-resummed GF $G_{CF}(d_{\sigma}|d_{\sigma}^{{\dagger}})_{\omega}$. Using NRG results as reference, we show that these modifications significantly improve the response of the impurity spin to the bath bias. Combined with DMFT, improved description is obtained for the antiferromagnetic phase of the half-filled Hubbard model on Bethe lattice. We also present an extended CF formalism for the GF of AIM, which may be the starting point of an alternative strong-coupling expansion theory.

6 Acknowledgements

This work is supported by 973 Program of China (2012CB921704), NSFC grant (11374362), Fundamental Research Funds for the Central Universities, and the Research Funds of Renmin University of China 15XNLQ03.

Appendix A Appendix: expressions for $W_{1}^{\sigma}\sim W_{5}^{\sigma}(\omega)$

In Appendix A, we give the explicit expression for $W_{1}^{\sigma}\sim W_{5}^{\sigma}(\omega)$ in Eq.(29). To obtain them, we follow Eqs.(84), (88)-(91) of I and insert the self-consistent calculation schemes Eq.(23) and (24) into them. We also multiplied $\left(\omega+E_{\gamma}-E_{\delta}\right)$ to $K^{\sigma}_{\alpha\beta,\gamma\delta}(\omega)$ to take care of the Lehmann representability of $W_{i}^{\sigma}(\omega)$ ($i=3,4,5$). At particle-hole symmetric point, we obtain for spin up,

[TABLE]

For spin down, we obtain

[TABLE]

The functions $K^{\sigma}_{\alpha\beta,\gamma\delta}(\omega)$ in the above equations are given by

[TABLE]

The functions $F^{\sigma}_{\alpha\beta,\gamma\delta}(\omega)$ are given by

[TABLE]

and

[TABLE]

Here, $I_{\alpha\beta}\equiv\langle A_{\alpha\alpha}\rangle+\langle A_{\beta\beta}\rangle$. The particle-hole symmetry gives the relations $\Phi_{13}^{\uparrow}(\omega)=\Phi_{41}^{\downarrow}(-\omega)$ and $\Phi_{42}^{\uparrow}(\omega)=-\Phi_{23}^{\downarrow}(-\omega)$.

Appendix B Appendix: extended CF formula for $G(d_{\sigma}|d_{\sigma}^{\dagger})_{\omega}$

In Appendix B, we derive the formally exact extended CF formula Eq.(34) for $G(d_{\sigma}|d_{\sigma}^{\dagger})_{\omega}$ of AIM, using the Mori-Zwanzig formalism.

We consider the Heisenberg equation of motion for the operator $A_{0}(t)=e^{iHt}A_{0}e^{-iHt}$ in the Heisenberg picture,

[TABLE]

Using the Liouville operator $LO\equiv[H,O]$, we get

[TABLE]

Here and below, $O(t)$, $O(z)$, and $O$ are used to denote the Heisenberg operator, its Laplace transformation, and $O(t=0)$, respectively. We introduce the inner product $(A|B)$ of the operator $A$ and $B$ and define the projection operator $P_{1}$ as

[TABLE]

$P_{1}$ is a superoperator that projects any operator into the one-dimensional subspace of $A_{0}$. The projection operator of the orthogonal subspace is $Q_{1}=1-P_{1}$. We have the relations $P_{1}^{2}=P_{1}$, $Q_{1}^{2}=Q_{1}$, and $P_{1}Q_{1}=Q_{1}P_{1}=0$. In this work, we use the inner product $(A|B)=\langle\{A^{\dagger},B\}\rangle$, which guarantees $P_{1}=P_{1}^{\dagger}$, $Q_{1}=Q_{1}^{\dagger}$, and $L=L^{\dagger}$. The generalized Langevin equation for $A_{0}(t)$ reads

[TABLE]

It is an alternative formulation of the Heisenberg equation of motion Eq.(45). $\Omega_{0}=(A_{0}|[A_{0},H])/(A_{0}|A_{0})$ is the frequency and the random force is given by $A_{1}(t)=e^{iQ_{1}Lt}A_{1}$ with $A_{1}=Q_{1}LA_{0}$. $A_{1}(t)$ is orthogonal to $A_{0}$. The memory function reads $K_{1}(t)=(A_{1}|A_{1}(t))/(A_{0}|A_{0})$. The derivation of Eq.(48) can be found in Ref.RefMZ:19 where the use of the following Dyson identity is made

[TABLE]

Applying the Laplace transformation $f(z)=\int_{0}^{\infty}f(t)exp(-zt)dt$ to Eq.(46) and Eq.(48), we obtain respectively the Heisenberg equation of motion and the generalized Langevin equation on the complex variable $z$ axis

[TABLE]

Projecting the second equation of Eq.(B) to $A_{0}$, we obtain

[TABLE]

The Fermion-type Matsubara Green’s function $G(A_{0}|A_{0}^{\dagger})_{i\omega_{n}}$ is obtained from $K_{0}(z)$ via $G(A_{0}|A_{0}^{\dagger})_{i\omega_{n}}=-iK_{0}(z=\omega_{n})$.

To calculate $K_{1}(z)$, one needs to obtain $A_{1}(t)$. It satisfies the equation of motion and the (equivalent) generalized Langevin equation,

[TABLE]

In Eq.(53), $\Omega_{1}=(A_{1}|[A_{1},H]))/(A_{1}|A_{1})$ and $A_{2}(t)=e^{iQ_{2}Q_{1}Lt}A_{2}$ with $A_{2}=Q_{2}Q_{1}LA_{1}$. The memory function reads $K_{2}(t)=(A_{2}|A_{2}(t))/(A_{1}|A_{1})$. The projection operator $Q_{2}=1-P_{2}$ and $P_{2}O=\left[(A_{1}|O)/(A_{1}|A_{1})\right]A_{1}$. Hence $A_{2}(t)$ is orthogonal to both $A_{1}$ and $A_{0}$. The Laplace transformation to Eqs.(52) and (53) reads

[TABLE]

and

[TABLE]

respectively. Projecting Eq.(55) to $A_{1}$ gives

[TABLE]

Repeating the above derivation, we will obtain the CF for $K_{0}(z)$, as did by MoriRefMZ:19 , LeeRefLee , TserkovnikovRefTserkovnikov , etc.

The specialty of the AIM is that some components in $A_{i}(z)$ ($i=1,2,...$), such as the bath operators $c_{k\sigma}(z)$, can be solved exactly from their Heisenberg EOM. A $z$-dependent term can thus be separated from the memory function, leading to an extended CF for $K_{0}(z)$. To employ this properties of the impurity model, we split the operator $A_{1}(z)$ as $A_{1}(z)=B_{1}(z)+C_{1}(z)$. Both $B_{1}(z)$ and $C_{1}(z)$ have their own equation of motion and the generalized Langevin equation, similar to those of $A_{1}(z)$ in Eqs.(54) and (55). As to be shown below, from $A_{1}(z)$ we can collect the bath operators $c_{k\sigma}(z)$ into $B_{1}(z)$ and solve its EOM exactly. It provides the frequency-dependent term in the extended CF. We need only to apply the generalized Langevin equation for $C_{1}(z)$ alone.

To guarantee the orthogonality between different hierarchy of operators, we also need $(A_{0}|B_{1})=(A_{0}|C_{1})=0$ to be fulfilled. The time-evolution is given by $B_{1}(t)=e^{iQ_{1}Lt}B_{1}$ and $C_{1}(t)=e^{iQ_{1}Lt}C_{1}$, similar to that of $A_{1}(t)$. We then have

[TABLE]

with

[TABLE]

Once $B_{1}(z)$ is solved exactly from its equation of motion $zB_{1}(z)=B_{1}+i\left[Q_{1}LB_{1}\right](z)$, it is easy to calculate $(B_{1}|B_{1}(z))$ and $(C_{1}|B_{1}(z))$. $(B_{1}|C_{1}(z))$ can be calculated from $(B_{1}|C_{1}(z))=\int_{0}^{\infty}dte^{-zt}(B_{1}(-t)|C_{1})$. So $K_{1B}(z)$ is easily calculated. For $K_{1C}(z)$, we write down the generalized Langevin equation for $C_{1}(z)$ as

[TABLE]

which gives

[TABLE]

Here, the frequency and memory functions are given by

[TABLE]

The second-order memory function $K_{2}(z)$ is defined by the random force $A_{2}(z)$ of $C_{1}(z)$. $A_{2}(z)$ is the Laplace transformation of $A_{2}(t)=e^{iQ_{2}Q_{1}Lt}A_{2}$, with $A_{2}=Q_{2}Q_{1}LC_{1}$. Here $Q_{2}=1-P_{2}$ and $P_{2}O=[(C_{1}|O)/(C_{1}|C_{1})]C_{1}$. Therefore, $A_{2}(z)$ is orthogonal to $C_{1}$.

The same procedure can be implemented for $A_{2}(z)$, i.e., from $A_{2}(z)$ we collect the bath operators into $B_{2}(z)$ and do the splitting $A_{2}(z)=B_{2}(z)+C_{2}(z)$. We apply the equation of motion and generalized Langevin equation to $B_{2}(z)$ and $C_{2}(z)$, respectively. Employing the exactly solved $B_{2}(z)$, we obtain the second-order memory function $K_{2}(z)=K_{2B}(z)+K_{2C}(z)$. $K_{2B}(z)$ is an exactly solvable frequency-dependent term

[TABLE]

and $K_{2C}(z)$ is an unknown memory function

[TABLE]

Here, we also require that $(C_{1}|B_{2})=(C_{1}|C_{2})=0$. The time evolution is given by $B_{2}(t)=e^{iQ_{2}Q_{1}Lt}B_{2}$ and $C_{2}(t)=e^{iQ_{2}Q_{1}Lt}C_{2}$. By repeating the above process, we can express $K_{0}(z)$ into an infinite extended CF.

For the Hamiltonian in Eq.(12), we start from $A_{0}=d_{\sigma}$ and obtain the following operators through straightforward calculation,

[TABLE]

Here $\tilde{\beta}_{\sigma}$ is given in Eq.(35). They fulfil the orthogonal requirements and the EOM of $B_{1}(z)$ and $B_{2}(z)$ can be solved exactly. From these operators, we obtain the projecting coefficients used in the extended CF

[TABLE]

Finally, putting Eq.(65) into Eqs.(57), (58), (60), (62), and (63), we obtain the two-level extended CF Eq.(34).