Are multi-quasiparticle interactions important in molecular ionization?

The quantitative understanding of electronic excitations in complex molecular extended systems is one of the most fundamental open challenges in modern theoretical chemistry. Perhaps the most direct experimental probe of the excited state manifold is given by photo-emission spectra (PES), which directly access individual electronic states^1–4 and, in principle, access information on single-particle orbitals.^5–7 Theoretical approaches are meant to provide an interpretative connection between measured spectral features and chemical concepts, thus helping design systems with tailored (opto)electronic properties. An accurate realization of this ideal often requires a treatment beyond effective one-body theories, such as Hartree–Fock (HF) or Kohn–Sham density-functional theory, to capture all interactions that renormalize single electron properties and lead to the formation of quasiparticles (QPs).

Traditional gold-standard wave function approaches, such as the configuration interaction or coupled cluster methods, can yield highly accurate predictions, and they have been extensively applied to describe the outer valence and core electron spectra of small molecules,^8–13 which fundamentally behave as individual QPs. However, the description of the inner valence region is more complicated: multiple excitation mechanisms are available, and ionization can be accompanied by simultaneous neutral excitation of the system. At the energy scales associated with inner valence PES, the single QP (SQP) picture breaks down, and the spectra exhibit a multitude of features,^14,15 which we refer to as the multi-QP (MQP) regime. In analyzing PES, these are typically referred to as shake-up satellites. Describing these ionized states requires accounting for the nontrivial interactions among many excited states. The inner valence excitations in the MQP regime thus represent an unambiguous measure of electron correlation.¹⁶ This correlation may be dynamic or static in nature and, depending on this distinction, it may be captured by different approximations: the dynamical correlation would correspond to a frequency dependent potential on top of an effective one-body theory and can (in principle) be described by perturbation theory. On the other hand, the static correlation requires an intrinsically many-body description, i.e., it cannot be qualitatively approximated by a single effective one-body Hamiltonian and its perturbation.

Ab initio wave function methods capture both types of correlations; they, in principle, describe the high energy holes but often require extensive exploration of the exponentially large Hilbert space, limiting their application to small molecular systems. However, while convenient from a computational perspective, the formulation of these methods in terms of determinants makes a definite distinction between the dynamic and static correlations complicated in most cases. In particular, it is unclear whether MQP features in PES belong to the former or latter class. On the one hand, they are intrinsically many-body effects that cannot be captured by a single Slater determinant. On the other hand, they often arise from neutral excitations weakly coupled to holes, which suggests that they should be amenable to a perturbative description.

In the pursuit of an accurate but scalable theoretical description of PES, many-body perturbation theory (MBPT) has proven to be a viable alternative for computing QP excitations. It offers a systematically improvable theoretical framework for capturing dynamical correlations. MBPT relies on physically motivated concepts such as screening, while retaining a polynomial computational scaling. In particular, the popular GW approximation is extensively applied to treat molecules and solids, and it predicts outer valence and core electron spectra with good accuracy.^17–22 With recent algorithmic developments, GW can be applied to systems with thousands of electrons.^20,23–28 Unfortunately, the GW framework is fundamentally limited to systems where classical electrodynamics dominates the electronic correlation. Hence, the description of the MQP regime is beyond the GW capabilities.^29–36 Within the MBPT framework, the missing higher-order quantum interactions are represented by the vertex term Γ, which encodes dynamical two-particle correlations. The inclusion of Γ leads to the GWΓ approach, which, for molecular systems, has been studied only on first ionization energies and electron affinities.^{30,31,37–40} However, the role of Γ on the entire valence states and the MQP excitations has not been explored up to now.

In this Communication, we demonstrate that MBPT with vertex terms successfully captures the correlations necessary for a qualitative description of the MQP regime of PES, suggesting them to be dynamic in nature. By comparison with explicitly correlated and numerically exact methods, such as the adaptive sampling configuration interaction (ASCI) approach,^41–44 we show that the GWΓ approach yields superior results compared to GW throughout all energy scales. Moreover, we discover that the non-local quantum interactions play a significant role in all valence electron excitations and correctly capture the breakdown of the SQP picture, eliminating spurious artifacts of GW in the shake-up region of PES. We choose, as a case study, the PES of selected closed shell molecules (NH₃, H₂O, CH₄, C₂H₂, and N₂), which, as we show, present rich MQP character in the inner valence region.

The spectral function A(ω) computed with ASCI and MBPT is the figure of merit. This corresponds to the physical observable (i.e., PES) and is given as the trace over the imaginary part of Green’s function (GF) G_i,j(ω),³² i.e., $A(ω)=1πTrImG(ω)$⁠. The i, j sub-indices correspond to a chosen single particle basis. The peaks in A(ω) correspond to the poles of G⁴⁵ and for a particular hole component,

where $|ΨmN$ is the mth eigenstate of the N particle system, with energy $EmN$⁠. Furthermore, c_i is the ith annihilation operator and η is an integrating factor. In the mean-field and SQP regimes, we expect only one non-vanishing overlap $ΨmN−1|ci|Ψ0N2$ per single-particle state i, corresponding to a dressed hole with energy $EmN−1−E0N$⁠. However, multiple states $|ΨmN−1$ with distinct $EmN−1$ exist in the MQP regime.

First, we evaluate Eq. (1) without further approximations via the ASCI algorithm,^41–44,46 which captures both SQP and MQP regimes, regardless of the degree of correlation of the excited states. ASCI provides accurate Green’s functions in model systems;^43,47 here, we compute Green’s functions of ab initio Hamiltonians using ASCI for the first time. The spectral features calculated from Eq. (1) represent a series of infinitely sharp peaks due to the use of finite atomic basis sets and a Hermitian Hamiltonian. As a result, no scattering states are considered and the finite lifetimes of individual excitations (arising due to the coupling to continuum) are neglected. While the MBPT calculations also discretize the state space, this is performed in the converged real space grid, which enables the description of the continuum states as well (for the details of the calculations, see the supplementary material). To facilitate comparison with the MBPT results, which thus include lifetimes due to electron–electron scattering, the ASCI-computed spectra have been artificially broadened using the η parameter in Eq. (1). For the heuristics and technical details of ASCI, including basis set extrapolation, see the supplementary material.

MBPT is based on the conceptually different approach of obtaining QP energies through the Dyson equation.³² This relates Green’s function of the fully interacting QPs G_i,j(ω) and a reference (mean-field) Green’s function G_0,i,j(ω), through the self-energy Σ_i,j(ω). The latter, in principle, contains all many-body effects. The poles of G_i,j(ω) are subject to the QP fixed point equation,

Here, $εj0$ is a pole of G₀ and Δ_j(ω) comprises the coupling due to the off-diagonal elements of Σ_i,j(ω). In this work, we adopt the common diagonal self-energy approximation, i.e., we assume Δ_j(ω) = 0 in Eq. (2). In the SQP regime, Σ(ω)_i,j merely shifts the poles of G with respect to G₀, and Eq. (2) only has one solution per orbital. For MQPs, the structure of the self-energy is necessarily more complex and, in principle, yields multiple solutions to Eq. (2). To capture such correlated states, Σ(ω)_i,j must account for the interactions between the particle–hole pairs and ionized holes (or injected electrons) in the system. It is often helpful to represent the solutions of Eq. (2) graphically as done in the discussion below.

In this work, we base the perturbation expansion on top of the Hartree–Fock (HF) Hamiltonian employing a one-step iteration in building the self-energy. In principle, the correlation self-energy, Σ_c, is derived from the equation of motion of QP Green’s function and the non-interacting reference (i.e., HF in this case). The full (self-consistently renormalized) expression reads³² 

where we employ a short-hand notation for space–time coordinates 1 ≡ (r₁, t₁) and the bar indicates a coordinate that is integrated over. Here, Σ_Hxc contains Hartree, exchange, and correlation interactions, and ν(1, 2) = δ(t₁ − t₂)/|r₁ − r₂| is the bare Coulomb interaction. Furthermore, we introduce a generalized response function connecting the induced Green’s function to a perturbing potential U(1) as ³χ(1, 2, 3) = −iδG(1, 2)/δU(3).³¹ To lower the computational cost, it is common to approximate δΣ_Hxc/δG, the two-particle interaction kernels.^32–36 In particular, taking δΣ_Hxc/δG ≈ δΣ_H/δG leads to the popular GW approximation,^48,49 which describes the correlations as induced time-dependent density fluctuations. In contrast, GWΓ evaluates Eq. (3) in full and captures the QP couplings. The additional terms (compared to GW) stemming from δΣ_xc/δG are referred to as vertex corrections, and they are responsible for mutual coupling of MQP interactions.^29–32,35 In particular, it has been shown that including vertex corrections is necessary for self-consistent renormalization⁵⁰ and to capture satellite spectral features for simple models of core electrons.⁵¹ While the computational cost of Γ is large,^30,52 a recent stochastic formalism introduced a linear-scaling algorithm, which we apply here³¹ to study the effect of vertex corrections on the satellite features in the valence spectrum.

In practice, we resort to the single step correction scheme for both GW and GWΓ. In this approximation, the two-particle kernel derives from the mean-field Hamiltonian, i.e., it becomes δΣ_H/δG for GW. For GWΓ, the kernel is δΣ_Hx/δG, i.e., it also includes the HF exchange term.³¹ Thus, the GWΓ implementation contains, on top of induced Hartree potentials, also the dynamical induced exchange interactions. The latter are introduced through an additional term in the self-energy (see the supplementary material for details) that is responsible for additional oscillations in Σ_c(ω) and, hence, for additional solutions to Eq. (2). Furthermore, as a usual process, we employ the random-phase approximation (RPA) in GW; however, it is avoided in GWΓ, where the time evolution of states is governed by the full Hamiltonian (by the time-dependent Hartree–Fock approximation in this particular case). Hence, the propagated states incorporate excitonic effects through the dynamical exchange potential.³³ Effectively, our GWΓ implementation corresponds to performing one full iteration through Hedin’s pentagon. An alternative formulation based on the T-matrix formalism allows us to incorporate other channels of correlations,^29,32,53 which would be accounted for if the Hedin’s pentagon [with Σ_xc defined by Eq. (3)] was run to self-consistency. The use of GW without RPA is possible, but it worsens the quality of A(ω).³¹ Furthermore, we deliberately choose HF as the starting point as it often lacks spurious multiple solutions to Eq. (3) ⁵⁴ and HF single-particle states are close to true Dyson orbitals.⁵⁵ 

The computed spectral functions with different theoretical approaches are illustrated in the top panel of Fig. 1 for the CH₄, NH₃, and H₂O molecules. The frontier orbitals appear at the lowest absolute frequencies, and thus we, henceforth, refer to the orbitals with energy close to zero frequency as low energy, or outer valence, and those with energy far away from zero as high energy, or inner valence. The distinction between the SQP and MQP regimes is evident: the highest occupied molecular orbitals (HOMO) and states energetically close to these are composed of a single sharp peak per orbital, consistent with the SQP picture. In contrast, excitations far away from the HOMO exhibit broader peaks; their spectral intensity is often redistributed to multiple satellite features. This is a signature of the MQP regime. These results alone, however, do not indicate whether the satellites represent a weak coupling of, in principle, distinguishable QPs, and to what extent the MQP couplings can be captured perturbatively using single particle states forming only one Slater determinant.

To answer this, we compare the ASCI and MBPT calculations. We show the QP energies for all valence excitations and all methods in Fig. 2. As expected,^{30,31,33–36,52} GWΓ is closest to ASCI results (compared to HF and GW) and we find the best agreement for the holes of the HOMO states, which are in the SQP regime. For HOMO, the self-energy is merely responsible for shifting the poles of Green’s function, but the presence of the vertex corrections is important, as illustrated in the inset of Fig. 2. The one-shot GW approach performs only slightly better than HF; the mean absolute deviation (MAD) with respect to ASCI is 1.0 and 0.8 eV for HF and GW, respectively. Upon inclusion of vertex terms, however, the MAD decreases to 0.3 eV. The presence of Γ is responsible for the excitonic effects in W, and these are likely non-negligible in small systems, where electrons and holes have large spatial overlaps. Furthermore, the vertex correction cancels, at least partially, the spurious self-polarization in GW.^56,57 The latter is possibly the major driving force of the improvement because the electron–hole interactions in the screening tend to have little or (surprisingly) negative impact on the QP energy predictions.^30,31,58

Unlike previous studies, we go beyond the HOMO and also examine higher energy excitations, which exhibit MQP behavior. Due to the lack of correlation, the HF ionization potentials deviate significantly from the ASCI QP energies (MAD ∼ 2 eV), as illustrated in Fig. 2. Even here, the ASCI spectral function retains a dominant peak (in Fig. 1 at −23 for CH₄, −27 for NH₃, and −33 eV for H₂O); yet, a fraction of the spectral weight is redistributed to satellite features. This departure from the SQP regime progresses with the increasing excitation energy. For the highest energy valence states of H₂O and N₂ (at −33 and −39 eV, respectively), A(ω) shows multitude of sizable satellites located at >10 eV away from the main peak. Unsurprisingly, the MBPT quasiparticle energies (QPEs) deviate more strongly from ASCI in this regime; see the inset of Fig. 2. Nevertheless, for these higher energy states, the vertex correction again significantly improves the QPEs of the most prominent peaks, bringing the MAD below 0.5 eV. Considering the main QP signatures alone, the vertex corrections seem necessary for an accurate description throughout all studied energy scales, suggesting that MQP effects are important even for the simple orbital ionization energies.

We now turn to address the satellites in the MQP regime. Processes at this energy scale are characterized by the presence of neutral (e.g., optical) excitations interacting with the ionized hole. However, the strength of their mutual “hybridization” [i.e., the degree of entanglement in Eq. (1)] is not encoded in A(ω). The spectral function only provides information about the probability of the given (SQP or MQP) excitation, as given by the overlaps $ΨmN−1|ci|Ψ0N$ in Eq. (1). In the rest of this paper, we explore how the MBPT methods may offer some clarification by investigating the properties of Σ(ω), since the poles of A(ω) are related to Σ(ω) through Eq. (2).

A graphical solution to Eq. (2) is illustrated in the lower panels of Fig. 1 for each orbital explicitly. We realize that the distinction between SQP and MQP regimes made in the upper panels is evident in the lower panels as well. Indeed, for the HOMO and two subsequent orbitals (labeled I, II, and III), the self-energy crosses the y = ω line exactly once, while for the inner valence orbital (IV), there may be multiple crossings, giving rise to satellite features. The GW A(ω) in the MQP regime shows a common simple structure: besides the main peak, there are two other smaller peaks accompanying it (see the arrows in the inset of Fig. 1). This three-peak feature corresponds to the self-energy with only a single true pole (middle arrow in the lower central panel of Fig. 1), but having two nearby frequencies where Eq. (2) is approximately fulfilled (left and right arrows). In almost all of the molecules, the main maximum corresponds to one of these two “pseudo-poles” (rightmost arrow). This same three-peak structure has been discussed in the context of solids as an artifact of GW.^59–61 Note that due to the finite real-time propagation employed in our stochastic implementation,^{20,23–25,31} some of the GW “pseudo-poles” may correspond to actual poles. Only for N₂ does Eq. (2) have two solutions (i.e., two poles). Here, the resulting A(ω) features only a single prominent peak associated with the low-energy QP; see the supplementary material. In general, neither the main QP peaks nor the “secondary” solutions in GW appear at the correct energy and the overall spectrum is on average shifted to too high energies. To summarize, ASCI shows a complex satellite structure, whereas GW only presents a regular three-peak spectra for the bottom valence states.

Clearly, GW can provide multiple solutions, in principle, though they do not match the fully correlated results.^32,61,62 The QPs predicted by GW in closed shell systems are consistent with a surmised “plasmaron,” representing a resonantly bound hole and a collective neutral excitation. In solids, this was interpreted as an electron–plasmon state. However, this was eventually identified as an artifact of GW.⁶³ In practice, the GW approximation spuriously substitutes the multiple satellites with a single secondary QP that does not correspond to a physical excited state. Nevertheless, for weakly interacting MQPs, the absence of satellites can be remedied by reconstructing Green’s function [and A(ω)] via the cumulant expansion technique.^25,61 This method reproduces MQP structures, but in its common linear order formulation assumes the presence of a distinguishable neutral excitation, which corresponds to a pole in the (classical) polarizability, which describes the charge density fluctuations. This in turn implies that $|ΨmN−1$ in Eq. (1) has a direct product state of an ionized hole and excited state determinant; the total energy of such a state is merely the sum of the ionization potential and the neutral excitation energy. If true, A(ω) would exhibit a regular satellite structure, in which peaks appear at energies corresponding to the multiples of the “plasmon” energy. In other words, the nth satellite maximum represents the energy of a single hole plus n excited plasmons. The weak coupling regime is justified for localized holes in the presence of delocalized neutral excitations.⁶³ However, this separation is hard to conceive in small finite systems and cannot be readily justified for molecules. Furthermore, the fully correlated ASCI results do not exhibit a regular pattern of satellite peaks, suggesting that the MQP regime comprises at least some entangled $|ΨmN−1$ states. A recent extension of the cumulant approximation^64,65 includes non-linear contributions and it is not limited to the regular satellite structure (important, e.g., for molecular core spectra). Unlike core holes, the energy scale of inner valence holes is comparable to that of the valence neutral excitations. As a result, a high degree of “hybridization” between the two QPs is expected, requiring an explicit inclusion of MQP interactions through the vertex Γ.

The two particle interactions present in GWΓ enable mutual couplings between holes and neutral excitations. Hence, if the MQP regime stems from the dynamical correlation, the frequency dependence of the self-energy should capture it. Indeed, the GWΓ yields multiple peaks in A(ω), recovering the same kind of rich satellite structure present in the ASCI spectra. In general, the maxima in A(ω) are close to the excitations predicted from ASCI. From the current results, it appears that there is no common tendency to shift the spectrum (toward neither higher nor lower frequency). To understand the origin of the multi-peak structure, we illustrate GWΓ self-energy in Fig. 1 for CH₄, NH₃, and H₂O. By inspecting the results closely, we see that the distant satellites appear either (i) due to the small denominator in Eq. (1),⁶⁶ or (ii) because Eq. (2) is satisfied at multiple frequencies.

Regardless of the system, the vertex correction is responsible for new features in Σ, especially in the inner-most valence energy regions. For NH₃, the double peak in the ∼−30 eV region is accompanied by small satellites at lower energies stemming from the rapid variation of Σ. The oscillatory behavior increases with the energy of the excited holes. This is clearly seen for the bottom valence region in H₂O and N₂; their inner valence states are at lower energies than NH₃ and the self-energy, indeed, exhibits stronger oscillations leading to a plethora of poles (more than ten) with energies almost 30 eV below the main QP peak (see Fig. 1 and the supplementary material).⁶⁷ Furthermore, we observe that the vertex corrections introduce “pseudo-poles” in the self-energy of the highest energy hole states that are resonant with the QP peaks of the lower energy holes, cf. the “pseudo-pole” at −11 eV in the highest energy hole state (IV) of NH₃, resonant with the HOMO (I) QPE. This seems to underline the inter-orbital couplings arising from the vertex correction, which may be further related to a mixed orbital character that appears in the ASCI treatment; see the supplementary material. As of now, it is not possible to verify the nature of the GWΓ poles as true poles or pseudo-poles, since, as mentioned above, ASCI Green’s functions are artificially broadened. A reliable classification of satellites in this regard would provide valuable information on the strength of the multi-quasiparticle interactions generating them.

In this Communication, we theoretically investigated the electronic correlation that affects valence ionization energies and is directly linked to photoemission experiments. We provided virtually exact ab initio valence spectra for small molecular systems and explored whether and how many-body perturbation theory captures the various excitations. Namely, we compared the spectral functions computed with one-shot GW and single-iteration GWΓ to adaptive-sampling CI. For the first time, we provide such a comparison for the entire range of valence excitations, i.e., the outer and inner valence holes and shake-up satellites.

We show that the neglect of explicit two-particle interactions in GW leads to substantial errors. In contrast, GWΓ results are close to the fully correlated calculations and exhibit a rich structure of multi-quasiparticle excitations stemming from the coupling between holes and optically excited states. While GW yields spurious solutions, the presence of vertex terms removes the artifacts and correctly reproduces the peak structure in A(ω).

The high-energy regime is typically associated with the breakdown of the quasiparticle picture. However, our results clearly show that the shake-up satellite features arise due to dynamical correlations. In other words, they can be described perturbatively using single-particle states of only one Slater determinant. Our findings should encourage the further development of perturbative methods that explicitly account for mutual multi-quasiparticle excitations via vertex terms. Beyond small molecular systems, this will have a decisive effect on the path toward a first-principles understanding of excited states, photoactivated chemical reactions, and quantum materials.

See the supplementary material for a detailed description of the ASCI, G₀W₀, and G₀W₀Γ implementation, completed with a convergence discussion. We further provide tables with the QPEs for all main peaks of the valence orbitals and details of the extrapolation to the complete basis set limit.

This work was supported by the NSF CAREER award through Grant No. DMR-1945098 (V.V.). In part, this work was supported by the NSF Quantum Foundry through QAMASE-i program Award No. DMR-1906325 (V.V.). This work was partially supported by an Obra Social “La Caixa” graduate fellowship (No. 100010434), with code LCF/BQ/AA16/11580047 (C.M.-Z.). N.M.T. and S.C. are grateful for support from the NASA Ames Research Center and support from the AFRL Information Directorate under Grant No. F4HBKC4162G001. The calculations were performed as part of XSEDE computational Project Nos. TG-CHE180051 and TG-MCA93S030. Computational facilities purchased with funds from the National Science Foundation (No. CNS-1725797) and administered by the Center for Scientific Computing (CSC) were used. The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF DMR 1720256) at UC Santa Barbara. We thank Mark Babin, Blake Erickson, and Diptarka Hait for helpful discussions.