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 He_{2} 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^{−9}*E*_{h} with only a few hundred expansion coefficients.

## I. INTRODUCTION

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 theory^{3} and its extensions,^{4–6} and diagrammatic Monte Carlo.^{7} *Ab 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 theory^{27–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 states^{39,40} is needed. Consequently, it is important to develop a compact representation that yields highly converged properties.

With the standard approach using equidistant Matsubara frequency^{41} 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 polynomials^{47,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 *N*_{L}, $\u223cO(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.

Domain . | Basis . | Convergence . | Compactness . | Dyson scaling . |
---|---|---|---|---|

Matsubara frequency | Finite frequency cutoff | $O(1)$ | Poor | $O(N)$ |

Tail correction, pth order^{42–44} | $O(N\u2212p)$ | Fair | … | |

Spline grid^{14} | … | Good | … | |

Both frequency and time | Sparse sampling^{52} | … | See Ref. 52 ^{a} | $O(N)$ |

Minimax isometry^{53} | … | See Ref. 53 | $O(N)$ | |

Imaginary time | Uniform mesh | $O(N\u22121)$ | Poor | $O(N3)$ |

Power mesh^{45,46,58} | … | Fair | $O(N3)$ | |

Orthogonal functions | Intermediate representation^{49–51} | $\u2272O(e\u2212N)$ | Excellent | No |

Chebyshev polynomials^{48} | $\u2272O(e\u2212N)$ | Very good | $O(N3)$ | |

Legendre polynomials (this work) | $\u2272O(e\u2212N)$^{47} | Very good | $O(N2)$ |

Domain . | Basis . | Convergence . | Compactness . | Dyson scaling . |
---|---|---|---|---|

Matsubara frequency | Finite frequency cutoff | $O(1)$ | Poor | $O(N)$ |

Tail correction, pth order^{42–44} | $O(N\u2212p)$ | Fair | … | |

Spline grid^{14} | … | Good | … | |

Both frequency and time | Sparse sampling^{52} | … | See Ref. 52 ^{a} | $O(N)$ |

Minimax isometry^{53} | … | See Ref. 53 | $O(N)$ | |

Imaginary time | Uniform mesh | $O(N\u22121)$ | Poor | $O(N3)$ |

Power mesh^{45,46,58} | … | Fair | $O(N3)$ | |

Orthogonal functions | Intermediate representation^{49–51} | $\u2272O(e\u2212N)$ | Excellent | No |

Chebyshev polynomials^{48} | $\u2272O(e\u2212N)$ | Very good | $O(N3)$ | |

Legendre polynomials (this work) | $\u2272O(e\u2212N)$^{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 $\u223cO(NL2)$, one order better than previous formulations.^{48} The super-exponential convergence allows us to achieve an energy accuracy of 10^{−9}*E*_{h} 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 He_{2} 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 He_{2}. In Sec. VI, we present conclusions.

## II. THEORY

The imaginary time single particle Green’s function *G* is defined on the interval *τ* ∈ [−*β*, *β*], *G* ≡ *G*(*τ*), 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*(*τ*) is^{54–57}

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 type^{59} imaginary time convolution is defined as $\Sigma *G\u2261\u222b0\beta d\tau \xaf\u2009\Sigma (\tau \u2212\tau \xaf)G(\tau \xaf)$.

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

where the Matsubara frequencies *iω*_{n} are given by $i\omega n\u2261i\pi \beta (2n+\eta )$ with *η* = (1 − *ξ*)/2 and *n* integers.^{54–57} In Matsubara frequency space, the Dyson equation [Eq. (1)] is diagonal,^{41}

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

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

where *L*_{n}[*x*] are Legendre polynomials defined on *x* ∈ [−1, 1] and $x(\tau )=2\tau \beta \u22121$. The Legendre coefficients *G*_{n} have a faster than exponential asymptotic decay^{47} [see Fig. 1(c)]. This also causes the finite *N*_{L} expansion at the right-hand side of Eq. (3) to converge faster than exponential $\u2272O(e\u2212NL)$ to the analytical *G*(*τ*).

## III. LEGENDRE SPECTRAL METHOD

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 *G*_{n} [Eq. (3)]. In the space of a finite Legendre expansion of order *N*_{L}, Eq. (1) is cast to a linear equation system,

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 $\u223cO(NL2)$, one order better than previous approaches.^{48}

Hence, the derivative matrix *D*_{kn} in Eq. (4) is given by

and is upper triangular (see Fig. 3). Using *L*_{n}(±1) = (±1)^{n}, the Dyson equation boundary condition can be written as

### A. Spectral convolution

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

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\u2276$, representing the two Volterra terms [Eq. (8)],

Stable recursion relations for $Bnk\u2276$ have been derived by Hale and Townsend^{62} 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,

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

The two first columns are given by the starting relations,

with the special case for *k* = 0, $B0,1\u2276=\u2213B1,0\u2276/3$, using the Legendre coefficients Σ_{n} of the self-energy Σ [cf. Eq. (3)].

### B. Convergence and scaling

Since each coefficient in $Bkn\u2276$ can be computed in $O(1)$ operations, the scaling of the convolution matrix construction is $\u223cO(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 *N*_{L} (see Fig. 2). This is in stark contrast to the polynomial convergence of the standard Matsubara tail approach^{42–44} (also shown in Fig. 2).

### C. Imaginary time transform

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,

To construct the well-conditioned transform matrices *S*_{ni} and *L*_{in}, we employ the Legendre quadrature and the Legendre–Gauss–Lobatto points $xi\u2208{x:(1\u2212x2)LNL(x)=0}$, *x*_{0} = −1, *x*_{N} = 1, re-scaled to the imaginary time interval [0, *β*], $\tau i=\beta xi+12$. Using *x*_{i}, the matrices *S*_{ni} and *L*_{in} can be directly constructed (avoiding matrix inversion),

## IV. APPLICATION (GF2)

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 applications^{8,15} and in quantum chemistry^{11,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

where *i*, *j*, *k* are orbital indices, *S*_{ij} is the overlap matrix, and *F*_{ij} is the so-called Fock matrix, $Fij\u2261hij+\Sigma ij(HF)$. The boundary condition for this equation is *∑*_{j}(*G*_{ij}(0) − *ξG*_{ij}(*β*)) · *S*_{jk} = −**1**_{ik}. Here, the single particle term *h*_{ij} accounts for electronic kinetic and nuclear-electronic matrix elements, and the Hartree–Fock self-energy $\Sigma ij(HF)$ is given by

where *P*_{ij} is the density matrix *P*_{ij} = −2*G*_{ij}(*β*) and $v$_{ijkl} 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

The evaluation of Σ^{(GF2)}(*τ*) for fixed *τ* scales as $\u223cO(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}.

The total energy *E* of the system is given by

## V. RESULTS

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 *D*_{e} in noble-gas dimers, where the weak bonding requires high precision calculations of total energies. The binding energy *D*_{e} is obtained from the minimum of the interaction energy *E*_{int}(*r*) as a function of atomic separation *r*,

where *r*_{e} is the equilibrium atomic distance. The interaction energy *E*_{int} is, in turn, given by

where $EA2$ is the total energy of the dimer and *E*_{A} is the total energy of the single atom (the monomer) evaluated using the standard counterpoise correction.^{67} In the noble gases, the total energies *E*_{A} and $EA2$ are of the order of Hartrees (∼*E*_{h} ≡ 1 Hartree), while the binding energy *D*_{e} is of the order of tens of micro-Hartrees (∼10*μE*_{h}), hence requiring high precision calculation of the total energies.

We use He_{2} as a prototype system since there exist published reference results for the binding energy *D*_{e} and equilibrium distance *r*_{e} 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 *E*_{int}(*r*) (and −*D*_{e}) of He_{2} computed with HF, MP2, and GF2 in the aug-cc-pvqz basis together with CCSD and CCSD(T) reference results on *D*_{e}.^{68} The GF2 results are obtained by fitting a 4th order polynomial to 21 *r*-points of *E*_{int}(*r*) computed in a 0.1 Bohr range centered around the minimum at *r*_{e}. 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, He_{2} 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.

### A. Complete basis set limit

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-pv*n*z) 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 *D*_{e} and equilibrium distance *r*_{e} of He_{2} 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 $\beta =50\u2009Eh\u22121$ 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 $\beta =50\u2009Eh\u22121$) 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.

D_{e} (μE_{h})
. | MP2 . | CCSD . | CCSD(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 |

D_{e} (μE_{h})
. | MP2 . | CCSD . | CCSD(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 |

r_{e} (Bohr)
. | MP2 . | CCSD . | CCSD(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 |

r_{e} (Bohr)
. | MP2 . | CCSD . | CCSD(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, $\u223cO(N\tau \u22c5N5)$, is comparable to the scaling of CCSD, $\u223cO(N6)$. As given in Table II, the accuracy of the GF2 result for *D*_{e} is comparable to that of CCSD when compared to that of CCSD(T), while the CCSD result for *r*_{e} is closer to the CCSD(T) result than that of GF2. This makes GF2 a considerable improvement over MP2.

With the systematic convergence of *D*_{e} and *r*_{e} as a function of basis set size *n*, it is possible to extrapolate to the complete basis limit *n* → *∞*.^{68} We extrapolate *D*_{e} and *r*_{e} using our GF2 aug-ccpv{t,q,5}z results by fitting the exponential model: *A* · *e*^{−B(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 *D*_{e} ≈ 29.67*μE*_{h} and *r*_{e} ≈ 5.680 *a*_{0} (see also Table II).

## VI. CONCLUSION AND OUTLOOK

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 *N*_{L} 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 $\u223c\u2009O(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 He_{2} using self-consistent second order perturbation theory (GF2). The exponential convergence of our algorithm allows us to reach the required 10^{−9}*E*_{h} zero temperature total-energy precision using only 100–200 Legendre coefficients in the Dunning basis series aug-ccpv*n*z.^{72–74}

## ACKNOWLEDGMENTS

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.

### APPENDIX: CONVOLUTION MATRIX

#### 1. Convolution and Fourier transform

The convolution of two continuous integrable functions is defined as^{62}

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

which satisfy the Fourier inversion theorem $F\u22121{F{f}}=f$ and convolution theorem,^{86}

#### 2. Legendre polynomials

The Legendre polynomials *P*_{n}(*x*) can be defined recursively using the three term recurrence relation,

They are orthogonal on [−1, 1],

and the derivatives satisfy the recurrence relation,

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

where *j*_{n}(*z*) is the *n*th spherical Bessel function and *P*_{n} = 0 outside [−1, 1].

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

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,

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 *f*_{M}(*x*) and *g*_{N}(*x*) of sufficiently high degree. With two Legendre series *f*_{M}(*x*) and *g*_{N}(*x*) supported on *x* ∈ [−1, 1],

Equation (A1) becomes

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* $\u2208$ [−2, 0]

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

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

collecting terms, we can write $\gamma k<=\u2211n=0NBk.n<\beta n$, where

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

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

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

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

##### b. Second interval *x* $\u2208$ [0, 2]

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

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

collecting terms, we can write $\gamma k>=\u2211n=0NBk.n>\beta n$, where

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

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,

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

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,

#### 3. Initial values $Bk,0\u2276$ and $Bk,1\u2276$

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,

for *a* = ±1, we get

where we have used *P*_{n}(±1) = (±1)^{n} to cancel the constant terms.

Returning to the convolution matrices, we have for $Bk,n<$ and *n* = 0, using *P*_{0}(*x*) = 1,

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

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

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

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 (*x* − *t* + 1), instead, of (*x* − *t* − 1) in Eq. (A35), so the sign before *B*_{k,0} is changed to +1. By using Eq. (A33), we obtain the recursion relation,

with the special case for *k* = 0, $B0,1\u2276=\u2213B1,0\u2276/3$.

## REFERENCES

_{0}W

_{0}for molecular systems

This follows from the generalized imaginary time trace $Tr[G]\u2261\u2212\xi \beta \u2211ab\u222c0\beta \u2009d\tau d\tau \u2032\u2009\delta a,b\delta (\tau \u2212\tau \u2032+0\u2212)Gab(\tau ,\tau \u2032)$ of a response function $Gab(\tau ,\tau \u2032)\u2261\u2212\u27e8Tca(\tau )cb\u2020(\tau \u2032)\u27e9$.

_{2}, Ne

_{2}and Ar

_{2}using correlation consistent basis sets through augmented sextuple zeta

_{2}, Ne

_{2}, and Ar

_{2}

_{2}and first row/second row AB diatomic molecules

_{2}–HF complex: Accurate structure and energetics

*ab initio*benchmark studies.

_{2}-F

_{2}