# A ‘moment-conserving’ reformulation of GW theory

Charles J. C. Scott<EMAIL_ADDRESS>Department of Physics, King’s
College London, Strand, London WC2R 2LS, U.K.  Oliver J. Backhouse
<EMAIL_ADDRESS>Department of Physics, King’s College London,
Strand, London WC2R 2LS, U.K.  George H. Booth<EMAIL_ADDRESS>Department of Physics, King’s College London, Strand, London WC2R 2LS, U.K.

###### Abstract

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 $\mathcal{O}[N^{4}]$ approach to build and solve
this Hamiltonian, with a proposal to reduce this further to
$\mathcal{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 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
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 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(\omega)$, and the screened Coulomb interaction $W_{p}(\omega)$,
obtained (in general) at the RPA level of theory. This provides the dynamical
part of the self-energy, $\Sigma(\omega)$, formally written as
$\Sigma(\omega)=(i/2\pi)\int
d\omega^{\prime}e^{i\eta\omega^{\prime}}G(\omega+\omega^{\prime})W_{p}(\omega^{\prime})$.
There are many different variants of $GW$ theory[21, 3, 13, 14, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43,
44], which primarily differ due to i) the choice (or absence) of self-
consistency conditions on $G(\omega)$ and/or $W_{p}(\omega)$[21], ii) the
approach to find the quasi-particle energies once $\Sigma(\omega)$ 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 $\mathcal{O}[N^{6}]$ scaling
step, required to find the entire set of poles in $W_{p}(\omega)$ from the
RPA. However, there are a number of techniques to approximate this frequency
integration which can reduce this scaling (generally down to
$\mathcal{O}[N^{3}-N^{4}]$) based on plasmon pole approximations, analytic
continuation, contour deformation and explicit grid resolved approaches for
the dynamics of these quantities, amongst 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 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$ self-energy. The dynamical information is therefore implicitly
recast into a small number of static matrices, each of which can be obtained
in $\mathcal{O}[N^{4}]$ time (with a proposed $\mathcal{O}[N^{3}]$ algorithm
also given). This removes the need for the definition of any frequency or time
grids in which to resolve dynamical quantities, spectral broadening, finite
temperatures, Fourier transforms or analytic continuation, with all dynamics
implicitly represented by this series of static quantities.

Furthermore, once these spectral moments of the self-energy are obtained, 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
[50, 45], 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 $\mathcal{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 $GW100$ test set frequently used to benchmark $GW$
implementations, demonstrating a rapid convergence of the moment expansion,
and accurate and efficient results across this test set for the $G_{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(\omega)$ and $W_{p}(\omega)$, can be formally expanded
as a sum over the $\mathcal{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 $\mathcal{O}[N]$ mean-field molecular orbital energies. This allows the
self-energy to be explicitly evaluated in the frequency-domain as

$\displaystyle\Sigma_{pq}$
$\displaystyle(\omega)=\sum_{\nu}\sum_{ia,jb,k}\frac{(pk|ia)(X_{ia}^{\nu}+Y_{ia}^{\nu})(X_{jb}^{\nu}+Y_{jb}^{\nu})(qk|jb)}{\omega-(\epsilon_{k}-\Omega_{\nu})-i0^{+}}$
$\displaystyle+\sum_{\nu}\sum_{ia,jb,c}\frac{(pc|ia)(X_{ia}^{\nu}+Y_{ia}^{\nu})(X_{jb}^{\nu}+Y_{jb}^{\nu})(qc|jb)}{\omega-(\epsilon_{c}+\Omega_{\nu})+i0^{+}}.$
(1)

In these expressions, $\Omega_{\nu}$ are the neutral excitation energies from
RPA theory, and $X_{ia}^{\nu}$ and $Y_{ia}^{\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 $\epsilon_{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 ($\Omega_{\nu}$) scales as
$\mathcal{O}[N^{6}]$, rendering it unsuitable for large-scale implementations.
However, in this work, we are only interested in evaluating the spectral
moments of the resulting self-energy, and finding a resulting compressed
representation of the self-energy with fewer poles, but which by construction
matches the spectral moments up to a desired order. These frequency-
independent spectral moments are defined separately for the greater and lesser
parts, and represent the $n^{\mathrm{th}}$-order moments of the resulting
dynamical distributions, as

$\displaystyle\Sigma^{(n,<)}_{pq}$
$\displaystyle=-\frac{1}{\pi}\int_{-\infty}^{\mu}\mathrm{Im}[\Sigma(\omega)_{pq}]\omega^{n}d\omega$
(2)
$\displaystyle=(-1)^{n}\left.\frac{d^{n}\Sigma(\tau)_{pq}}{d\tau^{n}}\right|_{\tau=0^{+}},$
(3)

and similarly

$\displaystyle\Sigma^{(n,>)}_{pq}$
$\displaystyle=\frac{1}{\pi}\int_{\mu}^{\infty}\mathrm{Im}[\Sigma(\omega)_{pq}]\omega^{n}d\omega$
(4)
$\displaystyle=(-1)^{n}\left.\frac{d^{n}\Sigma(\tau)_{pq}}{d\tau^{n}}\right|_{\tau=0^{-}},$
(5)

where $\mu$ 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

$\displaystyle\Sigma_{pq}^{(n,<)}=\sum_{\nu}$
$\displaystyle\sum_{ia,jb}\sum_{k}\left[(pk|ia)(X_{ia}^{\nu}+Y_{ia}^{\nu})\right.$
$\displaystyle\left.(\epsilon_{k}-\Omega_{\nu})^{n}(X_{jb}^{\nu}+Y_{jb}^{\nu})(qk|jb)\right],$
(6) $\displaystyle\Sigma_{pq}^{(n,>)}=\sum_{\nu}$
$\displaystyle\sum_{ia,jb}\sum_{c}\left[(pc|ia)(X_{ia}^{\nu}+Y_{ia}^{\nu})\right.$
$\displaystyle\left.(\epsilon_{c}+\Omega_{\nu})^{n}(X_{jb}^{\nu}+Y_{jb}^{\nu})(qc|jb)\right].$
(7)

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

$\eta^{(n)}_{ia,jb}=\sum_{\nu}(X_{ia}^{\nu}+Y_{ia}^{\nu})\Omega_{\nu}^{n}(X_{jb}^{\nu}+Y_{jb}^{\nu}),$
(8)

we can rewrite Eqs. 6 and 7 as

$\displaystyle\Sigma_{pq}^{(n,<)}$
$\displaystyle=\sum_{ia,jb,k}\sum_{t=0}^{n}\binom{n}{t}(-1)^{t}\epsilon_{k}^{n-t}(pk|ia)\eta_{ia,jb}^{(t)}(qk|jb)$
(9) $\displaystyle\Sigma_{pq}^{(n,>)}$
$\displaystyle=\sum_{ia,jb,c}\sum_{t=0}^{n}\binom{n}{t}\epsilon_{c}^{n-t}(pc|ia)\eta_{ia,jb}^{(t)}(qc|jb).$
(10)

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(\omega)$ and $W_{p}(\omega)$
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^{(n)}_{ia,jb}$.
These relations are subsequently used in Sec. IV to demonstrate how the self-
energy moments can be evaluated in $\mathcal{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, $\mathbf{\Sigma}_{\infty}$, which can be calculated as

$\mathbf{\Sigma}_{\infty}=\mathbf{K}[\mathbf{D}]-\mathbf{V}_{\mathrm{xc}},$
(11)

where $\mathbf{K}[\mathbf{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, $\mathbf{f}[\mathbf{D}]$, with $\mathbf{V}_{\mathrm{xc}}$ being
the exchange-correlation potential used in $\mathbf{f}$. Note that for a
Hartree–Fock reference, this static self-energy contribution is zero.

###  II.1 Full $GW$ spectrum from self-energy moments

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

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

$\mathbf{\tilde{H}}=\begin{bmatrix}\mathbf{f}+\mathbf{\Sigma}_{\infty}&\mathbf{\tilde{W}}\\\
\mathbf{\tilde{W}}^{\dagger}&\mathbf{\tilde{d}}\end{bmatrix},$ (12)

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

$\mathbf{G}(\omega)=\left(\omega\mathbf{I}-\mathbf{f}+\mathbf{\Sigma}_{\infty}-\mathbf{\tilde{W}}(\omega\mathbf{I}-\mathbf{\tilde{d}})^{-1}\mathbf{\tilde{W}}^{\dagger}\right)^{-1}.$
(13)

Such upfolded representations have been considered previously in diagrammatic
theories, in a recasting of GF2 theory in terms of its moments [57, 46, 47] as
well as more recently to $GW$ amongst others[58, 59, 60, 61, 62, 48, 63]. For
‘exact’ $G_{0}W_{0}$, this auxiliary space (i.e. the dimension of
$\mathbf{d}$) must scale as $\mathcal{O}[N^{3}]$ [61, 60]. However, in the
moment truncation, $\mathbf{\tilde{W}}$ and $\mathbf{\tilde{d}}$ can be
directly constructed such that their effect exactly matches that of a
truncated set of conserved $GW$ self-energy moments (separately in the
particle and hole sectors), yet rigorously scales in dimension as
$\mathcal{O}[nN]$, where $n$ is the number of conserved self-energy moments.
This allows for a complete diagonalization of $\mathbf{\tilde{H}}$, obtaining
all excitations in a single shot, and a reconstruction of the full $GW$
Green’s function from its Lehmann representation in $\mathcal{O}[(nN)^{3}]$
computational effort, avoiding the need for any grids or iterative solutions
once $\mathbf{\tilde{H}}$ is found.

To find this effective upfolded representation of the moment-conserving
dynamics, we modify the block Lanczos procedure to ensure the construction a
$\mathbf{\tilde{H}}$ 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 self-
energy representation[63], via inspection from Eq. 1, with

$\displaystyle\mathbf{W}_{p,k\nu}$
$\displaystyle=\sum_{ia}(pk|ia)(X_{ia}^{\nu}+Y_{ia}^{\nu})$ (14)
$\displaystyle\mathbf{d}_{k\nu,l\nu^{\prime}}$
$\displaystyle=(\epsilon_{k}-\Omega_{\nu})\delta_{k,l}\delta_{\nu,\nu^{\prime}},$
(15)

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

$\displaystyle\tilde{\mathbf{H}}_{\mathrm{tri}}$
$\displaystyle=\tilde{\mathbf{q}}^{(j),\dagger}\begin{bmatrix}\mathbf{f}+\boldsymbol{\Sigma}_{\infty}&\mathbf{W}\\\
\mathbf{W}^{\dagger}&\mathbf{d}\end{bmatrix}\tilde{\mathbf{q}}^{(j)}$
$\displaystyle=\begin{bmatrix}\mathbf{f}+\boldsymbol{\Sigma}_{\infty}&\mathbf{L}&&&&\mathbf{0}\\\
\mathbf{L}^{\dagger}&\mathbf{M}_{1}&\mathbf{C}_{1}&&&\\\
&\mathbf{C}_{1}^{\dagger}&\mathbf{M}_{2}&\mathbf{C}_{2}&&\\\
&&\mathbf{C}_{2}^{\dagger}&\mathbf{M}_{3}&\ddots&\\\
&&&\ddots&\ddots&\mathbf{C}_{j-1}\\\
\mathbf{0}&&&&\mathbf{C}_{j-1}^{\dagger}&\mathbf{M}_{j}\end{bmatrix},$ (16)

where we define $\tilde{\mathbf{q}}^{(j)}$ as

$\displaystyle\tilde{\mathbf{q}}^{(j)}$
$\displaystyle=\begin{bmatrix}\mathbf{I}&\mathbf{0}\\\
\mathbf{0}&\mathbf{q}^{(j)}\end{bmatrix}.$ (17)

The $\mathbf{q}^{(j)}$ are block Lanczos vectors of depth $j$, which form a
recursive Krylov space as
$\mathbf{q}^{(j)}=\begin{bmatrix}\mathbf{q}_{1}&\mathbf{q}_{2}&\cdots&\mathbf{q}_{j}\end{bmatrix}$,
ensuring that when taken together, they project to a block tridiagonal
representation of the upfolded self-energy with $j$ on-diagonal blocks as
shown. The action of this block Lanczos tridiagonalization of the upfolded
(hole or particle) self-energy is to exactly conserve these spectral moments
of the self-energy[64, 65]. This block tridiagonal representation is
equivalent to a truncated continued fraction form[50, 66], widely used in the
representation of dynamical quantities[67, 68], and even as an expansion
previously considered within $GW$ theory [45]. We therefore seek to
reformulate the Lanczos recursion in terms of just these moments, rather than
the action of the full upfolded Hamiltonian which we seek to avoid due to its
scaling.

The initial couplings $\mathbf{L}$ and Lanczos vectors $\mathbf{q}_{1}$ can be
found via a QR factorisation of the exact $GW$ couplings $\mathbf{W}$, as

$\displaystyle\mathbf{W}^{\dagger}$
$\displaystyle=\mathbf{q}_{1}\mathbf{L}^{\dagger}.$ (18)

However, this will scale poorly, and so we can rewrite this to directly
compute $\mathbf{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 $\mathbf{L}$ to the zeroth order self-energy
moment as

$\displaystyle\mathbf{L}^{\dagger}$
$\displaystyle=\left(\mathbf{W}\mathbf{W}^{\dagger}\right)^{\frac{1}{2}}=\left(\boldsymbol{\Sigma}^{(0)}\right)^{\frac{1}{2}},$
(19)

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

$\displaystyle\mathbf{q}_{i+1}\mathbf{C}_{i}^{\dagger}$
$\displaystyle=\left[\mathbf{d}\mathbf{q}_{i}-\mathbf{q}_{i}\mathbf{M}_{i}-\mathbf{q}_{i-1}\mathbf{C}_{i-1}\right]$
(20)

where the on-diagonal blocks are defined as

$\displaystyle\mathbf{M}_{i}$
$\displaystyle=\mathbf{q}_{i}^{\dagger}\mathbf{d}\mathbf{q}_{i}.$ (21)

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

$\displaystyle\mathbf{S}_{1,1}^{(n)}$
$\displaystyle=\mathbf{q}_{1}^{\dagger}\mathbf{d}^{n}\mathbf{q}_{1}$
$\displaystyle=\mathbf{L}^{-1}\boldsymbol{\Sigma}^{(n)}\mathbf{L}^{-1,\dagger},$
(22)

where the subscript indices on $\mathbf{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
$\mathbf{\Sigma}^{(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 ($\mathbf{d}$
or $\mathbf{W}$), as

$\displaystyle\mathbf{S}_{i+1,i}^{(n)}$
$\displaystyle=\mathbf{q}_{i+1}^{\dagger}\mathbf{d}^{n}\mathbf{q}_{i}=\mathbf{C}_{i}^{-1}\left[\mathbf{S}_{i,i}^{(n+1)}-\mathbf{M}_{i}\mathbf{S}_{i,i}^{(n)}-\mathbf{C}_{i-1}^{\dagger}\mathbf{S}_{i-1,i}^{(n)}\right],$
(23) $\displaystyle\mathbf{S}_{i+1,i+1}^{(n)}$
$\displaystyle=\mathbf{q}_{i+1}^{\dagger}\mathbf{d}^{n}\mathbf{q}_{i+1}$
$\displaystyle=\mathbf{C}_{i}^{-1}\left[\mathbf{S}_{i,i}^{(n+2)}+\mathbf{M}_{i}\mathbf{S}_{i,i}^{(n)}\mathbf{M}_{i}+\mathbf{C}_{i-1}^{\dagger}\mathbf{S}_{i-1,i-1}^{(n)}\mathbf{C}_{i-1}-P(\mathbf{S}_{i,i}^{(n+1)}\mathbf{M}_{i})+P(\mathbf{M}_{i}\mathbf{S}_{i,i-1}^{(n)}\mathbf{C}_{i-1})-P(\mathbf{S}_{i,i-1}^{(n+1)}\mathbf{C}_{i-1})\right]\mathbf{C}_{i}^{-1,\dagger},$
(24)

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

$\displaystyle\mathbf{C}_{i}^{2}=$
$\displaystyle\left[\mathbf{S}_{i,i}^{(2)}+\mathbf{M}_{i}^{2}+\mathbf{C}_{i-1}^{\dagger}\mathbf{C}_{i-1}\right.$
$\displaystyle-\left.P(\mathbf{S}_{i,i}^{(1)}\mathbf{M}_{i})-P(\mathbf{S}_{i,i-1}^{(1)}\mathbf{C}_{i-1})\right],$
(25)

and the on-diagonal $\mathbf{M}$ matrices can be found using Eq. 21

$\displaystyle\mathbf{M}_{i}$
$\displaystyle=\mathbf{q}_{i}^{\dagger}\mathbf{d}\mathbf{q}_{i}=\mathbf{S}_{i,i}^{(1)}.$
(26)

These recurrence relations allow the calculation of the on- and off-diagonal
blocks resulting in the block tridiagonal form of the Hamiltonian in Eq. 16.
Despite the apparent complexity of the recurrence relations, this algorithm
contains no step scaling greater than $\mathcal{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[61, 71, 72, 73, 74,
60, 75, 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

$\displaystyle\tilde{\mathbf{H}}$
$\displaystyle=\begin{bmatrix}\mathbf{f}+\boldsymbol{\Sigma}_{\infty}&\tilde{\mathbf{W}}\\\
\tilde{\mathbf{W}}^{\dagger}&\tilde{\mathbf{d}}\\\
\end{bmatrix}=\begin{bmatrix}\mathbf{f}+\boldsymbol{\Sigma}_{\infty}&\tilde{\mathbf{W}}^{<}&\tilde{\mathbf{W}}^{>}\\\
\tilde{\mathbf{W}}^{<,\dagger}&\tilde{\mathbf{d}}^{<}&\mathbf{0}\\\
\tilde{\mathbf{W}}^{>,\dagger}&\mathbf{0}&\tilde{\mathbf{d}}^{>}\\\
\end{bmatrix},$ (27)

where $\tilde{\mathbf{W}}$ are equal to the $\mathbf{L}$ matrix padded by
zeros

$\displaystyle\tilde{\mathbf{W}}^{\lessgtr}=\begin{bmatrix}\mathbf{L}^{\lessgtr}&\mathbf{0}&\cdots&\mathbf{0}\end{bmatrix},$
(28)

and $\tilde{\mathbf{d}}$ are defined by the block tridiagonal elements

$\displaystyle\tilde{\mathbf{d}}^{\lessgtr}$
$\displaystyle=\begin{bmatrix}\mathbf{M}_{1}^{\lessgtr}&\mathbf{C}_{1}^{\lessgtr}&&\mathbf{0}\\\
\mathbf{C}_{1}^{\lessgtr,\dagger}&\mathbf{M}_{2}^{\lessgtr}&\ddots\\\
&\ddots&\ddots&\mathbf{C}_{j-1}^{\lessgtr}\\\
\mathbf{0}&&\mathbf{C}_{j-1}^{\lessgtr,\dagger}&\mathbf{M}_{j}^{\lessgtr}\end{bmatrix}.$
(29)

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 $\tilde{\mathbf{d}}$ and appropriately rotating
$\tilde{\mathbf{W}}$ into this basis. The eigenvalues of $\tilde{\mathbf{H}}$
are moment-conserving approximations to those of the exact upfolded
Hamiltonian, and the corresponding eigenvectors $\mathbf{u}$ can be
transformed into Dyson orbitals via $\mathbf{L}\mathbf{P}\mathbf{u}$, where
$\mathbf{P}$ is a projection into the physical space, and the $\mathbf{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_{\mathrm{iter}}$ is used; in
this notation the $n_{\mathrm{iter}}=0$ calculation corresponds to the
inclusion of only a single on-diagonal block $\mathbf{M}_{1}$, with modified
couplings $\mathbf{L}$ to the physical space. As such, in this notation the
number of conserved moments equals $2n_{\mathrm{iter}}+2$, i.e. up to and
including the $2n_{\mathrm{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_{\mathrm{iter}}$ applications of this algorithm to both
the lesser and greater self-energy sectors, this results in
$N(2n_{\mathrm{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_{\mathrm{iter}}\sim
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 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 [78, 79, 80, 5, 63]. This
effectively ensures that all information required to build the 2-point RPA dd-
response is contained in the first two spectral moments, analogous to how all
the information on the density of states in mean-field theory (quadratic in
Fermionic operators) is contained in the first two Green’s function moments
(i.e. the one-body density matrix and Fock matrix).

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

$\begin{bmatrix}\mathbf{A}&\mathbf{B}&\\\
\mathbf{-B}&\mathbf{-A}&\end{bmatrix}\begin{bmatrix}\mathbf{X}&\mathbf{Y}\\\
\mathbf{Y}&\mathbf{X}\end{bmatrix}=\begin{bmatrix}\mathbf{X}&\mathbf{Y}\\\
\mathbf{Y}&\mathbf{X}\end{bmatrix}\begin{bmatrix}\mathbf{\Omega}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{\Omega}\end{bmatrix},$ (30)

where the left and right eigenvectors form the biorthogonal set as

$\begin{bmatrix}\mathbf{X}&\mathbf{Y}\\\
\mathbf{Y}&\mathbf{X}\end{bmatrix}^{T}\begin{bmatrix}\mathbf{X}&\mathbf{Y}\\\
\mathbf{-Y}&\mathbf{-X}\end{bmatrix}=\begin{bmatrix}\mathbf{I}&\mathbf{0}\\\
\mathbf{0}&\mathbf{-I}\end{bmatrix}.$ (31)

This biorthogonality ensures an inverse relationship between
$(\mathbf{X}+\mathbf{Y})$ and $(\mathbf{X}-\mathbf{Y})^{T}$, as

$(\mathbf{X}+\mathbf{Y})(\mathbf{X}-\mathbf{Y})^{T}=(\mathbf{X}+\mathbf{Y})^{T}(\mathbf{X}-\mathbf{Y})=\mathbf{I}.$
(32)

The $\mathbf{A}$ and $\mathbf{B}$ matrices are defined as

$\displaystyle A_{ia,jb}$
$\displaystyle=\left(\epsilon_{a}-\epsilon_{i}\right)\delta_{ij}\delta_{ab}+\mathcal{K}_{ia,bj}$
(33) $\displaystyle B_{ia,jb}$ $\displaystyle=\mathcal{K}_{ia,jb}.$ (34)

Here, $\mathcal{K}$ is an interaction kernel which couples particle-hole
excitations and de-excitations. In the traditional RPA (without second-order
exchange), this coupling is taken to be the same for excitations, de-
excitations and their coupling, given by the static, bare Coulomb interaction,
$\mathcal{K}_{ia,jb}=(ia|jb)=\mathcal{K}_{ia,bj}$. Hole and particle orbital
energies are given by $\epsilon_{i}$ and $\epsilon_{a}$ respectively, defining
the irreducible polarizability of the system from the reference state in
$\mathbf{A}$. Upon diagonalization, the eigenvectors defined by $X_{ia,\nu}$
and $Y_{ia,\nu}$ define the coefficients of the RPA excitations in the
particle-hole and hole-particle basis, with energies $\Omega_{\nu}$, with
$\mathbf{\Omega}$ 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

$\mathbf{\chi}(\omega)=\begin{bmatrix}\mathbf{X}&\mathbf{Y}\\\
\mathbf{Y}&\mathbf{X}\end{bmatrix}\begin{bmatrix}\omega\mathbf{I}-\mathbf{\Omega}\hskip
10.0pt&\mathbf{0}\\\ \mathbf{0}\hskip
10.0pt&-\omega\mathbf{I}-\mathbf{\Omega}\end{bmatrix}^{-1}\begin{bmatrix}\mathbf{X}&\mathbf{Y}\\\
\mathbf{Y}&\mathbf{X}\end{bmatrix}^{T}.$ (35)

Note that this matrix is formally equivalent to the alternative directly
dynamical construction from
$\mathbf{\chi}(\omega)=(\mathbf{P}(\omega)^{-1}-\mathcal{K})^{-1}$, where
$\mathbf{P}(\omega)$ 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

$\eta(\omega)=(\mathbf{X}+\mathbf{Y})(\omega\mathbf{I}-\mathbf{\Omega})^{-1}(\mathbf{X}+\mathbf{Y})^{T},$
(36)

which sums contributions from particle-hole and hole-particle fluctuations
together, and from which optical properties such as dynamic polarizabilities
can be computed [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

$\eta^{(n)}_{ia,jb}=-\frac{1}{\pi}\int_{0}^{\infty}\textrm{Im}[\eta_{ia,jb}(\omega)]\omega^{n}d\omega.$
(37)

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

$\mathbf{\eta}^{(n)}=(\mathbf{X}+\mathbf{Y})\mathbf{\Omega}^{n}(\mathbf{X}+\mathbf{Y})^{T},$
(38)

and constitutes a central object of interest in this work, required for the GW
self-energy moment construction of Eqs. 9-10. 111 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
$(\mathbf{X+Y})^{-1}=(\mathbf{X-Y})^{T}$, that is
$(\mathbf{\eta}^{(n)})^{-1}=(\mathbf{X}-\mathbf{Y})\mathbf{\Omega}^{-n}(\mathbf{X}-\mathbf{Y})^{T}$.

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

$\displaystyle\mathbf{\eta}^{(1)}$
$\displaystyle=(\mathbf{X}+\mathbf{Y})\mathbf{\Omega}(\mathbf{X}+\mathbf{Y})^{T}=(\mathbf{A}-\mathbf{B}),$
(39)
$\displaystyle=\mathbf{\eta}^{(0)}(\mathbf{A}+\mathbf{B})\mathbf{\eta}^{(0)},$
(40)

noting that
$(\mathbf{A}+\mathbf{B})=(\mathbf{X}-\mathbf{Y})\mathbf{\Omega}(\mathbf{X}-\mathbf{Y})^{T}$.
By taking the sum and difference of the two Casida equations of Eq. 30, we
also find

$\displaystyle\left(\mathbf{A}+\mathbf{B}\right)\left(\mathbf{X}+\mathbf{Y}\right)$
$\displaystyle=\left(\mathbf{X}-\mathbf{Y}\right)\mathbf{\Omega}$ (41)
$\displaystyle\left(\mathbf{A}-\mathbf{B}\right)\left(\mathbf{X}-\mathbf{Y}\right)$
$\displaystyle=\left(\mathbf{X}+\mathbf{Y}\right)\mathbf{\Omega},$ (42)

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

$\left(\mathbf{A}-\mathbf{B}\right)\left(\mathbf{A}+\mathbf{B}\right)\left(\mathbf{X}+\mathbf{Y}\right)=\left(\mathbf{X}+\mathbf{Y}\right)\mathbf{\Omega}^{2},$
(43)

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

$(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})\mathbf{\eta}^{(0)}=\mathbf{\eta}^{(2)}.$
(44)

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

$\displaystyle\mathbf{\eta}^{(m)}$
$\displaystyle=(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})\mathbf{\eta}^{(m-2)}$
(45)
$\displaystyle=[\mathbf{\eta}^{(0)}(\mathbf{A}+\mathbf{B})]^{m}\mathbf{\eta}^{(0)}.$
(46)

While these are important relations in themselves, they also illustrate that
all the RPA excitations and weights in the dd-response of Eq. 36 are
implicitly accessible without requiring the explicit solution to the Casida
equation ($\mathbf{X}$, $\mathbf{Y}$ and $\mathbf{\Omega}$ matrices). This
reformulation just requires knowledge of $\mathbf{\eta}^{(0)}$ as the central
variable (which can be defined independently of the original equations via Eq.
40), the $\mathbf{A}$ and $\mathbf{B}$ matrices defining the system and their
interactions, and the recursive relations of Eq. 46. As an aside, the
Tamm–Dancoff approximation (TDA) sets $\mathbf{B}=0$, which dramatically
simplifies the resulting expressions due to the lack of correlation in the
ground state, with $\eta^{(0)}=\mathbf{I}$, and $\eta^{(n)}=\mathbf{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

$E_{\mathrm{corr}}^{\mathrm{RPA}}=\frac{1}{2}\mathrm{Tr}[\mathbf{\eta}^{(0)}(\mathbf{A}+\mathbf{B})-\mathbf{A}].$
(47)

Related expressions for the RPA correlation energy can be found found in Eqs.
39-40 of Ref. 86 in terms of the quantities $\mathbf{\eta}^{(0)}$,
$\mathbf{A+B}$ and $\mathbf{A-B}$ (where there, these quantities are referred
to as $\mathbf{Q}^{\text{dRPA}}_{\text{I}}$, $\mathbf{\varepsilon}$ and
$\mathbf{\varepsilon}+2\mathbf{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

$\displaystyle\text{Tr}\left[\mathbf{\eta}^{(0)}(\mathbf{A+B})\right]$ (48)
$\displaystyle=$
$\displaystyle\text{Tr}\left[(\mathbf{\eta}^{(0)})^{-1}(\mathbf{A-B})\right]$
(49) $\displaystyle=$
$\displaystyle\text{Tr}\left[\left((\mathbf{A-B})^{\frac{1}{2}}(\mathbf{A+B})(\mathbf{A-B})^{\frac{1}{2}}\right)^{\frac{1}{2}}\right]=\text{Tr}\left[\mathbf{\Omega}\right].$
(50)

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

##  IV Efficient evaluation of self-energy and density response moments

Given the recasting of the RPA dd-response in Sec. III in terms of its lowest
order moment (Eq. 40) and recursion to access the higher moments via Eq. 46,
we now consider the efficient $\mathcal{O}[N^{4}]$ evaluation of these
quantities which are central to the approach in this work, thus avoiding their
formal $\mathcal{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

$(ia|jb)\simeq\sum_{P}V_{ia,P}V_{jb,P}=\mathbf{V}\mathbf{V}^{T},$ (51)

where we use $P,Q,\dots$ to index elements of this auxiliary (RI) basis, whose
dimension $N_{\mathrm{aux}}$ scales $\mathcal{O}[N]$ with system size. We
define an intermediate quantity

${\tilde{\eta}}^{(n)}_{ia,P}=\sum_{jb}\eta^{(n)}_{ia,jb}V_{jb,P}\quad.$ (52)

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

$\Sigma_{pq}^{(n,>)}=\sum_{t=0}^{n}\binom{n}{t}\left(\epsilon_{c}^{n-t}V_{pc,Q}\left(V_{qc,P}\left({\tilde{\eta}}_{ia,P}^{(t)}V_{ia,Q}\right)\right)\right),$
(53)

where the brackets indicate the order of contractions in order to preserve
$\mathcal{O}[N^{4}]$ scaling, and einstein summation is implied. The lesser
self-energy moment of Eq. 9 can be recast in an analogous fashion. 222We 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 $\mathcal{O}[N]$ (in
this case the auxiliary basis).

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 ${\tilde{\eta}}_{ia,P}^{(0)}$
and ${\tilde{\eta}}_{ia,P}^{(1)}$, via use of the recursive relationship
between the moments as given by Eq. 45, as for even moment orders,

$\tilde{\mathbf{\eta}}^{(n)}=[(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})]^{n/2}\tilde{\mathbf{\eta}}^{(0)}$
(54)

and for odd moment orders,

$\tilde{\mathbf{\eta}}^{(n)}=[(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})]^{(n-1)/2}\tilde{\mathbf{\eta}}^{(1)}$
(55)

where we have omitted explicit indices for brevity. Ensuring that an
$\mathcal{O}[N^{4}]$ scaling is retained in this recursion relies on
$(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{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-34),

$(\mathbf{A}-\mathbf{B})_{ia,jb}=(\epsilon_{a}-\epsilon_{i})\delta_{ij}\delta_{ab}=\mathbf{D}$
(56)

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

$(\mathbf{A}+\mathbf{B})=\mathbf{D}+2\mathcal{K}=\mathbf{D}+2\sum_{P}V_{ia,P}V_{jb,P}.$
(57)

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

$(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})=\mathbf{D}^{2}+\mathbf{S}_{L}\mathbf{S}_{R}^{T},$
(58)

where for the RPA, $\mathbf{D}$ is defined by Eq. 56,
$\mathbf{S}_{L}=\mathbf{D}\mathbf{V}$ and $\mathbf{S}_{R}=2\mathbf{V}$. Future
work will explore other analogous approaches where
$(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{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-55, a fixed number
of dd-moments of the form of Eq. 52 can be found in $\mathcal{O}[N^{4}]$ time,
provided the original $\tilde{\eta}^{(0)}$ and $\tilde{\eta}^{(1)}$ values are
known.

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

$\tilde{\eta}^{(1)}=\mathbf{D}\mathbf{V}=\mathbf{S}_{L}.$ (59)

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

$(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})=\left(\eta^{(0)}(\mathbf{A}+\mathbf{B})\right)^{2},$
(60)

from which we can find an expression for $\eta^{(0)}$ as

$\eta^{(0)}=[(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})]^{\frac{1}{2}}(\mathbf{A}+\mathbf{B})^{-1}.$
(61)

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

$\tilde{\eta}^{(0)}=(\mathbf{D}^{2}+\mathbf{S}_{L}\mathbf{S}_{R}^{T})^{\frac{1}{2}}(\mathbf{D}+2\mathbf{V}\mathbf{V}^{T})^{-1}\mathbf{V}.$
(62)

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

$\displaystyle\mathbf{T}$
$\displaystyle=(\mathbf{D}+2\mathbf{V}\mathbf{V}^{T})^{-1}\mathbf{V}$ (63)
$\displaystyle=\mathbf{D}^{-1}\mathbf{V}-2\mathbf{D}^{-1}\mathbf{V}(\mathbf{I}+2\mathbf{V}^{T}\mathbf{D}^{-1}\mathbf{V})^{-1}\mathbf{V}^{T}\mathbf{D}^{-1}\mathbf{V}.$
(64)

This now only requires the inversion of the diagonal matrix, $\mathbf{D}$, and
a matrix of dimension $N_{\mathrm{aux}}$, with the overall $ov\times
N_{\mathrm{aux}}$ matrix able to be constructed in
$\mathcal{O}[N_{\mathrm{aux}}^{3}+N_{\mathrm{aux}}^{2}ov]$ time.

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

$\mathbf{M}^{\frac{1}{2}}=\frac{1}{\pi}\int_{-\infty}^{\infty}\left(\mathbf{I}-z^{2}(\mathbf{M}+z^{2}\mathbf{I})^{-1}\right)dz.$
(65)

From Eq. 62, this results in

$\displaystyle\tilde{\mathbf{\eta}}^{(0)}$
$\displaystyle=(\mathbf{D}^{2}+\mathbf{S}_{L}\mathbf{S}_{R}^{T})^{\frac{1}{2}}\mathbf{T}$
(66)
$\displaystyle=\frac{1}{\pi}\int_{-\infty}^{\infty}\left(\mathbf{I}-z^{2}(\mathbf{D}^{2}+\mathbf{S}_{L}\mathbf{S}_{R}^{T}+z^{2}\mathbf{I})^{-1}\right)\mathbf{T}dz.$
(67)

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

$\displaystyle\mathbf{F}(z)$
$\displaystyle=(\mathbf{D}^{2}+z^{2}\mathbf{I})^{-1}$ (68)
$\displaystyle\mathbf{Q}(z)$
$\displaystyle=\mathbf{S}_{R}^{T}\mathbf{F}(z)\mathbf{S}_{L},$ (69)

where $\mathbf{F}(z)$ is a diagonal matrix in the $ov$ space, and
$\mathbf{Q}(z)$ is a $N_{\mathrm{aux}}\times N_{\mathrm{aux}}$ matrix which
can be constructed in $\mathcal{O}[N_{\mathrm{aux}}^{2}ov]$ time. This casts
Eq. 67 into the form

$\tilde{\eta}^{(0)}=\frac{1}{\pi}\int_{-\infty}^{\infty}\left[\mathbf{I}-z^{2}\mathbf{F}(z)\left(\mathbf{I}-\mathbf{S}_{L}(\mathbf{I}+\mathbf{Q}(z))^{-1}\mathbf{S}_{R}^{T}\mathbf{F}(z)\right)\right]\mathbf{T}dz.$
(70)

For each value of the integration variable $z$, the integrand is a matrix of
size $ov\times N_{\mathrm{aux}}$, which can be constructed in
$\mathcal{O}[N^{4}]$ scaling, rendering it efficient for numerical quadrature.
Along with the results of Eqs. 59, 58, 54 and 55, 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
$\mathcal{O}[N^{4}]$ scaling (and $\mathcal{O}[N^{3}]$ memory).

However, manipulations of the resulting integrand and choice of quadrature
points can further improve the efficiency of their construction by ensuring a
faster decay of the integrand and separating components which can be
analytically integrated. This derivation is given in Appendix A, and results
in a final $\mathcal{O}[N_{\mathrm{aux}}^{2}ov+N_{\mathrm{aux}}^{3}]$
expression to evaluate for the zeroth-order dd-moment as

$\displaystyle\tilde{\eta}^{(0)}$ $\displaystyle=\mathbf{D}\mathbf{T}$
$\displaystyle+\int_{0}^{\infty}e^{-t\mathbf{D}}\mathbf{S}_{L}\mathbf{S}_{R}^{T}e^{-t\mathbf{D}}\mathbf{T}dt$
$\displaystyle+\frac{1}{\pi}\int_{-\infty}^{\infty}z^{2}\mathbf{F}(z)\mathbf{S}_{L}\left((\mathbf{I}+\mathbf{Q}(z))^{-1}-\mathbf{I}\right)\mathbf{S}_{R}^{T}\mathbf{F}(z)\mathbf{T}dz.$
(71)

The first numerical integral in Eq. 71 (where the integrand decays
exponentially) is computed via Gauss-Laguerre quadrature, while the second
(where the integrand decays as $\mathcal{O}[z^{-4}]$) is evaluated via
Clenshaw–Curtis quadrature. A comparison of the decay of the original and
refined integrands is shown in the inset to Fig. 1.

Figure 1: Exponential convergence of the numerical integration (NI) error for
$\tilde{\mathbf{\eta}^{(0)}}$ in Eq. 71 with respect to integration points,
for the singlet oxygen dimer in a cc-pVTZ basis at equilibrium (1.207 Å) bond
length. Also included are the two error estimates of the true NI error, used
to check convergence and estimate the required number of points (see App. B,
Eq. B10 for the ‘Nested Fit’ and Eq. B5 for the ‘Lower Bound’ definitions).
Inset: The originally derived integrand (Eq. 70), and the form optimized for
efficient NI given in Eq. 71 and derived in App. A, showing the faster decay.

The scaling of the grid spacing of both numerical quadratures is optimised to
ensure exact integration of the trace of a diagonal approximation for the
integrand, analogous to the grid optimization discussed in Ref. 56. For
numerical robustness, we optimize the quadrature for evaluating
$\textrm{Tr}[(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})]^{\frac{1}{2}}$,
rather than the full integrand. We write a diagonal approximation to
$(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})$ as

$\mathbf{M}^{D}=\mathbf{D}^{2}+\mathbf{S}_{L}^{D}(\mathbf{S}_{R}^{D})^{T},$
(72)

with

$\displaystyle(S_{L}^{D})_{ia,R}$
$\displaystyle=(S_{L})_{ia,Q}(S_{R})_{ia,Q}\delta_{R,ia}$ (73)
$\displaystyle(S_{R}^{D})_{jb,R}$ $\displaystyle=\delta_{R,jb}.$ (74)

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 $\mathcal{O}[ov]$ computational time. We then
seek to ensure the trace of the difference between the exact and numerically
integrated estimate of ${\mathbf{M}^{D}}^{\frac{1}{2}}$ vanishes. This is
achieved for both integrals in Eq. (71), with the analytic form for the first
integral given by
$\frac{1}{2}\text{diag}(\mathbf{D}^{-1}\mathbf{S}_{L}\mathbf{S}_{R}^{T})=I^{\text{offset}}$
and the second numerical integral as
${\mathbf{M}^{D}}^{\frac{1}{2}}-\mathbf{D}-I^{\text{offset}}=I^{\text{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

$\displaystyle O^{\text{offset}}(a)$
$\displaystyle=I^{\text{offset}}-\frac{1}{\pi}\sum_{i}^{n_{p}}aw_{i}\text{Tr}(e^{-2\mathbf{D}az_{i}}\mathbf{S}_{L}\mathbf{S}_{R}^{T})$
(75) $\displaystyle O^{\text{int}}(a)$
$\displaystyle=I^{\text{int}}-\frac{1}{\pi}\sum_{i}^{n_{p}}aw_{i}\text{Tr}(I^{\text{D}}(az_{i}))$
(76) $\displaystyle I^{\text{D}}(z)$
$\displaystyle=z^{2}\mathbf{F}^{D}(z)\mathbf{S}^{D}_{L}\left((\mathbf{I}+\mathbf{Q}^{D}(z))^{-1}-\mathbf{I}\right)(\mathbf{S}^{D}_{R})^{T}\mathbf{F}^{D}(z),$
(77)

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 $\mathcal{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 lowest-scaling $GW$
approaches [92, 24, 93, 42]. 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

$(ia|jb)\simeq\sum_{P,Q}X_{iP}X_{aP}Z_{PQ}X_{jQ}X_{bQ},$ (78)

This factorizes the orbital product into separate terms, and is also known as
tensor hypercontraction or CP decomposition, used for various recent low-
scaling formulations of quantum chemical methods, where the dimension of the
$P$ and $Q$ space rigorously grow linearly with system size[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
$\mathcal{O}[N^{3}]$ 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 $\mathcal{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}\eta^{(0)}_{ia,jb}X_{jR}X_{bR}Z_{PR}$. Inspired by the low-
scaling 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 $\mathbf{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

$\displaystyle{F}_{ia,ia}(z)$
$\displaystyle=D_{ia,ia}\int_{p=0}^{\infty}\frac{\sin{zp}}{z}e^{-pD_{ia,ia}}dp,$
(79)

which allows the key matrix $\mathbf{Q}(z)$ of Eq. 69 to be obtained as

$\displaystyle Q_{PS}(z)$
$\displaystyle=2\int_{p=0}^{\infty}\frac{\sin{zp}}{z}Y_{PQ}(A_{QR}(p)-B_{QR}(p))Y_{SR}dp$
(80)

where both intermediates

$\displaystyle A_{QR}(p)$
$\displaystyle=X_{iQ}X_{aQ}\epsilon_{a}e^{-\epsilon_{a}p}e^{\epsilon_{i}p}X_{iR}X_{aR}$
(81) $\displaystyle B_{QR}(p)$
$\displaystyle=X_{iQ}X_{aQ}\epsilon_{i}e^{-\epsilon_{a}p}e^{\epsilon_{i}p}X_{iR}X_{aR}$
(82)

can be evaluated in $\mathcal{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 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 forthcoming work,
with numerical results in the rest of this work employing the quartic scaling
algorithm described in Sec. IV.

##  VI Results

###  VI.1 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, 103, 104, 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 reference [11]. In Fig. 2 we
consider the convergence of the first ionization potential (IP) of singlet O2
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 $\mathcal{O}[N^{6}]$), and one which
performs an analytic continuation (AC) of the imaginary frequency self-energy
to real-frequencies via a fit to Padé approximants in order to perform the
convolution (denoted ‘AC QP-$G_{0}W_{0}$’, which scales as
$\mathcal{O}[N^{4}]$) [44, 106, 105]. 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.

Figure 2: Convergence of the first $G_{0}W_{0}$ IP of singlet O2 in a cc-pVDZ
basis with respect to the number of conserved self-energy moments. Also shown
is the same quasiparticle state computed from traditional $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 $\Sigma$ moment’). All approaches find
an IP of $8.49$ eV to within $10$ meV, with the difference between full-
frequency and AC approaches $5$ meV, and the relaxation of diagonal
approximation also accounting for small $\sim 5$ meV differences, with the
reference Hartree–Fock IP for comparison being 9.73 eV.

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

Having demonstrated correctness compared to a high-scaling exact frequency
implementation, we can compare our results to AC-$G_{0}W_{0}$ across a larger
test set to consider the moment truncation convergence. We use the established
‘$GW100$’ benchmark test set, where many $GW$ (and other excited state method)
implementations have been rigorously compared [12, 107, 108, 47]. This
benchmark set contains 102 diverse systems, with the IP of the molecules in
the set ranging from $\sim 4-25$ eV, and featuring molecules with bound metal
atoms (including metal diatomics and clusters), strongly ionic bonding, and
molecules with a strongly 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.

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

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

Figure 4: Mean absolute errors (eV) for the IP ($x$-axis) and EA ($y$-axis)
across the $GW100$ benchmark set in a def2-TZVPP basis compared to accurate
CCSD(T) values. AC-$G_{0}W_{0}$ values, as well as moment-truncated results
are shown, with the number in each data point marker giving the number of
conserved moments. Extrapolation of individual data points from the
$9^{\mathrm{th}}$\- and $11^{\mathrm{th}}$-order conserved moment values are
also provided, denoted by the truncation label ‘$\infty$’.

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

###  VI.2 Full frequency spectra

Figure 5:  Spectral functions for guanine in a def2-TZVPP basis set, for HF,
AC-$G_{0}W_{0}$, and the moment-conserving $G_{0}W_{0}$ approach with zero to
five block Lanczos iterations, thereby conserving up to the
$11^{\mathrm{th}}$-order self-energy moment. The values indicated in the
spectra indicate the Wasserstein metric taken with 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.

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 $\mathcal{O}[N^{5}]$ 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 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 $\sim
7^{\mathrm{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 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 the next
section.

###  VI.3 Multiple solutions and additional spectral features

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

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, 110, 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
$Z_{p}=(1-\frac{\partial\Sigma_{pp}(\epsilon_{p})}{\partial\omega})^{-1}$,
which indicates the approximate weight in the quasi-particle solutions.

Since the moment-conserving $GW$ approach obtains all poles in one step
(including satellites and low-weighted features), all the excitations can be
characterized by their quasiparticle weight, and either selected as specific
excitation energies, or all excitations included in the spectrum to ensure
smooth changes with molecular geometry. Points of discontinuity will therefore
manifest as the presence of multiple lower-weighted solutions at a given
energy, giving a smooth change to a broadened feature in the spectrum near
self-energy poles. To demonstrate this, we consider the same simple system as
Ref. 111, observing the $GW$ quasiparticle energies as a function of the
inter-nuclear distance in the H2 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^{\mathrm{st}}$\- and $11^{\mathrm{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^{\mathrm{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 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 the 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 $\mathcal{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 $GW100$
test set, we find rapid convergence to established $GW$ methodology results
for both state-specific and full spectral properties, with errors due to the
incompleteness of the moment expansion many times smaller than the inherent
accuracy of the method. The formulation follows relatively closely from
previous low-scaling approaches to RPA correlation energies, enabling these
codes to be adapted to low-scaling $GW$ methods with relatively little effort.

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

## Acknowledgements

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 programme under grant agreement No.
759063. We are grateful to the UK Materials and Molecular Modelling Hub for
computational resources, which is partially funded by EPSRC (EP/P020194/1 and
EP/T022213/1). Further computational resources were awarded under the embedded
CSE programme of the ARCHER2 UK National Supercomputing Service
(http://www.archer2.ac.uk).

## References

  * Cohen, Mori-Sánchez, and Yang [2012] A. J. Cohen, P. Mori-Sánchez, and W. Yang, ‘‘Challenges for Density Functional Theory,’’ Chem. Rev. 112, 289--320 (2012).
  * Onida, Reining, and Rubio [2002] G. Onida, L. Reining, and A. Rubio, ‘‘Electronic excitations: density-functional versus many-body Green’s-function approaches,’’ Rev. Mod. Phys. 74, 601--659 (2002).
  * Hybertsen and Louie [1985] M. S. Hybertsen and S. G. Louie, ‘‘First-Principles Theory of Quasiparticles: Calculation of Band Gaps in Semiconductors and Insulators,’’ Phys. Rev. Lett. 55, 1418--1421 (1985).
  * Aryasetiawan and Gunnarsson [1998] F. Aryasetiawan and O. Gunnarsson, ‘‘The GW method,’’ Rep. Prog. Phys. 61, 237--312 (1998).
  * Hedin [1999] L. Hedin, ‘‘On correlation effects in electron spectroscopies and the $GW$ approximation,’’ J. Phys. Condens. Matter 11, R489--R528 (1999).
  * Aulbur, Jönsson, and Wilkins [2000] W. G. Aulbur, L. Jönsson, and J. W. Wilkins, ‘‘Quasiparticle Calculations in Solids,’’ in  _Solid State Physics_ (Elsevier, 2000) pp. 1--218.
  * Friedrich and Schindlmayr [2006] C. Friedrich and A. Schindlmayr, ‘‘Many-Body Perturbation Theory : The GW Approximation,’’ Computing 31, 335--355 (2006).
  * Kutepov, Savrasov, and Kotliar [2009] A. Kutepov, S. Y. Savrasov, and G. Kotliar, ‘‘Ground-state properties of simple elements from GW calculations,’’ Phys. Rev. B 80, 1--4 (2009).
  * Ke [2011] S. H. Ke, ‘‘All-electron GW methods implemented in molecular orbital space: Ionization energy and electron affinity of conjugated molecules,’’ Phys. Rev. B 84, 1--6 (2011).
  * Bruneval [2012] F. Bruneval, ‘‘Ionization energy of atoms obtained from GW self-energy or from random phase approximation total energies,’’ J. Chem. Phys. 136 (2012).
  * van Setten, Weigend, and Evers [2013] M. J. van Setten, F. Weigend, and F. Evers, ‘‘The GW-method for quantum chemistry applications: Theory and implementation,’’ J. Chem. Theory Comput. 9, 232--246 (2013).
  * van Setten  _et al._ [2015] M. J. van Setten, F. Caruso, S. Sharifzadeh, X. Ren, M. Scheffler, F. Liu, J. Lischner, L. Lin, J. R. Deslippe, S. G. Louie, C. Yang, F. Weigend, J. B. Neaton, F. Evers, and P. Rinke, ‘‘GW100: Benchmarking G0W0 for Molecular Systems,’’ J. Chem. Theory Comput. 11, 5665--5687 (2015).
  * Reining [2017] L. Reining, ‘‘The GW approximation: content, successes and limitations,’’ Wiley Interdiscip. Rev. Comput. Mol. Sci. 8 (2017).
  * Golze, Dvorak, and Rinke [2019] D. Golze, M. Dvorak, and P. Rinke, ‘‘The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy,’’ Front. Chem. 7 (2019).
  * Hedin [1965] L. Hedin, ‘‘New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem,’’ Phys. Rev. 139, A796--A823 (1965).
  * Hedin and Lundqvist [1970] L. Hedin and S. Lundqvist, ‘‘Effects of Electron-Electron and Electron-Phonon Interactions on the One-Electron States of Solids,’’ in  _Solid State Physics_ (Elsevier, 1970) pp. 1--181.
  * Langreth and Perdew [1977] D. C. Langreth and J. P. Perdew, ‘‘Exchange-correlation energy of a metallic surface: Wave-vector analysis,’’ Phys. Rev. B 15, 2884--2901 (1977).
  * Chong [1995] D. P. Chong,  _Recent Advances in Density Functional Methods_ (World Scientific, 1995).
  * Heßelmann and Görling [2010] A. Heßelmann and A. Görling, ‘‘Random phase approximation correlation energies with exact Kohn-Sham exchange,’’ Mol. Phys. 108, 359--372 (2010).
  * Ren  _et al._ [2012a] X. Ren, P. Rinke, C. Joas, and M. Scheffler, ‘‘Random-phase approximation and its applications in computational chemistry and materials science,’’ J. Mater. Sci. 47, 7447--7471 (2012a).
  * Stan, Dahlen, and van Leeuwen [2009] A. Stan, N. E. Dahlen, and R. van Leeuwen, ‘‘Levels of self-consistency in the GW approximation,’’ J. Chem. Phys. 130, 114105 (2009).
  * Foerster, Koval, and Snchez-Portal [2011] D. Foerster, P. Koval, and D. Snchez-Portal, ‘‘An O(N3) implementation of Hedins GW approximation for molecules,’’ J. Chem. Phys. 135, 1--19 (2011).
  * Bruneval and Marques [2013] F. Bruneval and M. A. L. Marques, ‘‘Benchmarking the Starting Points of the GW Approximation for Molecules,’’ J. Chem. Theory Comput. 9, 324--329 (2013).
  * Liu  _et al._ [2016] P. Liu, M. Kaltak, J. Klimeš, and G. Kresse, ‘‘Cubic scaling $GW$: Towards fast quasiparticle calculations,’’ Phys. Rev. B 94, 165109 (2016).
  * Knight  _et al._ [2016] J. W. Knight, X. Wang, L. Gallandi, O. Dolgounitcheva, X. Ren, J. V. Ortiz, P. Rinke, T. Körzdörfer, and N. Marom, ‘‘Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules III: A Benchmark of GW Methods,’’ J. Chem. Theory Comput. 12, 615--626 (2016).
  * Maggio  _et al._ [2017] E. Maggio, P. Liu, M. J. van Setten, and G. Kresse, ‘‘$GW_{1}00$: A Plane Wave Perspective for Small Molecules,’’ J. Chem. Theory Comput. 13, 635--648 (2017).
  * Bruneval, Dattani, and van Setten [2021] F. Bruneval, N. Dattani, and M. J. van Setten, ‘‘The GW Miracle in Many-Body Perturbation Theory for the Ionization Potential of Molecules,’’ Front. Chem. 9 (2021).
  * von Barth and Holm [1996] U. von Barth and B. Holm, ‘‘Self-consistent $\mathit{GW}_{0}$ results for the electron gas: Fixed screened potential $\mathit{W}_{0}$ within the random-phase approximation,’’ Phys. Rev. B 54, 8411--8419 (1996).
  * Holm and von Barth [1998] B. Holm and U. von Barth, ‘‘Fully self-consistent $\mathrm{GW}$ self-energy of the electron gas,’’ Phys. Rev. B 57, 2108--2117 (1998).
  * Schöne and Eguiluz [1998] W.-D. Schöne and A. G. Eguiluz, ‘‘Self-Consistent Calculations of Quasiparticle States in Metals and Semiconductors,’’ Phys. Rev. Lett. 81, 1662--1665 (1998).
  * García-González and Godby [2001] P. García-González and R. W. Godby, ‘‘Self-consistent calculation of total energies of the electron gas using many-body perturbation theory,’’ Phys. Rev. B 63 (2001).
  * Faleev, van Schilfgaarde, and Kotani [2004] S. V. Faleev, M. van Schilfgaarde, and T. Kotani, ‘‘All-electron self-consistent GW approximation: Application to Si, MnO, and NiO,’’ Phys. Rev. Lett. 93, 12--15 (2004).
  * van Schilfgaarde, Kotani, and Faleev [2006] M. van Schilfgaarde, T. Kotani, and S. Faleev, ‘‘Quasiparticle Self-Consistent $GW$ Theory,’’ Phys. Rev. Lett. 96, 226402 (2006).
  * Stan, Dahlen, and van Leeuwen [2006] A. Stan, N. E. Dahlen, and R. van Leeuwen, ‘‘Fully self-consistent GW calculations for atoms and molecules,’’ Europhys. Lett. 76, 298--304 (2006).
  * Kotani, van Schilfgaarde, and Faleev [2007] T. Kotani, M. van Schilfgaarde, and S. V. Faleev, ‘‘Quasiparticle self-consistent GW method: A basis for the independent-particle approximation,’’ Phys. Rev. B 76, 1--24 (2007).
  * Shishkin and Kresse [2007] M. Shishkin and G. Kresse, ‘‘Self-consistent $\mathit{GW}$ calculations for semiconductors and insulators,’’ Phys. Rev. B 75 (2007).
  * Caruso  _et al._ [2012] F. Caruso, P. Rinke, X. Ren, M. Scheffler, and A. Rubio, ‘‘Unified description of ground and excited states of finite systems: The self-consistent $GW$ approach,’’ Phys. Rev. B 86, 81102 (2012).
  * Bruneval  _et al._ [2016] F. Bruneval, T. Rangel, S. M. Hamed, M. Shao, C. Yang, and J. B. Neaton, ‘‘MOLGW 1: Many-body perturbation theory software for atoms, molecules, and clusters,’’ Comput. Phys. Commun. 208, 149--161 (2016).
  * Kaplan  _et al._ [2016] F. Kaplan, M. E. Harding, C. Seiler, F. Weigend, F. Evers, and M. J. van Setten, ‘‘Quasi-Particle Self-Consistent GW for Molecules,’’ J. Chem. Theory Comput. 12, 2528--2541 (2016).
  * Jin, Su, and Yang [2019] Y. Jin, N. Q. Su, and W. Yang, ‘‘Renormalized Singles Green’s Function for Quasi-Particle Calculations beyond the G0W0 Approximation,’’ J. Phys. Chem. Lett. 10, 447--452 (2019).
  * Duchemin and Blase [2020] I. Duchemin and X. Blase, ‘‘Robust Analytic-Continuation Approach to Many-Body GW Calculations,’’ J. Chem. Theory Comput. 16, 1742--1756 (2020).
  * Duchemin and Blase [2021] I. Duchemin and X. Blase, ‘‘Cubic-Scaling All-Electron GW Calculations with a Separable Density-Fitting Space-Time Approach,’’ J. Chem. Theory Comput. 17, 2383--2393 (2021).
  * Yeh  _et al._ [2022] C.-N. Yeh, S. Iskakov, D. Zgid, and E. Gull, ‘‘Fully self-consistent finite-temperature $GW$ in Gaussian Bloch orbitals for solids,’’ Phys. Rev. B 106, 235104 (2022).
  * Ren  _et al._ [2012b] X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo, K. Reuter, and M. Scheffler, ‘‘Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions,’’ New J. Phys. 14, 053020 (2012b).
  * Engel  _et al._ [1991] G. E. Engel, B. Farid, C. M. M. Nex, and N. H. March, ‘‘Calculation of the GW self-energy in semiconducting crystals,’’ Phys. Rev. B 44, 13356--13373 (1991).
  * Backhouse and Booth [2020] O. J. Backhouse and G. H. Booth, ‘‘Efficient Excitations and Spectra within a Perturbative Renormalization Approach,’’ J. Chem. Theory Comput. 16, 6294--6304 (2020).
  * Backhouse, Santana-Bonilla, and Booth [2021] O. J. Backhouse, A. Santana-Bonilla, and G. H. Booth, ‘‘Scalable and Predictive Spectra of Correlated Molecules with Moment Truncated Iterated Perturbation Theory,’’ J. Phys. Chem. Lett. 12, 7650--7658 (2021).
  * Backhouse and Booth [2022] O. J. Backhouse and G. H. Booth, ‘‘Constructing “Full-Frequency”Spectra via Moment Constraints for Coupled Cluster Green’s Functions,’’ J. Chem. Theory Comput. 18, 6622--6636 (2022).
  * Sriluckshmy  _et al._ [2021] P. V. Sriluckshmy, M. Nusspickel, E. Fertitta, and G. H. Booth, ‘‘Fully algebraic and self-consistent effective dynamics in a static quantum embedding,’’ Phys. Rev. B 103, 085131 (2021).
  * Haydock [1980] R. Haydock, ‘‘The recursive solution of the Schrödinger equation,’’ Comp. Phys. Comms. 20, 11--16 (1980).
  * Adachi and Lipparini [1988] S. Adachi and E. Lipparini, ‘‘Sum rules in extended RPA theories,’’ Nuc. Phys. A 489, 445--460 (1988).
  * Karlsson and van Leeuwen [2016] D. Karlsson and R. A. van Leeuwen, ‘‘Partial self-consistency and analyticity in many-body perturbation theory: Particle number conservation and a generalized sum rule,’’ Physical Review B 94, 125124 (2016).
  * Sabin  _et al._ [2010] J. R. Sabin, J. Oddershede, R. Cabrera-Trujillo, S. P. Sauer, E. Deumens, and Y. Öhrn, ‘‘Stopping power of molecules for fast ions,’’ Mol. Phys. 108, 2891--2897 (2010), https://doi.org/10.1080/00268976.2010.508753 .
  * Van Caillie and Amos [2000] C. Van Caillie and R. D. Amos, ‘‘Static and dynamic polarisabilities, Cauchy coefficients and their anisotropies: an evaluation of DFT functionals,’’ Chem. Phys. Lett. 328, 446--452 (2000).
  * Kalugina and Thakkar [2016] Y. N. Kalugina and A. J. Thakkar, ‘‘Ab initio calculations of static dipole polarizabilities and Cauchy moments for the halomethanes,’’ Chem. Phys. Lett. 644, 20--24 (2016).
  * Eshuis, Yarkony, and Furche [2010] H. Eshuis, J. Yarkony, and F. Furche, ‘‘Fast computation of molecular random phase approximation correlation energies using resolution of the identity and imaginary frequency integration,’’ J. Chem. Phys. 132, 234114 (2010).
  * Backhouse, Nusspickel, and Booth [2020] O. J. Backhouse, M. Nusspickel, and G. H. Booth, ‘‘Wave Function Perspective and Efficient Truncation of Renormalized Second-Order Perturbation Theory,’’ J. Chem. Theory Comput. 16, 1090--1104 (2020).
  * Rebolini and Toulouse [2016] E. Rebolini and J. Toulouse, ‘‘Range-separated time-dependent density-functional theory with a frequency-dependent second-order Bethe-Salpeter correlation kernel,’’ J. Chem. Phys. 144 (2016).
  * Véril  _et al._ [2018] M. Véril, P. Romaniello, J. A. Berger, and P.-F. Loos, ‘‘Unphysical Discontinuities in GW Methods,’’ J. Chem. Theory Comput. 14, 5220--5228 (2018).
  * Quintero-Monsebaiz  _et al._ [2022] R. Quintero-Monsebaiz, E. Monino, A. Marie, and P.-F. Loos, ‘‘Connections between many-body perturbation and coupled-cluster theories,’’ J. Chem. Phys. 157, 231102 (2022), https://doi.org/10.1063/5.0130837 .
  * Bintrim and Berkelbach [2021] S. J. Bintrim and T. C. Berkelbach, ‘‘Full-frequency GW without frequency,’’ J. Chem. Phys. 154, 041101 (2021).
  * Bintrim and Berkelbach [2022] S. J. Bintrim and T. C. Berkelbach, ‘‘Full-frequency dynamical Bethe--Salpeter equation without frequency and a study of double excitations,’’ J. Chem. Phys. 156, 044114 (2022).
  * Tölle and Chan [2022] J. Tölle and G. K.-L. Chan, ‘‘Exact relationships between the GW approximation and equation-of-motion coupled-cluster theories through the quasi-boson formalism,’’  (2022).
  * Meyer and Pal [1989] H.-D. Meyer and S. Pal, ‘‘A band-Lanczos method for computing matrix elements of a resolvent,’’ J. Chem. Phys. 91, 6195--6204 (1989).
  * Weikert  _et al._ [1996] H.-G. Weikert, H.-D. Meyer, L. S. Cederbaum, and F. Tarantelli, ‘‘Block Lanczos and many-body theory: Application to the one-particle Green’s function,’’ J. Chem. Phys. 104, 7122--7138 (1996).
  * Haydock and Nex [1985] R. Haydock and C. M. M. Nex, ‘‘A general terminator for the recursion method,’’ J. Phys. C: Solid State Physics 18, 2235 (1985).
  * Engel and Farid [1992] G. E. Engel and B. Farid, ‘‘Calculation of the dielectric properties of semiconductors,’’ Phys. Rev. B 46, 15812--15827 (1992).
  * Farid [2002] B. Farid, ‘‘Dynamical correlation functions expressed in terms of many-particle ground-state wavefunction; the dynamical self-energy operator,’’ Phil. Mag. B 82, 1413--1610 (2002), https://doi.org/10.1080/13642810208222682 .
  * Fukaya  _et al._ [2014] T. Fukaya, Y. Nakatsukasa, Y. Yanagisawa, and Y. Yamamoto, ‘‘CholeskyQR2: A Simple and Communication-Avoiding Algorithm for Computing a Tall-Skinny QR Factorization on a Large-Scale Parallel System,’’ in  _2014 5th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems_ (2014) pp. 31--38.
  * Fukaya  _et al._ [2020] T. Fukaya, R. Kannan, Y. Nakatsukasa, Y. Yamamoto, and Y. Yanagisawa, ‘‘Shifted Cholesky QR for Computing the QR Factorization of Ill-Conditioned Matrices,’’ SIAM J. Sci. Comput. 42, A477--a503 (2020).
  * Schirmer, Trofimov, and Stelter [1998] J. Schirmer, A. B. Trofimov, and G. Stelter, ‘‘A non-Dyson third-order approximation scheme for the electron propagator,’’ J. Chem. Phys. 109, 4734--4744 (1998).
  * Schirmer [2018] J. Schirmer,  _Many-Body Methods for Atoms, Molecules and Clusters_ (Springer Chem, Switzerland, 2018).
  * Berkelbach [2018] T. C. Berkelbach, ‘‘Communication: Random-phase approximation excitation energies from approximate equation-of-motion coupled-cluster doubles,’’ J. Chem. Phys. 149, 041103 (2018), https://doi.org/10.1063/1.5032314 .
  * Rishi, Perera, and Bartlett [2020] V. Rishi, A. Perera, and R. J. Bartlett, ‘‘A route to improving rpa excitation energies through its connection to equation-of-motion coupled cluster theory,’’ The Journal of Chemical Physics 153, 234101 (2020), https://doi.org/10.1063/5.0023862 .
  * Dempwolff  _et al._ [2019] A. L. Dempwolff, M. Schneider, M. Hodecker, and A. Dreuw, ‘‘Efficient implementation of the non-Dyson third-order algebraic diagrammatic construction approximation for the electron propagator for closed- and open-shell molecules,’’ J. Chem. Phys. 150 (2019).
  * Trofimov and Schirmer [2005] A. B. Trofimov and J. Schirmer, ‘‘Molecular ionization energies and ground- and ionic-state properties using a non-Dyson electron propagator approach,’’ J. Chem. Phys. 123, 144115 (2005).
  * Scott and Booth [2021] C. J. C. Scott and G. H. Booth, ‘‘Extending density matrix embedding: A static two-particle theory,’’ Phys. Rev. B 104, 245114 (2021).
  * Langreth [1970] D. C. Langreth, ‘‘Singularities in the x-ray spectra of metals,’’ Phys. Rev. B 1, 471--477 (1970).
  * Nozières and De Dominicis [1969] P. Nozières and C. T. De Dominicis, ‘‘Singularities in the x-ray absorption and emission of metals. iii. one-body theory exact solution,’’ Phys. Rev. 178, 1097--1107 (1969).
  * Lundqvist and Samathiyakanit [1969] B. Lundqvist and V. Samathiyakanit, ‘‘Single-particle spectrum of the degenerate electron gas iv. ground state energy.’’ Phys. kondens. Materie. , 231–235 (1969).
  * Heßelmann and Görling [2011] A. Heßelmann and A. Görling, ‘‘Random-phase approximation correlation methods for molecules and solids,’’ Mol. Phys. 109, 2473--2500 (2011).
  * Chen  _et al._ [2017] G. P. Chen, V. K. Voora, M. M. Agee, S. G. Balasubramani, and F. Furche, ‘‘Random-Phase Approximation Methods,’’ Annu. Rev. Phys. Chem. 68, 421--445 (2017).
  * Parkinson and Zerner [1989] W. A. Parkinson and M. C. Zerner, ‘‘The calculation of dynamic molecular polarizability,’’ J. Chem. Phys. 90, 5606--5611 (1989).
  * Note [1] 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 $(\mathbf{X+Y})^{-1}=(\mathbf{X-Y})^{T}$, that is $(\mathbf{\eta}^{(n)})^{-1}=(\mathbf{X}-\mathbf{Y})\mathbf{\Omega}^{-n}(\mathbf{X}-\mathbf{Y})^{T}$.
  * Furche [2001] F. Furche, ‘‘On the density matrix based approach to time-dependent density functional response theory,’’ J. Chem. Phys. 114, 5982--5992 (2001), https://doi.org/10.1063/1.1353585 .
  * Ángyán  _et al._ [2011] J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, ‘‘Correlation energy expressions from the adiabatic-connection fluctuation--dissipation theorem approach,’’ J. Chem. Theory Comput. 7, 3116--3130 (2011), https://doi.org/10.1021/ct200501r .
  * Furche [2008] F. Furche, ‘‘Developing the random phase approximation into a practical post-Kohn--Sham correlation model,’’ J. Chem. Phys. 129, 114105 (2008).
  * Note [2] 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 $\mathcal{O}[N]$ (in this case the auxiliary basis).
  * Hummel  _et al._ [2019] F. Hummel, A. Grüneis, G. Kresse, and P. Ziesche, ‘‘Screened Exchange Corrections to the Random Phase Approximation from Many-Body Perturbation Theory,’’ J. Chem. Theory Comput. 15, 3223--3236 (2019).
  * Förster [2022] A. Förster, ‘‘Assessment of the Second-Order Statically Screened Exchange Correction to the Random Phase Approximation for Correlation Energies,’’ J. Chem. Theory Comput. 18, 5948--5965 (2022).
  * Hale, Higham, and Trefethen [2008] N. Hale, N. J. Higham, and L. N. Trefethen, ‘‘Computing $A^{\alpha},\log(A)$, and Related Matrix Functions by Contour Integrals,’’ SIAM J. Numer. Anal. 46, 2505--2523 (2008).
  * Wilhelm  _et al._ [2016] J. Wilhelm, P. Seewald, M. Del Ben, and J. Hutter, ‘‘Large-Scale Cubic-Scaling Random Phase Approximation Correlation Energy Calculations Using a Gaussian Basis,’’ J. Chem. Theory Comput. 12, 5851--5859 (2016), pMID: 27779863.
  * Förster and Visscher [2020] A. Förster and L. Visscher, ‘‘Low-Order Scaling G0W0 by Pair Atomic Density Fitting,’’ J. Chem. Theory Comput. 16, 7381--7399 (2020), pMID: 33174743.
  * Parrish  _et al._ [2012] R. M. Parrish, E. G. Hohenstein, T. J. Martínez, and C. D. Sherrill, ‘‘Tensor hypercontraction. II. Least-squares renormalization,’’ J. Chem. Phys. 137, 224106 (2012).
  * Lu and Thicke [2017] J. Lu and K. Thicke, ‘‘Cubic scaling algorithms for RPA correlation using interpolative separable density fitting,’’ J. Comput. Phys. 351, 187--202 (2017).
  * Shenvi  _et al._ [2014] N. Shenvi, H. van Aggelen, Y. Yang, and W. Yang, ‘‘Tensor hypercontracted ppRPA: Reducing the cost of the particle-particle random phase approximation from O($r^{6}$) to O($r^{4}$),’’ J. Chem. Phys. 141, 024119 (2014).
  * Lu and Ying [2015] J. Lu and L. Ying, ‘‘Compression of the electron repulsion integral tensor in tensor hypercontraction format with cubic scaling cost,’’ J. Comput. Phys 302, 329--335 (2015).
  * Pierce and Valeev [2023] K. Pierce and E. F. Valeev, ‘‘Efficient Construction of Canonical Polyadic Approximations of Tensor Networks,’’ J. Chem. Theory Comput. 19, 71--81 (2023), https://doi.org/10.1021/acs.jctc.2c00861 .
  * Schurkus and Ochsenfeld [2016] H. F. Schurkus and C. Ochsenfeld, ‘‘Communication: An effective linear-scaling atomic-orbital reformulation of the random-phase approximation using a contracted double-Laplace transformation,’’ J. Chem. Phys. 144, 31101 (2016).
  * Helmich-Paris and Visscher [2016] B. Helmich-Paris and L. Visscher, ‘‘Improvements on the minimax algorithm for the Laplace transformation of orbital energy denominators,’’ J. Comput. Phys 321, 927--931 (2016).
  * Graf  _et al._ [2018] D. Graf, M. Beuerle, H. F. Schurkus, A. Luenser, G. Savasci, and C. Ochsenfeld, ‘‘Accurate and Efficient Parallel Implementation of an Effective Linear-Scaling Direct Random Phase Approximation Method,’’ J. Chem. Theory Comput. 14, 2505--2515 (2018).
  * Sun  _et al._ [2017] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, ‘‘PySCF: the Python-based simulations of chemistry framework,’’ Wiley Interdiscip. Rev. Comput. Mol. Sci. 8, e1340 (2017).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, ‘‘Recent developments in the PySCF program package,’’ J. Chem. Phys. 153, 024109 (2020).
  * Zhu and Chan [2021a] T. Zhu and G. K.-L. Chan, ‘‘Ab Initio Full Cell $GW+\mathrm{DMFT}$ for Correlated Materials,’’ Phys. Rev. X 11, 021006 (2021a).
  * Zhu and Chan [2021b] T. Zhu and G. K.-L. Chan, ‘‘All-Electron Gaussian-Based G0W0 for Valence and Core Excitation Energies of Periodic Systems,’’ J. Chem. Theory Comput. 17, 727--741 (2021b).
  * Wilhelm, Del Ben, and Hutter [2016] J. Wilhelm, M. Del Ben, and J. Hutter, ‘‘GW in the Gaussian and Plane Waves Scheme with Application to Linear Acenes,’’ J. Chem. Theory Comput. 12, 3623--3635 (2016).
  * Caruso  _et al._ [2016] F. Caruso, M. Dauth, M. J. van Setten, and P. Rinke, ‘‘Benchmark of GW Approaches for the GW100 Test Set,’’ J. Chem. Theory Comput. 12, 5076--5087 (2016).
  * Lange and Berkelbach [2018] M. F. Lange and T. C. Berkelbach, ‘‘On the Relation between Equation-of-Motion Coupled-Cluster Theory and the GW Approximation,’’ J. Chem. Theory Comput. 14, 4224--4236 (2018).
  * Loos, Romaniello, and Berger [2018] P.-F. Loos, P. Romaniello, and J. A. Berger, ‘‘Green functions and self-consistency: Insights from the spherium model,’’ J. Chem. Theory Comput. 14, 3071--3082 (2018), pMID: 29746773, https://doi.org/10.1021/acs.jctc.8b00260 .
  * Di Sabatino, Loos, and Romaniello [2021] S. Di Sabatino, P.-F. Loos, and P. Romaniello, ‘‘Scrutinizing gw-based methods using the hubbard dimer,’’ Front. Chem. 9 (2021), 10.3389/fchem.2021.751054.
  * Monino and Loos [2022] E. Monino and P.-F. Loos, ‘‘Unphysical discontinuities, intruder states and regularization in $GW$ methods,’’ J. Chem. Phys. 156, 231101 (2022).
  * Weiße  _et al._ [2006] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, ‘‘The kernel polynomial method,’’ Rev. Mod. Phys. 78, 275--306 (2006).
  * Bhatia and Rosenthal [1997] R. Bhatia and P. Rosenthal, ‘‘How and Why to Solve the Operator Equation AX-XB = Y,’’ Bulletin of the London Mathematical Society 29, 1--21 (1997).
  * Häser and Almlöf [1992] M. Häser and J. Almlöf, ‘‘Laplace transform techniques in Møller–Plesset perturbation theory,’’ J. Chem. Phys. 96, 489--494 (1992).

## Appendix A: Improved Numerical Quadrature

While the integral expression for $\tilde{\eta}^{(0)}$ derived in Eq. 70 in
the main text is sufficient for numerical quadrature in $\mathcal{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
$\tilde{\eta}^{(0)}$, i.e. setting the Coulomb interaction to zero in the
expression for
$[(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})]^{\frac{1}{2}}$. This just
leaves a contribution from the irreducible polarizability, which can be
analytically integrated, as

$\displaystyle\mathbf{D}\mathbf{T}$
$\displaystyle=(\mathbf{D}^{2})^{\frac{1}{2}}\mathbf{T}$
$\displaystyle=\frac{1}{\pi}\int_{-\infty}^{\infty}\left(\mathbf{I}-z^{2}(\mathbf{D}^{2}+z^{2}\mathbf{I})^{-1}\right)\mathbf{T}dz,$
(A1)

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

$\tilde{\eta}^{(0)}=\mathbf{D}\mathbf{T}+\frac{1}{\pi}\int_{-\infty}^{\infty}z^{2}\mathbf{F}(z)\mathbf{S}_{L}(\mathbf{I}+\mathbf{Q}(z))^{-1}\mathbf{S}_{R}^{T}\mathbf{F}(z)\mathbf{T}dz,$
(A2)

which substantially reduces the magnitude of the numerically integrated
component of $\tilde{\eta}^{(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 $\mathbf{F}(z)$ and $\mathbf{Q}(z)$
decay as $z^{-2}$ in the large-$z$ limit. Series expanding
$[\mathbf{I}+\mathbf{Q}(z)]^{-1}\sim\mathbf{I}-\mathbf{Q}(z)+\mathbf{Q}^{2}(z)-\dots$,
we find that the leading order decay of the integrand in Eq. A2 is
$\mathcal{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 $\mathcal{O}[z^{-4}]$
rate.

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

$\mathbf{C}_{2}=\frac{1}{\pi}\int_{-\infty}^{\infty}z^{2}\mathbf{F}(z)\mathbf{S}_{L}\mathbf{S}_{R}^{T}\mathbf{F}(z)\mathbf{T}dz.$
(A3)

This can be analytically integrated via the residue theorem, to give

$\mathbf{C}_{2}=\left((\mathbf{S}_{L}\mathbf{S}_{R}^{T})\circ\mathbf{E}\right)\mathbf{T},$
(A4)

where $\circ$ denotes the Hadamard or element-wise product, and we define
$\mathbf{E}$ as

$\mathbf{E}_{ia,jb}=(\mathbf{D}_{ia,ia}+\mathbf{D}_{jb,jb})^{-1}.$ (A5)

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

$(\mathbf{C}_{2})_{ia,P}=\sum_{jb}\frac{\sum_{Q}(\mathbf{S}_{L})_{ia,Q}(\mathbf{S}_{R})_{jb,Q}}{\epsilon_{a}+\epsilon_{b}-\epsilon_{i}-\epsilon_{j}}\mathbf{T}_{jb,P}.$
(A6)

However, in this form it is unfortunately unable to be efficiently evaluated
in $\mathcal{O}[N^{4}]$ time (instead scaling as $\mathcal{O}[N^{5}]$). We
note that a similar subtraction of this leading-order term in the computation
of the RPA correlation energy is in contrast efficiently computable, due to
the presence of a cyclically-invariant trace operation [56].

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

$((\mathbf{S}_{L}\mathbf{S}_{R}^{T})\circ\mathbf{E})\mathbf{D}+\mathbf{D}((\mathbf{S}_{L}\mathbf{S}_{R}^{T})\circ\mathbf{E})=\mathbf{S}_{L}\mathbf{S}_{R}^{T},$
(A7)

which can be verified by substitution of Eq. A5. Writing the matrix of
interest we wish to efficiently compute as
$\mathbf{N}=(\mathbf{S}_{L}\mathbf{S}_{R}^{T})\circ\mathbf{E}$, this Sylvester
equation has an integral solution of the form [113]

$\mathbf{N}=\int_{0}^{\infty}e^{-t\mathbf{D}}\mathbf{S}_{L}\mathbf{S}_{R}^{T}e^{-t\mathbf{D}}dt.$
(A8)

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

$\displaystyle\mathbf{D}\mathbf{N}+\mathbf{N}\mathbf{D}$
$\displaystyle=\int_{0}^{\infty}\left(\mathbf{D}e^{-t\mathbf{D}}\mathbf{S}_{L}\mathbf{S}_{R}^{T}e^{-t\mathbf{D}}\right.$
$\displaystyle\qquad+\left.e^{-t\mathbf{D}}\mathbf{S}_{L}\mathbf{S}_{R}^{T}e^{-t\mathbf{D}}\mathbf{D}\right)dt$
(A9)
$\displaystyle=-\int_{0}^{\infty}\frac{d}{dt}\left(e^{-t\mathbf{D}}\mathbf{S}_{L}\mathbf{S}_{R}^{T}e^{-t\mathbf{D}}\right)dt$
(A10)
$\displaystyle=-\left[e^{-t\mathbf{D}}\mathbf{S}_{L}\mathbf{S}_{R}^{T}e^{-t\mathbf{D}}\right]_{0}^{\infty}$
(A11) $\displaystyle=\mathbf{S}_{L}\mathbf{S}_{R}^{T}.$ (A12)

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

$\displaystyle\tilde{\eta}^{(0)}$ $\displaystyle=\mathbf{D}\mathbf{T}$
$\displaystyle+\int_{0}^{\infty}e^{-t\mathbf{D}}\mathbf{S}_{L}\mathbf{S}_{R}^{T}e^{-t\mathbf{D}}\mathbf{T}dt$
$\displaystyle+\frac{1}{\pi}\int_{-\infty}^{\infty}z^{2}\mathbf{F}(z)\mathbf{S}_{L}\left((\mathbf{I}+\mathbf{Q}(z))^{-1}-\mathbf{I}\right)\mathbf{S}_{R}^{T}\mathbf{F}(z)\mathbf{T}dz.$
(A13)

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

## Appendix B: Numerical Integration Error Estimates

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

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

$(\tilde{\eta}^{(0)}_{\infty})^{T}(\mathbf{A}+\mathbf{B})\tilde{\eta}^{(0)}_{\infty}=\mathbf{V}^{T}(\mathbf{A}-\mathbf{B})\mathbf{V}.$
(B1)

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

$\displaystyle X(n_{p})=$
$\displaystyle(\tilde{\eta}^{(0)}_{n_{p}})^{T}(\mathbf{A}+\mathbf{B})\tilde{\eta}^{(0)}_{n_{p}}-\mathbf{V}^{T}(\mathbf{A}-\mathbf{B})\mathbf{V}$
(B2) $\displaystyle=$
$\displaystyle(\Delta\tilde{\eta}^{(0)}_{n_{p}})^{T}(\mathbf{A}+\mathbf{B})\tilde{\eta}^{(0)}_{n_{p}}+(\tilde{\eta}^{(0)}_{n_{p}})^{T}(\mathbf{A}+\mathbf{B})\Delta\tilde{\eta}^{(0)}_{n_{p}}$
$\displaystyle-(\Delta\tilde{\eta}^{(0)}_{n_{p}})^{T}(\mathbf{A}+\mathbf{B})\Delta\tilde{\eta}^{(0)}_{n_{p}}.$
(B3)

We can then explicitly evaluate $||X(n_{p})||_{2}$ in $\mathcal{O}[N^{4}]$
time under the same conditions as the rest of this manuscript, and apply the
properties of the Frobenius norm to write

$\displaystyle||X(n_{p})||_{2}\leq$ $\displaystyle
2||\Delta\tilde{\eta}^{(0)}_{n_{p}}||_{2}||(\mathbf{A}+\mathbf{B})\tilde{\eta}^{(0)}_{n_{p}}||_{2}$
(B4)
$\displaystyle-||\Delta\tilde{\eta}^{(0)}_{n_{p}}||_{2}^{2}||(\mathbf{A}+\mathbf{B})||_{2}.$

Writing $||\Delta\tilde{\eta}^{(0)}_{n_{p}}||_{2}=x$, the error is therefore
rigorously bound by the solutions satisfying the quadratic inequality,

$x^{2}-2\frac{||(\mathbf{A}+\mathbf{B})\tilde{\eta}^{(0)}_{n_{p}}||_{2}}{||(\mathbf{A}+\mathbf{B})||_{2}}x+\frac{||X(n_{p})||_{2}}{||(\mathbf{A}+\mathbf{B})||_{2}}\leq
0.$ (B5)

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 $\mathcal{O}[N^{4}]$ scaling of evaluating $X(n_{p})$ this makes this
estimate of limited practical use, and we seek an improved estimate, that
provides a tighter upper bound on the true error.

To do this we take advantage of the nested nature of the Clenshaw--Curtis
quadrature within the integral providing the dominant error contribution.
Using this, we can obtain estimates at both one half and one quarter of the
current number of integration points with no additional cost. We show how this
can give an estimate of the error in
$\tilde{\eta}^{(0)}_{n_{p}}\left(\mathbf{A}+\mathbf{B}\right)$, which we can
relate to the error in $\tilde{\eta}^{(0)}_{n_{p}}$. Parameterising the
exponential convergence of the $L_{2}$ norm error as
$||\Delta\tilde{\eta}^{(0)}_{n_{p}}||_{2}=\alpha e^{-\beta n_{p}}$, we can
rigorously relate the differences between estimates at different levels of
sampling, as

$\displaystyle\Delta_{n_{1},n_{2}}$
$\displaystyle=||\tilde{\eta}^{(0)}_{n_{1}}-\tilde{\eta}^{(0)}_{n_{2}}||_{2}$
(B6)
$\displaystyle\leq||\Delta\tilde{\eta}^{(0)}_{n_{1}}||_{2}+||\Delta\tilde{\eta}^{(0)}_{n_{2}}||_{2}$
(B7) $\displaystyle\approx\alpha\left(e^{-\beta n_{1}}+e^{-\beta
n_{2}}\right).$ (B8)

We seek to estimate the value of $||\Delta\tilde{\eta}^{(0)}_{4n_{p}}||_{2}$
given the values $\Delta_{4n_{p},2n_{p}}$ and $\Delta_{4n_{p},n_{p}}$ which
are computed as intermediates in the nested quadrature of $4n_{p}$ points.

Setting $x=e^{-\beta n_{p}}$, and using the approximation defined in Eq. B8
for $\Delta_{4n_{p},2n_{p}}=\alpha(x^{4}+x^{2})$ and
$\Delta_{4n_{p},n_{p}}=\alpha(x^{4}+x)$, eliminating the factor of $\alpha$ we
find

$(\Delta_{4n_{p},n_{p}}-\Delta_{4n_{p},2n_{p}})x^{3}+\Delta_{4n_{p},n_{p}}x-\Delta_{4n_{p},2n_{p}}=0.$
(B9)

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

$||\Delta\tilde{\eta}^{(0)}_{4n_{p}}||_{2}\approx\Delta_{4n_{p},2n_{p}}(1+x^{-2})^{-1}.$
(B10)

This estimator provides a consistent and systematic overestimate of the error,
allowing us to use this as a reliable upper bound to the true error of the
quadrature. This property, combined with only requiring the quantities
$\Delta_{4n_{p},2n_{p}}$ and $\Delta_{4n_{p},n_{p}}$, which are unavoidable
intermediates of the nested integration, makes this essentially a
computationally free estimate of the error for a given number of grid points,
and therefore as a reliable check for convergence. If the resulting error
estimate is too large, the number of grid points can be increased by a factor
of two, and all points already evaluated can be reused in the next estimate.
The convergence of both error estimates can be seen in Fig. 1 of the main
text, with the result of Eq. B10 denoted ‘Nested Fit’, while the lower bound
from the solution of Eq. B5 denoted ‘Lower Bound’.