Approximate energy functionals for one-body reduced density matrix functional theory from many-body perturbation theory

We develop a systematic approach to construct energy functionals of the one-particle reduced density matrix (1RDM) for equilibrium systems at finite temperature. The starting point of our formulation is the grand potential Ω[G] regarded as variational functional of the Green’s function G based on diagrammatic many-body perturbation theory and for which we consider either the Klein or Luttinger–Ward form. By restricting the input Green’s function to be one-to-one related to a set on one-particle reduced density matrices (1RDM) this functional becomes a functional of the 1RDM. To establish the one-to-one mapping we use that, at any finite temperature and for a given 1RDM γ in a finite basis, there exists a non-interacting system with a spatially non-local potential v[γ] which reproduces the given 1RDM. The corresponding set of non-interacting Green’s functions defines the variational domain of the functional Ω. In the zero temperature limit we obtain an energy functional E[γ] which by minimisation yields an approximate ground state 1RDM and energy. As an application of the formalism we use the Klein and Luttinger–Ward functionals in the GW-approximation to compute the binding curve of a model hydrogen molecule using an extended Hubbard Hamiltonian. We compare further to the case in which we evaluate the functionals on a Hartree–Fock and a Kohn–Sham Green’s function. We find that the Luttinger–Ward version of the functionals performs the best and is able to reproduce energies close to the GW energy which corresponds to the stationary point.


Introduction
One-body reduced density matrix (1RDM) functional theory provides an interesting alternative to density function theory (DFT), which holds the promise to alleviate many of the current practical failures of DFT. Especially the inherent better capability of 1RDM functional theory to deal with strongly correlated systems, e.g. the breaking of chemical bonds [1][2][3][4], Mott insulator transitions [5,6] and the fractional quantum Hall effect [7], are an encouragement to invest even further in the development of 1RDM functionals.
Currently, two main strategies are followed in the development of new 1RDM functionals. The first strategy uses the fact that the exact 1RDM functional for twoelectron systems is available which is used as a paradigm to device general N -electron functionals [1,2,[8][9][10]. The other main approach is to reconstruct the two-body reduced density matrix (2RDM) from the 1RDM guided by N -representability conditions [3,[11][12][13][14][15][16][17]. A more extensive overview of different approaches can be found in reference [18].
A disadvantage of the current approaches is that they do not provide a systematic route towards improved 1RDM functionals. In many-body perturbation theory (MBPT), however, such a route is available by including an increasing amount of terms in the perturbation expansion [19,20]. These terms in the perturbation expansion are conveniently represented by Feynman diagrams. It would therefore be advantageous to connect the MBPT framework to 1RDM functional theory, since this would allow one to use the MBPT functionals for the construction of increasingly more accurate 1RDM functionals in a systematic manner. This approach has already been heavily pursued in DFT with some considerable success [21][22][23][24][25][26][27][28][29][30], so the additional flexibility of 1RDM functional theory would lead to a further improvement of the results, especially when strong correlations play an important role. Due to the peculiar properties of 1RDM functional theory, however, the connection to MBPT is not so trivial as for DFT. In DFT one can use the Kohn-Sham (KS) system [31], from which a non-interacting Green's function can be extracted to be inserted into the MBPT functional. The construction of a similar KS system in 1RDM functional theory -a non-interacting system with the same 1RDM as the interacting system -leads to a completely degenerate orbital spectrum, since the orbitals are fractionally occupied [32][33][34]. For Coulomb systems we know at least that infinitely many orbitals are fractionally occupied [35] and there is strong evidence that they all are [36][37][38]. The corresponding Green's function of the non-interacting system in 1RDM functional theory has therefore all its poles located at the same position. Such an unphysical Green's function is bound to lead to nonsensical results when inserted directly into an MBPT functional.
The degeneracy of the orbital energies in 1RDM functional theory seems to pose a serious problem if one wants to use MBPT for the construction of functionals. However, the massive degeneracy of the orbital energies is caused by the fact that we evaluate the system at zero temperature and can be lifted by elevating the temperature to any finite value, T > 0. The degeneracy problem can therefore circumvented by performing the calculation at a finite temperature. Within a finite basis set, 1RDM functional theory allows for a completely rigorous foundation in which all required functionals are nicely differentiable and no v-representability issues arise [39]. The zero temperature limit, T → 0, can be taken at the end of the calculation. An additional advantage is that also the finite-temperature MBPT formalism is more reliable, since finite-temperature MBPT does not rely on adiabatically turning on the interaction, which leads to problems for zero-temperature MBPT in the case of level crossings [19,20]. We will elaborate later in this article more on this extraction process of the non-interacting Green's function from the 1RDM.
The extracted Green's function can now be deployed in the MBPT framework in several ways. The most straightforward option is to use the extracted Green's function as a zeroth order Green's function in the MBPT formalism. Unfortunately, this leads to complicated expressions: the effective potential needed to force the non-interacting system to have the same 1RDM as the interacting system needs to be subtracted from all occurrences of the selfenergy [40][41][42][43][44]. This means that also the energy expression depends on the effective potential itself, so becomes a severe complication. More importantly, following this approach one only obtains a reformulation of MBPT in terms of a different (less convenient) zeroth order Green's function. The results will therefore be the same as for an MBPT calculation with the same diagrams, so the only achievement would be a more complicated formalism, if one aims for fully self-consistent solutions.
A more viable approach is to use variational MBPT functionals Ω[G] for the grand potential [22,23,[45][46][47]. The total energy from these functionals are relatively stable with respect to variations in the Green's function, due their variational property δΩ/δG = 0 when G is obtained self-consistently form the Dyson equation. This has the advantage that one does not need to solve the full Dyson equation, but that it is sufficient to use a reasonable approximate Green's function to obtain a reasonably accurate value for the total energy. Such an approximate Green's function can be the KS Green's function or the non-interacting Green's function extracted from the 1RDM.
We have different variational functionals at our disposal and popular forms are the Luttinger-Ward [48] and Klein [49] functionals. Especially the Klein functional is quite popular, as it yields the simplest expression for the energy [22,23]. The Klein functional can be reformulated in terms of the exchange-correlation kernel to yield the adiabatic connection fluctuation-dissipation expression [26,27,50,51]. Though the Luttinger-Ward functional results in a more complicated expression for the energy than the Klein functional, it has the advantage of superior variational properties. We will take both functionals into consideration in this article, but limit ourselves to one of the simplest perturbative expansions: the GW approximation. The performance of both the Klein and Luttinger-Ward functionals in the GW approximation will be tested on an extensions of the two-site Hubbard model which is capable to give a good representation of the ground state of the hydrogen molecule. To achieve this, the extended two-site Hubbard model also takes the intersite interactions into account and the interaction between the electrons and the nuclei.
The paper is divided as follows. In Section 2 we present the general theory of the variational functionals that we use. In Section 3 we show how to construct a noninteracting Green's function G s [γ] that yields a prescribed 1RDM γ. In Section 4 we present the details for the molecular model that we use to test our functionals. In Section 5 we describe the various input Green's functions that we use as input for the variational functionals. In Section 6 we perform variational calculations on the molecular model and finally, in Section 7 we present our conclusions.

Variational functionals
The required theoretical framework for the use of variational MBPT functionals has basically already been presented in the framework of DFT in reference [22]. The only difference is that in the framework, we now use non-local potentials to keep the complete 1RDM fixed at all interaction strengths. The procedure is very analogous to the derivation of the Luttinger-Ward functional [20,Sect. 11.4]. For completeness, we will provide here a short overview. Consider the Hamiltonian where k ij and w ijkl are are one-and two-body matrix elements. The term v λ is a potential to keep the 1RDM identical at all interaction strengths λ to the fully interacting 1RDM, γ λ = γ 1 = γ. The potential at λ = 1 is therefore the real external potential v 1 = v ext . In the framework of KS-DFT the potential v λ would be local, i.e. diagonal in the coordinate representation or site basis, and only the density would be kept fixed at all interaction strengths, n λ = n. At λ = 0, the Hamiltonian becomes non-interactingĤ s :=Ĥ 0 and the its corresponding Green's function G s has a particularly simple form (see Sect. 3). The corresponding potential to keep the 1RDM or density of the non-interacting system identical to the one of the non-interacting system is denoted by v s := v 0 . At general λ, the Green's function is related to the noninteracting Green's function via the Dyson equation To relate the grand potentials of the interacting and non-interacting system, we will use the fundamental theorem of calculus. So we first work out out the derivative of the grand potential with respect to the interaction strength where Tr indicates a summation over all indices as well as a summation over the Matsubara frequencies. We now use the Φ-functional which can be constructed by summing over all irreducible self-energy diagrams closed with an additional Green's function [48]. Each of these closed diagrams is multiplied by the pre-factor 1/(2n), where n denotes the number of interaction lines The self-energy derived from the Φ functional as leads to a conserving approximation for Σ [52,53]. The derivative of Φ respect to the coupling strength is readily worked out as which allows us to rewrite the derivative of the grand potential as The last term is readily verified by working out the derivative of Tr ln 1 − G sΣ λ [20]. Integrating from the non-interacting system to the fully interacting one, we find that their grand potentials are related as whereΣ and tr only sums over matrix indices. By regarding the grand potential as a functional of the one-body Green's function we have retrieved the Luttinger-Ward (LW) functional Ω LW [G], for an arbitrary non-interacting reference Green's function G s . Note that this expression assumes that we are keeping the complete 1RDM constant for varying interaction strength. If we were only to keep the density fixed, γ should be replaced by n and the potential v ext and v s would be local. Note that apart from the allowance of an arbitrary spatially non-local potential instead of the true external potential to be able to work with more general G s , the result in (8) is identical to that obtained originally by Luttinger and Ward [48]. When integrating over the interaction strength, we have assumed that the chemical potential is also kept fixed at its interacting value. This condition fixes the constant in v λ , i.e. tr v λ , as it should be chosen such that the particle number remains constant at varying interaction strengths. This is readily guaranteed by assuring that the Luttinger-Ward functional yields correct number of particles, N = −∂Ω LW /∂µ. This is particularly convenient as the term µN can now easily be eliminated to obtain the total energy. The Luttinger-Ward functional is variational in the sense that when G solves the Dyson equation (2) with Σ = δΦ/δG, perturbations in Ω LW vanish to first order This does not guarantee that the Luttinger-Ward functional achieves its minimum at the solution of the Dyson equation, but only that it is stationary. One can construct infinitely many different expressions for the grand potential with a different functional dependence of the Green's function which both yield the correct value of the grand potential and are stationary at the solution of the Dyson equation [20]. A popular alternative variational functional is the one due to Klein [49] The Klein functional is simpler to evaluate, as there is no explicit reference to the self-energy. The price we pay for this simplification is that the Klein functional is less stable to perturbations than the Luttinger-Ward functional.
Though the first order variations vanish at the solution of the Dyson equation, the higher order terms typically have larger amplitudes. This higher sensitivity to the input Green's function of the Klein functional has been demonstrated for atoms [22] and diatomic molecules [23]. The simplicity of the Klein functional is especially apparent if we insert the non-interacting Green's function G s into the Klein functional. All the terms within Tr{·} disappear and only the Φ-functional remains as a non-trivial part. Finally, we address some practical applications of the variational functionals. First of all let us consider that the stationary point of the functional is attained for a Green's function G but that we insert a different input Green's functionG. Due to the stationarity property (10) (and similarly for the Klein functional) we have, with ∆G =G − G that where the second order derivative in the integrand on the right hand side is a tensor contracted with terms ∆G and we further imply a double integration of imaginary times. We therefore see that if we make an error ∆G in the Green's function then the error in Ω is, due to its variational property, only quadratic in this error. This is a very useful property as it ensure that we may obtain good values of the grand potential from rather simple input Green's functions, provided they are close enough in some sense to the stationary point. The size of the error also depends on the value of second derivative in the integral, which is different for the LW and Klein forms of the energy functional and in practice the variational property of the LW form is found to be superior [22,23]. Another way to use the energy functionals is by looking for the stationary point on a restricted domain. For example by inserting Green's functions from a general non-interacting system with a local potential we obtain a density functional [22,23] since the potential is in one-to-one correspondence with a density by means of the Hohenberg-Kohn (or Mermin) theorem. In this work we extend this to non-local potentials since at finite temperature and in a finite basis it can be proven that any ensemble N -representable 1RDM is v-representable by a non-local potential [39]. Using this property, we can consider the variational domain of Ω[G] to be the set of Green's functions of non-interacting systems at a finite temperature with a non-local potential. By varying over the non-local potentials we vary over a domain of 1RDMs and we are effectively using a 1RDM functional, whose stationary point yields an approximate 1RDM. This approach will be used in the present work. Note that by restricting the variational domain we will not recover the exact Green's function anymore at the stationary point of the variational equations. This would even be the case had we used the exact Φ-functional. We further remark that a different way to construct density matrix functionals on the basis of Green's functions is given in reference [44] and the proposed approximate scheme in Section IV C of reference [44] coincides with our use of the Luttinger-Ward functional. The method that we present here is close in spirit to the one presented in reference [54].

Construction of G s [γ]
In this section we outline how we can construct a noninteracting system at finite temperature that produces a given 1RDM. The Green's function G s [γ] of this system produces functionals of the 1RDM when inserted in the variational Luttinger-Ward and Klein functionals. These 1RDM functionals will be studied in greater detail in later sections. The Matsubara Green's function of a non-interacting system consisting of bosonic/fermionic particles can be shown to be of the form [19,20] where f k := (e β M k ∓ 1) −1 is the Bose/Fermi distribution function,f k := 1 ± f k and β := 1/T is the inverse temperature. The Heaviside step function is defined as The Matsubara energies are obtained by subtracting the chemical potential from the eigenvalues of the one-body Hamiltonian, M k := k − µ, where the chemical potential can be adjusted such that the system has the desired number of particles The corresponding non-interacting 1RDM is readily extracted from the Green's function as Since the 1RDM is diagonal in the eigenbasis of the onebody Hamiltonian, the Bose/Fermi distribution functions are the natural occupation numbers. Now let us limit the discussion to electrons, i.e. fermions. It is clear from the form of the Fermi function that it satisfies the strict inequalities 0 < f k < 1 at finite temperatures, so the occupation numbers are purely fractional at T > 0. In the limit T → 0, however, the Fermi function collapses to a Heaviside function and the occupation numbers become integers 0 or 1 depending on the sign of M k . For M k = 0 the value of the occupation number is undetermined at T = 0 and can be anything between zero and one, 0 ≤ f k ≤ 1, in principle. However, by considering how the M k → 0 in the zero temperature limit, the occupation numbers will have a well defined value in the T → 0 limit even if M k → 0. Suppose we have a given 1RDM and want to find the Matsubara Hamiltonian,Ĥ M :=Ĥ − µN , which yields this 1RDM. Since the eigenstates of the one-body Hamiltonian and the 1RDM coincide, it is clear that the eigenfunctions of the one-body Hamiltonian should be the natural orbitals (eigenfunctions of the 1RDM). The corresponding eigenvalues are readily obtained by inverting the Fermi distribution From this expression it is clear that if the Matsubara energies in the β → ∞ limit decay as M k ∼ 1/β that the occupation number will have a definite value 0 < f k < 1. If its decay is slower it will end up at one of the integer values depending on its sign. If the Matsubara energy decays faster than 1/β, the corresponding occupation number converges to f k = 1 /2. The behaviour of the occupation number in the T → 0 for these various decay rates are depicted in Figure 1. The collapse of the Matsubara energy spectrum for fractional occupation numbers in the T → 0 limit is demonstrated in Figure 2.

Two-orbital model for H 2
To keep the calculations as simpel as possible, we limit ourselves to a two-orbital model for the hydrogen molecule. Let us start from the full interacting where the one-electron matrix elements, contain the kinetic energy and external potential. The two-electron matrix elements, w, describe the Coulomb interaction between the electrons and are given as In this equation we also used the quantum chemical notation [55] which regards the two-electron integral as a weighted overlap between charge distribution ϕ * i (r)ϕ l (r) and ϕ * j (r )ϕ k (r ). When we work with orthonormal orbitals, we can make the tight-binding approximation to the two-electron integrals The underlying idea is that the charge distributions ϕ * i (r)ϕ l (r) integrates to zero for i = l and therefore leads to a much smaller value of the two-electron integral than the diagonal term i = l. If we further only use one orbital per hydrogen atom, the Hamiltonian simplifies in the tight-binding approximation tô wheren iσ :=ĉ † iσĉiσ ,n i :=n i↑ +n i↓ andN =n 1 +n 2 . This Hamiltonian can be regarded as an extension of the Hubbard model for two sites. The normal Hubbard Hamiltonian only includes the hopping term between the two sites, whose strength is governed by t, and the on-site interaction, whose strength is set by U . We additionally include a term depending on the number of particles with strength α. This term does not have an influence the eigenstates, since they are eigenfunctions of the total number of particles, but does change the corresponding eigenvalues, however, so is important to recover the correct electronic energy. Apart from the on-site interaction, we also include the interaction between the electrons if they reside on different sites and its strength is determined by the parameter w.
The last term in the extended Hubbard Hamiltonian is a self-interaction term and does not contribute in principle. However, when we construct energy functionals from the perturbation expansion, the expansion for the Φ functional cannot always directly be associated to a proper expansion of an anti-symmetric 2-body Green's function. In that case the energy functional is not guaranteed to be self-interaction free andŪ typically does make a contribution. The GW (RPA) approximation considered in this article is an example of such an expansion which is not self-interaction free, hence the final result depends on U . In this work, we will consider two different values for U . The valueŪ = U corresponds to a spin-independent interaction, such as the non-relativistic Coulomb interaction used in molecular Hamiltonians. The other sensible value to use isŪ = 0, which explicitly eliminates the selfinteraction at the level of the Hamiltonian. This is the usual choice in the Hubbard model.
To determine the values of the extended Hubbard parameters, we choose normalised 1s orbitals located at each hydrogen atom to construct our basis where the normalization factor is given by N ζ = ζ 3 /π. The exponent ζ will be variationally optimized to obtain the lowest energy. In the dissociation limit R := |R A − R B | → ∞ we have that ζ = 1. A localised orthonormal basis is readily construct by Löwdin orthogonalization [56], ϕ = χS − 1 /2 , where S ij = χ i |χ j is the overlap matrix of the non-orthogonal basis. Using this Löwdin orthogonalised basis, the one-electron matrix elements in the extended Hubbard Hamiltonian can be determined as where s := χ 1 |χ 2 denotes the overlap. For more details one these expressions and explicit forms in terms of the orbital exponent ζ and the internuclear distance R, consult Appendix A.
The two-electron matrix elements can be expressed as where we need the two-electron matrix elements in the non-orthogonal 1s-basis These two-electron matrix elements for the 1s-basis are worked out explicitly in terms of ζ and R in Appendix A. Due to the limited dimension of the Hilbert space, the extended two-site Hubbard model is easy to solve exactly. In the two-particle sector we obtain the following triplet solutionŝ with the eigenvalue E = 2α + w. The singlet states can be subdivided according to their parity. There is only one ungerade singlet state, 1  eigenvalue E = 2α + U . There are two gerade singlet states which can be expressed as where the angles α ± can be determined from the equation The corresponding electronic energies of the gerade singlet states are The energy E − is the ground state energy which is minimized by optimizing the value for the orbital exponent ζ. The optimal value, ζ opt , is shown in Figure 3. The value of the exponent goes to 1 when stretching the bond. This is expected, since in the dissociation limit the system consists of two separated hydrogen atoms. The asymptotic value is approached from below, since the 1s orbital needs to become more diffuse to facilitate binding. Binding would be more efficiently achieved by allowing the 1s orbital to polarize towards the other hydrogen atom, for example by mixing in the 2p z orbital. By allowing for polarization, the reduction in the exponent would be less. The values for the total energy from our extended two-site Hubbard model are compared to the accurate values obtained for the non-relativistic hydrogen molecule in the Born-Oppenheimer approximation by Ko los and Wolniewicz [57][58][59][60]. Since the free parameter ζ is optimized for the ground state, X 1 Σ + g , it is reproduced quite accurately. The energy of the triplet b 3 Σ + u state is of less quality around the equilibrium bond length R e = 1.4 Bohr, but behaves very well upon dissociation. The other excited states are upshifted and there minima are located at too large bond lengths. Nevertheless the overall shape of the B 1 Σ + u state is reasonable and   correctly becomes degenerate with the EF 1 Σ + g state in the dissociation limit. The double-well structure of the EF 1 Σ + g is completely absent, which is no surprise, since the higher lying GK 1 Σ + g required for the necessary avoided crossing is not present in our simple two-orbital model.

Input Green's functions
In the following sections we will assess the performance of the Klein and the Luttinger-Ward (LW) functionals using the Green's functions as input: the Kohn-Sham (KS), the Hartree-Fock (HF) and 1RDM Green's function. These non-interacting Green's functions are based on a one-body Hamiltonian. Since we only use two basis functions in the model, the symmetry of the system dictates that the eigenfunctions of the one-body Hamiltonian need to be So the Matsubara non-interacting Green's function (13) are also diagonal in this symmetry adapted basis. As we restrict the optimisation to symmetry adapted solutions, we do not allow for broken symmetry states to handle strong correlation effects. The full burden is placed on the functional. We will limit the discussion to neutral H 2 , so N = 2(f g + f u ) = 2, which immediately implies that the chemical potential should be chosen such that M g + M u = 0, so the chemical potential in the non-interacting system is required to be (32)

Hartree-Fock approximation
One of the simplest perturbation expansions involves only the Hartree-Fock (HF) diagrams whereñ := f g + f u = N/2 = 1 andf := f g − f u denotes the difference in occupation of the gerade and ungerade orbital. The corresponding HF self-energy is particularly simple, since it is local in time where the HF potential for our simple model system is diagonal in the symmetry-adapted basis and the nonvanishing elements are (see Supplementary Material) The lowest energy is obtained by fully occupying the σ g orbitals,f = 1. As the HF potential differs by w between the gerade and ungrade orbitals, so the HF gap is increased by w with respect to the gap of the non-interacting system. Hence, we have As the HF potential is completely fixed by (35), the chemical potential in the HF approximation becomes Later, we will demonstrate that the correlation part of the self-energy does not introduce any additional contributions to the constant in the potential, so µ = µ HF . Hence, the effective potential in (8) becomes

The KS Green's function
The KS system has a particularly simple realisation in this two-orbital model. As the KS potential is local, its can only have diagonal elements in the site basis. As the KS should not break the symmetry of the system, it needs to be equal on both sites, so it can only be a constant shift. As the only constant contributions come from the external potential and the Hartree part, we need to set in order to have µ KS = µ. Since t < 0 for any finite R H-H , occupying the σ g orbital will always lead to the lowest KS energy. The gap in the KS system is now completely fixed by hopping matrix elements

The 1RDM Green's function
The 1RDM Green's function is defined in a similar manner as the KS Green's function, except that we do not constrain the effective potential to be local anymore. Allowing for non-local potential, gives more flexibility, but as the potential should retain the symmetry of the system, the non-local 1RDM potential remains diagonal in the symmetry adapted basis of our model system. As the trace of the potential is fixed by the particle number, the nonlocal potential can only adjust the gap in the following manner v 1RDM So in the 1RDM setting, we have one variable to optimise: the gap . At finite inverse temperature, the relation between the gap and the difference in occupation numbers is one-to-one and explicitly given bỹ In the zero-temperature limit (β → ∞), this function collapses into a step function as illustrated in Figure 5. This demonstrates explicitly that the occupation numbers can only be fractional at zero temperature, if the gap closes [32- 34,61]. In particular, we havẽ As we loose the one-to-one relation between the gap and the occupation number difference in the zero-temperature limit, we will need to optimise over the separate pieces of the β = ∞ curve: { ≤ 0,f = −1}, { = 0, −1 ≤f ≤ 1} and { ≥ 0,f = 1}. for H 2 , but all three possibilities can readily be used in more general settings. The simplest case is the use of the HF Green's function, since one only needs to do a (restricted) Hartree-Fock calculation and insert the resulting Green's function into either the Klein (11) or Luttinger-Ward functional (8) with some appropriate approximation for Φ[G]. As we will use exclusively the GW approximation for Φ in this work, we will denote these calculations by K-GW@HF and LW-GW@HF respectively. The insertion of Kohn-Sham type Green's function results in some variational freedom, as the Kohn-Sham type Green's function encompasses all non-interacting Green's functions which can be generated via a local potential, i.e. potentials diagonal in spatial representation or site basis. So in a general setting, we can use the freedom of the local potential optimise the energy expression which results from inserting the KS type Green's function into the Klein or LW functional. Within the GW approximation, we will denote such calculations as K-GW@KS and LW-GW@KS respectively. As discussed before in Section 5.2, the symmetry constraints in the limited two-orbital model (homogeneous density) effectively leaves no variational freedom in the KS potential. This is clearly an artefact of our limited setting.

Considerations for the general case
The space of non-interacting Green's functions to be searched over can be enlarged by allowing the potential to be non-local, i.e. we only require the potential to be hermitian (v † s = v s ). This class of potentials is exactly the type of potentials used in 1RDM functional theory. The insertion of these more general non-interacting Green's functions in the energy functionals within the GW approximation and subsequent optimisation, will therefore be denoted as K-GW@1RDM and LW-GW@1RDM respectively.

Variational calculations
In this section we consider the evaluation of the LW and Klein functionals for various input Green's functions.
To benchmark the results we have also solved the Dyson equation in the GW approximation self-consistently to obtain the true stationary Green's function of these functionals. The code has been implemented for spinindependent interactions, so we only have these results for the spin-independent interaction,Ū = U . For bond distances larger than 6.4 Bohr, we had difficulties converging the results. We therefore only include the results up to 6.4 Bohr, as only those are useful for comparison.

Klein functional
The Hartree and exchange part of the Φ functional were already evaluated in the previous section (Hartree-Fock), so we only need to evaluate the correlation part. The correlation part of the Φ-functional in the GW approximation can be written as (41) where the trace Tr is both over matrix entries as well as the Matsubara time/frequencies. 1 The bar over the interaction indicates that its indices are in chemical ordening,w ijkl := [il|jk] = w iklj , and the polarization bubble is defined as in frequency space. Using that tr ln(M ) = ln |M |, i.e. the logarithm of the determinant, the expression for the correlation part of the Φ-functional in the GW approximation can be further simplified to where the trace, tr, is now only over matrix entries. For our simple non-interacting Green's function, the determinant can be worked out as (see Supplementary Material) where := u − g denotes the gap and we introduced the short-hand notations u := 1 2 Ū − w and v := 1 2 U − w . Evaluating the remaining sum over the bosonic frequencies (see Supplementary Material) yields Φ GW c,s = ln sinh where Taking the T → 0 limit, we find that the Klein functional gives the following approximation for the correlation energy forñ = 1 Since the K-GW correlation energy (47) depends explicitly on the self-interaction termŪ , we should investigate different values. As mentioned in the introduction, we will investigate a spin-independent interaction (Ū = U ) and the explicitly self-interaction free model (Ū = 0). The total K-GW energy is plotted as a function of the gap/ occupation numbers for these two cases in Figure 6 for R H-H = 1.4 Bohr along the 3 segments of the β = ∞-path as depicted in Figure 5. In the case of a spin-independent interaction, the K-GW energy can be calculated for all values of the gap as is shown by the continuous line. If the self-interaction term is set to zero,Ū = 0, small values of the gap yield an imaginary part to the energy, since we have u − v = −U/2 < 0 in the second square root ζ − of the K-GW correlation energy (47). Only when | | ≥ U or = 0, we obtain a real value for the K-GW energy as is apparent from the discontinuous curve in Figure 6.
The relative positioning of the minima of the total GW energy as depicted in Figure 6 turns out to be the same at all bond distances. The minimum is always located at a fully occupied gerade orbitalf = 1 for both choices of U . In the case of a spin-independent interaction,Ū = U , the optimal gap is always located at = 0. Excluding selfinteraction term,Ū = 0, always yields the global minimum at the smallest strictly positive gap, = U . Different methods to use the K-GW energy expression (47) in the spin-independent interaction case are compared to the exact ground state in Figure 7. The 1RDM  version of the K-GW functional (K-GW@1RDM) is consistently too low, but becomes exact in the dissociation limit. Restricting the potential to be local (K-GW@KS) raises the total energy and it becomes quite accurate around the equilibrium distance, though still slightly too low. In the dissociation regime, the upshift is too high and creates the infamous RPA bump. Since the KS gap also  vanishes in the dissociation limit, the correct dissociation limit is retained [23,62]. Alternatively, one could use the HF gap as input (K-GW@HF). The HF gap is larger than the KS gap due to the additional non-local exchange in the potential (36a). This leads to a higher total energy, since the K-GW correlation energy increases monotonically in for positive gaps as can be seen from (47). The K-GW@HF energy gives an improvement over the KS gap result around the equilibrium distance, in the sense that the total energy agrees even better with the exact result (see Fig. 7). Upon dissociation however, the artificial bump is raised to even higher values, leading to a worse performance than the KS gap. Though not so clear from Figure 7, the dissociation limit remains correct, since also the HF gap goes to zero when R H-H → ∞. Albeit, the closure rate is much slower. Due to the additional exchange contribution to the gap, the HF gap only closes as 1/R H-H , whereas the KS gap closes exponentially.
Since the GW approximation is not self-interaction free (depends onŪ ), one would expect that the results would improve if the self-interaction term would be explicitly excluded from the Hamiltonian. The results for the total energy forŪ = 0 are shown in the upper panel of Figure 8. In fact, elimination of the self-interaction does actually not lead to any improvement at all. Using the K-GW functional as a 1RDM functional leads to even lower energies around the equilibrium distance and upon dissociation the energy becomes higher than the exact results. Even the dissociation limit is not correct anymore, since the gap needs to remain finite, = U , to prevent the energy from having a complex part, as mentioned before in connection with Figure 6. The results are even more disastrous when the KS or HF gaps are used. Since both gaps decrease when stretching the bond, there is always a point from where they become smaller than U and the K-GW correlation energy becomes complex. This is illustrated in the lower panel in Figure 8, where the different gaps of each calculation are plotted. Since the 1RDM gap is always the minimum gap, the failure of the K-GW functional with the KS/HF gap exactly occurs when it crosses the 1RDM gap = U . At these points, the K-GW energy from the KS/HF gap becomes exactly equal to the K-GW@1RDM functional.

Luttinger-Ward functional
In an attempt to improve the results, one can use the Luttinger-Ward (LW) functional instead of the Klein functional, as the LW functional has superior variational properties [22,23]. The superior variational properties come at the cost of a more complicated functional which requires the evaluation of the following additional terms where the modified self-energyΣ was defined in (9). Since the HF part has already been evaluated before, cf. (34) and (35), we only need the correlation part of the self-energy in the GW approximation evaluated at G s . The evaluation is somewhat simplified by the fact that the self-energy is also diagonal in the symmetry-adapted basis. Its diagonal elements can be evaluated to be where we introduced As in that case also M g = − M u , the components of the self-energy are simply related as This relation is convenient, as it implies that the number of particles is not affected by the correlation part forñ = 1 (see Appendix B), proving the validity of the constant in (trace of) the potentials used to generate the input Green's functions (38). The linear term of the LW correction in (48) can now be evaluated to be Unfortunately, the logarithmic term in (48) cannot be evaluated analytically. The integrand decays as 1/ω 2 , so numerical evaluation is not a problem. In the gapless limit, → 0, however, an analytical solution is in reach.
In the limit of vanishing gap, the self-energy simplifies considerably to In the β → ∞ limit its contribution will therefore vanish, so only the contribution of the correlation potential in the logarithmic term remains. Evaluation of the logarithmic term with only the correlation potential is rather straightforward and gives where The values for the total energies evaluated with the LW functional in the GW approximation are plotted as a function of the gap and the occupation number difference in Figure 9 for a bond distance of R H-H = 1.4 Bohr. In the case of a spin-independent interaction,Ū = U , the total energy can be calculated for any gap and the minimum is now located at a positive gap. In the case of the selfinteraction free interaction,Ū = 0, the energy is only real when the gap is either sufficiently large or zero. This is exactly the same situation as described before for the Klein functional. For short bond distances the minimum is located atf = 1 and = 0, but at larger bond distances the minimum shifts to a positive gap. This relocation of the minimum is demonstrated in Figure 10 for R H-H = 5.0 Bohr, and occurs around R H-H = 4.25 Bohr. The relocation of the minimum only occurs for self-interaction free case and not for the spin-independent interaction case.
In the top panel of Figure 11 we show the results for the total energies from the Luttinger-Ward functional for the spin-free Hamiltonian,Ū = U . Compared to the Klein functional we see a huge improvement of the total energy when the LW functional is used as a 1RDM functional, but unfortunately, the dissociation limit is not correct anymore. As a matter of fact, the optimised 1RDM results are almost identical to the results when the HF Green's function is used as input. Only in the dissociation limit the optimisation of the non-local potential yields a somewhat lower result. This trend is corroborated by comparing the HF gap and the 1RDM gap visualised in the lower panel of Figure 11. At equilibrium distance the gaps are nearly identical and they start to deviate somewhat when the interatomic distance is increased. Somewhat surprisingly, increasing the gap slightly from its HF value yields a lower energy.
We have also included the results when the KS Green's function is used as input for the LW functional. The KS Green's function now gives the worst result of all Greens functions, which is somewhat surprising, as it performed so well in the Klein functional. The results for the Klein functional are mainly due to error cancelations between the inflexibility of the KS Green's function and the bad variational properties of the Klein functional. A similar inferior behaviour of the KS Green's function as input for the LW functional has been obtained for the H 2 molecule in a large Slater basis [23]. Now let us consider the results in the self-interaction free case,Ū = 0. The results for the total energies are depicted in Figure 12. Most notable are similar catastrophes for the KS and HF input Green's functions as for the Klein functional, cf. Figure 8. When the bond is sufficiently stretched, both the KS and HF Green's function yield from some distance onwards a complex energy, due to the square root in ζ − . These catastrophes therefore occur at exactly the same bond distances as for the Klein   function the catastrophe for the LW functional occurs in a different manner than for the Klein functional. In the case of the KS functional the behaviour of the two functionals is nearly identical. When using the full flexibility of a non-local potential to generate a trial Green's function, the LW functional does not perform better than the Klein functional in theŪ = 0 case. Around the equilibrium distance the total energy is somewhat lower compared to the Klein functional, so the LW functional yields worse results. Approaching the dissociation limit, R H-H = 10.0 Bohr, the results are somewhat improved. Arguably, the most interesting feature when using the LW as a 1RDM functional is the kink at R H-H = 4.25 Bohr. This is the location where the global minimum jumps from the = 0 solution to > U solution (see Figs. 9 and 10). When one of these solutions becomes higher in energy, we have continued plotting the solution with a dotted line in Figure 12. If now the = 0 solution would be disregarded altogether, we see that the finite gap solution actually yields a quite reasonable result. We could not think of a proper physical argument to discard the gapless solution, so we do not promote this procedure. All these issues remain a topic for future research.

Conclusion
In this work we considered two different energy functionals of the 1RDM. These were obtained by restricting the domain of the LW and Klein functionals to non-interacting Green's function for a Hamiltonian with a spatially non-local potential. These functionals were evaluated within the GW approximation for two different cases of self-interaction for an extended Hubbard Hamiltonian representing a model system for the hydrogen molecule. We compared the binding curves of this system with exact results as well as with energy functionals with HF and KS input Green's functions. We found that the LW functional outperforms the Klein functional and gives results very close to the optimal selfconsistent GW Green's function which shows that there is room for future improvement of the method. Especially the GW approximation does not correctly treat the selfinteraction correctly and a probable improvement can be attained by including T-matrix and ladder diagram expansions to introduce the necessary exchange-type diagrams. The system that we considered is very rigid due to its low dimensionality and systems in larger basis sets are likely to increase the variational freedom to attain more accurate results. We further studied the system in the zero-temperature limit which led to discontinuities in the energy landscape. It is therefore worthwhile to further explore systems at finite temperature to smoothen these discontinuities.
KJHG thanks B.C.E. Giesbertz for useful suggestions and gratefully acknowledges a VENI grant by the Netherlands Foundation for Research NWO (722.012.013).

Author contribution statement
The idea came forth in a discussion between KJHG and RvL. The integrals for the model parameters have been evaluated by RvL. All calculations on the variational functionals have be done by KJHG. AMU has performed the benchmark calculation. The paper has been written by KJHG and RvL.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Appendix A: Derivation of matrix elements of the two-orbital model for H 2 Defining ρ := ζR, the one-electron integrals in the nonorthogonal 1s-basis become s := χ 1 |χ 2 = 1 + ρ + ρ 2 3 e −ρ , The unique two-electron integrals in the non-orthogonal 1s-basis are + 16 γ + ln(ρ) 9 + 18ρ + 15ρ 2 + 6ρ 3 + ρ 4 − 32 9 − 3ρ 2 + ρ 4 e 2ρ Ei(−2ρ) where γ ≈ 0.5772 is the Euler-Mascheroni constant and we used the exponential integral The formula (A.2d) has been derived using the procedure presented in [63] and agrees to all digits with the numerical results presented in [64]. In the calculations we use the Mulliken approximation to the two-electron integrals [65] (il|jk) ≈ S il S jk 4 (ii|jj) + (ll|jj) + (ii|kk) + (ll|kk) , where S ij := χ i |χ j is the overlap. That means that the following integrals are approximated as The Mulliken approximation has been verified numerically as a function of the bond length combined with the optimized exponent. The error of the Mulliken approximation is shown in Figure A.1 and is of the order of mH in around the equilibrium distance and decreases quickly for the (12|12) integral upon dissociation. When transforming to the Löwdin orthogonalised basis we obtain exactly the one-electron elements used in the extended Hubbard model (24). The transformed twoelectron integrals can be worked out as

Appendix B: Particle number consistency
To have a consistent number of particles in the various approximations to the grand potential of the interacting system, we need to guarantee that we do not change the particle number when going from Ω s to Ω, so we require Let us first consider the HF approximation. The HF selfenergy is local in time, so G s can be replaced by γ.
The derivative of the 1RDM with respect to the chemical potential is readily worked out as ∂γ/∂µ = βf g f u 1, so the particle conservation condition (B.2) puts the following condition on the trace of the effective potential of the non-interacting reference system tr v s = tr Σ HF + v ext = 4 α + U + 2w 2 , (B.3) where we used the matrix elements of the HF selfenergy (35) and have setñ = 1. This exactly agrees with the chemical potential in the HF approximation (37), as this leads to M g + M u = 0. Including the correlation part of the GW functional does not lead to a change in the trace of v s . This is most easily established by considering the derivative of Φ GW We find therefore, that the trace condition on v s remains the same as in (B.3) when using the GW approximation in the Klein functional.
In the Luttinger-Ward functional we have some additional terms (48). For the derivative of the linear term with respect to the chemical potential we find directly from (52)

∂ ∂µ
Tr Σ G s which has the following useful property ∂G s,gg (ω) ∂µ = ∂G s,uu (−ω) ∂µ . (B.8) Likewise, the derivative of the GW self-energy with G s inserted, can be calculated directly from its explicit expression in (49 where we used ω k = ω − M k as an abbreviation. Due to the relation between the GW self-energy in (51) forñ = 1, the derivative with respect to the chemical potential obeys the following relation where we used that the direction of summation over the Matsubara frequencies ω p is immaterial. Note that we could drop the convergence factor e ηωp as the integrand already converges without it. So also for the Luttinger-Ward functional (B.3) is the appropriate trace of the single particle potential to keep N = 2.