Quantum many-body systems in thermal equilibrium can be described by the imaginary time Green’s function formalism. However, the treatment of large molecular or solid ab initio problems with a fully realistic Hamiltonian in large basis sets is hampered by the storage of the Green’s function and the precision of the solution of the Dyson equation. We present a Legendre-spectral algorithm for solving the Dyson equation that addresses both of these issues. By formulating the algorithm in Legendre coefficient space, our method inherits the known faster-than-exponential convergence of the Green’s function’s Legendre series expansion. In this basis, the fast recursive method for Legendre polynomial convolution enables us to develop a Dyson equation solver with quadratic scaling. We present benchmarks of the algorithm by computing the dissociation energy of the helium dimer He2 within dressed second-order perturbation theory. For this system, the application of the Legendre spectral algorithm allows us to achieve an energy accuracy of 10−9Eh with only a few hundred expansion coefficients.

The equilibrium properties of many-body quantum systems can be described by the finite temperature imaginary-time Green’s function formalism,1 which is widely applicable to condensed matter physics, quantum chemistry, and materials science. Applications include numerical methods for low energy effective model Hamiltonians such as lattice Monte Carlo,2 dynamical mean field theory3 and its extensions,4–6 and diagrammatic Monte Carlo.7Ab initio calculations using the random phase approximation,8 self-consistent second order perturbation theory,9–17 Hedin’s GW approach,18–26 and self-energy embedding theory27–33 can also be formulated in imaginary time.

While the finite temperature Green’s function formalism is very successful in applications to model Hamiltonians, its applicability to quantum chemistry and materials science remains limited to simple molecular and periodic problems. This is due to the necessity of simultaneously describing both the core and valence orbitals, which results in an energy scale that is difficult to describe by a single imaginary time/frequency grid. A simple equidistant Matsubara grid would contain millions of points, thus making the storage and manipulation of the Green’s functions computationally costly. In contrast, a grid with only a small number of equidistant points will result in a poorly converged energy or density matrix, making calculations with μHartree precision challenging. Such precision is necessary in applications where the evaluation of interaction energies,34–37 energies of conformers,38 or energies of high-lying excited states39,40 is needed. Consequently, it is important to develop a compact representation that yields highly converged properties.

With the standard approach using equidistant Matsubara frequency41 grids with finite frequency cutoff, the imaginary time Green’s function only converges to the analytical result linearly in the number of Matsubara frequencies. Amending the representation with a low order high frequency expansion results in polynomial convergence.42–44 In practice, this is problematic, since for systems with a wide range of energy scales, the number of coefficients is controlled by the largest energy scale.14 Alternatives such as uniform power meshes have had some success.45,46 However, the most compact representations are achieved using a set of (orthogonal) continuous basis functions directly in imaginary time, such as orthogonal polynomials47,48 or numerical basis functions.49–53 The convergence of such a representation is faster than exponential,47,48 and asymptotically superior to any polynomially converging representation.

In all imaginary time methods, a central step besides the solution of the impurity problem is the solution of the Dyson equation for the single particle Green’s function.54–57 In the Matsubara frequency representation,41 the Dyson equation is diagonal and can be readily solved. However, the solution is plagued by the polynomial convergence with respect to the number of frequency coefficients used. In imaginary time, the Dyson equation is a non-trivial integro-differential equation with a mixed boundary condition. Recently, an algorithm for solving the Dyson equation in imaginary time using the Chebyshev polynomials has been presented.48 This algorithm preserves the exponential convergence of the orthogonal polynomial expansion.47 However, the central convolution step has a cubic scaling in the expansion order NL, O(NL3), which limits the applicability of the algorithm.

The development of compact representations and algorithms for solving the Dyson equation is an active field of research (see Table I for an overview of the state-of-the-art methods). For a recent development, see Ref. 53.

TABLE I.

Overview of Green’s function representation approaches in both Matsubara frequency space and imaginary time combined with the scaling of solvers for the Dyson equation. Where no convergence is listed, the scaling either involves additional parameters or is unknown.

DomainBasisConvergenceCompactnessDyson scaling
Matsubara frequency Finite frequency cutoff O(1) Poor O(N) 
 Tail correction, pth order42–44  O(Np) Fair … 
 Spline grid14  … Good … 
Both frequency and time Sparse sampling52  … See Ref. 52a O(N) 
 Minimax isometry53  … See Ref. 53  O(N) 
Imaginary time Uniform mesh O(N1) Poor O(N3) 
 Power mesh45,46,58 … Fair O(N3) 
Orthogonal functions Intermediate representation49–51  O(eN) Excellent No 
 Chebyshev polynomials48  O(eN) Very good O(N3) 
 Legendre polynomials (this workO(eN)47  Very good O(N2) 
DomainBasisConvergenceCompactnessDyson scaling
Matsubara frequency Finite frequency cutoff O(1) Poor O(N) 
 Tail correction, pth order42–44  O(Np) Fair … 
 Spline grid14  … Good … 
Both frequency and time Sparse sampling52  … See Ref. 52a O(N) 
 Minimax isometry53  … See Ref. 53  O(N) 
Imaginary time Uniform mesh O(N1) Poor O(N3) 
 Power mesh45,46,58 … Fair O(N3) 
Orthogonal functions Intermediate representation49–51  O(eN) Excellent No 
 Chebyshev polynomials48  O(eN) Very good O(N3) 
 Legendre polynomials (this workO(eN)47  Very good O(N2) 
a

The compactness of the sparse sampling approach depends on the real-time basis employed.

In this paper, we present a Legendre spectral method for solving the Dyson equation with super-exponential convergence and a convolution that scales quadratically O(NL2), one order better than previous formulations.48 The super-exponential convergence allows us to achieve an energy accuracy of 10−9Eh in a realistic quantum chemistry system with a few hundred expansion coefficients. We show this in a proof-of concept benchmark: computing the dissociation energy of He2 using self-consistent second-order perturbation theory, taking both the zero temperature and the complete basis limit.

This paper is organized as follows: In Sec. II, we introduce the Dyson equation. In Sec. III, we present our Legendre spectral method. In Secs. IV and V, we apply our method to a realistic quantum chemistry problem, the dissociation energy of the noble gas He2. In Sec. VI, we present conclusions.

The imaginary time single particle Green’s function G is defined on the interval τ ∈ [−β, β], GG(τ), where β is the inverse temperature β = 1/T. It obeys the periodicity condition G(−τ) = ξG(βτ), with ξ = + 1 (−1) for bosons (fermions), making it an (anti-)periodic function with a step discontinuity at τ = 0 [see Fig. 1(a)]. The imaginary time Dyson equation for G(τ) is54–57 

[τh]G(τ)Σ*G=0,
(1)

where h is the single particle energy and Σ is the self-energy, which accounts for all many-body interactions. We note in passing that Σ(τ) has the same periodicity as G(τ). The boundary condition for Eq. (1) is G(0) − ξG(β) = −1, and the Fredholm type59 imaginary time convolution is defined as Σ*G0βdτ¯Σ(ττ¯)G(τ¯).

FIG. 1.

Single particle Green’s function in (a) imaginary time G(τ), (b) Matsubara frequency G(n) [with (iωn)1 black line], and (c) Legendre expansion coefficients Gn for site one in the fermionic two level system with the second quatization Hamiltonian H=μc1c1+V(c1c2+c2c1)+ϵc2c2 at inverse temperature β = 1, where ci creates and ci annihilates a fermion at site i and μ = −3, ϵ = 3.3, and V = 4.

FIG. 1.

Single particle Green’s function in (a) imaginary time G(τ), (b) Matsubara frequency G(n) [with (iωn)1 black line], and (c) Legendre expansion coefficients Gn for site one in the fermionic two level system with the second quatization Hamiltonian H=μc1c1+V(c1c2+c2c1)+ϵc2c2 at inverse temperature β = 1, where ci creates and ci annihilates a fermion at site i and μ = −3, ϵ = 3.3, and V = 4.

Close modal

Analytically, the Dyson equation [Eq. (1)] can be solved using the Fourier series expansion,

G(τ)=1βn=eiωnτG(iωn),G(iωn)=0βdτeiωnτG(τ),

where the Matsubara frequencies n are given by iωniπβ(2n+η) with η = (1 − ξ)/2 and n integers.54–57 In Matsubara frequency space, the Dyson equation [Eq. (1)] is diagonal,41 

[iωnhΣ(iωn)]G(iωn)=1.
(2)

Numerically, however, the discontinuity at τ = 0 results in a slow asymptotic decay G(iωn)(iωn)1 as n → ±i∞ [see Fig. 1(b)]. This prevents a naïve finite frequency |n| < Nω approximation G(τ)1β|n|<NωeiωnτG(iωn) from converging in Nω [the maximal error in G(τ) scales as O(Nω0)=O(1)]. The standard solution to this problem is to use a finite number p of high-frequency “tail” coefficients G¯k to approximate G(iωn)k=1pG¯k/(iωn)k for |n| > Nω, where the known asymptotic decay implies G¯1 = 1. This type of tail correction procedure gives polynomial convergence in G(τ) with the power determined by the order p of the tail expansion O(Nωp) [see, e.g., Refs. 42–44]. In Fig. 2, this is shown for the case of p = 3 using the TRIQS library.60 

FIG. 2.

Error in density Δn as a function of Legendre expansion order NL and number of Matsubara frequencies Nω for the same system as in Fig. 1.

FIG. 2.

Error in density Δn as a function of Legendre expansion order NL and number of Matsubara frequencies Nω for the same system as in Fig. 1.

Close modal

Since G(τ) is continuous on τ ∈ [0, β], it can be much more efficiently represented by a finite orthogonal polynomial expansion,

G(τ)n=0NLGnLn[x(τ)],
(3)

where Ln[x] are Legendre polynomials defined on x ∈ [−1, 1] and x(τ)=2τβ1. The Legendre coefficients Gn have a faster than exponential asymptotic decay47 [see Fig. 1(c)]. This also causes the finite NL expansion at the right-hand side of Eq. (3) to converge faster than exponential O(eNL) to the analytical G(τ).

Here, we develop a Legendre spectral method for solving the Dyson equation [Eq. (1)], reformulating the integro-differential equation in the space of Legendre coefficients Gn [Eq. (3)]. In the space of a finite Legendre expansion of order NL, Eq. (1) is cast to a linear equation system,

n=0NLDknh1knΣ*knGn=0k,
(4)

where terms with one and two indices are vectors and matrices in Legendre coefficient space. The last row of the left-hand side of the matrix is modified to enforce the boundary condition of Eq. (1). The resulting method has faster than exponential convergence and quadratic scaling O(NL2), one order better than previous approaches.48 

The differential operator τ in Eq. (1) acting on the Legendre polynomials takes the form61 

τLn[x(τ)]=2βxLn(x)=2βk=0,k+noddn1(2k+1)Lk(x)=kDknLk(x).
(5)

Hence, the derivative matrix Dkn in Eq. (4) is given by

β2Dkn2k+1,0kn,k+nodd,0,elsewhere
(6)

and is upper triangular (see Fig. 3). Using Ln(±1) = (±1)n, the Dyson equation boundary condition can be written as

1=G(0)ξG(β)=n((1)nξ)Gn.
(7)
FIG. 3.

Matrix structure of the spectral derivative operator Dkn and the convolution operator [Σ*]kn for ϵ = 1 and Σ(τ)=eϵτ(ξeϵβ1)1 at β = 1 and ξ = −1 (fermions).

FIG. 3.

Matrix structure of the spectral derivative operator Dkn and the convolution operator [Σ*]kn for ϵ = 1 and Σ(τ)=eϵτ(ξeϵβ1)1 at β = 1 and ξ = −1 (fermions).

Close modal

The imaginary time convolution [Σ * G] in the Dyson equation [Eq. (1)] can be separated into the two terms of Volterra type,

[Σ*G](τ)=0βdτΣ(ττ)G(τ)=0τdτΣ(ττ)G(τ)+τβdτξΣ(β+ττ)G(τ),
(8)

using the periodicity property Σ(−τ) = ξΣ(βτ). In Eq. (8), Σ(τ) is only evaluated for τ ∈ [0, β], avoiding the discontinuity at τ = 0.

In Legendre coefficient space, the convolution operator [Σ *] can be written as a sum of two matrices Bkn, representing the two Volterra terms [Eq. (8)],

[Σ*]knBkn<+ξBkn>.
(9)

Stable recursion relations for Bnk have been derived by Hale and Townsend62 using the Fourier connection between Legendre polynomials and spherical Bessel functions. Since the derivation is detailed in Ref. 62, we only state the result specialized to the imaginary time convolution in Eq. (8) here and provide a derivation in the  Appendix.

The coefficients are related by the recursion relation,

Bk,n+1=2n+12k+3Bk+1,n+2n+12k1Bk1,n+Bk,n1,
(10)

which for each column require two previous columns to be known. The recursion is only stable for the lower triangular coefficients in Bkn. The upper triangular coefficients are computed using the transpose relation,

Bk,n=(1)n+k2k+12n+1Bn,k.
(11)

The two first columns are given by the starting relations,

Bk,0=Σ0±Σ13,k=0,±(Σk12k1Σk+12k+3),k1,Bk,1=Bk,0+Bk1,02k1Bk+1,02k+3,k1,
(12)

with the special case for k = 0, B0,1=B1,0/3, using the Legendre coefficients Σn of the self-energy Σ [cf. Eq. (3)].

Since each coefficient in Bkn can be computed in O(1) operations, the scaling of the convolution matrix construction is O(NL2). The self-energy Σ(τ) is a smooth function with asymptotic exponentially decaying Legendre coefficients, which causes the entries of the dominantly diagonal spectral convolution operator [Σ *]kn to decay exponentially both along and away from the diagonal (see Fig. 3).

The numerical solution of G(τ) from the Dyson equation constructed in terms of the linear system in Eq. (4) converges faster than exponentially to the analytical solution, with the increased number of Legendre coefficients NL (see Fig. 2). This is in stark contrast to the polynomial convergence of the standard Matsubara tail approach42–44 (also shown in Fig. 2).

To retain the high accuracy of the Legendre spectral Dyson solver, the method has to be complemented with stable transforms between Legendre coefficients and imaginary time,

Gn=i=0NLSniG(τi),G(τi)=n=0NLLinGn.
(13)

To construct the well-conditioned transform matrices Sni and Lin, we employ the Legendre quadrature and the Legendre–Gauss–Lobatto points xi{x:(1x2)LNL(x)=0}, x0 = −1, xN = 1, re-scaled to the imaginary time interval [0, β], τi=βxi+12. Using xi, the matrices Sni and Lin can be directly constructed (avoiding matrix inversion),

Lin=Ln(x(τi)),Sni=β2WnωiLn(x(τi)),
(14)

where 11dxLn(x)Lm(x)=δnm22n+1δnmWn and ωi=2N(N+1)1LNL(xi)2 (see, e.g., Refs. 61 and 63).

As a proof of concept application of the Legendre spectral Dyson solver developed in this paper, we employ the solver in a quantum chemistry setting using a Gaussian basis set. We will employ self-consistent second order perturbation theory, also known as GF2,9–17 which has seen a revival in recent years, both in ab initio condensed matter applications8,15 and in quantum chemistry11,16,58,64 in combination with embedding methods.29 Our implementation is built on the Coulomb integrals of the pyscf library.65 

In the resulting non-orthogonal basis set, the Dyson equation takes the form

j[Sij(τμ)+Fij+Σij*]Gjk(τ)=0,
(15)

where i, j, k are orbital indices, Sij is the overlap matrix, and Fij is the so-called Fock matrix, Fijhij+Σij(HF). The boundary condition for this equation is j(Gij(0) − ξGij(β)) · Sjk = −1ik. Here, the single particle term hij accounts for electronic kinetic and nuclear-electronic matrix elements, and the Hartree–Fock self-energy Σij(HF) is given by

Σij(HF)=klPkl(vijklvilkj/2),
(16)

where Pij is the density matrix Pij = −2Gij(β) and vijkl is the electron–electron Coulomb repulsion integral.

In GF2, the imaginary-time-dependent part of the self-energy Σ(τ) is approximated with the second order self-energy diagram using the full electron Greens function G, Σ ≈ Σ(GF2)[G], where

Σij(GF2)(τ)=klmnpqGkl(τ)Gmn(τ)Gpq(βτ)×vimpk(2vjnlqvjlnq).
(17)

The evaluation of Σ(GF2)(τ) for fixed τ scales as O(N5),64 where N is the number of atomic orbitals.

Solving for the GF2 Green’s function, G amounts to solving Eqs. (15)–(17), which is a highly non-linear problem. To find the solution, we perform self-consistent iterations (see Fig. 4 for a schematic picture). The inner loop solves the Dyson equation [Eq. (15)] and updates the Hartree–Fock self-energy Σ(HF) [Eq. (16)] until convergence (in the Fock-matrix F). At convergence in F, one step of the outer loop is performed by re-evaluating the GF2 self-energy Σ(GF2) [Eq. (17)] and computing the relative change in total energy E. If the change is above a fixed threshold, the inner loop is started again. To compute the inter molecular energies, which is an energy difference, we need a threshold of 10−10.

FIG. 4.

Schematic of the GF2-self-consistency loop.

FIG. 4.

Schematic of the GF2-self-consistency loop.

Close modal

The total energy E of the system is given by

E=12Tr[(h+F)P]+Tr[Σ*G]+Enn,
(18)

where Enn is the nuclei–nuclei Coulomb energy. The imaginary time trace Tr[·] is defined as Tr[A] ≡ −iAii(β),44,66 and the Σ * G convolution is computed with the spectral Legendre convolution as in Eq. (4).

The faster than exponential convergence of the Legendre spectral Dyson solver [Eq. (4)] is particularly suited for high precision calculations. A prime example is the computation of the binding energy De in noble-gas dimers, where the weak bonding requires high precision calculations of total energies. The binding energy De is obtained from the minimum of the interaction energy Eint(r) as a function of atomic separation r,

DeEint(re)minrEint(r),
(19)

where re is the equilibrium atomic distance. The interaction energy Eint is, in turn, given by

Eint(r)EA2(r)2EA(r),
(20)

where EA2 is the total energy of the dimer and EA is the total energy of the single atom (the monomer) evaluated using the standard counterpoise correction.67 In the noble gases, the total energies EA and EA2 are of the order of Hartrees (∼Eh ≡ 1 Hartree), while the binding energy De is of the order of tens of micro-Hartrees (∼10μEh), hence requiring high precision calculation of the total energies.

We use He2 as a prototype system since there exist published reference results for the binding energy De and equilibrium distance re calculated with Hartree–Fock (HF), second-order Moller–Plesset perturbation theory (MP2), coupled cluster singles doubles (CCSD) theory, and coupled cluster singles doubles and non-iterative perturbative triples [CCSD(T)] theory.68 The MP2 method is closely related to GF2 and uses the second order self-energy [Eq. (17)] evaluated at the HF Green’s function G(HF), Σ(MP2) ≡ Σ(GF2)[G(HF)]. However, note that the prefactors in the total energy differ.10,69

Figure 5 shows Eint(r) (and −De) of He2 computed with HF, MP2, and GF2 in the aug-cc-pvqz basis together with CCSD and CCSD(T) reference results on De.68 The GF2 results are obtained by fitting a 4th order polynomial to 21 r-points of Eint(r) computed in a 0.1 Bohr range centered around the minimum at re. The GF2 results are obtained using the Legendre spectral Dyson solver, while HF and MP2 are computed using pyscf.65 As shown in Fig. 5, He2 does not bind within the Hartree–Fock approximation, which gives a strictly positive interaction energy. Compared to MP2 our GF2 results are a considerable improvement using the coupled cluster CCSD and CCSD(T) as a reference.

FIG. 5.

Interaction energy Eint as a function of atomic distance r of He2 with basis aug-ccpvqz using HF, MP2, and GF2. The HF and MP2 results are computed with pyscf;65 the GF2 results are computed using β=50Eh1, Nτ = 192. The CCSD and CCSD(T) results are from Ref. 68.

FIG. 5.

Interaction energy Eint as a function of atomic distance r of He2 with basis aug-ccpvqz using HF, MP2, and GF2. The HF and MP2 results are computed with pyscf;65 the GF2 results are computed using β=50Eh1, Nτ = 192. The CCSD and CCSD(T) results are from Ref. 68.

Close modal

In order to extrapolate the results to the complete basis set (CBS) limit,70,71 we repeat the calculations using the augmented correlation consistent (aug-cc-pvnz) basis set series, with n = d, t, q, 5 (i.e., n = 2, 3, 4, 5).72–74 This series has been shown to enable accurate extrapolation of a number of properties due to its systematic convergence in n.70,75–85

In Table II, we summarize the binding energy De and equilibrium distance re of He2 computed by MP2, CCSD, CCSD(T), and GF2 using the aug-cc-pv{d,t,q,5}z basis sets. The aug-cc-pv{d,t,q,5}z GF2 energies are computed at β=50Eh1 using Nτ = 128, 160, 192, and 250 τ-points, respectively. The convergence in Nτ is imposed so that the absolute values of the elements in highest Legendre coefficient matrix are smaller than 10−10. The zero temperature convergence (at β=50Eh1) is ensured by requiring that the finite temperature MP2 total energy differs by less than 0.1 nHartree compared to the zero temperature MP2 total energy from pyscf.

TABLE II.

Dissociation energies De (top) and equilibrium distances re computed by MP2, CCSD, CCSD(T), and GF2 with the basis sets aug-cc-pvnz, with n = d, t, q, 5. The MP2, CCSD, and CCSD(T) results are from Ref. 68.

De (μEh)MP2CCSDCCSD(T)GF2
aug-ccpvdz 12.69 16.78 18.57 18.17 
aug-ccpvtz 17.97 23.77 27.10 24.63 
aug-ccpvqz 19.66 25.79 29.64 26.59 
aug-ccpv5z 20.71 27.09 31.25 27.79 
CBS 22.98 30.06 34.70 29.67 
 
De (μEh)MP2CCSDCCSD(T)GF2
aug-ccpvdz 12.69 16.78 18.57 18.17 
aug-ccpvtz 17.97 23.77 27.10 24.63 
aug-ccpvqz 19.66 25.79 29.64 26.59 
aug-ccpv5z 20.71 27.09 31.25 27.79 
CBS 22.98 30.06 34.70 29.67 
 
re (Bohr)MP2CCSDCCSD(T)GF2
aug-ccpvdz 6.1680 6.0580 6.0086 6.0547 
aug-ccpvtz 5.9175 5.8060 5.7452 5.8244 
aug-ccpvqz 5.8606 5.7546 5.6891 5.7722 
aug-ccpv5z 5.8244 5.7210 5.6537 5.7388 
CBS 5.769 5.672 5.607 5.680 
re (Bohr)MP2CCSDCCSD(T)GF2
aug-ccpvdz 6.1680 6.0580 6.0086 6.0547 
aug-ccpvtz 5.9175 5.8060 5.7452 5.8244 
aug-ccpvqz 5.8606 5.7546 5.6891 5.7722 
aug-ccpv5z 5.8244 5.7210 5.6537 5.7388 
CBS 5.769 5.672 5.607 5.680 

We note that the number of τ-points Nτ used for the aug-cc-pv{d,t,q,5}z basis sets is of the same order as the number of atomic orbitals N. Hence, the scaling of GF2, O(NτN5), is comparable to the scaling of CCSD, O(N6). As given in Table II, the accuracy of the GF2 result for De is comparable to that of CCSD when compared to that of CCSD(T), while the CCSD result for re is closer to the CCSD(T) result than that of GF2. This makes GF2 a considerable improvement over MP2.

With the systematic convergence of De and re as a function of basis set size n, it is possible to extrapolate to the complete basis limit n.68 We extrapolate De and re using our GF2 aug-ccpv{t,q,5}z results by fitting the exponential model: A · eB(n−2) + C, proposed in Ref. 68, where A, B, and C are parameters. The applicability of the model is checked by a logarithmic plot (see Fig. 6). The resulting CBS limit of our GF2 results is De ≈ 29.67μEh and re ≈ 5.680 a0 (see also Table II).

FIG. 6.

Basis extrapolation of equilibrium distance (re) and dissociation energy (De) He2 with basis aug-ccpvnz with n = 2, 3, 4, 5. Left panels: raw data and fitting. Right panels: check of fitting results.

FIG. 6.

Basis extrapolation of equilibrium distance (re) and dissociation energy (De) He2 with basis aug-ccpvnz with n = 2, 3, 4, 5. Left panels: raw data and fitting. Right panels: check of fitting results.

Close modal

We introduce a Legendre-spectral algorithm for solving the Dyson equation in Legendre coefficient space. By staying in Legendre-coefficient space, the algorithm converges super exponentially with respect to the number of Legendre coefficients NL used to represent the imaginary time Green’s function.47 This is in stark contrast to the Matsubara frequency space based approach with only polynomial convergence.42–44 The exponential convergence is shared with a recently presented Chebyshev polynomial based algorithm,48 where the convolution scales as O(NL3). Currently, there is no known algorithm for Chebyshev series that can compute the convolution term with the same efficiency as in the Legendre series.62 Our work goes beyond this, employing a Legendre convolution with O(NL2) scaling and enabling the application to larger ab initio systems.

To benchmark the algorithm, we apply it to the quantum chemistry computation of the dissociation energy of the noble gas He2 using self-consistent second order perturbation theory (GF2). The exponential convergence of our algorithm allows us to reach the required 10−9Eh zero temperature total-energy precision using only 100–200 Legendre coefficients in the Dunning basis series aug-ccpvnz.72–74 

The algorithm is also relevant for condensed matter ab initio applications in periodic systems that require high precision, such as GF29–17 and Hedin’s GW.18–26 This is a promising venue for future research.

H.U.R.S. would like to acknowledge helpful discussions and support from (i) Alex Barnett and Manas Rachh on the fundamentals of spectral methods, (ii) Keaton Burns for pointing out Ref. 62, (iii) Sergei Iskakov for providing independent reference results for testing the pyscf-GF2 implementation and useful discussions, (iv) Lewin Boehnke, Andreas Herrmann, Philipp Werner, and Hiroshi Shinaoka, who took part in early discussion on Green’s function representations, which guided the development of the method, and (v) Antoine Georges, Olivier Parcollet, Manuel Zingl, Alexandru Georgescu, Igor Krivenko, and Nils Wentzell, who contributed through discussions and contributions to the TRIQS project. E.G. and X.D. were supported by the Simons Collaboration on the many-electron problem. D.Z. acknowledges support from Grant No. NSFCHE-1453894. The Flatiron Institute is a division of the Simons Foundation.

In this appendix, we derive Eqs. (10)–(12) in the main text. The derivation follows Ref. 62 but with more details for both integrals in Eq. (8).

1. Convolution and Fourier transform

The convolution of two continuous integrable functions is defined as62 

h(x)=(f*g)(x)dtf(t)g(xt).
(A1)

With the assumption f and g being periodic functions, their Fourier transform can be written as

F{f}(ω)=dxeiωxf(x),
(A2)
F1{f}(x)=12πdxeiωxf(x),
(A3)

which satisfy the Fourier inversion theorem F1{F{f}}=f and convolution theorem,86 

F{f*g}=F{f}F{g}.
(A4)

2. Legendre polynomials

The Legendre polynomials Pn(x) can be defined recursively using the three term recurrence relation,

P0(x)=1,P1(x)=x,(n+1)Pn+1(x)=(2n+1)xPn(x)nPn1(x).
(A5)

They are orthogonal on [−1, 1],

11dxPm(x)Pn(x)=δm,n22n+1,
(A6)

and the derivatives satisfy the recurrence relation,

(2n+1)Pn(x)=ddxPn+1(x)Pn1(x).
(A7)

The Fourier transform and inverse Fourier transform of the Legendre polynomials can be expressed in terms of Bessel functions of the first kind,

F{Pn}=11dxeiωxPn(x)=2(i)njn(ω),
(A8)
F1{Pn}=11dxeiωxPn(x)=2injn(ω),
(A9)

where jn(z) is the nth spherical Bessel function and Pn = 0 outside [−1, 1].

By combining Eqs. (A4) and (A8), the convolution of Legendre polynomials can be expressed in terms of Bessel functions,

(Pm*Pn)(x)=2(i)m+nπdωeiωxjm(ω)jn(ω).
(A10)

This is the central observation of Ref. 62 that enables the derivation of recursion relations for the Legendre polynomial convolution.

The main property of spherical Bessel functions used is the three term recurrence relation,

j1(z)=coszz,j0(z)=sinzz,jn+1(z)=2n+1zjn(z)jn1(z),n0.
(A11)

The convolution equation [Eq. (A1)] can be computed by replacing the two continuous function f(x) and g(x) on the bounded interval with polynomial approximates fM(x) and gN(x) of sufficiently high degree. With two Legendre series fM(x) and gN(x) supported on x ∈ [−1, 1],

fM(x)=m=0MαmPm(x),gN(x)=n=0NβnPn(x).
(A12)

Equation (A1) becomes

h(x)=(fM*gN)(x)=max(1,x1)min(1,x+1)dtfM(t)gN(xt)=1x+1dtfM(t)gN(xt)+x11dtfM(t)gN(xt),
(A13)

which can be computed separately in two integration domains x ∈ [−2, 0] and x ∈ [0, 2] (see Fig. 4.1 in Ref. 62).

a. First interval x [−2, 0]

For x ∈ [−2, 0], we have h(x) = h<(x), where

h<(x)=1x+1dtfM(t)gN(xt)=k=0M+N+1γk<Pk(x+1).
(A14)

Using the orthogonality of Legendre polynomials [Eq. (A6)], we have

γk<=2k+1220dxPk(x+1)1x+1dtfM(t)gN(xt)=n=0Nβn2k+12m=0Mαm20dxPk(x+1)×1x+1dtPm(t)Pn(xt),
(A15)

collecting terms, we can write γk<=n=0NBk.n<βn, where

Bk,n<=2k+12m=0Mαm20dxPk(x+1)×1x+1dtPm(t)Pn(xt)=2k+12m=0Mαm20dxPk(x+1)(Pm*Pn)(x)=2k+12m=0Mαm11dsPk(s)(Pm*Pn)(s1).
(A16)

Using the Fourier expression for the Legendre convolution [Eq. (A10)], Bk,n< can be expressed in terms of spherical Bessel functions,

Bk,n<=2k+1πm=0M(i)m+nαm11dsPk(s)×dωeiω(s1)jm(ω)jn(ω).
(A17)

Consider the Bk,n+1< term, changing the order of integration and Fourier transforming the remaining Legendre polynomial give

Bk,n+1<=2(2k+1)πm=0M(i)m+n+1ikαm×dωjk(ω)jm(ω)jn+1(ω)eiω.
(A18)

Applying the recursion relation of the spherical Bessel functions [Eq. (A11)] on n and k, we have

(i)m+n+1ikjk(ω)jm(ω)jn+1(ω)  =(i)m+n+1ikjk(ω)jm(ω)2n+1ωjn(ω)jn1(ω)  =2n+12k+1(i)m+n+1ikjk+1(ω)+jk1(ω)jm(ω)jn(ω)   +(i)m+n1ikjk(ω)jm(ω)jn1(ω).
(A19)

Back insertion in Eq. (A18) and simplifying prefactors in k give

Bk,n+1<=2n+12k+3Bk+1,n<+2n+12k1Bk1,n<+Bk,n1<.
(A20)
b. Second interval x [0, 2]

For x ∈ [0, 2], we have h(x) = h>(x), where

h>(x)=x11dtfM(t)gN(xt)=k=0M+N+1γk>Pk(x1).
(A21)

γk> can be computed in the same way as γk< [see Eq. (A15)],

γk>=2k+1202dxPk(x1)x11dtfM(t)gN(xt)=n=0Nβn2k+12m=0Mαm02dxPk(x1)×x11dtPm(t)Pn(xt),
(A22)

collecting terms, we can write γk>=n=0NBk.n>βn, where

Bk,n>=2k+12m=0Mαm02dxPk(x1)×x11dtPm(t)Pn(xt)=2k+12m=0Mαm02dxPk(x1)(Pm*Pn)(x)=2k+12m=0Mαm11dsPk(s)(Pm*Pn)(s+1).
(A23)

Using the Fourier expression for the Legendre convolution [Eq. (A10)] gives

Bk,n>=2k+1πm=0M(i)m+nαm11dsPk(s)×dωeiω(s+1)jm(ω)jn(ω).
(A24)

Since the exponent in the integral is unchanged when applying the recursion relations of the spherical Bessel functions, we conclude that B> obeys the same recursion relation as B<, albeit with a different starting point since the “seeding” integrals have a different sign in the exponent.

c. Summary

The convolution matrices for both intervals can be expressed as the integral sums,

Bk,n=2(2k+1)πm=0M(i)m+nikαm×dωjk(ω)jm(ω)jn(ω)eiω,
(A25)

differing only in the sign in the exponent. The coefficients are related by the recursion relation,

Bk,n+1=2n+12k+3Bk+1,n+2n+12k1Bk1,n+Bk,n1.
(A26)

In practice this recursion relation is only stable below the diagonal with k > n. To get entries above diagonal, the transpose relation that can be derived from the integral expression [Eq. (A18)], is used,

Bk,n=(1)n+k2k+12n+1Bn,k.
(A27)

3. Initial values Bk,0 and Bk,1

To start the recursion, the initial values for n = 0 and 1 are needed. To derive explicit expressions for these terms, we repeatedly use the Volterra integral formula for Legendre polynomials from Ref. 87,

Sa,n(x)=axdtPn(t),
(A28)
Sa,0(x)=xa,
(A29)
Sa,n(x)=12n+1Pn+1(t)Pn1(t)ax,
(A30)

for a = ±1, we get

S±1,0(x)=x1=P1(x)P0(x),
(A31)
S±1,n(x)=12n+1Pn+1(x)Pn1(x),
(A32)

where we have used Pn(±1) = (±1)n to cancel the constant terms.

Returning to the convolution matrices, we have for Bk,n< and n = 0, using P0(x) = 1,

Bk,0=±2k+12m=0Mαm11dxPk(x)1xdtPm(t)=±2k+12m=0Mαm11dxPk(x)S1,m(x)=±2k+12m=0Mαm2m+111dxPk(x)Pm+1(x)Pm1(x),
(A33)

repeatedly using the Legendre orthogonality relation [Eq. (A6)] gives

Bk,0=α0α13,k=0,±(αk12k1αk+12k+3),k1.
(A34)

For the second column with n = 1, we detail the derivation of Bk,1<, the other case Bk,1> can be done analogously. Using P1(x) = x, we get

Bk,1<=2k+12m=0Mαm20dxPk(x+1)1x+1dtPm(t)P1(xt)=2k+12m=0Mαm11dxPk(x)1xdtPm(t)(xt1)=Bk,0<+2k+12m=0Mαm11dxPk(x)1xdtPm(t)txds=Bk,0<+2k+12m=0Mαm11dxPk(x)r1xds1sdtPm(t),
(A35)

where the last step is obtained by changing the order of integration. The last integral relation is a double Volterra integral and can, hence, be written using S−1,m(x) as

Bk,1<=Bk,0<+2k+12m=0Mαm11dxPk(x)1xdsS1,m(s)=Bk,0<+12m=0Mαm11dx[Pk1(x)Pk+1(x)]S1,m(x),
(A36)

where we, in the second step, have used partial integration and the Legendre derivative relation [Eq. (A7)].

For the second case Bk,1>, the only difference is when we change the integration variable, we get (xt + 1), instead, of (xt − 1) in Eq. (A35), so the sign before Bk,0 is changed to +1. By using Eq. (A33), we obtain the recursion relation,

Bk,1=Bk,0+Bk1,02k1Bk+1,02k+3,k1,
(A37)

with the special case for k = 0, B0,1=B1,0/3.

1.
A. A.
Abrikosov
,
I.
Dzyaloshinskii
,
L. P.
Gorkov
, and
R. A.
Silverman
,
Methods of Quantum Field Theory in Statistical Physics
(
Dover
,
New York, NY
,
1975
).
2.
R.
Blankenbecler
,
D. J.
Scalapino
, and
R. L.
Sugar
, “
Monte Carlo calculations of coupled boson-fermion systems. I
,”
Phys. Rev. D
24
,
2278
2286
(
1981
).
3.
A.
Georges
,
G.
Kotliar
,
W.
Krauth
, and
M. J.
Rozenberg
, “
Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions
,”
Rev. Mod. Phys.
68
,
13
125
(
1996
).
4.
A.
Toschi
,
A. A.
Katanin
, and
K.
Held
, “
Dynamical vertex approximation: A step beyond dynamical mean-field theory
,”
Phys. Rev. B
75
,
045118
(
2007
).
5.
A. N.
Rubtsov
,
M. I.
Katsnelson
, and
A. I.
Lichtenstein
, “
Dual fermion approach to nonlocal correlations in the Hubbard model
,”
Phys. Rev. B
77
,
033101
(
2008
).
6.
T.
Maier
,
M.
Jarrell
,
T.
Pruschke
, and
M. H.
Hettler
, “
Quantum cluster theories
,”
Rev. Mod. Phys.
77
,
1027
1080
(
2005
).
7.
N.
Prokof’ev
and
B.
Svistunov
, “
Bold diagrammatic Monte Carlo technique: When the sign problem is welcome
,”
Phys. Rev. Lett.
99
,
250201
(
2007
).
8.
M.
Kaltak
,
J.
Klimeš
, and
G.
Kresse
, “
Low scaling algorithms for the random phase approximation: Imaginary time and Laplace transformations
,”
J. Chem. Theory Comput.
10
,
2498
2507
(
2014
).
9.
P.
García-González
and
R. W.
Godby
, “
Self-consistent calculation of total energies of the electron gas using many-body perturbation theory
,”
Phys. Rev. B
63
,
075112
(
2001
).
10.
N. E.
Dahlen
and
R.
van Leeuwen
, “
Self-consistent solution of the Dyson equation for atoms and molecules within a conserving approximation
,”
J. Chem. Phys.
122
,
164102
(
2005
).
11.
J. J.
Phillips
and
D.
Zgid
, “
Communication: The description of strong correlation within self-consistent Green’s function second-order perturbation theory
,”
J. Chem. Phys.
140
,
241101
(
2014
).
12.
J. J.
Phillips
,
A. A.
Kananenka
, and
D.
Zgid
, “
Fractional charge and spin errors in self-consistent Green’s function theory
,”
J. Chem. Phys.
142
,
194108
(
2015
).
13.
A. A.
Kananenka
,
J. J.
Phillips
, and
D.
Zgid
, “
Efficient temperature-dependent Green’s functions methods for realistic systems: Compact grids for orthogonal polynomial transforms
,”
J. Chem. Theory Comput.
12
,
564
571
(
2016
).
14.
A. A.
Kananenka
,
A. R.
Welden
,
T. N.
Lan
,
E.
Gull
, and
D.
Zgid
, “
Efficient temperature-dependent Green’s function methods for realistic systems: Using cubic spline interpolation to approximate Matsubara Green’s functions
,”
J. Chem. Theory Comput.
12
,
2250
2259
(
2016
).
15.
A. A.
Rusakov
and
D.
Zgid
, “
Self-consistent second-order Green’s function perturbation theory for periodic systems
,”
J. Chem. Phys.
144
,
054106
(
2016
).
16.
A. R.
Welden
,
A. A.
Rusakov
, and
D.
Zgid
, “
Exploring connections between statistical mechanics and Green’s functions for realistic systems: Temperature dependent electronic entropy and internal energy from a self-consistent second-order Green’s function
,”
J. Chem. Phys.
145
,
204106
(
2016
).
17.
S.
Iskakov
,
A. A.
Rusakov
,
D.
Zgid
, and
E.
Gull
, “
Effect of propagator renormalization on the band gap of insulating solids
,”
Phys. Rev. B
100
,
085112
(
2019
).
18.
L.
Hedin
, “
New method for calculating the one-particle Green’s function with application to the electron-gas problem
,”
Phys. Rev.
139
,
A796
A823
(
1965
).
19.
F.
Aryasetiawan
and
O.
Gunnarsson
, “
The GW method
,”
Rep. Prog. Phys.
61
,
237
312
(
1998
).
20.
A.
Stan
,
N. E.
Dahlen
, and
R.
van Leeuwen
, “
Levels of self-consistency in the GW approximation
,”
J. Chem. Phys.
130
,
114105
(
2009
).
21.
A.
Kutepov
,
S. Y.
Savrasov
, and
G.
Kotliar
, “
Ground-state properties of simple elements from GW calculations
,”
Phys. Rev. B
80
,
041103
(
2009
).
22.
M. J.
van Setten
,
F.
Caruso
,
S.
Sharifzadeh
,
X.
Ren
,
M.
Scheffler
,
F.
Liu
,
J.
Lischner
,
L.
Lin
,
J. R.
Deslippe
,
S. G.
Louie
,
C.
Yang
,
F.
Weigend
,
J. B.
Neaton
,
F.
Evers
, and
P.
Rinke
, “
GW100: Benchmarking G0W0 for molecular systems
,”
J. Chem. Theory Comput.
11
,
5665
5687
(
2015
).
23.
E.
Maggio
,
P.
Liu
,
M. J.
van Setten
, and
G.
Kresse
, “
GW100: A plane wave perspective for small molecules
,”
J. Chem. Theory Comput.
13
,
635
648
(
2017
).
24.
M.
Grumet
,
P.
Liu
,
M.
Kaltak
,
J.
Klimeš
, and
G.
Kresse
, “
Beyond the quasiparticle approximation: Fully self-consistent GW calculations
,”
Phys. Rev. B
98
,
155143
(
2018
).
25.
A. L.
Kutepov
, “
Electronic structure of Na, K, Si, and LiF from self-consistent solution of Hedin’s equations including vertex corrections
,”
Phys. Rev. B
94
,
155101
(
2016
).
26.
A. L.
Kutepov
, “
Self-consistent solution of Hedin’s equations: Semiconductors and insulators
,”
Phys. Rev. B
95
,
195120
(
2017
).
27.
A. A.
Kananenka
,
E.
Gull
, and
D.
Zgid
, “
Systematically improvable multiscale solver for correlated electron systems
,”
Phys. Rev. B
91
,
121111
(
2015
).
28.
T. N.
Lan
,
A. A.
Kananenka
, and
D.
Zgid
, “
Communication: Towards ab initio self-energy embedding theory in quantum chemistry
,”
J. Chem. Phys.
143
,
241102
(
2015
).
29.
D.
Zgid
and
E.
Gull
, “
Finite temperature quantum embedding theories for correlated systems
,”
New J. Phys.
19
,
023047
(
2017
).
30.
T. N.
Lan
and
D.
Zgid
, “
Generalized self-energy embedding theory
,”
J. Phys. Chem. Lett.
8
,
2200
2205
(
2017
).
31.
T. N.
Lan
,
A.
Shee
,
J.
Li
,
E.
Gull
, and
D.
Zgid
, “
Testing self-energy embedding theory in combination with GW
,”
Phys. Rev. B
96
,
155106
(
2017
).
32.
L. N.
Tran
,
S.
Iskakov
, and
D.
Zgid
, “
Spin-unrestricted self-energy embedding theory
,”
J. Phys. Chem. Lett.
9
,
4444
4450
(
2018
).
33.
A. A.
Rusakov
,
S.
Iskakov
,
L. N.
Tran
, and
D.
Zgid
, “
Self-energy embedding theory (SEET) for periodic systems
,”
J. Chem. Theory Comput.
15
,
229
240
(
2019
).
34.
D. E.
Taylor
,
J. G.
Ángyán
,
G.
Galli
,
C.
Zhang
,
F.
Gygi
,
K.
Hirao
,
J. W.
Song
,
K.
Rahul
,
O.
Anatole von Lilienfeld
,
R.
Podeszwa
,
I. W.
Bulik
,
T. M.
Henderson
,
G. E.
Scuseria
,
J.
Toulouse
,
R.
Peverati
,
D. G.
Truhlar
, and
K.
Szalewicz
, “
Blind test of density-functional-based methods on intermolecular interaction energies
,”
J. Chem. Phys.
145
,
124105
(
2016
).
35.
G.
Chalasinski
and
M. M.
Szczesniak
, “
Origins of structure and energetics of van der Waals clusters from ab initio calculations
,”
Chem. Rev.
94
,
1723
1765
(
1994
).
36.
G.
Chalasinski
and
M. M.
Szczesniak
, “
State of the art and challenges of the ab initio theory of intermolecular interactions
,”
Chem. Rev.
100
,
4227
4252
(
2000
).
37.
G.
Chalasinski
and
M.
Gutowski
, “
Weak interactions between small systems. models for studying the nature of intermolecular forces and challenging problems for ab initio calculations
,”
Chem. Rev.
88
,
943
962
(
1988
).
38.
R.
Podeszwa
,
R.
Bukowski
, and
K.
Szalewicz
, “
Potential energy surface for the benzene dimer and perturbational analysis of π−π interactions
,”
J. Phys. Chem. A
110
,
10345
10354
(
2006
).
39.
S.
Sharma
,
K.
Sivalingam
,
F.
Neese
, and
G. K.-L.
Chan
, “
Low-energy spectrum of iron-sulfur clusters directly from many-particle quantum mechanics
,”
Nat. Chem.
6
,
927
933
(
2014
).
40.
Z.
Li
,
S.
Guo
,
Q.
Sun
, and
G. K.-L.
Chan
, “
Electronic landscape of the P-cluster of nitrogenase as revealed through many-electron quantum wavefunction simulations
,”
Nat. Chem.
11
,
1026
1033
(
2019
).
41.
T.
Matsubara
, “
A new approach to quantum-statistical mechanics
,”
Prog. Theor. Phys.
14
,
351
378
(
1955
).
42.
N.
Blümer
, “
Mott-Hubbard metal-insulator transition and optical conductivity in high dimensions
,” Ph.D. thesis,
Universität Augsburg
,
2002
.
43.
A.-B.
Comanac
, “
Dynamical mean field theory of correlated electron systems: New algorithms and applications to local observables
,” Ph.D. thesis,
Columbia University
,
2007
.
44.
D.
Hügel
,
P.
Werner
,
L.
Pollet
, and
H. U. R.
Strand
, “
Bosonic self-energy functional theory
,”
Phys. Rev. B
94
,
195119
(
2016
).
45.
W.
Ku
and
A. G.
Eguiluz
, “
Band-gap problem in semiconductors revisited: Effects of core states and many-body self-consistency
,”
Phys. Rev. Lett.
89
,
126401
(
2002
).
46.
W.
Ku
, “
Electronic excitations in metals and semiconductors: Ab initio studies of realistic many-particle systems
,” Ph.D. thesis,
University of Tennessee
,
2000
.
47.
L.
Boehnke
,
H.
Hafermann
,
M.
Ferrero
,
F.
Lechermann
, and
O.
Parcollet
, “
Orthogonal polynomial representation of imaginary-time Green’s functions
,”
Phys. Rev. B
84
,
075145
(
2011
).
48.
E.
Gull
,
S.
Iskakov
,
I.
Krivenko
,
A. A.
Rusakov
, and
D.
Zgid
, “
Chebyshev polynomial representation of imaginary-time response functions
,”
Phys. Rev. B
98
,
075127
(
2018
).
49.
H.
Shinaoka
,
J.
Otsuki
,
M.
Ohzeki
, and
K.
Yoshimi
, “
Compressing Green’s function using intermediate representation between imaginary-time and real-frequency domains
,”
Phys. Rev. B
96
,
035147
(
2017
).
50.
N.
Chikano
,
J.
Otsuki
, and
H.
Shinaoka
, “
Performance analysis of a physically constructed orthogonal representation of imaginary-time Green’s function
,”
Phys. Rev. B
98
,
035104
(
2018
).
51.
N.
Chikano
,
K.
Yoshimi
,
J.
Otsuki
, and
H.
Shinaoka
, “
irbasis: Open-source database and software for intermediate-representation basis functions of imaginary-time Green’s function
,”
Comput. Phys. Commun.
240
,
181
188
(
2019
).
52.
J.
Li
,
M.
Wallerberger
,
N.
Chikano
,
C.-N.
Yeh
,
E.
Gull
, and
H.
Shinaoka
, “
Sparse sampling approach to efficient ab initio calculations at finite temperature
,”
Phys. Rev. B
101
,
035144
(
2020
).
53.
M.
Kaltak
and
G.
Kresse
, “
Minimax isometry method
,” arXiv:1909.01740 [cond-mat.mtrl-sci] (
2019
).
54.
J. W.
Negele
and
H.
Orland
,
Quantum Many-Particle Systems
(
Westview Press
,
1998
).
55.
A. L.
Fetter
and
J. D.
Walecka
,
Quantum Theory of Many-Particle Systems
(
Dover Publications, Inc.
,
Mineloa, NY
,
2003
).
56.
A.
Altland
and
B.
Simons
,
Condensed Matter Field Theory
, 2nd ed. (
Cambridge University Press
,
Cambridge, UK
,
2010
).
57.
G.
Stefanucci
and
R.
van Leeuwen
,
Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction
(
Cambridge University Press
,
2013
).
58.
M.
Schüler
and
Y.
Pavlyukh
, “
Spectral properties from Matsubara Green’s function approach: Application to molecules
,”
Phys. Rev. B
97
,
115164
(
2018
).
59.
P. D.
Lax
,
Functional Analysis
(
John Wiley & Sons, Inc.
,
2002
).
60.
O.
Parcollet
,
M.
Ferrero
,
T.
Ayral
,
H.
Hafermann
,
I.
Krivenko
,
L.
Messio
, and
P.
Seth
, “
TRIQS: A toolbox for research on interacting quantum systems
,”
Comput. Phys. Commun.
196
,
398
415
(
2015
).
61.
L.-L. W.
Jie Shen
and
T.
Tang
,
Spectral Methods Algorithms, Analysis and Applications
, Springer Series in Computational Mathematics Vol. 41 (
Springer
,
2011
).
62.
N.
Hale
and
A.
Townsend
, “
An algorithm for the convolution of Legendre series
,”
SIAM J. Sci. Comput.
36
,
A1207
A1220
(
2014
).
63.
T. T.
Jie Shen
, in
Spectral and High-Order Methods With Applications
, Mathematics Monograph Series 3, edited by
C.
Yuzhuo
(
Science Press
,
Beijing
,
2006
).
64.
D.
Neuhauser
,
R.
Baer
, and
D.
Zgid
, “
Stochastic self-consistent second-order Green’s function method for correlation energies of large electronic systems
,”
J. Chem. Theory Comput.
13
,
5396
5403
(
2017
).
65.
Q.
Sun
,
T. C.
Berkelbach
,
N. S.
Blunt
,
G. H.
Booth
,
S.
Guo
,
Z.
Li
,
J.
Liu
,
J. D.
McClain
,
E. R.
Sayfutyarova
,
S.
Sharma
,
S.
Wouters
, and
G. K. L.
Chan
, “
PySCF: The python-based simulations of chemistry framework
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1340
(
2018
).
66.

This follows from the generalized imaginary time trace Tr[G]ξβab0βdτdτδa,bδ(ττ+0)Gab(τ,τ) of a response function Gab(τ,τ)Tca(τ)cb(τ).

67.
S.
Boys
and
F.
Bernardi
, “
The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors
,”
Mol. Phys.
19
,
553
566
(
1970
).
68.
T.
Van Mourik
,
A. K.
Wilson
, and
T. H.
Dunning
, “
Benchmark calculations with correlated molecular wavefunctions. XIII. Potential energy curves for He2, Ne2 and Ar2 using correlation consistent basis sets through augmented sextuple zeta
,”
Mol. Phys.
96
,
529
547
(
1999
).
69.
L. J.
Holleboom
and
J. G.
Snijders
, “
A comparison between the Møller–Plesset and Green’s function perturbative approaches to the calculation of the correlation energy in the many-electron problem
,”
J. Chem. Phys.
93
,
5826
5837
(
1990
).
70.
D.
Feller
, “
Application of systematic sequences of wave functions to the water dimer
,”
J. Chem. Phys.
96
,
6104
6114
(
1992
).
71.
T.
Helgaker
,
W.
Klopper
,
H.
Koch
, and
J.
Noga
, “
Basis-set convergence of correlated calculations on water
,”
J. Chem. Phys.
106
,
9639
9646
(
1997
).
72.
R. A.
Kendall
,
T. H.
Dunning
, and
R. J.
Harrison
, “
Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions
,”
J. Chem. Phys.
96
,
6796
6806
(
1992
).
73.
D. E.
Woon
and
T. H.
Dunning
, “
Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon
,”
J. Chem. Phys.
98
,
1358
1371
(
1993
).
74.
D. E.
Woon
and
T. H.
Dunning
, “
Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response properties
,”
J. Chem. Phys.
100
,
2975
2988
(
1994
).
75.
D. E.
Woon
and
T. H.
Dunning
, “
Benchmark calculations with correlated molecular wave functions. I. Multireference configuration interaction calculations for the second row diatomic hydrides
,”
J. Chem. Phys.
99
,
1914
1929
(
1993
).
76.
D. E.
Woon
and
T. H.
Dunning
, “
Calculation of the electron affinities of the second row atoms: Al–Cl
,”
J. Chem. Phys.
99
,
3730
3737
(
1993
).
77.
D. E.
Woon
, “
Accurate modeling of intermolecular forces: A systematic Møller-Plesset study of the argon dimer using correlation consistent basis sets
,”
Chem. Phys. Lett.
204
,
29
35
(
1993
).
78.
K. A.
Peterson
,
R. A.
Kendall
, and
T. H.
Dunning
, “
Benchmark calculations with correlated molecular wave functions. II. Configuration interaction calculations on first row diatomic hydrides
,”
J. Chem. Phys.
99
,
1930
1944
(
1993
).
79.
S. S.
Xantheas
and
T. H.
Dunning
, “
Theoretical studies of sulfurous species of importance in atmospheric chemistry. 1. Characterization of the mercaptooxy (HSO) and hydroxythio (SOH) isomers
,”
J. Phys. Chem.
97
,
6616
6627
(
1993
).
80.
D. E.
Woon
, “
Benchmark calculations with correlated molecular wave functions. V. The determination of accurate ab initio intermolecular potentials for He2, Ne2, and Ar2
,”
J. Chem. Phys.
100
,
2838
2850
(
1994
).
81.
D. E.
Woon
and
T. H.
Dunning
, “
Benchmark calculations with correlated molecular wave functions. VI. Second row A2 and first row/second row AB diatomic molecules
,”
J. Chem. Phys.
101
,
8877
8893
(
1994
).
82.
D. E.
Woon
,
T. H.
Dunning
, and
K. A.
Peterson
, “
Ab initio investigation of the N2–HF complex: Accurate structure and energetics
,”
J. Chem. Phys.
104
,
5883
5891
(
1996
).
83.
T.
van Mourik
and
T. H.
Dunning
, “
Ab initio characterization of the structure and energetics of the ArHF complex
,”
J. Chem. Phys.
107
,
2451
2462
(
1997
).
84.
K. A.
Peterson
and
T. H.
Dunning
, “
The CO molecule: The role of basis set and correlation treatment in the calculation of molecular properties
,”
J. Mol. Struct.: THEOCHEM
400
,
93
117
(
1997
), ab initio benchmark studies.
85.
K. A.
Peterson
,
A. K.
Wilson
,
D. E.
Woon
, and
T. H.
Dunning
, Jr.
, “
Benchmark calculations with correlated molecular wave functions XII. Core correlation effects on the homonuclear diatomic molecules B2-F2
,”
Theor. Chem. Acc.
97
,
251
259
(
1997
).
86.
Y.
Katznelson
,
Fundamental Principles of Optical Lithography
(
Dover
,
1976
).
87.
A. R.
DiDonato
, “
Recurrence relations for the indefinite integrals of the associated Legendre functions
,”
Math. Comput.
38
,
547
551
(
1982
).