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[N4]$ approach to build and solve this Hamiltonian, with a proposal to reduce this further to $O[N3]$. 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 *GW*100 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.^{1} 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}(*ω*),^{21} (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[N6]$ 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[N3\u2212N4]$) 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[N4]$ time (with a proposed $O[N3]$ 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[N4]$ cost, based on ideas from the seminal work in 2010 by Furche and collaborators on low-scaling approaches for the RPA correlation energy.^{56} 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 *GW*100 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*_{0}*W*_{0} 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[N2]$ 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*_{0}*W*_{0}” 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

In these expressions, Ω_{ν} are the neutral excitation energies from RPA theory, and $Xia\nu $ and $Yia\nu $ 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. (1) therefore represents the “lesser” part, and the second term the “greater” part of the full *G*_{0}*W*_{0} self-energy.

Exact evaluation of the RPA excitations (Ω_{ν}) scales as $O[N6]$, 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 *n*th-order moments of the resulting dynamical distributions, as

and similarly

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. (1), the moments can be constructed as

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 *n*th-order spectral moments of the RPA density-response, summed over both particle–hole excitation and de-excitation components, as

Evaluating the self-energy spectral moments of Eqs. (9) and (10) 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, $\eta ia,jb(n)$. These relations are subsequently used in Sec. IV to demonstrate how the self-energy moments can be evaluated in $O[N4]$ 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

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.^{12} 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,

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

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*_{0}*W*_{0}, this auxiliary space (i.e., the dimension of **d**) must scale as $O[N3]$.^{60,61} However, in the moment truncation, $W\u0303$ and $d\u0303$ 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[nN]$, where *n* is the number of conserved self-energy moments. This allows for a complete diagonalization of $H\u0303$, obtaining all excitations in a single shot, and a reconstruction of the full *GW* Green’s function from its Lehmann representation in $O[(nN)3]$ computational effort, avoiding the need for any grids or iterative solutions once $H\u0303$ is found.

To find this effective upfolded representation of the moment-conserving dynamics, we modify the block Lanczos procedure to ensure the construction a $H\u0303$ of minimal size, whose effective hole and particle self-energy moments exactly match the ones from Eqs. (9) and (10). 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,^{63} via inspection from Eq. (1), with

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

where we define $q\u0303(j)$ as

The **q**^{(j)} are block Lanczos vectors of depth *j*, which form a recursive Krylov space as $q(j)=q1q2\cdots qj$, 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.^{45} 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**_{1} can be found via a QR factorization of the exact *GW* couplings **W**, as

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

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**_{1} = **W**^{†}**L**^{−1,†}. Subsequent block Lanczos vectors can then be defined according to the standard three-term recurrence

where the on-diagonal blocks are defined as

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. (18) permits the definition

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 **Σ**^{(n)} matrices are therefore the input to the procedure, defined by Eqs. (9) and (10). One can then use the definition of the three-term recurrence in Eq. (20) to express all definitions in terms of these moments, without formal reference to the large auxiliary space quantities (**d** or **W**), as

where the permutation operator *P* is defined as *P*(**Z**) = **Z** + **Z**^{†}. For a Hermitian theory, we can write $Si,j=Sj,i\u2020$ and assume the zeroth order Lanczos vector to be zero, i.e., $Si,0(n)=S0,j(n)=S0,0(n)=0$. Additionally, the orthogonality of the Lanczos vectors requires that $Si,j(0)=\delta ijI$. By considering again the Cholesky QR algorithm, the off-diagonal **C** matrices can therefore be computed as

and the on-diagonal **M** matrices can be found using Eq. (21),

These recurrence relations allow the calculation of the on- and off-diagonal blocks resulting in the block tridiagonal form of the Hamiltonian in Eq. (16). Despite the apparent complexity of the recurrence relations, this algorithm contains no step scaling greater than $O[N3]$, 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. (14) 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. (12), requires the combination of the block tridiagonal Hamiltonian in Eq. (16) resulting from both the hole and particle self-energy moments. This is constructed as

where $W\u0303$ are equal to the **L** matrix padded by zeros

and $d\u0303$ are defined by the block tridiagonal elements

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 $d\u0303$ and appropriately rotating $W\u0303$ into this basis. The eigenvalues of $H\u0303$ 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 2*j* 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**_{1}, with modified couplings **L** to the physical space. As such, in this notation, the number of conserved moments equals 2*n*_{iter} + 2, i.e., up to and including the 2*n*_{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*(2*n*_{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*^{2}.

## 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.^{77} 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

where the left and right eigenvectors form the biorthogonal set as

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

The **A** and **B** matrices are defined as

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, $Kia,jb=(ia|jb)=Kia,bj$. 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

Note that this matrix is formally equivalent to the alternative directly dynamical construction from $\chi (\omega )=(P(\omega )\u22121\u2212K)\u22121$, 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. (35)], we can write a more compact form of the dd-response as

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

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. (8), which can be written more compactly in matrix form as

and constitutes a central object of interest in this work, required for the GW self-energy moment construction of Eqs. (9) and (10).^{84}

We now show that the RPA can be entirely reformulated in terms of the dd-moments [Eq. (38)] 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. (30), along with insertion of a resolution of the identity via Eq. (32), we find a relation between the first two dd-moments as

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

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

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

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

While these are important relations in themselves, they also illustrate that all the RPA excitations and weights in the dd-response of Eq. (36) are implicitly accessible without requiring the explicit solution to the Casida equation (**X**, **Y**, and **Ω** matrices). This reformulation just requires knowledge of *η*^{(0)} as the central variable [which can be defined independently of the original equations via Eq. (40)], the **A** and **B** matrices defining the system and their interactions, and the recursive relations of Eq. (46). 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 *η*^{(0)} = **I**, and *η*^{(n)} = **A**^{n}. 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

Related expressions for the RPA correlation energy can be found in Eqs. (39) and (40) of Ref. 86 in terms of the quantities *η*^{(0)}, **A** + **B** and **A** − **B** (where there, these quantities are referred to as $QIdRPA$, ** ɛ**, and

**+ 2**

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

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[N4]$ evaluation of these quantities that are central to the approach in this work, thus avoiding their formal $O[N6]$ 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,^{56} 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

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

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

where the brackets indicate the order of contractions in order to preserve $O[N4]$ scaling, and Einstein summation is implied. The lesser self-energy moment of Eq. (9) can be recast in an analogous fashion.^{88}

Obtaining all dd-moments of the form of Eq. (52) up to order *n* can be simply reduced to knowledge of the first two moments $\eta \u0303ia,P(0)$ and $\eta \u0303ia,P(1)$, via use of the recursive relationship between the moments as given by Eq. (45), as for even moment orders,

and for odd moment orders,

where we have omitted explicit indices for brevity. Ensuring that an $O[N4]$ 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. (33) and (34)]

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

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

Where, for the RPA, **D** is defined by Eq. (56), **S**_{L} = **DV** and **S**_{R} = 2**V**. 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. (54) and (55), a fixed number of dd-moments of the form of Eq. (52) can be found in $O[N4]$ time, provided the original $\eta \u0303(0)$ and $\eta \u0303(1)$ values are known.

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

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

from which we can find an expression for *η*^{(0)} as

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

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

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[Naux3+Naux2ov]$ time.

Having constructed **T**, we can complete the evaluation of $\eta \u0303(0)$ using the definition of the matrix square-root as an integration in the complex plane,^{91}

From Eq. (62), this results in

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,

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[Naux2ov]$ time. This casts Eq. (67) into the form

For each value of the integration variable *z*, the integrand is a matrix of size *ov* × *N*_{aux}, which can be constructed in $O[N4]$ scaling, rendering it efficient for numerical quadrature. Along with the results of Eqs. (54), (55), (58), and (59), 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. (52), in no more than $O[N4]$ scaling (and $O[N3]$ 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[Naux2ov+Naux3]$ expression to evaluate for the zeroth-order dd-moment as

The first numerical integral in Eq. (71) (where the integrand decays exponentially) is computed via Gauss–Laguerre quadrature, while the second (where the integrand decays as $O[z\u22124]$) 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.

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\u2212B)(A+B)]12$, rather than the full integrand. We write a diagonal approximation to (**A** − **B**)(**A** + **B**) as

with

This is the same form as Eq. (58), 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[ov]$ computational time. We then seek to ensure the trace of the difference between the exact and numerically integrated estimate of $MD12$ vanishes. This is achieved for both integrals in Eq. (71), with the analytic form for the first integral given by $12\pi diag(D\u22121SLSRT)=Ioffset$ and the second numerical integral as $MD12\u2212D\u2212Ioffset=Iint$.

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

where Eq. (75) is minimized to optimize the grid for the first integral of Eq. (71), and Eq. (76) 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[ov]$ 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

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.^{94} 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[N3]$ time.^{97,98}

Once found, the *Z*_{PQ} can be symmetrically decomposed as *Z*_{PQ} = *Y*_{PR}*Y*_{QR}. By replacing the density-fitted integral tensor *V*_{iaR} in the above expressions with the fully factorized form *X*_{iP}*X*_{aP}*Y*_{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[N3]$ with formation of appropriate intermediates. We also require the numerical computation of the partially transformed dd-moment, $ZQSXaSXiS\eta ia,jb(0)XjRXbRZPR$. 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*.^{93} While this becomes a two-dimensional numerical integral, optimal quadrature grids can be calculated in a minimax sense.^{100,101}

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

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

where both intermediates

can be evaluated in $O[N3]$ cost. Further contractions in the evaluation of Eq. (71) 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*_{0}*W*_{0} 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.^{11} In Fig. 2, we consider the convergence of the first ionization potential (IP) of singlet O_{2} with conserved self-energy moment order. We compare this to two *G*_{0}*W*_{0} implementations, one of which performs an exact frequency integration (denoted “full QP-*G*_{0}*W*_{0},” which scales as $O[N6]$), 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*_{0}*W*_{0},” which scales as $O[N4]$).^{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.

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[N4]$, 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*_{0}*W*_{0} across a larger test set to consider the moment truncation convergence. We use the established “*GW*100” 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 $\u223c4\u221225$ 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*_{0}*W*_{0}. 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*_{0}*W*_{0} 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*_{0}*W*_{0} 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*_{0}*W*_{0} limit, the general trend, convergence, and level of accuracy, which can be reached with moment order is likely to be similar.

It is important to put the scale of the moment truncation convergence in the context of the overall accuracy of the *G*_{0}*W*_{0} method for these systems. In Fig. 4, we therefore compare the aggregated mean absolute error (MAE) in the moment-truncated *G*_{0}*W*_{0} IP and EA values over this *GW*100 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*_{0}*W*_{0} 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).^{46} 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*_{0}*W*_{0} comparison, albeit noting the other potential sources of discrepancy between these values discussed earlier.

### 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*_{0}*W*_{0} 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[N5]$, in general, to obtain the full spectrum.^{41} 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 *GW*100 benchmark set, and compared to the AC-*G*_{0}*W*_{0} full-frequency spectrum.

The convergence of the spectrum is shown for a series of conserved *G*_{0}*W*_{0} moment orders, from the HF level up to the AC-*G*_{0}*W*_{0} spectrum. The AC-*G*_{0}*W*_{0} 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*_{0}*W*_{0} 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*_{0}*W*_{0} spectrum with increasing numbers of included moments.

This Wasserstein metric convergence plateaus at the $\u223c7th$-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*_{0}*W*_{0}, 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*_{0}*W*_{0} 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 $Zp=(1\u2212\u2202\Sigma pp(\u03f5p)\u2202\omega )\u22121$, 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_{2} dimer in a 6-31G basis set. Figure 6 shows quasiparticle energies, self-energies, and quasiparticle renormalization factors for AC-*G*_{0}*W*_{0}, and the moment conserving *G*_{0}*W*_{0} 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.

The figure shows the HOMO and first three unoccupied states found in this system, with the AC-*G*_{0}*W*_{0} (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*_{0}*W*_{0} 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*_{0}*W*_{0} 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*_{0}*W*_{0}; 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[N4]$ 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 *GW*100 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*_{0}*W*_{0} 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.^{112}

## 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 *GW*100 benchmark.

### APPENDIX A: IMPROVED NUMERICAL QUADRATURE

While the integral expression for $\eta \u0303(0)$ derived in Eq. (70) in the main text is sufficient for numerical quadrature in $O[N4]$ 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;^{56} however, there are also important differences.

We can first consider a mean-field contribution to the integral component of $\eta \u0303(0)$, i.e., setting the Coulomb interaction to zero in the expression for $[(A\u2212B)(A+B)]12$. This just leaves a contribution from the irreducible polarizability, which can be analytically integrated as

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

which substantially reduces the magnitude of the numerically integrated component of $\eta \u0303(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. (68) and (69) that both **F**(*z*) and **Q**(*z*) decay as *z*^{−2} in the large-*z* limit. Series expanding [**I** + **Q**(*z*)]^{−1} ∼ **I** − **Q**(*z*) + **Q**^{2}(*z*) − …, we find that the leading order decay of the integrand in Eq. (A2) is $O[z\u22122]$. 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\u22124]$ rate.

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

This can be analytically integrated via the residue theorem to give

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

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

However, in this form, it is unfortunately unable to be efficiently evaluated in $O[N4]$ time (instead scaling as $O[N5]$). 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.^{56}

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

which can be verified by substitution of Eq. (A5). Writing the matrix of interest we wish to efficiently compute as $N=(SLSRT)\u25e6E$, this Sylvester equation has an integral solution of the form^{113}

This can be verified via substitution into Eq. (A7), as

The numerical integration of Eq. (A8) can therefore be used to compute the leading order term in the numerical integration of Eq. (A2). 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. (A5).^{114} It may appear as though we have swapped one numerical integral for another; however, the transformation to the form given in Eq. (A8) 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

The integrand of the first integral decays exponentially and is efficiently computed with Gauss–Laguerre quadrature, with the second integral (decaying as $O[z\u22124]$) is evaluated with Clenshaw–Curtis quadrature.

### APPENDIX B: NUMERICAL INTEGRATION ERROR ESTIMATES

We make use of two separate approaches to estimate the error in $\eta \u0303(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 $\eta \u0303np(0)=\eta \u0303\u221e(0)+\Delta \eta \u0303np(0)$.

For our first estimate, we begin by noting from Eq. (40) that the exact $\eta \u0303(0)$ will satisfy the relation

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

We can then explicitly evaluate ‖*X*(*n*_{p})‖_{2} in $O[N4]$ time under the same conditions as the rest of this manuscript and apply the properties of the Frobenius norm to write

Writing $\Vert \Delta \eta \u0303np(0)\Vert 2=x$, the error is therefore rigorously bound by the solutions satisfying the quadratic inequality,

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[N4]$ 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 $\eta \u0303np(0)A+B$, which we can relate to the error in $\eta \u0303np(0)$. Parameterizing the exponential convergence of the *L*_{2} norm error as $\Vert \Delta \eta \u0303np(0)\Vert 2=\alpha e\u2212\beta np$, we can rigorously relate the differences between estimates at different levels of sampling, as

We seek to estimate the value of $\Vert \Delta \eta \u03034np(0)\Vert 2$ given the values $\Delta 4np,2np$ and $\Delta 4np,np$, which are computed as intermediates in the nested quadrature of 4*n*_{p} points.

Setting $x=e\u2212\beta np$, using the approximation defined in Eq. (B7) for $\Delta 4np,2np=\alpha (x4+x2)$ and $\Delta 4np,np=\alpha (x4+x)$, and eliminating the factor of *α*, we find

The smallest positive, real solution to this is used to estimate the error (given as *αx*^{4}) as

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 $\Delta 4np,2np$ and $\Delta 4np,np$, 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. (B9) denoted “nested fit,” while the lower bound from the solution of Eq. (B4) denoted “lower bound.”

## REFERENCES

^{3}) implementation of Hedin’s GW approximation for molecules

_{1}00: A plane wave perspective for small molecules

_{0}results for the electron gas: Fixed screened potential W

_{0}within the random-phase approximation

_{m}ClnF

_{4−}m−n

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 $(\eta (n))\u22121=(X\u2212Y)\Omega \u2212n(X\u2212Y)T$.

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

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

^{6}) to O(r

^{4})