A “moment-conserving” reformulation of GW theory

Scott, Charles J. C.; Backhouse, Oliver J.; Booth, George H.

doi:10.1063/5.0143291

We show how to construct an effective Hamiltonian whose dimension scales linearly with system size, and whose eigenvalues systematically approximate the excitation energies of GW theory. This is achieved by rigorously expanding the self-energy in order to exactly conserve a desired number of frequency-independent moments of the self-energy dynamics. Recasting GW in this way admits a low-scaling $O [N^{4}]$ approach to build and solve this Hamiltonian, with a proposal to reduce this further to $O [N^{3}]$ ⁠. This relies on exposing a novel recursive framework for the density response moments of the random phase approximation, where the efficient calculation of its starting point mirrors the low-scaling approaches to compute RPA correlation energies. The frequency integration of GW, which distinguishes so many different GW variants, can be performed without approximation directly in this moment representation. Furthermore, the solution to the Dyson equation can be performed exactly, avoiding analytic continuation, diagonal approximations, or iterative solutions to the quasiparticle equation, with the full-frequency spectrum obtained from the complete solution of this effective static Hamiltonian. We show how this approach converges rapidly with respect to the order of the conserved self-energy moments and is applied across the GW100 benchmark dataset to obtain accurate GW spectra in comparison to traditional implementations. We also show the ability to systematically converge all-electron full-frequency spectra and high-energy features beyond frontier excitations, as well as avoiding discontinuities in the spectrum, which afflict many other GW approaches.

I. INTRODUCTION

Despite the phenomenal success of density functional theory (DFT) in electronic structure, its standard approach is both conceptually and (often) practically ill-suited for an accurate description of the energy levels in a material or chemical system.¹ These quantities are however essential for predictions of fundamental bandgaps and other charged excitation properties that govern the photo-dynamics, transport, and response properties of a system. Into this, GW theory has grown in popularity, first for materials and more recently for molecular systems, as a post-mean-field approach to obtain charged excitation spectra in a principled diagrammatic fashion, free from empiricism.^2–14

The GW approach is based on Hedin’s equations^15,16 and in its most common formulation builds a self-energy to dress a reference description of the quasi-particles of a system [generally from DFT or Hartree–Fock (HF)] with an infinite resummation of all “bubble” diagrams. These diagrams make up the random phase approximation (RPA)^17–20 and physically describes all collective quantum charge fluctuations in the electron density from the reference state arising from their correlated mutual Coulomb repulsion. This dynamically screens the effective interaction between the constituent quasi-particles of the system, whose physics generally dominates in small gapped semi-conducting systems. The use of this RPA screened interaction in GW has therefore become widespread, correcting many of the failures of DFT for spectral properties.

At the core of GW theory is a convolution, between the Green’s function of the system, G(ω), and the screened Coulomb interaction W_p(ω), obtained (in general) at the RPA level of theory. This provides the dynamical part of the self-energy, Σ(ω), formally written as Σ(ω) = (i/2π)∫dω′e^iηω′G(ω + ω′)W_p(ω′). There are many different variants of GW theory,^{3,13,14,21–44} which primarily differ due to (i) the choice (or absence) of self-consistency conditions on G(ω) and/or W_p(ω),²¹ (ii) the approach to find the quasi-particle energies once Σ(ω) is obtained (i.e., application of Dyson’s equation), and (iii) different approximations to perform the frequency integration in the convolution itself. A numerically exact formulation of this convolution entails an $O [N^{6}]$ scaling step, required to find the entire set of poles in W_p(ω) from the RPA. However, there are a number of techniques to approximate this frequency integration that can reduce this scaling (generally down to $O [N^{3} - N^{4}]$ ⁠) based on plasmon pole approximations, analytic continuation, contour deformation, and explicit grid resolved approaches for the dynamics of these quantities, among others. All of these approaches compress and approximate the dynamical resolution of the key quantities in order to simplify the resulting convolutional integral.

Another key difference between the approaches rests on how the quasi-particle energies are updated from their mean-field molecular orbital (MO) energy starting point, once the self-energy has been constructed. This is formally an application of Dyson’s equation but is commonly approximated via a self-consistent solution to a quasi-particle (QP) equation entailing a diagonal approximation to the self-energy. This is valid when the quasiparticle energies are far from the self-energy poles, thereby asserting that the GW largely just provides a shift in the original MO energies, rather than introducing significant quasiparticle renormalization, additional satellite peaks from state splitting, or relaxation of the mean-field electron density. These assumptions can however break down (especially in more correlated systems), while the numerical solution of the QP equation can also converge to different (“spurious”) results based on the specifics of how it is solved. A thorough study of the discrepancies due to different approaches to both the QP equation solution and the frequency integral in the convolution can be found in Ref. 12.

In this work, we introduce a new approach to this frequency integration and reformulate GW theory with a number of desirable properties, while retaining a low scaling. The key step is that the self-energy is not represented as an explicitly dynamical quantity, but instead in terms of a series of static matrices representing the moments of its frequency-dependence up to a given order. These can be directly obtained, and from them a compressed representation of the full self-energy can be algebraically constructed, which only has a number of poles that scales linearly with system size, but nonetheless exactly preserves the moment distribution of the exact GW self-energy dynamics up to the desired order.^45–49 This order can be systematically increased to more finely resolve the full dynamical dependence of the GW self-energy. The dynamical information is therefore implicitly recast into a small number of static matrices, each of which can be obtained in $O [N^{4}]$ time (with a proposed $O [N^{3}]$ algorithm also given). This removes the need for the definition of any frequency or time grids in which to resolve dynamical quantities, spectral broadening, finite temperatures, Fourier transforms, or analytic continuation, with all dynamics implicitly represented by this series of static quantities.

Furthermore, once these spectral moments of the self-energy are obtained, and the QP equation and restrictions to a diagonal self-energy representation can be entirely removed, with an exact application of Dyson’s equation possible in this “moment” representation. This leads to the self-energy represented as a small number of explicit poles at specific energies which taken together have exactly the moment distribution in their dynamics as described. This allows for a simple construction of the full frequency dependence of the resulting quasi-particle spectrum (including any additional emergent satellite structures from the correlations) via diagonalization in an “upfolded” and explicit Hamiltonian representation. Moment expansions have a long history in the representation of dynamical quantities, with the use in numerical approaches,^45,50 characterizing sum rules,^51,52 and as physical observables in their own right.^53–55

In Sec. II, we show how these spectral moments of the self-energy can be directly constructed from moments of Green’s function and the two-point density–density response function from RPA, without any further approximation. We show how these can be used to directly obtain the full-frequency GW spectrum without the requirement of an explicit grid. In Sec. III, we show how the RPA can be fully reformulated as a series expansion of moments of the density–density (dd) response, and in Sec. IV, we show how they can be efficiently obtained in $O [N^{4}]$ cost, based on ideas from the seminal work in 2010 by Furche and collaborators on low-scaling approaches for the RPA correlation energy.⁵⁶ Furthermore, in Sec. V, we propose an approach to further reduce the scaling of the whole algorithm directly (rather than asymptotically) to cubic cost with system size, without invoking screening or locality assumptions. We then apply the approach in Sec. VI to the commonly used molecular GW100 test set frequently used to benchmark GW implementations, demonstrating a rapid convergence of the moment expansion, and accurate and efficient results across this test set for the G₀W₀ level of theory.

II. MOMENT-TRUNCATED GW THEORY

In GW theory, the dynamical part of the self-energy, obtained as the convolution of the G(ω) and W_p(ω), can be formally expanded as a sum over the

O [N^{2}]

neutral excitations of RPA theory (representing the poles of the screened Coulomb interaction) and the charged excitations of the reference Green’s function. In the absence of self-consistency (the most common “G₀W₀” formulation of the method, which we exclusively consider in this work), the Green’s function is just given from the

O [N]

mean-field molecular orbital energies. This allows the self-energy to be explicitly evaluated in the frequency-domain as

\begin{align} Σ_{p q} (ω) & = \sum_{ν} \sum_{i a, j b, k} \frac{(p k | i a) (X_{i a}^{ν} + Y_{i a}^{ν}) (X_{j b}^{ν} + Y_{j b}^{ν}) (q k | j b)}{ω - (ϵ_{k} - Ω_{ν}) - i 0^{+}} \\ + \sum_{ν} \sum_{i a, j b, c} \frac{(p c | i a) (X_{i a}^{ν} + Y_{i a}^{ν}) (X_{j b}^{ν} + Y_{j b}^{ν}) (q c | j b)}{ω - (ϵ_{c} + Ω_{ν}) + i 0^{+}} . \end{align}

(1)

In these expressions, Ω_ν are the neutral excitation energies from RPA theory, and

X_{i a}^{ν}

and

Y_{i a}^{ν}

are the excitation and de-excitation amplitude components associated with this excitation. These are expanded in the basis of hole and particle spin-orbitals of the reference mean-field, represented by the indices i, j, k (a, b, c) of dimension o (v), respectively, and with orbital energies denoted by ϵ_x. The bare two-electron integrals are denoted by (pk|ia) in standard Mulliken (“chemists”) notation, which are therefore screened by the RPA reducible density response. More details on these quantities are given in Sec. III. The first term in Eq. therefore represents the “lesser” part, and the second term the “greater” part of the full G₀W₀ self-energy.

Exact evaluation of the RPA excitations (Ω_ν) scales as

O [N^{6}]

⁠, rendering it unsuitable for large-scale implementations. However, in this work, we are only interested in evaluating the spectral moments of the resulting self-energy and finding a resulting compressed representation of the self-energy with fewer poles, but which by construction matches the spectral moments up to a desired order. These frequency-independent spectral moments are defined separately for the greater and lesser parts and represent the nth-order moments of the resulting dynamical distributions, as

Σ_{p q}^{(n, <)} = - \frac{1}{π} \int_{- \infty}^{μ} I m [Σ {(ω)}_{p q}] ω^{n} d ω,

(2)

= {(- 1)}^{n} {\frac{d^{n} Σ {(τ)}_{p q}}{d τ^{n}}|}_{τ = 0^{+}},

(3)

and similarly

Σ_{p q}^{(n, >)} = \frac{1}{π} \int_{μ}^{\infty} I m [Σ {(ω)}_{p q}] ω^{n} d ω

(4)

= {(- 1)}^{n} {\frac{d^{n} Σ {(τ)}_{p q}}{d τ^{n}}|}_{τ = 0^{-}},

(5)

where μ represents the chemical potential of the system. This exposes the relationship of these spectral moments to a Taylor expansion of the short-time dynamics of the greater and lesser parts of the self-energy, with the moments defining the integrated weight, mean, variance, skew, and higher-order moments of the dynamical distribution of each element of the self-energy in the frequency domain.

Applied to the GW self-energy of Eq. , the moments can be constructed as

\begin{align} Σ_{p q}^{(n, <)} & = \sum_{ν} \sum_{i a, j b} \sum_{k} [(p k | i a) (X_{i a}^{ν} + Y_{i a}^{ν}) \\ \times {(ϵ_{k} - Ω_{ν})}^{n} (X_{j b}^{ν} + Y_{j b}^{ν}) (q k | j b)], \end{align}

(6)

\begin{align} Σ_{p q}^{(n, >)} & = \sum_{ν} \sum_{i a, j b} \sum_{c} [(p c | i a) (X_{i a}^{ν} + Y_{i a}^{ν}) \\ \times {(ϵ_{c} + Ω_{ν})}^{n} (X_{j b}^{ν} + Y_{j b}^{ν}) (q c | j b)] . \end{align}

(7)

The moment distribution of a convolution of two quantities can be expressed via the binomial theorem as a sum of products of the moments of the individual quantities. This enables us to split apart the expressions above into products of the individual Green’s function and density–density response moments. Defining the nth-order spectral moments of the RPA density-response, summed over both particle–hole excitation and de-excitation components, as

η_{i a, j b}^{(n)} = \sum_{ν} (X_{i a}^{ν} + Y_{i a}^{ν}) Ω_{ν}^{n} (X_{j b}^{ν} + Y_{j b}^{ν}),

(8)

we can rewrite Eqs. and as

Σ_{p q}^{(n, <)} = \sum_{i a, j b, k} \sum_{t = 0}^{n} (\binom{n}{t}) {(- 1)}^{t} ϵ_{k}^{n - t} (p k | i a) η_{i a, j b}^{(t)} (q k | j b),

(9)

Σ_{p q}^{(n, >)} = \sum_{i a, j b, c} \sum_{t = 0}^{n} (\binom{n}{t}) ϵ_{c}^{n - t} (p c | i a) η_{i a, j b}^{(t)} (q c | j b) .

(10)

Evaluating the self-energy spectral moments of Eqs. and up to a desired order n, represents the central step of the proposed “moment-conserving” GW formulation, defining the convolution between G(ω) and W_p(ω) in this moment expansion of the dynamics. In Sec. III, we show how the RPA can be reformulated to define specific constraints on the relations between different orders of the RPA density-response moments,

η_{i a, j b}^{(n)}

⁠. These relations are subsequently used in Sec. IV to demonstrate how the self-energy moments can be evaluated in

O [N^{4}]

scaling, with Sec. V going further to propose a cubic scaling algorithm for their evaluation (and therefore full GW algorithm). In addition to these moments representing the dynamical part of the self-energy, we also require a static (exchange) part of the self-energy, Σ_∞, which can be calculated as

Σ_{\infty} = K [D] - V_{xc},

(11)

where K[D] is the exchange matrix evaluated with the reference density matrix. This reference density matrix is found from a prior mean-field calculation via a self-consistent Fock or Kohn–Sham single-particle Hamiltonian, f[D], with V_xc being the exchange-correlation potential used in f. Note that for a Hartree–Fock reference, this static self-energy contribution is zero.

A. Full GW spectrum from self-energy moments

Once the moments of the self-energy are found, it is necessary to obtain the resulting dressed GW excitations and spectrum. While this is formally an application of Dyson’s equation, the most common approach is to find each GW excitation explicitly via a self-consistent solution (or linearized approximation) of the quasiparticle equation, while assuming a diagonal self-energy in the MO basis.¹² This assumption neglects physical effects due to electron density relaxation and mixing or splitting of quasiparticle states in more strongly correlated systems. In this work, we allow for an exact invocation of Dyson’s equation, which can be achieved straightforwardly in this moment domain of the effective dynamics, allowing extraction of quasi-particle weights associated with transitions, and a full matrix-valued form of the resulting GW Green’s function over all frequencies, with all poles obtained analytically without artificial broadening.

This is achieved by constructing an “upfolded” representation of an effective Hamiltonian, consisting of coupling between a physical and “auxiliary” space (with the latter describing the effect of the moment-truncated self-energy). Specifically, we seek an effective static Hamiltonian,

\tilde{H} = [\begin{matrix} f + Σ_{\infty} & \tilde{W} \\ {\tilde{W}}^{†} & \tilde{d} \end{matrix}],

(12)

whose eigenvalues are the charged excitation energies at the level of the moment-truncated GW, with quasiparticle weights and Dyson orbitals explicitly obtained from the projection of the corresponding eigenvectors into the physical (MO) space. The full Green’s function can therefore be constructed as

G (ω) = {(ω I - f + Σ_{\infty} - \tilde{W} {(ω I - \tilde{d})}^{- 1} {\tilde{W}}^{†})}^{- 1} .

(13)

Such upfolded representations have been considered previously in diagrammatic theories, in a recasting of GF2 theory in terms of its moments^46,47,57 as well as more recently to GW among others.^48,58–63 For “exact” G₀W₀, this auxiliary space (i.e., the dimension of d) must scale as

O [N^{3}]

⁠.^60,61 However, in the moment truncation,

\tilde{W}

and

\tilde{d}

can be directly constructed such that their effect exactly matches that of a truncated set of conserved GW self-energy moments (separately in the particle and hole sectors), yet rigorously scales in dimension as

O [n N]

⁠, where n is the number of conserved self-energy moments. This allows for a complete diagonalization of

\tilde{H}

⁠, obtaining all excitations in a single shot, and a reconstruction of the full GW Green’s function from its Lehmann representation in

O [{(n N)}^{3}]

computational effort, avoiding the need for any grids or iterative solutions once

\tilde{H}

is found.

To find this effective upfolded representation of the moment-conserving dynamics, we modify the block Lanczos procedure to ensure the construction a

\tilde{H}

of minimal size, whose effective hole and particle self-energy moments exactly match the ones from Eqs. and . We first proceed by splitting the auxiliary space into a space denoting the effect of the hole (lesser) and particle (greater) self-energy, and consider each in turn. Focusing on the lesser self-energy, we can construct an exact upfolded self-energy representation,⁶³ via inspection from Eq. , with

W_{p, k ν} = \sum_{i a} (p k | i a) (X_{i a}^{ν} + Y_{i a}^{ν}),

(14)

d_{k ν, l ν^{'}} = (ϵ_{k} - Ω_{ν}) δ_{k, l} δ_{ν, ν^{'}},

(15)

where we remove the tilde above upfolded auxiliary quantities when denoting the exact upfolded GW self-energy components. We now consider the projection of the exact GW upfolded matrix representation into a truncated block tridiagonal form, as

\begin{align} {\tilde{H}}_{tri} & = {\tilde{q}}^{(j), †} [\begin{matrix} f + Σ_{\infty} & W \\ W^{†} & d \end{matrix}] {\tilde{q}}^{(j)} \\ = [\begin{matrix} f + Σ_{\infty} & L & 0 \\ L^{†} & M_{1} & C_{1} \\ C_{1}^{†} & M_{2} & C_{2} \\ C_{2}^{†} & M_{3} & ⋱ \\ ⋱ & ⋱ & C_{j - 1} \\ 0 & C_{j - 1}^{†} & M_{j} \end{matrix}], \end{align}

(16)

where we define

{\tilde{q}}^{(j)}

as

{\tilde{q}}^{(j)} = [\begin{matrix} I & 0 \\ 0 & q^{(j)} \end{matrix}] .

(17)

The q^(j) are block Lanczos vectors of depth j, which form a recursive Krylov space as

q^{(j)} = [\begin{matrix} q_{1} & q_{2} & \dots & q_{j} \end{matrix}]

⁠, ensuring that when taken together, they project to a block tridiagonal representation of the upfolded self-energy with j on-diagonal blocks as shown. The action of this block Lanczos tridiagonalization of the upfolded (hole or particle) self-energy is to exactly conserve these spectral moments of the self-energy.^64,65 This block tridiagonal representation is equivalent to a truncated continued fraction form,^50,66 widely used in the representation of dynamical quantities,^67,68 and even as an expansion previously considered within GW theory.⁴⁵ We therefore seek to reformulate the Lanczos recursion in terms of just these moments, rather than the action of the full upfolded Hamiltonian, which we seek to avoid due to its scaling.

The initial couplings L and Lanczos vectors q₁ can be found via a QR factorization of the exact GW couplings W, as

W^{†} = q_{1} L^{†} .

(18)

However, this will scale poorly, and so we can rewrite this to directly compute L from the computed self-energy moments, rather than requiring manipulations of the full auxiliary space. Via the Cholesky QR algorithm,^69,70 we can relate L to the zeroth order self-energy moment as

L^{†} = {(W W^{†})}^{\frac{1}{2}} = {(Σ^{(0)})}^{\frac{1}{2}},

(19)

where the indication of the sector of the self-energy has been dropped, with this process considered independently for the hole and particle (lesser and greater, respectively) parts of the self-energy. The initial Lanczos vector can then be computed as q₁ = W^†L^−1,†. Subsequent block Lanczos vectors can then be defined according to the standard three-term recurrence

q_{i + 1} C_{i}^{†} = [d q_{i} - q_{i} M_{i} - q_{i - 1} C_{i - 1}],

(20)

where the on-diagonal blocks are defined as

M_{i} = q_{i}^{†} d q_{i} .

(21)

In order to recast this process in terms of the self-energy moments directly, we wish to express the block Lanczos recurrence in terms of the inner space of the Lanczos vectors rather than spanning the large auxiliary space. The choice of initial vectors in Eq. permits the definition

\begin{align} S_{1, 1}^{(n)} & = q_{1}^{†} d^{n} q_{1} \\ = L^{- 1} Σ^{(n)} L^{- 1, †}, \end{align}

(22)

where the subscript indices on S indicate the projection into the block Lanczos space on the left and the right, respectively. These provide the initialization of the recurrence terms, and have a dimension that scales linearly with system size (the same as the input self-energy moments). These Σ⁽ⁿ⁾ matrices are therefore the input to the procedure, defined by Eqs. and . One can then use the definition of the three-term recurrence in Eq. to express all definitions in terms of these moments, without formal reference to the large auxiliary space quantities (d or W), as

S_{i + 1, i}^{(n)} = q_{i + 1}^{†} d^{n} q_{i} = C_{i}^{- 1} [S_{i, i}^{(n + 1)} - M_{i} S_{i, i}^{(n)} - C_{i - 1}^{†} S_{i - 1, i}^{(n)}],

(23)

\begin{align} S_{i + 1, i + 1}^{(n)} & = q_{i + 1}^{†} d^{n} q_{i + 1} \\ = C_{i}^{- 1} [S_{i, i}^{(n + 2)} + M_{i} S_{i, i}^{(n)} M_{i} + C_{i - 1}^{†} S_{i - 1, i - 1}^{(n)} C_{i - 1} \\ - P (S_{i, i}^{(n + 1)} M_{i}) + P (M_{i} S_{i, i - 1}^{(n)} C_{i - 1}) - P (S_{i, i - 1}^{(n + 1)} C_{i - 1})] C_{i}^{- 1, †}, \end{align}

(24)

where the permutation operator P is defined as P(Z) = Z + Z^†. For a Hermitian theory, we can write

S_{i, j} = S_{j, i}^{†}

and assume the zeroth order Lanczos vector to be zero, i.e.,

S_{i, 0}^{(n)} = S_{0, j}^{(n)} = S_{0, 0}^{(n)} = 0

⁠. Additionally, the orthogonality of the Lanczos vectors requires that

S_{i, j}^{(0)} = δ_{i j} I

⁠. By considering again the Cholesky QR algorithm, the off-diagonal C matrices can therefore be computed as

C_{i}^{2} = [S_{i, i}^{(2)} + M_{i}^{2} + C_{i - 1}^{†} C_{i - 1} - P (S_{i, i}^{(1)} M_{i}) - P (S_{i, i - 1}^{(1)} C_{i - 1})],

(25)

and the on-diagonal M matrices can be found using Eq. ,

M_{i} = q_{i}^{†} d q_{i} = S_{i, i}^{(1)} .

(26)

These recurrence relations allow the calculation of the on- and off-diagonal blocks resulting in the block tridiagonal form of the Hamiltonian in Eq. . Despite the apparent complexity of the recurrence relations, this algorithm contains no step scaling greater than

O [N^{3}]

⁠, by eliminating the explicit reference to the full upfolded Hamiltonian, while still conserving the exact self-energy moments by construction.

It can be seen from Eq. that the auxiliary space formally couples to both the hole and particle physical sectors (both occupied and virtual MOs). To allow for this coupling (which is critical for generating effective higher order diagrams and to avoid a “non-Dyson” approximation^{60,61,71–76}), the solution of the Dyson equation using the compressed self-energy, i.e., diagonalization of Eq. , requires the combination of the block tridiagonal Hamiltonian in Eq. resulting from both the hole and particle self-energy moments. This is constructed as

\tilde{H} = [\begin{matrix} f + Σ_{\infty} & \tilde{W} \\ {\tilde{W}}^{†} & \tilde{d} \end{matrix}] = [\begin{matrix} f + Σ_{\infty} & {\tilde{W}}^{<} & {\tilde{W}}^{>} \\ {\tilde{W}}^{<, †} & {\tilde{d}}^{<} & 0 \\ {\tilde{W}}^{>, †} & 0 & {\tilde{d}}^{>} \end{matrix}],

(27)

where

\tilde{W}

are equal to the L matrix padded by zeros

{\tilde{W}}^{≶} = [\begin{matrix} L^{≶} & 0 & \dots & 0 \end{matrix}],

(28)

and

\tilde{d}

are defined by the block tridiagonal elements

{\tilde{d}}^{≶} = [\begin{matrix} M_{1}^{≶} & C_{1}^{≶} & 0 \\ C_{1}^{≶, †} & M_{2}^{≶} & ⋱ \\ ⋱ & ⋱ & C_{j - 1}^{≶} \\ 0 & C_{j - 1}^{≶, †} & M_{j}^{≶} \end{matrix}] .

(29)

This ensures the conservation of both the separate hole and particle moments of the self-energy, as well as conservation in the central moments according to their sum.

The compressed Hamiltonian can be returned to a diagonal representation of the self-energy by diagonalizing $\tilde{d}$ and appropriately rotating $\tilde{W}$ into this basis. The eigenvalues of $\tilde{H}$ are moment-conserving approximations to those of the exact upfolded Hamiltonian, and the corresponding eigenvectors u can be transformed into Dyson orbitals via LPu, where P is a projection into the physical space, and the L is required to transform the physical component of the eigenvectors back to the MO representation. This process conserves exactly the first 2j hole and particle self-energy moments. Commonly, a notation referring to the number of iterations of the block Lanczos recurrence n_iter is used; in this notation, the n_iter = 0 calculation corresponds to the inclusion of only a single on-diagonal block M₁, with modified couplings L to the physical space. As such, in this notation, the number of conserved moments equals 2n_iter + 2, i.e., up to and including the 2n_iter + 1 order moment. This is the same number of moments as required as input to the recurrence relations, and therefore the algorithm conserves all the moments used as input, which should be up to an odd order. After n_iter applications of this algorithm to both the lesser and greater self-energy sectors, this results in N(2n_iter + 3) quasiparticle states (demonstrating the potential to capture satellite features with these additional poles). As such, application of this algorithm becomes theoretically equivalent to a full diagonalization of the exact upfolded Hamiltonian in the limit of n_iter ∼ N².

III. DENSITY RESPONSE MOMENTS IN THE RPA

Having described the overall approach in Sec. II, what remains for a practical implementation is to ensure that the GW self-energy moments described in Eqs. (9) and (10) can be computed efficiently. As a first step toward this, in this section, we show how the RPA can be motivated from the perspective of the two-point density–density (dd) response moments of Eq. (8), which are central quantities to obtain in this approach to GW theory. We find that we can reformulate RPA entirely in terms of these dd-moments of the system and a strict recursive form for their inter-relation.⁷⁷ This recursive relation between the moments is a direct result of the fact that the RPA can be written as a quadratic Hamiltonian in Bosonic operators.^5,63,78–80 This effectively ensures that all information required to build the 2-point RPA dd-response is contained in the first two spectral moments, analogous to how all the information on the density of states in mean-field theory (quadratic in Fermionic operators) is contained in the first two Green’s function moments (i.e., the one-body density matrix and Fock matrix).

We start from the Casida formulation of RPA,^81,82 as a generalized eigenvalue decomposition

[\begin{matrix} A & B \\ - B & - A \end{matrix}] [\begin{matrix} X & Y \\ Y & X \end{matrix}] = [\begin{matrix} X & Y \\ Y & X \end{matrix}] [\begin{matrix} Ω & 0 \\ 0 & - Ω \end{matrix}],

(30)

where the left and right eigenvectors form the biorthogonal set as

{[\begin{matrix} X & Y \\ Y & X \end{matrix}]}^{T} [\begin{matrix} X & Y \\ - Y & - X \end{matrix}] = [\begin{matrix} I & 0 \\ 0 & - I \end{matrix}] .

(31)

This biorthogonality ensures an inverse relationship between (X + Y) and (X − Y)^T, as

(X + Y) {(X - Y)}^{T} = {(X + Y)}^{T} (X - Y) = I .

(32)

The A and B matrices are defined as

A_{i a, j b} = (ϵ_{a} - ϵ_{i}) δ_{i j} δ_{a b} + K_{i a, b j},

(33)

B_{i a, j b} = K_{i a, j b} .

(34)

Here,

K

is an interaction kernel that couples particle–hole excitations and de-excitations. In the traditional RPA (without second-order exchange), this coupling is taken to be the same for excitations, de-excitations, and their coupling, given by the static, bare Coulomb interaction,

K_{i a, j b} = (i a | j b) = K_{i a, b j}

⁠. Hole and particle orbital energies are given by ϵ_i and ϵ_a, respectively, defining the irreducible polarizability of the system from the reference state in A. Upon diagonalization, the eigenvectors defined by X_ia,ν and Y_ia,ν define the coefficients of the RPA excitations in the particle–hole and hole–particle basis, with energies Ω_ν, with Ω therefore a diagonal matrix of the positive (neutral) RPA excitation energies.

These neutral excitations define the poles of the full RPA reducible density–density (dd) response function, which can be constructed as

χ (ω) = [\begin{matrix} X & Y \\ Y & X \end{matrix}] {[\begin{matrix} ω I - Ω & 0 \\ 0 & - ω I - Ω \end{matrix}]}^{- 1} {[\begin{matrix} X & Y \\ Y & X \end{matrix}]}^{T} .

(35)

Note that this matrix is formally equivalent to the alternative directly dynamical construction from

χ (ω) = {(P {(ω)}^{- 1} - K)}^{- 1}

⁠, where P(ω) is the irreducible polarizability of the reference state. Considering the positive-frequency part of the dd-response [noting that the negative frequency part is symmetric due to the bosonic-like symmetry of Eq. ], we can write a more compact form of the dd-response as

η (ω) = (X + Y) {(ω I - Ω)}^{- 1} {(X + Y)}^{T},

(36)

which sums contributions from particle–hole and hole–particle fluctuations together, and from which optical properties, such as dynamic polarizabilities, can be computed.⁸³ However, in this work, we are interested in the order-by-order moments of the spectral distribution of Eq. over all RPA excitation energies, which is given as

η_{i a, j b}^{(n)} = - \frac{1}{π} \int_{0}^{\infty} I m [η_{i a, j b} (ω)] ω^{n} d ω .

(37)

The non-negative integer index n denotes the order of this static dd spectral moment information. Performing this integration results in the direct construction of the dd-moments as defined in Eq. , which can be written more compactly in matrix form as

η^{(n)} = (X + Y) Ω^{n} {(X + Y)}^{T},

(38)

and constitutes a central object of interest in this work, required for the GW self-energy moment construction of Eqs. and .⁸⁴

We now show that the RPA can be entirely reformulated in terms of the dd-moments [Eq. ] without loss of information, and expose constraints on the relationship between the moments of different order at the RPA level. First, we note that from the definition of the original eigendecomposition of Eq. , along with insertion of a resolution of the identity via Eq. , we find a relation between the first two dd-moments as

η^{(1)} = (X + Y) Ω {(X + Y)}^{T} = (A - B),

(39)

= η^{(0)} (A + B) η^{(0)},

(40)

noting that (A + B) = (X − Y)Ω(X − Y)^T. By taking the sum and difference of the two Casida equations of Eq. , we also find

(A + B) (X + Y) = (X - Y) Ω,

(41)

(A - B) (X - Y) = (X + Y) Ω,

(42)

from which an equation for the square of the RPA excitations can be found as

(A - B) (A + B) (X + Y) = (X + Y) Ω^{2},

(43)

which has previously appeared in the RPA literature.⁸⁵ By right-multiplying by (X + Y)^T, this leads to a relation between the zeroth and second dd-moment, as

(A - B) (A + B) η^{(0)} = η^{(2)} .

(44)

The above equations can be further generalized as a recursive relation to generate all higher-order moments from lower-order ones, as

η^{(m)} = (A - B) (A + B) η^{(m - 2)}

(45)

= {[η^{(0)} (A + B)]}^{m} η^{(0)} .

(46)

While these are important relations in themselves, they also illustrate that all the RPA excitations and weights in the dd-response of Eq. are implicitly accessible without requiring the explicit solution to the Casida equation (X, Y, and Ω matrices). This reformulation just requires knowledge of η⁽⁰⁾ as the central variable [which can be defined independently of the original equations via Eq. ], the A and B matrices defining the system and their interactions, and the recursive relations of Eq. . As an aside, the Tamm–Dancoff approximation (TDA) sets B = 0, which dramatically simplifies the resulting expressions due to the lack of correlation in the ground state, with η⁽⁰⁾ = I, and η⁽ⁿ⁾ = Aⁿ. This reflects that the 2-RDM in the TDA is equivalent to that of mean-field theory. Finally, we note in passing that the ground state correlation energy from the RPA can be similarly formulated in terms of the zeroth-order dd-moment, as

E_{corr}^{RPA} = \frac{1}{2} Tr [η^{(0)} (A + B) - A] .

(47)

Related expressions for the RPA correlation energy can be found in Eqs. (39) and (40) of Ref. 86 in terms of the quantities η⁽⁰⁾, A + B and A − B (where there, these quantities are referred to as

Q_{I}^{dRPA}

⁠, ɛ, and ɛ + 2K). Equivalence between these expressions (as well as the more commonly used expression found in Ref. 87) can be seen by noting [using Eqs. , , and ] that

Tr [η^{(0)} (A + B)]

(48)

= Tr [{(η^{(0)})}^{- 1} (A - B)]

(49)

= Tr [{({(A - B)}^{\frac{1}{2}} (A + B) {(A - B)}^{\frac{1}{2}})}^{\frac{1}{2}}] = Tr [Ω] .

(50)

Overall, this perspective on the RPA in terms of dd-moments is key to open new avenues, such as the ones explored in this work.

IV. EFFICIENT EVALUATION OF SELF-ENERGY AND DENSITY RESPONSE MOMENTS

Given the recasting of the RPA dd-response in Sec. III in terms of its lowest order moment [Eq. (40)] and recursion to access the higher moments via Eq. (46), we now consider the efficient $O [N^{4}]$ evaluation of these quantities that are central to the approach in this work, thus avoiding their formal $O [N^{6}]$ construction via Eq. (38). The derivation here is heavily inspired by the seminal resolution of the identity (RI) approach of Furche to compute RPA correlation energies,⁵⁶ with key adaptations to target these dd-moments to arbitrary order, rather than the correlation energy.

We first employ a standard low-rank decomposition of the two-electron repulsion integrals (e.g., via density fitting or Cholesky decomposition) as

(i a | j b) ≃ \sum_{P} V_{i a, P} V_{j b, P} = V V^{T},

(51)

where we use P, Q, … to index elements of this auxiliary (RI) basis, whose dimension N_aux scales

O [N]

with system size. We define an intermediate quantity

{\tilde{η}}_{i a, P}^{(n)} = \sum_{j b} η_{i a, j b}^{(n)} V_{j b, P} .

(52)

If this intermediate can be efficiently found, then the greater self-energy moment of Eq. can be rewritten as

Σ_{p q}^{(n, >)} = \sum_{t = 0}^{n} (\binom{n}{t}) (ϵ_{c}^{n - t} V_{p c, Q} (V_{q c, P} ({\tilde{η}}_{i a, P}^{(t)} V_{i a, Q}))),

(53)

where the brackets indicate the order of contractions in order to preserve

O [N^{4}]

scaling, and Einstein summation is implied. The lesser self-energy moment of Eq. can be recast in an analogous fashion.⁸⁸

Obtaining all dd-moments of the form of Eq. up to order n can be simply reduced to knowledge of the first two moments

{\tilde{η}}_{i a, P}^{(0)}

and

{\tilde{η}}_{i a, P}^{(1)}

⁠, via use of the recursive relationship between the moments as given by Eq. , as for even moment orders,

{\tilde{η}}^{(n)} = {[(A - B) (A + B)]}^{n / 2} {\tilde{η}}^{(0)},

(54)

and for odd moment orders,

{\tilde{η}}^{(n)} = {[(A - B) (A + B)]}^{(n - 1) / 2} {\tilde{η}}^{(1)},

(55)

where we have omitted explicit indices for brevity. Ensuring that an

O [N^{4}]

scaling is retained in this recursion relies on (A − B)(A + B) admitting a form where it can be written as a diagonal plus low-rank correction. For the RPA, this is true since [from Eqs. and ]

{(A - B)}_{i a, j b} = (ϵ_{a} - ϵ_{i}) δ_{i j} δ_{a b} = D

(56)

is a purely diagonal matrix, while using Eq. we can cast (A + B) into an appropriate form as

(A + B) = D + 2 K = D + 2 \sum_{P} V_{i a, P} V_{j b, P} .

(57)

We therefore express the low-rank asymmetrically decomposed form of (A − B)(A + B) in a general fashion as a diagonal plus asymmetric low-rank part, as

(A - B) (A + B) = D^{2} + S_{L} S_{R}^{T},

(58)

Where, for the RPA, D is defined by Eq. , S_L = DV and S_R = 2V. Future work will explore other analogous approaches where (A − B)(A + B) can be decomposed in this way, for applicability to, e.g., the Bethe–Salpeter equation or other RPA variants with (screened) exchange contributions.^62,89,90 From this low-rank decomposition and the recursive definition of Eqs. and , a fixed number of dd-moments of the form of Eq. can be found in

O [N^{4}]

time, provided the original

{\tilde{η}}^{(0)}

and

{\tilde{η}}^{(1)}

values are known.

We now consider how to obtain these initial low-order dd-moments efficiently. From the definitions of Eqs. and , and specifying the standard RPA definition of Eq. , we find that it is straightforward to efficiently construct the first moment, as

{\tilde{η}}^{(1)} = D V = S_{L} .

(59)

The zeroth-order dd-moment can be constructed via a rapidly convergent numerical integration, as we will show. From Eq. , we can write

(A - B) (A + B) = {(η^{(0)} (A + B))}^{2},

(60)

from which we can find an expression for η⁽⁰⁾ as

η^{(0)} = {[(A - B) (A + B)]}^{\frac{1}{2}} {(A + B)}^{- 1} .

(61)

We note that, for RPA to be well-defined with positive excitation energies, (A − B) and (A + B) must both be positive-definite matrices.⁸⁵ Using Eqs. and , we can write the low-rank RPA form of this as

{\tilde{η}}^{(0)} = {(D^{2} + S_{L} S_{R}^{T})}^{\frac{1}{2}} {(D + 2 V V^{T})}^{- 1} V .

(62)

We first consider the evaluation of the second half of this expression, which we denote as T. We can use the Woodbury matrix identity to rewrite it as

T = {(D + 2 V V^{T})}^{- 1} V

(63)

= D^{- 1} V - 2 D^{- 1} V {(I + 2 V^{T} D^{- 1} V)}^{- 1} V^{T} D^{- 1} V .

(64)

This now only requires the inversion of the diagonal matrix, D, and a matrix of dimension N_aux, with the overall ov × N_aux matrix able to be constructed in

O [N_{aux}^{3} + N_{aux}^{2} o v]

time.

Having constructed T, we can complete the evaluation of

{\tilde{η}}^{(0)}

using the definition of the matrix square-root as an integration in the complex plane,⁹¹

M^{\frac{1}{2}} = \frac{1}{π} \int_{- \infty}^{\infty} (I - z^{2} {(M + z^{2} I)}^{- 1}) d z .

(65)

From Eq. , this results in

{\tilde{η}}^{(0)} = {(D^{2} + S_{L} S_{R}^{T})}^{\frac{1}{2}} T

(66)

= \frac{1}{π} \int_{- \infty}^{\infty} (I - z^{2} {(D^{2} + S_{L} S_{R}^{T} + z^{2} I)}^{- 1}) T d z .

(67)

We can modify this integrand into one more efficient for numerical integration, via another application of the Woodbury matrix identity to reduce the scaling of the matrix inverse. We also simplify the notation by introducing the intermediates,

F (z) = {(D^{2} + z^{2} I)}^{- 1},

(68)

Q (z) = S_{R}^{T} F (z) S_{L},

(69)

where F(z) is a diagonal matrix in the ov space, and Q(z) is a N_aux × N_aux matrix that can be constructed in

O [N_{aux}^{2} o v]

time. This casts Eq. into the form

{\tilde{η}}^{(0)} = \frac{1}{π} \int_{- \infty}^{\infty} [I - z^{2} F (z) (I - S_{L} {(I + Q (z))}^{- 1} S_{R}^{T} F (z))] T d z .

(70)

For each value of the integration variable z, the integrand is a matrix of size ov × N_aux, which can be constructed in

O [N^{4}]

scaling, rendering it efficient for numerical quadrature. Along with the results of Eqs. , , , and , this therefore completes the ambition of constructing a fixed number of dd-response moments needed for the moment-truncated GW method as defined in Eq. , in no more than

O [N^{4}]

scaling (and

O [N^{3}]

memory).

However, manipulations of the resulting integrand and choice of quadrature points can further improve the efficiency of their construction by ensuring a faster decay of the integrand and separating components that can be analytically integrated. This derivation is given in Appendix A, and results in a final

O [N_{aux}^{2} o v + N_{aux}^{3}]

expression to evaluate for the zeroth-order dd-moment as

\begin{align} {\tilde{η}}^{(0)} & = D T + \frac{1}{π} \int_{0}^{\infty} e^{- t D} S_{L} S_{R}^{T} e^{- t D} T d t \\ + \frac{1}{π} \int_{- \infty}^{\infty} z^{2} F (z) S_{L} ({(I + Q (z))}^{- 1} - I) S_{R}^{T} F (z) T d z . \end{align}

(71)

The first numerical integral in Eq. (where the integrand decays exponentially) is computed via Gauss–Laguerre quadrature, while the second (where the integrand decays as

O [z^{- 4}]

⁠) is evaluated via Clenshaw–Curtis quadrature. A comparison of the decay of the original and refined integrands is shown in the inset to Fig. 1.

FIG. 1.

Exponential convergence of the numerical integration (NI) error for η(0)̃ in Eq. (71) with respect to integration points, for the singlet oxygen dimer in a cc-pVTZ basis at equilibrium (1.207 Å) bond length. Also included are the two error estimates of the true NI error, used to check convergence and estimate the required number of points [see Appendix B, Eq. (B9) for the “nested fit” and Eq. (B4) for the “lower bound” definitions]. Inset: The originally derived integrand [Eq. (70)], and the form optimized for efficient NI given in Eq. (71) and derived in Appendix A, showing the faster decay.

View large Download slide

Exponential convergence of the numerical integration (NI) error for $\tilde{η^{(0)}}$ in Eq. (71) with respect to integration points, for the singlet oxygen dimer in a cc-pVTZ basis at equilibrium (1.207 Å) bond length. Also included are the two error estimates of the true NI error, used to check convergence and estimate the required number of points [see Appendix B, Eq. (B9) for the “nested fit” and Eq. (B4) for the “lower bound” definitions]. Inset: The originally derived integrand [Eq. (70)], and the form optimized for efficient NI given in Eq. (71) and derived in Appendix A, showing the faster decay.

The scaling of the grid spacing of both numerical quadratures is optimized to ensure exact integration of the trace of a diagonal approximation for the integrand, analogous to the grid optimization discussed in Ref. 56. For numerical robustness, we optimize the quadrature for evaluating

Tr {[(A - B) (A + B)]}^{\frac{1}{2}}

⁠, rather than the full integrand. We write a diagonal approximation to (A − B)(A + B) as

M^{D} = D^{2} + S_{L}^{D} {(S_{R}^{D})}^{T}

(72)

with

{(S_{L}^{D})}_{i a, R} = {(S_{L})}_{i a, Q} {(S_{R})}_{i a, Q} δ_{R, i a},

(73)

{(S_{R}^{D})}_{j b, R} = δ_{R, j b} .

(74)

This is the same form as Eq. , but contains only diagonal matrices (denoted by subscript “D” labels) and has an auxiliary space of size ov. As such, the exact square root and all quantities within the numerical integrations can be obtained in

O [o v]

computational time. We then seek to ensure the trace of the difference between the exact and numerically integrated estimate of

{M^{D}}^{\frac{1}{2}}

vanishes. This is achieved for both integrals in Eq. , with the analytic form for the first integral given by

\frac{1}{2 π} diag (D^{- 1} S_{L} S_{R}^{T}) = I^{offset}

and the second numerical integral as

{M^{D}}^{\frac{1}{2}} - D - I^{offset} = I^{int}

⁠.

Writing this explicitly, given quadrature points and weights {z_i, w_i} for a n_p-point infinite or semi-infinite quadrature, we seek to scale our points by a factor a, which is a root the objective functions

O^{offset} (a) = I^{offset} - \frac{1}{π} \sum_{i}^{n_{p}} a w_{i} Tr (e^{- 2 D a z_{i}} S_{L} S_{R}^{T}),

(75)

O^{int} (a) = I^{int} - \frac{1}{π} \sum_{i}^{n_{p}} a w_{i} Tr (I^{D} (a z_{i})),

(76)

I^{D} (z) = z^{2} F^{D} (z) S_{L}^{D} ({(I + Q^{D} (z))}^{- 1} - I) {(S_{R}^{D})}^{T} F^{D} (z),

(77)

where Eq. is minimized to optimize the grid for the first integral of Eq. , and Eq. minimized for the second integral. This can be done via either simple root-finding or minimization and gives a robust optimization of the integration grids in

O [o v]

computational cost. The resulting exponential convergence of the zeroth dd-moment estimate with number of quadrature points, along with the error estimates derived in Appendix B, are shown in Fig. 1 for both numerical integrands. We find that as few as 12 quadrature points are sufficient for high accuracy in the results of this work, while the number of points is expected to increase for systems with a small or vanishing spectral gap.

V. REDUCTION TO CUBIC-SCALING GW

With this reformulation of GW in terms of the moments of the self-energy, it is possible to further reduce the scaling to cubic in time and quadratic in memory with respect to system size, in common with the lowest-scaling GW approaches.^24,42,92,93 We stress that this is not an asymptotic scaling after exploiting screening and locality arguments, but rather a formal scaling exploiting further rank reduction of quantities. To do this, we employ a double factorization of the Coulomb integrals, allowing them to be written as a product of five rank-2 tensors, as

(i a | j b) ≃ \sum_{P, Q} X_{i P} X_{a P} Z_{P Q} X_{j Q} X_{b Q} .

(78)

This factorizes the orbital product into separate terms, and is also known as tensor hypercontraction or CP decomposition, used for various recent low-scaling formulations of quantum chemical methods, where the dimension of the P and Q space rigorously grow linearly with system size.⁹⁴ The use of this doubly factorized form has also been previously suggested in the use of a reduced-scaling RPA and particle–particle RPA schemes.^95,96 This form for the integrals can be directly constructed with controllable errors in

O [N^{3}]

time.^97,98

Once found, the Z_PQ can be symmetrically decomposed as Z_PQ = Y_PRY_QR. By replacing the density-fitted integral tensor V_iaR in the above expressions with the fully factorized form X_iPX_aPY_RP, the contractions to form the moments of the GW self-energy, and hence the Green’s function and quasi-particle spectrum naturally follow as $O [N^{3}]$ with formation of appropriate intermediates. We also require the numerical computation of the partially transformed dd-moment, $Z_{Q S} X_{a S} X_{i S} η_{i a, j b}^{(0)} X_{j R} X_{b R} Z_{P R}$ ⁠. Inspired by the low-scaling approach taken to RPA correlation energies in the work of Ref. 99, this can be achieved with the use of a contracted double-Laplace transform in the place of the original numerical integration procedure. This factorizes the squared energy denominator F(z), allowing the occupied and unoccupied indices to be contracted independently, similar to the space–time approaches to GW.⁹³ While this becomes a two-dimensional numerical integral, optimal quadrature grids can be calculated in a minimax sense.^100,101

Applied to Eq. , this contracted double-Laplace transform takes the form

F_{i a, i a} (z) = D_{i a, i a} \int_{p = 0}^{\infty} \frac{\sin z p}{z} e^{- p D_{i a, i a}} d p,

(79)

which allows the key matrix Q(z) of Eq. to be obtained as

Q_{P S} (z) = 2 \int_{p = 0}^{\infty} \frac{\sin z p}{z} Y_{P Q} (A_{Q R} (p) - B_{Q R} (p)) Y_{S R} d p,

(80)

where both intermediates

A_{Q R} (p) = X_{i Q} X_{a Q} ϵ_{a} e^{- ϵ_{a} p} e^{ϵ_{i} p} X_{i R} X_{a R},

(81)

B_{Q R} (p) = X_{i Q} X_{a Q} ϵ_{i} e^{- ϵ_{a} p} e^{ϵ_{i} p} X_{i R} X_{a R}

(82)

can be evaluated in

O [N^{3}]

cost. Further contractions in the evaluation of Eq. also follow naturally with cubic scaling.

An alternative approach to reduce the scaling (to asymptotically linear) without requiring the doubly factorized integrals is to screen the atomic orbital density-fitted integral contributions (constructed with the overlap metric), along with the double-Laplace transform, exploiting locality as has been recently performed for the RPA correlation energy.^99,101 An explicit numerical demonstration of this reduction to cubic cost via the double factorization of the Coulomb tensor of Eq. (78) will follow in a forthcoming work, with numerical results in the rest of this work employing the quartic scaling algorithm described in Sec. IV.

VI. RESULTS

A. Comparison to quasiparticle GW approaches

We first consider the convergence of the moment-truncated G₀W₀ algorithm compared to more traditional implementations, as found in the PySCF software package.^102–105 As found to be more effective for molecular systems due to the importance of exact exchange, we perform all calculations on a restricted Hartree–Fock.¹¹ In Fig. 2, we consider the convergence of the first ionization potential (IP) of singlet O₂ with conserved self-energy moment order. We compare this to two G₀W₀ implementations, one of which performs an exact frequency integration (denoted “full QP-G₀W₀,” which scales as $O [N^{6}]$ ⁠), and one which performs an analytic continuation (AC) of the imaginary frequency self-energy to real-frequencies via a fit to Padé approximants in order to perform the convolution (denoted “AC QP-G₀W₀,” which scales as $O [N^{4}]$ ⁠).^44,105,106 However, both of these two implementations also solve the diagonal approximation to the quasiparticle equation in solving for each state, effectively imposing a diagonal approximation to the self-energy in the MO basis. This is avoided in our work; however, we can constrain a similar diagonal approximation by simply removing the off-diagonal components of our computed self-energy moments. This does not result in a significant computational saving in our approach and, therefore, is only relevant for comparison purposes when considering the effect of this neglected off-diagonal part of the self-energy.

FIG. 2.

Convergence of the first G0W0 IP of singlet O2 in a cc-pVDZ basis with respect to the number of conserved self-energy moments. Also shown is the same quasiparticle state computed from traditional G0W0 via an exact frequency integration (“full”) and analytic continuation approach to G0W0 (“AC”), albeit both imposing a diagonal approximation in their solution to the QP equation. In the moment expansion, we also consider a diagonal approximation, by explicitly removing the non-diagonal parts of our computed self-energy moments (“diagonal Σ moment”). All approaches find an IP of 8.49 eV to within 10 meV, with the difference between full-frequency and AC approaches 5 meV, and the relaxation of diagonal approximation also accounting for small ∼5 meV differences, with the reference Hartree–Fock IP for comparison being 9.73 eV.

View large Download slide

Convergence of the first G₀W₀ IP of singlet O₂ in a cc-pVDZ basis with respect to the number of conserved self-energy moments. Also shown is the same quasiparticle state computed from traditional G₀W₀ via an exact frequency integration (“full”) and analytic continuation approach to G₀W₀ (“AC”), albeit both imposing a diagonal approximation in their solution to the QP equation. In the moment expansion, we also consider a diagonal approximation, by explicitly removing the non-diagonal parts of our computed self-energy moments (“diagonal Σ moment”). All approaches find an IP of 8.49 eV to within 10 meV, with the difference between full-frequency and AC approaches 5 meV, and the relaxation of diagonal approximation also accounting for small $\sim 5$ meV differences, with the reference Hartree–Fock IP for comparison being 9.73 eV.

As can be seen in Fig. 2, the IP converges rapidly with moment order, with the full self-energy moments converging slightly faster and more systematically without the diagonal approximation (something also observed in other applications). The diagonal approximation converges to the “exact” frequency integration as expected, with our more complete (non-diagonal) self-energy moment approach only very slightly different, indicating the relative unimportance of the non-diagonal self-energy components in this system and lack of significant correlation-induced coupling between the mean-field quasiparticle states. Furthermore, the analytic continuation approach is also highly accurate for this system, introducing an error of only 5 meV compared to the full frequency integration. We have furthermore numerically verified that our approach scales as $O [N^{4}]$ ⁠, with computational cost comparable to the analytic continuation approach.

Having demonstrated correctness compared to a high-scaling exact frequency implementation, we can compare our results to AC-G₀W₀ across a larger test set to consider the moment truncation convergence. We use the established “GW100” benchmark test set, where many GW (and other excited state method) implementations have been rigorously compared.^{12,47,107,108} This benchmark set contains 102 diverse systems, with the IP of the molecules in the set ranging from $\sim 4 - 25$ eV, and featuring molecules with bound metal atoms (including metal diatomics and clusters), strongly ionic bonding, and molecules with a strongly delocalized electronic structure. The molecules range in size from simple atomic systems to the five canonical nucleobases and large aliphatic and aromatic hydrocarbons, providing a suitable range in system size.

In Fig. 3, we consider the discrepancy in the first IP, electron affinity (EA), and quasiparticle gap across all systems as the order of conserved self-energy moments increases in a realistic def2-TZVPP basis, again compared to AC-G₀W₀. Errors in individual systems, along with the mean signed error (MSE) across the set (white circle), is shown for each conserved moment order. This MSE for the first IP decreases from −0.142 eV for the lowest order moment conservation, to −11 meV when up to the 11th-order moment is conserved, with the EA errors generally a little larger. Similarly, the gap calculations converge to a MSE of −34.8 meV, with standard deviation about the AC-G₀W₀ result of only 91 meV across all systems. We note that there may also be small differences arising due to the comparison with AC-G₀W₀ due to the approximations of the frequency integration via analytically continued quantities, as well as the differences in whether off-diagonal parts of the self-energy are included. These will contribute to the discrepancy between the methods at each order; however, the comparison is not strictly equivalent, and therefore these errors will be overestimated compared to an exact frequency and non-diagonal G₀W₀ limit, the general trend, convergence, and level of accuracy, which can be reached with moment order is likely to be similar.

FIG. 3.

Errors in the IP, EA, and gap for each system in the GW100 benchmark set, compared to AC-G0W0 in a def2-TZVPP basis, for each order of self-energy moment conservation. White circles show the mean signed error (MSE) aggregated over the test set for the given moment truncation in each quantity.

View large Download slide

Errors in the IP, EA, and gap for each system in the GW100 benchmark set, compared to AC-G₀W₀ in a def2-TZVPP basis, for each order of self-energy moment conservation. White circles show the mean signed error (MSE) aggregated over the test set for the given moment truncation in each quantity.

It is important to put the scale of the moment truncation convergence in the context of the overall accuracy of the G₀W₀ method for these systems. In Fig. 4, we therefore compare the aggregated mean absolute error (MAE) in the moment-truncated G₀W₀ IP and EA values over this GW100 test set to highly accurate coupled-cluster with singles, doubles and perturbative triple excitation [CCSD(T)] calculations on the separate charged and neutral systems, which is often used as a more faithful benchmark to compare against than experiment. We can therefore see the convergence of the moment truncation to the AC-G₀W₀ values compared to the intrinsic error in the method. This intrinsic error is found to dominate over the error due to the self-energy moment truncation for higher numbers of conserved moments. We note, however, that the mean error compared to CCSD decreases systematically with increasing numbers of moments for these systems, which contrasts with observed behavior for moment-truncated GF2 theory (where lower order self-energy truncations were found to give rise to a more accurate overall excitations).⁴⁶ It is natural to, therefore, also consider whether a simple extrapolation can improve the moment-truncated results to an infinite moment limit. We therefore apply a linear extrapolation of the excitation energies to the infinite moment limit from the two most complete moment calculations of each system. We can see from Fig. 4 that these results continue the trend of the MAE across the test set, slightly overshooting the AC-G₀W₀ comparison, albeit noting the other potential sources of discrepancy between these values discussed earlier.

FIG. 4.

Mean absolute errors (eV) for the IP (x-axis) and EA (y-axis) across the GW100 benchmark set in a def2-TZVPP basis compared to accurate CCSD(T) values. AC-G0W0 values, as well as moment-truncated results are shown, with the number in each data point marker giving the number of conserved moments. Extrapolation of individual data points from the 9th- and 11th-order conserved moment values are also provided, denoted by the truncation label “∞.”

View large Download slide

Mean absolute errors (eV) for the IP (x-axis) and EA (y-axis) across the GW100 benchmark set in a def2-TZVPP basis compared to accurate CCSD(T) values. AC-G₀W₀ values, as well as moment-truncated results are shown, with the number in each data point marker giving the number of conserved moments. Extrapolation of individual data points from the 9th- and 11th-order conserved moment values are also provided, denoted by the truncation label “∞.”

B. Full frequency spectra

One of the strengths of the moment-conserving approach to GW of this work is the ability to obtain all excitations from a given order of truncation in a single complete diagonalization of the effective Hamiltonian of Eq. (27). This allows full frequency spectra to be obtained, with the approximation not expected to bias significantly toward accuracy in any particular energy range, making it suitable for G₀W₀ excitations beyond frontier excitations. The description of low-lying states is a particular challenge for many other GW approaches, with analytic continuation becoming less reliable, and alternatives like contour-deformation scaling as $O [N^{5}]$ ⁠, in general, to obtain the full spectrum.⁴¹ In Fig. 5, we therefore show a series of spectra plotting on the real frequency axis for the guanine nucleobase in a def2-TZVPP basis set over a 100 eV energy window about the quasiparticle gap, taken from the GW100 benchmark set, and compared to the AC-G₀W₀ full-frequency spectrum.

FIG. 5.

Spectral functions for guanine in a def2-TZVPP basis set, for HF, AC-G0W0, and the moment-conserving G0W0 approach with zero to five block Lanczos iterations, thereby conserving up to the 11th-order self-energy moment. The values indicated in the spectra indicate the Wasserstein metric taken with respect to the AC-G0W0 spectrum, quantifying the difference between the spectral distributions. The AC-G0W0 spectrum is indicated transparently behind the other spectra to ease visualization of the convergence, and was calculated using an iterative diagonal approximation to the quasiparticle equation.

View large Download slide

Spectral functions for guanine in a def2-TZVPP basis set, for HF, AC-G₀W₀, and the moment-conserving G₀W₀ approach with zero to five block Lanczos iterations, thereby conserving up to the 11th-order self-energy moment. The values indicated in the spectra indicate the Wasserstein metric taken with respect to the AC-G₀W₀ spectrum, quantifying the difference between the spectral distributions. The AC-G₀W₀ spectrum is indicated transparently behind the other spectra to ease visualization of the convergence, and was calculated using an iterative diagonal approximation to the quasiparticle equation.

The convergence of the spectrum is shown for a series of conserved G₀W₀ moment orders, from the HF level up to the AC-G₀W₀ spectrum. The AC-G₀W₀ spectrum is also shown “behind” the other spectra to allow the deviations to be observed for each moment. It can be seen that the full-frequency spectrum rapidly converges with conserved self-energy moment order, even for high-energy states where the HF approximation is poor. The similarity of each spectrum to the AC-G₀W₀ result across the full frequency range can be rigorously quantified via the Wasserstein or “earth-mover” metric, which describes the similarity between probability distributions. This metric is shown as the value inside of each plot, indicating a rapid and robust convergence of the spectral features from the mean-field to the full G₀W₀ spectrum with increasing numbers of included moments.

This Wasserstein metric convergence plateaus at the $\sim 7 th$ -order conserved moments, with further orders not improving this metric further. This could be due to numerical precision of the algorithm or fundamental approximations in the AC-G₀W₀, such as the precision of the analytic continuation, or the diagonal approximation to the self-energy. Furthermore, it should be noted that the moment-conserving GW approximation will rigorously have a larger number of poles included in its spectrum compared to those GW approaches that rely on an iterative solution to the QP equation, which considers the change to a single MO at a time. These additional peaks are likely very low weighted for this weakly correlated system, yet could be contributing to this discrepancy with AC-G₀W₀ in the Wasserstein metric. We consider this point in more detail in Sec. VI C.

C. Multiple solutions and additional spectral features

This fact that many low-scaling GW implementations rely on an iterative solution to solve the quasiparticle equation can be a source of error and loss of robustness. This is because when self-energy poles are found near quasiparticle energies, the GW poles can split into multiple peaks, where the final excitation energy converged to can depend sensitively on the specifics of the root-finding algorithm used to solve the QP equation. This was highlighted in the Ref. 12 as a significant source of error, where a number of simple systems were found to exhibit a number of poles close to the HOMO energy level (with these solutions spanning a range of up to 1 eV). The specifics of which pole is converged to (with undesired solutions called “spurious”) then depended on initial conditions, choices of optimization method, and specifics of the self-energy construction, or linearization of the QP equation. The requirement to select one of these multiple solutions to the QP equation can also manifest in undesirable discontinuities in the excitation energies as, e.g., the molecular geometry or correlation strength changes, as described in Refs. 109–111. An indicator for the presence of these “spurious” poles and multiple solutions is the magnitude of the quasi-particle weight or renormalization factor evaluated at the quasi-particle energy, defined as $Z_{p} = {(1 - \frac{\partial Σ_{p p} (ϵ_{p})}{\partial ω})}^{- 1}$ ⁠, which indicates the approximate weight in the quasi-particle solutions.

Since the moment-conserving GW approach obtains all poles in one step (including satellites and low-weighted features), all the excitations can be characterized by their quasiparticle weight, and either selected as specific excitation energies, or all excitations included in the spectrum to ensure smooth changes with molecular geometry. Points of discontinuity will therefore manifest as the presence of multiple lower-weighted solutions at a given energy, giving a smooth change to a broadened feature in the spectrum near self-energy poles. To demonstrate this, we consider the same simple system as Ref. 111, observing the GW quasiparticle energies as a function of the inter-nuclear distance in the H₂ dimer in a 6-31G basis set. Figure 6 shows quasiparticle energies, self-energies, and quasiparticle renormalization factors for AC-G₀W₀, and the moment conserving G₀W₀ approach with both up to the 1st- and 11th-order conserved self-energy moments in each sector. The self-energies plotted are the diagonal elements corresponding to the particular MO, evaluated at the respective quasiparticle energy.

FIG. 6.

Quasiparticle energies, self-energies, and renormalization factors for the H2 dimer in a 6-31G basis set with varying bond lengths. Shown are results for AC-G0W0 and moment-conserved G0W0 with zero and five iterations, thereby conserving up to the 1st- and 11th-order moments, respectively. The self-energy corresponds to the diagonal element corresponding to the particular orbital, evaluated at the quasiparticle energy. The renormalization factor corresponds to Zp=(1−∂Σpp(ϵp)∂ω)−1, with larger values indicating a more quasiparticle-like excitation. The transparent lines in the lower panel show the existence of multiple solutions, broadening spectral features near the self-energy poles, and where a single dominant solution can be chosen (defined here by maximum overlap with the corresponding MO) denoted by the thicker line.

View large Download slide

Quasiparticle energies, self-energies, and renormalization factors for the H₂ dimer in a 6-31G basis set with varying bond lengths. Shown are results for AC-G₀W₀ and moment-conserved G₀W₀ with zero and five iterations, thereby conserving up to the 1st- and 11th-order moments, respectively. The self-energy corresponds to the diagonal element corresponding to the particular orbital, evaluated at the quasiparticle energy. The renormalization factor corresponds to $Z_{p} = {(1 - \frac{\partial Σ_{p p} (ϵ_{p})}{\partial ω})}^{- 1}$ ⁠, with larger values indicating a more quasiparticle-like excitation. The transparent lines in the lower panel show the existence of multiple solutions, broadening spectral features near the self-energy poles, and where a single dominant solution can be chosen (defined here by maximum overlap with the corresponding MO) denoted by the thicker line.

The figure shows the HOMO and first three unoccupied states found in this system, with the AC-G₀W₀ (first row) exhibiting discontinuities in the LUMO+2 state at slightly compressed geometries, and discontinuities in the LUMO+1 state at slightly stretched geometries. As discussed, these changing solutions arise from the specifics of the root-finding in the solution to the QP equation, which will generally (but not always) converge to the solution with largest quasiparticle weight between the multiple options, indicated by the renormalization factor. These discontinuous changes between states are also shown to coincide with poles in the self-energy in the second column at the MO energies for a given separation. These AC-G₀W₀ results are essentially the same as those found in Ref. 111.

When only the first two self-energy moments are conserved in each sector (second row), the approximation to the self-energy renders its pole structure sufficiently sparse such that their poles are pushed far from the MO energies at all geometries for these states. While this regularization removes the discontinuities, this significant approximation renders the renormalization factors close to one at all points, indicating only small changes from the original MOs. The final row represents a G₀W₀ calculation with up to the 11th-order conserved moments. With this finer resolution of the self-energy dynamics, the structure of the self-energy closely matches the one from AC-G₀W₀; however, the multiple solutions are all found simultaneously, with their changing quasiparticle weights shown. The points of discontinuity are replaced by the presence of multiple solutions at similar energies and with competing quasiparticle weights, providing broad spectral features at those points.

If a single solution is required, the specific excitation can be selected from the manifold based rigorously on, e.g., the maximum overlap with the MO of interest (shown as thicker lines in the plot) or largest quasiparticle weight, all of which can easily be selected. This removes the uncertainty in the energies of the states based on the unphysical specifics of the QP solution algorithm, without incurring additional complexity or cost in the moment-conserving algorithm. Furthermore, the relaxation of the diagonal approximation of the self-energy in this approach is expected to be more significant at these points of multiple solution, where mixing between the different MOs is expected to be more pronounced.

VII. CONCLUSIONS AND OUTLOOK

In this work, we present a reformulation of the GW theory of quasiparticle excitations, based around a systematic expansion and conservation of the spectral moments of the self-energy. This contrasts with other approaches designed to approximate the central frequency integration of GW theory, which use, e.g., grid expansions, analytic continuation, or contour deformation in order to affect a scaling reduction from the exact theory. The moment expansion presented in this work has appealing features arising from the avoidance of an iterative solution to the quasiparticle equation for each state (avoiding “spurious” solutions), diagonal self-energy approximations, or requirements for analytic continuation of dynamical quantities. It allows for all excitation energies and weights to be obtained directly in a non-iterative single diagonalization of a small effective Hamiltonian, controlled by a single parameter governing the number of conserved self-energy moments.

Full RPA screening and particle–hole coupling in the self-energy is included, which is captured with $O [N^{4}]$ computational scaling via a numerical one-dimensional integration, with a reduction to cubic scaling also proposed. This approach is enabled by a recasting of the RPA in terms of the moments of the density–density response function. Applied across GW100 test set, we find rapid convergence to established GW methodology results for both state-specific and full spectral properties, with errors due to the incompleteness of the moment expansion many times smaller than the inherent accuracy of the method. The formulation follows relatively closely from previous low-scaling approaches to RPA correlation energies, enabling these codes to be adapted to low-scaling GW methods with relatively little effort.

Going forward, we will aim to test the limits of the moment-truncated GW formulation, pushing it to larger systems, including the solid-state and different reference states, lower-scaling variants, and the inclusion of various self-consistent flavors of the theory (beyond the G₀W₀ implementation here). The reformulation of RPA in terms of a recursive moment expansion also lends itself to a low-scaling implementation of the Bethe–Salpeter equation for neutral excitations, which we will explore in the future, as well as other beyond-RPA approaches. Finally, we will also explore the connections of this moment expansion to kernel polynomial approaches, which expand spectral quantities in terms of Chebyshev and other orthogonal polynomial expansions.¹¹²

ACKNOWLEDGMENTS

The authors thank Filipp Furche and Johannes Tölle for useful feedback on the manuscript. G.H.B. gratefully acknowledges support from the Royal Society via a University Research Fellowship, as well as funding from the European Union’s Horizon 2020 Research and Innovation Program under Grant Agreement No. 759063. We are grateful to the UK Materials and Molecular Modeling Hub for computational resources, which is partially funded by EPSRC (Grant Nos. EP/P020194/1 and EP/T022213/1). Further computational resources were awarded under the embedded CSE program of the ARCHER2 UK National Supercomputing Service (http://www.archer2.ac.uk).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Charles J. C. Scott: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Oliver J. Backhouse: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). George H. Booth: Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

Open-source code for reproduction of all results in this paper, along with examples, can be found at https://github.com/BoothGroup/momentGW. The repository also includes the data used in this paper relating to the GW100 benchmark.

APPENDIX A: IMPROVED NUMERICAL QUADRATURE

While the integral expression for ${\tilde{η}}^{(0)}$ derived in Eq. (70) in the main text is sufficient for numerical quadrature in $O [N^{4}]$ scaling, it can be refined by deducting large and/or slowly decaying contributions that can be efficiently integrated separately, in order to minimize the number of quadrature points that are required to obtain a given overall accuracy. This again follows many of the developments in the numerical integration (NI)-RPA approaches for the correlation energy in the literature due to the commonalities in their forms;⁵⁶ however, there are also important differences.

We can first consider a mean-field contribution to the integral component of

{\tilde{η}}^{(0)}

⁠, i.e., setting the Coulomb interaction to zero in the expression for

{[(A - B) (A + B)]}^{\frac{1}{2}}

⁠. This just leaves a contribution from the irreducible polarizability, which can be analytically integrated as

D T = {(D^{2})}^{\frac{1}{2}} T = \frac{1}{π} \int_{- \infty}^{\infty} (I - z^{2} {(D^{2} + z^{2} I)}^{- 1}) T d z,

(A1)

where we have used the integral form of the matrix square root again [Eq. ]. We can add and subtract these different forms within Eq. to obtain

{\tilde{η}}^{(0)} = D T + \frac{1}{π} \int_{- \infty}^{\infty} z^{2} F (z) S_{L} {(I + Q (z))}^{- 1} S_{R}^{T} F (z) T d z,

(A2)

which substantially reduces the magnitude of the numerically integrated component of

{\tilde{η}}^{(0)}

⁠. However, we can further improve on this, by increasing the rate of decay of the remaining integrand with respect to z. It can be seen from Eqs. and that both F(z) and Q(z) decay as z⁻² in the large-z limit. Series expanding [I + Q(z)]⁻¹ ∼ I − Q(z) + Q²(z) − …, we find that the leading order decay of the integrand in Eq. is

O [z^{- 2}]

⁠. We consider this contribution in isolation, noting that if it can be removed from the numerical integration, the next leading order will result in the integrand decaying at an improved

O [z^{- 4}]

rate.

This leading-order contribution to the integral of Eq. can be written as

C_{2} = \frac{1}{π} \int_{- \infty}^{\infty} z^{2} F (z) S_{L} S_{R}^{T} F (z) T d z .

(A3)

This can be analytically integrated via the residue theorem to give

C_{2} = \frac{1}{π} ((S_{L} S_{R}^{T}) ◦ E) T,

(A4)

where ◦ denotes the Hadamard or element-wise product, and we define E as

E_{i a, j b} = {(D_{i a, i a} + D_{j b, j b})}^{- 1} .

(A5)

Writing Eq. with explicit indices for clarity it can be seen as a second-order direct-MP2-like contribution of

{(C_{2})}_{i a, P} = \frac{1}{π} \sum_{j b} \frac{\sum_{Q} {(S_{L})}_{i a, Q} {(S_{R})}_{j b, Q}}{ϵ_{a} + ϵ_{b} - ϵ_{i} - ϵ_{j}} T_{j b, P} .

(A6)

However, in this form, it is unfortunately unable to be efficiently evaluated in

O [N^{4}]

time (instead scaling as

O [N^{5}]

⁠). We note that a similar subtraction of this leading-order term in the computation of the RPA correlation energy is in contrast efficiently computable, due to the presence of a cyclically invariant trace operation.⁵⁶

We therefore take a different approach, by noting that we can form an algebraic Sylvester equation, as

((S_{L} S_{R}^{T}) ◦ E) D + D ((S_{L} S_{R}^{T}) ◦ E) = S_{L} S_{R}^{T},

(A7)

which can be verified by substitution of Eq. . Writing the matrix of interest we wish to efficiently compute as

N = (S_{L} S_{R}^{T}) ◦ E

⁠, this Sylvester equation has an integral solution of the form¹¹³

N = \int_{0}^{\infty} e^{- t D} S_{L} S_{R}^{T} e^{- t D} d t .

(A8)

This can be verified via substitution into Eq. , as

D N + N D = \int_{0}^{\infty} (D e^{- t D} S_{L} S_{R}^{T} e^{- t D} + e^{- t D} S_{L} S_{R}^{T} e^{- t D} D) d t

(A9)

= - \int_{0}^{\infty} \frac{d}{d t} (e^{- t D} S_{L} S_{R}^{T} e^{- t D}) d t

(A10)

= - {[e^{- t D} S_{L} S_{R}^{T} e^{- t D}]}_{0}^{\infty}

(A11)

= S_{L} S_{R}^{T} .

(A12)

The numerical integration of Eq. can therefore be used to compute the leading order term in the numerical integration of Eq. . We note that this is essentially analogous to a matrix-generalization of the Laplace transform to the direct second-order diagram, which in frequency space involves energy denominators of the type given in Eq. .¹¹⁴ It may appear as though we have swapped one numerical integral for another; however, the transformation to the form given in Eq. gives an integrand, which is exponentially rather than quadratically decaying, leading to a form that can be very effectively integrated via Gauss–Laguerre quadrature with exponential convergence and only a small number of samples, in practice.

The final expression for the numerical integration of the RPA zeroth-order dd-response moment is

\begin{align} {\tilde{η}}^{(0)} & = D T + \frac{1}{π} \int_{0}^{\infty} e^{- t D} S_{L} S_{R}^{T} e^{- t D} T d t \\ + \frac{1}{π} \int_{- \infty}^{\infty} z^{2} F (z) S_{L} ({(I + Q (z))}^{- 1} - I) S_{R}^{T} F (z) T d z . \end{align}

(A13)

The integrand of the first integral decays exponentially and is efficiently computed with Gauss–Laguerre quadrature, with the second integral (decaying as

O [z^{- 4}]

⁠) is evaluated with Clenshaw–Curtis quadrature.

APPENDIX B: NUMERICAL INTEGRATION ERROR ESTIMATES

We make use of two separate approaches to estimate the error in ${\tilde{η}}^{(0)}$ evaluated through numerical integration, which we can use as a check for convergence of this key intermediate. These arise from quite different considerations and provide complementary estimates. In the following, we will denote the moment estimate resulting from a numerical integration with n_p points as ${\tilde{η}}_{n_{p}}^{(0)} = {\tilde{η}}_{\infty}^{(0)} + Δ {\tilde{η}}_{n_{p}}^{(0)}$ ⁠.

For our first estimate, we begin by noting from Eq. that the exact

{\tilde{η}}^{(0)}

will satisfy the relation

{({\tilde{η}}_{\infty}^{(0)})}^{T} (A + B) {\tilde{η}}_{\infty}^{(0)} = V^{T} (A - B) V .

(B1)

We can then take the deviation from this relation at finite integration samples as an error estimate, as X(n_p). Expanding the estimated quantity in terms of the exact result and error gives

X (n_{p}) = {({\tilde{η}}_{n_{p}}^{(0)})}^{T} (A + B) {\tilde{η}}_{n_{p}}^{(0)} - V^{T} (A - B) V

(B2)

\begin{align} = {(Δ {\tilde{η}}_{n_{p}}^{(0)})}^{T} (A + B) {\tilde{η}}_{n_{p}}^{(0)} + {({\tilde{η}}_{n_{p}}^{(0)})}^{T} (A + B) Δ {\tilde{η}}_{n_{p}}^{(0)} \\ - {(Δ {\tilde{η}}_{n_{p}}^{(0)})}^{T} (A + B) Δ {\tilde{η}}_{n_{p}}^{(0)} . \end{align}

(B3)

We can then explicitly evaluate ‖X(n_p)‖₂ in

O [N^{4}]

time under the same conditions as the rest of this manuscript and apply the properties of the Frobenius norm to write

‖ X (n_{p}) ‖_{2} \leq 2 ‖ Δ {\tilde{η}}_{n_{p}}^{(0)} ‖_{2} ‖ (A + B) {\tilde{η}}_{n_{p}}^{(0)} ‖_{2} - ‖ Δ {\tilde{η}}_{n_{p}}^{(0)} ‖_{2}^{2} ‖ (A + B) ‖_{2} .

Writing

‖ Δ {\tilde{η}}_{n_{p}}^{(0)} ‖_{2} = x

⁠, the error is therefore rigorously bound by the solutions satisfying the quadratic inequality,

x^{2} - 2 \frac{‖ (A + B) {\tilde{η}}_{n_{p}}^{(0)} ‖_{2}}{‖ (A + B) ‖_{2}} x + \frac{‖ X (n_{p}) ‖_{2}}{‖ (A + B) ‖_{2}} \leq 0 .

(B4)

In practice, while we find that the two roots to this equation rigorously bound the error in the quadrature, their values are not a tight bound on the true error (in particular, the upper bound is excessively large to be used as an estimate, though the lower bound on the error that it provides is reasonable, shown in Fig. 1 of the main text as “lower bound”). Combined with the

O [N^{4}]

scaling of evaluating X(n_p) this makes this estimate of limited practical use, and we seek an improved estimate that provides a tighter upper bound on the true error.

To do this we take advantage of the nested nature of the Clenshaw–Curtis quadrature within the integral providing the dominant error contribution. Using this, we can obtain estimates at both one-half and one quarter of the current number of integration points with no additional cost. We show how this can give an estimate of the error in

{\tilde{η}}_{n_{p}}^{(0)} (A + B)

⁠, which we can relate to the error in

{\tilde{η}}_{n_{p}}^{(0)}

⁠. Parameterizing the exponential convergence of the L₂ norm error as

‖ Δ {\tilde{η}}_{n_{p}}^{(0)} ‖_{2} = α e^{- β n_{p}}

⁠, we can rigorously relate the differences between estimates at different levels of sampling, as

Δ_{n_{1}, n_{2}} = ‖ {\tilde{η}}_{n_{1}}^{(0)} - {\tilde{η}}_{n_{2}}^{(0)} ‖_{2}

(B5)

\leq ‖ Δ {\tilde{η}}_{n_{1}}^{(0)} ‖_{2} + ‖ Δ {\tilde{η}}_{n_{2}}^{(0)} ‖_{2}

(B6)

\approx α (e^{- β n_{1}} + e^{- β n_{2}}) .

(B7)

We seek to estimate the value of

‖ Δ {\tilde{η}}_{4 n_{p}}^{(0)} ‖_{2}

given the values

Δ_{4 n_{p}, 2 n_{p}}

and

Δ_{4 n_{p}, n_{p}}

⁠, which are computed as intermediates in the nested quadrature of 4n_p points.

Setting

x = e^{- β n_{p}}

⁠, using the approximation defined in Eq. for

Δ_{4 n_{p}, 2 n_{p}} = α (x^{4} + x^{2})

and

Δ_{4 n_{p}, n_{p}} = α (x^{4} + x)

⁠, and eliminating the factor of α, we find

(Δ_{4 n_{p}, n_{p}} - Δ_{4 n_{p}, 2 n_{p}}) x^{3} + Δ_{4 n_{p}, n_{p}} x - Δ_{4 n_{p}, 2 n_{p}} = 0 .

(B8)

The smallest positive, real solution to this is used to estimate the error (given as αx⁴) as

‖ Δ {\tilde{η}}_{4 n_{p}}^{(0)} ‖_{2} \approx Δ_{4 n_{p}, 2 n_{p}} {(1 + x^{- 2})}^{- 1} .

(B9)

This estimator provides a consistent and systematic overestimate of the error, allowing us to use this as a reliable upper bound to the true error of the quadrature. This property, combined with only requiring the quantities

Δ_{4 n_{p}, 2 n_{p}}

and

Δ_{4 n_{p}, n_{p}}

⁠, which are unavoidable intermediates of the nested integration, makes this essentially a computationally free estimate of the error for a given number of grid points and, therefore, as a reliable check for convergence. If the resulting error estimate is too large, the number of grid points can be increased by a factor of two, and all points already evaluated can be reused in the next estimate. The convergence of both error estimates can be seen in Fig. 1 of the main text, with the result of Eq. denoted “nested fit,” while the lower bound from the solution of Eq. denoted “lower bound.”

REFERENCES

1.

A. J.

Cohen

,

P.

Mori-Sánchez

, and

W.

Yang

, “

Challenges for density functional theory

,”

Chem. Rev.

112

,

289

–

320

(

2012

).

https://doi.org/10.1021/cr200107z

Google Scholar

Crossref

PubMed

2.

G.

Onida

,

L.

Reining

, and

A.

Rubio

, “

Electronic excitations: Density-functional versus many-body Green’s-function approaches

,”

Rev. Mod. Phys.

74

,

601

–

659

(

2002

).

https://doi.org/10.1103/revmodphys.74.601

Google Scholar

Crossref

3.

M. S.

Hybertsen

and

S. G.

Louie

, “

First-principles theory of quasiparticles: Calculation of band gaps in semiconductors and insulators

,”

Phys. Rev. Lett.

55

,

1418

–

1421

(

1985

).

https://doi.org/10.1103/physrevlett.55.1418

Google Scholar

Crossref

PubMed

4.

F.

Aryasetiawan

and

O.

Gunnarsson

, “

The GW method

,”

Rep. Prog. Phys.

61

,

237

–

312

(

1998

).

https://doi.org/10.1088/0034-4885/61/3/002

Google Scholar

Crossref

5.

L.

Hedin

, “

On correlation effects in electron spectroscopies and the GW approximation

,”

J. Phys.: Condens. Matter

11

,

R489

–

R528

(

1999

).

https://doi.org/10.1088/0953-8984/11/42/201

Google Scholar

Crossref

6.

W. G.

Aulbur

,

L.

Jönsson

, and

J. W.

Wilkins

, “

Quasiparticle calculations in solids

,” in

Solid State Physics

(

Elsevier

,

2000

), pp.

1

–

218

.

Google Scholar

Crossref

7.

C.

Friedrich

and

A.

Schindlmayr

, “

Many-body perturbation theory: The GW approximation

,” in

Computational Nanoscience: Do It Yourself!

, edited by J. Grotendorst, S. Bluegel, and D. Marx (NIC, 2006), NIC Series Vol. 31, pp. 335-355, ISBN: 3-00-017350-1.

8.

A.

Kutepov

,

S. Y.

Savrasov

, and

G.

Kotliar

, “

Ground-state properties of simple elements from GW calculations

,”

Phys. Rev. B

80

,

041103

(

2009

).

https://doi.org/10.1103/physrevb.80.041103

Google Scholar

Crossref

9.

S. H.

Ke

, “

All-electron GW methods implemented in molecular orbital space: Ionization energy and electron affinity of conjugated molecules

,”

Phys. Rev. B

84

,

205415

(

2011

).

https://doi.org/10.1103/physrevb.84.205415

Google Scholar

Crossref

10.

F.

Bruneval

, “

Ionization energy of atoms obtained from GW self-energy or from random phase approximation total energies

,”

J. Chem. Phys.

136

,

194107

(

2012

).

https://doi.org/10.1063/1.4718428

Google Scholar

Crossref

PubMed

11.

M. J.

van Setten

,

F.

Weigend

, and

F.

Evers

, “

The GW-method for quantum chemistry applications: Theory and implementation

,”

J. Chem. Theory Comput.

9

,

232

–

246

(

2013

).

https://doi.org/10.1021/ct300648t

Google Scholar

Crossref

PubMed

12.

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

).

https://doi.org/10.1021/acs.jctc.5b00453

Google Scholar

Crossref

PubMed

13.

L.

Reining

, “

The GW approximation: Content, successes and limitations

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

8

,

e1344

(

2017

).

https://doi.org/10.1002/wcms.1344

Google Scholar

Crossref

14.

D.

Golze

,

M.

Dvorak

, and

P.

Rinke

, “

The GW compendium: A practical guide to theoretical photoemission spectroscopy

,”

Front. Chem.

7

,

377

(

2019

).

https://doi.org/10.3389/fchem.2019.00377

Google Scholar

Crossref

PubMed

15.

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

).

https://doi.org/10.1103/physrev.139.a796

Google Scholar

Crossref

16.

L.

Hedin

and

S.

Lundqvist

, “

Effects of electron-electron and electron-phonon interactions on the one-electron states of solids

,” in

Solid State Physics

(

Elsevier

,

1970

), pp.

1

–

181

.

Google Scholar

Crossref

17.

D. C.

Langreth

and

J. P.

Perdew

, “

Exchange-correlation energy of a metallic surface: Wave-vector analysis

,”

Phys. Rev. B

15

,

2884

–

2901

(

1977

).

https://doi.org/10.1103/physrevb.15.2884

Google Scholar

Crossref

18.

D. P.

Chong

,

Recent Advances in Density Functional Methods

(

World Scientific

,

1995

).

Google Scholar

Crossref

19.

A.

Heßelmann

and

A.

Görling

, “

Random phase approximation correlation energies with exact Kohn–Sham exchange

,”

Mol. Phys.

108

,

359

–

372

(

2010

).

https://doi.org/10.1080/00268970903476662

Google Scholar

Crossref

20.

X.

Ren

,

P.

Rinke

,

C.

Joas

, and

M.

Scheffler

, “

Random-phase approximation and its applications in computational chemistry and materials science

,”

J. Mater. Sci.

47

,

7447

–

7471

(

2012

).

https://doi.org/10.1007/s10853-012-6570-4

Google Scholar

Crossref

21.

A.

Stan

,

N. E.

Dahlen

, and

R.

van Leeuwen

, “

Levels of self-consistency in the GW approximation

,”

J. Chem. Phys.

130

,

114105

(

2009

).

https://doi.org/10.1063/1.3089567

Google Scholar

Crossref

PubMed

22.

D.

Foerster

,

P.

Koval

, and

D.

Sánchez-Portal

, “

An O(N³) implementation of Hedin’s GW approximation for molecules

,”

J. Chem. Phys.

135

,

074105

–

074119

(

2011

).

https://doi.org/10.1063/1.3624731

Google Scholar

Crossref

PubMed

23.

F.

Bruneval

and

M. A. L.

Marques

, “

Benchmarking the starting points of the GW approximation for molecules

,”

J. Chem. Theory Comput.

9

,

324

–

329

(

2013

).

https://doi.org/10.1021/ct300835h

Google Scholar

Crossref

PubMed

24.

P.

Liu

,

M.

Kaltak

,

J.

Klimeš

, and

G.

Kresse

, “

Cubic scaling GW: Towards fast quasiparticle calculations

,”

Phys. Rev. B

94

,

165109

(

2016

).

https://doi.org/10.1103/physrevb.94.165109

Google Scholar

Crossref

25.

J. W.

Knight

,

X.

Wang

,

L.

Gallandi

,

O.

Dolgounitcheva

,

X.

Ren

,

J. V.

Ortiz

,

P.

Rinke

,

T.

Körzdörfer

, and

N.

Marom

, “

Accurate ionization potentials and electron affinities of acceptor molecules III: A benchmark of GW methods

,”

J. Chem. Theory Comput.

12

,

615

–

626

(

2016

).

https://doi.org/10.1021/acs.jctc.5b00871

Google Scholar

Crossref

PubMed

26.

E.

Maggio

,

P.

Liu

,

M. J.

van Setten

, and

G.

Kresse

, “

GW₁00: A plane wave perspective for small molecules

,”

J. Chem. Theory Comput.

13

,

635

–

648

(

2017

).

https://doi.org/10.1021/acs.jctc.6b01150

Google Scholar

Crossref

PubMed

27.

F.

Bruneval

,

N.

Dattani

, and

M. J.

van Setten

, “

The GW miracle in many-body perturbation theory for the ionization potential of molecules

,”

Front. Chem.

9

,

749779

(

2021

).

https://doi.org/10.3389/fchem.2021.749779

Google Scholar

Crossref

PubMed

28.

U.

von Barth

and

B.

Holm

, “

Self-consistent GW₀ results for the electron gas: Fixed screened potential W₀ within the random-phase approximation

,”

Phys. Rev. B

54

,

8411

–

8419

(

1996

).

https://doi.org/10.1103/physrevb.54.8411

Google Scholar

Crossref

29.

B.

Holm

and

U.

von Barth

, “

Fully self-consistent GW self-energy of the electron gas

,”

Phys. Rev. B

57

,

2108

–

2117

(

1998

).

https://doi.org/10.1103/physrevb.57.2108

Google Scholar

Crossref

30.

W.-D.

Schöne

and

A. G.

Eguiluz

, “

Self-consistent calculations of quasiparticle states in metals and semiconductors

,”

Phys. Rev. Lett.

81

,

1662

–

1665

(

1998

).

https://doi.org/10.1103/physrevlett.81.1662

Google Scholar

Crossref

31.

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

).

https://doi.org/10.1103/physrevb.63.075112

Google Scholar

Crossref

32.

S. V.

Faleev

,

M.

van Schilfgaarde

, and

T.

Kotani

, “

All-electron self-consistent GW approximation: Application to Si, MnO, and NiO

,”

Phys. Rev. Lett.

93

,

126406

–

126415

(

2004

).

https://doi.org/10.1103/PhysRevLett.93.126406

Google Scholar

Crossref

PubMed

33.

M.

van Schilfgaarde

,

T.

Kotani

, and

S.

Faleev

, “

Quasiparticle self-consistent GW theory

,”

Phys. Rev. Lett.

96

,

226402

(

2006

).

https://doi.org/10.1103/PhysRevLett.96.226402

Google Scholar

Crossref

PubMed

34.

A.

Stan

,

N. E.

Dahlen

, and

R. v.

Leeuwen

, “

Fully self-consistent GW calculations for atoms and molecules

,”

Europhys. Lett.

76

,

298

–

304

(

2006

).

https://doi.org/10.1209/epl/i2006-10266-6

Google Scholar

Crossref

35.

T.

Kotani

,

M.

van Schilfgaarde

, and

S. V.

Faleev

, “

Quasiparticle self-consistent GW method: A basis for the independent-particle approximation

,”

Phys. Rev. B

76

,

165106

(

2007

).

https://doi.org/10.1103/physrevb.76.165106

Google Scholar

Crossref

36.

M.

Shishkin

and

G.

Kresse

, “

Self-consistent GW calculations for semiconductors and insulators

,”

Phys. Rev. B

75

,

235102

(

2007

).

https://doi.org/10.1103/physrevb.75.235102

Google Scholar

Crossref

37.

F.

Caruso

,

P.

Rinke

,

X.

Ren

,

M.

Scheffler

, and

A.

Rubio

, “

Unified description of ground and excited states of finite systems: The self-consistent GW approach

,”

Phys. Rev. B

86

,

81102

(

2012

).

https://doi.org/10.1103/physrevb.86.081102

Google Scholar

Crossref

38.

F.

Bruneval

,

T.

Rangel

,

S. M.

Hamed

,

M.

Shao

,

C.

Yang

, and

J. B.

Neaton

, “

MOLGW 1: Many-body perturbation theory software for atoms, molecules, and clusters

,”

Comput. Phys. Commun.

208

,

149

–

161

(

2016

).

https://doi.org/10.1016/j.cpc.2016.06.019

Google Scholar

Crossref

39.

F.

Kaplan

,

M. E.

Harding

,

C.

Seiler

,

F.

Weigend

,

F.

Evers

, and

M. J.

van Setten

, “

Quasi-particle self-consistent GW for molecules

,”

J. Chem. Theory Comput.

12

,

2528

–

2541

(

2016

).

https://doi.org/10.1021/acs.jctc.5b01238

Google Scholar

Crossref

PubMed

40.

Y.

Jin

,

N. Q.

Su

, and

W.

Yang

, “

Renormalized singles Green’s function for quasi-particle calculations beyond the G0W0 approximation

,”

J. Phys. Chem. Lett.

10

,

447

–

452

(

2019

).

https://doi.org/10.1021/acs.jpclett.8b03337

Google Scholar

Crossref

PubMed

41.

I.

Duchemin

and

X.

Blase

, “

Robust analytic-continuation approach to many-body GW calculations

,”

J. Chem. Theory Comput.

16

,

1742

–

1756

(

2020

).

https://doi.org/10.1021/acs.jctc.9b01235

Google Scholar

Crossref

PubMed

42.

I.

Duchemin

and

X.

Blase

, “

Cubic-scaling all-electron GW calculations with a separable density-fitting space–time approach

,”

J. Chem. Theory Comput.

17

,

2383

–

2393

(

2021

).

https://doi.org/10.1021/acs.jctc.1c00101

Google Scholar

Crossref

PubMed

43.

C.-N.

Yeh

,

S.

Iskakov

,

D.

Zgid

, and

E.

Gull

, “

Fully self-consistent finite-temperature GW in Gaussian Bloch orbitals for solids

,”

Phys. Rev. B

106

,

235104

(

2022

).

https://doi.org/10.1103/physrevb.106.235104

Google Scholar

Crossref

44.

X.

Ren

,

P.

Rinke

,

V.

Blum

,

J.

Wieferink

,

A.

Tkatchenko

,

A.

Sanfilippo

,

K.

Reuter

, and

M.

Scheffler

, “

Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions

,”

New J. Phys.

14

,

053020

(

2012

).

https://doi.org/10.1088/1367-2630/14/5/053020

Google Scholar

Crossref

45.

G. E.

Engel

,

B.

Farid

,

C. M. M.

Nex

, and

N. H.

March

, “

Calculation of the GW self-energy in semiconducting crystals

,”

Phys. Rev. B

44

,

13356

–

13373

(

1991

).

https://doi.org/10.1103/physrevb.44.13356

Google Scholar

Crossref

46.

O. J.

Backhouse

and

G. H.

Booth

, “

Efficient excitations and spectra within a perturbative renormalization approach

,”

J. Chem. Theory Comput.

16

,

6294

–

6304

(

2020

).

https://doi.org/10.1021/acs.jctc.0c00701

Google Scholar

Crossref

PubMed

47.

O. J.

Backhouse

,

A.

Santana-Bonilla

, and

G. H.

Booth

, “

Scalable and predictive spectra of correlated molecules with moment truncated iterated perturbation theory

,”

J. Phys. Chem. Lett.

12

,

7650

–

7658

(

2021

).

https://doi.org/10.1021/acs.jpclett.1c02383

Google Scholar

Crossref

PubMed

48.

O. J.

Backhouse

and

G. H.

Booth

, “

Constructing “full-Frequency” spectra via moment constraints for coupled cluster Green’s functions

,”

J. Chem. Theory Comput.

18

,

6622

–

6636

(

2022

).

https://doi.org/10.1021/acs.jctc.2c00670

Google Scholar

Crossref

PubMed

49.

P. V.

Sriluckshmy

,

M.

Nusspickel

,

E.

Fertitta

, and

G. H.

Booth

, “

Fully algebraic and self-consistent effective dynamics in a static quantum embedding

,”

Phys. Rev. B

103

,

085131

(

2021

).

https://doi.org/10.1103/physrevb.103.085131

Google Scholar

Crossref

50.

R.

Haydock

, “

The recursive solution of the Schrödinger equation

,”

Comput. Phys. Commun.

20

,

11

–

16

(

1980

).

https://doi.org/10.1016/0010-4655(80)90101-0

Google Scholar

Crossref

51.

S.

Adachi

and

E.

Lipparini

, “

Sum rules in extended RPA theories

,”

Nucl. Phys. A

489

,

445

–

460

(

1988

).

https://doi.org/10.1016/0375-9474(88)90006-1

Google Scholar

Crossref

52.

D.

Karlsson

and

R.

van Leeuwen

, “

Partial self-consistency and analyticity in many-body perturbation theory: Particle number conservation and a generalized sum rule

,”

Phys. Rev. B

94

,

125124

(

2016

).

https://doi.org/10.1103/physrevb.94.125124

Google Scholar

Crossref

53.

J. R.

Sabin

,

J.

Oddershede

,

R.

Cabrera-Trujillo

,

S. P. A.

Sauer

,

E.

Deumens

, and

Y.

Öhrn

, “

Stopping power of molecules for fast ions

,”

Mol. Phys.

108

,

2891

–

2897

(

2010

).

https://doi.org/10.1080/00268976.2010.508753

Google Scholar

Crossref

54.

C.

Van Caillie

and

R. D.

Amos

, “

Static and dynamic polarisabilities, Cauchy coefficients and their anisotropies: An evaluation of DFT functionals

,”

Chem. Phys. Lett.

328

,

446

–

452

(

2000

).

https://doi.org/10.1016/s0009-2614(00)00942-8

Google Scholar

Crossref

55.

Y. N.

Kalugina

and

A. J.

Thakkar

, “

Ab initio calculations of static dipole polarizabilities and Cauchy moments for the halomethanes CH_mClnF₄₋m−n

,”

Chem. Phys. Lett.

644

,

20

–

24

(

2016

).

https://doi.org/10.1016/j.cplett.2015.11.044

Google Scholar

Crossref

56.

H.

Eshuis

,

J.

Yarkony

, and

F.

Furche

, “

Fast computation of molecular random phase approximation correlation energies using resolution of the identity and imaginary frequency integration

,”

J. Chem. Phys.

132

,

234114

(

2010

).

https://doi.org/10.1063/1.3442749

Google Scholar

Crossref

PubMed

57.

O. J.

Backhouse

,

M.

Nusspickel

, and

G. H.

Booth

, “

Wave function perspective and efficient truncation of renormalized second-order perturbation theory

,”

J. Chem. Theory Comput.

16

,

1090

–

1104

(

2020

).

https://doi.org/10.1021/acs.jctc.9b01182

Google Scholar

Crossref

PubMed

58.

E.

Rebolini

and

J.

Toulouse

, “

Range-separated time-dependent density-functional theory with a frequency-dependent second-order Bethe-Salpeter correlation kernel

,”

J. Chem. Phys.

144

,

094107

(

2016

).

https://doi.org/10.1063/1.4943003

Google Scholar

Crossref

PubMed

59.

M.

Véril

,

P.

Romaniello

,

J. A.

Berger

, and

P.-F.

Loos

, “

Unphysical discontinuities in GW methods

,”

J. Chem. Theory Comput.

14

,

5220

–

5228

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00745

Google Scholar

Crossref

PubMed

60.

R.

Quintero-Monsebaiz

,

E.

Monino

,

A.

Marie

, and

P.-F.

Loos

, “

Connections between many-body perturbation and coupled-cluster theories

,”

J. Chem. Phys.

157

,

231102

(

2022

).

https://doi.org/10.1063/5.0130837

Google Scholar

Crossref

PubMed

61.

S. J.

Bintrim

and

T. C.

Berkelbach

, “

Full-frequency GW without frequency

,”

J. Chem. Phys.

154

,

041101

(

2021

).

https://doi.org/10.1063/5.0035141

Google Scholar

Crossref

PubMed

62.

S. J.

Bintrim

and

T. C.

Berkelbach

, “

Full-frequency dynamical Bethe–Salpeter equation without frequency and a study of double excitations

,”

J. Chem. Phys.

156

,

044114

(

2022

).

https://doi.org/10.1063/5.0074434

Google Scholar

Crossref

PubMed

63.

J.

Tölle

and

G. K.-L.

Chan

, “

Exact relationships between the GW approximation and equation-of-motion coupled-cluster theories through the quasi-boson formalism

,” arXiv:2212.08982 (

2022

).

Google Scholar

64.

H. D.

Meyer

and

S.

Pal

, “

A band-Lanczos method for computing matrix elements of a resolvent

,”

J. Chem. Phys.

91

,

6195

–

6204

(

1989

).

https://doi.org/10.1063/1.457438

Google Scholar

Crossref

65.

H. G.

Weikert

,

H. D.

Meyer

,

L. S.

Cederbaum

, and

F.

Tarantelli

, “

Block Lanczos and many-body theory: Application to the one-particle Green’s function

,”

J. Chem. Phys.

104

,

7122

–

7138

(

1996

).

https://doi.org/10.1063/1.471429

Google Scholar

Crossref

66.

R.

Haydock

and

C. M. M.

Nex

, “

A general terminator for the recursion method

,”

J. Phys. C: Solid State Phys.

18

,

2235

(

1985

).

https://doi.org/10.1088/0022-3719/18/11/007

Google Scholar

Crossref

67.

G. E.

Engel

and

B.

Farid

, “

Calculation of the dielectric properties of semiconductors

,”

Phys. Rev. B

46

,

15812

–

15827

(

1992

).

https://doi.org/10.1103/physrevb.46.15812

Google Scholar

Crossref

68.

B.

Farid

, “

Dynamical correlation functions expressed in terms of many-particle ground-state wavefunction; the dynamical self-energy operator

,”

Philos. Mag. B

82

,

1413

–

1610

(

2002

).

https://doi.org/10.1080/13642810208222682

Google Scholar

Crossref

69.

T.

Fukaya

,

Y.

Nakatsukasa

,

Y.

Yanagisawa

, and

Y.

Yamamoto

, “

CholeskyQR2: A simple and communication-avoiding algorithm for computing a tall-skinny QR factorization on a large-scale parallel system

,” in

5th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems

(

IEEE

,

2014

), pp.

31

–

38

.

Google Scholar

Crossref

70.

T.

Fukaya

,

R.

Kannan

,

Y.

Nakatsukasa

,

Y.

Yamamoto

, and

Y.

Yanagisawa

, “

Shifted Cholesky QR for computing the QR factorization of ill-conditioned matrices

,”

SIAM J. Sci. Comput.

42

,

A477

–

A503

(

2020

).

https://doi.org/10.1137/18m1218212

Google Scholar

Crossref

71.

J.

Schirmer

,

A. B.

Trofimov

, and

G.

Stelter

, “

A non-Dyson third-order approximation scheme for the electron propagator

,”

J. Chem. Phys.

109

,

4734

–

4744

(

1998

).

https://doi.org/10.1063/1.477085

Google Scholar

Crossref

72.

J.

Schirmer

,

Many-Body Methods for Atoms, Molecules and Clusters

(

Springer Chem

,

Switzerland

,

2018

).

Google Scholar

Crossref

73.

T. C.

Berkelbach

, “

Communication: Random-phase approximation excitation energies from approximate equation-of-motion coupled-cluster doubles

,”

J. Chem. Phys.

149

,

041103

(

2018

).

https://doi.org/10.1063/1.5032314

Google Scholar

Crossref

PubMed

74.

V.

Rishi

,

A.

Perera

, and

R. J.

Bartlett

, “

A route to improving RPA excitation energies through its connection to equation-of-motion coupled cluster theory

,”

J. Chem. Phys.

153

,

234101

(

2020

).

https://doi.org/10.1063/5.0023862

Google Scholar

Crossref

PubMed

75.

A. L.

Dempwolff

,

M.

Schneider

,

M.

Hodecker

, and

A.

Dreuw

, “

Efficient implementation of the non-Dyson third-order algebraic diagrammatic construction approximation for the electron propagator for closed- and open-shell molecules

,”

J. Chem. Phys.

150

,

064108

(

2019

).

https://doi.org/10.1063/1.5081674

Google Scholar

Crossref

PubMed

76.

A. B.

Trofimov

and

J.

Schirmer

, “

Molecular ionization energies and ground- and ionic-state properties using a non-Dyson electron propagator approach

,”

J. Chem. Phys.

123

,

144115

(

2005

).

https://doi.org/10.1063/1.2047550

Google Scholar

Crossref

PubMed

77.

C. J. C.

Scott

and

G. H.

Booth

, “

Extending density matrix embedding: A static two-particle theory

,”

Phys. Rev. B

104

,

245114

(

2021

).

https://doi.org/10.1103/physrevb.104.245114

Google Scholar

Crossref

78.

D. C.

Langreth

, “

Singularities in the x-ray spectra of metals

,”

Phys. Rev. B

1

,

471

–

477

(

1970

).

https://doi.org/10.1103/physrevb.1.471

Google Scholar

Crossref

79.

P.

Nozières

and

C. T.

De Dominicis

, “

Singularities in the x-ray absorption and emission of metals. III. One-body theory exact solution

,”

Phys. Rev.

178

,

1097

–

1107

(

1969

).

https://doi.org/10.1103/physrev.178.1097

Google Scholar

Crossref

80.

B. I.

Lundqvist

and

V.

Samathiyakanit

, “

Single-particle spectrum of the degenerate electron gas IV. Ground state energy

,”

Phys. Kondens. Materie.

,

9

,

231

–

235

(

1969

).

https://doi.org/10.1007/bf02422566

Google Scholar

Crossref

81.

A.

Heßelmann

and

A.

Görling

, “

Random-phase approximation correlation methods for molecules and solids

,”

Mol. Phys.

109

,

2473

–

2500

(

2011

).

https://doi.org/10.1080/00268976.2011.614282

Google Scholar

Crossref

82.

G. P.

Chen

,

V. K.

Voora

,

M. M.

Agee

,

S. G.

Balasubramani

, and

F.

Furche

, “

Random-phase approximation methods

,”

Annu. Rev. Phys. Chem.

68

,

421

–

445

(

2017

).

https://doi.org/10.1146/annurev-physchem-040215-112308

Google Scholar

Crossref

PubMed

83.

W. A.

Parkinson

and

M. C.

Zerner

, “

The calculation of dynamic molecular polarizability

,”

J. Chem. Phys.

90

,

5606

–

5611

(

1989

).

https://doi.org/10.1063/1.456413

Google Scholar

Crossref

84.

We also note that a related expansion and recursive relation can also be constructed in terms of the inverse of these moments courtesy of the equivalence (X + Y)−1 = (X − Y)T, that is ${(η^{(n)})}^{- 1} = (X - Y) Ω^{- n} {(X - Y)}^{T}$ ⁠.

85.

F.

Furche

, “

On the density matrix based approach to time-dependent density functional response theory

,”

J. Chem. Phys.

114

,

5982

–

5992

(

2001

).

https://doi.org/10.1063/1.1353585

Google Scholar

Crossref

86.

J. G.

Ángyán

,

R.-F.

Liu

,

J.

Toulouse

, and

G.

Jansen

, “

Correlation energy expressions from the adiabatic-connection fluctuation–dissipation theorem approach

,”

J. Chem. Theory Comput.

7

,

3116

–

3130

(

2011

).

https://doi.org/10.1021/ct200501r

Google Scholar

Crossref

PubMed

87.

F.

Furche

, “

Developing the random phase approximation into a practical post-Kohn–Sham correlation model

,”

J. Chem. Phys.

129

,

114105

(

2008

).

https://doi.org/10.1063/1.2977789

Google Scholar

Crossref

PubMed

88.

We note that the derivation in this section does not rely on the V_jb,P tensor in Eq. 52 arising from this Coulomb form specifically, but rather that it involves a more general linear transformation from a space in the particle-hole product basis to a space which scales no more than $O [N]$ ⁠, (in this case the auxiliary basis).

89.

F.

Hummel

,

A.

Grüneis

,

G.

Kresse

, and

P.

Ziesche

, “

Screened exchange corrections to the random phase approximation from many-body perturbation theory

,”

J. Chem. Theory Comput.

15

,

3223

–

3236

(

2019

).

https://doi.org/10.1021/acs.jctc.8b01247

Google Scholar

Crossref

PubMed

90.

A.

Förster

, “

Assessment of the second-order statically screened exchange correction to the random phase approximation for correlation energies

,”

J. Chem. Theory Comput.

18

,

5948

–

5965

(

2022

).

https://doi.org/10.1021/acs.jctc.2c00366

Google Scholar

Crossref

PubMed

91.

N.

Hale

,

N. J.

Higham

, and

L. N.

Trefethen

, “

Computing A^α, log(A), and related matrix functions by contour integrals

,”

SIAM J. Numer. Anal.

46

,

2505

–

2523

(

2008

).

https://doi.org/10.1137/070700607

Google Scholar

Crossref

92.

J.

Wilhelm

,

P.

Seewald

,

M.

Del Ben

, and

J.

Hutter

, “

Large-scale cubic-scaling random phase approximation correlation energy calculations using a Gaussian basis

,”

J. Chem. Theory Comput.

12

,

5851

–

5859

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00840

Google Scholar

Crossref

PubMed

93.

A.

Förster

and

L.

Visscher

, “

Low-order scaling G0W0 by pair atomic density fitting

,”

J. Chem. Theory Comput.

16

,

7381

–

7399

(

2020

).

https://doi.org/10.1021/acs.jctc.0c00693

Google Scholar

Crossref

PubMed

94.

R. M.

Parrish

,

E. G.

Hohenstein

,

T. J.

Martínez

, and

C. D.

Sherrill

, “

Tensor hypercontraction. II. Least-squares renormalization

,”

J. Chem. Phys.

137

,

224106

(

2012

).

https://doi.org/10.1063/1.4768233

Google Scholar

Crossref

PubMed

95.

J.

Lu

and

K.

Thicke

, “

Cubic scaling algorithms for RPA correlation using interpolative separable density fitting

,”

J. Comput. Phys.

351

,

187

–

202

(

2017

).

https://doi.org/10.1016/j.jcp.2017.09.012

Google Scholar

Crossref

96.

N.

Shenvi

,

H.

van Aggelen

,

Y.

Yang

, and

W.

Yang

, “

Tensor hypercontracted ppRPA: Reducing the cost of the particle-particle random phase approximation from O(r⁶) to O(r⁴)

,”

J. Chem. Phys.

141

,

024119

(

2014

).

https://doi.org/10.1063/1.4886584

Google Scholar

Crossref

PubMed

97.

J.

Lu

and

L.

Ying

, “

Compression of the electron repulsion integral tensor in tensor hypercontraction format with cubic scaling cost

,”

J. Comput. Phys.

302

,

329

–

335

(

2015

).

https://doi.org/10.1016/j.jcp.2015.09.014

Google Scholar

Crossref

98.

K.

Pierce

and

E. F.

Valeev

, “

Efficient construction of canonical polyadic approximations of tensor networks

,”

J. Chem. Theory Comput.

19

,

71

–

81

(

2023

).

https://doi.org/10.1021/acs.jctc.2c00861

Google Scholar

Crossref

PubMed

99.

H. F.

Schurkus

and

C.

Ochsenfeld

, “

Communication: An effective linear-scaling atomic-orbital reformulation of the random-phase approximation using a contracted double-Laplace transformation

,”

J. Chem. Phys.

144

,

031101

(

2016

).

https://doi.org/10.1063/1.4939841

Google Scholar

Crossref

PubMed

100.

B.

Helmich-Paris

and

L.

Visscher

, “

Improvements on the minimax algorithm for the Laplace transformation of orbital energy denominators

,”

J. Comput. Phys.

321

,

927

–

931

(

2016

).

https://doi.org/10.1016/j.jcp.2016.06.011

Google Scholar

Crossref

101.

D.

Graf

,

M.

Beuerle

,

H. F.

Schurkus

,

A.

Luenser

,

G.

Savasci

, and

C.

Ochsenfeld

, “

Accurate and efficient parallel implementation of an effective linear-scaling direct random phase approximation method

,”

J. Chem. Theory Comput.

14

,

2505

–

2515

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00177

Google Scholar

Crossref

PubMed

102.

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

(

2017

).

https://doi.org/10.1002/wcms.1340

Google Scholar

Crossref

103.

Q.

Sun

,

X.

Zhang

,

S.

Banerjee

,

P.

Bao

,

M.

Barbry

,

N. S.

Blunt

,

N. A.

Bogdanov

,

G. H.

Booth

,

J.

Chen

,

Z.-H.

Cui

,

J. J.

Eriksen

,

Y.

Gao

,

S.

Guo

,

J.

Hermann

,

M. R.

Hermes

,

K.

Koh

,

P.

Koval

,

S.

Lehtola

,

Z.

Li

,

J.

Liu

,

N.

Mardirossian

,

J. D.

McClain

,

M.

Motta

,

B.

Mussard

,

H. Q.

Pham

,

A.

Pulkin

,

W.

Purwanto

,

P. J.

Robinson

,

E.

Ronca

,

E. R.

Sayfutyarova

,

M.

Scheurer

,

H. F.

Schurkus

,

J. E. T.

Smith

,

C.

Sun

,

S.-N.

Sun

,

S.

Upadhyay

,

L. K.

Wagner

,

X.

Wang

,

A.

White

,

J. D.

Whitfield

,

M. J.

Williamson

,

S.

Wouters

,

J.

Yang

,

J. M.

Yu

,

T.

Zhu

,

T. C.

Berkelbach

,

S.

Sharma

,

A. Y.

Sokolov

, and

G. K.-L.

Chan

, “

Recent developments in the PySCF program package

,”

J. Chem. Phys.

153

,

024109

(

2020

).

https://doi.org/10.1063/5.0006074

Google Scholar

Crossref

PubMed

104.

T.

Zhu

and

G. K.-L.

Chan

, “

Ab initio full cell GW + DMFT for correlated materials

,”

Phys. Rev. X

11

,

021006

(

2021

).

https://doi.org/10.1103/physrevx.11.021006

Google Scholar

Crossref

105.

T.

Zhu

and

G. K.-L.

Chan

, “

All-electron Gaussian-based G0W0 for valence and core excitation energies of periodic systems

,”

J. Chem. Theory Comput.

17

,

727

–

741

(

2021

).

https://doi.org/10.1021/acs.jctc.0c00704

Google Scholar

Crossref

PubMed

106.

J.

Wilhelm

,

M.

Del Ben

, and

J.

Hutter

, “

GW in the Gaussian and plane waves scheme with application to linear acenes

,”

J. Chem. Theory Comput.

12

,

3623

–

3635

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00380

Google Scholar

Crossref

PubMed

107.

F.

Caruso

,

M.

Dauth

,

M. J.

van Setten

, and

P.

Rinke

, “

Benchmark of GW approaches for the GW100 test set

,”

J. Chem. Theory Comput.

12

,

5076

–

5087

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00774

Google Scholar

Crossref

PubMed

108.

M. F.

Lange

and

T. C.

Berkelbach

, “

On the relation between equation-of-motion coupled-cluster theory and the GW approximation

,”

J. Chem. Theory Comput.

14

,

4224

–

4236

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00455

Google Scholar

Crossref

PubMed

109.

P.-F.

Loos

,

P.

Romaniello

, and

J. A.

Berger

, “

Green functions and self-consistency: Insights from the spherium model

,”

J. Chem. Theory Comput.

14

,

3071

–

3082

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00260

Google Scholar

Crossref

PubMed

110.

S.

Di Sabatino

,

P.-F.

Loos

, and

P.

Romaniello

, “

Scrutinizing GW-based methods using the Hubbard dimer

,”

Front. Chem.

9

,

751054

(

2021

).

https://doi.org/10.3389/fchem.2021.751054

Google Scholar

Crossref

PubMed

111.

E.

Monino

and

P.-F.

Loos

, “

Unphysical discontinuities, intruder states and regularization in GW methods

,”

J. Chem. Phys.

156

,

231101

(

2022

).

https://doi.org/10.1063/5.0089317

Google Scholar

Crossref

PubMed

112.

A.

Weiße

,

G.

Wellein

,

A.

Alvermann

, and

H.

Fehske

, “

The Kernel polynomial method

,”

Rev. Mod. Phys.

78

,

275

–

306

(

2006

).

https://doi.org/10.1103/revmodphys.78.275

Google Scholar

Crossref

113.

R.

Bhatia

and

P.

Rosenthal

, “

How and why to solve the operator equation AX − XB = Y

,”

Bull. London Math. Soc.

29

,

1

–

21

(

1997

).

https://doi.org/10.1112/s0024609396001828

Google Scholar

Crossref

114.

M.

Häser

and

J.

Almlöf

, “

Laplace transform techniques in Møller–Plesset perturbation theory

,”

J. Chem. Phys.

96

,

489

–

494

(

1992

).

https://doi.org/10.1063/1.462485

Google Scholar

Crossref

Published open access through an agreement with JISC Collections

A “moment-conserving” reformulation of GW theory

I. INTRODUCTION

II. MOMENT-TRUNCATED GW THEORY

A. Full GW spectrum from self-energy moments

III. DENSITY RESPONSE MOMENTS IN THE RPA

IV. EFFICIENT EVALUATION OF SELF-ENERGY AND DENSITY RESPONSE MOMENTS

V. REDUCTION TO CUBIC-SCALING GW

VI. RESULTS

A. Comparison to quasiparticle GW approaches

B. Full frequency spectra

C. Multiple solutions and additional spectral features

VII. CONCLUSIONS AND OUTLOOK

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: IMPROVED NUMERICAL QUADRATURE

APPENDIX B: NUMERICAL INTEGRATION ERROR ESTIMATES

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

A “moment-conserving” reformulation of GW theory

I. INTRODUCTION

II. MOMENT-TRUNCATED GW THEORY

A. Full GW spectrum from self-energy moments

III. DENSITY RESPONSE MOMENTS IN THE RPA

IV. EFFICIENT EVALUATION OF SELF-ENERGY AND DENSITY RESPONSE MOMENTS

V. REDUCTION TO CUBIC-SCALING GW

VI. RESULTS

A. Comparison to quasiparticle GW approaches

B. Full frequency spectra

C. Multiple solutions and additional spectral features

VII. CONCLUSIONS AND OUTLOOK

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: IMPROVED NUMERICAL QUADRATURE

APPENDIX B: NUMERICAL INTEGRATION ERROR ESTIMATES

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only