A unified Green's function approach for spectral and thermodynamic properties from algorithmic inversion of dynamical potentials

Dynamical potentials appear in many advanced electronic-structure methods, including self-energies from many-body perturbation theory, dynamical mean-field theory, electronic-transport formulations, and many embedding approaches. Here, we propose a novel treatment for the frequency dependence, introducing an algorithmic inversion method that can be applied to dynamical potentials expanded as sum over poles. This approach allows for an exact solution of Dyson-like equations at all frequencies via a mapping to a matrix diagonalization, and provides simultaneously frequency-dependent (spectral) and frequency-integrated (thermodynamic) properties of the Dyson-inverted propagators. The transformation to a sum over poles is performed introducing $n$-th order generalized Lorentzians as an improved basis set to represent the spectral function of a propagator, and using analytic expressions to recover the sum-over-poles form. Numerical results for the homogeneous electron gas at the $G_0W_0$ level are provided to argue for the accuracy and efficiency of such unified approach.


I. INTRODUCTION
Electronic-structure calculations have been and remain a powerful and ever expanding field of research to understand and predict materials properties [1]. The development of methods, algorithms, and hardware brings in continuous progress, allowing for computational materials discovery [2][3][4], accurate comparison with experiments [5,6], and even hybrid quantum-computation algorithms [7,8].
Due to the interaction between the electrons in a system, solving the many-body quantum problem is often at the core of many approaches. Focusing on condensedmatter systems, density-functional theory (DFT) has been one of the most used and successful methods so far [9]. The possibility to map exactly the groundstate solution of the N -body problem to the minimization of a density functional for the energy [10] offers great computational simplifications, allowing the accurate computation of ground-state quantities for most materials. Although mathematically well-defined [11] and computationally inexpensive, it remains challenging to improve the approximate functionals [12,13] -often resulting in incorrect predictions for complex or stronglycorrelated systems [14] -or to address spectroscopic properties [15,16].
Dynamical (i.e. frequency-dependent) theories like many-body perturbation theory, dynamical mean-field theory, and in general embedding theories offer the flexibility to overcome these limitations of DFT. While the type of embedding differs in different approaches, a common element is the appearance of dynamical potentials. * corresponding author: tommaso.chiarotti@epfl.ch As an example, many-body perturbation theory (MBPT) reduces the multi-particle electronic degrees of freedom to one via frequency embedding [17]. Dynamical meanfield theory (DMFT) couples a real-space impurity with the rest of the system, requiring self-consistency between the two self-energies acting on the impurity and on the bath [18]. Self-energy-embedding theory (SEET) calculates exactly the frequency-dependent self-energy of strongly-correlated manifolds in solids, and applies it to the remaining weakly interacting orbitals [19,20]. Coherent electronic-transport theories use a Green's function embedding to calculate the electronic conductance of e.g. a conductor between two semi-infinite leads, coupling the three systems dynamically [21,22]. Clearly, handling properly frequency-dependent potentials is of central interest in the field.
Using here MBPT as a paradigmatic example, we highlight that the difficulty in treating dynamical quantities has often led to different methodological approaches when calculating spectral or thermodynamic quantities (such as energies, number of particles, chemical potentials). Real-axis calculations are commonly performed to compute the frequency-dependent spectral properties [23][24][25], while the frequency-integrated thermodynamic properties are typically calculated using an imaginary-axis formalism [26][27][28][29][30][31]. In a series of papers [32][33][34][35] von Barth and coworkers have proposed a formalism partially able to tackle spectra and thermodynamics together for the homogeneous electron gas [36], by modelling the spectral function in frequencymomentum space using Gaussians with k-parametrized centers (quasi-particle energies), broadening (weights), and satellites. Due to its model nature, the approach does not easily offer the flexibility to target realistic systems and in general extend to embedding problems.
Here we introduce a novel approach, termed arXiv:2109.07972v1 [cond-mat.mtrl-sci] 16 Sep 2021 algorithmic-inversion method, applied on sum-over-pole expansions (AIM-SOP), to address the simultaneous calculation of accurate spectral and thermodynamic quantities. Within AIM-SOP, dynamical (frequencydependent) self-energies are expanded on sum over poles, and the exact solution -at all frequencies -of the Dyson equation is found via a matrix diagonalization. The transformation of a frequency-dependent propagator into a SOP via a representation of its spectral function on a target basis set is greatly improved with the introduction of n-th order generalized Lorentzians as a basis with improved decay properties. The SOP form allows one to compute analytically convolutions and moments of propagators for the calculation of spectral, and thermodynamic properties. Owing to the fulfillment of all sum rules implied by the Dyson equation, we show that the AIM-SOP method becomes essential to have accurate frequency-integrated quantities in a real-axis (thus, spectral oriented) formalism. As a case study, we consider the paradigmatic case of the homogeneous electron gas (HEG), for r s from 1 to 10, treated at the G 0 W 0 level [37][38][39].
The paper is organized as follows: In Sec. II we introduce the AIM-SOP approach, discussing its main goal and the SOP form for propagators and self-energies. In Sec. II A we provide an overview of the connection between a propagator and its spectral function, first for a continuum and then extending it to treat spectral functions represented on discrete basis sets, as will be used in this work. Then, we consider different basis sets to represent the spectral function and obtain a SOP representation introducing n-th order Lorentzians. In Sec. II B we provide the numerical procedure to transform a propagator sampled on a grid to a SOP representation, and viceversa. In Sec. II C we introduce several useful expressions when dealing with propagators on SOP, such as analytic convolutions and moments, and in Sec. II D we show with a numerical example the representation on SOP for a test propagator. Finally, in Sec. II E we present the algorithmic-inversion method on sum over poles to obtain exact solutions on SOP of any Dyson-like equation. We first provide a mathematical proof for the case of a self-energy on SOP, then we discuss the case of the polarizability inversion, providing a numerical example as proof-of-concept for the procedure. In Sec. III we discuss the application of the method to the test case of the homogeneous electron gas. In Sec. IV we discuss the results obtained applying AIM-SOP to the homogeneous electron gas at the G 0 W 0 level, first discussing the r s = 4 case in detail and then presenting results for r s from 1 to 10. Finally, in Sec. V we draw the conclusions for the paper. Technical aspects of the method are further presented in the Appendices.

II. METHOD: AIM-SOP FOR DYNAMICAL POTENTIALS
In this Section we introduce the algorithmic-inversion method to treat dynamical (frequency-dependent) potentials. The crucial goal for AIM-SOP is to solve exactly and at all frequencies Dyson-like equations for dynamical potential expressed as sum over poles. For this purpose we express frequency-dependent propagators and self-energies (or, say, polarizabilities or screened Coulomb interactions) in a SOP form: where the constant term A 0 may be present for selfenergies and potentials. Generally, we consider here having complex residues A i and poles In order to provide the correct analytical structure respecting time-ordering, δ i ≷ 0 when i ≶ µ, where µ is the effective chemical potential of the propagator (for a Green's function µ is the Fermi energy of the system, for a polarizability or a screened potential µ = 0). Throughout this work we will use as case of study the homogeneous electron gas (HEG) also in view of to the extensive algorithmic and numerical results in the literature. In the HEG, due to translational symmetry, the two-point operators (including Green's functions, self-energies, polarizabilities) are diagonal on the plane-wave basis, but the AIM-SOP can be generalized to non-homogeneous systems, as will be discussed in future work.
A. Spectral representations Following Ref. [33], we consider the spectral representation of a propagator (here the Green's function for simplicity), where G is expressed in terms of its spectral function A, by performing a time-ordered Hilbert transform (TOHT), where C is a time-ordered contour which is shifted above/below the real axis for ω ≶ µ, and where the shift is sent to zero after the integral is computed. Accordingly, the inverse relation to go from G to A is given by: the last expression being valid for a scalar Green's function, as is the case for the HEG. Representing the spectral function on a (finite) basis set {b j (ω)}, (4) with b j (ω) centred on j and positive (negative) for j ≶ µ, respectively, and a j > 0, we also induce a representation of G. This is achieved by introducing a discrete time-ordered Hilbert transform (D-TOHT) as where the sign chosen for b j in Eq. (4) gives by construction the time-ordered analytic structure of the Green's function. In the case of all δ j → 0 with the number of b j becoming infinite (continuum representation limit), Eq. (5) becomes the standard TOHT of Eq. (2) (with C shifted by ±i0 + ).
A natural choice is to use a basis of Lorentzian functions centered at different frequencies j , according to: for which the D-TOHT for the single element is analytical, yielding a pole function 1/(ω − z j ) with z j = j + iδ j , with the sign convention for δ j defined as discussed above according to time ordering. Thus, choosing b j as in Eq. (6) induces a SOP representation for G according to Eq. (1), with A i = a i ∈ R. Once the SOP representation of G is known, i.e. poles and amplitudes are known, the grid evaluation (inverse of the above) is trivial and amounts to performing the finite sum in Eq. (1). This approach ensures a full-frequency treatment of the propagator (approaching the continuum representation limit where Lorentzians becomes delta functions), while preserving an explicit knowledge of the analytical structure and continuation of G.
The main drawback of using Lorentzians to represent G is related to the slowly decaying tails (1/ω 2 for ω → ∞) induced in the spectral function when using finite broadening values δ j . In order to improve on this, we introduce here n-th order generalized Lorentzians to obtain fast-decay basis functions. These are defined as where N n = n sin π 2n −1 is the normalization factor (see Appendix A). The D-TOHT of L n δ remains analytic and still yields a SOP representation for G (see Appendix A): with residues α m and poles ζ j,m given by Importantly, α m are complex (and so become the residues A i = a j α m in the SOP representation, i being a combined index). Thus, the spectral function of this SOP has contribution by both the real and the imaginary part of each Lorentzian pole 1/(ω − ζ j,m ), resulting in a overall faster decay than each single Lorentzian. Also, it is worth noting that, as for standard Lorentzians, a normalized n-th-order-Lorentzian approaches a Dirac delta for δ j → 0 + . Owing to their fast decay and to this last property, using a SOP for G 0 in term of n-th Lorentzians provides a faster convergence for δ → 0 + , in comparison with a SOP representation built on ordinary Lorentzians, as will be also shown later. While the use of n-th order generalized Lorentzians to represent the spectral function A(ω) provides a faster decay in the imaginary-part of the propagator, it results in a multiplication of the number of poles in the SOP for G (by the degree of the Lorentzian), and in having complex residues. As it will be shown in Sec. II C, the decay properties are fundamental for evaluating the moments of a SOP representation, assuring absolute convergence up to order 2(n − 1). Also, the use of faster decay basis elements when representing the spectral function improves on the stability of the representation procedure, reducing the off-diagonal elements of the overlap matrix of the basis (see Sec. II B).
Alternatively to n-th order Lorentzians, one could consider e.g. using Gaussian functions to represent A(ω), and consequently G(ω), as done in Refs. [32,33]. Gaussians also allow for an analytical expression of the D-TOHT, at the price, though, of invoking the Dawson [40] or Faddeeva [41] functions to evaluate the real part of the propagator. Because of this, SOP expressions are not available, and basic operations involving propagators (such as those described in Sec. II C) cannot be evaluated analytically and need to be worked out in other ways, e.g. numerically or recasting the expressions in terms of propagator spectral functions [32].

B. Transform to a sum over poles
Once the SOP representation has been introduced, the next important step is to determine numerically the SOP coefficients A i in Eq. (1), given an evaluation of G on a frequency grid. According to the discussion of Sec. II A, the SOP representation can be seen equivalently as a representation for the Green's function G or for the spectral function A.
As a first case, we consider representing A(ω) according to Eq. (4), and we do it using the basis of nth-order generalized Lorentzians introduced in Eq. (7). First, we obtain the coefficients a j of the representa-tion by performing a non-negative-least-square (NNLS) fit [41,42], thus assuring the positivity of all a j . Then, we use Eqs. (9)(10) to get the SOP representation for the propagator. While the position and broadening ( j , δ j ) of the n-th order Lorentzians could also be optimized by means of a non-linear NNLS fit, here we consider them centred at j = 1 2 (ω j + ω j−1 ) and broadened with δ j = |ω j − ω j−1 |, and we just linearly optimize a j . Also, for numerical reasons we prefer to work with the bare imaginary part of G, i.e. without imposing the sign factor of Eq. (3), since this function is smoother than the actual spectral function A(ω) close to the Fermi level.
Alternatively, one could consider the basis representation induced on G via Eq. (1) in order to directly obtain the A i and z i coefficients (residues and poles). As for A(ω), this can be achieved by a linear or non-linear LS fit (or interpolation) taking advantage of the knowledge of the whole G(ω) on a frequency grid (and not just of A). Interestingly, the SOP representation in Eq. (1) is a special case of a Padè approximant, written as the ratio of polynomials of order N − 1 and N , respectively. Because of this, one can exploit Padè-specific approaches to determine (A i ,z i ), such as, for instance, Thiele's recursive scheme [43]. We found that this leads to a very efficient method when few tens of poles are considered, becoming numerically unstable beyond. Moreover, since the residues are not constrained to be real and positive (A i are actually complex), there is no control over the timeordered position of the poles, and the procedure is nontrivial to extend to the case of n-th order Lorentzians. For the above reasons, in the present work we adopt the first approach, based on the representation of A(ω).

C. Analytical expressions
Once the SOP representation of a dynamical propagator is available, a number of analytical expressions hold. For instance, the convolution of propagators, such as those involved in the evaluation of the independentparticle polarizabilities in terms of the Green's functions, can be evaluated using Cauchy's residue theorem: Using the the SOP for G, the following integrals can also be computed explicitly: where we refer to the term E m [G] as the m-th (regularized) moment of G. We restrict the discussion to the first m = 0 and m = 1 moments, since those are of interest for calculating the number of particles and the total energy in MBPT (see Sec. III B for details). Higher order moments would require a stronger regularization factor in Eq. (12) than e iω0 + . We underline that if one uses an n-th-order Lorentzian basis to represent G(ω) on SOP, the first 2(n − 1) moments coincide with the moments of its occupied spectral function µ −∞ dω ω 2(n−1) A(ω). This is shown in Appendix B.

D. Numerical validation
In the following we highlight numerically some properties of the SOP representation. To this aim, we consider a propagator G obtained as the Hilbert transform (HT) of a Gaussian, analytically expressed via the Faddeeva [41] function (black curves in the top panels of Fig. 1). Here we have assumed the Fermi level to be far enough from the imaginary part of G such that the retarded HT can be used. The objective of the validation is to transform the Faddeeva Green's function sampled on a finite grid to a SOP representation. Following Sec. II B, we represent the Gaussian spectral function on first (ordinary) and second-order Lorentzians, use Eqs. (9) and (10) to get the resulting SOP representation, and then compare the results -orange and green lines for 1 st and 2 nd order basis, respectively-with the starting Faddeeva Green's function calculated on a much finer grid. We choose to feed the NNLS fitting algorithm with 10 sampling points (for the imaginary part of G) and to use 9 basis functions centered at the midpoints of the grid, and broadened with the width of the interval (in order to ensure an exhaustive cover of the domain). Then we use Eqs. (8)(9)(10) to obtain the SOP representation of G (in the case of nth order Lorentzians) from the output of the NNLS procedure.
Underscoring the quality of the representation, the upper panels of Fig. 1 show that using faster-decay 2 nd order Lorentzians provides a more accurate result for both the real (left) and imaginary (right) part of G. In the lower panel we also compare the moments of the Faddeeva function with those obtained from Eq. (12) and Eq. (B1). Since absolute convergence for the first and second moments is not ensured for 1 st order Lorentzians, meaning that Eq. (B1) does not hold, E n>1 [G] has to be calculated according to Eq. (12) and is in general complex. In order to obtain a meaningful result we take its real part, and consider the imaginary one as an error that must be controlled by extending the basis set towards completeness.

E. Algorithmic inversion on SOP
As anticipated in the introduction, within the SOP approach the exact solution at all frequencies of the Dyson equation can be remapped into the diagonalization of a static effective Hamiltonian (Hermitian only under special conditions), a procedure that we refer to as "algorithmic-inversion method on sum over poles" (AIM-SOP); this is a central result for the present work. As mentioned we will use the HEG as a paradigmatic test case, leaving the treatment of the non-homogeneous case to later work. Suppressing then the k momentum index for simplicity, let us suppose to have the SOP representation of the self-energy Σ(ω) and the non-interacting Green's function G 0 (ω) given by Taking advantage of these expressions, the Dyson equation can be rewritten as in which the N + 1 roots of the polynomial are the N + 1 poles of the Green's function (as expected when the self-energy has N poles). Then, the key statement of this Section is that the roots of T N can be obtained as the eigenvalues of the (N + 1)×(N + 1) matrix We prove this statement by observing that the characteristic polynomial of H AIM is T N (ω), and we proceed by induction. Since the N = 1 case is trivial we move to the N -th case: using the Laplace expansion on the last line, the characteristic polynomial of the N -th case can be written as where we have used the induction hypotheses in the first term of the rhs. Applying the same procedure to the last column of the second term we obtain which completes the proof.
Calling i the poles of G we calculate the residues by equating and performing the limit lim ω→zi (ω − z i ) on both sides (Heaviside cover-up method [44]), obtaining: We have thus proven that by knowing Σ represented on SOP, the SOP expression of G can be found by the diagonalization of the AIM-SOP matrix H AIM followed by the evaluation of the residues using Eq. (20).
It is worth noting that the H AIM matrix becomes (or can be made) Hermitian under special conditions. This happens when the self-energy residues Γ i are real and positive, and the self-energy poles have all the same imaginary part iδ, with the usual time-ordered convention according to Eq. (1), also equal to the broadening assumed for the G 0 pole. Then it is possible to include the imaginary part of the poles in the frequency variable ω, and invert G(ω ∈ γ) on this time-ordered-complex path, in order to have H AIM with only the real part of the poles along the diagonal. Finally, in order to have G(ω ∈ R) we analytically continue the solution to the real axis, obtaining Im{ i } = δ i for the SOP of the Green's function.
We also stress that, given a self-energy represented on SOP, the solution provided by the algorithmic-inversion procedure is exact at all frequencies. This ensures the Green's function fulfills all the sum-rules implied by the Dyson equation, including e.g. the normalization of the spectral weight, and the first and second moments sum rules of the spectral function derived in Ref. [32]. This result is crucial when evaluating frequency-integrated quantities of a Green's function, such as the number of particles or the total energy (see Sec. III B).
Besides the solution of the Dyson equation for G, the AIM-SOP can also be used to solve the Dyson equation for the screened Coulomb interaction W (ω), i.e. to compute the SOP representation of W (ω) once a SOP for the irreducible polarizability P (ω) is provided.
Here v c is the Coulomb potential (recalling that we are suppressing the momentum dependence for simplicity). By letting P (ω) = i Si ω−gi , we can write: and following Eq. (21) (multiplied by ω ω ) we have for which the AIM-SOP matrix can be used to find the poles of −1 (ω) and W (ω). The amplitudes of W are easily found using Eq. (20). Note that by multiplying For validation, we apply the AIM-SOP approach to the paradigmatic case of the homogeneous electron gas (HEG), treated at the G 0 W 0 level of theory [6,17,45]. Since we calculate propagators on the real axis we can easily access spectral (frequency-dependent) properties. The calculation of frequency-integrated groundstate quantities (occupation numbers, total energies, and thermodynamic quantities in general) can be obtained directly from the SOP representation of the spectral quantities computed in the procedure. We stress that usually [26,46,47] thermodynamic properties are obtained via additional calculations of propagators (e.g. on the imaginary axis), while in this work spectral properties and integrated quantities are obtained simultaneously using the SOP representation of propagators computed on the real axis.
While some quantities computed using the freepropagator G 0 have known analytical expressions, as is the case for the irreducible polarizability P 0 expressed via the Lindhard function [48,49], here we recompute explicitly all the propagators needed to evaluate the GW self-energy, making the treatment suitable also for selfconsistent calculations. Therefore in the following the only assumption we make is to consider the Green's function as represented on SOP.

A. HEG propagators on the real frequency axis
In order to solve a one-shot G 0 W 0 cycle for the spin-unpolarized HEG, we first need to compute the irreducible polarizability at the independent-particle (or RPA) level, according to the integral where k = |k| and q = |q| are the moduli of the electron and transferred quasi-momenta, respectively. To compute Eq. (24), the frequency integral (convolution) is performed analytically according to Eq. (11). Then we integrate numerically in spherical coordinates by performing the variable change x = |k + q| on the azimuthal angle of k, which allows for the pre-calculation of the analytical convolutions on the two-dimensional (x, k) grid, instead of on the three-dimensional (k, q, θ) space. Exploiting the par-ity of P (ω), it is also possible to limit the k integration to the occupied states (see Appendix C). The numerical integration on the momentum is performed using the trapezoidal rule, which ensures exponential convergence for decaying functions [50]. In order to have a SOP representation for the screened potential W we transform the polarizability calculated on a frequency grid (at fixed momentum q) to a SOP performing a NNLS fitting, following the procedure of Sec. II B. We then solve the Dyson equation using the algorithmic inversion for the polarizability (see Sec. II E) to obtain a SOP for W , and use it for the GW integral. An alternative possibility would be to solve the Dyson equation on a grid (which, due to homogeneity, is an algebraic inversion), and then transform W to a SOP representation. Even admitting for an exact interpolation for the SOP of W on the calculated frequencies (where the Dyson equation is solved on grid), this SOP would suffer from not having solved the Dyson equation for all other frequencies. Very differently, the SOP obtained from the algorithmic inversion provides for an exact solution of the Dyson equation at all frequencies (see Sec. II E). Thus, the sum rules implied by the Dyson equation (moments of the spectral function) are all obeyed by the SOP obtained from the algorithmic inversion, being the exact solution at all frequencies. Conversely, this is not true for the grid inversion where the solution is exact only for isolated frequencies.
Concerning the self-energy integral where W corr = W − v c , we can still use Eq. (11) since we have the SOP representation of W . Again, in Eq. (26) we perform the x = |k + q| change of variable obtaining which allows for fewer convolutions (as for the polarizability integral), and use trapezoidal weights as in Eq. (25) for the momentum integration. The solution of the Dyson equation for the Green's function using the algorithmic inversion, and the calculation of frequencyintegrated (thermodynamic) quantities, are discussed in the next section. In Fig. 3 we show the overall flow chart describing the process of going from the knowledge of the initial Green's function to the calculation of the corresponding self-energy (for the HEG in the GW approximation), as implemented in the heg sgm.x program of the AGWX suite [51], by means of the SOP approach. As opposed to the path in red, where the Dyson equations are solved on grids, in the green path we highlight the protocol followed in the present work. The crucial difference between the two approaches is the use of the algorithmic-inversion method in order to solve exactly the Dyson equation, providing a SOP for W obeying all sum rules implied by the Dyson equation, as previously discussed in this Section.

B. Frequency-integrated quantities and thermodynamics
Having obtained the self-energy on a frequency grid following the procedure described in Sec. III A, we evaluate the Green's function together with some related frequency-integrated quantities. As mentioned, the SOP approach plays here a central role, enabling the possibility of performing analytical integrals for the moments of G, as those involved in the Galitskii-Migdal expression for the total energy [see Eq. (28) below], and thus to have accurate thermodynamic (frequency-integrated) quantities. Moreover, the use of the algorithmic inversion allows for the exact solution the Dyson equation for the Green's function at all frequencies. The conservation of all sum rules (implied by the Dyson equation, see Sec. II) guaranteed by the AIM-SOP is fundamental when calculating the occupied moments of the spectral function. As an example, the normalization condition of the spectral function is automatically satisfied when G on SOP is obtained using the algorithmic inversion, and allows for not having fitting constrains which would be required, e.g., if we were to use a grid inversion.
In order to exploit the AIM-SOP to get G on a SOP, we obtain the SOP representation of the self-energy by performing a NNLS fitting of ImΣ(ω) (see Sec.II B). Then, in order to compute the total energy from the knowledge of the Green's function G, we use the Galitski-Migdal expression [17,48], here in Hartree units, where V is the volume of the periodic cell of the electron gas. In this expression, the frequency integrals are performed using the SOP for G, and exploiting Eq. (12) with m = 1 and m = 0 for the first and second terms, respectively. Here n k is the k-resolved occupation function, which sums to the total number of particles when integrated over momentum, and k is the occupied band, i.e. the first momentum of the occupied spectral function. For both m = 0 and m = 1 moments, the equality between the moments of the Green's function and the moments of the occupied spectral function, Eq. (12) and Eq. (B1), is assured by having used the algorithmic inversion when obtaining the SOP for the Green's function. Indeed, the knowledge of the selfenergy on SOP and the use of the algorithmic inversion for solving exactly the Dyson equation ensures that the spectral function decays at least as ImΣ , thereby making the first two occupied moments (see Sec. II C) converge.
Similarly to the discussion in Sec. III A, the SOP approach combined with the algorithmic inversion allows one to follow the workflow highlighted by the green path in Fig. 4. Overall, the results presented Sec. IV are obtained using an implementation of the above approach in the heg sgm.x program of the AGWX suite [51].

C. Numerical details
Here we discuss and report the parameters that control the numerical accuracy of the quantities (polariz- FIG. 4. Flow chart representing different strategies for the calculation of the total-energy given a self-energy Σ, or a spectral function A, on a frequency grid as input for the heg sgm.x code. The strategy used in this article is highlighted with green lines. FIG. 5. Convergence study for the correlation energy per particle Ecorr, obtained with the Galitzki-Migdal formula, and using a Green function from a G0W0 calculation for the HEG at rs = 4. The parameters to converge are explained in Sec. III C. We choose to converge Ecorr for each parameter taking all the others fixed at the converged (second to last point) value. For each different parameter, we increase stepby-step its value by 20% in the convergent direction. ability, self-energy, total energy) computed by means of Eqs. (25), (27) and (28). In practice, this corresponds to going from left to right in the flow diagram of Fig. 3 following the green path, performing all calculations mentioned in the boxes. The first quantity to be computed is the polarizability P (q, ω). For each momentum q and frequency ω, we perform the integral of Eq. (25). As k in the integral is limited by k f (see Sec. III A), the discretization of the k-and x-grids, ∆k P and ∆x P , has to be converged to the zero spacing limit. Also, it is necessary to converge to zero the spacing of the momentum and frequency points of the polarizability-(q, ω) grid, controlled by ∆q and ∆ω P , along with the grid-upper limits (to infinity) q max and ω max P .
Moving to the central part of the flow chart in Fig. 3, the SOP representation of the polarizability is obtained following the method of Sec. II B, and placing the center of the 2 nd order Lorentzians on the mid points of the frequency grid, which improves the accuracy of the fit as ∆ω P → 0. Next, we employ the algorithmic-inversion method to go from the SOP representation of the po-larizability to the SOP of the screened-potential W (exact to machine precision, see Sec. II E). Using the SOP representation of W (and of G), the self-energy integral (right part of Fig. 3), Eq. (27), is formally identical to the integral in Eq. (25) for the polarizability. Therefore, the remaining parameters to converge are ∆x Σ , ∆k, ∆ω Σ , k max , and ω max Σ (using the same notation adopted above). As for the screened-potential W , we obtain the SOP representation of the self-energy following Sec. II B, and placing 2 nd -order Lorentzians on the mid points of the frequency grid. Finally, we obtain the SOP representation of the Green's function employing the algorithmicinversion method.
In principle, for each computed quantity which depends on the Green's function G, e.g. via the spectral function or its integrals, we should study the numerical stability of the computational procedure with respect to all the above parameters. Our numerical approach allows for the evaluation of the Green's function and the related spectral quantities on the real-axis, which are then used for the computation of thermodynamic quantities. In this work, we choose to converge the total energy (as obtained in Sec. III C), which is sensitive enough to guarantee a reasonable convergence for the other (spectral) properties of interest here. By changing individually each parameter (increase or decrease by 20% of its value towards convergence), we study the stability of the total energy against the selected parameter, keeping the values of all the others fixed at a reference point (baseline calculation of Fig. 5). Each target parameter is then converged separately until a plateau for the subsequent values of the computed quantity is observed. We evaluate the error on the result considering the two most distant values among those in the plateau. Importantly, it is possible to reduce the number of parameters to converge from 13 to 5, by linking all the gridspacing and broadening parameters together into a single variable, ∆, which ensures convergence for ∆ → 0 + . Specifically, we bind those parameters together by setting ∆ = ∆k P = 5∆x P = 1 6 ∆ω W = 1 9 ∆q = 1 25 ∆ω Σ = ∆x Σ = 1 3 ∆k Σ = 5 4 δ P = 1 100 δ Σ . Together with ∆, the grid-limit parameters are converged separately, following the strategy designed above. The converged values obtained for all the calculated densities are: ∆ = 0.004 k f ,

IV. RESULTS
In this Section we discuss the results obtained applying the SOP approach to the case of the one-shot G 0 W 0 calculation in the HEG. First we extensively discuss the r s = 4 case, also one of the most studied in the literature, then in Sec. IV C we provide the results for densities ranging from r s = 1 to r s = 10.

A. Spectral propagators on the real axis
We start by considering the independent particle polarizability P 0 (q, ω) computed at the G 0 level. In Fig. 6 we compare the imaginary part of P 0 , calculated using Eq. (25) and represented on SOP (fitted to 2 nd order Lorentzians with NNLS and then evaluated on a frequency grid, see Secs. II B and II D), with its analytic expression [52] (note that this is the only analytic result we use as a check -all others are evaluated numerically). The δ → 0 + broadening used in G 0 in order to converge the momentum integration does not sensibly affect the calculations. It is worth noting that the use of 2 nd order Lorentzians with respect to simple Lorentzians eases this convergence, providing for the same δ and k-grid spacing better agreement with the analytic result at δ = 0 (thermodynamic limit, see Sec. II A). From the plot comparison we can qualitatively infer that the SOP approach, together with its numerical implementation, is working Next, we look at the self-energy numerical procedures by examining directly the G 0 W 0 spectral function as shown in Fig. 7. This is obtained evaluating Eq. (27), representing the self-energy on SOP with 2 nd order Lorentzians, using the algorithmic inversion for the self-energy, and then evaluating the Green's function on a frequency grid. Focusing the attention on the lower satellite as well as on the quasi-particle band, we can see that Fig. 7 compares well with [30,53] (note that, at variance with [30], we use a logarithmic scale to represent the intensity of the spectral function, in order to highlight its structure). The plasmaron peak [53] is very visible for small momenta where the quasi-particle band broadens, while the satellite band in the occupied-frequency range (ω < µ) is sharper. As k approaches k f , the plasmaron disappears and the quasi-particle band becomes more peaked. At k = k f the spectral function presents the typical metallic divergence along the quasi-particle band, and occupied and empty satellites are almost of the same weight, in agreement with Ref. [32]. For k > k f satellites coming from empty states (ω > µ) become dominant along with the quasi-particle band, and no structure resembling a plasmaron hole appears.

B. Frequency integrated quantities and thermodynamics
We now study convergence and stability of the total energies. Following the prescription of Sec. III, we use the spectral function on SOP obtained in Sec. IV A, Eq. (12) to get analytically the occupation number n k and the occupied-band energy k (see Sec. III B), and finally numerically integrate the momenta of Eq. (28) to obtain the . total energy. To perform the convergence study on the total energy, we follow the approach described in Sec. III C which consists in converging all parameters for the calculation separately. Being the HEG a metal, the use of the algorithmic-inversion method to get a spectral function that obeys all sum rules (implied by the Dyson equation, see Sec. II E for details), including the normalization condition for the spectral function, is crucial for obtaining well-converged results. Indeed, the Luttinger discontinuity of n k makes the value of the total energy from the Galitzki-Migdal very sensitive to the converging parameters.
In Fig. 5 we show the convergence study for the correlation energy per particle (total energy minus Fockexchange): the convergence value for r s = 4 is 0.0381 ± 0.0003 Ha in agreement with Refs. [26] (with a difference of 0.0003 Ha), where calculations were done along the imaginary axis. In panel a) of Fig. 8 we plot n k , and in panel b) k (as defined in Sec. III B). The occupation number n k presents a sharp Luttinger discontinuity, which indicates that the broadening used in Eq. (27) is well-controlled and does not spoil the quality of the results. In panel c) of Fig. 8 we plot the total-energy resolved over k-contributions e k [rhs of Eq. (28)]. As previously mentioned, due to the presence of the Luttinger discontinuity, this function is sharp and thus difficult to integrate, at variance, e.g. with the RPA-Klein-energy functional, which is expected to be smoother [54].  [55]), using the data of Table I, and the covariance matrix of the fit. The fitted function is plotted in Fig. 11.

C. G0W0 for a broad range of HEG densities
In this Section we report results for the HEG with r s ranging from 1 to 10 studied at the G 0 W 0 level, following the same approach used for r s = 4. In Fig. 9 we show the computed data for the spectral function obtained with the AIM-SOP approach. In the chosen units ( f for the energy and k f for the momentum) the spectral function for increasing r s shows an increase in the separation between the quasi-particle band and the satellite occupied and empty bands. Indeed, in these units r s controls the interaction strength -see Eq. (3.24) of Ref. [52]-with the limits of the non-interacting gas obtained for r s → 0 and the strongly interacting gas corresponding to r s → ∞. Accordingly, the plasmaron peak of the satellite band at small momenta is weakened for smaller r s . The same behaviours can be observed for the occupation factor of Fig. 10 for the different densities. For r s → 0 the HEG approaches the non-interacting limit and the occupation number drops from 1 to 0 for increasing k/k f . Going toward r s = 10 the jump becomes smaller, since the quasi-particle is reduced due to the more evident satellite bands, as it can be seen from Fig. 9.
In Table I we report the corresponding total energies computed at the different densities, together with the available (to our knowledge) results in the literature. Since the calculations of Ref. [26] were done on the imaginary axis, we shall consider those as the most accurate for the comparison. We refer to the Supplemental Material [57] for the convergence studies of the total energies for the different densities. We find at r s = 1 the largest discrepancy (0.0059 Ha) with respect to the data of Ref. [26]. This can be rationalized by noting, e.g., that n k is a steeper function, thereby enhancing the numerical issues of the Galitzki-Migdal expression discussed in Sec. III B. To deepen the understanding of this numerical discrepancy, aside the convergence study of Fig. 1 provided in the Supplemental Material [57], we performed an additional calculation increasing the refinement parameter ∆ by 20%, aiming at increasing the accuracy in the integral grids, to target the steeper character of r s = 1. The result, 0.0736 Ha against 0.0749 Ha of Table I, is acceptable considering the error of 0.0015 Ha of Table I. Most importantly we stress that at variance FIG. 11. Correlation energy in Hartree units for several densities within the G0W0 approximation in the HEG. In green we show the results found in this work (see Sec. IV C) compared with those found in [35] in blue, and with [26] in orange. In green we plot the correlation energy fit of Eq. (30) (same functional form as in [55]) on the present (green) data. For reference, in dashed grey we add also the Quantum Monte Carlo data obtained by Ceperley and Alder [56] in the fit made by Perdew and Zunger [55].
with Ref. [26], our procedure provides not only accurate frequency-integrated quantities (e.g. the total energy), but also precise spectral properties on the real axis (key quantities for spectroscopy).
In Fig. 11 we plot the correlation energy of Table I as a function of r s , including the Perdez-Zunger (PZ) fit of the Quantum Monte Carlo (QMC) Ceperley Alder data as a reference [55,56]. We also exploit the same functional form of PZ to fit our data, providing in Table II γ, β 1 , and β 2 for the fitting function for the correlation energy of the HEG (in Hartree): together with the covariance matrix of the fit. In Fig. 11 we plot the result of the fit as a green line.

V. CONCLUSIONS
In this work we introduce the novel algorithmicinversion method on sum over poles (AIM-SOP) to handle frequency-dependent quantities in dynamical theories. Specializing to the case of many-body perturbation theory, we show that the AIM-SOP is able to provide a unified formalism for spectral and thermodynamic properties of an interacting-electron system. Expanding all frequency-dependent quantities on SOP, we use AIM-SOP to solve exactly and at all frequencies Dyson-like equations, getting analytic frequency-dependent (spectral) and frequency-integrated (thermodynamic) properties. This is allowed by the mapping of the Dyson equation to an effective Hamiltonian of dimension controlled by the number of poles in the SOP of the self-energy (see Sec. II E). The transformation of frequency-dependent quantities into SOP is performed exploiting the representation of their spectral functions on different basis sets: aside from the standard choice of a basis of Lorentzians, we introduce n-th order generalized Lorentzian basis elements (see Sec. II A) with improved decay properties. This allows for better numerical stability when transforming a propagator to SOP (see Sec. II B), improved analytic properties for calculating the thermodynamic quantities (see Sec. II C), and an acceleration of convergence to the thermodynamic limit (zero broadening and infinite k-space sampling). Also, once the SOP representation of a propagator is known, we use the Cauchy residue theorem to calculate convolutions and (occupied) moments, accessing both spectral and thermodynamic quantities (see Sec. II C).
In order to have a working example of the AIM-SOP approach, we apply it to the paradigmatic case of manybody perturbation theory at the G 0 W 0 level for the HEG at several densities (r s from 1 to 10). Using AIM-SOP, we are able to provide accurate spectra simultaneously with precise frequency-integrated quantities (e.g. occupation numbers and total energies). At the available densities, we find very good agreement with Refs. [30,53] for the spectral function. Moving to the total energy, we provide an in depth study of the stability and convergence of our results, finding quantitative agreement with Ref. [26] for the available r s , where calculations are performed on the imaginary axis.
Although in this article we study a homogeneous system as test case, the AIM-SOP approach aims to treat realistic non-homogeneous systems in the more general framework of dynamical embedding theories, for a full-frequency representation of potentials and propagators, the flexibility for self-consistent calculations, and the exact solution of Dyson-like equations. In this Appendix we obtain the SOP representation of a Green's function from an n-th order Lorentzian spectral function. Recalling Sec. II A, the discrete time-ordered Hilbert transform (Eq. (5)) of a (not normalized) n-th order Lorentzian,

VI. ACKNOWLEDGMENTS
induces a SOP representation for the Green's function, see Sec. II A. The expression in Eq. A1 can be computed using the residue theorem. Closing the contour in the upper/lower plane for ω ≶ µ, the poles of the integrand ζ j,m = j + e i π 2n (1+2m) δ j come only from the spectral function A. Using L'Hôpital's rule, the residues of the integrand are reduced to . (A2) Thus, taking the limit for C on the real-axis, poles and residues of the SOP for G are those in Eqs. (10) and (9). The normalization of the n-th order Lorentzian is given by summing α m of Eq. (9) and using the geometric sum, (A3) Appendix B: Moments of a propagator and occupied moments of its spectral function In this Section we discuss the equality between the (regularized) moments of a propagator Eq. (12) and the occupied moments of its spectral function. For simplicity of notation we restrict to the case of a single n-th Lorentzian L n δ , as defined in Eq. (7), and focus on the m = 2(n − 1) case [again here we suppose the integral in Eq. (12) converges which is assured by m = 0 and m = 1, but must be stronger regularized for higher degrees]: where A(ω) is the spectral function of G. To go from the second to the third line, we used the 1/ω 2n decay of the n-th Lorentzian, and applied the dominated convergence theorem which allows for the 0 + limit to be performed inside the integral. The same derivation holds for lower degree moments. For the higher order moments, m > 2(n − 1) it is not possible to discard the e iω0 + factor in the integral, thus E m>2(n−1) [G] becomes complex. The equality between E m>2(n−1) [G] [first and second line of Eq. (B1)] and the occupied moments of A [third line of (B1)] is lost, with the integral for the occupied moments of A diverging. The divergence happens because we cannot exchange the limit of the finite representation (controlled by δ i ) and the lower bound a → −∞ of the integral µ a dω ω 2(n−1) A(ω). Numerically, this translates into performing the two limits in order, i.e. fixing the lower bound of the integral and controlling the integral stability for δ i → 0, then lower a and again convergence the result for δ i → 0, and repeat until both convergences are achieved. In this continuous limit for the represen-

Supplemental Material
A unified Green's function approach for spectral and thermodynamic properties from algorithmic inversion of dynamical potentials This manuscript contains the supplemental material for the paper: "A unified Green's function approach for spectral and thermodynamic properties from algorithmic inversion of dynamical potentials".
FIG. 1. Correlation energy Ecorr convergence study obtained with the Galitzki-Migdal formula using a Green function from a G0W0 calculation for the HEG at rs = 1. See Fig. 5 of [1] for further reference.

I. CONVERGENCE STUDIES AT SEVERAL DENSITIES OF THE HEG
In this section we show the convergence studies at all densities of [1]. Following the method presented in Secs. III, and particularly III C, of [1] at r s = 4, we study stability and convergence of the correlation (total minus Fock) energy per particle at r s from 1 to 10. In Figs. 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 we plot the convergence study at the different densities.