Peng Fan(范鵬)and Ning-Hua Tong(同寧華)

Department of Physics,Renmin University of China,Beijing 100872,China

1.Introduction

Green’s function(GF)is a widely used tool in the study of quantum many-body physics.Among many methods for calculating GF,the equation of motion(EOM)approach to the two-time GF is based on the Heisenberg equation of motion of operators.[1–4]For a given interacting Hamiltonian,the EOM of a given GF contains higher order GFs and repeatedly applying the EOM generates a chain of GFs.In the conventional Tyablikov-type truncation approximation,[3]the GF of certain order is approximated as a linear combination of the lower order GFs in the frequency domain.This leads to a set of closed but approximate algebraic equations for the GFs,which can be solved to obtain the desired GFs.

Decades of experience on this practice for various Hamiltonians shows that naive truncation of the EOM has some drawbacks.First,causality of the GFs is not guaranteed.The truncation may destroy the correct analytical structure of GF,i.e.,GF containing only real simple poles.Second,for a given higher order GF,the truncation scheme is not unique.Different truncations may lead to drastically different results.Since there is no transparent clue for the optimal truncation scheme,in practice,the truncation depends heavily on experience.Due to these drawbacks,the EOM truncation approach is usually regarded as an uncontrolled approximation.Its application in the modern study of quantum many-body physics is therefore severely limited.

There are efforts to overcome the drawbacks of the EOM truncation approach.Using the idea of operator projection,Mori[5,6]and Zwanzig[7]developed the generalized Langevin equation formalism for the operator EOM,which can be used to calculate GFs.In particular,an elegant and exact continued fraction formalism was proposed to express GF or time correlation functions in terms of projecting coefficients in the operator space.[6,8,9]Similar theories have been developed by Tserkovnikov[10,11]and applied by many other researchers under different names,including the two pole approximation,[12,13]composite operator approach,[14,15]projection operator approach,[16–20]irreducible GF method,[21]and many others.[22,23]These closely related theories employ the operator projection idea to truncate the EOM.They have the advantage that the causality of GF is guaranteed by the formalism.Also,the time translation invariance of the equilibrium state is strictly obeyed by the GF.This is embodied by the fact that?G(t?t′)/?t and ?G(t?t′)/?t′give equivalent formula.Therefore,the first drawback of the EOM approach is removed.

However,in these theories,except for special cases,[24]the projection coefficients cannot be calculated without introducing additional approximations.The second drawback of the EOM truncation approach,i.e.,the arbitrariness in the truncation,is still present in these theories.In our recent work,[25]we proposed a practical and systematic projective truncation approximation(PTA)for the EOM of GFs.For a selected set of operator basis{A1,A2,...,An},our theory is equivalent to the matrix form of the Mori–Zwanzig formula for GFs,with the memory function matrix neglected.By using a partial projection approximation,we reduce the calculation of projecting coefficients into that of two matrices,the inner product matrix III and the natural closure matrix MMM.They are calculated self-consistently by the fluctuation-dissipation theorem(for III )and through the commutators[Ai,H](i=1,2,...,n)(for MMM),respectively.The arbitrariness in the truncation is thus removed.For the Anderson impurity model,our PTA at the same truncation level is superior to the conventional Lacroix approximation[26]and the results are in quantitative agreement with those of the numerical renormalization group(NRG).[27]

In this paper, we study the convergence properties of PTA,by comparing the results of PTA on three successively larger bases.Using the Anderson impurity model and NRG results as reference,we examine whether the PTA results are improved with enlarging basis size and converge towards the exact ones.The positive results of this convergence check establishes PTA as a method of controlled precision for quantum many-body systems.

2.The projective truncation approximation for EOM

The projective truncation approximation for the EOM of GFs was developed in Ref.[25].In this section,for the sake of completeness,we overview the general formalism.

For a given Hamiltonian H,we choose n linearly independent operators{A1,A2,...,An}to span a subspace of the full Liouville space.To truncate the higher order operators generated by EOM,we project them into this subspace and neglect the component orthogonal to it.Therefore,the basis set should contain the most important excitations of the system.Such projective truncation becomes exact if the basis is complete.

For the column vector AAA formed by basis operators,the retarded GF matrix is defined as

where θ(t?t′)is the Heaviside step function and AAA(t)is the vector of basis operators in Heisenberg picture.In this paper we only consider the Fermion-type GF and the curly bracket in the above equation denotes anti-commutator.Below we take the natural unit and drop?.

The equation of motion for the GF matrix in the frequency domain reads

We define the inner product of operators X and Y as

where〈〉=Tr(ρ)and ρ=e?βH/Tr(e?βH)is the equilibrium density operator of H at temperature T.Note that in our theory,the Hamiltonian or density operator is not projected.All the averages are to be calculated with a Gibbs distribution in the full Hilbert space.Equation(4)fulfils the standard requirements for the inner product in a linear space.Other definitions of inner product can be found in the literature as well.[5,6,22]

We write the commutator between the basis operators and the Hamiltonian as

The first term on the right-hand side includes all the basis operators that naturally appear in the commutator.MMM is called the natural closure matrix.Biis the newly generated operator outside the basis(In certain situations,Bimay include basis operators for symmetry reasons,see below.).We further decompose

where δBiis the component orthogonal to the subspace of basis operators.That is,(Ak|δBi)=0 for i,k=1,2,...,n. NNN is obtained by projecting Eq.(6)onto the basis operators and solving the obtained equation

with PPP and III defined as PPPij≡(Ai|Bj)and IIIij≡(Ai|Aj),respectively.

Putting Eq.(6)into Eq.(5)and projecting it onto the basis operators,we obtain

Here,MMMt=MMM+ NNN is the total closure matrix.The Liouville matrix LLL is defined as LLLij≡(Ai|[Aj,H]).Note that III is Hermitian and positive definite. LLL is Hermitian under the inner product Eq.(4),which guarantees the causality of GF.

The general idea of projective truncation[5–10,12,13,22]is to neglect the orthogonal component δBiin Eq.(6),i.e.,

Putting Eqs.(5),(6),and(9)into the EOMs of GF,Eqs.(2)and(3),we obtain

Note that the left side and the right side time derivatives of GF,Eqs.(2)and(3),produce the same equation,respecting the time translational invariance of the equilibrium state.Equation(10)is equivalent to the matrix form of the Mori–Zwanzig equation with the memory function matrix neglected.It becomes exact as the operator basis covers the complete Liouville space.

For given matrices MMMtand III,the GF matrix in Eq.(10)can be calculated directly by matrix inversionas done in previous analytical studies,[12,13]or by numerically solving the generalized eigen-value problem of the pair of matrices( LLL, III ),

Here Λ =diag{λ1,λ2,...,λn}is a real diagonal matrix. UUU is the generalized eigen vector matrix which diagonalizes MMMt,UUU?1MMMtUUU= Λ.It fulfills the generalized orthonormal relation UUU?III UUU =1.Equation(10)can be reformulated in terms of UUU and Λ as

The corresponding spectral function reads

The Hermitian matrices L and I contain the average of operators on the state defined by the density matrix ρ in Eq.(4).The calculation of them usually relies on additional approximations which cause the arbitrariness.In Ref.[25],we proposed a practical and systematic method to calculate the matrices I and L.The averages of the kind 〈A?jAi〉can be obtained from the corresponding GF via the spectral theorem as

For the average of the type〈Ai〉(is an operator outside the basis set{Ak}),we use

〉can then be calculated self-consistently from III, UUU,and Λ,provided that the averages〈{Ai,O?}〉(i=1,2,...,n)are linear combinations of{〈A?kAp〉}.[12,13]

Usually equations(14)and(15)are sufficient to produce III but not L.Therefore,we introduce the following partial projection approximation for L.We first classify the basis operators into two groups,The superscripts 1 and 2 denote subset-1 and 2,respectively.Subset-1 is composed of basis operators whose commutators with H close automatically.Subset-2 contains the rest basis operators.That is,we have=0(i=1,2,...,m)and=0(j=m+1,m+2,...,n).Associated with this split of basis,the matrices I ,M, P, N,and L all become 2×2 block matrices.In particular,where PPP12and PPP22are the projection matrices from BBB(2)to AAA(1)and to AAA(2),respectively.Similarly,we have

Employing the Hermiticity of L,we proposed the following partial projection approximation to L,

Here,

and we use the exact expression for PPP12,

Under this approximation,the input of the calculation are MMM and II

I matrices only.The precision of results is determined only by the selection of basis operators.Below we will show that the precision is improved systematically with enlarged basis size.

3.Application to Anderson impurity model:formalism

In this section,we apply the PTA to the Anderson impurity model(AIM).Taking the NRG results as a reference,we compare the results from three successively larger basis sets:the basis at the level of Hubbard-I approximation[28,29](HIA basis),at the level of alloy analogy approximation[30](AAA basis),and at the level of Lacroix approximation[26](Lacroix basis).These bases form a chain of sets:HIA basis?AAA basis?Lacroix basis,so that we can speak of enlarging the basis.Our aim is to study how the results depend on the basis size and to observe the convergence of results to the exact ones in the large basis limit.The Hamiltonian of the AIM that we will study reads

The denotations are standard.For the hybridization function?σ(ω)≡ ∑kV2kσδ(ω ?εkσ),we use the Lorentzian type,

Here,δω is a parameter used to tune the band energy of conduction electrons with different spins.It is equivalent to applying an external magnetic field to the conduction electrons.Due the spin dependence of the band energy,the spin of the conduction electrons could be polarized.This can simulate a quantum dot coupled to ferromagnetic leads and can be used to study the influence of the ferromagnetism on the Kondo effect.[31]Besides,in the dynamical mean- field theory,[32,33]the hybridization function will become spin dependent in the magnetic phase.Equation(22)is used to test the reliability of our impurity solver for the magnetic case,although in the real dynamical mean- field calculation for magnetic phase,the spin dependence of hybridization is more complicated.ωc=1.0 is set as the energy unit.In accordance with?ˉσ(ω)=?σ(?ω),we assume the following constraints in the parameters of H,

The particle–hole symmetry of H is realized at the parameter point εd= ?U/2 and μ =0.

For HIA and AAA bases,the projection matrix PPP is calculated analytically without using the partial projection approximation. NNN and G( AAA| AAA?)ωare obtained by analytically solving the linear equations Eqs.(7)and(10).For AAA basis,additional decoupling approximations are used to calculate some of the averages in PPP.The particle–hole symmetry is fulfilled automatically in these approximations.The results for Lacroix basis are taken from Ref.[25],where we used a particle–hole symmetric form of the partial projection approximation and solved the GF numerically via the generalized eigen-value formalism Eq.(11)on a discretized bath.

3.1.HIA Basis

The conventional HIA for AIM is obtained by truncating the EOM at the second order.The involved operators are selected here to form the HIA basis,

In this basis set,k goes through all nkwave vectors of the conduction electron,giving the HIA basis a dimension of d=2+nk.Due to the conservation of total number of electronsand the z component of total spin Sz,the basis operators are confined to the type that annihilates an electron with spin σ.Here we did not use the full spin SU(2)symmetry of AIM.The inner product matrix of the basis operators is

In Eq.(25),the rank of the matrix corresponds to?the sequenc?e{A1,A2k,A3}and the column corresponds toA1,A2p,A3.For simplicity,the full d×d matrix is abbreviated as a 3×3 matrix,in which the sub-matrix involving bath operators is abbreviated to a k-dependent number.For an example,δkpis used to represent the unity sub-matrix(A2k|A2p)(k,p=1,2,...,nk).Below,our matrix expression will always use this abbreviation convention.

The commutators of the basis operators with H are summarized in Appendix A.The generated new operators are B1=B2k=0 and

The matrices MMM and PPP are written as

and

whereNote that we have included an additional terminto B3to make it particle–hole symmetric.The averages are taken as real numbers,i.e.,〉=〉.

For this simple case,the full projective truncation can be carried out,i.e., PPP (i.e.,〈nσˉ〉and βσ)can be calculated selfconsistently without further approximation.Solving Eqs.(7)and(10)analytically,we obtain the impurity self-energy as

Here=βσ/[〈nσˉ〉(1?〈nσˉ〉)].Equation(29)has the form of atomic limit,same as the conventional HIA,[34]but with an additional spin-dependent shiftβ?σof the impurity level.It is exactly the extended continued fraction expression with the memory function omitted at this level.[35]The impurity GF is obtained fromthe Dyson equationWe also obtainThe non-interacting impurity GF is given by

with Γσ(ω)=R+∞?∞?σ(ε)/(ω?ε)dε.For the averages,〈nˉσ〉can be calculated from G(dˉσ|d?ˉσ)ω.?βσneeds to be calculated from the right-hand side EOMs

For a paramagnetic bath and at the particle–hole symmetric point,=0.Equation(29)recovers that of the conventional HIA.Away from the particle–hole symmetry or for δ ω =0,=0 and equation(29)differs from the conventional HIA.

3.2.AAA basis

The AAA basis is obtained by adding the operators A4k=nˉσckσ(k=1,2,...,nk)into the HIA basis,

The total dimension is d=2+2nk.The Tyablikov-type decoupling of EOM at this level results in an approximation which,when combined with dynamical mean- field theory for Hubbard model,gives the conventional AAA.[29,34]Similar calculation is carried out as for the HIA basis.The inner product matrix reads

Using the commutators summarized in Appendix A,we obtain the matrix MMM as

Commutators of A3and A4kwith H generate the new operators)and

The projection matrix PPP in the block form of Eq.(16)has submatrices PPP11= PPP21= PPP12=0,and

For the AAA basis,it is still possible to analytically solve the linear equations Eqs.(7)and(10)for the local GFs G(dσ|dσ?)ωand G(nσˉdσ|dσ?)ω.However,the averages in( PPP22)3,4pand( PPP22)4k,4pcannot be written in the form 〈A?iAj〉.If we use the partial projection approximation Eq.(18)for LLL,due to PPP12=0,we recover the conventional AAA which is equivalent to setting B(32)=B(42k)≈0.To go beyond the conventional AAA, we use a simple decoupling approximation for the elements of PPP22,Since this approximation keeps PPP22symmetric and IIIM is also symmetric,the Hermiticity of LLL is conserved.The particle–hole symmetry is also fulfilled.

The key problem for the projective truncation of equation of motion of Green’s function is how to evaluate the projection coefficients,such as those in Eq.(36).Besides the partial projection approximation that we introduced in Ref.[25],one could calculate the averages in Eq.(36)from new Green’s functions G(X|Y)ωwhere X and/or Y are certain operators outside the basis.The equation of motion for G(X|Y)ωcan be closed by similar projection truncation.However,except for special situations,[12,13]this process will generate new averages which do not form closed set of equations,as for the case of Eq.(36).In Eq.(37),we therefore perform a decoupling approximation due to Lacroix.[30]It takes into account the hybridization between the conduction electrons and the impurity and captures the main aspects of the Kondo effect.Therefore,for this decoupling approximation the Kondo physics should be captured and the Kondo peak is expected to occur at the Fermi energy in the spectral function(see Fig.5).

Solving Eqs.(7)and(10)analytically,we obtain the selfenergy for AAA basis as

Compared to the self-energy of HIA basis,a frequency dependent shift and broadening of the impurity level appears as

with

Neglecting Bpσand?βσ,equation(38)recovers the conventional AAA.[34]Equation(38)is also consistent with the form of extended continued fraction,[35]but with an approximate expression for the memory function.Using EOM of GF G(dσ|c?pσ)ω,we reduce Eq.(39)to

Here,GRσ(ω)=.The symbol〈g(ω)〉represents(?1/π))dω with η being an infinitesimal positive number.In Eq.(39),the shift ofVkσis generated by projecting B3to A4kand B4kto A3.It contains the spin exchange between impurity and bath electrons.As will be shown below,this renormalization of hybridization produces improved description of the Kondo peak at low temperatures compared to conventional AAA.

3.3.Lacroix basis

In the work of Lacroix,[30]the GFs generated by the commutator of H and nˉσdσare kept and the truncation is done in the next order EOM.Here we take the corresponding operators to form the Lacroix basis{A1,A2k,A3,A4k,A5k,A6k}(k=1,2,...,nk),with

The superscripts 1 and 2 denote the grouping of basis operators according to the closure properties of their commutators with H:B(i1)=0 and B(i2)=0.This basis has a dimension 4nk+2.For this basis,we directly take the results from Ref.[25],where we used the particle–hole symmetric partial projection truncation.We numerically solved the PTA equations for a linearly-discretized bath with nk=401 bath sites(i.e.,the total basis size 1606),which already represents the continuous bath satisfactorily.The half band width is D=5.0.The δpeaks in the LDOS were broadened with η =0.01～ 0.02.

4.Numerical results and comparison

Using the formalism in previous sections,we obtain numerical results for HIA,AAA,as well as Lacroix bases.Below,these approximations are called projective-HIA(pHIA),projective-AAA(pAAA),and projective-Lacroix(pLacroix),respectively.The NRG results,used as a reference,are obtained from the full density matrix algorithm.[36,37]For the local density of states(LDOS),we use the self-energy trick[38]and average on Nz=8 interleaved discretizations.[39]The logarithmic discretization parameter is Λ=2.0 and we keep Ms=350～380 states.Though not extrapolated to the exact limit Λ =1 and Ms= ∞,[40]we have checked that the uncertainties in NRG results are much smaller than the difference between NRG and all the approximate results.For the results below,we fixμ=0.0 and?=0.1.

We first study the impurity electron occupation 〈nσ〉and the double occupancy 〈n↑n↓〉as functions of εd,δω,U,and T.They describe the static magnetic and the charge response of the impurity to external parameters.In Fig.1,we plot〈n↑〉(Fig.1(a))and 〈n↑n↓〉(Fig.1(b))as functions of εd,for a nonmagnetic bath at U=2.0 and T=0.1.The curves of pHIA,pAAA,and pLacroix are compared with those of NRG.As the basis is enlarged from HIA to Lacroix,both quantities shift towards NRG results in the whole εdregime,with slight overshooting in the pLacroix results.At εd= ?1.0,all methods give 〈n↑〉=0.5 due to the particle–hole symmetry.The significant improvement of pLacroix over pHIA and pAAA shows that the hybridization effect lacking in the HIA and AAA bases is important for quantitative accuracy.

Fig.1.(a)〈n↑〉and(b)〈n↑n↓〉as functions of εd.The parameters are U=2.0,T=0.1,δω =0.0.

Fig.2.(a)〈n↑〉and(b)〈n↑n↓〉as functions of δω.The parameters are U=2.0,T=0.1,εd=?U/2.

In Fig.2,we plot〈n↑〉and 〈n↑n↓〉as functions of δω for U=2.0,T=0.1,and εd= ?U/2.The particle–hole symmetry at this parameter set assures the exact〈n↑〉=0.5 at δω =0.Away from δω =0,the deviation 〈n↑〉begins to increase for all bases,but pLacroix gives the smallest deviation.In particular,pLacroix gives the correct sign in the impurity spin response to the bath bias.The double occupancy shown in Fig.2(b)has a weak δω dependency.Again,we observe that the results from projective truncations tend to those of NRG systematically with increasing basis size.

Figure 3 shows the same averages as functions of U for T=0.1,εd= ?U/2,and a negative bath bias δω = ?0.2.At U=0.0,all projection truncations give exact results for〈n↑〉and 〈n↑n↓〉.In Fig.3(a),as U increases from zero,the agreement with the NRG curve is maintained to larger U values for larger basis,up to U=1.5 for pLacroix.In Fig.3(b),pHIA gives significant deviation in 〈n↑n↓〉as soon as U>0.pAAA gives good agreement up toU=1.0.The result from pLacroix is not as accurate as pAAA in the small U regime but the overall agreement, especially in the intermediate to large U regime,is much better.

Fig.3.(a)〈n↑〉and(b)〈n↑n↓〉as functions of U.The parameters are T=0.1,δω =?0.2,and εd=?U/2.

The temperature dependence of the same quantities are shown in Fig.4 for U=2.0,εd= ?U/2,δω = ?0.2.In Fig.4(a),using the particle–hole symmetry properties 〈n↑〉+〈n↓〉=1,we can deduce that the impurity spin polarization is zero at T=∞and it increases as temperature is lowered for all approximations.While pHIA and pAAA give 〈n↑〉<0.5 which has the wrong sign of spin polarization,pLacroix gives correct sign and quantitative agreement with NRG for all temperatures.In Fig.4(b),〈n↑n↓〉decreases with decreasing temperature but a slight increase is observed in both pAAA and pLacroix curves(below T=0.05 for pAAA and T=0.15 for pLacroix).This upturn of double occupancy, also seen in NRG result,reflects the screening of local moment and forming of the Fermi liquid state below the Kondo temperature.It is notable that the double occupancy from pLacroix is very accurate at T=0.Similar pattern of convergence is observed in these data as the basis is enlarged from HIA to Lacroix.

Fig.4.(a)〈n↑〉and(b)〈n↑n↓〉as functions of temperature.The parameters areU=2.0,δω =?0.2,εd=?U/2.

Fig.5.Impurity density of states calculated at U=2.0,T=0.001,δω =0.0,εd=?0.7.

Figure 5 shows the LDOS obtained from different bases at a particle–hole asymmetric point εd= ?0.7,U=2.0,and at low temperature T=0.001.They are compared with NRG result.The LDOS from various PTAs correctly contain the upper and lower Hubbard peaks.As the basis is enlarged from HIA to Lacroix,both the position and the weight of the Hubbard peaks tend to those of NRG systematically.At ω=0,pHIA does not produce a Kondo peak while pAAA and pLacroix produce a sharp Kondo peak.Quantitatively comparing the weight and the shape of Kondo peaks,pAAA gives a too much sharper peak with small weight while pLacroix produces a slightly broader peak with closer weight to NRG.The overall tendency of the convergence in LDOS is apparent.

In PTA,using the inner product,we can quantify the truncation error.Here we take the HIA basis as an example,in which the only truncation approximation is δB3≈ 0.For pHIA,we propose the following quantity to measure the truncation error,

where|X|≡(X|X)1/2is the norm of operator X under the inner product Eq.(4).B3and δB3are given by Eq.(26)and Eq.(6),respectively.Extension of the defination of eTto more complicated truncations is possible.The Pythagorean theorem implies that 0≤eT≤1.At finite temperature,eT=0 is equivalent to δB3=0.In this case no approximation is made and PTA becomes exact.In the other limit,eT=1 means B3≈0,a complete negligence of the new operators produced by EOM.Therefore,eTis a quantitative gauge of the truncation error in PTA.In Fig.6(a),we plot eTas a function of εdfor U=2.0,T=0.1,and δω =0.0.We use pLacroix to calculate all the inner product involved in Eq.(43), including the projecting matrix N of Eq.(6).For comparison,in Fig.6(b),we show the particle–hole symmetrized double occupancies from NRG and pHIA as functions of εd.

Fig.6.(a)relative truncation error eTin pHIA,(b)symmetrized double occupancy,as functions of εd.Parameters are U=2.0,T=0.1,and δω=0.0.

It is seen from Fig.6 that eTqualitatively reflects the errors in the physical quantities.At the particle–hole symmetric point εd= ?1.0,eT=1.0,being consistent with the fact that at this point,pHIA recovers the conventional HIA which amounts to B3≈0.In the regime?2.0≤εd≤0.0,eTstays close to 1.0,showing that pHIA has the largest truncation error in this regime.Correspondingly,the discrepancy in the symmetrized double occupancy between NRG and pHIA is significant in this regime.Further away from this regime,eTquickly decreases and in Fig.6(b),the double occupancies from pHIA and NRG merge.In the limits εd= ±∞,〈nσ〉=0 or 1,we expect eT=0.pHIA will become exact since no electron correlation is present in those limits.Figure 6 shows a qualitative correlation between eTand the error in the double occupancy of pHIA.We observe similar correlations in other quantities as well,but there is no strict monotonous correspondence between eTand the errors.This is because in PTA,the physical quantities depend on the truncation error non-linearly.Therefore,we conclude that eTcan be used to gauge the overall level of approximation of PTA.It is especially suitable for comparing the error among different parameter regimes.

5.Discussion and summary

First,let us discuss the scaling of error with the basis size.We did not study this scaling quantitatively because for AIM with a continuous bath,it is difficult to quantify the size of basis.In fact,for pHIA and pAAA,the dimension of the basis is infinity if we count the number of linearly independent operators in the basis,due to infinitely many bath modes in AIM.Therefore,below we only give a qualitative discussion of this issue.If we regard truncating the EOM of GFs as such a problem:for an original operator A in the full operator spaceS,look for the operator A′in a subspace S′?S and require that A′is as close to A as possible.The solution A′will be the projection of A intoS′.In this sense,for a given inner product,projective truncation is the optimal way of truncating the chain of EOMs.This is why the PTA could be superior to conventional decoupling in accuracy.However,in PTA,the inner product(X|Y)is not calculate exactly(except for those trivial ones such as(dσ|dσ))but self-consistently from the GFs obtained in this theory.As a result,even without partial projection approximation,PTA has two sources of error,Liouville space truncation and the approximate evaluation of projection.As the basis is enlarged,both the space truncation and the precision of projection are improved.Thus we expect that the accuracy in the final result improves beyond linear fashion with the basis size.

With enlarging basis,the rate of convergence in results depends crucially on the basis selection method.Here,we simply collect those separate operators appearing in the successive EOMs,[A1,H],[[A1,H],H],etc.There are other ways of selecting basis operators.[41]Especially,the selection of orthogonal basis operators in the Krylov subspace produces a continued fraction form for the local GF.[8,9]The self-energy functions from pHIA and pAAA in this work have the extended continued fraction form[35]used to do resummation for the strong-coupling series expansions of GF.[42]The efficiency of basis,measured by the accuracy versus basis size scaling,depends on how fast the key excitations are taken into account as the basis is enlarged.Finding the optimal basis selection procedure is an important research topic for the future.

Up to now,the best results that we obtain for AIM is from pLacroix.From Fig.2 and Fig.3,it is clear that even for pLacroix,the accuracy in the impurity spin polarization under the bath bias is not satisfactory,especially in the large U and large δω regime which is important for describing the antiferromagnetic phase in Hubbard model within DMFT.Therefore,it is natural to go beyond Lacroix basis.However,we find that it is not easy to maintain the positive definiteness of the inner product matrix III for larger basis set,possibly due to the inclusion of basis operators which has very small norm.Method such as singular value decomposition is being considered to removed those excitation modes of tiny norm.

In summary,in this paper we compare the PTA results obtained from three successively larger bases with those from NRG.The results improve systematically with increasing basis size and a clear tendency of convergence to NRG results is observed.We also propose a quantity to gauge the truncation error in PTA and demonstrate its usefulness in pHIA.Our results confirm that the PTA is a computational method of controllable precision for quantum many-body systems.

Appendix A:Commutators of[AAii,,HH]]

In this appendix,we summarize the commutators between the basis operators and the AIM Hamiltonian H Eq.(21).Below,we give the commutators for HIA and AAA bases.For Lacroix basis,see Appendix B of Ref.[25].

[1]Martin P C and Schwinger J 1959 Phys.Rev.115 1342

[2]Bogolyubov N and Tyablikov S V 1959 Doklady.Akad.Nauk SSSR 126 53

[3]Tyablikov S V,1959 Vkrain.Mat.Zhur.11 287

[4]Zubarev D N 1960 Usp.Fiz.Nauk 71 71

[5]Mori H 1965 Prog.Theor.Phys.33 423

[6]Mori H 1965 Prog.Theor.Phys.34 399

[7]Zwanzig R 1961 Lectures in Theoretical Physics,Vol.3(New York:Interscience)

[8]Lee M H 1982 Phys.Rev.Lett.49 1072

[9]Lee M H 2000 Phys.Rev.E 62 1769

[10]Tserkovnikov Yu A 1981 Theor.Math.Phys.49 993

[11]Tserkovnikov Yu A 1999 Theor.Math.Phys.118 85

[12]Roth L M 1968 Phys.Rev.Lett.20 1431

[13]Roth L M 1969 Phys.Rev.184 451

[14]Mancini F and Avella A 2004 Adv.Phys.53 537

[15]Avella A 2014 Adv.Conden.Matter Phys.2014 515698

[16]Kakehashi Y and Fulde P 2004 Phys.Rev.B 70 195102

[17]Fulde P 1995 Electron Correlations in Molecules and Solids,3rd edn.(Berlin,Heidelberg,New York:Springer-Verlag)

[18]Onoda S and Imdada M 2001 J.Phys.Soc.Jpn.70 632

[19]Onoda S and Imdada M 2001 J.Phys.Soc.Jpn.70 3398

[20]Onoda S and Imada M 2003 Phys.Rev.B 67 161102(R)

[21]Kuzemsky A L 2002 Rivista Nuovo Cimento 25 1

[22]Rowe D J 1968 Rev.Mod.Phys.40 153

[23]Plakida N M 2012 “Strongly Correlated Systems”ed.Avella A and Mancini F,Springer Series in Solid-State Sciences,Vol.171(Berlin,Heidelberg:Springer)

[24]Lee M H,Hong J,and Florencio J 1987 Phys.Scr.T19 498

[25]Fan P,Yang K,Ma K H and Tong N H 2018 Phys.Rev.B 97 165140

[26]Lacroix C 1981 J.Phys.F:Metal Phys.11 2389

[27]Wilson K G 1975 Rev.Mod.Phys.47 773

[28]Hubbard J 1963 Proc.Roy.Soc.London Ser.A 276 238

[29]Hubbard J 1963 Proc.Roy.Soc.London Ser.A 277 237

[30]Hubbard J 1964 Proc.Roy.Soc.London Ser.A 281 401

[31]Sindel M,Borda L,Martinek J,Bulla R,K?nig J,Sch?n G,Maekawa S and von Delft J 2007 Phys.Rev.B 76 045321

[32]Metzner W and Vollhardt D 1989 Phys.Rev.Lett.62 324

[33]Georges A,Kotliar G,Krauth W and Rozenberg M J 1996 Rev.Mod.Phys.68 13

[34]Gebhard F 1997 The Mott Metal-Insulator Transition(Berlin,Heidelberg,New York:Springer-Verlag)

[35]Ma K H and Tong N H arXiv:1901.11471(unpublished)

[36]Weichselbaum A and von Delft J 2007 Phys.Rev.Lett.99 076402

[37]Peters R,Pruschke T and Anders F B 2006 Phys.Rev.B 74 245114

[38]Bulla R,Hewson A C and Pruschke Th 1998 J.Phys.:Condens.Matter 10 8365

[39]Yoshida M,Whitaker M A and Oliveira L N 1990 Phys.Rev.B 41 9403

[40]The convergence of the NRG result for LDOS of AIM with respect to Λ and Mswas discussed in Supplementary Materials of Li Z H et al.2012 Phys.Rev.Lett.109 266403

[41]Ciolo A D and Avella A 2018 Physica B 536 687

[42]Tong N H 2015 Phys.Rev.B 92 165126