Exact Generalized Kohn-Sham Theory for Hybrid Functionals

Hybrid functionals have proven to be of immense practical value in density functional theory calculations. While they are often thought to be a heuristic construct, it has been established that this is in fact not the case. Here, we present a rigorous and formally exact analysis of generalized Kohn-Sham (GKS) density functional theory of hybrid functionals, in which exact remainder exchange-correlation potentials combine with a fraction of Fock exchange to produce the correct ground state density. First, we extend formal GKS theory by proving a generalized adiabatic connection theorem. We then use this extension to derive two different deﬁnitions for a rigorous distinction between multiplicative exchange and correlation components - one new and one previously postulated. We examine their density-scaling behavior and discuss their similarities and differences. We then present a new algorithm for obtaining exact GKS potentials by inversion of accurate reference electron densities and employ this algorithm to obtain exact potentials for simple atoms and ions. We establish that an equivalent description of the many-electron problem is indeed obtained with any arbitrary global fraction of Fock exchange and we rationalize the Fock-fraction dependence of the computed remainder exchange-correlation potentials in terms of the new formal theory. Finally, we use the exact theoretical framework and numerical results to shed light on the exchange-correlation potential used in approximate hybrid functional calculations and to assess the consequences of different choices of fractional exchange.


I. INTRODUCTION
Density functional theory (DFT) [1] is an approach to the many-electron problem, in which the fundamental variable is the electron density, rather than the many-body wavefunction. DFT has become the method of choice for first principles electronic structure calculations theory of almost every conceivable property of an exceptionally wide range of molecules and materials [2][3][4][5][6][7][8]. DFT is in principle exact, but almost always approximate in practice, as the general mapping between wave-function and density is unknown. Early practical applications of DFT were all based on the Kohn-Sham (KS) formalism [9], which maps the original interacting-electron problem to a non-interacting one, in the sense that the fictitious system retains the same ground state density as the real one. In this approach, the quantum properties of the original system are reflected in the Kohn-Sham system by means of a multiplicative potential known as the exchange-correlation potential, V xc , which is a functional of the density.
Early approximations to V xc , notably the local density approximations (LDA) [9] and the generalized gradient approximation (GGA) [10], provided approximate functionals that were explicitly density-dependent. In the early 1990s, Becke [11,12] showed that so-called "hybrid" functionals, in which LDA or GGA exchange is mixed with the non-multiplicative Fock exchange operator, offered further and significant improvements in accuracy. Such hybrid functionals, and many variants thereof, have since proven to be indispensable for a great number of molecular systems (see, e.g., Refs. [13][14][15][16]). More recently they have been gaining increasing popularity in condensed phase calculations too (see, e.g., Refs. [17][18][19][20]). * t.gould@griffith.edu.au ; leeor.kronik@weizmann.ac.il Originally, the introduction of non-multiplicative Fock exchange was thought to be a further, uncontrolled approximation to KS theory. Within a few years, however, it was shown that it can also be viewed rigorously [21], in terms of a generalized Kohn-Sham (GKS) theory [22]. In this theory, one maps the original interacting system into a partially interacting system that is represented by a single Slater determinant, again such that the ground state density is retained. GKS theory then provides a conceptual advantage over KS theory (which itself emerges as a special case of GKS theory) in that it allows more freedom in functional construction. This is because while there is only one exact KS map to the non-interacting system, there are infinitely many partially-interacting systems to which an exact GKS map may be formed.
Just like KS theory, GKS theory is also exact in principle but approximate in practice. In fact, seeking new approximate density functionals that offer better accuracy and/or a reduction in computational cost is an ongoing effort [23]. In KS theory, such development has been greatly advanced by studies of exact definitions and relations, e.g., the adiabatic connection theorem [24] and precise distinction between exchange and correlation [10] . Furthermore, for KS theory there is a long and distinguished history of aiding such exact analysis by examining the properties of exact potentials for practical systems, obtained through inversion of accurate reference densities [25][26][27][28][29][30][31][32][33][34][35][36]. To the best of our knowledge, neither analytical nor numerical analysis of this kind has been performed for GKS theory. The purpose of this article is to fill that gap.
The article is arranged as follows. In Section II, we present a brief reminder of GKS formalism, followed by proof of a generalized adiabatic connection theorem for exact hybrid functionals. We use the generalization to derive rigorously two different definitions of multiplicative exchange and correlation components -one new, one previously postulated -and examine their differences and similarities, also in terms of density-scaling. In Section III, we present a new algorithm for obtaining exact GKS potentials by inversion of accurate reference electron densities and employ this algorithm to obtain exact potentials for simple atoms and ions. We then use the results to demonstrate and assess analytical theory. In Section IV, we use the exact theoretical framework together with numerical results to shed light on the exchange-correlation potential used in approximate hybrid functional calculations and to assess the consequences of different choices of fractional exchange.

A. A brief overview of GKS theory
The KS equation [9], which as explained in the introduction maps the fully-interacting problem to a non-interacting one, is given by: where ε i and φ i (r) are Kohn-Sham eigenvalues and orbitals, respectively, and V KS ([n]; r) is a multiplicative potential which is a functional of the electron density, n(r), at each point in space r. V KS ([n]; r) is typically expressed as the sum of the external ionic potential term, V ext (r), the classical electron-electron repulsion (Hartree) term, V H ([n]; r) = n(r')/|r − r'|dr', and the "exchangecorrelation" term, V xc ([n]; r), which accounts for all manybody quantum effects beyond classical repulsion via a (typically unknown) functional of the density: In GKS theory [22], mapping to a partially interacting model system that can still be represented by a single Slater determinant is achieved in practice by defining an energy functional of the Slater determinant, S[Φ], or equivalently of the orbitals that comprise it, S[{φ j }], and an associated energy density functional, F s [n], obtained from the Slater determinant that minimizes S[·] while yielding a density n(r), i.e., The minimizing orbitals, {φ j }, play a role similar to that of KS orbitals. According to the Hohenberg-Kohn theorem [1], the total energy, E tot , can be expressed as a sum of the potential energy owing to the external potential and an energy contribution that is a universal functional of the density, F HK [n] (defined in the next Section in terms of many-electron quantities), such that We can then therefore define a "remainder energy" functional, R s [n], by the difference between F HK [n] and F s [n], which naturally depends on the initial choice of S[·]. The total energy is then trivially given by: Eq. (5) can be thought of as an orbital-dependent expression for the total energy [14]. If one chooses to treat the orbitals as implicit density functionals and then take a variational derivative with respect to the density, one obtains a multiplicative KS potential. However, this requires the solution of a nontrivial integro-differential equation, known for historical reasons as the optimized effective potential (OEP) [14,[37][38][39]. If one instead minimizes the energy with respect to the underlying orbitals, subject to the constraint that they integrate to the density n(r), one directly obtains the generalized KS equation: where is a multiplicative "remainder potential" andÔ S [{φ j }] is a generally non-multiplicative operator, obtained from taking functional derivatives with respect to the orbitals S[{φ j }]. This operator obviously depends on the choice of S[·], but does not depend on V ext (r) or V R (r). Importantly, the GKS equation (6) is as rigorous as the KS equation (1) in retaining the original interacting-electron ground state density. As explained in the introduction, one important aspect of GKS theory is that it provides a rigorous framework for the concept of the hybrid functional, which uses the nonmultiplicative Fock operator. This is achieved by choosing [21,40,41]: whereT = − 1 2 ∑ i ∇ 2 i is the many-electron kinetic energy operator,Ŵ = ∑ i< j 1 |r i −r j | is the many-electron Coulomb operator, 0 ≤ α ≤ 1 is the fraction of Coulomb repulsion used in the choice of S[Φ α ], and Φ α is the Slater determinant that minimizes F S α [n], given by and therefore: whereV F is the Fock operator, defined bŷ We emphasize that Eq. (11) is exact for any choice of α. We note that the original KS equation [Eq.
(1)] is obtained as a special case of Eq. (11) when choosing α = 0. Another special case of note is obtained for α = 1, for whichÔ S becomes the sum of the kinetic energy operator and the full Hartree-Fock operator. The exact remainder potential in this case transforms the Hartree-Fock approximation into an exact construct, sometimes referred to as the Hartree-Fock-Kohn-Sham equation [2]. Finally, we note the connection between exact GKS theory and approximate hybrid functionals. For a general 0 < α < 1, if one subsequently chooses to approximate V α R (r) as i.e., to introduce semilocal approximations (denoted by the 'sl' subscript) for the exchange [V x,sl (r)] and correlation [V c,sl (r)] potentials, then Eq. (11) becomes the standard equation for an approximate global hybrid functional, given by [12,14,42]: B. The GKS adiabatic connection theorem An adiabatic connection formula [14,24] continuously connects properties of a fully interacting quantum system to properties of a reference system, with the aim of making formal analysis and/or approximations easier. In fact, the original introduction of fractional Fock exchange as an ingredient in approximate exchange-correlation energy expression, i.e., what we now call approximate hybrid functional theory, was motivated by adiabatic connection considerations within KS theory [11,43].
In regular KS theory, an adiabatic connection is achieved by defining the many-body operatorĤ λ ≡T +V ext,λ + λŴ , where 0 ≤ λ ≤ 1 is a parameter andV ext,λ is a λ -dependent external potential operator.V ext,λ is adjusted such that the ground-state many-electron wave-function, |Ψ λ , corresponding toĤ λ , integrates to the original many-electron density, n, for all λ .
In particular, when λ = 1 all interactions are included.V ext,1 is simply the external potential of the true many-electron,Ĥ 1 is the true many-electron Hamiltonian, and |Ψ 1 is the original many-electron wave-function, |Ψ . Similarly,V ext,0 andĤ 0 are the many-electron external potential operator and Hamiltonian corresponding to the Kohn-Sham potential, and |Ψ 0 is the Kohn-Sham determinant, |Φ S . The advantage, then, is that mapping from |Ψ to |Φ S is no longer "abrupt" but rather can be obtained adiabatically by slowly "turning off" the electron-electron repulsion while retaining the overall density.
The many-electron and KS limits can be further connected by defining a λ -dependent universal functional, For λ = 1, F λ [n] becomes the universal Hohenberg-Kohn functional for the fully-interacting system, namely For λ = 0, F λ [n] becomes the non-interacting kinetic energy function, T s [n], of Kohn-Sham theory, i.e., Formulating an adiabatic connection theorem within the GKS framework is more challenging and less direct. This is because GKS theory involves an ansatz for the partially interacting system [Eq. (8)], which does not lend itself to continuous switching from off to on, as done for the interaction in Eq. (15). We must therefore define an appropriate adiabatic connection density functional, F λ α [n], that can accommodate the ansatz. This functional must satisfy several criteria: (i) it must reduce to the GKS result, F S α [n] [Eq. (9)] when λ = 0, for all α ∈ [0, 1]; (ii) it must reduce to the fully interacting result, F HK [n] [Eq. (16)] when λ = 1, for all α ∈ [0, 1]; (iii) It must reduce to the adiabatic connection of KS theory when α = 0, for all λ ∈ [0, 1]. As explained in detail in the next sub-section, in order to be useful for defining exchange and correlation components, we shall also ask that it (iv) allow the attainment of the Hartree-exchange energy from its right derivative with respect to λ at λ = 0 [35], for all α ∈ [0, 1]. This means that the expression should involve only Slater determinants at λ = 0 + .
In keeping with the spirit of Eq. (15) as closely as possible, we seek an expression that involves only terms that are explicitly linear in λ , which we refer to henceforth as a linear adiabatic connection (LAC). This expression is: and |Ψ λ =0 ≡ |Φ S , it is straightforward to verify that F 0 α = F S α and F 1 α = F 1 , so that (18) satisfies conditions (i) and (ii). Furthermore, Eq. (18) reduces to the Kohn-Sham adiabatic connection equation (15) for α = 0 and any λ , because |Φ α=0 ≡ |Φ S , thereby fulfulling condition (iii). Finally, condition (iv) is also trivially obeyed.
Importantly, the LAC of Eq. (18) is a simple adiabatic connection, but it is not the only one possible. Another example of an adiabatic connection. which is shown below to be of importance, is that of a quadratic adiabatic connection (QAC) of GKS theory, defined by: The quadratic term vanishes for λ = 0, λ = 1, or α = 0, so that conditions (i)-(ii), respectively, are still obeyed. Furthermore, the quadratic term is phrased entirely in terms of Slater determinants, so that condition (iv) is also still obeyed. Both definitions are equally exact. A possible advantage of the LAC definition over the QAC one is that the linearity in the coupling parameter λ makes it more easily amenable to Görling-Levy perturbation theory analysis [44,45] and to future generalization, e.g., to GKS maps beyond global hybrid functionals and/or to degenerate ground states and other ensemble DFT scenarios [35,[46][47][48][49][50].

C. Exact remainder Hartree-exchange and correlation components
In KS theory, dividing energies and potentials into Hartree (H), exchange (x) and correlation (c) components has historically been an important aid in interpreting and approximating density functionals. In the absence of degeneracies, which is the case throughout this article, such division is straightforward [10]. In GKS theory it is more complicated.
The Hartree energy is an explicit and known functional of the density. The Hartree-exchange (Hx) energy of KS-DFT is, thereby defining its x component as the usual exchangeenergy expression of Hartree-Fock theory, but evaluated using the KS orbitals. The correlation energy is then given by . In GKS theory, the Hartree energy is defined just as in KS theory, but breaking down the remainder potential rigorously into Hx and c components is no longer obvious, because some of the Hx is necessarily incorporated into the non-multiplicative operator,Ô S . Görling and Levy (GL) [21] suggested that this quandary can be overcome by using an ansatz that defines the remainder Hx energy in terms of the KS Hx energy, namely: where here and throughout bars indicate quantities relevant to the GL definition. Using Eqs. (4), (5), (9), and (21) the remainder correlation energy is then given by: The above definitions of remainder x and c are intuitive and add up to the correct remainder xc energy by consruction. Nevertheless, with a view towards both present analysis and future generalization, it is important to examine whether partition between remainder exchange and correlation can be derived generally and without an ansatz. Our starting point for such analysis are cases where it is already known that the conventional definition of the Hx energy of Eq. (20) breaks down, notably where |Φ S is degenerate, so that use of Eq. (20) no longer produces a unique value -the "non-uniqueness disaster" of Ref. [35]. Standard application of theory then leads to spurious "ghost interactions" [51] in open shells [52] and dissociation [35]. In such cases, the right derivative with respect to λ of the adiabatically connected universal functional F λ provides a more rigorous way to define E Hx , that is consistent with conventional theories [35]. We therefore adopt this broader definition in our analysis and apply it to exact GKS theory by defining the remainder Hartree-exchange energy of a GKS system as: Using the above definition with the LAC of Eq. (18) yields: which yields: with the exchange part then given by . Specifically, using Eqs. (9), (20), and (17), we can write: The remainder correlation energy is then defined as Using Eq. (25) in Eq. (27) lets us obtain: The above LAC expressions for the remainder Hartreeexchange and correlation, namely Eqs. (26) and (28), respectively, provide a partitioning of the Hx and c energies that is different from that of the GL expressions for the same quantities, given by Eqs. (21) and (22). Straightforward algebra shows that as required, but where for all α>0. Some further and equally straightforeward algebra shows that if the remainder Hx definition of Eq. (23) is used on the QAC of Eq. (19), rather than on the LAC of Eq. (18), then the quantities of GL theory are obtained. One could argue that this result is somewhat contrived, because it is precisely the quantity ∆F S α [n], which distinguishes between the LAC and GL definitions, that we chose to add as the multiplicative factor of the λ -quadratic term in the QAC expression of Eq. (19). Nevertheless, the result does show that GL quantities can also be reconciled with the general relation of Eq. (23) between the Hartree-exchange remainder energy and the adiabatic connection, if one is willing to sacrifice the linearity of the latter.
Physically, the difference between the LAC and GL expressions is that in the former definition the exchange energy changes with α and the correlation is α-independent, whereas in the latter definition the opposite is true. Clearly, this difference must be inherited by the corresponding potentials, V R,Hx = δ E R,Hx /δ n and V R,c = δ E R,c /δ n, which are at the focus of this manuscript. The difference between the two definitions is most easily rationalized by considering the limiting case of α = 1, i.e., Hartree-Fock with exact correlation. From the point of view of GL theory, exact exchange is by definition already expressed fully in the Fock operator and therefore the remainder potential must be entirely due to correlation. However, from the point of view of LAC theory, there is in fact a remainder exchange component, given by Eq. (25), which is non-zero even for α = 1. The physical origin of the remainder exchange is the difference between the KS and GKS Slater determinants, |Φ S ≡ |Φ 0 and |Φ α , i.e., the fact that even when the energy expression is the same in KS and GKS theory, the underlying orbitals are not, owing to the use of a multplicative or a non-multiplicative potential, respectively.

D. Behavior under uniform coordinate-scaling of densities
An important property of exact density functionals is their behavior under a uniform coordinate-scaling of the density, i.e. upon varying γ in n γ (r) ≡ γ 3 n(γr). Such scaling arguments have been enforced as part of the construction of important non-empirical density functional approximations, notably the Perdew-Burke-Ernzerhof (PBE) [53] GGA functional and the strongly constrained and appropriately normed (SCAN) [54] meta-GGA functional. It is therefore of immediate interest to examine exact scaling behavior from the perspective of exact hybrid GKS theory.
Under uniform coordinate-scaling, the following three exact relations in KS theory have been established for the kinetic, Hx, and correlation terms, respectively [10,55,56]: The first two relationships involve direct scaling of the functional. The last relationship involves the definition namely the correlation energy at interaction strength λ (note, E c ≡ E 1 c ). Here, Ψ λ [n] is the wavefunction at scaled Coulomb interaction,Ŵ → λŴ , defined as The direct scaling relationships for T s and E Hx follow from the scaling ofT andŴ , because ∇ 2 r → γ 2 ∇ 2 γr and 1 |r| → γ 1 |γr| . To derive Eq. (35), consider that under coordinate scaling, Ψ λ [n γ ] → Ψ λ /γ [n] because of the relative scaling ofT and W in (37). Use of Eq. (36) then leads to One γ pre-factor in (38) arises from the Coulomb potential.
The second γ pre-factor and 1/γ limit in the integrand of (38) arises from the variable change (λ /γ → λ ), which accommo- For the above GKS theory, the only additional scaling that needs to be investigated is that of ∆F S α [n], given in Eq. (32), which appears as an additional ingredient in either the exchange (LAC theory) or the correlation (GL theory) component. By applying the same arguments used to derive Eq. (35) to the constrained minimization over Slater-determinants used to derive ∆F S α , we obtain where Φ λ [n] = arg min Φ→n Φ|T + λŴ |Φ involves a restricted minimization over Slater determinants. Then, because the Slater determinants are unique, we can use the same arguments behind Eq. (38) to show that, Clearly, the scaling behavior of ∆F S α is similar to that of KS correlation (though not identical, except for α = 1) and different from that of KS exchange. From the LAC point of view, this is simply another aspect of the difference between GKS exchange and KS exchange. From a GL theory point of view, this is simply because ∆F S α is already part of the correlation.

III. NUMERICAL INVESTIGATION OF EXACT REMAINDER POTENTIALS
A. An inversion algorithm We now wish to examine some of the theoretical concepts of the preceding section in numerical practice. As mentioned in the introduction, exact KS theory is often examined by studying the behavior of exact KS potentials, obtained from inversion of the KS equation based on an accurate reference density. This has been achieved using many different approaches [25]. Several iterative algorithms are based on the linear response of the potential to the difference between the density at a given iteration and the reference density [26,27]. Other algorithms include in the potential a Coulomb term that enforces the correct asymptotic decay of the density [28][29][30]. Yet other approaches employ orbital-by-orbital corrections to the potential, in the spirit of orbital shifts introduced in the Krieger-Li-Iafrate (KLI) [57] approximation to OEP [31]. More recently, Peirs et al. included correction terms that combine the local linear response and the enforcement of the asymptotic decay [32]. Even more recently, it has been suggested that these seemingly different methods can be derived from one general algorithm [33]. Inversion has also been used to find an effective local potential from Hartree-Fock densities and compare it to the KS potential obtained from exact exchange OEP calculations [34].
To the best of our knowledge, an inversion algorithm that accounts for a non-multiplicative operator has not been implemented and exact GKS remainder potentials have been neither obtained nor investigated. As a first step, we therefore devise an inversion algorithm, shown schematically in Fig. 1, that can work within a GKS formalism.
Due to the orbital nature of the GKS problem, it makes sense to adopt an orbital-based solution for inversion. Our goal is to approach the given reference density n(r) (to within a specified numerical precision) by iterating on a density n φ (r) = ∑ i |φ i (r)| 2 , determined from GKS orbitals. To that end, we rewrite the exact hybrid functional equation (11) as: where V (r) is the full multiplicative potential acting on the GKS orbitals, φ i (r). Upon an initial guess for V (r), any of the obtained GKS orbitals can be used to calculate a current potential, V φ (r), via As n φ → n, it follows from uniqeness (and assumed smoothness) that V φ → V . However, when iterations are not yet converged, n φ = n and therefore V φ = V . Our approach thus seeks to renormalize the current orbitals, φ i (r), at each iteration, so as to provide "pseudo-orbitals", ψ i (r), which are guaranteed to give the correct density n(r). That is, we define: such that ∑ i |ψ i (r)| 2 = n(r) by definition. This transformation comes at a cost, however, as the pseudo-orbitals are not orthogonal, nor are they solutions of Eq. (41). Two adjustments to (42) are therefore needed in order to iterate to the correct potential: 1) the energy ε i must be re-evaluated based on ψ i (r); 2) as the pseudo-orbitals are neither orthogonal nor obtained from the same multiplicative potential, their corresponding effective potentials, V i (r) = ε i + [∇ 2 ψ i (r) − αV F ψ i (r)]/ψ i (r), will generally differ and therefore must be averaged in some way. The first issue can be dealt with by applying perturbation theory to the kinetic energy, to give Formally, one should do the same for the non-local Fock operator, but in practice we did not find this to be necessary. The second issue can be dealt with by taking a weighted average of the potentials, i.e., by setting Combining everything, we can thus construct a new potential by using: where ε i are the current orbital energies, andV F is the Fock operator based on the current orbitals. In practice, we typically mix the potential obtained from (44) with that of the previous step to aid in convergence. Finally, based on Eq. (11) we subtract V ext (r) and αV H (r) to obtain the exact remainder potential V R (r).

B. Exact Remainder potentials
The algorithm described above can in principle be implemented in any code, albeit with varying levels of expected success owing to different numerical issues. In order to be confident in the accuracy of the inversion, we wish to remove additional numerical approximations that are almost always used in DFT, i.e., the use of a basis set and/or pseudopotentials. To this end, we have implemented the algorithm described above in DARSEC [58], an all electron fully numerical (i.e. basisset free) code. DARSEC uses prolate spheroidal coordinates [59] to calculate atomic and diatomic systems. It has been specifically designed to be useful in endeavors of functional development [60] and has already been used successfully in cases where use of pseudopotentials or basis sets was a specific cause of concern [61,62]. For performing hybrid functional calculations, we have implemented the Fock operator in DARSEC using the self-consistent scheme described in ref. [63]. We note that efficient projection approximations for the Fock exchange operator have been given [63][64][65][66], but for the small systems studied here use of the full (unprojected) Fock operator was found to be sufficient.
Using these tools, we calculated the exact GKS remainder potentials corresponding to various fractions of exact exchange for four systems: Li + , Li − , Be, and F − . These systems were chosen due to their simplicity, spherical symmetry, and the availability of high-quality reference data. Here we used reference densities obtained by fitting Quantum Monte Carlo calculations using known exact constraints [67]. First, we note that in the limit α = 0, in which exact KS map is retrieved, the obtained V (r) agree well with exact KS inversion results reported previously for the same systems using the same reference densities [67] (see Fig. S1 in the Supplementary Material (SM)). We also note that as a second test, we have successfully replicated the exact exchange potential obtained in Ref. [33] from inversion of Hartree-Fock orbitals for the Ar atom.
With the accuracy of our method established, we now turn to analyzing the results. We find that in all four systems the remainder potential V R (r) changes significantly with α. However, a large component of V R (r) is clearly given by the complementary fraction of the Hartree potential which is not included inÔ S , namely (1−α)V H (r) (cf. Eq. (13)). The Hartree potential is necessarily identical for all choices of α because it is an explicit functional of the exact density and therefore its inclusion, scaled by (1 − α), obscures the observation of trends with α. For this reason, we subtract this Hartree component from the remainder potential and define the result as the remainder exchange-correlation potential, Trends in V α R,xc (r) as a function of α for all four systems are shown in Fig. 2. In the figure, spurious noise arising from convergence difficulties at the empty focus of the DARSEC prolate spheroidal grid has been removed for clarity. For completeness, the unedited figure is given as Fig. S2 in the SM. Figure 2 demonstrates the logic behind obtaining different yet equally exact GKS maps. For every choice of α, the nonmultiplicative operator changes, but the remainder exchangecorrelation potential adjusts itself such that the same ground state density is obtained. This confirms that in exact GKS theory there is no a priori reason to prefer one value of α over another.

C. Partition between exchange and correlation potentials
We now turn to the partitioning of the remainder potential into exchange and correlation potentials. In particular, it is important to examine whether the term ∆F S α [n] of Eq. (32) is meaningful in practice, regardless of whether it is viewed as contributing to exchange (LAC theory) or to correlation (GL theory). This is because although in general |Φ S ≡ |Φ 0 and |Φ α need not be the same (as they are only constrained to yield the same overall density), in practice they can often be close, which is indeed the case for the orbitals found by our inversion process (see Fig. S3 in the SM for an example).
If we assume that the difference is small enough to neglect ∆F S α [n], then from the perspective of LAC theory, this means that we neglect the remainder exchange component for α = 1, i.e., approximate the α-independent V R,c by V α=1 R,xc , which is known from inversion with α = 1. Similarly, from the GL perspective this would mean that we neglect the α-dependence of V α R,c and useV α=1 R,c throughout instead. In either case, we can then approximate the remainder exchange potential for GKS maps with α < 1 as: Anticipating, based on Eq. (21), that a major contribution to the remainder exchange potential is its scaling by 1 − α, i.e., the fraction complementary to the one weighing the Fock exchange inÔ S , we then define a scaled remainder exchange potential:ṽ We have applied Eq. (47) to calculateṽ α R,x from the remainder potentials given in Fig. 2. The results are shown in Fig.  3. Clearly,ṽ α R,x is found to be essentially independent of α, i.e., the multiplicative exchange component of a given map simply scales with 1 − α, which means that ∆F S α [n] is of negligible consequence, to within our numerical accuracy. This nearly-linear scaling relationship holds even in Li − , where one could have expected that the weakly-bound outermost electron would be more sensitive to small changes in the potential.
While the above examples obviously do not prove that ∆F S α [n] would always be negligible, the fact that it often is can still be rationalized. To that end, we return to the analytical form of E α R,Hx . For α = 0 we obtain the usual KS solution. Perturbing the α > 0 GKS wave-functions around α = 0 Using Eq. (48) in Eq. (25), we obtain: We can further show that the second term above, involving |∆Φ 0 , is proportional to α 2 . This involves recognizing that |Φ 0 ≡ |Φ S minimizesT , so that Φ 0 |T |∆Φ 0 = 0 for any perturbation. Therefore, −2αℜ Φ 0 |(T + αŴ )|∆Φ 0 = −2α 2 ℜ Φ 0 |Ŵ |∆Φ 0 , and we obtain: The above result can also be interpreted as a GKS viewpoint of the proof given by Kümmel and Perdew [68], which showed that KS and GKS maps are identical to first order. Up to order α 2 , then, the only α-dependence of the exact Hartreeexchange energy is that it scales directly with 1-α. Any missing contributions that are linear in α, if at all, must come from our approximation that V R,c ≈ V α=1 R,xc . Clearly, this is an excellent approximation in the cases studied here.
Finally, we note that as we are unaware of attempts to explicitly approximate ∆F S α [n] or to set bounds on its magnitude. The most popular global hybrid approximations, notably B3LYP [69] and PBE0 [43]) tacitly set it to zero. More heavily parameterized global hybrid functionals may implicitly incorporate some of it, albeit not in an intentional manner.

IV. COMPARISON OF EXACT AND APPROXIMATE HYBRID FUNCTIONALS
Throughout the previous sections, we have emphasized that exact GKS maps based on any choice of α map the exact density equally well. At the same time, we know that practical application using approximate hybrid functionals can be very sensitive to the choice of α. For example, the generally recommended values of α for the popular B3LYP [69] and PBE0 [43,70,71] hybrid functionals are 20% and 25%, respectively. Furthermore, cases involving more significant delocalization and/or self-interaction errors, or larger static correlation contributions, may call for larger [72][73][74] or smaller [75,76] values of α, respectively. The purpose of this Section is to examine whether the GKS remainder potential plays a role in this optimal choice, and to gain insight into scenarios and properties where the transition from the KS to the GKS picture may result in qualitative differences.

A. Remainder exchange-correlation potentials
We choose to focus our comparisons on the PBE0 hybrid functional, which is based on non-empirical semi-local exchange and correlation components. In Figure 4, we compare the exact, inversion-based V α R,xc with the corresponding approximate V α R,xc from a self-consistent PBE0-based calculation with the same value of α (which we denote as PBE0 α ). We use Be as an example, with similar trends observed for all other systems -see Figs. S4-S6 of the SM.
We focus first on the case of α = 1 (panel (e) in the figure). Neglecting the small difference between the correlation in the Hartree-Fock sense and in GKS sense, as in the previous Section, we can compare the remainder exchange-correlation potential in this case directly to the PBE0 correlation potential. Clearly, the two are very different and in fact disagree even in the sign of the potential. This is to be expected. It is wellknown that semi-local correlation is fundamentally incompatible with exact exchange, because the exact-exchange hole is delocalized, whereas the semi-local hole is localized, resulting in an overall non-physical delocalized exchange-correlation hole [14,77]. The present results demonstrate just how severe this incompatibility is.
For the opposite extreme fraction of Fock exchange, α = 0, PBE0 α reduces to the PBE functional. It has been noted previously [78,79]  results even when PBE correlation is opposite in sign. This is confirmed in panel (a) of Fig. 4, where the PBE exchangecorrelation potential is a much better approximation of the exact one than the PBE exchange or correlation potentials individually. Clearly, some differences between the PBE and exact potential remain, notably a small spurious peak, as well as quantitative differences.
The above analysis therefore suggests that some finite fraction of α may be best in balancing error cancellation between semi-local exchange and semi-local correlation while enjoying some of the benefits of exact exchange, e.g., reduction of self-interaction errors. Panels (b-d) of Fig. 4 suggest that this is likely the case, but the differences in remainder exchangecorrelation potentials for various choices of α, in the absence of energy considerations, do not establish an optimal value of α unequivocally. In particular the results shown in the Figure do establish that use of the default value of α = 0.25 is superior to use of α = 0, i.e., PBE, they do not identify this fraction as optimal in mimicking the exact potential. Rather, it emerges as one useful value out of a fairly broad range of possible Fock exchange fractions.
To understand this behavior, we note that the overall difference between the self-consistent PBE0 α and the exact potential is given by: Kim et al. suggested that this difference can be decomposed into a functional-driven (FD) and a density-driven (DD) component, ∆V tot (r) = ∆V DD (r) + ∆V FD (r) [80], where the FD difference is given by: (53) while the DD difference is: where n exact (r) is the exact density and n PBE0 α (r) is the selfconsistent PBE0α density. The results of this decomposition, for the calculations presented in Fig. 4, are given in Fig. S7 of the SM. We find that for all systems studied in this work, the differences are overwhelmingly dominated by the functional-driven component. Therefore, in the systems studied the choice of an optimal fraction is strongly driven by exchange-correlation energy, rather than potential.

B. Steps in remainder exchange potentials
Given the conclusions of the preceding sub-section, which did not reveal strong trends with α, one may think that using a GKS approach over a KS is merely a computational convenience. However, one issue where the potential itself plays a clear role is the step structure at electronic shell closings. Within the KS picture, it has long been known that this step is related to orbital-dependent terms in the optimized effective potential equation and therefore difficult to mimic with semi-local functionals [81,82].
To illustrate this, we return to the Be atom and consider the difference between the remainder exchange potential and the appropriately scaled Slater potential [14], given by using orbitals from α = 1 (varying α leads to differences of less than 1 mHa, see Figure S3, which are indistinguishable from noise in the remainder potential). We adopt this approach because the Slater potential is known to serve as a rough approximation for the KS exchange potential corresponding to Hartree-Fock theory and the step structure has been traced back to terms beyond the Slater potential [14,81,82].
The result is given in Fig. 5. For α = 0, i.e., KS theory, as expected the step is vital. Furthermore, PBE calculations can mimic some rough features of the needed step (and certainly do a better job at that than calculations based on LDA, in which the "step" is not as sharp and attains a higher maximum). Still, PBE calculations struggle to capture essential features of the step, such as its abruptness and its overall height. As one increases α, i.e. uses GKS, the step diminishes due to scaling and so the error is less significant (compare the PBE0 result to the exact result with α = 0.25). In the limit of α=1, the non-local features of the Fock operator make this step unimportant.
Using increasing values of α requires more sophisticated approximate correlation expressions that limit their use. However, the above observations are still important, because a similar step structure in the exchange potential has long been known to arise between dissociating atoms [83,84]. We do not pursue the issue numerically here for lack of suitable exact reference densities. However, we note that the fact that this step need not be mimicked for α = 1 is surely a major factor in the success of asymptotically correct functionals at treating charge transfer scenarios, which are very difficult to handle with KS theory (See, e.g., Ref. [40], and references therein).
Finally, we note that in the presence of strong correlation one may also observe abrupt correlation features in the spindependent KS potential [85]. These would not be captured in the spin-resticted formalism employed throughout this article, but the elimination or at least mitigation of the exchange step is already a major step forward.

C. Orbital energies
An even greater difference between different (exact or approximate) GKS maps has to do with orbital energies. Throughout our discussion, we have emphasized the densityequivalence of all exact GKS maps. However, the orbital energies, i.e., the GKS eigenvalues obtained from different exact maps may still change without violating this equivalence.
Strictly speaking, the above formalism does not attach physical meaning to individual GKS orbitals and orbital energies, except inasmuch as that they conspire to produce the correct density. The only exception to this statement is for the energy of the highest occupied orbital, for which the ionization potential theorem, originally proven for exact KS theory [86,87] but also valid for exact GKS theory [40,88,89], establishes that it is equal and opposite to the ionization energy. Understanding orbital energy trends is still interesting because, although they do not correspond exactly to quasiparticle excitation energies [90], they are often used as useful approximations thereof [91,92].
Differences between the energy of the highest occupied orbital and the lower occupied orbital(s), for Be, Li − , and F − , obtained from both the exact and the PBE0 α hybrid functional, as a function of α, are given in Fig. 6. Interestingly, the orbital energies obtained from both exact and approximate calculations change significantly and mostly linearly with α, albeit not necessarily with the same slope, despite the fact that (as discussed in the previous Section) the orbitals themselves change very little throughout. This shows that even in exact hybrid theory the α-scaling of the exchange components in the GKS Hamiltonian is clearly manifested in the eigenvalues, in agreement with past results from approximate hybrid functionals [92].
Finally, the strong α-dependence of the orbital energies obtained even in exact GKS theory raises an interesting question for further research: how to choose, from otherwise equivalent maps, the one best-suited for approximating quasi-particle excitation energies.

V. CONCLUSIONS
In conclusion, we have presented a rigorous and formally exact generalized Kohn-Sham (GKS) density functional theory of hybrid functionals. First, we provided a brief overview of how GKS theory defines an exact remainder potential that combines with a fraction of Fock exchange to produce the correct ground state density. We then generalized the adiabatic connection theorem to the case of exact hybrid functional theory and showed that this can be done in different ways. We used this to derive two different yet equally rigorous distinctions between remainder (multiplicative) exchange and correlation components, one of which has been previously given as an ansatz, which we compared and contrasted.
We examined the formal theory numerically by developing a novel algorithm for inverting accurate reference electron densities to obtain exact remainder potentials. We used this algorithm for selected atoms and ions to demonstrate how an equivalent description of the many-electron problem is obtained with an arbitrary fraction of Fock exchange. We further found that, at least for the systems studied, the formal component that sets apart the above two definitions for the partition of exchange and correlation is negligibly small, and rationalized this behavior using perturbation theory arguments.
We then compared the behavior of exact GKS calcluations to those obtained from approcimate ones based on the PBE0 α functional. We showed that the normally recommended value of α = 0.25 is governed by energetics rather than by potential. We established that the choice of α may play a far greate role when there is a need to approximate accurately a step structure. Finally, we showed that different, equally exact GKS maps can feature significantly different eigenvalue spectra and discussed how this may affect the approximate interpretation of eigenvalues as quasi-particle excitation energies.
We sincerely hope that the general theory provided in this article will be used to guide further development efforts of hybrid density functionals and that the generalized Kohn-Sham inversion procedures and analysis given here may be used in the future to understand and overcome limitations of existing state-of-art approximate hybrid density functionals.