Pierre-François Loos

Antoine Marie

Abdallah Ammar

Cumulant Green's function methods for molecules

The cumulant expansion of the Green's function is a computationally efficient beyond-GW approach renowned for its significant enhancement of satellite features in materials. In contrast to the ubiquitous GW approximation of many-body perturbation theory, ab initio cumulant expansions performed on top of GW (GW+C) have demonstrated the capability to handle multi-particle processes by incorporating higher-order correlation effects or vertex corrections, yielding better agreements between experiment and theory for satellite structures. While widely employed in condensed matter physics, very few applications of GW+C have been published on molecular systems. Here, we assess the performance of this scheme on a series of 10-electron molecular systems (Ne, HF, H 2 O, NH 3 , and CH 4 ) where full configuration interaction estimates of the outer-valence quasiparticle and satellite energies are available.

I. INTRODUCTION

The cumulant expansion is a versatile mathematical and theoretical tool that finds applications in various fields of physics. [START_REF] Kubo | Generalized cumulant expansion method[END_REF] From a mathematical point of view, the cumulant expansion is an alternative to the moments for characterizing a probability distribution function. It is often employed to obtain fluctuation and/or correlation functions beyond the mean-field approximation and is particularly valuable in dealing with many-body phenomena via the inclusion of higher-order effects. It has been shown to be useful for understanding a wide range of physical phenomena across different energy scales, from the quantum realm to the cosmological scales. For example, cumulant-based approaches have been developed and applied in quantum field theory, 2 statistical physics, [START_REF] Evans | Cumulant expansion of a Fokker-Planck equation: Rotational and translational motion in dense fluids[END_REF] plasma physics, [START_REF] Cherif | Higher order cumulants in colorless partonic plasma[END_REF] optics, [START_REF] Schack | Moment hierarchies and cumulants in quantum optics[END_REF] photonics, [START_REF] Sánchez-Barquilla | Cumulant expansion for the treatment of light-matter interactions in arbitrary material structures[END_REF] cosmology, [START_REF] Garny | Perturbation theory with dispersion and higher cumulants: Framework and linear theory[END_REF] and many others.

In condensed matter physics, the cumulant expansion has been particularly fruitful and is often used in the framework of Green's function theory. [START_REF] Hedin | On correlation effects in electron spectroscopies and the GW approximation[END_REF][START_REF] Onida | Electronic excitations: Densityfunctional versus many-body green's function approaches[END_REF][START_REF] Martin | Interacting Electrons: Theory and Computational Approaches[END_REF] Cumulant-based Green's function methods have been instrumental in providing a better description of satellite peaks in the context of photoemission spectroscopy of materials. [START_REF] Aryasetiawan | Multiple Plasmon Satellites in Na and Al Spectral Functions from Ab Initio Cumulant Expansion[END_REF][START_REF] Aryasetiawan | The gw method[END_REF][START_REF] Vos | Determination of the energy-momentum densities of aluminium by electron momentum spectroscopy[END_REF][START_REF] Vos | Electron correlation effects in the spectral momentum density of graphite[END_REF][START_REF] Vos | Quantitative measurement of the spectral function of aluminum and lithium by electron momentum spectroscopy[END_REF][START_REF] Kheifets | Spectral properties of quasiparticles in silicon: A test of many-body theory[END_REF][START_REF] Guzzo | Valence Electron Photoemission Spectrum of Semiconductors: Ab Initio Description of Multiple Satellites[END_REF][START_REF] Lischner | Physical Origin of Satellites in Photoemission of Doped Graphene: An Ab Initio G W Plus Cumulant Study[END_REF][START_REF] Gatti | Dynamical screening in correlated metals: Spectral properties of srvo 3 in the gw approximation and beyond[END_REF][START_REF] Guzzo | Multiple satellites in materials with complex plasmon spectra: From graphite to graphene[END_REF][START_REF] Kas | Cumulant expansion of the retarded one-electron Green function[END_REF][START_REF] Caruso | Band Structures of Plasmonic Polarons[END_REF][START_REF] Zhou | Dynamical effects in electron spectroscopy[END_REF][START_REF] Lischner | Satellite band structure in silicon caused by electron-plasmon coupling[END_REF][START_REF] Nakamura | Ab initio gw plus cumulant calculation for isolated band systems: Application to organic conductor (TMTSF) 2 pf 6 and transition-metal oxide srvo 3[END_REF][START_REF] Gumhalter | On the combined use of GW approximation and cumulant expansion in the calculations of quasiparticle spectra: The paradigm of Si valence bands[END_REF][START_REF] Kas | Particle-hole cumulant approach for inelastic losses in x-ray spectra[END_REF][START_REF] Verdi | Origin of the crossover from polarons to fermi liquids in transition metal oxides[END_REF][START_REF] Zhou | Cumulant Green's function calculations of plasmon satellites in bulk sodium: Influence of screening and the crystal environment[END_REF] Compared to the GW approximation of many-body perturbation theory, [START_REF] Aryasetiawan | The gw method[END_REF][START_REF] Hedin | New method for calculating the one-particle Green's function with application to the electron-gas problem[END_REF][START_REF] Reining | The GW approximation: Content, successes and limitations: The GW approximation[END_REF][START_REF] Golze | The gw compendium: A practical guide to theoretical photoemission spectroscopy[END_REF][START_REF] Marie | The GW Approximation: A Quantum Chemistry Perspective[END_REF] the cumulant form yields better agreements between experimental observations a) Electronic mail: loos@irsamc.ups-tlse.fr b) Electronic mail: amarie@irsamc.ups-tlse.fr c) Electronic mail: aammar@irsamc.ups-tlse.fr and theoretical predictions for satellite structures, effectively reproducing the series of multiple satellites observed, for instance, in the photoemission spectrum of sodium [START_REF] Aryasetiawan | Multiple Plasmon Satellites in Na and Al Spectral Functions from Ab Initio Cumulant Expansion[END_REF][START_REF] Zhou | Dynamical effects in electron spectroscopy[END_REF][START_REF] Zhou | Cumulant Green's function calculations of plasmon satellites in bulk sodium: Influence of screening and the crystal environment[END_REF] and silicon. [START_REF] Kheifets | Spectral properties of quasiparticles in silicon: A test of many-body theory[END_REF][START_REF] Guzzo | Valence Electron Photoemission Spectrum of Semiconductors: Ab Initio Description of Multiple Satellites[END_REF][START_REF] Caruso | Band Structures of Plasmonic Polarons[END_REF][START_REF] Lischner | Satellite band structure in silicon caused by electron-plasmon coupling[END_REF][START_REF] Gumhalter | On the combined use of GW approximation and cumulant expansion in the calculations of quasiparticle spectra: The paradigm of Si valence bands[END_REF][START_REF] Caruso | Spectral fingerprints of electron-plasmon coupling[END_REF][START_REF] Vlček | Quasiparticle spectra from molecules to bulk[END_REF] The cumulant expansion has also been succesfully employed to model x-ray photoemission spectra that probe core electrons. [START_REF] Nozières | Singularities in the X-Ray Absorption and Emission of Metals. III. One-Body Theory Exact Solution[END_REF][START_REF] Becrstedt | On the Theory of Plasmon Satellite Structures in the Photoelectron Spectra of Non-Metallic Solids I. Soluble Model[END_REF][START_REF] Bechstedt | Theory of Core Excitons in Semiconductors[END_REF][START_REF] Bechstedt | Electronic Relaxation Effects in Core Level Spectra of Solids[END_REF][START_REF] Kas | Green's function methods for excited states and x-ray spectra of functional materials[END_REF] This success can be understood thanks to the close connection of the cumulant ansatz with electron-boson models (see below). [START_REF] Hedin | On correlation effects in electron spectroscopies and the GW approximation[END_REF][START_REF] Lundqvist | Characteristic structure in core electron spectra of metals due to the electron-plasmon coupling[END_REF][START_REF] Langreth | Singularities in the X-Ray Spectra of Metals[END_REF][START_REF] Gunnarsson | Corrections to Migdal's theorem for spectral functions: A cumulant treatment of the time-dependent Green's function[END_REF][START_REF] Hedin | Effects of Recoil on Shake-Up Spectra in Metals[END_REF] Indeed, satellites are many-body electronic excitations that go beyond the singleparticle picture and the inclusion of higher-order correlation effects is required. [START_REF] Cederbaum | Application of Green's functions to excitations accompanying photoionization in atoms and molecules[END_REF][START_REF] Mejuto-Zaera | Are multi-quasiparticle interactions important in molecular ionization?[END_REF] In this context, cumulants are employed to approximate the higher-order terms (or vertex corrections) in the expansion of the one-body Green's function G. Developments around ab initio cumulant expansions are still ongoing and constitute an area of active research, [START_REF] Tzavala | Nonlinear response in the cumulant expansion for core-level photoemission[END_REF][START_REF] Cudazzo | Correlation satellites in optical and loss spectra[END_REF][START_REF] Cudazzo | First-principles description of the exciton-phonon interaction: A cumulant approach[END_REF][START_REF] Kas | Ab Initio Multiplet-Plus-Cumulant Approach for Correlation Effects in X-Ray Photoelectron Spectroscopy[END_REF] especially for extensions to strongly correlated materials. [START_REF] Biermann | Electronic polarons, cumulants and doubly dynamical mean field theory: Theoretical spectroscopy for correlated and less correlated materials[END_REF] The fundamental concept behind the cumulant expansion is an exponential ansatz for the one-body Green's function in the time domain

G(t) = G 0 (t)e C(t) (1)
where G 0 (t) represents a reference Green's function and C(t) is the so-called cumulant. This exponential ansatz shares obvious similarities with coupled-cluster (CC) theory [START_REF] Shavitt | Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory[END_REF] and quantum Monte Carlo. [START_REF] Foulkes | Quantum monte carlo simulations of solids[END_REF][START_REF] Austin | Quantum monte carlo and related approaches[END_REF] The cumulant expansion generates a Poisson series of satellites in the spectral function A(ω) = π -1 |Im G(ω)|, which establishes a direct link to experimental photoemission spectra via its connection to the photocurrent. [START_REF] Hedin | On correlation effects in electron spectroscopies and the GW approximation[END_REF][START_REF] Guzzo | Valence Electron Photoemission Spectrum of Semiconductors: Ab Initio Description of Multiple Satellites[END_REF][START_REF] Guzzo | Multiple satellites in materials with complex plasmon spectra: From graphite to graphene[END_REF][START_REF] Zhou | Dynamical effects in electron spectroscopy[END_REF][START_REF] Zhou | Cumulant Green's function calculations of plasmon satellites in bulk sodium: Influence of screening and the crystal environment[END_REF][START_REF] Kas | Green's function methods for excited states and x-ray spectra of functional materials[END_REF][START_REF] Almbladh | Handbook on synchrotron radiation[END_REF] In practice, the central component of the cumulant is the GW self-energy, aligning the procedure and computational cost of a cumulant calculation with that of GW. The GW+cumulant approach (henceforth GW+C) can thus be seen as an economical post-treatment that goes beyond the GW approximation.

The cumulant has been also employed in realistic molecular systems [START_REF] Vlček | Quasiparticle spectra from molecules to bulk[END_REF][START_REF] Vila | Real-Time Coupled-Cluster Approach for the Cumulant Green's Function[END_REF][START_REF] Vila | Equation-of-Motion Coupled-Cluster Cumulant Green's Function for Excited States and X-Ray Spectra[END_REF] but much less than in solids. Our goal here is to assess how this approach performs in the context of molecular systems, where highly-accurate reference data for outervalence quasiparticle and satellite energies are available. [START_REF] Marie | Reference energies for valence ionizations and satellite transitions[END_REF] Our manuscript is organized as follows. In Sec. II, we discuss several key developments and applications of the cumulant expansion in condensed matter physics before delving into the definition of the Green's function in Sec. III and the derivation of the cumulant expansion in Sec. IV. Section V reports a detailed derivation of the cumulant expansion based on the GW self-energy. In Sec. VI, we compute quasiparticle and satellite energies on a series of 10-electron molecular systems (Ne, HF, H 2 O, NH 3 , and CH 4 ), as well as their corresponding spectral functions. Various comparisons are proposed to gauge the quality of these physical quantities. Our conclusions are drawn in Sec. VII

II. SHORT REVIEW

Fundamentally, the cumulant approach is rooted in electronboson (or polaron) models where one or more fermions are coupled to bosons. In 1970, Langreth, [START_REF] Langreth | Singularities in the X-Ray Spectra of Metals[END_REF] building upon the work of Nozieres and De Dominicis, [START_REF] Nozières | Singularities in the X-Ray Absorption and Emission of Metals. III. One-Body Theory Exact Solution[END_REF] studied a simple electron-boson model where a deep core electron is coupled with bosonic excitations. He successfully provided the exact solution to this model, the spectral function revealing a quasiparticle peak accompanied by a series of satellites following a Poissonian distribution, the corresponding Green's function having the form of Eq. ( 1). [START_REF] Langreth | Singularities in the X-Ray Spectra of Metals[END_REF] This forms the basis of the cumulant ansatz.

One year before this, Lundqvist had shown that the approximate solution of this model, correct up to second order, is closely linked to the GW approximation. [START_REF] Lundqvist | Characteristic structure in core electron spectra of metals due to the electron-plasmon coupling[END_REF] The exact solution clearly evidences that the GW approximation provides a good description of the quasiparticle peak but a poor description of the satellite region, GW predicting one broad, wrongly-placed peak for the incoherent part of the spectrum. This so-called plasmaron, initially thought to be a novel type of quasiparticle excitation resulting from strong coupling between electrons and plasmons, [START_REF] Lundqvist | Characteristic structure in core electron spectra of metals due to the electron-plasmon coupling[END_REF] was later attributed to a spurious solution of the Dyson equation or an artifact introduced by GW which disappears at higher levels of theory.

The plasmaron was actually first observed in the uniform electron gas (UEG) by Hedin, Lundqvist, and coworkers, [START_REF] Hedin | Many-body effects in soft x-ray emission in metals[END_REF][START_REF] Lundqvist | Single-particle spectrum of the degenerate electron gas: I. the structure of the spectral weight function[END_REF][START_REF] Lundqvist | Single-particle spectrum of the degenerate electron gas: Ii. numerical results for electrons coupled to plasmons[END_REF][START_REF] Lundqvist | Single-particle spectrum of the degenerate electron gas: Iii. numerical results in the random phase approximation[END_REF][START_REF] Lundqvist | Single-particle spectrum of the degenerate electron gas iv. ground state energy[END_REF] before being identified in the spectrum of core electrons. [START_REF] Lundqvist | Characteristic structure in core electron spectra of metals due to the electron-plasmon coupling[END_REF] Despite claims and reports asserting the observation of the plasmaron, [START_REF] Bostwick | Observation of Plasmarons in Quasi-Freestanding Doped Graphene[END_REF][START_REF] Dial | Observations of plasmarons in a two-dimensional system: Tunneling measurements using time-domain capacitance spectroscopy[END_REF] recent consensus suggests its non-existence in materials. The spurious prediction of the plasmaron highlights GW's limitations in describing plasmon satellites. We refer the interested reader to Ref. 66 for an exhaustive discussion about the plasmaron.

Langreth's polaron model provided crucial theoretical insights and found early applications in core-level spectroscopy. Building upon Langreth's work, Hedin 8,[START_REF] Hedin | Effects of Recoil on Shake-Up Spectra in Metals[END_REF][START_REF] Hedin | Properties of electron self-energies and their role in electron spectroscopies[END_REF] went on to gener-alize the electron-boson model, occasionally termed a quasiboson model, to the valence region. This extension includes dispersion or recoil effects, enhancing the model's applicability. Interestingly, Hedin's model shares close similarities with the electron-boson Hamiltonian used by Tölle and Chan [START_REF] Tölle | Exact relationships between the GW approximation and equation-of-motion coupled-cluster theories through the quasi-boson formalism[END_REF] to highlight the connections between GW and CC theory (see also Refs. 69 and 70).

In 1996, Aryasetiawan et al. [START_REF] Aryasetiawan | Multiple Plasmon Satellites in Na and Al Spectral Functions from Ab Initio Cumulant Expansion[END_REF] made a convincing comparison between the experimental valence photoemission spectra of sodium and aluminum and the spectra simulated via the ab initio cumulant expansion. The shortcomings of GW became evident, generating only one plasmon satellite at a considerably larger binding energy. In contrast, the cumulant expansion remarkably enhanced the spectral function, notably revealing multiple plasmon satellite structures and improving the satellite positions to align closely with experimental observations. However, the improvement in intensities was comparatively modest. One reason behind this observation is the absence of considerations for extrinsic and interference effects, which predominantly impact intensities. This study likely marks the first exploration of first principles cumulant-based calculations that specifically address satellite features, moving beyond the paradigmatic UEG model. Despite the cumulant's exactness for core ionizations, [START_REF] Langreth | Singularities in the X-Ray Spectra of Metals[END_REF][START_REF] Almbladh | Handbook on synchrotron radiation[END_REF] the study reveals its surprising effectiveness for valence electrons.

One year later, Holm and Aryasetiawan conducted a thorough study on the impact of self-consistency within the cumulant approach applied to the UEG. [START_REF] Holm | Self-consistent cumulant expansion for the electron gas[END_REF] Comparing the selfconsistent (sc) versions of GW and GW+C revealed that (i) scGW and scGW+C yield similar quasiparticle energies with a small impact of the self-consistency, (ii) while the effects on satellite positions is marginal at the scGW+C level, there are notable modifications in intensity, and (iii) the improvement brought by self-consistency is more significant at the GW level, the satellite structure becoming more realistic. In short, despite the small impact on satellites, the study highlights the nuanced effects of self-consistency within the cumulant approach, especially concerning the enhancement of plasmon satellite structure.

Following a period of relative quiet, Lucia Reining's and John Rehr's groups revived the general interest in GW+C, applying this approach to semiconductors, specifically silicon. [START_REF] Guzzo | Valence Electron Photoemission Spectrum of Semiconductors: Ab Initio Description of Multiple Satellites[END_REF] The spectral function produced by the cumulant expansion revealed multiple satellites in the valence band photoemission spectrum that eluded accurate description by GW alone. This study reports the estimation of extrinsic, intrinsic, and interference effects, showcasing impressive agreement between theory and experiment. Indeed, to compare theoretical calculations with experimental data, the intrinsic spectral function alone is insufficient. Extrinsic losses, stemming from the scattering of the outgoing electron on its way to the detector, and interference effects between extrinsic and intrinsic contributions have to be considered. Notably, Ref. 17 highlights the shortcomings of GW that were attributed to the plasmaron (see above). They also suggested that GW might be more effective in cases where a sharp plasmaron is not formed, such as in the context of local plasmon structures within strongly correlated materials.

In a subsequent investigation, detailed in Ref. 20, the group extended their work to graphite, maintaining a similar approach to Ref. 17. Their data revealed multiple satellite replicas of intrinsic origin, augmented by extrinsic losses. Note that, in this material possessing more than one significant plasmon, there was no manifestation of a plasmaron at the GW level.

In lower-dimensional materials such a doped graphene [START_REF] Lischner | Physical Origin of Satellites in Photoemission of Doped Graphene: An Ab Initio G W Plus Cumulant Study[END_REF] and the two-dimensional version of the UEG, [START_REF] Lischner | Satellite structures in the spectral functions of the two-dimensional electron gas in semiconductor quantum wells: A G W plus cumulant study[END_REF] Steven Louie's group showed that GW also predicts a plasmaron while GW+C demonstrates commendable accuracy in reproducing satellite features.

In 2014, Kas, Rehr, and Reining introduced a variation of the cumulant expansion employing the retarded Green's function instead of the time-ordered Green's function, with a first application on the UEG. [START_REF] Kas | Cumulant expansion of the retarded one-electron Green function[END_REF] In the time-ordered ansatz of the cumulant, there is a decoupling between electron and hole branches, yielding "asymmetric" spectral functions. The separation of the electron and hole branches proves justified for core levels. However, this becomes a significant limitation near the Fermi level. Overcoming this limitation is challenging, as attempts to surpass the basic cumulant often result in issues, such as negative spectral functions. [START_REF] Gunnarsson | Corrections to Migdal's theorem for spectral functions: A cumulant treatment of the time-dependent Green's function[END_REF] Another interesting paper from Reining's and Rehr's groups was published one year later. [START_REF] Zhou | Dynamical effects in electron spectroscopy[END_REF] It deals with dynamical effects from a general perspective, focusing on the generation of new structures, such as satellites, arising from Coulomb interactions resulting from the coupling of excitations. The study presents a unified derivation of GW and GW+C based on the equation-ofmotion (EOM) formalism and reports a specific examination of bulk sodium in both valence and core regions, emphasizing the crucial role of self-consistency, particularly in the core region. The study demonstrates good agreement with experimental results when the intrinsic spectral function is adjusted for extrinsic and interference effects (see above). Although suitable for electron-hole satellites, the authors also point out the limitations GW+C for hole-hole satellites, as observed in nickel. [START_REF] Springer | First-principles T-matrix theory with application to the 6 ev satellite in ni[END_REF] Interestingly, they outlined how one can apply the cumulant expansion to the two-body Green's function. This extension of the cumulant approach was further investigated and developed by Cudazzo and Reining to describe phenomena like double plasmon satellites or exciton-exciton coupling within the Bethe-Salpeter formalism. [START_REF] Cudazzo | Correlation satellites in optical and loss spectra[END_REF][START_REF] Cudazzo | First-principles description of the exciton-phonon interaction: A cumulant approach[END_REF] In 2016, McClain et al. [START_REF] Mcclain | Spectral functions of the uniform electron gas via coupled-cluster theory and comparison to the G W and related approximations[END_REF] presented a notable study comparing the spectral function obtained at the CC level on finite-size UEGs with results from GW and GW+C. One of the interesting points of this paper is the comparison of GW+C with another state-of-the-art method, CC with single and double excitations (CCSD). For the 14-electron UEG at Wigner-Seitz radius r s = 4 (which corresponds approximately to the density of the valence electron in metallic sodium), they found, based on additional density-matrix renormalization group [START_REF] Chan | The density matrix renormalization group in quantum chemistry[END_REF][START_REF] Baiardi | The density matrix renormalization group in chemistry and molecular physics: Recent developments and new challenges[END_REF] (DMRG) calculations and the inclusion of the triple excitations at the CC level, that CCSD performs better than both GW and GW+C. Large errors were imputed to the underlying Hartree-Fock (HF) starting point which is known to be grossly inaccurate for metallic systems. Much better results were obtained by relying on a local-density approximation (LDA) starting point. Because, by construction, GW+C produces a plasmon-replica satellite structure even for finite systems, sev-eral satellite peaks with incorrect energies and underestimated peak heights were found. The picture is quite different for the 114-electron system, where GW produces a single satellite peak too high in energy (the infamous spurious plasmaron), while GW+C@LDA and CCSD spectra were found to be qualitatively similar. However, the CCSD spectral function has a stronger quasiparticle peak with a larger spectral width, and more fine structure overall than the GW+C spectral function.

In the same timeframe, several notable developments emerged: (i) Vigil-Fowler et al. [START_REF] Vigil-Fowler | Dispersion and line shape of plasmon satellites in one, two, and three dimensions[END_REF] explored the dispersion and line shape of plasmon satellites across various systems, employing the retarded GW+C method to investigate systems with variable dimensionality (the one-dimensional UEG, twodimensional doped graphene, the three-dimensional UEG, and silicon); (ii) Reichman's group [START_REF] Mayers | Description of quasiparticle and satellite properties via cumulant expansions of the retarded one-particle Green's function[END_REF] proposed to use the improper self-energy instead of the usual proper self-energy in the framework of the retarded GW+C method; and (iii) Caruso and Giustino conducted a comprehensive analysis of the spectral function of the UEG, employing both GW and GW+C methodologies, with a specific emphasis on angle-resolved spectral functions. [START_REF] Caruso | Band Structures of Plasmonic Polarons[END_REF][START_REF] Caruso | Spectral fingerprints of electron-plasmon coupling[END_REF][START_REF] Caruso | The GW plus cumulant method and plasmonic polarons: Application to the homogeneous electron gas*[END_REF] As a follow-up of their 2015 paper, [START_REF] Zhou | Dynamical effects in electron spectroscopy[END_REF] Zhou et al. provided a few years later an insightful comparison between the timeordered and retarded cumulant expansions, with a special focus on the dispersion and intensity of the plasmon satellites in bulk sodium and the UEG. [START_REF] Zhou | Cumulant Green's function calculations of plasmon satellites in bulk sodium: Influence of screening and the crystal environment[END_REF] Although both approaches are exact for deep core electrons, the investigation reveals that small yet significant changes accumulate due to variations in the ansatz details, causing a significant shift in satellite positions. This evidences the intricate nature of satellite calculations, which were shown to be much more sensitive to those details than the quasiparticle energies. Factors such as the linear response approximation, the level of theory employed for computing the dynamical screening W, and the diagonal approximation of the self-energy were identified as particularly influential. The validity of the linear response approximation was later investigated by Tzavala et al. [START_REF] Tzavala | Nonlinear response in the cumulant expansion for core-level photoemission[END_REF] , who demonstrated how to derive GW+C beyond linear response, with the Kadanoff-Baym functional serving as the starting point for this derivation.

In 2018, Vleck presented, using a stochastic GW approach, [START_REF] Neuhauser | Breaking the Theoretical Scaling Limit for Predicting Quasiparticle Energies: The Stochastic G W Approach[END_REF][START_REF] Vlček | Stochastic GW Calculations for Molecules[END_REF] one of the very few papers applying GW+C to small molecular systems (NH 3 , PH 3 , and C 2 H 2 ), although his primary focus remains on bulk-like (i.e. large) silicon nanocrystals. [START_REF] Vlček | Quasiparticle spectra from molecules to bulk[END_REF] The comparison of the GW+C spectral functions with experimental photoemission spectra, as well as SAC-CI calculations, [START_REF] Wasada | Computational studies of satellite peaks of the inner-valence ionization of C2H4, C2H2 and H2S using the SAC CI method[END_REF][START_REF] Ishida | Outer-and inner-valence ionization spectra of NH3, PH3, and AsH3: Symmetry-adapted cluster configuration interaction general-R study[END_REF] shows qualitative agreement with a significant weight transfer from the quasiparticle peak to the satellite region.

Finally, in 2020, Kowalski, Peng, and coworkers initiated developments of a real-time EOM-CC approach for cumulants, [START_REF] Vila | Real-Time Coupled-Cluster Approach for the Cumulant Green's Function[END_REF][START_REF] Vila | Equation-of-Motion Coupled-Cluster Cumulant Green's Function for Excited States and X-Ray Spectra[END_REF][START_REF] Rehr | Equation of motion coupled-cluster cumulant approach for intrinsic losses in x-ray spectra[END_REF][START_REF] Vila | Real-time equation-of-motion CC cumulant and CC Green's function simulations of photoemission spectra of water and water dimer[END_REF][START_REF] Vila | Real-Time Equation-of-Motion CCSD Cumulant Green's Function[END_REF][START_REF] Pathak | Real-time equation-of-motion coupled-cluster cumulant green's function method: Heterogeneous parallel implementation based on the tensor algebra for many-body methods infrastructure[END_REF] with applications to core excitations in small molecular systems. It is worth mentioning that their approach goes beyond the usual linear response of the cumulant expansion and yields, thanks to the CC exponential ansatz, a nonperturbative expression for the cumulant. They found that the nonlinear terms significantly improved the results, yielding accurate core binding energies as well as a satisfactory treatment of the absolute positions of the satellites and the overall shape of their feature, especially when double excitations are included.

III. ONE-BODY GREEN'S FUNCTION

We denote as G T and G R the time-ordered and retarded Green's functions, [START_REF] Zhou | Cumulant Green's function calculations of plasmon satellites in bulk sodium: Influence of screening and the crystal environment[END_REF] and their matrix elements in the spinorbital basis are respectively defined as

G T pq (t) = (-i) ⟨Ψ N 0 | T [a p (t)a † q (0)]|Ψ N 0 ⟩ (2a) G R pq (t) = (-i)Θ(t) ⟨Ψ N 0 |{a p (t), a † q (0)}|Ψ N 0 ⟩ (2b)
where T and {•, •} are the time-ordering and anticommutation operators, respectively, Θ(t) is the Heaviside function, and |Ψ N 0 ⟩ is the exact N-electron ground-state wave function. In the following, the indices i, j, k, and l are occupied (hole) spinorbitals; a, b, c, and d are unoccupied (particle) spinorbitals; p, q, r, and s indicate arbitrary spinorbitals; and ν labels single (de)excitations. Here a p (t) and a † p (t) are annihilation and creation operators in the Heisenberg picture, i.e.,

a p (t) = e i Ĥt a p e -i Ĥt a † p (t) = e i Ĥt a † p e -i Ĥt (3)
Introducing the greater and lesser components of the Green's function,

G > pq (t) = -i ⟨Ψ N 0 |a p (t)a † q (0)|Ψ N 0 ⟩ (4a) G < pq (t) = +i ⟨Ψ N 0 |a † p (t)a q (0)|Ψ N 0 ⟩ (4b)
the time-ordered and retarded Green's functions are expressed as

G T pq (t) = Θ(+t)G > pq (t) + Θ(-t)G < pq (t) (5a) G R pq (t) = Θ(+t)G > pq (t) -Θ(+t)G < pq (t) (5b)
In frequency space, the corresponding elements of the spectral function, which is directly related to the photoemission spectrum of the system (see Sec. II), are given by

A pq (ω) = i 2π G > pq (ω) -G < pq (ω) = 1 π Im G T pq (ω) = - 1 π Im G R pq (ω) (6)
The greater and lesser components of G can be combined in various ways which are equivalent if observables are calculated consistently. [START_REF] Zhou | Cumulant Green's function calculations of plasmon satellites in bulk sodium: Influence of screening and the crystal environment[END_REF] In the time-ordered and retarded versions of the cumulant expansion, the diagonal matrix elements of G ≶ (t) share the same form

G ≶ pp (t) ∝ e -iϵ p t e C ≶ pp (t) (7)
where the ϵ p 's are the one-body energies employed to build G 0 (t) [see Eq. ( 1)] and C ≶ pp (t) are the cumulant matrix elements to be determined. However, in the time-ordered formulation, because only hole (particle) states contribute to the electron removal (addition) spectrum, we have

G > ii (t) = 0 [G < aa (t) = 0], which yields G T ii (t) = iΘ(-t)e -iϵ i t e C T ii (t) (8a) G T aa (t) = iΘ(+t)e -iϵ a t e C T aa (t) (8b)
For the retarded Green's function, hole and particle states contribute to both the electron addition and removal spectra, which explains the more symmetric definition of the following retarded Green's function

G R pp (t) = -iΘ(t)e -iϵ p t e C R pp (t) (9)
As mentioned in Sec. II, in the time-ordered formulation of the cumulant, there is a separation between the electron and hole branches, as readily seen in Eqs. (8a) and (8b), leading to spectral functions that exhibit asymmetry. While this electron-hole branch decoupling is appropriate for core levels, it poses a notable limitation when approaching the Fermi level. Therefore, in the following, our derivation is based on the retarded Green's function. As explained in Chapter 5 of Ref. 10, the poles of retarded quantities must be shifted to the lower part of the complex plane by an infinitesimal amount due to the requirement of causality. This is a major difference from the more conventional time-ordered quantities where the poles are displaced in the upper or lower plane depending on their relative position with respect to the chemical potential.

IV. THE CUMULANT EXPANSION

The derivation of the cumulant expansion can be approached through different methods: diagrammatic techniques, [START_REF] Nozières | Singularities in the X-Ray Absorption and Emission of Metals. III. One-Body Theory Exact Solution[END_REF][START_REF] Hedin | Effects of Recoil on Shake-Up Spectra in Metals[END_REF] the EOM formalism, [START_REF] Almbladh | Handbook on synchrotron radiation[END_REF] the Baym-Kadanoff equation, [START_REF] Guzzo | Valence Electron Photoemission Spectrum of Semiconductors: Ab Initio Description of Multiple Satellites[END_REF] or by imposing its form and identifying the cumulant with the Dyson equation linking the Green's function and the self-energy. [START_REF] Aryasetiawan | Multiple Plasmon Satellites in Na and Al Spectral Functions from Ab Initio Cumulant Expansion[END_REF][START_REF] Gunnarsson | Corrections to Migdal's theorem for spectral functions: A cumulant treatment of the time-dependent Green's function[END_REF] The latter stands out as the simplest and we shall follow this procedure in our derivation. Additional details are provided in the supporting information. It is also worth mentioning that, because we deal with small molecular systems, we rely on a HF reference Green's function.

By definition, the cumulant ansatz of the Green's function is

G(t) = G HF (t)e C(t) (10)
where, at the HF level, we have

G HF pq (t) = -iΘ(t)e -iϵ HF p t δ pq (11)
and the ϵ HF p 's are HF one-electron orbital energies. We start from the Dyson equation that links the Green's function to the HF Green's function via the correlation part of the self-energy

G(x 1 x 1 ′ ; t) = G HF (x 1 x 1 ′ ; t) + d(t 1 t 2 ) d(x 2 x 3 ) G HF (x 1 x 2 ; t-t 1 )Σ c (x 2 x 3 ; t 1 -t 2 )G(x 3 x 1 ′ ; t 2 ) (12)
where we have taken into account the time-translation invariance of the Hamiltonian and x is a composite variable gathering spin and space.

Expanding Eq. ( 10) as a function of C(t) and Eq. ( 12) with respect to Σ c yields

G(t) = G HF (t) + G HF (t)C(t) + • • • (13a) G(t) = G HF (t) + dt 1 dt 2 G HF (t -t 1 )Σ c (t 1 -t 2 )G HF (t 2 ) + • • • (13b)
where we have omitted, for the sake of simplicity, the spinspatial coordinates. In order to determine C(t) up to first order in W, the second terms of the right-hand side of the two previous equations are equated, which gives

G HF (t)C(t) = dt 1 dt 2 G HF (t -t 1 )Σ c (
C pq (t) = i dω 2π e -i(ω-ϵ HF p )t G HF pp (ω)Σ c pq (ω)G HF qq (ω) (16)
This last identity represents the general expression of the matrix elements of the cumulant within the linear response approximation.

Using the expression of the elements of the HF retarded Green's function in frequency space

G HF pp (ω) = lim δ→0 + 1 ω -ϵ HF p -iδ (17)
enforcing the diagonal approximation for the self-energy, i.e. Σ pq (ω) ≈ δ pq Σ pp (ω), and applying a frequency shift, one gets

C pp (t) = i lim δ→0 + dω 2π Σ c pp ω + ϵ HF p (ω + iδ) 2 e -iωt (18)
Substituting Eq. ( 18) into the cumulant ansatz gives the diagonal elements of the Green's function in the time domain

G pp (t) = -iΘ(t)e -iϵ HF p t+C pp (t) ( 19
)
The cumulant is often expressed as a function of the cumulant kernel

β p (ω) = - 1 π Im Σ c pp (ω) (20)
as

C pp (t) = dω β p ω + ϵ HF p ω 2 e -iωt + iωt -1 ( 21
)
This expression is known as the Landau form of the cumulant (see the supporting information for a detailed derivation). 88

V. GW-BASED CUMULANT EXPANSION

To derive the GW+C expression of the Green's function, let us now consider the diagonal elements of the GW retarded self-energy, which reads

Σ c pp (ω) = iν M 2 piν ω -ϵ i + Ω ν + iη + aν M 2 paν ω -ϵ a -Ω ν + iη (22)
where the ϵ p 's are potentially quasiparticle energies depending on the level of self-consistency. The transition densities are given by

M pqν = jb ⟨p j|qb⟩ X jb,ν + Y jb,ν (23)
where ⟨pq|rs⟩ are the usual two-electron integrals in Dirac notation and the matrices X and Y gather the eigenvectors of the direct RPA problem [START_REF] Ring | The Nuclear Many-Body Problem[END_REF][START_REF] Chen | Random-phase approximation methods[END_REF][START_REF] Ren | Beyond the G W approximation: A second-order screened exchange correction[END_REF] obtained via the following linear eigenvalue problem

A B -B -A • X Y Y X = X Y Y X • Ω 0 0 -Ω (24)
where the diagonal matrix Ω contains the positive RPA eigenvalues Ω ν . The matrix elements of the (anti)resonant block A and the coupling block B read

A ia, jb = (ϵ a -ϵ i )δ i j δ ab + ⟨ib|a j⟩ (25a) B ia, jb = ⟨i j|ab⟩ (25b)
Substituting the GW retarded self-energy into Eq. ( 18) and performing the frequency integration yields

C pp (t) = iν ζ piν e -i∆ piν t + i∆ piν t -1 + aν
ζ paν e -i∆ paν t + i∆ paν t -1 (26) with

ζ piν = M piν ∆ piν 2 ζ paν = M paν ∆ paν 2 ( 27
)
and

∆ piν = ϵ i -ϵ HF p -Ω ν -iη (28a) ∆ paν = ϵ a -ϵ HF p + Ω ν -iη (28b)
Therefore, defining

ϵ QP p = ϵ HF p + ∆ϵ QP p ( 29
)
where

∆ϵ QP p = - iν ∆ piν ζ piν - aν ∆ paν ζ paν = Σ c pp ω = ϵ HF p ( 30
)
is the quasiparticle shift, the diagonal elements of the Green's function in the frequency domain are given by

G pp (ω) = -iZ p dt Θ(t)e i ω-ϵ QP p t e iν ζ pqν e -i∆ piν t + aν ζ paν e -i∆paν t (31)
where

Z QP p = exp        - iν ζ piν - aν ζ paν        = exp        ∂Σ c pp (ω) ∂ω ω=ϵ HF p        (32)
is the quasiparticle weight (or renormalization factor) and the last term of Eq. ( 31) containing the double exponential is responsible for the appearance of satellites. In the case of a G 0 W 0 calculation where one linearizes the quasiparticle equation to obtain the quasiparticle energies, the GW+C renormalization factor associated with the quasiparticle peak and its GW counterpart

Z GW p = 1 1 - ∂Σ c pp (ω) ∂ω ω=ϵ HF p (33)
agree up to first order, as readily seen by comparing Eqs. ( 32) and (33). Moreover, in this very specific case, it is easy to show that Re Z QP p ≤ Re Z GW p , which evidences that the cumulant expansion systematically implies a redistribution of weights from the quasiparticle peak to the satellite structure. However, this is not always true in the general case.

Expanding the last term of Eq. ( 31) to first order, one obtains the following expression for the diagonal elements of the spectral function [see Eq. ( 6) for its definition]

A pp (ω) ≈ Z p δ ω -ϵ QP p + iν Z sat piν δ ω -ϵ sat piν + aν Z sat paν δ ω -ϵ sat paν + • • • (34)
which features two sets of satellites at energies

ϵ sat piν = ϵ QP p + ∆ piν = ∆ϵ QP p + ϵ i -Ω ν (35)
ϵ sat paν = ϵ QP p + ∆ paν = ∆ϵ QP p + ϵ a + Ω ν (36)
each located on a different branch and associated with the respective weights

Z sat piν = Z QP p ζ piν Z sat paν = Z QP p ζ paν (37)
where one can readily see that they are directly proportional to the quasiparticle spectral weight. Here, we limit our analysis to these two sets of satellite peaks (especially the satellite peaks on the hole branch) as expanding to second order would produce satellites with even smaller weights and further away from the quasiparticle peak.

VI. SPECTRAL FUNCTIONS

The GW+C scheme has been implemented in quack, [START_REF] Loos | QuAcK: a software for emerging quantum electronic structure methods[END_REF] an open-source electronic structure package for emerging methods. The present implementation of GW is described in Ref. 33. To compute the quasiparticle energies at the one-shot G 0 W 0 level, we do not linearize the self-energy, solving instead the frequency-dependent quasiparticle equation using Newton's method:

ω -ϵ HF p -Re Σ c pp (ω) = 0 (38)
In the G 0 W 0 calculations, we set η = 0.001 E h in the expression of the self-energy [see Eq. ( 22)] unless we plot the spectral function in which case we set η = 0.01 E h to broaden the peaks further (see below). Note also that, at the G 0 W 0 level, we have ϵ p = ϵ HF p in the various expressions reported in Sec. IV. The HF starting point quantities are computed in the restricted formalism. The qsGW calculations are performed with the regularized scheme based on the similarity renormalization group approach, as described in Ref. 93. A flow parameter of s = 500 is employed. In the case of self-consistent calculations, we have ϵ p = ϵ GW p . All (hole and particle) states are corrected. We systematically employed Dunning's aug-cc-pVDZ basis set for all calculations and consider below the well-known 10-electron series of molecular systems, namely, Ne, HF, H 2 O, NH 3 , and CH 4 . As reference for the outer-valence quasiparticle and satellite energies, we rely on the full configuration interaction (FCI) values reported in Ref. 58, from where the geometries have also been extracted.

At the G 0 W 0 level, satellite transition energies have been obtained via the linear version of the GW equations as described in Refs. 94-97. Within this scheme, the so-called "upfolded" matrix is diagonalized in a larger space that includes the 2h1p and 2p1h configurations, the satellite energies being obtained as higher/lower roots. Although this procedure is more computationally expensive than the usual "downfolded" version where one solves the non-linear equation (38), it is still technically feasible for the present systems and ease the obtention of satellite energies which are extremely challenging to get as solutions of the non-linear equation. Satellite energies are not computed at the qsGW level due to the static nature of this approximation.

At the GW level, the spectral function is

A GW p (ω) = - 1 π Im Σ c pp (ω) ω -ϵ HF p -Re Σ c pp (ω) 2 + Im Σ c pp (ω)
2 (39) where the expression of the elements Σ c pp (ω) is given in Eq. (22). At the GW+C level, the spectral function reads (see supporting information)

A GW+C p (ω) = A QP p (ω) + iν A C piν (ω) + aν A C paν (ω) + . . . ( 40
)
where the quasiparticle part is pressions:

A QP p (ω) = -
A C piν (ω) = - 1 π Re Z sat piν Im ϵ sat piν + Im Z sat piν ω -Re ϵ sat piν ω -Re ϵ sat piν 2 + Im ϵ sat piν 2 (42a) A C paν (ω) = - 1 π Re Z sat paν Im ϵ sat paν + Im Z sat paν ω -Re ϵ sat paν ω -Re ϵ sat paν 2 + Im ϵ sat paν 2 (42b)
The outer-valence spectral function A(ω) = p A p (ω) (where the sum is performed on the outer-valence orbitals) computed at the G 0 W 0 and G 0 W 0 +C levels for each 10-electron molecular system is represented in Fig. 1. For example, the spectral function of water features three well-defined quasipar-

G 0 W 0 qsGW G 0 W 0 +C qsGW +C -3 -2 -1 0 1 FIG. 2.
Percentage of error with respect to FCI for the outer-valence ionization energies reported in Table I computed at the G 0 W 0 , qsGW, G 0 W 0 +C, qsGW+C levels. All calculations are performed with the aug-cc-pVDZ basis and η = 0.001 E h . II computed at the G 0 W 0 , G 0 W 0 +C, and qsGW+C levels. All calculations are performed with the aug-cc-pVDZ basis and η = 0.001 E h ticle peaks around 13, 15, and 19 eV and satellite structures higher in energy. Figure 1 evidences that the main difference between the G 0 W 0 and G 0 W 0 +C spectral functions is a global shift to lower energy of the satellite features, the quasiparticle peaks being much less affected. In addition, the intensities of the satellite peaks are also slightly affected but the overall structure of the satellite band is not modified by the cumulant expansion. On the other hand, vertex corrections to the G 0 W 0 scheme have been shown to create new satellite features in ionization spectra (see Ref. 46). Tables I and II report the outer-valence ionization energies and satellite transition energies, respectively, for the present set of molecules. Most of these satellites have vanishing weights. The percentage of error with respect to the reference FCI values is represented in Figs. 2 and3. The first ionization of water corresponds to an electron removal from the orbital 1b 1 , and this process is denoted as (1b 1 ) -1 . The lowest satellite transition in H 2 O corresponds to the process (3a 1 ) -1 (1b 1 ) -1 (4a 1 ) [START_REF] Kubo | Generalized cumulant expansion method[END_REF] , where one electron is removed from orbitals 3a 1 and 1b 1 and one is attached in orbital 4a 1 . The symmetry labels of these charged excited states are also reported in Tables I andII.

G 0 W 0 G 0 W 0 +C qsGW +C -15 -
Figure 2 shows that the quasiparticle energies are not improved by the cumulant expansion. They are actually slightly worse. It is particularly true at the qsGW level where the error is increased by 1 or 2% when one considers the cumulant scheme. The spectral weights of the quasiparticle peaks reported in Table I show that the cumulant induces a small redistribution of spectral weight from the quasiparticle to the satellites but this effect is quite subtle.

Concerning the satellites (see Fig. 3), we find that G 0 W 0 +C can significantly reduce the error of G 0 W 0 in certain situations. This is the case for the 2 2 ∆ satellite of HF (from 2 eV to 0.05 eV), the 1 2 A 1 and 3 2 B 1 satellites of H 2 O (from 1.6 eV to 0.2 eV and from 2 eV to 0.6 eV), and to a lesser extent, the 2 2 A 1 satellite of NH 3 (from 0.6 eV to -0.3 eV). On the contrary, the qsGW+C scheme tends to overcorrect the satellite energies. The 2 2 P satellite of Ne, the 2 2 B 1 satellite of H 2 O, the 2 2 E satellite of NH 3 , and the 2 2 T 2 satellite of CH 4 are much harder to describe at the GW level and errors remain larger even within the cumulant scheme.

VII. CONCLUDING REMARKS

In the present article, we have reviewed the literature on cumulant Green's function methods, provided a detailed derivation of the associated equations, and investigated the performance of this scheme in the context of molecular systems where solid reference data for satellite transitions are now available. In particular, we have compared the satellite transition energies obtained via the GW+C scheme and the "upfolded" version of the GW equations. In a nutshell, G 0 W 0 +C does sometimes improve upon G 0 W 0 but it is far from being systematic. The cumulant version of qsGW has been found to usually overcorrect satellite energies. However, the cumulant approach allows us to estimate satellite energies without solving the dynamical GW equations or its larger (hence more expensive) upfolded version. These observations show that the development of new Green's function methods capable of describing accurately satellite states in molecular systems would be extremely useful.

CONFLICTS OF INTEREST

There are no conflicts to declare.

1 π 2 + 2 ( 41 )

 12241 Re Z QP p Im ϵ QP p + Im Z QP p ω -Re ϵ QP p ω -Re ϵ QP p Im ϵ QPp and the satellite contributions are given by the following ex-

8 FIG. 1 .

 81 FIG.1. Outer-valence spectral functions of the 10-electron molecular series computed at the G 0 W 0 (black) and G 0 W 0 +C (red) levels with the aug-cc-pVDZ basis and η = 0.01 E h

15 FIG. 3 .

 153 FIG.3. Percentage of error with respect to FCI for the satellite transition energies reported in Table II computed at the G 0 W 0 , G 0 W 0 +C, and qsGW+C levels. All calculations are performed with the aug-cc-pVDZ basis and η = 0.001 E h

  t 1t 2 )G HF (t 2 )(14)

	Projecting this expression in the spinorbital basis yields
	G HF pr (t)C rq (t) =	dt 1 dt 2 G HF pr (t-t 1 )Σ c rs (t 1 -t 2 )G HF sq (t 2 )
	r	rs
		(15)
	or, by moving to frequency space,

TABLE I .

 I Selection of quasiparticle energies of the 10-electron series computed at various levels of theory with the aug-cc-pVDZ basis and η = 0.001 E h . The spectral weight is reported in parenthesis.

	System	State	Process	G 0 W 0	qsGW	G 0 W 0 +C	qsGW+C	FCI
	Ne	1 2 P	(2p) -1	21.104(0.947)	21.435(0.937)	20.983(0.942)	20.733(0.930)	21.426
	HF	1 2 Σ +	(3σ) -1	15.868(0.937)	16.144(0.924)	15.740(0.931)	15.510(0.916)	16.059
		1 2 Π	(2a 1 ) -1	19.812(0.942)	20.084(0.931)	19.740(0.938)	19.497(0.926)	20.043
	H 2 O	1 2 B 1	(1b 1 ) -1	12.485(0.933)	12.640(0.920)	12.384(0.927)	12.228(0.912)	12.540
		1 2 A 1	(3a 1 ) -1	14.781(0.935)	14.932(0.921)	14.698(0.929)	14.466(0.914)	14.829
		1 2 B 2	(1b 2 ) -1	18.865(0.941)	19.069(0.931)	18.822(0.938)	18.706(0.928)	18.995
	NH 3	1 2 A 1	(3a 1 ) -1	10.837(0.933)	10.870(0.922)	10.776(0.928)	10.663(0.915)	10.762
		1 2 E 1	(1e g ) -1	16.578(0.940)	16.655(0.930)	16.544(0.936)	16.461(0.926)	16.534
	CH 4	1 2 T 2	(1t 2 ) -1	14.466(0.943)	14.446(0.936)	14.445(0.940)	14.406(0.933)	14.285

TABLE II .

 II Selection of satellite transition energies of the 10-electron series computed at various levels of theory with the aug-cc-pVDZ basis and η = 0.001 E h . The IP-EOM-CCSDT and IP-EOM-CCSDTQ results from Ref. 58 are reported for comparison purposes.

								IP-EOM
	System	State	Process	G 0 W 0	G 0 W 0 +C	qsGW+C	CCSDT	CCSDTQ	FCI
	Ne	2 2 P	(2p) -2 (3s) 1	54.398	52.168	48.259			49.349
	HF	2 2 ∆	(1π) -2 (4σ) 1	36.453	34.492	31.058	34.885	34.403	34.445
	H 2 O	2 2 B 1	(3a 1 ) -1 (1b 1 ) -1 (4a 1 ) 1	30.846	29.370	26.868	27.694	27.049	27.065
		2 2 A 1	(3a 1 ) -2 (4a 1 ) 1	28.770	27.293	24.577	27.476	27.104	27.131
		3 2 B 1	(3a 1 ) -1 (1b 1 ) -1 (4a 1 ) 1	30.867	29.387	26.881	29.129	28.729	28.754
	NH 3	2 2 A 1	(3a 1 ) -2 (4a 1 ) 1	24.410	23.510	21.657	24.101	23.818	23.829
		2 2 E	(3a 1 ) -2 (3e) 1	24.997	24.098	22.317	25.882	25.648	25.655
	CH 4	2 2 T 2	(1t 1 ) -2 (3a 1 ) 1	30.681	30.317	29.438	28.415	28.123	28.108

ACKNOWLEDGEMENTS

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 863481). Additionally, it was supported by the European Centre of Excellence in Exascale Computing (TREX), and has received funding from the European Union's Horizon 2020 -Research and Innovation program -under grant agreement no. 952165.