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

*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

_{ν}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.

_{ν}) 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

*μ*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.

*GW*self-energy of Eq. (1), the moments can be constructed as

*n*th-order spectral moments of the RPA density-response, summed over both particle–hole excitation and de-excitation components, as

*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

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

*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

^{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.

*exact*upfolded self-energy representation,

^{63}via inspection from Eq. (1), with

*GW*self-energy components. We now consider the projection of the exact

*GW*upfolded matrix representation into a truncated block tridiagonal form, as

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

**L**and Lanczos vectors

**q**

_{1}can be found via a QR factorization of the exact

*GW*couplings

**W**, as

**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

**q**

_{1}=

**W**

^{†}

**L**

^{−1,†}. Subsequent block Lanczos vectors can then be defined according to the standard three-term recurrence

**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

*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

**M**matrices can be found using Eq. (21),

^{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

**L**matrix padded by zeros

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

^{81,82}as a generalized eigenvalue decomposition

**X**+

**Y**) and (

**X**−

**Y**)

^{T}, as

**A**and

**B**matrices are defined as

*ϵ*

_{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.

**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

^{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

*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

^{84}

**A**+

**B**) = (

**X**−

**Y**)

**Ω**(

**X**−

**Y**)

^{T}. By taking the sum and difference of the two Casida equations of Eq. (30), we also find

^{85}By right-multiplying by (

**X**+

**Y**)

^{T}, this leads to a relation between the zeroth and second dd-moment, as

**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

*η*^{(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

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

*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

^{88}

*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,

**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)]

**A**+

**B**) into an appropriate form as

**A**−

**B**)(

**A**+

**B**) in a general fashion as a diagonal plus asymmetric low-rank part, as

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

*η*

^{(0)}as

**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

**T**. We can use the Woodbury matrix identity to rewrite it as

**D**, and a matrix of dimension

*N*

_{aux}, with the overall

*ov*×

*N*

_{aux}matrix able to be constructed in $O[Naux3+Naux2ov]$ time.

**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}

**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

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

**A**−

**B**)(

**A**+

**B**) as

*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$.

*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

## V. REDUCTION TO CUBIC-SCALING GW

*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

*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}

**Q**(

*z*) of Eq. (69) to be obtained as

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.

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

**E**as

^{56}

^{113}

^{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.

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

*X*(

*n*

_{p}). Expanding the estimated quantity in terms of the exact result and error gives

*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

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

*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

*n*

_{p}points.

*α*, we find

*αx*

^{4}) as

## REFERENCES

*GW*approximation

*Solid State Physics*

*GW*approximation

*Computational Nanoscience: Do It Yourself!*

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

*GW*self-energy or from random phase approximation total energies

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

*GW*100: Benchmarking

*G*0

*W*0 for molecular systems

*GW*compendium: A practical guide to theoretical photoemission spectroscopy

*Solid State Physics*

*Recent Advances in Density Functional Methods*

*O*(

*N*

^{3}) implementation of Hedin’s

*GW*approximation for molecules

*GW*approximation for molecules

*GW*: Towards fast quasiparticle calculations

*GW*methods

*GW*

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

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

*GW*

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

*W*

_{0}within the random-phase approximation

*GW*approximation: Application to Si, MnO, and NiO

*GW*theory

*GW*calculations for atoms and molecules

*GW*method: A basis for the independent-particle approximation

*GW*calculations for semiconductors and insulators

*GW*approach

*GW*for molecules

*G*0

*W*0 approximation

*GW*calculations

*GW*calculations with a separable density-fitting space–time approach

*GW*in Gaussian Bloch orbitals for solids

*GW*with numeric atom-centered orbital basis functions

_{m}Cl

*n*F

_{4−}

*m*−

*n*

*GW*methods

*Many-Body Methods for Atoms, Molecules and Clusters*

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

*A*

^{α}, log(

*A*), and related matrix functions by contour integrals

*G*0

*W*0 by pair atomic density fitting

*O*(

*r*

^{6}) to

*O*(

*r*

^{4})

*Ab initio*full cell

*GW*+ DMFT for correlated materials

*G*0

*W*0 for valence and core excitation energies of periodic systems

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

*GW*approaches for the

*GW*100 test set

*GW*approximation

*GW*-based methods using the Hubbard dimer

*GW*methods

*AX*−

*XB*=

*Y*