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 and 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 , 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 and the natural closure matrix . They are calculated self-consistently by the fluctuation-dissipation theorem (for ) and through the commutators ( ) (for ), 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.

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 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 formed by basis operators, the retarded GF matrix is defined as

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

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

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

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

For given matrices and , the GF matrix in Eq. (10) can be calculated directly by matrix inversion , as done in previous analytical studies,^{[12, 13]} or by numerically solving the generalized eigen-value problem of the pair of matrices ,

The Hermitian matrices and 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 and . The averages of the kind can be obtained from the corresponding GF via the spectral theorem as

Usually equations (14) and (15) are sufficient to produce but not . Therefore, we introduce the following partial projection approximation for . 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 ( ) and ( ). Associated with this split of basis, the matrices , , , , and all become 2 × 2 block matrices. In particular,

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

For HIA and AAA bases, the projection matrix is calculated analytically without using the partial projection approximation. and 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 . 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 n_k wave vectors of the conduction electron, giving the HIA basis a dimension of . Due to the conservation of total number of electrons and the z component of total spin S_z, 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 isIn Eq. (25), the rank of the matrix corresponds to the sequence and the column corresponds to . 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, is used to represent the unity sub-matrix ( ). 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 and

The matrices and are written asandwhere . Note that we have included an additional term into B₃ to 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., (i.e., and ) can be calculated self-consistently without further approximation. Solving Eqs. (7) and (10) analytically, we obtain the impurity self-energy as

Here . 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 from the Dyson equation . We also obtain . The non-interacting impurity GF is given bywith . For the averages, can be calculated from . needs to be calculated from the right-hand side EOMsFor a paramagnetic bath and at the particle–hole symmetric point, . Equation (29) recovers that of the conventional HIA. Away from the particle–hole symmetry or for , and equation (29) differs from the conventional HIA.

3.2. AAA basis

The AAA basis is obtained by adding the operators ( ) into the HIA basis,

The total dimension is . 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 readsUsing the commutators summarized in Appendix A, we obtain the matrix asCommutators of A₃ and with H generate the new operators and ,The projection matrix in the block form of Eq. (16) has sub-matrices , and

For the AAA basis, it is still possible to analytically solve the linear equations Eqs. (7) and (10) for the local GFs and . However, the averages in and cannot be written in the form . If we use the partial projection approximation Eq. (18) for , due to , we recover the conventional AAA which is equivalent to setting . To go beyond the conventional AAA, we use a simple decoupling approximation for the elements of ,

Since this approximation keeps symmetric and is also symmetric, the Hermiticity of 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 where X and/or Y are certain operators outside the basis. The equation of motion for 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 self-energy for AAA basis as

Compared to the self-energy of HIA basis, a frequency-dependent shift and broadening of the impurity level appears aswithNeglecting 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 , we reduce Eq. (39) toHere, . The symbol represents with (η being an infinitesimal positive number. In Eq. (39), the shift of is generated by projecting B₃ to and to A₃. 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 are kept and the truncation is done in the next order EOM. Here we take the corresponding operators to form the Lacroix basis ( ), with

The superscripts 1 and 2 denote the grouping of basis operators according to the closure properties of their commutators with H: and . This basis has a dimension . 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 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 .

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 interleaved discretizations.^[39] The logarithmic discretization parameter is Λ=2.0 and we keep states. Though not extrapolated to the exact limit Λ=1 and ,^[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 and the double occupancy as functions of , , U, and T. They describe the static magnetic and the charge response of the impurity to external parameters. In Fig. 1, we plot (Fig. 1(a)) and (Fig. 1(b)) as functions of , for a non-magnetic 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 regime, with slight overshooting in the pLacroix results. At , all methods give 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.

	Figure Option ViewDownloadNew Window
	Fig. 1. (a) and (b) as functions of . The parameters are U=2.0, T=0.1, δ ω = 0.0.

In Fig. 2, we plot and as functions of for U=2.0, T=0.1, and . The particle–hole symmetry at this parameter set assures the exact at . Away from , the deviation 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.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (a) and (b) as functions of . The parameters are U=2.0, T=0.1, .

Figure 3 shows the same averages as functions of U for T=0.1, , and a negative bath bias δω= −0.2. At U=0.0, all projection truncations give exact results for and . 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 as soon as . pAAA gives good agreement up to U=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.

	Figure Option ViewDownloadNew Window
	Fig. 3. (a) and (b) as functions of U. The parameters are T=0.1, δ ω = −0.2, and .

The temperature dependence of the same quantities are shown in Fig. 4 for U=2.0, , δω= −0.2. In Fig. 4(a), using the particle–hole symmetry properties , we can deduce that the impurity spin polarization is zero at and it increases as temperature is lowered for all approximations. While pHIA and pAAA give which has the wrong sign of spin polarization, pLacroix gives correct sign and quantitative agreement with NRG for all temperatures. In Fig. 4(b), 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.

	Figure Option ViewDownloadNew Window
	Fig. 4. (a) and (b) as functions of temperature. The parameters are U=2.0, δ ω = −0.2, .

Figure 5 shows the LDOS obtained from different bases at a particle–hole asymmetric point , 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 . For pHIA, we propose the following quantity to measure the truncation error,

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Impurity density of states calculated at U=2.0, T=0.001, δ ω =0.0, .

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (a) relative truncation error eT in pHIA, (b) symmetrized double occupancy, as functions of . Parameters are U=2.0, T=0.1, and δ ω =0.0.

It is seen from Fig. 6 that e_T qualitatively reflects the errors in the physical quantities. At the particle–hole symmetric point , e_T=1.0, being consistent with the fact that at this point, pHIA recovers the conventional HIA which amounts to . In the regime , e_T stays 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, e_T quickly decreases and in Fig. 6(b), the double occupancies from pHIA and NRG merge. In the limits , or 1, we expect e_T=0. pHIA will become exact since no electron correlation is present in those limits. Figure 6 shows a qualitative correlation between e_T and 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 e_T and the errors. This is because in PTA, the physical quantities depend on the truncation error non-linearly. Therefore, we conclude that e_T can be used to gauge the overall level of approximation of PTA. It is especially suitable for comparing the error among different parameter regimes.

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 space , look for the operator in a subspace and require that is as close to A as possible. The solution will be the projection of A into . 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 is not calculate exactly (except for those trivial ones such as ) 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, , , 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 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.

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].