Wave-function inspired density functional applied to the H$_2$/H$_2^+$ challenge

We start from the Bethe-Goldstone equation (BGE) to derive a simple orbital-dependent correlation functional -- BGE2 -- which terminates the BGE expansion at the second-order, but retains the self-consistent coupling of electron-pair orrelations. We demonstrate that BGE2 is size consistent and one-electron"self-correlation"free. The electron-pair correlation coupling ensures the correct H$_2$ dissociation limit and gives a finite correlation energy for any system even if it has a no energy gap. BGE2 provides a good description of both H$_2$ and H$_2^+$ dissociation, which is regarded as a great challenge in density functional theory (DFT). We illustrate the behavior of BGE2 analytically by considering H$_2$ in a minimal basis. Our analysis shows that BGE2 captures essential features of the adiabatic connection path that current state-of-the-art DFT approximations do not.

A viable approach in functional development is to learn from wave-function theory in constructing nonlocal correlation functionals that involve unoccupied KS orbitals and thus stand on the fifth (and currently highest) rung of Perdew's Jacob's ladder [7]. A prominent example is second-order Görling-Levy perturbation theory (GL2) [20] that is closely related to secondorder Møllet-Plesset perturbation theory (MP2) [21][22][23] in wave-function theory [24]. In fact, even MP2 can be viewed as an implicit density functional by means of the adiabatic connection approach [20,25,26]. An important feature of 2nd-order perturbation theories (PT2) such as GL2 and MP2 is that they are one-electron "self-correlation" free [1,5], i.e. the correlation is zero for one-electron systems, and that they are size consistent (i.e. if the system is fragmented into two parts, the total energy becomes the sum of the two fragments) * zhang@fhi-berlin.mpg.de [27,28]. Furthermore, PT2 is fully non-local and captures the correct R −6 decay behaviour at long distances, which is essential to provide an accurate description of weak interactions. Therefore, PT2 correlation is an ideal building block for fifth-level density functionals, following Perdew's nomenclature. This feature has been exploited in a range of promising double-hybrid functionals [29][30][31][32][33][34][35][36][37], which linearly mix generalized gradient approximations (GGAs), e.g. BLYP [38,39] or PBE [40], with both exact exchange and PT2 correlation. The admixture of semi-local exchange and correlation can be viewed as an efficient yet semi-empirical way to take higher-order perturbative contributions into account that would go beyond PT2 [41,42]. These double-hybrid functionals provide a satisfactory accuracy for various chemical interactions, but they fail for systems with small KS energy gaps, e.g. heavily-stretched H 2 and metallic systems. In such systems, two or more determinants become degenerate in energy and mean-field theories that rely on a single reference determinant such as HF or KS-DFT break down. Also, perturbation theory diverges at any order, making it essential to find an appropriate resummation. Seidl, Perdew and Kurth suggested an empirical adiabatic-connection (AC) model [20,25,26], namely interaction-strength interpolation (ISI) [7], to implicitly resum the perturbation expansion by using only exact exchange, PT2 correlation and an explicit density functional derived from the "point charge plus continuum" model in the strong-interaction limit where the coupling constant parameter goes to infinity [43]. Frequently, the analytic dependence on the coupling constant is approximated by a Padé formula, whose parameters can be determined either empirically on theoretical grounds [6,10]. Recently, an explicit density functional in the stronginteraction limit was suggested. It can be constructed from so-called co-motion functions and captures fully non-local effects in the strong-interaction limit [44].
An example of a successful resummation is the particle-arXiv:1604.03929v1 [cond-mat.mtrl-sci] 13 Apr 2016 hole random-phase approximation (RPA), in which an infinite number of "ring diagrams" is summed [8,[45][46][47][48][49][50]. This has been widely recognized as key to make RPA applicable to small-gap or metallic systems. The RPA method provides the correct H 2 dissociation limit [2,15]. Unfortunately, it suffers from a heavy "self-correlation" error for one-electron systems [1][2][3][4][5], and thus yields an even worse H + 2 dissociation behaviour than conventional density functionals. In addition, RPA exhibits an incorrect repulsive "bump" at intermediate H 2 bond distances. These deficiencies have in the past been attributed to a lack of self-consistency in RPA [2,5], which was disproved by actual self-consistent calculations [15,51]. Another widely accepted hypothesis attributes these deficiencies to the lack of higher order diagrams and spurred considerable beyond-RPA developments in the past few years [3,16,18,19,[52][53][54][55][56][57]. Other interesting developments in this realm include reduced density matrix theory [58] or self-consistent Green's function frameworks [15,59]. A successful beyond-RPA method is renormalized second-order perturbation theory (rPT2), which adds an infinite summation of the second-order exchange diagram of PT2 (termed second order screened exchange (SOSEX)) [4,[60][61][62] and renormalized singleexcitation (rSE) [52,63] diagrams on top of RPA [16]. rPT2 does not diverge for small-gap systems and is free of one-electron "self-correlation". It thus provides the correct description of one-electron systems including H + 2 and individual H atoms, but fails for H 2 dissociation if breaking spin symmetry is not allowed [Refs. 3, 16, and 53 and also below]. Further improvements have been stipulated in the context of the couple-cluster (CC) theory [18,[53][54][55][56][57] or the Bethe-Salpeter equation [3,19,22]. These methods, although proposed from different perspectives, can all be interpreted as attempts to explicitly introduce more "selective summations to infinite order" in the density functional perturbation framework. Even though these methods improve over standard RPA schemes with varying degrees of success for the H + 2 /H 2 dissociation problem, no improvement to date removes the "bump" while simultaneously yielding the correct dissociation limit for H + 2 and H 2 , indicating the difficulty of understanding this problem in any perturbative framework.
In this paper, we lay the ground for an efficient orbitaldependent correlation functional based on the Bethe-Goldstone equation (BGE) [64], which is derived from the correlation of two particles [22]. As the BGE is the simplest approximation which provides the exact solution for one-and two-electron systems, it is a good starting point to understand the aforementioned H + 2 /H 2 dissociation challenge. In Sec. II, we formulate the BGE in the context of DFT through the adiabatic connection approach. In contrast to the normal resummation strategy in density functional perturbation theory, we propose a new correlation functional by terminating the BGE expansion at the second order (BGE2). As shown in Sec. III, this BGE2 approximation gives a good description of both H + 2 and H 2 dissociations, without requiring any higher order connected Goldstone diagrams that are commonly believed to be necessary. We further analyse BGE2 analytically in the minimal basis H 2 model. We show that BGE2 is size-extensive and free of one-electron "self-correlation".

II. BETHE-GOLDSTONE EQUATION IN DENSITY FUNCTIONAL THEORY
In the adiabatic-connection (AC) approach of density functional theory (DFT) [20,25,26], the non-interacting KS system is connected to the physical system by an adiabatic path. The density n along the path is fixed to the exact ground-state density. The Hamiltonian for a family of partially interacting N -electron systems in this path is controlled by a coupling-constant parameter λ, (atomic units are used hereafter unless stated otherwise): Here,Ĥ s is the Hamiltonian of the non-interacting KS systemĤ where v s (r) is a multiplicative one-electron potential comprising the external potential (v ext ) arising from the Coulomb interaction between the electrons and the nuclei, the Hartree potential (v H ), and the exchange (v x ) and correlation (v c ) potential. The operatorv λ is also multiplicative and constrained to satisfyv 0 = 0 and v 1 =v H +v xc . Thus,Ĥ λ=0 =Ĥ s , whileĤ 1 is the Hamiltonian of the fully interacting system. From the perspective of many-body perturbation theory, ∆ λ =V ee −v λ /λ is a perturbation of the non-interaction KS Hamiltonian, which does not change the ground-state density n. By using coordinate scaling [20,26,65], it was shown that where n α (x, y, z) = α 3 n(αx, αy, αz), and v c (r i , α) is the correlation potential of the scaled correlation energy with respect to the normal density n.
In the AC framework [20,25,26], the xc functional can be interpreted as the coupling-constant integration, where E λ xc [n] is the xc functional for a given couplingconstant λ [66] V λ xc is the corresponding xc potential for a given coupling constant λ. We can further define the exchange E λ x [n] and correlation E λ c [n] components separately Here Ψ λ n is the ground-state wave-function on the AC path with the coupling constant λ, which gives the same ground-state density n as the physical system (λ = 1). Ψ 0 n = Φ n is thus the ground-state wave-function of the non-interacting KS system. E H [n] is the Hartree energy which is an explicit functional of the density. Immediately, we have V λ , as the exchange density functional E λ x [n] defined in this manner is linear in the coupling constant λ. And the corresponding Hartree potential v H is written as with the definition of the two-electron four-center integral as In contrast, it is in general not possible to obtain the exact E λ xc [n] for any λ = 0, since the electron-electron repulsion operatorV ee appears explicitly in the Hamiltonian, and the ground-state wave-function Ψ λ n cannot be obtained exactly. This is also true for two-electron systems, although the ground-state wave-function is now just a simple electron-pair function Ψ λ n = Ψ ab .
As one of the motivations in this paper is to construct a functional which can provide an accurate description for both H 2 and H + 2 dissociations, we start from the Bethe-Goldstone equation (BGE) ofĤ λ [22], which is derived from the correlation of two particles, and is thus the exact solution for one-and two-electron systems. The BGE explicitly solves the Schrödinger equation for each electron pair ab interacting through a perturbationĤ 1 (λ) Here, we consider the electron-electron interactionV ee of electron pair ab explicitly, while leaving the interaction with the other N -2 electrons on the mean field levelv MF ab . For two electrons we trivially havev MF ab = 0. However, for more than two electrons we would have to make this approximation explicitly. The resulting perturbation operator isĤ (13) with the definition ofv MF ab as where {φ i } are the KS orbitals and The KS orbitals can be used to generate an antisymmetric non-interacting KS electronpair function Here the numbers 1 and 2 are a short-hand notation for the tuple of space and spin variables of the first and second electron, respectively. The corresponding noninteracting Green's function G 0 (1, 2; 1 , 2 ; E ab ) for this electron pair ab is where { i } are the KS eigenvalues AsĤ s (Eq. 2) is a one-electron operator, it is also possible to reorganize the eigenvalues in terms of each electron pair ab which could be considered as the zero-order approximation of the electron pair energy E ab . Now we introduce the first approximation to the BGE. We neglect the single excitation contribution, i.e., the first sum in Eq. 16.
As will be discussed in Sec. IV A, this approximation, together with the other two approximations we will make later, is essential for achieving an efficient correlation functional which, however, keeps the exact solution for the H 2 dissociation limit in the minimal basis.
With this approximation, the electron-pair function Ψ ab can be written as It is convenient to introduce intermediate normalization Φ ab |Ψ ab = 1. Together with the expression of the expectation value of the perturbation energy the BGE electron pair function Ψ ab (1, 2) becomes Since the BGE electron pair function Ψ ab appears on both sides of this equation, both Ψ ab and E ab (eq. 23) contain an infinite sequence of Goldstone diagrams [22,23], as one can easily see by inserting Eq. 22 into Eq. 21 We now consider the different terms step by step. First, we expand the e 1st ab (λ) term on the right-hand side It is the first-order correction to the non-interaction electron pair energy ab defined in Eq. 18. Utilizing the definitions ofv λ (Eq. 4) andv MF ab (Eq. 14) for the fullyinteracting system (λ = 1), we have where 1st a (1) is the corresponding first-order correction of the non-interaction electron energy a for the fullyinteracting system. Andv HF x is the Hartree-Fock like exact exchange operator defined as Together with the the non-interaction electron pair energy ab , this leads to the electron-pair total energy at the first-order many-body perturbation level, which contains only the exact exchange. Next we define the BGE electron-pair correlation energy e BGE ab (λ) as For two electrons, e BGE ab (λ) is the total correlation energy E BGE c [n](λ). For more than two electrons we have to sum up the the correlation energies of all electron pairs: Finally, the BGE total energy for the fully-interacting system (λ = 1) becomes where E EX tot is the exact-exchange total energy in the KS-DFT framework. BGE is thus exact for one-and twoelectron systems, but approximate for more electrons, because interaction terms between three or more electrons are missing.
With the definition of e BGE ab (λ), eq. 23 becomes The second term on the third line of eq. 30 emerges when we replace Ψ ab by eq. 22. This expansion reveals that the BGE correlation energy contains an infinite summation of a sequence of Goldstone pair diagrams [22,23]. Therefore, this sequence of Goldstone diagrams contains only two hole lines, representing the electron pair ab, and the infinite summation goes through all the ladder diagrams over two particle lines (see fig. 1). In other words, the intermediate pairs always propagate as electrons [22]. Conversely, eq. 30 has to be solved self-consistently, as the electron-pair energy E ab depends on itself.
As alluded to before, the BGE accounts for the correlation of two electrons and thus provides the exact solution for one-and two-electron systems. However, eq. 30 is not exact, because we omitted the single excitations at the very beginning. Nonetheless, our approximation should still be able to capture the subtle (near)-degeneracy static correlation effects at H 2 dissociation and eliminate the one-electron "self-correlation" error as does the configuration-interaction method with FIG. 1. The Goldstone diagrams in the BGE are an infinite sequence of particle-particle ladder diagrams (pp-ladder) [22]. The squiggly lines refer to the bare Coulomb interaction.
double excitations (CID) or the coupled-cluster doubles (CCD) method. Nesbet [67] has demonstrated that BGE is equivalent to the so-called independent electron-pair approximation (IEPA) in quantum chemistry [23], which can be considered as an intermediate approximation between CID and MP2. However, we prefer to keep the BGE acronym (eq. 30) as it is more compact and easily linked to an expansion of Goldstone diagrams [23,68,69] which is helpful for further discussions ( fig. 1).

A. Derivation of BGE2
The BGE electron-pair correlation energy e BGE ab in eq. 30 contains an infinite sequence of particle-particle ladder diagrams (pp-ladder resummation), as shown in fig. 1, because the BGE wave function Ψ ab also appears on the right-hand side of the equation. In addition, eq. 30 should be solved iteratively as e BGE ab appears on both sides of the equation (e ab -coupling effect). These two mechanisms cooperate to deliver an accurate description of exchange and correlation in one-and two-electron systems. It has been argued that an explicit (or implicit) resummation of a selected series of diagrams (e.g., the pp-ladder resummation shown in fig. 1) is necessary to remove the divergence at degeneracies of any finite-order perturbation theory [7,8,16,19,48]. However, in this work we will show that the same effect can be achieved by the e ab -coupling effect at finite orders of perturbation theory. This allows us to terminate the BGE expansion at the second order, as long as we retain the e ab -coupling effect. This second-order BGE (BGE2) is the second approximation we make: Here we have utilized the fact that bothv MF ab andv λ are one-electron operators, which do not contribute to the expectation value between the ground state and a double excitation.
We will show in Sec. IV A that BGE2 only dissociates H 2 in a minimal basis correctly, if we remove e 1st ab from the denominator. So in the following we drop e 1st ab . This is our final approximation: e BGE2 ab (λ) now appears as a simple sum-over-state formula that is similar to standard PT2 and thus exhibits the same computational scaling. We will again need to sum all electron pairs to obtain the full BGE2 correlation energy E BGE2 c (λ) for systems with more than two electrons A distinct advantage of eq. 32 is that the dependence on the coupling constant λ is now simple and well-defined. Following eq. 5 we can easily obtain the BGE2 electron- (34) and then for a many electron system The BGE2 correlation potential also has to be solved iteratively, which prevents us from making further analytical manipulations. However, for H 2 in a minimal basis, the BGE2 correlation energy and potential can be solved analytically, which will give us more insight into BGE2. This will be discussed later in Sec. IV A. At this point, we will recap the approximations made in the derivation of the BGE2 correlation functional: 1. From the outset we chose a pair theory. The full BGE (eqs. 11 and 12) explicitly treats interactions in one electron pair and is exact for one-and twoelectron systems. For more than two electrons, the interaction from other electrons can be taken into account in a mean field fashion (eq. 14). Then the correlation energy sums up the correlations of all electron pairs (eqs. 28 and 33).
2. The BGE can be solved by means of Green's functions (eq. 20). Here we omit the single excitation contribution in the construction of the noninteracting Green's function G 0 (eq. 19). We argue that this approximation is justified and captures the subtle exchange-correlation effects in one-and two-electron systems as does CID or CCD.
3. The next approximation is to terminate the BGE expansion at the second order (eq. 31). This implies that we remove the infinite summation of particleparticle ladder diagrams ( fig. 1). The resulting BGE2 approximation still retains the e ab -coupling effect. We will show later that the e ab -coupling at second-order in perturbation theory is sufficient to capture correlations that emerge from (near)degeneracies and that higher-order connected Goldstone diagrams are then not needed. However, if we are far from degeneracies (e.g., H 2 in the middle of dissociation) BGE2 alone is not sufficient and higher order diagrams would be required.
4. The final approximation (eq. 32) removes the firstorder perturbation term e 1st ab from the denominator of eq. 31. On the one hand, this omission removes the difficulty of having to consider the unknown density adaptive operatorv λ explicitly along the AC path (see eqs. 4 and 24). On the other hand, we will show in Sec. IV A that together with the other two approximations, this approximation is necessary to deliver an accurate description of H 2 dissociation in a minimal basis. Comparing to the standard PT2 expression, the only difference in the BGE2 correlation expression (eqs. 32 and 33) is that the BGE2 electron-pair correlation e BGE2 ab (λ) itself appears in the denominator. e BGE2 ab (λ) should acquire a finite value to prevent the numerical divergence for small-gap systems where ∆ rs ab → 0. To distinguish BGE2 from PT2, modified particle and hole lines (double line) are introduced in fig. 2 to represent the e ab -coupling effect in the BGE2 method which corrects the double-excitation energies ∆ rs ab and should be solved iteratively. We will demonstrate the accuracy of the e ab -coupling effect later both numerically (Sec. III C) and analytically (Sec. IV A). In this section, we will provide a many-body perturbation theory perspective of the e ab -coupling effect.
In quantum chemistry, the configuration interaction equation with singles and doubles (CISD) is usually solved with iterative techniques [70], to avoid a direct diagonalization of the large Hamiltonian matrices in configuration space. Pople et. al. [71] demonstrated that such iterative algorithms are more than just a technical trick. From a many-body perturbation theory viewpoint, each iteration introduces higher order terms. For example, after the second iteration, the second-and third-order terms emerge in the CISD energy expression but with a scaled weight [71,72]. This would also be true for the e abcoupling effect in the complete BGE expansion (eq. 30). But can we write down a perturbative expansion for the second-order BGE expression (eq. 32 and fig. 2)? In other words, does the e ab -coupling effect introduce higher order perturbation terms during the iterative procedure? In this section, we will answer these questions step-bystep.
To demonstrate the behavior of BGE2 for (near)degeneracies, we analyse its limit as ∆ rs ab goes to zero. To do so we need to introduce a level-shift (L) into the expression of the BGE2 correlation energy (eq. 32) where To be able to expand eq. 36 into a geometric series we require This leads to the following constraint for L By definition we have e BGE2 ab ≤ 0 and ∆ rs ab ≥ 0. In addition, ∆ rs ab > |e BGE2 ab (λ)| holds for insulators, most semiconductors and even for most of the double excitations in small-gap systems (excluding cases where ab refers to the highest occupied molecular orbital (HOMO)). We can then always choose L=0.
However, we cannot take L=0 when the energy gap of double excitations from the HOMO to the LUMO (∆ rs ab ) tends to zero and ∆ rs ab ≤ |e BGE2 ab (λ)|. Then only a positive level shift L > − 1 2 (∆ rs ab + e BGE2 ab (λ)) will guarantee a convergent geometric expansion. We will discuss the positive level shift later and first analyse the L=0 case. Expanding eq. 36 into a geometric series yields which for L = 0 becomes This expression allows us to analyze the BGE2 correlation energy from a many-body perspective by iterating the right-hand side. Here, we examine the simplest two terms. The first term in the BGE2 expansion (n = 0) is nothing but the standard PT2 correlation energy e BGE2,1st The Goldstone diagrams of the PT2 correlation are shown in fig. 2). The second term becomes e BGE2,2nd where S is the normalization of the first-order perturbative wave-function Φ 1st ab of the electron-pair ab ab is a fourth-order perturbation in terms of the coupling constant λ. This expansion includes only even powers of the perturbation. On the other hand, e BGE2,2nd ab can be interpreted as 32 quadruple-excitation Goldstone diagrams which, however, are both disconnected. The e ab -coupling effect in BGE2 does therefore not produce higher-order connected Goldstone diagrams. Using the so-called Hugenholtz diagram rule [73], these 32 quadruple-excitations can be represented by two Hugenholtz diagrams, which are shown in fig. 3. Recently, a self-consistent Green's function schem was proposed at 2nd arder as well [59]. This self-consistent second-order self-energy method also exhibits promising performance for systems with strong correlation. It would be very interesting to compare the diagrams of the Green's function theory with BGE2 in the future. Now we consider the case of ∆ rs ab ≤ |e BGE2 ab (λ)| where the (near)-degeneracy effects are dominant. As mentioned above, to guarantee a convergent geometric expansion, a positive level-shift L > − 1 2 (∆ rs ab + e BGE2 ab (λ)) is required. By inserting the binomial series e BGE2 We note that the L=0 case in eq. 41 is a special case of eq. 45. For positive level shifts, it can be easily proven that the first two expansion terms of eq. 45 are the same as in eqs. 42 and 43, only that the level-shift L appears in the denominator. We plot the perturbative expansion of the BGE2 correlation in fig. 4. For the L=0 case, the BGE2 correlation can be interpreted based on the standard perturbative expansion. However, if ∆ rs ab ≤ |e BGE2 ab (λ)|, a positive level shift is required to guarantee a well-defined perturbative expansion of the BGE2 correlation. In fig. 4, we show the firstand second-order geometric expansion of the direct term of the BGE2 correlation. The exchange part can be expanded in the same way. Thick lines are introduced to represent a positive level-shift L to the double excitation energies ∆ rs ab . For both cases, higher order excitations are not involved, but only rescale weights of the existing contributions. The S is the modified S with a positive level-shift L: The divergence of Møller-Plesset (MP) and Görling-Levy (GL) perturbation theories has been widely discussed in quantum chemistry [74][75][76][77][78]. Small-gap systems with strong (near)-degeneracy effects are, of course, one kind of failure, as the perturbation energy diverges at any finite order. For non-degenerate systems where ∆ rs ab is not exactly zero, the MP (or GL) perturbation expansion does not diverge. However, Leininger et al. found that the perturbation expansion also does not converge toward the exact solution even for very simple systems such as Ne, F − , and Cl − [75]. The individual terms in the perturbative expansion of these systems do not diverge, but exhibit an oscillatory divergence with increasing order [74,75]. From a mathematical point of view the MP expansion fails, if the single reference, e.g., HF or KS, is far from the exact ground state [76][77][78]. However, there is no simple diagnostic tool to determine when a multi-reference problem breaks the MP expansion.
In this work we propose to use the condition ∆ rs ab ≤ |e BGE2 ab (λ)| as a simple criterion to judge the divergence of a single-reference perturbation method. If the HOMO-LUMO gap of a given single reference is smaller than the absolute value of the corresponding BGE2 electron-pair correlation energy, it is not advisable to use a perturbative method based on this single reference, because the BGE2 correlation cannot be expanded directly without a proper level-shift L (eq. 36). The value of the level shift L then quantifies the multi-reference nature of each electron pair in the investigated systems.

C. H2 and H + 2 dissociation
Our BGE2 xc functional encompasses the exact exchange and the BGE2 correlation term (eqs. 32 and 33) It has been implemented in the Fritz Haber Institute ab initio molecular simulations (FHI-aims) code package [79,80]. Due to its simple sum-over-state formula, BGE2 has the same computational scaling as standard PT2 in terms of both time and memory. Although the e abcoupling requires an iterative solution, convergence is fast in our experience, and an accuracy of 10 −8 Hartree can be achieved within a few iterations.
In fig. 5 we plot the H 2 and H + 2 dissociation curves for various methods (BGE2, PBE, PBE0, PT2, RPA, and rPT2). All results are obtained with input KS orbitals from a PBE0 calculation [81][82][83]. All calculations, including the CISD reference, are carried out with FHIaims using the NAO-VCC-5Z basis set [80,84]. In fig. 6 = + ×S + · · ·, if ∆ rs ab > |e BGE2 ab | and L = 0 H2 (A) and H + 2 (B) dissociation curves without breaking spin symmetry. All calculations, including the configuration interaction method with singles and doubles (CISD), have been carried out with FHI-aims [79] and the NAO-VCC-5Z basis set [80]. For one and two electron systems, CISD provides the exact curves, which are thus denoted as CI. The HF, PBE and PBE0 results are obtained self-consistently. The HF orbitals are employed to evaluate the CISD results, and the PT2, RPA, rPT2, and BGE2 calculations are on top of a PBE0 reference. In the smaller panel for H2 dissociation, the total energies of two isolated Hydrogen atoms are plotted for each method. And for H + 2 dissociation, the smaller panel shows the total energies of one isolated Hydrogen atom.
we show the same curves for RPA, rPT2 and BGE2 for different starting points (PBE, PBE0 and HF).
In the following we will analyse the performance of the different approaches shown in fig. 5 class by class. In the dissociation of H + 2 only non-local exact exchange is required and correlation is absent, while dissociated H 2 contains strong (near)-degeneracy static correlation, which current DFT methods typically underestimate. As illustrated in fig. 5, the PBE0 functional fails in both cases, yielding a heavy one-electron "selfcorrelation" error (around 66 mHartree in the H + 2 dissociation limit) and a significant underestimation of the (near)-degeneracy static correlation limit (around 119 The H2 dissociations calculated by RPA (A), rPT2 (B) and BGE2 (C). The nomenclature adopted here: F@SC is the advanced functionals (F), i.e. RPA, rPT2, and BGE2, respectively, evaluated with the orbitals of different self-consistent (SC) schemes, i.e. HF, PBE0, and PBE. mHartree in the H 2 dissociation limit).
Fifth-level correlation functionals, e.g., PT2, RPA and rPT2, are fully non-local. Fig. 5 reveals that PT2 is exact for one-electron systems such as H + 2 dissociation, in agreement with previous investigations [1][2][3][4][5], but completely fails for heavily stretched H 2 , where static correlation becomes dominant. RPA exhibits the opposite behavior. It is accurate in the H 2 dissociation limit, but its severe one-electron "self-interaction" or delocalization error [85,86] affects H + 2 dissociation adversely. Adding the SOSEX term to RPA removes the one-electron selfinteraction error again such that rPT2 dissociates H + 2 correctly, but simultaneously the performance for H 2 deteriorates [15]. Henderson et al. have ascribed this behavior to a removal of static correlation by the SOSEX term [4]. Fig. 5 reveals that the BGE2 functional is free of oneelectron "self-correlation" and thus dissociates H + 2 correctly. It also delivers the correct H 2 dissociation limit and reduces the incorrect repulsive "bump" that RPA exhibits at intermediate bond lengths. We attribute this consistent improvement to three factors. 1) The oneelectron "self-correlation" error is removed at the PT2 level. 2) (Near)-degeneracy static correlation (for two electrons) is incorporated by the e ab -coupling mechanism at the second-order perturbation level without having to invoke higher-order Goldstone diagrams. 3) Due to the systematic nature of the approximations we made, we can trace the repulsive "bump" back to the 2nd order approximation, when we went from the full BGE to BGE2. In other words higher-order pp-ladder diagrams will alleviate the "bump" [19].
In fig. 6 we illustrate the starting-point dependence of RPA, rPT2 and BGE2 by evaluating all three approaches for a PBE, a PBE0 and a HF reference. The startingpoint dependence in RPA is quite pronounced, which is akin to the much investigated starting-point dependence in the corresponding GW approach for excitation spectra [87][88][89][90]. Judging by fig. 6, the starting-point dependence of rPT2 and BGE2 appears to be as pronounced as for RPA. Although it remains to be seen in the future, if this statement can be generalized.

IV. PROMISING PROPERTIES OF THE BGE2 CORRELATION FUNCTIONAL
In wave-function theory, there is a systematic way to improve the theoretical approach for the electronic correlation energy [23,67,91,92]. However, it is challenging to design methods that fulfil a certain number of exact conditions and constraints [92]. One requirement is that the solution for one-and two-electron systems such as H 2 /H + 2 is exact. Another requirement is size consistency [23,27,28], i.e., the ground-state total energy itself is an extensive quality, which should be asymptotically proportional to the system size [23]. These two criteria are widely used to judge the universal applicability of a given theoretical method from small isolated molecules to extensive solids. We will analyse analytically, how well BGE2 fares for H 2 in a minimal basis.

A. H2 in minimal basis
The BGE2 approximation is not exact for two-electron systems. However, in Sec. III C, we reported a significant improvement of BGE2 over PT2 and RPA for H 2 dissociation. In this section, we analyze the BGE2 correlation functional (eq. 32) and its correlation potential (eq. 34) for H 2 in a minimal basis. We demonstrate that the BGE2 approximation captures essential features of the adiabatic connection path that current state-of-the-art approximations do not.
The minimal basis of H 2 consists of one bonding (ψ + ) and one anti-bonding (ψ − ) KS orbital determined by the where R 1 and R 2 are the two atomic positions, φ 1s (r − R 1/2 ) is the normalized 1s atomic orbital located at the each hydrogen atom, and S R1,R2 is the overlap integral of the two atomic orbitals. In this minimal basis representation, the BGE2 correlation energy (eq. 32) and its potential with respect to the coupling constant λ (eq. 34) can be derived analytically where A = | ψ + ψ + |ψ − ψ − | 2 and B = ∆ −− ++ . Now we see that several important features of the AC path are captured by the BGE2 approximation: i) When λ → 0, we have which shows that the initial slope of the BGE2 correlation potential is twice the energy in PT2. Comparing to the exact condition, i.e. V c (0) = 2E GL2 c , the singleexcitation contribution is ignored during the first approximation in this work (eqs. 16 and 19).
ii) In the H 2 dissociation limit (|R 1 − R 2 | → ∞ and B = ∆ −− ++ → 0) we find for the initial slope V BGE2 The BGE2 correlation potential itself tends to a constant value that is independent of λ thanks to the e ab -coupling effect: The coupling-constant integration (eq. 5) is trivial to carry out and the correlation energy E BGE2 c has the same value as V BGE2 c . This is the exact correlation energy in the minimal basis [23]. In conjunction with the exact exchange energy, it completely cancels out the undesired error originating from the Hartree approximation, and thus guarantees the correct dissociation limit.
If we do not make the fourth approximation, i.e. to remove the first-order perturbation term e 1st ab in the denominator of eq. 31, the parameter B now equals ∆ −− ++ + e 1st ++ (λ). By using the definitions ofv λ and e 1st ab (eqs. 4 and 24) we have e 1st ++ (λ) in the H 2 dissociation limit Here, for simplicity, we neglect the exchange v x (r i ) and the scaled correlation potential v c (r i , α) (eq. 4), since they are small compared to the Hartree energy E H . The resulting correlation energy E c = − 2 √ 5 φ 1s φ 1s |φ 1s φ 1s recovers only 90% of the exact value in the minimal basis (eq. 51). This motivates a posteriori why we made the second approximation, i.e. omitted the single excitation contribution. In practice this single excitation contribution is non-zero in a DFT framework. However, our minimal basis consideration shows that if we were to include it, H 2 would no longer dissociate correctly, unless we would also include higher-order particle-particle ladder diagrams (the third approximation), which would significantly increase the computational cost. The BGE2 correlation functional proposed in this work is thus the simplest approximation that provides the exact H 2 dissociation in the minimal basis. Any further improvement over BGE2 has to simultaneously deal with all approximations in a proper and balanced way.
v) The λ-dependent AC path is closely connected to the behavior under uniform density scaling [94], for which the low-density limit is related to the strong correlation limit λ → ∞. The exact correlation functional should reach a finite value in the strong correlation limit [6,93,94], which is satisfied by the BGE2 model according to eq. 49 (V c → − √ A when λ → ∞). Figure 7 presents the correlation potentials of RPA and BGE2 for H 2 dissociation in the minimal basis. At the equilibrium geometry (R=0.8Å) where the energy gap is large, the BGE2 model exhibits a quasi-linear behavior. For stretched geometries (R= 3.0 and 6.0 Å), the e abcoupling then bends the correlation potential, whereas RPA overestimates the initial slopes. The BGE2 correlation potentials are similar to those of configuration interaction calculations with a quadrupole-ζ basis set [93].
We also include the common [1/1]-Padé AC model [6,10] V Padé The two parameters are determined by fixing the initial slope to 2E PT2 c and the correlation energy to E BGE2 c . As illustrated in fig. 7, the [1/1]-Padé formula describes the AC path at the equilibrium geometry well, but exhibits a tendency to underestimate the curvature and overestimate the correlation potential at stretched geometries. A better agreement can be expected by using a more sophisticated [2/2]-Padé formula [6,10].

B. Size consistency for the ground-state energy calculation
Recently, the application of advanced correlation methods, e.g., MP2 [95][96][97], RPA [61,63,97,98], and coupled-cluster theories [99][100][101], in materials science has attracted increased attention. In this paper, we adopt the K-dependence criterion [28,[102][103][104][105] to demonstrate the size consistency of the BGE2 correlation functional and its applicability to complex extended materials. The advantage and usage of the K-dependence criterion has been discussed comprehensively in Ref. 28. In short, the number (K) of k-points in the Brillouin zone is a direct measure of system size in periodic boundary conditions. Then a size-consistent method must have an asymptotic K 1 dependence.
Before we turn to BGE2, we first discuss Brillouin-Wigner second-order perturbation theory (BW2) [103,106]. The size consistency of BW2 has been disproved in Ref. 28. This helps us to better understand the size consistency of BGE2 which shares a very similar sumover-state formula as BW2.
In periodic boundary conditions, the BW2 correlation is given by Here, φ iki is a canonical HF or KS spin-orbital in the ith band with wave vector k i . Due to momentum conservation, the summation only goes over three wave vectors (k b , k r , k s ), giving rise to a factor K 3 . It can be proven that | φ aka φ bk b ||φ rkr φ sks | 2 exhibits an asymptotic K −2 dependence. E EX x and ∆ rkrsks akabk b scale as K 1 and K 0 , respectively [28]. If the denominator scales the same as E EX x (the first line in eq. 56), the overall scaling of E BW2 c becomes K 0 , which would not be size consistent. Removing E EX x from the denominator (second line in eq. 56) changes the scaling of E BW2 c to K 1/2 , which is still not size consistent. These non-physical size dependences make the second-order energy per unit cell, E BW2 c /K, go to zero as K → ∞. It has been argued that the presence of an extensive quantity, E BW2 c + E EX x , in the denominator is responsible for the lack of size consistency [28,106].
The BGE2 correlation energy (eqs. 32 and 33) in periodic boundary conditions takes the form e BGE2 ab (λ) = r<s krks Compared to the BW2 correlation energy, the only difference in BGE2 is the appearance of the correlation coupling for each electron pair ab, i.e. the e ab -coupling. Following a similar approach as for BW2, it is easy to prove that the electron-pair correlation term e BGE2 ab scales as K 0 [28]. Due to momentum conservation the summation in the BGE2 correlation energy E BGE2 c only runs over one wave vectors k b (see eq. 57). Therefore, E BGE2 c scales as K 1 and thus is size consistent.

C. Orbital invariance
It should be stated that an electron-pair approximation such as BGE2 does not have a derived wave function. It thus breaks another important feature: the orbital invariance [23], i.e. the total energies cannot be determined uniquely with respect to rotations among occupied and/or unoccupied orbitals. However, the simple sum-over-state formula of BGE2 (eqs. 32 and 57) clearly suggests that the e ab -coupling decays very quickly to standard PT2 when the energy difference becomes larger. Therefore, we expect that this orbital invariance deficiency does not affect real applications. We leave a detailed examination of this issue to the future work.

V. CONCLUSION
In this work, we present a new insight into the H 2 /H + 2 challenge in DFT. We establish the Bethe-Goldstone equation in the context of DFT through the adiabaticconnection approach. BGE is the simplest approximation to provide the exact solution for one-and two-electron systems. We propose a simple orbital-dependent correlation functional, BGE2, by terminating the BGE expansion at the second order, but reversing the e ab -coupling effect in BGE. BGE2 has a similar sum-over-state formula as the standard PT2, thus sharing the same computational scaling as PT2 in terms of both time and memory. We demonstrates that the e ab -coupling iteration procedure at the second-order expansion does not invoke higher-order connected Goldstone diagrams, but partially captures the multicenter character of each electron pair, especially in heavily stretched H 2 . A remarkable improvement of BGE2 over PT2 and RPA in the H 2 /H + 2 challenge can be observed, which suggests that the one-electron "self-correlation" and the two-electron (near)-degeneracy static correlation can be included simultaneously well at the second-order perturbation level in conjunction with a proper treatment of the multireference contributions of each electron pair. In addition of the size consistency, the advantage of the BGE2 correlation functional has been further demonstrated using H 2 in minimal basis.
However, for systems with large energy gaps or more electrons, BGE2 reduces to standard PT2 since the effect of the e ab -coupling is nearly damped out. Further development on top of BGE2 could proceed as follows: 1) From the semi-empirical double hybrid perspective, the BGE2 correlation formula could be a promising substitute of normal PT2, which opens an opportunity to extend the double-hybrid scheme into the realm of transition metals while keeping the accuracy achieved for main group elements [31]. 2) From the AC modeling perspective, we can improve the BGE2 model by satisfying more physical constraints. For example, in the strong correlation limit λ → ∞, the BGE2 model depends asymptotically on λ −1 rather than the correct λ −1/2 [7,93]. Promisingly, the plain sum-over-state PT2-like formula makes it easy to fulfill the correct asymptotic behavior by introducing an additional term that depends on λ −3/2 in the denominator of the BGE2 formula (Eq. 31). 3) From the many-body perturbation theory perspective, it would be appealing to renormalize the BGE2 scheme in rPT2 [16], which would be an alternative to construct advanced orbital-dependent functionals systematically.