Many-body theory of positron binding to polyatomic molecules

Positron binding to molecules is key to extremely enhanced positron annihilation and positron-based molecular spectroscopy1. Although positron binding energies have been measured for about 90 polyatomic molecules1–6, an accurate ab initio theoretical description of positron–molecule binding has remained elusive. Of the molecules studied experimentally, ab initio calculations exist for only six; these calculations agree with experiments on polar molecules to at best 25 per cent accuracy and fail to predict binding in nonpolar molecules. The theoretical challenge stems from the need to accurately describe the strong many-body correlations including polarization of the electron cloud, screening of the electron–positron Coulomb interaction and the unique process of virtual-positronium formation (in which a molecular electron temporarily tunnels to the positron)1. Here we develop a many-body theory of positron–molecule interactions that achieves excellent agreement with experiment (to within 1 per cent in cases) and predicts binding in formamide and nucleobases. Our framework quantitatively captures the role of many-body correlations and shows their crucial effect on enhancing binding in polar molecules, enabling binding in nonpolar molecules, and increasing annihilation rates by 2 to 3 orders of magnitude. Our many-body approach can be extended to positron scattering and annihilation γ-ray spectra in molecules and condensed matter, to provide the fundamental insight and predictive capability required to improve materials science diagnostics7,8, develop antimatter-based technologies (including positron traps, beams and positron emission tomography)8–10, and understand positrons in the Galaxy11.

Pioneering technological developments have enabled the trapping, accumulation and delivery [8][9][10] of positrons for study of their fundamental interactions with atoms and molecules, 1,12 and the formation, exploitation and interrogation of positronium (Ps) 13,14 and antihydrogen. 15The ability of positrons to annihilate with atomic electrons forming characteristic γ rays makes them a unique probe over vast length scales, giving them important use in e.g., materials science for ultra-sensitive diagnostics of industrially important materials, 7,8 medical imaging (positron emission tomography), 16 and in astrophysics. 11roper interpretation of the materials science techniques, and the development of next-generation antimatter-based technologies rely on an accurate understanding of the fundamental interactions of positrons with atoms and molecules.Yet, many basic aspects are poorly understood theoretically.A striking example is the open fundamental problem of positron binding to molecules.Positrons can readily attach to molecules that bind them via vibrational Feshbach resonances, leading to spectacularly enhanced annihilation spectra that exhibit pronounced resonances downshifted from the vibrational mode energies by the positron-molecule binding energy. 1 Observation of such energy-resolved annihilation spectra have enabled measurement of positron binding energies (ranging from a few to a few hundred meV) for more than 90 molecules.The majority of these (∼60) are nonpolar or weakly polar species, such as alkanes, aromatics, partially halogenated hydrocarbons, alcohols, formates, and acetates.Ab initio calculations have been performed predominantly for strongly polar molecules, see e.g., Refs.[17-24] (we b c a Fig. 1.Main contributions to the positron-molecule self energy.a, the 'GW ' contribution, which involves the positron Green's function Gν and the (dynamic part of the) screened Coulomb interaction W .It describes the positroninduced polarization of the molecular electron cloud and corrections to it due to screening of the electron-positron Coulomb interaction by molecular electrons, and electron-hole attractions.b, the virtual-positronium contribution Σ Γ , which includes the summed infinite ladder series ('Γ-block') of screened electron-positron interactions.c, the positron-hole ladder series (the 'Λ block') contribution Σ Λ .Lines directed to the right (left) represent particles (holes) propagating on the N -electron ground-state molecule: red lines labelled ε represent the external positron state; other red (blue) lines represent positron (excited electron or hole) intermediate states that are summed over; single (double) wavy lines represent bare (screened) Coulomb interactions.See 'Methods' and Extended Data Fig. 1 for details of their calculation via Bethe-Salpeter equations.
The theoretical difficulty lies in the need to identify and accurately describe the strong many-body correlations that characterise the positron-molecule system.Here, we develop the ab initio many-body theory of positron interactions with polyatomic molecules, which enables the natural, explicit, and systematic account of important correlations including polarization of the molecular electron cloud, screening of the electron-positron interaction, and the unique process of virtual-positronium formation (where a molecular electron temporarily tunnels to the positron).We focus its application, via a state-of-the-art computational implementation, to calculations of positron binding and annihilation for the six molecules for which both previous theory and measurements exist, and additionally formamide, CSe 2 , benzene and A, C, T, G and U nucleobases.The positron binding energy ε and bound-state wavefunction ψ ε are found by solving the Dyson equation 29 where H (0) is the Hamiltonian of the positron in the Hartree-Fock (HF) field of the ground-state molecule, and Σ E is a non-local, energy-dependent correlation potential (irreducible self energy of the positron).It acts as an integral operator Σ E ψ(r) ≡ Σ E (r, r )ψ(r )dr and encapsulates the full complexity of the many-body problem.We calculate Σ via its expansion in residual electron-electron and electron-positron interactions, see Fig. 1.Diagram (a), the 'GW ' self energy Σ GW , describes the positron-induced polarization of the molecular electron cloud, and corrections to it due to screening of the electron-positron Coulomb interaction by the molecular electrons, and electron-hole attractions (the Bethe-Salpeter Equation approximation, 'GW @BSE').Diagram (b), denoted by Σ Γ , which involves the summed infinite ladder series of (screened) electron-positron interactions (the 'Γ-block', see Extended Data Fig. 1), represents virtual-positronium formation. 30,31 ts importance is unique to the positron problem because successive diagrams in this series contribute to the positron-molecule self energy with the same sign, whereas for all-electron systems the series is sign alternating and gives a small overall contribution.We also consider the smaller positron-hole ladder series contribution, Σ Λ , Fig. 1 c.'Methods' details the state-of-the-art construction of Σ and solution of the Dyson equation.
Results and discussion: positron binding energies and lifetimes.
Table 1 shows our calculated binding energies at successively more sophisticated approximations to the correlation potential: HF; Σ (2) (bare polarisation); Σ GW (polarisation including electron screening and screened electron-hole interactions (Fig. 1 a); Σ GW +Γ (Fig. 1 a+b); Σ GW +Γ+Λ (Fig. 1 a+b+c): for this the first (second) number is the result using bare (dressed) Coulomb interactions in the ladders, whilst the third (our most sophisticated, in bold) is that using dressed interactions and energies.Also see Fig. 2 for a graphical comparison of theory and experiment, and Extended Data Table 2 for more details.  1 for anisotropic polarizabilities).Binding energy calculations are presented in the Hartree-Fock, Σ (2) (bare-polarisation) and GW @BSE (bare-polarisation plus screening and electron-hole corrections) approximations, and additionally including virtual-positronium formation Σ GW +Γ , and the positron-hole ladder contribution Σ GW +Γ+Λ : † for the latter, the first (second) number is that using bare (dressed) Coulomb interactions in the Γ and Λ blocks, and the third (our most sophisticated calculation, in bold) additionally uses GW energies in the diagrams.Their difference gives a measure of the theoretical uncertainty; ‡ Experiment from; 1, 4, 5 * except that for formamide, which is preliminary and unpublished. 37ther calculations: ECG (explicitly correlated Gaussian), CI (configuration interaction), and any-particle-molecular orbital (APMO, 'REN-PP3').Also see Fig. 2 for a graphical comparison of experiment with present and previous theory., the electron HOMO wavefunction density isosurface (blue lobe is negative region at 40% of maximum, and brown is positive region at 10% of the maximum).Also shown is the positron wavefunction calculated along the molecular axis in the static HF approximation (black curve) and at the Σ GW +Γ+Λ level of many-body theory (red curve).b, acetonitrile; c, propionitrile; d, acetone; e, propanal; f, acetaldehyde; and g, formamide; show the positron wavefunction density isosurface at 80% of the maximum.For the nonpolars: h, CS 2 ; i, CSe 2 ; and j, benzene, the isosurfaces are at 90% of the maximum.
Benchmarking and general trends.-Webenchmark our approach against a highly-accurate explicitly correlated gaussian (ECG) calculation (ε b = 1043 meV) 21 for the strongly polar molecule LiH.
The results demonstrate the general trends seen in all the molecules considered.The HF binding en-ergy (ε b = 130 meV) is severely deficient.Including the bare polarization attraction Σ (2) significantly increases the binding energy (to ε b = 434 meV).The addition of short-range screening corrections reduces the polarizability and binding energy (to ε b = 336 meV, see Extended Data Table 2), but this is compensated by the inclusion of the electron-hole attractions (Σ GW : ε b = 518 meV).This is still, however, less than half of the ECG result.The previous CI calculation 18 is similarly deficient.Strikingly, however, including the virtual-positronium formation correlation potential (Σ GW +Γ ) strongly enhances the binding, more than doubling it (to ε b = 1291 meV).Including the positron-hole ladder (Σ GW +Γ+Λ ) slightly reduces binding (to ε b = 1106 meV); using screened interactions in the ladders reduces it slightly further (ε b = 1038 meV); and additionally using the dressed energies in the diagram construction gives ε b = 1060 meV, agreeing with the ECG result to ∼1%.As for all the polar molecules, the maximum of the positron wavefunction density (Fig. 2) in LiH is highly localized at the negative end of the molecule, but overall the wavefunction is quite diffuse, asymptotically taking the form ψ ∼ e −κr where κ = √ 2ε b .We also calculate the positron Dyson wavefunction renormalization constants a (see 'Methods Eqn.7 and Extended Data Table 2).These represent the contribution of the 'positron plus molecule in the ground state' component to the positron-molecule bound state.Their closeness to unity suggests the picture of a positron bound to the neutral molecule (rather than a Ps atom orbiting a molecular cation). 38mparison with experiment and previous theory .-Thebest prior agreement between theory and experiment for any molecule was for acetonitrile ( 25%).Considering the polar molecules first (Table 1 and Fig. 2), we immediately see that our full many-body theory (Σ GW +Γ+Λ ) is much superior, giving near exact agreement ( 1% level) with experiment for propionitrile, propanal, acetaldehyde and formamide, and within 10% for acetonitrile and acetone (including the experimental error).(Overall we find excellent convergence in our calculation: see 'Methods' and Extended Data Fig. 2).For all the polar molecules, the HF and bare (Σ (2) ) and dressed (GW ) polarization results significantly underestimate binding.The effect of virtual-positronium is crucial: it enhances the binding energy by a factor of ∼2 and is essential to bring theory into agreement with experiment.We note that the previous configuration interaction and any-particle-molecular orbital (renormalized PP3 "REN-PP3", which employs a diagonal approximation and does not explicitly account for virtual-positronium formation) calculations are severely deficient.The ECG approach is not easily scalable to these molecules. 21or the nonpolars, binding is exclusively enabled by correlations.For CS 2 a considerable binding energy of 75 meV has been measured, whilst the CI calculation failed to predict binding. 35We find that polarization correlations (GW ) alone are insufficient to support binding.Spectacularly, however, including the virtual-positronium contribution results in significant binding: our Σ GW +Γ+Λ result of ε b = 46-68 meV is close to experiment: the range is larger than in the calculations for the polars as the delocalization of the positron wavefunction (Fig. 2 h-j) makes accurately describing virtual-positronium more demanding.For CSe 2 and benzene, in contrast to the molecules already considered, we have not performed any optimization of the bases (due to the delocalised positron wavefunction, these molecules require computational resources beyond our disposal) and our calculated ε b should be considered as lower bounds.Nevertheless, the results further elucidate the essential role of virtual-positronium formation in enabling (significant) binding, and the positron wavefunctions also provide fundamental insight that may prove instructive to refine ab initio and model calculations (see below).Prediction for formamide.-Forthis, the archetypal molecule for the investigation of protein and peptide chemistry, we are unaware of any prior calculation.We predict binding (ε b ∼ 189 meV).Preliminary experiments see evidence of ε b ∼200 meV, although a final value has yet to be determined. 37lecular orbital contributions to binding: anisotropy and strength of correlations.-Atthe static HF level, we find ε b to be (monotonically and non-linearly) related to the permanent dipole moment (expected from the dipole-potential model 40 ).Ultimately the correlation potential is anisotropic Fig. 3. Molecular orbital contributions to binding, and scaling formula fo to the dimensionless strength of the virtual-positronium formation correlation pot polarization S (2+ ) (squares) and the ratio g ⌘ S (2+ ) /S (2) (crosses, where S (2) = where I is the MO ionization energy and EPs = 6.8 eV is the ground state en orbitals with I < 15 eV, and inset MO plots show HOMO ( type) and next HO and blue: positive and negative regions of electronic MO; red wireframe: positron lines show the fits S ⇡ ae bI + cI d with a = 0.67 (2.57); b = 0.121 (0.092); c = 2. b, the positron wavefunction density in acetaldehyde (in the plane containing th plane), which protrudes along the ⇡ bond.c, comparison of binding energies f ⌃ = ⌃ GW + +⇤ and accounting for (the computationally demanding to calculat (circles) and 1.5 (squares) (see text).Table 2, " b in nucleobases, calculated with (see Extended Data Table 1 for calculated anisotropic polarizabilities)  the ratio g ⌘ S (2+ ) /S (2) , where where S (2) = S (2+ ) S ( ) .Both S 149 Ps-formation threshold to higher ionization energies: it is more di cu (see Extended Data Table 1 for calculated anisotropic polarizabilities), and depends non-linearly on the polarizabilities and ionization energies of the individual MOs.Moreover, the binding energy depends nonlinearly on the correlation potential (e.g., see Extended Data Fig. 3).The ordering of ε b with respect to dipole moment persists to the full Σ (2+Γ+Λ) calculation, with the exception of acetaldehyde and propanal, and we note that for acetone correlations considerably enhance ε b .It is instructive to consider the dimensionless quantity S = − ν>0 ε −1 ν ν|Σ|ν 41 (where the sum is over excited HF positron basis states of energy ε ν , see 'Methods'), which gives an effective measure of the strength of the correlation potential Σ.The magnitudes of the strength of Σ (2) , S (2) ranges from 4-15 (see Extended Data Table 2), and follows the ordering of the isotropic polarizability, with the exception of acetone and propanal (acetone has a larger polarizability and smaller ionization energy than propanal), and benzene and CSe 2 (owing to benzene's π bonds, see below).This suggests that (the short range contributions to) Σ (2) cannot be parametrized solely by the polarizability.Similarly, the magnitudes of S (Γ) (ranging from 2-5) do not strictly follow the ordering of the ionization energies.To illuminate this, note that at the bare-polarization, Σ (2) , and polarization plus virtual-positronium formation approximations, Σ (2+Γ) = Σ (2) + Σ (Γ) , we can delineate the contribution of individual MOs to positron binding.Fig. 3 a shows the partial S (Γ) and S (2+Γ) for individual occupied MOs against their respective ionization energies, and the ratio g ≡ S (2+Γ) /S (2) , where S (2) = S (2+Γ) − S (Γ) .Both S (Γ) and S (2+Γ) decrease from the Ps-formation threshold to higher ionization energies: it is more difficult to perturb more tightly bound electrons.However, the decrease is not monotonic: we see that despite having larger ionization energies, π-type electronic MOs below the HOMO can contribute significantly more than a σ-type HOMO to S (Γ) and S (2+Γ) , e.g., in acetone, propanal, and acetaldehyde, the strength of the π-type (H-1)OMO is larger than the σ-type HOMO, and in propanal, the (H-3)OMO of π type contributes more strongly than the (H-2)OMO etc.It was previously speculated in Ref. [3] that π bonds were important due to the ability of the positron to more easily access electron density that is delocalized from (repulsive) nuclei.This is borne out by our calculations, and we see that in       Table : Positron lifetimes with respect to annihilation.τ (0) : lifetime calculated in the Hartree-Fock independent particle approximation excluding the vertex enhancement factors and using a positron wavefunction normalised to unity; τ GW and τ : lifetime calculated using the Dyson positron wavefunction at the Σ GW and Σ GW +Γ+Λ levels including vertex enhancement factors and renormalization constants (Eqns.6 and 7 of 'Methods').
strength.For acetonitrile this results in a larger strength parameter than formamide.
Predicting positron binding in larger molecules: nucleobases as an example.-Theratio g ≡ S (2+Γ) /S (2) is weakly dependent on the ionization energy, with a value of ∼ 1.4-1.5 for the HOMOs (I ∼ 10 eV).We propose that binding energies of large molecules (e.g., 15-100 atoms, for which a converged calculation of the virtual-positronium diagram Fig. 2 c is too computationally demanding) can be calculated by approximating Σ ≈ gΣ (2) +Σ Λ .As well as accounting for virtual-positronium formation, this model potential reflects the anisotropy of the true interactions.For the molecules considered in Table 1, this works well (see Fig. 3 c and also Extended Data Fig. 3).Using this approximation, we calculate the positron binding energy in the five nucleobases (Table 2).Our results are larger than the previous APMO calculations, mirroring the results for the molecules in Table 1.We predict binding in adenine.
Positron lifetimes.Correlations and the relative MO symmetries also dictate the lifetime τ [ns] ≈ 0.02 δ −1 ep [a.u.] of the bound-state positron against annihilation, where δ ep is the electron-positron contact density [see 'Methods' Eqn. ( 6) and Extended Data Table 3 for values of δ ep ].We also calculate the annihilation lifetime of the bound positron (see Figs. 4, 'Methods' and Extended Data Fig. 4), finding that the correlations reduce it by two orders of magnitude to τ ∼ 1 ns, and find that the contribution of electron MOs to annihilation depends strongly on their symmetry relative to that of the positron MO, with the HOMO not necessarily dominating, e.g., in acetonitrile and formamide.

Conclusions and future perspectives.
Many-body theory has uncovered the mechanisms of positron binding to molecules.Binding is governed by the interplay of the static (dipole) interaction with the strong (and non-local) correlation potential, to which numerous MOs contribute significantly.The interactions are inherently anisotropic, and the binding energy depends non-linearly on the correlation potential.In particular, we uncovered the key role of virtual-positronium formation in significantly enhancing and enabling binding, and quantified the importance of π bonds.Overall, the ab initio approach, which takes proper account of the distinct correlations and the ansiotropy of the interactions gives binding energies in the best agreement with experiment to date (highlighting the severe deficiency of previous quantum chemistry calculations).We also predicted binding in formamide and nucleobases.
The present work directly supports resonant annihilation experiments and related model theory that requires positron binding energies as parameters. 1 Moreover, the many-body theory approach can be naturally extended to describe positron scattering and (non-resonant) annihilation rates and γ spectra in molecules and condensed matter.Calculations of γ spectra can provide insight on energy deposition and molecular fragmentation, [42][43][44][45] clusters and cluster surfaces, 46 and the determination of electronmomentum densities from annihilation γ spectra measurements. 47As well as novel positron-based molecular spectroscopy, more generally, such predictive capability could provide fundamental insight required to develop positron traps, accumulators and high-energy-resolution beams 9,46 (for e.g., shorter temporal pulses in lifetime spectroscopy, colder positrons for antihydrogen formation, and higher densities for Ps BEC production 48 ), to develop next-generation positron emission tomography (including understanding positron emitting tracers and positron-induced DNA damage, and new spectroscopic PET 48 ), and proper interpretation of positron-based materials science diagnostics. 7,8 urthermore, the formation and interrogation of antihydrogen 15 gives promise for more complicated anti-atoms/molecules: study of their structure, interactions and chemistry will require theoretical understanding.Finally, the many-body formalism provides a natural foundation to develop ab initio descriptions of positron-induced (excited-state) molecular processes, such as positron-induced interatomic Coulomb decay or positron capture, 49 charge migration and luminescence, 50 and (whilst extremely challenging) the inclusion of phonons and photons, which may enable ab initio description of vibrational Feshbach resonant annihilation spectra via coupling of the nuclear and electronic degrees of freedom, or modelling of positron pump-probe experiments. 7

Methods
Solving the Dyson equation for the positron binding energy and wavefunctions using a Gaussian basis.
We calculate positron-molecule binding energies ε and quasiparticle wavefunction ψε by solving the Dyson equation, main text Eqn.(1).We take the zeroth-order Hamiltonian H (0) to be that of the positron in the Hartree-Fock field of the frozen-target N -electron ground-state molecule.The self-energy diagrams thus begin at second order in the Coulomb interaction.Rather than computing the self energy Σ(r, r ) in the coordinate basis, it is more convenient to work with its matrix elements in the Hartree-Fock basis.Specifically, we expand the electron (-) and positron (+) Hartree-Fock MOs ϕ ± a (r) in distinct Gaussian basis sets as ϕ ± a (r) = , where A labels the N ± c basis centres, k labels the N ± A different Gaussians on centre A, each taken to be of Cartesian type with angular momentum l x + l y + l z , viz., , where NA k is a normalisation constant, and C are the expansion coefficients to be determined (see below).
For both electrons and positrons, we use the diffuse-function-augmented correlation-consistent polarised augcc-pVXZ (X=T or Q) Dunning basis sets centred on all atomic nuclei of the molecule, which enables accurate determination of the electronic structure including cusps 51 and expulsion of the positron density from the nuclei.To capture the long-range correlation effects, for the positron we also additionally include at least one large eventempered set at the molecular centre or region of maximum positron density of the form N s(N −1)p(N −2)d(N − 3)f (N − 4)g with N ∼ 10 − 15 (where it should be understood that the full degenerate set of non-zero angular momentum functions is used) and exponents ζA k = ζA 1 β k−1 , k = 1, . . ., N , for each angular momentum type, where ζA 1 > 0 and β > 1 are parameters.The value of ζA 1 is important as the bound positron wavefunction behaves asymptotically as ψ ∝ e −κr , where κ = √ 2ε b .Thus, to ensure that the expansion describes the wavefunction well at r ∼ 1/κ, i.e., that the broadest Gaussian covers the extent of the positron-wavefunction, one must have ζA 1 κ 2 = 2ε b .In practice we performed binding energy calculations for a range of ζA 1 and β for each molecule, finding that there are broad ranges of stability.The optimal ζA 1 was typically found to be in the range of 10 −4 − 10 −3 for s-and p-type Gaussians and 10 −3 − 10 −2 for d-and f -type Gaussians, whilst g-type Gaussian exponents usually had ζA 1 = 10 −1 (atomic units are assumed throughout unless otherwise specified).The optimal β ranges from 2.2 to 3.0 depending on the number of functions N in a given shell.Finally, to improve the description of the virtual-Ps formation process, which occurs several atomic units away from the molecule and requires large angular momenta, additional (aug-cc-pVXZ, X=T,Q) basis sets are strategically placed at 'ghost' centres close to the regions of maximum positron density.To check convergence with respect to the number and location of these ghost centres, for each molecule we performed calculations including TZ or QZ bases on a successively increasing number of ghosts centres in different arrangements until the increase in binding energy fell below a few percent.We found that including ghosts can increase binding energies by ∼ 10% in the polar molecules, and easily by ∼ 30% for the nonpolar ones, e.g., for CS2 we obtained ε b = 39 meV at GW @BSE+Γ+Λ level with no ghosts, rising to ε b = 68 meV with 16 additional ghosts.The use of higher angular momenta and more ghosts would further increase the binding energies.We also investigated the difference of using aug-cc-pVXZ for X=T,Q in the atomic centred and ghost bases, higher angular momenta in the even tempered basis.Some improvement was noted moving from X=T to Q, and from including g states in addition to f , to a level of a 5-10% in polar molecules, and 10-30% in nonpolars.Overall, good convergence with respect to both the electron and positron bases was observed (see e.g., Extended Data Fig. 2).
The coefficients C in the expansion of the positron wavefunction in Gaussians are found by solving the Roothaan equations F ± C ± = S ± C ± ε ± , where F ± is the Fock matrix and S is the overlap matrix.The one-body and two-body Coulomb integrals of the Fock matrix are calculated using the McMurchie Davidson algorithm. 52We eliminate linearly-dependent states by excluding eigenvalues < 10 −5 of the overlap matrices (typically 5% of the states).In practice, to minimise the basis dimensions we transform all quantities to a spherical harmonic Gaussian basis (for a given angular momentum, the number of Cartesian Gaussians is greater than or equal to the number of spherical Harmonic Gaussians). 53Solution of the Roothaan equations yield bases of electron and positron Hartree-Fock MOs {ϕ ± α (r)} (which include ground and other negative energy states, and discretized continuum states) with which the self-energy diagrams can be constructed (see below for details).
Expanding the positron Dyson wavefunction (see Eqn. 1 of main text) in the positron HF MO basis as ψε(r) = ν D ε ν ϕ + ν (r) transforms the Dyson equation to the linear matrix equation HD = εD, where ν1|H|ν2 = εν 1 δν 1 ν 2 + ν1|Σε|ν2 .Note that we calculate the full self-energy matrix including off-diagonal terms.Such a nonperturbative approach is essential for nonpolar molecules, where binding is enabled exclusively by correlations.In practice, to obtain the self-consistent solution to the Dyson equation, we calculate the self energy at a number of distinct energies Ei spanning the true binding energy ε b , with the latter determined from the intersection of the ε b (Ei) data with the line ε b (E) = E.

Constructing the positron-molecule self energy via solution of the BSE equations.
As discussed in the main text (Fig. 1), we consider three contributions to the irreducible self energy of the positron in the field of the molecule: Σ GW (which describes polarization, screening and electron-hole interactions); Σ Γ (which describes the non-perturbative process of virtual-positronium formation); and Σ Λ (which includes the infinite ladder series of positron-hole interactions).In practice, we construct the individual contributions by first solving the respective Bethe-Salpeter equations (see Extended Data Fig. 1) for the electron-hole polarization propagator Π, the two-particle positron-electron propagator G ep II and the positron-hole two-'particle' propagator G ph II . 29Their general form is L(ω) = L (0) (ω) + L (0) (ω)KL(ω) where the L (0) are non-interacting two-body propagators and K are the interaction kernels 29,54,55 [e.g., see Extended Data Fig. 1 e for the BSE for the electron-hole polarization propagator Π].In the excitation space of pair product HF orbitals , where the pair transition amplitudes ξ are the solutions of the pseudo-Hermitian linearresponse generalised eigenvalue equations [55][56][57] Hξ = CξΩ, ξ † Cξ = C, where for excitation energies Ω α + and Ω α − , which are labelled by α = 1, . . ., dim(A).Here the A and B matrices depend on the particular two-particle propagator L under consideration and the approximation used for it, (see Extended Table 4 for the explicit matrix elements): note that B = 0 for the two-particle propagators involving the positron, since the vacuum state for the diagrammatic expansion is that of the N -electron molecule, and thus there are no positron holes and only time-forward positron propagators.To determine the amplitudes, we employ the parallel diagonalisation algorithm of Shao, 58 which exploits a similarity transform that gives the eigenvalues of C −1 H as the square roots of the eigenvalues of (A + B)(A − B) (thus requiring matrices of dimension of the A block, i.e., half of the full BSE matrix dimension) to obtain X = 1 2 (L2U + L1V ) , via the Cholesky decompositions A + B = L1L T 1 and A − B = L2L T 2 , and the singular value decomposition L2L T 1 = U ΩV T .The positron-molecule self-energy matrix elements can then be written as where ν1, ν2 and ν3 denote positron indices and µ and n denote electron excited states and holes respectively, and Σ (2) -which results from the Π (0) contribution to Σ GW and is present in both G ep II and G ph II -is subtracted to prevent double counting, and The total self energy is thus calculated as Σ = Σ GW + Σ Γ + Σ Λ .][63] We implement the above in the massively-parallelised EXCITON+ code developed by us, heavily adapting the EXCITON code [64][65][66] to include positrons and the many-body theory capability (calculation of the self energy and solution of the Dyson equation).We employ density fitting [66][67][68][69][70][71] (of the electronic density) to calculate the Coulomb integrals in the matrix elements of A and B, in w Π , w Γ and w Λ , and positron-electron contact density, via a parallel implementation that assigns matrix elements involving auxiliary basis functions on distinct atomic centres to distinct processors, similar to that used in the MolGW program. 72The employment of density fitting reduces four-centre Coulomb integrals to products of three-centre Coulomb integrals and matrix elements of the Coulomb operator between atomic orbital basis functions.Thus, the memory scaling is ∼ N 2 − M−, where N− is the total number of electron basis functions, and M− 3N− is the number of electron auxiliary basis functions.The most computationally demanding part of our approach is in the calculation of the virtual-Ps self-energy contribution Σ Γ .For this, dim A = dim X Γ = Nν × Nµ, the product of total number of positron MOs and excited electron MOs.For the calculations considered here, Nν ranged from 400-500 and Nµ from 300-400, resulting in dim X Γ = 120, 000-200,000, and thus diagonalising the matrix of (dim X Γ ) 2 elements demanded between ∼100 GB and 1.5 TB of RAM.The calculations were performed on two AMD EPYC 128 CPU @ 2 GHz, 768GB RAM nodes of the United Kingdom Tier-2 supercomputer 'Kelvin-2' at Queen's University Belfast.In contrast, the GW calculations involve dim A = dim X Π ≤ Nν × Nn, i.e., a maximum equal to the product of the number of occupied and excited electron MOs: in practice not all occupied orbitals need to be included because the tightly bound LOMOs are less susceptible to perturbation by the positron and have negligible contribution to the self energy.Thus, since Nn Nµ < Nν , ab initio GW @RPA/TDHF/BSE calculations are considerably less computationally expensive, and can be performed for molecules or clusters with ∼ 100 atoms, providing at least lower bounds on the positron binding energies.Moreover, as discussed in the main text (see Fig. 3 c of main text and Extended Data Fig. 3) and demonstrated for nucelobases (Table 2 of main text), the virtual-Ps formation contribution can be approximated by scaling the Σ (2) self energy by the strength parameter ratio g ≡ S (2+Γ) /S (2) , viz.Σ ≈ gΣ (2) +Σ Λ , thus enabling computationally relatively inexpensive binding-energy calculations that account for virtual-Ps formation for molecules of ∼ 100 atoms.Ab initio calculations for larger molecules including the virtual-positronium self energy will be feasible with additional computational resources, as would calculations using different truncated product spaces of excited electron and positron MOs and extrapolating to the basis set limit.

Improving the accuracy of calculations
As mentioned in the previous section, the computationally intensive calculations presented here were performed using relatively modest computational resources.Access to national supercomputing facilities would enable more complete basis sets and further exploration of the effect of ghost basis centres.Numerical accuracy can also be systematically improved in a number of ways.Exploiting the molecular point group symmetry via symmetry adapted bases and employing integral screening techniques would improve the efficiency of the calculations, enabling more complete basis sets to be used.This would ultimately improve the description of the correlations (particularly in generating higher angular momenta for improved description of the virtual-positronium formation process).The calculation of the positron-molecule self energy can be improved by implementing a self-consistent diagram approach in which the positron-molecule self energy is built from GW calculated electron and positron Dyson orbitals rather than HF ones, 29,73 and/or by coupling the three self-energy channels Σ GW , Σ Γ and Σ Λ by approximating the three-particle propagators via the Fadeev, 74 parquet 61 or ADC(3) methods 75 (expected to be computationally feasible for small molecules using national supercomputing facilities).Moreover, the diagrammatic series should be amenable to a diagrammatic Monte Carlo 76,77 prescription, a stochastic simulation method that enables the effective summation of many more (classes of) diagrams than considered here.

Positron annihilation rate in the bound state
Here the sum is over all occupied electron MOs with wavefunctions ϕn, and γn are MO dependent enhancement factors that account for the short-range electron-positron attraction. 78,79 ecent many-body calculations for atoms by one of us determined them to follow a physically motivated scaling with the ionization energy 78, 79 γn = 1 + 1.31/|εn| + (0.834/|εn|) 2.15 , which we assume to hold here.The positron Dyson wavefunction is a quasiparticle wavefunction that is the overlap of the wavefunction of the N -electron ground state molecule with the fully-correlated wavefunction of the positron plus N -electron molecule system. 29It is normalised as which estimates the contribution of the 'positron plus molecule in the ground state' component to the positronmolecule bound state wavefunction, i.e., the degree to which the positron-molecule bound state is a single-particle state, with smaller values of a signifying a more strongly-correlated state.Setting K = v, the direct part of the Coulomb interaction only, gives the 'random phase approximation' (GW @RPA).Setting K = v−v exch , i.e., including exchange which gives rise to interactions within the bubbles, yields the 'time-dependent Hartree-Fock' approximation (GW @TDHF).Using screened Coulomb interactions in the exchange term is 'Bethe-Salpeter' approximation (GW @BSE).g, the summed infinite ladder series of screened electron-positron interactions ('Γ-block').The Λ-block (diagram d) is the ladder series of positron-hole interactions, and it satisfies a linear integral equation of the same form as that shown in g.
Extended Data Table 1: Calculated polarizabilities (in Å3 ) and ionization energies (in eV).[80].Molecules are oriented such that the main axis of symmetry, or the main bond (C-O, C-N), above which the positron density is localised, is along z.The isotropic value is given by a sum of xx, yy, and zz terms multiplied by 2/3.Note that the zz components have larger differences between molecules than isotropic polarizabilities, e.g., for propionitrile, acetone and propanal the isotropic polarisabilities are within 1% of each other, whereas the zz components differ by ∼ 15%.Ionisation energies calculated at the GW @RPA level were performed using the diagonal approximation for the electron-molecule self energy Σ (−) , i.e., εµ = εµ + Z µ|Σ ep , using HF positron wavefunction; δ GW ep , using the Dyson wavefunction calculated with the GW @BSE self energy; δ GW +Γ ep , using the Dyson wavefunction calculated with the GW @BSE plus virtual-positronium self energy; δ GW +Γ+Λ ep , using the Dyson wavefunction calculated with the GW @BSE plus virtual-positronium plus positron-hole self energy with unscreened Coulomb interactions in the Γ and Λ ladders; δ GW + Γ+ Λ ep , using screened Coulomb interactions in the ladders and † additionally using GW energies in the diagrams.Numbers in brackets indicate powers of 10.Hyphens denote approximations in which the positron does not bind.Elements of the A and B blocks of the linear-response Hamiltonian matrices that result from the BSE equations ('Methods' Eqn. 1) for the electron-hole propagator, the positron-electron propagator, and the positron-hole propagator.Chemists' notation for Coulomb matrix elements in the MO basis is used (νµ|ν µ ) = drdr ϕ * ν (r)ϕµ(r)v(r, r )ϕ * µ (r )ϕν (r ), where (µ and µ ), (n and m) and (ν and ν ) denote electron particles, electron holes and positron particles respectively.Factors of two arise from summation over spin, and tildes on energy eigenvalues for BSE denote that these are calculated at the level of GW @RPA.For the virtual-Ps and positron-hole matrices, B = 0 because there are no positron 'holes' in the N -electron ground-state molecule vacuum-state, and thus only time-forward diagrams are present in the positron single-particle propagator and two-particle propagators (i.e., here the 'Tamm-Dancoff approximation' is exact for positrons).Matrix elements of the dressed Coulomb interaction W = v + W d (Extended Data   Extended Data Figure 3. Non-linearity of the binding energy and strength of correlation potential.Binding energy calculated approximating the positron self energy Σ as Σ ≈ gΣ (2) + Σ Λ as a function of the scaling parameter g ≡ S 2+Γ /S (2) (see main text for more details) (circles).Experiment (squares) is from [4, 5]; for formamide preliminary measurements find a binding energy of ε b ∼ 200 meV, but a final result is yet to be determined.See also Fig. 3 c 3 and Fig. 4 of main text.b, Contact density calculated at the HF (circles) and various levels of many-body theory (diamonds: GW @BSE; squares: GW @BSE+Γ+Λ) against the square root of the binding energy.

Fig. 2 .
Fig.2.Positron-molecule binding energies and bound-state Dyson wavefunction densities.Graph shows the comparison of the present many-body calculations (red circles) with experiment.Also shown are the CI and APMO calculations (blue square and crosses, respectively).Positron wavefunction densities: a, LiH, with Li atom at origin and H at ∼3 a.u.along the molecular axis, showing the positron wavefunction density isosurface at 70% of the maximum (red lobe), the electron HOMO wavefunction density isosurface (blue lobe is negative region at 40% of maximum, and brown is positive region at 10% of the maximum).Also shown is the positron wavefunction calculated along the molecular axis in the static HF approximation (black curve) and at the Σ GW +Γ+Λ level of many-body theory (red curve).b, acetonitrile; c, propionitrile; d, acetone; e, propanal; f, acetaldehyde; and g, formamide; show the positron wavefunction density isosurface at 80% of the maximum.For the nonpolars: h, CS 2 ; i, CSe 2 ; and j, benzene, the isosurfaces are at 90% of the maximum.

133
polarizabilities and ionization energies of the individual MOs.Moreover 134 linearly on the correlation potential (e.g., see Extended Data Fig.3).

135
to dipole moment persists to the full ⌃ (2+ +⇤) calculation, with th 136 propanal, and we note that for acetone correlations considerably enhan 137 the dimensionless quantity S = P ⌫>0 " 1 ⌫ h⌫|⌃|⌫i 41 (where the sum 138 states of energy " ⌫ , see 'Methods'), which gives an e↵ective measure 139 potential ⌃.The magnitudes of the strength of ⌃ (2) , S (2) ranges from 140 2), and follows the ordering of the isotropic polarizability, with the e 141 (acetone has a larger polarizability and smaller ionization energy than 142 (owing to benzene's ⇡ bonds, see below).This suggests that (the s 143 cannot be parametrized solely by the polarizability.Similarly, the m 144 2-5) do not strictly follow the ordering of the ionization energies.To 145 bare-polarization, ⌃ (2) , and polarization plus virtual-positronium form 146 ⌃ (2) + ⌃ ( ) , we can delineate the contribution of individual MOs to 147 the partial S ( ) and S (2+ ) for individual occupied MOs against their r 148
Fig. 3 b considerable positron density protrudes into the region of the π bond.Acetonitrile and propionitrile have a doubly-degenerate π HOMO of large of (screened) electron-positron interac It is unique to the positron problem because success 80 molecule self energy with the same sign, whereas f 81 We also consider the smaller positron-hole ladder se 82 details of the state-of-the-art construction of the se

90
Benchmarking and the general trends.-Wefir 91 explicitly correlated gaussian (ECG) calculation ( 92 LiH.The results demonstrate the general trend s 93 calculated binding energy (" b = 130 meV) severely 94 correlations.Including the bare polarisation attra 95 energy (to " b = 434 meV).The addition of short-r 96 polarisability and binding energy (to " b = 336 me 97 electron-hole attractions (GW @TDHF: " b = 542 m 98 half of the ECG result.The previous CI calculation 1 99 the virtual-positronium formation correlation poten 100 more than doubling it (to " b = 1290 meV).The ad 101 causes a relatively small reduction (to " b = 1095), a 102 ladder series of (screened) electron-positron interactions, r 79 It is unique to the positron problem because successive diag 80 molecule self energy with the same sign, whereas for all-e 81 We also consider the smaller positron-hole ladder series con 82 details of the state-of-the-art construction of the self-energ

Fig. 2 s 90 1 92
Benchmarking and the general trends.-Wefirst benc 91 explicitly correlated gaussian (ECG) calculation (" b = LiH.The results demonstrate the general trend seen in 93 calculated binding energy (" b = 130 meV) severely underes 94 correlations.Including the bare polarisation attraction (G 95 energy (to " b = 434 meV).The addition of short-range sc 96 polarisability and binding energy (to " b = 336 meV), but 97 electron-hole attractions (GW @TDHF: " b = 542 meV; and 98 half of the ECG result.The previous CI calculation 18 is sim 99 the virtual-positronium formation correlation potential (GW 100 more than doubling it (to " b = 1290 meV).The additional 101 causes a relatively small reduction (to " b = 1095), and inclu 102 3 ladder series of (screened) electron-positron interact 79 It is unique to the positron problem because successi 80 molecule self energy with the same sign, whereas fo 81 We also consider the smaller positron-hole ladder ser 82 details of the state-of-the-art construction of the sel

93 3 Fig. 4 .
Fig. 4. Calculated electron-positron contact densities and positron lifetimes with respect to annihilation.a, fractional contribution of individual MOs to the total electron-positron contact density (Eqn.6 in 'Methods).b and c, the electron-positron contact density (magenta) at the Σ GW +Γ+Λ level for the HOMO and (H-1)OMO in acetonitrile (blue and brown show positive and negative electron wavefunction regions, respectively), c.f., Fig. 2 (b) (positron density).Table:Positron lifetimes with respect to annihilation.τ (0) : lifetime calculated in the Hartree-Fock independent particle approximation excluding the vertex enhancement factors and using a positron wavefunction normalised to unity; τ GW and τ : lifetime calculated using the Dyson positron wavefunction at the Σ GW and Σ GW +Γ+Λ levels including vertex enhancement factors and renormalization constants (Eqns.6 and 7 of 'Methods').
Solution of the Dyson equation also yields the positron-bound state wavefunction ψε.Using it, the 2γ annihilation rate in the bound state Γ = πr 2 0 cδep (Γ[ns −1 ] = 50.47δep[a.u.]), whose inverse is the lifetime of the positronmolecule complex with respect to annihilation, can be calculated.Here r0 is the classical electron radius, c is the speed of light and δep is the electron-positron contact density δep = Ne n=1 γn |ϕn(r)| 2 |ψε(r)| 2 dr.

Fig. 4 and 1 .
Extended Data Figs. 4 present contact density data.Extended Data Fig. 4 a shows the individual MO contribution to the contact density as a function of the MO ionisation energy: similar to Fig. 3 of the main text (contribution of strength parameters from individual MOs), overall the contact density increases as the ionisation energy decreases: the positron overlap is greater with the more diffuse electronic HOMOs.However, MOs below the HOMO can in fact dominate, e.g., acetonitrile, as shown further in Fig 4 a, b and c.The main contributions to the positron-molecule self energy, including the twoparticle propagators.a, the GW diagram involves the positron Green's function Gν and the dynamic part (due to the absence of an electron-positron exchange interaction) of the screened Coulomb interaction W d = vΠv (bold denotes operator form), where Π is the electron-hole polarization propagator (see b).It satisfies the Bethe-Salpeter equation (diagram e) with kernel K = v − W RPA (diagram f), where W RPA = v + W d,RPA is the screened electron-hole Coulomb interaction calculated in the random phase approximation.Setting K = 0 results in the bare polarisation entering W only, and gives the Σ (2) approximation, so-called as it is a second-order diagram in the electron-positron Coulomb interaction.

Extended Data Figure 2 .
Fig.1 b), where W d = vΠ RPA v is the dynamic part determined from the polarization propagator in the random phase approximation, are determined asW µn,µ m = (µn|µ m)+ α w α µn w α µ m (ω − Ω α + + iη) −1 − (ω + Ω α + − iη) −1 .Convergence of positron binding energies in acetonitrile and CSe 2 with respect to electron and positron basis size.Positron binding energy calculated using the Σ GW @BSE+Γ+Λ self energy for varying number of electron (positron) HF MOs (whose energies are shown as blue and red crosses, respectively) included in the basis.For acetonitrile, the varying electron (positron) MO calculations included all positron (electron) MOs.For CSe 2 , the varying electron MO calculations included all positron MOs, whilst the varying positron MO calculation included 113 electron MOs (indicated by the lowest blue circle).The binding energy reaches convergence when the electronic orbital with energies up to ∼150-200 eV are included.Similar behaviour was also observed for the other molecules considered.ε b (meV), using Σ ≈ gΣ(2)   + Σ (Λ) g ≡ S (2+Γ) /S(2)

Extended Data Figure 4 .
Calculated electron-positron contact density.a, contact density for individual electronic MOs as a function of their ionisation energy, calculated including vertex enhancement factors and renormalization coefficients (see 'Methods' Eqns.6 and 7).Red dashed line: positronium ground state energy at |E Ps | = 6.8 eV.Grey line: δep = 0.008/(I − |E Ps |) (for a guide).Also see Extended Data Table

Table 1 :
Calculated positron-molecule binding energies (meV) Dipole moment µ from[36]; isotropic polarizabilities α and ionization energies I calculated at the GW level (see Extended Data Table

Table 3 :
Positron lifetimes (ns) 84Positron binding in polar and non-polar mol 85 86 Table 1 shows our calculated binding energies at su 87 correlation potential (⌃ GW in its various levels, ⌃ G 88 interactions in the ladders denoted by tildes).Exte 89 gies, and isotropic and anisotropic polarisabilities.

Table 3 :
Positron lifetimes (ns) 84Positron binding in polar and non-polar molecules.85 86 Table 1 shows our calculated binding energies at successiv 87 correlation potential (⌃ GW in its various levels, ⌃ GW + ⌃ 88 interactions in the ladders denoted by tildes).Extended D 89 gies, and isotropic and anisotropic polarisabilities.

Table 3 :
Positron lifetimes (ns) 84Positron binding in polar and non-polar mole 85 86Table 1 shows our calculated binding energies at suc 87 correlation potential (⌃ GW in its various levels, ⌃ GW 88 interactions in the ladders denoted by tildes).Exten 89 gies, and isotropic and anisotropic polarisabilities.F 90 Benchmarking and the general trends.-Wefirs 91 explicitly correlated gaussian (ECG) calculation (" 92 LiH.The results demonstrate the general trend se

Table 2 :
(2)E /∂E) −1 |ε µ .61ExtendedDataPositronbindingenergies in the GW approximation (meV), dimensionless correlation-potential strength parameters and Dyson wavefunction renormalisation constants a Positron binding energies (in meV, complementary data to Table1of main text) calculated at the HF and various levels of the GW approximation (see Extended Data Fig.1): Σ(2)(bare polarisation propagator); RPA (random phase approximation); TDHF (time-dependent Hartree Fock approximation); and BSE (Bethe-Salpeter equation).Dimensionless strength parameter of the correlation potential (defined in the main text) in different approximations to the positron-molecule self energy (see Fig.3of main text).Positive (negative) strength parameters denote attractive (repulsive) positron-molecule interactions.The final column gives the calculated positron Dyson wavefunction renormalization constants a for the Σ GW +Γ+Λ calculation (see 'Methods' Eqn. 7).