A Regularized Second-Order Correlation Method from Green’s Function Theory

We present a scalable single-particle framework to treat electronic correlation in molecules and materials motivated by Green’s function theory. We derive a size-extensive Brillouin-Wigner perturbation theory from the single-particle Green’s function by introducing the Goldstone self-energy. This new ground state correlation energy, referred to as Quasi-Particle MP2 theory (QPMP2), avoids the characteristic divergences present in both second-order Møller–Plesset perturbation theory and Coupled Cluster Singles and Doubles within the strongly correlated regime. We show that the exact ground state energy and properties of the Hubbard dimer are reproduced by QPMP2 and demonstrate the advantages of the approach for larger Hubbard models where the metal-to-insulator transition is qualitatively reproduced, contrasting with the complete failure of traditional methods. We apply this formalism to characteristic strongly correlated molecular systems and show that QPMP2 provides an efficient, size-consistent regularization of MP2.


INTRODUCTION
The characteristic divergences of second-order Møller−Plesset perturbation theory (MP2) and CCSD for typical condensedmatter and molecular systems represent a significant unsolved problem in electronic structure theory. 1−3 A scalable, quantitatively accurate, and nondivergent electronic structure theory for strong correlation remains elusive. While CCSD is a successful many-body theory of electronic correlation, the origin of its divergences stems from its dependence on the Hartree−Fock (HF) reference determinant. 4 The divergences of MP2 theory for insulators arise predominantly due to the denominator of the correlation energy. At internuclear separations where the HF singleparticle energies are similar and the coulomb matrix elements remain finite, the theory will clearly diverge. Historically, the issue of this divergence has been handled by introducing various regularization schemes. The κ-regularization scheme of Lee and Head-Gordon has attracted a lot of interest  (1) where {ϵ p } are the Fock energies, i and a index occupied and unoccupied Fock orbital states, ⟨ij∥ab⟩ is the antisymmetrized Coulomb integral, and κ is empirically determined though benchmarking data sets. 5 This regularization scheme is useful in curing the divergent behavior of MP2 when the divergences are caused by the vanishing denominator Δ ij ab . However, the κregularization scheme diverges for the Hubbard model and therefore cannot account for the qualitative metal-to-insulator transition observed as U t increases (Section 6). 6 Another approach to regularization takes the form of Brillouin-Wigner perturbation theory (BWPT) which is exact for two-level systems and avoids the divergence of the MP2 denominator. 7 However, BWPT is not size-extensive and gives multiple solutions which are nondesirable characteristics.
The Driven Similarity Renormalization Group (DSGR) method of Evangelista and co-workers can also be viewed as another size-consistent regularization of second-order perturbation theory. 8−11 At second-order in perturbation theory of the similarity transformed Hamiltonian obtained from the DSRG transformation, the correlation energy is given by 8 where s ∈ [0, ∞) is the flow parameter which is related to an energy cutoff, = s 1/2 . The flow parameter essentially decouples states of the Hamiltonian which are separated by an energy difference, Δ ij ab larger than the energy cutoff, thereby alleviating the intruder state problem. 12−15 The advantages of the DSRG method stem from the fact that it can be formulated as a modified unitary coupled cluster theory which allows for iterative solution of the flow equations thereby including interactions beyond second-order perturbation theory. 8 Additionally, the method can also be formulated based on a generalized reference wave function referred to as DSRG-MRPT. 11 However, while the regularized second-order perturbation theory in expression 2 (DSRG-PT2) cures the divergent behavior of MP2 theory caused by the vanishing denominator, it also diverges akin to the κ regularization scheme for the metal-to-insulator transition present in the Hubbard model as U t increases. Motivated by these issues, we introduce a size-consistent and natural second-order regularization scheme based on the single-particle Green's function. This implementation is rooted in the second-order Goldstone self-energy which we refer to as quasi-particle MP2 theory (QPMP2). Importantly, this provides us with a first-principles regularization which recovers the performance of MP2 in the weakly correlated regime, while extending its domain of applicability to the strongly correlated regime. Our results show that QPMP2 recovers the metal-toinsulator phase behavior of various Hubbard models and extends the domain of MP2 theory for strongly correlated molecular systems. However, we find that for Green's functions methods in general to be successful in quantum chemistry requires higher-order interaction terms to be taken into account.
The Green's function formalism represents an effective Hamiltonian theory resulting from a field theoretic approach to the electronic structure problem. 16−18 By working directly with the second-quantized electron field operators themselves, the Green's function appears naturally when computing expectation values of many observables. The advent of quantum field theory and quantum electrodynamics introduced many sophisticated methods for evaluating the single-particle Green's function based on time-dependent perturbation theory. 19−23 By folding in the large numbers of degrees of freedom into the dynamical self-energy, a single-particle potential which contains all the many-body interactions present in the system, the many-body problem is expressed entirely in terms of a "one-particle" theory. Of course, there are many other effective Hamiltonian approaches operating in both the single-particle and many-body formalism, which have been explored extensively in the literature. 24−29 In this work, we focus on the second-order Green's function perturbation theory approach, where the dynamical self-energy is approximated to second-order in the interaction, as a means of regularizing the MP2 correlation energy. The second-order, self-consistent Green's function perturbation theory (GF2) has also been shown to cure the divergences of MP2 in a size-consistent, single-particle framework and has been extremely successful in providing a self-consistent treatment of strongly correlated systems. However, this approach is implemented within a temperature dependent Matsubara frequency formulation and is a much more involved numerical implementation deriving from statistical mechanics. 30−34 The ground state correlation energy evaluated within the GF2 implementation of Zgid and co-workers is explicitly temperature dependent and obtained by evaluation of the Matsubara frequency summation of a modified Galitskii-Migdal (GM) formula dependent on both the real and imaginary parts of the Green's function and self-energy. 35 The implementation is performed entirely within the atomic orbital basis and by building the self-energy on the imaginary time grid which scales with the number of imaginary time grid points used. The subsequent Fourier transformations and selfconsistency requirements result in an increase in the formal scaling of the GF2 implementation as n n ( ) where n τ is the number of imaginary time grid points and n at is the number of atomic orbitals. The auxiliary GF2 method provides an alternative approach to the renormalized second-order self-energy by exploiting the fact that the effects of the self-energy can be expressed by coupling the system to a set of fictitious external degrees of freedom spanned by the so-called "auxiliary" states. 36 This approach reformulates the Dyson equation in the language of wave functions and avoids the use of frequency grids and numerical Fourier transforms. 36,37 The ground state correlation energy is self-consistently evaluated using a modified Galitskii-Migdal formula at zero temperature to generate a renormalized MP2-like correlation energy derived from the auxiliary degrees of freedom. The regularized MP2 correlation energy introduced in auxiliary-GF2 is obtained by identifying that directly applying the GM formula gives twice the MP2 energy for the noninteracting Green's function. 36 However, this regularized correlation energy is not directly obtained from the self-consistent Green's function formalism, but is rather obtained by halving the two-body correlation energy obtained from the GM formula using the "interacting' Green"s function containing the auxiliaries. Purely using Green's function auxiliary compression means that each iteration, until selfconsistency is obtained, would scale as N ( ) phys 9 where N phys is the number of spin−orbitals in the physical subspace. With a two-step auxiliary compression scheme (which introduces a small error in the true correlation energy obtained from the full procedure) this scaling can be formally reduced to N ( ) phys 5 , that is, the scaling of regular MP2 theory at each iteration until self-consistency is obtained.
Our contribution connects BWPT and MP2 regularization methods with the second-order frequency dependent selfenergy in a clear and simple way. By employing the quasiparticle approximation for the Green's function and introducing the forward time second-order Goldstone self-energy, we are able to show directly how the regularized MP2 correlation energy naturally appears. This sheds light on how second-order approximations for the self-energy act to regularize the MP2 correlation energy and allows us to investigate the effects of the quasi-particle approximation on ground state correlation energies. Our approach considers the effects of the dynamical second-order self-energy only on the occupied states by noting the similarity of the ground state density obtained from the quasi-particle Green's function with that of Kohn−Sham Density Functional Theory (KS-DFT). Our QPMP2 method takes the quasi-particle solution from the second-order Goldstone self-energy and uses these renormalized singleparticle energies to calculate the electronic correlation from the associated Green's function. From this approach, we are able to cure the divergences present in MP2 while remaining within a scaling of N ( ) 5 , where N is the number of spin−orbitals. Our regularization scheme is particularly adapt for Hubbard model systems, where both the κ and DSRG-PT2 regularization schemes as well as CCSD diverge. This is a dramatic result which demonstrates the regularization of second-order perturbation theory afforded by QPMP2 across the entire correlation regime. Our results for molecular systems show that this simple method provides similar performance to MP2 within the weakly correlated regime, but extends its domain of applicability when the MP2 solution diverges in the strongly correlated regime. To the best of our knowledge, this approach to regularization has not yet been explored, but is closely related to GF2 and other second-order regularization schemes. 6,38 The structure of this paper is as follows. Section 2 provides background on the single-particle Green's function and its relationship to the ground state energy. Section 3 outlines the structure of quasi-particle Green's function theory and the solutions of the quasi-particle equations. Section 4 provides the definition of the Goldstone self-energy which contains the interactions encoded by QPMP2 theory. Section 5 outlines how our approach can be viewed as a size-extensive Brillouin-Wigner perturbation theory. In Sections 6 and 7 we apply QPMP2 to Hubbard models and molecular systems, comparing with other electronic structure theories. Finally in Section 8, we summarize the results of our method and outlook for future scalable electronic structure theories.

SINGLE-PARTICLE GREEN'S FUNCTION AND GROUND STATE ENERGY
The second-quantized, nonrelativistic electronic structure Hamiltonian in atomic units and within the Born− Oppenheimer approximation is given by 16 where x 1 ≡ (r 1 , σ 1 ) is the composite spin-space coordinate and ψ † (x)/ψ(x) are the field operators that create/annihilate electrons at the composite spin-space coordinate x. The sets {r i } and {R I } refer to the electron and nuclear coordinates, respectively. {Z I } is the set of nuclear charges associated with the nuclei. The single particle Green's function is defined by the following time-ordered expectation value where |N, 0⟩ is the normalized, exact N-electron ground state and T is the time-ordering operator which places operators with the larger time argument to the left. The field operators are defined within the Heisenberg picture. From the Green's function, it is possible to find the ground state expectation value of any single-particle operator where η tends to zero from above. From the equation of motion for the Green's function it is also possible to extract the total ground state energy In general, the Green's function can be constructed from its Lehmann representation as where ψ j (x) = ⟨N − 1, j|ψ(x)|N, 0⟩ and ψ a (x) = ⟨N, 0|ψ(x)|N + 1, a⟩ are the "single-particle" Dyson orbital states.
are the exact electron addition and removal energies, respectively. From the Lehmann representation, the Green's function reduces to the particle density via Importantly, this relation is noninvertible and thus it is impossible to obtain the Green's function from the density. It is clear from above that the ground state energy of a system of N-interacting electrons can be expressed in terms of the Dyson orbitals as 39 In general, the nonorthogonal and overcomplete set of Dysonorbitals are normalized by their pole strength or probability factors, P i The probability factors are constrained between zero and one and the summation requires that the minimum number of ionization poles any Green's function possesses must be equal to the number of occupied spin−orbitals. In order to capture all the many-body correlation effects within a single-particle framework, it is necessary for the dimension of the Dyson spin−orbital basis to be larger than that of the spin−orbital one due to the one-to-one mapping between ionization processes and Dyson orbitals. The probability factors are the residues of the poles of the single-particle Green's function. In the HF approximation, these residues are all normalized to one for every pole of the Green's function and the number of ionization poles is equal to the total number of electrons. Under conditions where there is a dominant solution of the Dyson equation, one with the largest residue, the Green's function can be approximated by simply just taking this pole.
In this case, we refer to any such Green's function as a quasiparticle Green's function. From eq 8, the ground state electronic density is a functional only of the "occupied" Dyson orbitals {ψ i }, which is analogous to the expression for the density in KS-DFT. When working with the quasi-particle Green's function, the expression for the density is exactly that obtained from KS-DFT (a sum over N orthonormal spin− orbitals).

QUASI-PARTICLE GREEN'S FUNCTION THEORY
The equation of motion for the Green's function can be transformed into the Dyson equation. Written in the spin− orbital basis and frequency domain this takes the form of where Σ(ω) is the dynamical self-energy containing all the many-body interactions experienced by a single-particle and G 0 (ω) is a reference Green's function which we take here to be the HF propagator. This equation of motion can be recast as the quasi-particle equation where f is the Fock operator. The poles of the interacting Green's function are given by the excitation energies ε p which are the energies at which the inverse single-particle Green's function vanishes. As a result, the solutions of eq 12 yields the exact poles and residues of the single-particle Green's function, meaning that we need only find these solutions to completely specify the single-particle Green's function (eq 7). To find the quasi-particle solutions ({ϵ p QP }), which are those of the largest residue, we can express the quasi-particle equation in terms of a perturbation theory with respect to the Fock orbitals For each Fock wave function we define two projectors where f |ϕ p ⟩ = ϵ p |ϕ p ⟩. Resolving the identity onto the quasiparticle wave function, we have We can normalize the Fock single-particle state with the quasiparticle state by letting ⟨ϕ p |ψ p ⟩ = 1 as |ϕ p ⟩ and |ψ p ⟩ correspond to each other. Treating the self-energy as a perturbation to the Fock operator, we apply Q p to both sides of eq 13 and commuting it past f gives (16) where λ ∈ [0,1] is the usual perturbation parameter. Therefore, from eq 16, we have the Dyson equation for the quasi-particle wave function The resulting equation for the quasi-particle energy is given by which, using the normalization above gives The self-energy can also be expanded perturbatively in the order of interaction (Feynman diagram series) as follows where Σ (1) = 0 as the zeroth-order Green's function corresponds to the Fock operator. Inserting this expression into eq 19 gives ... is simply a restatement of the Dyson geometric series summation of eq 11 as a function of the self-energy perturbation. As the self-energy contains only connected diagrams, no reducible diagram contributions appear. Expanding the self-energy to second-order in the interaction we have which, by setting λ = μ = 1 gives We have converted the quasi-particle equation into a Brillouin-Wigner perturbation theory for the corrected quasi-particle energies from the Fock reference states. The second-order quasi-particle equations are written as which is equivalent to the diagonal approximation of the Dyson equation in eq 11. The diagonal approximation is regularly used in GW band structure calculations to compute the corrected addition and removal energies from a DFT calculation. 41,42 In this case, the self-energy is given by the GW approximation with the polarization calculated within the Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article random phase approximation (RPA). 43−45 However, usually eq 24 is solved in a linearized form by Taylor expansion of the GW self-energy about the KS-DFT energy which also introduces renormalization factors. 46,47 Additionally, these calculations are not concerned with the implications for the ground state correlation energy. 48 The solutions of eq 24 within the GW approximation have been explored extensively in the form of evGW and have been shown to be comparable to those of the fully self-consistent qsGW. 47 Importantly, the self-energy appearing in eq 24 is expanded as the Feynman- , not in terms of the renormalized self-energy functional appearing in Hedin's equations (Σ[G, W]). Our choice to normalize the HF single-particle state with the quasi-particle state means we have made a quasi-particle transformation where there is a one-to-one correspondence between the HF orbitals and the quasi-particle orbitals. Taylor expanding the self-energy about the HF energy gives (2) QP p (25) and taking the zeroth order term gives We shall see that this expression simply gives the MP2 correlation energy with the Goldstone self-energy (Section 4). However, working with expression 24 will give us a welldefined, simple, and size consistent expression for the ground state correlation energy: QPMP2. In practice, solving for the quasi-particle energy in eq 24 simply amounts to finding the largest root. By taking the quasi-particle solutions only, the Green's function represents a canonical transformation. This means that we need only find the root with the largest residue which can easily be done via an iterative solver. This substantially simplifies the issue of having to calculate all the solutions of this nonlinear equation and significantly improves the scaling. Inclusion of all the satellite solutions of the quasi-particle equation corresponding to less significant many-body interactions and their residues constitutes a significant increase in computational effort. We show from our Hubbard model results that inclusion of these satellites, required to construct the full Green's function, makes little to no impact on the ground state correlation energy within the second-order approximation.

THE GOLDSTONE EXCHANGE-CORRELATION SELF-ENERGY
In this section we define the second-order self-energy used to solve the quasi-particle eq 24. Our choice of self-energy provides a straightforward connection to MP2 regularization methods and size-consistent Brillouin-Wigner perturbation theory. This leads us to the definition of the second-order selfenergy in eq 27 which enters the QPMP2 correlation energy derived in Section 5. We can generate a perturbation expansion for the self-energy of occupied states by expanding the exchange-correlation energy from the infinite Goldstone diagram series. This is analogous to the Luttinger-Ward functional introduced in the temperature dependent formalism. 49,50 The Goldstone dia-grams are defined with a prespecified time direction with arrows moving upward corresponding to virtual states (particles), while those moving downward are occupied states (holes). 16 We define the self-energy as the functional derivative of the Goldstone correlation energy with respect to the hole Green's function. The Goldstone self-energy is identical to the forward time Feynman-Dyson self-energy up to third order. 51 The self-energy of the occupied states simply corresponds to cutting the downward arrows of the Goldstone diagrams. This definition means that we can regain the Møller−Plesset perturbation theory expansion by simply evaluating the trace of the Goldstone self-energy at the single-particle HF energies. 51 At second-order, the Goldstone self-energy is given by cutting the downward arrows of the direct and exchange diagrams to yield the expression (Figure 1 2 (27) This also ensures a real-valued self-energy with no lifetime effects. From this expression, we directly see that the selfenergy is a frequency dependent single-particle potential containing many-body interactions present in the system. We regard this definition of the self-energy as being the most useful in the context of ground state electronic correlation. In Feynman-Dyson perturbation theory, the self-energy is a complex time-ordered quantity with real and imaginary parts. 52 This is the quantity that is approximated in the GF2 and auxiliary GF2 methods. However, the presence of a complex self-energy no longer guarantees real energy solutions of the resulting quasi-particle equations. The real part of these solutions is usually interpreted as the renormalized singleparticle energies while the imaginary part is interpreted as an exponential lifetime of the state. 45,46 However, using the Goldstone self-energy, we have a zero imaginary part� precisely the conditions under which Brillouin-Wigner perturbation theory is valid. 17

QUASI-PARTICLE MP2 AS A SIZE-EXTENSIVE BRILLOUIN-WIGNER PERTURBATION THEORY
In this section, we derive the QPMP2 correlation energy from the quasi-particle Green's function and outline how QPMP2 can be interpreted as a size-consistent formulation of secondorder BWPT. The central result of this section is the QPMP2 and iQPMP2 correlation energy derived in eqs 35 and 37.
In an orthornormal spin−orbital basis, the exact Green's function is defined as  where C p i = ⟨ϕ p |ψ i ⟩ and D p i = ⟨ϕ p |ψ a ⟩ are the overlap between the Dyson orbitals and the occupied/virtual HF states and ε i/a are the exact electron removal and addition energies defined in Section 2. The ground state energy is given by the contour integral in eq 6 which evaluates to The residues/probability factors are simply given by By choosing the quasi-particle normalization condition, we have that as the set of quasi-particle states now constitute a complete set and the transformation is unitary. These states are obtained as the single largest solution of the quasi-particle equation as outlined in Section 3. Therefore, the quasi-particle GF energy is given by Using eq 24 as a very good first approximation, we see that the quasi-particle energy is given by Therefore, the quasi-particle GF energy reduces to From Section 4, the exchange-correlation self-energy can be systematically approximated according to the Goldstone diagram series (eq 27). The second-order expression of eq 27 gives the single-shot quasi-particle MP2 correlation energy (QPMP2) as This expression is clearly size extensive and represents a regularization of the MP2 expansion based on a Green's function approach. Although not pursued here, this expression can be trivially extended to finite temperature by working in the Matsubara formalism which introduces the relevant Fermi−Dirac occupation factors. Expression 35 is symmetric with respect to particle exchange as the summation runs over both occupied dummy indices i, j. It is cast in the canonical Hartree−Fock orbital basis and while an orbitally invariant formulation is possible, the resulting equations are much more involved. The MP2 correlation energy corresponds to taking the quasi-particle energy from the Møller−Plesset perturbation expression in eq 26. Therefore, for our definition of the selfenergy, all higher order MPn correlation energies can also be found by evaluating the trace of the n th -order Goldstone selfenergy at the HF energies.
Furthermore, one may iterate the self-energy as to include the effects casued by the renormalized occupied states. 16−18 Updating the second-order MP self-energy to account for the quasi-particle renormalization gives the expression for the interacting QPMP2 self-energy This renormalization is represented in Figure 2, where the hole lines are now dressed. By the same procedure, we may take the trace of the updated self-energy evaluated at the quasi-particle energy to give the interacting QPMP2 (iQPMP2) correlation energy as This represents no increase in scaling as the quasi-particle energies are already found by solution of the quasi-particle eq 33. As will be shown, this expression substantially improves the dissociation curves obtained for molecular systems, but reduces the performance of QPMP2 for Hubbard models. Many-Body Brillouin-Wigner perturbation theory is constructed by defining the following projector where E 0 is the exact ground state energy. The perturbative expression for the ground state energy within this formalism is given by Expanding self-consistently to second-order for the correlation energy gives This is the many-body equivalent of the Dyson series, but it is not clear what the energy difference of the denominator should be in order to decouple these equations. Additionally, the expression is not size-extensive and is evidently problematic for extended systems. 6 There have been a few attempts to introduce a size-extensive implementation of second-order BW perturbation theory via ad-hoc procedures of defining this energy difference as a correlation energy per electron. From the size-extensive, interacting quasi-particle MP2 correlation energy, we can view the combination of all the energy  This is exactly the expression for the interacting quasiparticle MP2 correlation energy. Hence we can view QPMP2 as a size extensive Brillouin-Wigner perturbation theory.

RESULTS FOR HUBBARD MODELS
The Hubbard Hamiltonian has been used for decades as a conceptually simple, yet insightful model to probe the nature strong electronic correlation and the metal-to-insulator transition. 53,54 The 1D periodic Hubbard model consists of a lattice of sites, each of which support a maximum of two electronic states (one spin up, one spin down). The model consists of a "hopping" term which couples nearest neighbor sites (favoring electron delocalization) and a repulsive interaction when two electrons are on the same site (favoring electron localization). The two terms provide an intuitive model of the behavior of electrons in molecules and materials. In the site basis, the Hamiltonian is (42) where imposing periodic boundary condition means c N+1σ † = c 1σ † . The Hamiltonian can be parametrized by the ratio U t , which provides a universal measure of the "correlation strength". As U t increases, the on-site repulsion increases relative to the "hopping" term, favoring configurations where the electron density is localized about each site. By neglecting the correlated fluctuations in the electron density about each site, we obtain the HF Hamiltonian for the Hubbard model (43) where = N n i i is the total number operator, and N s is the number of sites. This will always be a linear function of U t .

The Hubbard Dimer.
The half-filled Hubbard dimer provides a simplified model of the hydrogen molecule in a minimal basis. It is an insightful system to understand the nature of strong correlation present in the chemical bond and there have been numerous studies of the Hubbard dimer as a benchmark for electronic correlation methods. 55 To find the exact Green's function requires the diagonalization of the single and three particle Hamiltonians as well. Using these states, it is simple to construct the exact Green's function and polarization from their Lehmann representations. As the Hamiltonian is spin-independent, all expressions are spindiagonal and the spin-index can be dropped. From knowledge of the exact and HF Green's functions, we can invert Dyson's equation to find the exact exchange-correlation self-energy. In the diagonal plane-wave basis, the exact Feynman-Dyson selfenergy is given by Solving the quasi-particle equations for this self-energy and taking the solutions with the largest weight gives the QPMP2 quasi-particle energies The corresponding QPMP2 ground state energy is given by which is also exact for the Hubbard dimer as Therefore, QPMP2 also reproduces the exact ground state energy of the Hubbard dimer. The MP2 correlation energy can be found by simply evaluating the second-order Goldstone selfenergy at the Fock single-particle energy and diverges due to the multireference nature of the strongly correlated solution. 35 The CCSD equations can be solved analytically, reproducing the exact ground state wave function and energy over the entire correlation regime as it is exact for two electron systems.

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article Figure 3(a) depicts the ground state energy of the Hubbard dimer as a function of U t for a variety of different electronic structure methods. We compare QPMP2 to regular GF2, RPA, GW, GW-BSE, CCSD, MP2, and HF theories. These results provide a qualitative insight into the relative behavior of these methods when accounting for bond dissociation processes.
The RPA energy provides a slight improvement over the HF approximation for values of < 1.8

U t
, but the presence of complex excitation energies for 2 U t results in a discontinuous curve. 1,56 The sudden onset of complex excitation energies as the correlation strength increases can be related to the change in the character of the HF Hessian matrix. 61 The RPA matrix is the HF Hessian matrix and if this matrix is positive-definite, then all the excitation energy eigenvalues will be real. Therefore, the presence of complex excitations is a consequence of the fact that the HF solution is no longer at a stable minimum due to the multireference character of the correlated ground state and is related to the singlet−triplet instability. 62 The divergent behavior of RPA when accounting for ground state electronic correlation is a consequence of the fact that artificial boson commutation relations are imposed on the Fermion electron−hole pair operators.
The GW approximation sees a significant improvement over HF, by correlating single-particle states via the "first-order" renormalized self-energy of Hedin. 43,63,64 The GW energy is reasonable within the weak correlation regime (up to values of 2.6 U t ) but fails in the strong-correlation limit, displaying an approximately linear correlation energy. 56 Figure 3(b) displays how the weight of the quasi-particle and satellite solutions vary as the correlation strength increases. It is clear that the weight of the quasi-particle state is dominant over the entire regime for both GW and the exact system. Figure 3(c) shows the variation in the quasi-particle energies with the correlation strength. From this plot we see that the band gap increases as U t increases for both QPMP2 and GW. The "exact" band-gaps change much more significantly as a function of U t when compared to both the GW and HF approximations. We also see that the variation in the energy of the satellite solutions across the entire correlation regime is much less significant than for that of the quasi-particle solutions. The entire spectrum of single-particle states is

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article symmetric about the chemical potential due to particle-hole symmetry. Similar to the relationship between RPA and HF, the GW-BSE method is dependent on the nature of the GW solution. It shows a marked improvement over the GW approximation, matching the exact energy over the region where GW performs well, and is accurate up to < 7 U t . However, in the strong correlation regime, the GW-BSE excitation energies become complex as the GW solution becomes unstable. The RPA and GW-BSE approaches give complex excitation energies in the strong correlation limit where the HF/GW starting point is no longer at a stable minimum. 65,66 Therefore, we suspect that the GW-BSE method applied to real molecular systems will likely exhibit similar divergences to those of RPA for the ground state energy.
The GF2 and QPMP2 approximations are exact for the Hubbard dimer. The second-order approximation corresponds to an infinite order resummation of the MP2 diagrams via the Dyson equation. The fact that QPMP2 also reproduces the exact ground state energy is even more intriguing as it suggests that the Hubbard dimer quasi-particle solutions can be considered as "exact" single-particle states, and that the satellites do not contribute to the energy. Interestingly, the interacting QPMP2 expression is no longer exact for the Hubbard dimer system and substantially decreases the accuracy of the approach, resulting in a correlation energy that is similar to that obtained from the GW approximation.
The κ-regularization scheme also diverges akin to MP2 for the Hubbard dimer. Figure 3(a) displays the divergence of κ-MP2 for a range of κ values from 0.8 to 1.4. From eq 1, we see that as κ tends to zero we recover the HF energy and as κ → ∞ we obtain the MP2 correlation energy. We see that no typical value of κ will regularize MP2 for the Hubbard dimer and all regularization schemes based on this approach fail for larger Hubbard models (Section 6.2). This is because the divergence of MP2 for the Hubbard model is actually due to the numerator rather than the denominator as Δ ij ab is constant over the entire range of U t . Similar results are observed for the second-order perturbative DSRG expression of eq 2. Therefore, one must be careful when using Hubbard models to understand electronic correlation in molecules and materials.

Periodic Six-Site Hubbard Model.
The half-filled periodic six-site Hubbard model is a closed shell system where a single HF reference provides a reasonable starting point. , where the electronic motion is no longer weakly correlated, and drastically fails at larger correlation strength. Of the traditional post HF electronic structure methods widely used, it would appear that MP2 performs best. However, the MP2 energy diverges, due to the highly multireference nature of the FCI wave function and begins to diverge for 4 U t . This is verified in Figure 4(b), as the quasi-particle weights from the GF2 calculation also begin to change significantly around 4 U t . Again, κ-MP2 also diverges for the six-site model and does not provide a regularization of MP2. The curves show that ultimately the κ-MP2 regularization scheme provides little to no improvement over MP2 for the Hubbard model. Again, this divergence is due to the numerator rather than the denominator of the MP2 expression. Hence, the Hubbard model represents a signifi-cantly challenging condensed matter physics system for benchmarking regularization schemes and traditional electronic structure theories. 68 CCSD performs well in the weak correlation regime for values of 4 U t . However, it also diverges in the intermediatestrong correlation limit as the weight of higher excitations in the FCI begins to dominate the doubles amplitude. This is a spectacular failure as CCSD does not provide a qualitatively accurate description of the metal-to-insulator transition present in the FCI. 1,35,69 The GF2 and QPMP2 methods give the same ground state energy for the six-site model across the entire correlation regime. Both perform similarly to CCSD and MP2 in the weak correlation regime as all methods begin to move away from the FCI energy at 4 U t . However, the ultimate difference is that neither the GF2 or QPMP2 energy diverges in the strongly correlated regime and therefore, they provide a qualitatively accurate description of the strongly correlated system. As a result, the metal-to-insulator transition is captured and the divergences of MP2 are cured. Again, it is noted that the interacting QPMP2 correlation energy expression 37 decreases the accuracy of the single-shot QPMP2 expression for the sixsite Hubbard model. Figure 4(b) contains the variation of the GF2 quasi-particle weights as a function of the correlation strength. Through the variations in the quasi-particle weights as a function of correlation strength, we can understand at the nature of the different quasi-particle solutions and the QPMP2 correlation energy. As can be seen, the weights of the satellite solutions begin to become increasingly important for 8 U t . This is the hallmark of strong correlation, where many-body effects are increasingly present. However, the weight of the first satellite solutions is effectively zero across the entire correlation regime and can be effectively ignored altogether. The quasi-particle energies in Figure 4(c) are again symmetric about the chemical potential as the self-energy preserves the particle-hole symmetry. Their variation with the correlation strength again highlights the increase in the band gap as U t increases. It is the correlation embedded in these quasi-particle states which leads to the QPMP2 regularization of MP2 for the Hubbard model.
Analogous behavior is observed for the linear eight-and tensite Hubbard models which also display the characteristic divergences of CCSD, MP2, κ-MP2 and DSRG-PT2. Again, QPMP2 theory qualitatively reproduces the metal-to-insulator transition present in the FCI solution and is nondivergent over the entire range of U t . The performance of QPMP2 for the Hubbard dimer and larger Hubbard models indicate that this regularization is sufficient to cure the divergences of MP2 within the strongly correlated regime, providing a qualitatively accurate description of the metal-to-insulator transition.

RESULTS FOR MOLECULAR SYSTEMS
In this section, we apply QPMP2 to a series of molecular systems which demonstrate the transition from weak to strong correlation present in quantum chemistry. All dissociation curves level off at large distances unless stated otherwise. Figure 5(a) shows the ground state potential energy surface for H 2 in the cc-pVDZ basis set. Clearly CCSD is exact for any two-electron system and recovers the FCI dissociation curve and exact ground state wave function.
The well-known divergence of MP2 for the simplest molecular system is a result of energy denominator Δ ij ab tending to zero as the distance between the hydrogen nuclei increases. The dissociation curve for QPMP2 cures the divergence of the MP2 solution while showing an identical energy surface around the equilibrium configuration. The single-shot calculation is enough to open up the band gap and regularize the MP2 energy in the strongly correlated regime.
However, the method substantially underestimates the correlation energy in the dissociation limit for H 2 as it simply regularizes the MP2 correlation energy. This is in stark contrast to the case for the Hubbard dimer, where both QPMP2 and GF2 are exact over the entire range of U t . This again reiterates that using the Hubbard dimer as a simple benchmark system for electronic structure methods may not always be appropriate. The nonzero slope at dissociation of the QPMP2 dissociation curve (which levels off at large bond distances) is completely removed in the interacting QPMP2 correlation energy expression. The interacting QPMP2 energy introduces no additional scaling as all the quasi-particle energies have already been calculated from the quasi-particle equation and is as a result of "iterating" the self-energy once to include the effects of the renormalized hole states (Section 5). It has been shown in Section 6 that the interacting QPMP2 approach for Hubbard model systems reduces the performance of the QPMP2 correlation energy. However, we find this update of the self-energy to be necessary to correctly describe bond dissociation processes. The observed improvement resulting from iteration of the self-energy should be noted as a fundamental difference between the electronic correlation in real molecules and Hubbard models. We would like to comment on the similarity between our interacting QPMP2 expression 37 and the expression obtained in the recent work of Carter-Fenk and Head-Gordon on a size consistent Brillouin-Wigner correlation energy expression involving dressed occupied orbital eigenvalues. 70 In addition, the κ-MP2 ground state energy of H 2 is shown for a typical value of κ = 1.6 reported in the literature. 5,71 As we can see QPMP2, κ-MP2 and MP2 are essentially identical near the equilibrium nuclei configuration. However, as the correlation strength increases toward dissociation, both QPMP2 methods recover a much larger correlation energy. The divergences of κ-MP2 seen for the Hubbard models do not appear for H 2 as the divergence of MP2 in this case is due to the denominator Δ ij ab rather than the numerator. We also benchmark our method on characteristic systems for which CCSD and CCSD(T) diverge. The linear H 6 and H 10 chains as well as the cyclic H 6 ring are such systems where effects of strong correlation present difficulties. The spectacular failures of CCSD and CCSD(T) for these systems is concerning and highlights the current challenges strong correlation presents for scalable electronic structure theories.
For the linear H 6 system ( Figure 5(c)) we see the same divergent behavior of MP2 due to the vanishing band gap that is renormalized by QPMP2 theory. However, again the correlation energy is overestimated at "dissociation" and QPMP2 simply acts as a regularization of MP2. The interacting QPMP2 energy removes the nonzero slope of the QPMP2 energy at dissociation. Again, κ-MP2 also regularizes MP2 in contrast with the Hubbard model results of Section 6, but does not recover the correlation energy in the strongly correlated regime. CCSD diverges for both H 6 systems and linear H 10 due to the onset of strong correlation as the H−H distance increases. Figure 5(d) shows the ground state energy of the H 6 ring as a function of the breathing mode within the cc-pVDZ basis. We see similar quantitative results as for linear H 6 with the regularization scheme again providing a qualitatively accurate energy at "dissociation". However, again we see no improvement over MP2 about the equilibrium configuration as the curves are essentially identical. Results for linear H 10 (not shown here) also display the same behavior as the linear and cyclic H 6 systems. We note that the quantitative accuracy of CCSD for these systems is unrivalled within the weakintermediate correlation regime. Similar qualitative results for H 6 , H 12 , and H 32 lattices have been demonstrated by the GF2 method of Zigid et al. 31,32,35,72 However, their implementation is based on the finite temperature formalism, where the Green's function and selfenergy are constructed in the imaginary frequency and imaginary time domains.
At finite temperature, the Green's function is antiperiodic in imaginary time with a period of the inverse temperature β and can be expanded as a Fourier series over discrete Matsubara Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article frequencies. As a result, the ground state energy is expressed as a Matsubara frequency summation and their method utilizes several properties of the Green's function which are not present at zero temperature. 16 Additionally, their implementation is based on the renormalized Feynman-Dyson perturbation theory to perform self-consistent calculations which are much more computationally involved. Finally, we discuss our regularization for the ground state energy of N 2 which presents a significant challenge at dissociation due to the highly correlated nature of the triple bond. In Figure 5(b), we see that CCSD performs well until it diverges in the dissociation limit. MP2 theory fails completely for the N 2 molecule as it barely reproduces the minimum in the energy surface before it diverges. In contrast, κ-MP2 provides an efficient regularization of MP2, but is less accurate near the equilibrium configuration.
For N 2 , we see the biggest difference between the interacting and single-shot QPMP2 correlation energy. Single-shot QPMP2 substantially improves the qualitative performance of MP2, clearly displaying a minimum. However, it also diverges beyond 2 au as the nature of the electronic correlation exhibited by N 2 near dissociation is not described within this scheme. In stark contrast, iQMP2 completely cures the divergences of second-order perturbation theory for N 2 and provides an accurate potential energy surface which dissociates to the correct limit. From Figure 5(b), the dissociation curve of iQMP2 matches the FCI solution at around 1.9 au and outperforms CCSD within the weak-intermediate correlation regime while remaining nondivergent at dissociation. However, the potential energy of iQMP2 at dissociation is again overestimated. This is a spectacular result that demonstrates the power of our second-order regularization scheme for molecular systems.

CONCLUSIONS AND OUTLOOK
We have introduced a scalable, size-consistent regularized second-order correlation method based on the single-particle Green's function which can be cast as a size-extensive Brillouin-Wigner perturbation theory. We have also outlined its connections to different regularized expressions that have been explored in the literature, with the additional advantage that QPMP2 represents a physically meaningful regularization scheme. We have shown that QPMP2 is capable of curing the divergences of MP2 in the strongly correlated regime for a wide range of cases and that the method is capable of qualitatively describing the metal-to-insulator transition. However, gaining a quantitatively accurate description of strong correlation within the Green's function formalism is not straightforward and requires the rigorous inclusion of threebody interactions which will likely involve higher order Green's functions. Future work must be carried out to understand whether the quantitative accuracy required for quantum chemistry can be achieved within the Green's function formalism. There is currently no Green's function based approach in the literature that can compete with the quantitative accuracy of CCSD or CCSD(T) for ground state energies and dissociation curves. A major open question is whether the divergences of CCSD(T) can be cured by exploring the connections between the Green's function and Equation-of-Motion Coupled Cluster theory. 73 ■ AUTHOR INFORMATION