A 'moment-conserving' reformulation of GW theory

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


I. INTRODUCTION
Despite the phenomenal success of density functional theory (DFT) in electronic structure, its standard approach is both conceptually and (often) practically illsuited 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 which 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][3][4][5][6][7][8][9][10][11][12][13][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][18][19][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 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 which 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][46][47][48][49] . This order can be systematically increased to more finely resolve the full dynamical dependence of the GW selfenergy. The dynamical information is therefore implicitly recast into a small number of static matrices, each of which can be obtained in O[N 4 ] time (with a proposed O[N 3 ] algorithm also given). This removes the need for the definition of any frequency or time grids in which to resolve dynamical quantities, spectral broadening, finite temperatures, Fourier transforms or analytic continuation, with all dynamics implicitly represented by this series of static quantities.
Furthermore, once these spectral moments of the selfenergy are obtained, 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 use in numerical approaches 45,50 , characterizing sum rules 51,52 , and as physical observables in their own right [53][54][55] .
In Sec. II we show how these spectral moments of the self-energy can be directly constructed from moments of the 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 show how they can be efficiently obtained in O[N 4 ] cost, based on ideas from the seminal work in 2010 by Furche and collaborators on low-scaling approaches for the RPA correlation energy 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[N 2 ] neutral excitations of RPA theory (representing the poles of the screened Coulomb interaction) and the charged excitations of the reference Green's function. In the absence of self-consistency (the most common 'G 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 frequencydomain as (1) In these expressions, Ω ν are the neutral excitation energies from RPA theory, and X ν ia and Y ν ia 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 meanfield, 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[N 6 ], rendering it unsuitable for large-scale implementations. However, in this work, we are only interested in evaluating the spectral moments of the resulting selfenergy, 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 densitydensity 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 we can rewrite Eqs. 6 and 7 as ia,jb (qk|jb) ia,jb (qc|jb).
Evaluating the self-energy spectral moments of Eqs. 9-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, η (n) ia,jb . These relations are subsequently used in Sec. IV to demonstrate how the self-energy moments can be evaluated in O[N 4 ] scaling, with Sec. V going further to propose a cubic scaling algorithm for their evaluation (and therefore full GW algorithm). In addition to these moments representing the dynamical part of the self-energy, we also require a static (exchange) part of the self-energy, Σ ∞ , which can be calculated as 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 selfconsistent 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 selfenergy). Specifically, we seek an effective static Hamiltonian,H 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 To find this effective upfolded representation of the moment-conserving dynamics, we modify the block Lanczos procedure to ensure the construction aH of minimal size, whose effective hole and particle self-energy moments exactly match the ones from Eqs. 9-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 selfenergy 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 defineq (j) as The q (j) are block Lanczos vectors of depth j, which form a recursive Krylov space as q (j) = q 1 q 2 · · · q j , ensuring that when taken together, they project to a block tridiagonal representation of the upfolded selfenergy 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 factorisation 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 initialisation of the recurrence terms, and have a dimension which 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-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 S i,j = S † j,i , and assume the zeroth order Lanczos vector to be zero i.e. S i,j = δ ij I. By considering again the Cholesky QR algorithm, the offdiagonal 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 onand 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[N 3 ], by eliminating the explicit reference to the full upfolded Hamiltonian, whilst 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. diagonalisation of Eq. 12, requires the combination of the the block tridiagonal Hamiltonian in Eq. 16 resulting from both the hole and particle self-energy moments. This is constructed as whereW are equal to the L matrix padded by zeros andd are defined by the block tridiagonal elements This ensures 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 diagonalising d and appropriately rotatingW into this basis. The eigenvalues ofH are moment-conserving approximations to those of the exact upfolded Hamiltonian, and the corresponding eigenvectors u can be transformed into Dyson orbitals via LPu, where P is a projection into the physical space, and the L is required to transform the physical component of the eigenvectors back to the MO representation. This process conserves exactly the first 2j hole and particle self-energy moments. Commonly, a notation referring to the number of iterations of the block Lanczos recurrence n iter is used; in this notation the n iter = 0 calculation corresponds to the inclusion of only a single on-diagonal block M 1 , with modified couplings L to the physical space. As such, in this notation the number of conserved moments equals 2n iter + 2, i.e. up to and including the 2n iter + 1 order moment. This is the same number of moments as required as input to the recurrence relations, and therefore the algorithm conserves all the moments used as input, which should be up to an odd order. After n iter applications of this algorithm to both the lesser and greater self-energy sectors, this results in N (2n iter + 3) quasiparticle states (demonstrating the potential to capture satellite features with these additional poles). As such, application of this algorithm becomes theoretically equivalent to a full diagonalisation 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-10 can be computed efficiently. As a first step towards 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 ddmoments 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][79][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 which couples particlehole excitations and de-excitations. In the traditional RPA (without second-order exchange), this coupling is taken to be the same for excitations, de-excitations and their coupling, given by the static, bare Coulomb interaction, K ia,jb = (ia|jb) = K ia,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 χ(ω) = (P(ω) −1 − K) −1 , where P(ω) is the irreducible polarizability of the reference state. Considering the positivefrequency 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 holeparticle 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-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. Firstly, 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 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 meanfield 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 Overall, this perspective on the RPA in terms of ddmoments is key to open new avenues such as the ones explored in this work.

IV. EFFICIENT EVALUATION OF SELF-ENERGY AND DENSITY RESPONSE MOMENTS
Given the recasting of the RPA dd-response in Sec. III in terms of its lowest order moment (Eq. 40) and recursion to access the higher moments via Eq. 46, we now consider the efficient O[N 4 ] evaluation of these quantities which are central to the approach in this work, thus avoiding their formal O[N 6 ] construction via Eq. 38. The derivation here is heavily inspired by the seminal 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 quantitỹ 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[N 4 ] 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η ia,P , via use of the recursive relationship between the moments as given by Eq. 45, as for even moment orders, and for odd moment orders, 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 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[N 3 aux + N 2 aux ov] time. Having constructed T, we can complete the evaluation ofη (0) using the definition of the matrix square-root as an integration in the complex plane 91 , From Eq. 62, this results iñ 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 which can be constructed in O[N 2 aux ov] time. This casts Eq. 67 into the form 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 which can be analytically integrated. This derivation is given in Appendix A, and results in a final O[N 2 aux ov + N 3 aux ] expression to evaluate for the zeroth-order dd-moment as The 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 M D 1 2 vanishes. This is achieved for both integrals in Eq. (71), with the analytic form for the first integral given by 1 2 diag(D −1 S L S T R ) = I offset and the second numerical integral as M D 1 2 − D − I offset = I int . Writing this explicitly, given quadrature points and weights {z i , w i } for a n p -point infinite or semi-infinite quadrature, we seek to scale our points by a factor a, which is a root the objective functions 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 App. 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 lowestscaling 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 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 . Use of this doubly-factorized form has also been previously suggested in the use of a reduced-scaling RPA and particle-particle RPA schemes 95,96 . This form for the integrals can be directly constructed with controllable errors in O[N 3 ] time 97,98 . Once found, the Z P Q can be symmetrically decomposed as Z P Q = Y P R Y QR . By replacing the densityfitted 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[N 3 ] with formation of appropriate intermediates. We also require the numerical computation of the partially transformed dd-moment, Z QS X aS X iS η (0) ia,jb X jR X bR Z P R . Inspired by the lowscaling approach taken to RPA correlation energies in the work of Refs. 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 B QR (p) = X iQ X aQ i e − a p e ip X iR X aR (82) can be evaluated in O[N 3 ] 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 doublyfactorized 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 forthcoming work, with numerical results in the rest of this work employing the quartic scaling algorithm described in Sec. IV.

A. Comparison to quasiparticle GW approaches
We first consider the convergence of the momenttruncated G 0 W 0 algorithm compared to more traditional implementations, as found in the PySCF software package [102][103][104][105] . As found to be more effective for molecular systems due to the importance of exact exchange,  2: Convergence of the first G 0 W 0 IP of singlet O 2 in a cc-pVDZ basis with respect to the number of conserved self-energy moments. Also shown is the same quasiparticle state computed from traditional G 0 W 0 via an exact frequency integration ('Full') and analytic continuation approach to G 0 W 0 ('AC'), albeit both imposing a diagonal approximation in their solution to the QP equation. In the moment expansion, we also consider a diagonal approximation, by explicitly removing the non-diagonal parts of our computed self-energy moments ('Diagonal Σ moment'). All approaches find an IP of 8.49 eV to within 10 meV, with the difference between full-frequency and AC approaches 5 meV, and the relaxation of diagonal approximation also accounting for small ∼ 5 meV differences, with the reference Hartree-Fock IP for comparison being 9.73 eV.
we perform all calculations on a restricted Hartree-Fock reference 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[N 6 ]), and one which performs an analytic continuation (AC) of the imaginary frequency selfenergy to real-frequencies via a fit to Padé approximants in order to perform the convolution (denoted 'AC QP- 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) selfenergy moment approach only very slightly different, indicating the relative unimportance of the non-diagonal self-energy components in this system and lack of significant correlation-induced coupling between the mean-field quasiparticle states. Furthermore, the analytic continuation approach is also highly accurate for this system, introducing an error of only 5 meV compared to the full frequency integration. We have furthermore numerically verified that our approach scales as O[N 4 ], with computational cost comparable to the analytic continuation approach.
Having demonstrated correctness compared to a highscaling 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 ∼ 4 − 25 eV, and featuring molecules with bound metal atoms (including metal diatomics and clusters), strongly ionic bonding, and molecules with a strongly delocalised 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 11 th -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 while 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 trunca- tion 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 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 behaviour for moment-truncated GF2 theory (where lower order selfenergy 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 towards accuracy in any particular energy range, making it suitable for G 0 W 0 excitations beyond frontier excitations. Description of low-lying states is a particular challenge for many other GW approaches, with analytic continuation becoming less reliable, and alternatives like contour-deformation scaling as O[N 5 ] in general to obtain the full spectrum 41 . In Fig. 5 we therefore show a respect to the AC-G 0 W 0 spectrum, quantifying the difference between the spectral distributions. The AC-G 0 W 0 spectrum is indicated transparently behind the other spectra to ease visualisation of the convergence, and was calculated using an iterative diagonal approximation to the quasiparticle equation.
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 accross 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 ∼ 7 th -order conserved moments, with further orders not improving this metric further. This could be due to numerical precision of the algorithm, or fundamental approximations in the AC-G 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 which 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 weaklycorrelated 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 the next section.

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 rootfinding 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 renomalization factor evaluated at the quasi-particle energy, defined as ) −1 , which indicates the approximate weight in the quasi-particle solutions.
Since the moment-conserving GW approach obtains all poles in one step (including satellites and low-weighted features), all the excitations can be characterized by their quasiparticle weight, and either selected as specific excitation energies, or all excitations included in the spectrum to ensure smooth changes with molecular geome-try. 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 renormalisation factors for AC-G 0 W 0 , and the moment conserving G 0 W 0 approach with both up to the 1 st -and 11 th -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 11 th -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 momentconserving 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 the an iterative solution to the quasiparticle equation for each state (avoiding 'spurious' solutions), diagonal self-energy ap-proximations, or requirements for analytic continuation of dynamical quantities. It allows for all excitation energies and weights to be obtained directly in a non-iterative single diagonalization of a small effective Hamiltonian, controlled by a single parameter governing the number of conserved self-energy moments. Full RPA screening and particle-hole coupling in the self-energy is included, which is captured with O[N 4 ] computational scaling via a numerical one-dimensional integration, with a reduction to cubic scaling also proposed. This approach is enabled by a recasting of the RPA in terms of the moments of the density-density response function. Applied across we 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 forwards, 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 .

CODE AND 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η (0) derived in Eq. 70 in the main text is sufficient for numerical quadrature in O[N 4 ] scaling, it can be refined by deducting large and/or slowly-decaying contributions which can be efficiently integrated separately, in order to minimize the number of quadrature points which are required to obtain a given overall accuracy. This again follows many of the developments in the 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η  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 obtaiñ which substantially reduces the magnitude of the numerically integrated component ofη (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-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 −2 ]. We consider this contribution in isolation, noting that if it can be removed from the numerical integration, the next leading order will result in the integrand decaying at an improved O[z −4 ] rate.
This leading-order contribution to the integral of Eq. 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 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 which 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 −4 ]) is evaluated with Clenshaw-Curtis quadrature.

APPENDIX B: NUMERICAL INTEGRATION ERROR ESTIMATES
We make use of two separate approaches to estimate the error inη (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 η Writing ||∆η (0) np || 2 = x, the error is therefore rigorously bound by the solutions satisfying the quadratic inequality, np || 2 ||(A + B)|| 2 x + ||X(n p )|| 2 ||(A + B)|| 2 ≤ 0. (B4) In practice, while we find that the two roots to this equation rigorously bound the error in the quadrature, their values are not a tight bound on the true error (in particular, the upper bound is excessively large to be used as an estimate, though the lower bound on the error that it provides is reasonable, shown in Fig. 1 of the main text as 'Lower Bound'). Combined with the O[N 4 ] scaling of evaluating X(n p ) this makes this estimate of limited practical use, and we seek an improved estimate, that provides a tighter upper bound on the true error.
To do this we take advantage of the nested nature of the Clenshaw-Curtis quadrature within the integral providing the dominant error contribution. Using this, we can obtain estimates at both one half and one quarter of the current number of integration points with no additional cost. We show how this can give an estimate of the error inη np || 2 = αe −βnp , we can rigorously relate the differences between estimates at different levels of sampling, as We seek to estimate the value of ||∆η (0) 4np || 2 given the values ∆ 4np,2np and ∆ 4np,np which are computed as intermediates in the nested quadrature of 4n p points. Setting x = e −βnp , and using the approximation defined in Eq. B7 for ∆ 4np,2np = α(x 4 + x 2 ) and ∆ 4np,np = α(x 4 + x), 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 ∆ 4np,2np and ∆ 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'.