{'fraction_non_alphanumeric': 0.03381619864672698, 'fraction_numerical': 0.007844608004000438, 'mean_word_length': 4.463502091692991, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Strategies to build functionals of the density, or functionals of Green's functions: what can we learn?

The many-body problem can in general not be solved exactly, and one of the most prominent approximations is to build perturbation expansions. A huge variety of expansions is possible, which differ by the quantity to be expanded, the expansion variable, the starting point, and ideas how to resum or terminate the series. Although much has been discussed and much has been done, some choices were made for historical reasons, in particular, limited computation or storage capacities. The present work aims at examining the justifications for different choices made in different contexts, by comparing ingredients of functionals based on \GF s on one side, and on the charge density on the other side. Of particular interest will be the question of how to build an optimal starting point for the approximation of non-local quantities, making use of near- or far-sightedness, and daring to consider models beyond the homogeneous electron gas. This will include the use of connector approximations. We will also discuss why it is a good idea to build functionals of the density.


Introduction
The general aim of solving the many-body Schrödinger equation is to calculate an observable O. One could, in principle, obtain O as expectation value involving many-body wavefunctions, but in order to describe realistic materials, one has to find a simpler route. One prominent way to go is to express the expectation value as functional of a quantity Q that is simpler than the manybody wavefunction, O[Q]. This raises two problems: first, it is generally not easy to find the functional O[Q], or good approximations to it, and second, Q itself may be an expectation value that is not known a priori. Only if Q is the external potential do we not have to bother about the second question. Otherwise, in the frameworks used here one can obtain Q from an auxiliary system, with an auxiliary potential v aux . Again, however, this auxiliary potential itself is unknown, but can in principle be formulated as functional of Q, which makes the calculations self-consistent. The search for O[Q] has then be replaced by the search for v aux [Q]. This is sometimes equivalent to searching for an energy functional E[Q] and using a variational principle.
In the present work, we will compare two cases for the choice of Q: the density n, and the one-body Green's function G. This means that we compare two frameworks, namely Density Functional Theory 1-5 (DFT) and Green's function functional theory [6][7][8][9][10][11] (GFFT). These may seem to be two distinct worlds -it is often said that DFT is computationally efficient, but lacks systematic approximations 12,13 , whereas functionals of Green's functions are built with systematic perturbation theory, but require more computational effort 14 . Here we will explore the reasons underlying these statements, and discuss possible evolutions.
To make the discussion more focused, we will concentrate on one line of approximations. that powerful approximations for a real system may be obtained if one can profit from knowledge obtained from a simpler system, a model. This can happen in two ways: either one can directly replace the observable or auxiliary potential in the real system with the one of the model system, whose parameters are specified in an appropriate way. This is easy if O [Q] or v aux [Q] do not depend on all details of Q, but, for example, only on its value in a given point x 0 : if O[Q] is in this sense near-sighted, one can choose the parameters of the model system such that the quantity in the model Q m has the same value in x 0 , while the rest of Q m can be different, and therefore simpler. Another straightforward case would be a perfectly far-sighted situation, where O[Q] depends only on an average of Q, and the model system can be chosen such that the average of Q m equals the average of Q.
Usually, life is not that easy, but at least the appropriately specified model system may be a good starting point for a perturbation expansion. The hope would be that the resulting expansion converges rapidly. One has a tendency to consider expansions to be systematic approximations, having in mind converging series, where approximations can be made better and better.
Green's functions functional theory is often thought to be systematic, because one associates it with Many-Body Perturbation Theory 8,9 (MBPT): this implies a specific choice, namely, an expansion in the Coulomb interaction, starting from a model of noninteracting electrons subject to the external potential of the real system. It is mandatory to expand beyond this model, because the model cannot be made good enough to yield a satisfactory prescription of many observables that are at the center of interest of GFFT, such as spectral functions. The convergence problems of the perturbation expansion, though, make that it does not always lead to systematic improvements.
In the framework of DFT the typical model system is the homogeneous electron gas (HEG): this model has an external potential that is very much simplified with respect to the one of the real system, but the model is fully interacting 15 . It is mostly used to approximate directly the auxiliary potential, in many cases successfully; therefore, perturbation expansions are less prominent 2 Density Functionals DFT is a very general theory. In the present work, we focus on DFT for interacting electrons. Moreover, we do not consider problems linked to the spin of the electrons, and always refer to the total charge density n(r).

The density
In DFT, the charge density plays the role of the quantity Q in terms of which functionals are built. In second quantization, at zero temperature and for fixed electron number the density is where |N is the N-body ground state andΨ(r) is the field operator that destroys an electron in point r. This expectation value can be evaluated, for example, by using Quantum Monte Carlo 11,17,18 (QMC). Here we are interested in functionals, and indeed, if one knew the total energy as functional of the density, the latter could be obtained by minimizing the energy.

Auxiliary system
Kohn and Sham 19 have shown that the charge density of an interacting system can be obtained as solution of a non-interacting auxiliary system with an appropriately chosen potential, the Kohn-Sham (KS) potential. This potential can in principle be obtained from a functional derivative of the ground state energy, using the variational principle. However, since the density func-tional of the energy is unknown, so is the Kohn-Sham potential. Moreover, the kinetic energy in the Kohn-Sham system is expressed in terms of its orbitals, because even the non-interacting kinetic energy is not known as explicit functional of the density.

Observables
Only the expectation values of local one-body operatorsÔ 1 are easily expressed as explicit density functionals, O 1 = dr O 1 (r)n(r). Other functionals are unknown and must be approximated.

Approximations
A major effort is made to approximate the exchange-correlation part of the total energy as functional of the density, whereas other observables, in particular spectral functions that are also in principle functionals of the ground state density, seem to be out of reach. For the lack of better approximations, observables are sometimes replaced by the corresponding expectation values calculated in the Kohn-Sham auxiliary system. The such obtained approximation is a density functional, and it may be considered to be an approximation to the true functional, but it is often a poor approximation, except for the Kohn-Sham density itself, which is in principle exact.

Model systems
DFT has been successful early on because Kohn and Sham 19 suggested an efficient approximation for the auxiliary potential, the Local Density Approximation (LDA). This approximation is based on the hypothesis of nearsightedness 20,21 : whereas the unknown exchange-correlation contribution v xc (r; [n]) depends in principle on the density in the entire system, it is dominated by the density close to the point r and can therefore be replaced by the exchange-correlation potential v xc (r; [n m ]) of a model, provided n m (r) = n(r). This allows one to use as model a different HEG for each point r, with the corresponding homogeneous density n m r = n(r). In this way, DFT could profit from the fact that calculations in the HEG could be carried out using QMC 22,23 . The importance of this idea cannot be stressed enough: the QMC results were carried out once and forever and made available to all DFT practitioners, who therefore never had to solve the interacting electron problem. This aspect certainly explains much of the "computational efficiency" of DFT.
The LDA is not exact, and the search for better approximations is ongoing. The present article is not meant to be a review paper, and we mention only the two lines of research that fit to our discussion: first, efforts have been made for an improved use of the HEG, such as the weighted density approximation or the average density approximation [24][25][26][27][28][29] , where the requirement for perfect nearsightedness is dropped by taking the HEG at an appropriately averaged density. Second, the most natural idea was to use the LDA as starting point for further expansions, as discussed next.

Expansions
Starting from the homogeneous density, the logical next step is to consider density variations through gradient expansions. This leads to explicit functionals of the density and its gradients. Since a straightforward gradient expansion 1,19,30,31 has convergence problems, successful generalized gradient approximations [32][33][34][35][36][37] (GGA) make use of additional knowledge, such as exact constraints, or parameters fitted to prototype systems.
There have also been attempts to use straightforward Taylor expansion around the LDA 19,38 . For v xc such an expansion reads v xc (r, [n]) = v xc (r, n h r ) + dr where n h r is a homogeneous density distribution with value n(r). There has been discussion concerning the density distribution n 0 where the derivative should be taken 28,38 ; strictly speaking, the Taylor expansion would prescribe n 0 = n h r , but other choices may improve the result when one stops at first order. No conclusive results were obtained, and the attempt was limited to an expansion of the total energy or v xc , but, as we will discuss later, the general idea of expanding in the density around an appropriate starting point may still show its power.
Also expansions in the Coulomb interaction have been proposed 39 . These do however lead to orbital dependent, not to explicit density functionals. While being systematic, the such obtained approximations lose the extreme computational simplicity of LDA or GGAs. Finally, much effort was devoted to semiclassical expansions, see 2 .

GF functionals
Like DFT, working with functionals of Green's functions is a very general topic. Here we will focus on functionals of the one-body Green's function for interacting electron systems, and, as above, do not consider spin explicitly.

The one-body Green's function
In Green's functions functional theory, most often the one-body Green's function plays the role of Q. In second quantization, for a static system at zero temperature and for fixed electron number the time-ordered one-body Green's function is defined as where T indicates the time-ordering of the operators. The diagonal of this Green's function, taken as r = r and t → t + , equals the density Eq. 1. The density matrix is obtained in the same time limit, considering diagonal and off-diagonal elements in space.
Since the external potential is supposed to be static, the Green's function depends only on the time difference t −t or, equivalently, on one frequency G(r, r , ω). The poles of G(ω) are electron addition and removal energies. Therefore its spectral function, which is proportional to the imaginary part, is a good approximation for direct and inverse photoemission spectra 40 . The smallest differences between addition and removal poles correspond to the fundamental gaps.
As in the case of the density, one would like to avoid calculating the Green's function through a direct evaluation of the expectation value Eq. 3. QMC methods are almost exclusively limited to calculations on the imaginary energy axis, and the analytic con-tinuation to the real axis, which would yield spectra, is notoriously difficult 41 . The way to proceed is therefore analogous to the DFT approach: one can build an energy functional [42][43][44][45] of G, and the physical Green's function is obtained at its extrema, which however, contrary to DFT, do not have to be minima.
If one evaluates Eq. 3 for the Kohn-Sham auxiliary system, one obtains where φ v and ε v (φ c and ε c ) are orbitals and eigenvalues of occupied (unoccupied) Kohn-Sham states. The Green's function of the exact KS system is different from the exact G; the latter has in general a much richer spectral function, and the fundamental gaps differ because exact KS eigenvalues are not electron addition or removal energies. Nevertheless, G KS is often used as approximation, or as valid starting point for further calculations.

Auxiliary system
The variational approach to the energy as functional of the Green's function yields the auxiliary system, similar to DFT. It can be formulated as a Dyson equation, where G 0 is the non-interacting Green's function of the system, the compact arguments 1 ≡ (r 1 , σ 1 ,t 1 ) denote space, spin and time, and the bar indicates integration: The self-energy Σ(3, 4) consists of the Hartree potential δ (3 − 4)v H (r 3 ) and the remainder Σ xc , which contains all exchangecorrelation effects and plays the role of a space, spin and time non-local auxiliary potential. When Σ xc is replaced by the exchange-correlation part of the KS potential v xc , the Dyson equation yields the KS Green's function.

Observables
Expectation values of all space local and non-local one-body operators are simple functionals of the Green's function, since its time diagonal is the one-body density matrix. Moreover, the spectral function is where denotes matrix elements in a basis. Finally, the total energy can be expressed in terms of G using the Galitzkii-Migdal expression 46 . Therefore the main problem for many interesting applications, such as the calculation of photoemission spectra, is not to find the Green's function-functional for the observable, but to find G itself. Other functionals for observables related to N-body operators with N ≥ 2, such as absorption spectra, have instead no simple expression in terms of the one-body G, and must be approximated 14,47,48 . This also holds when one wishes to use variational energy functionals, which yield the same result as the Galitzkii-Migdal expression when G is the solution of the Dyson equation, but may be better for an approximate G that is not calculated self-consistently from the Dyson equation.

Approximations
To find G one in general uses the Dyson Eq. 5, and analogous to the case of DFT, the problem is then to find good approximations for the auxiliary potential, i.e. the self-energy, as functional of G. By the way of contrast, there are also attempts to approximate directly the Green's function, without passing through an auxiliary quantity 49,50 . As we will see later, this is however difficult when one is interested in the poles of G.

Model systems
In view of the parallels between DFT and GFFT, one might expect similar approximation strategies. In particular, since the LDA is one of the cornerstones for the success of DFT, the question is whether one can use a model system to simulate a real system also in the case of Green's functions. Indeed, the closest approach to the LDA is probably dynamical mean field theory 51 (DMFT) in the single-site approximation: In this approach, one is not interested in the full Green's function, but only in its diagonal in a site basis, G (ω). The auxiliary system to produce this diagonal, in principle exactly, has a local but frequency-dependent auxiliary potential, i.e. a local self-energyΣ (ω); note that it should not be confused with the diagonal element of the full self-energy, i.e.
Σ (ω) =Σ (ω) in general. As in the LDA, the hypothesis is nearsightedness, i.e. one assumes thatΣ jj (ω) on a site j depends only on G jj (ω), and not on G (ω) on sites = j elsewhere. This allows one to importΣ jj (ω) from a model system that has the same G jj (ω) on that site, and can be different elsewhere. The model chosen by single-site DMFT is the Anderson Impurity Model 52 (AIM), which can be tuned to yield the desired G jj (ω) by tuning a bath Green's function. This leads to a self-consistent calculation ofΣ (ω) and G jj (ω). The method is hence strictly analogous to the LDA, with three main differences: (i) the choice of the model system, which in the case of DMFT is meant to be closer to the target applications containing localized electrons, (ii) the fact that even in the model the interacting and non-interacting local Green's function are different, which is not the case for the density in the HEG, and (iii) the fact that the HEG DFT results have been calculated once and forever, whereas the solutions of the AIM are not tabulated and have to be produced for each calculation. This means that DMFT does not benefit from already existing calculations of interaction effects in the model, which makes it computationally much heavier. A good reason for this is of course the fact that the HEG is parametrized by just a number, i.e. its density, whereas the AIM has a much larger parameter space (orbitals, frequency). Nevertheless, one might hope hat with increased storage capacities and machine learning interpolations, this redundancy will be removed in the future. It is also worthwhile to mention a much simpler and approximate approach, which attempts to use the HEG as model system for Green's functions. The connection between the real and the model system is made through the density; indeed, like every expectation value also the self-energy can be formulated as a func-tional of the density. The approximation termed Quasi-Particle Local Density Approximation 53,54 (QPLDA) is meant to describe the quasi-particle states in the spectrum close to the fermi level, and it makes again use of the hypothesis of nearsightedness. As we will discuss later, nearsightedness is easy to define (though not necessarily valid) only for an object that depends on one space coordinate, whereas the self-energy depends on r and r . Therefore, in the QPLDA approach the self-energy is first converted into an effective local potential.
The most important model for GFFT, however, is that of a system of non-interacting electrons, which gives rise to G 0 in the Dyson equation. In this model the external potential equals the one of the real system, but the interaction is put to zero. Contrary to a model such as the interacting HEG, this non-interacting model is not supposed to be tuned: it can merely be used as starting point for further refinements which, as we will discuss in the following, imply expansions around G 0 . One can, of course, imagine alternatives, which could take the form of a simplified but tunable interaction, instead of a rigid putting the interaction to zero, and inspirations might be obtained e.g. by work on rangeseparation 55,56 . Still, there is no established method yet in this sense.
As a final remark, it should also be noted that it would probably be meaningless to tabulate the non-interacting model, since it depends on all details of the real system, through the external potential, although interpolations might be tempted. It is, however, comparably straightforward to solve. Finally, we should add that most often the model is not a strictly non-interacting one, but rather an effective non-interacting system, such as the KS one, which has proven to be a good starting point for further developments of GFFT.

Expansions
Starting from the non-interacting model the most natural thing to do is an expansion in the Coulomb interaction v c (|r − r |) around v c = 0. This is the heart of many-body perturbation theory 8,57 . It gives rise to the diagrammatic expansions in terms of v c and G 0 that are schematically v c + . . .. As we will see later, it is difficult to use expansions when one is interested in functions with poles. Therefore, except for some quantum chemistry work (see e.g. 49,50,58 ), rarely the Green's function itself is expanded, and one rather seeks to approximate the self-energy, i.e. the auxiliary system, as in DFT. This pushes the problem to the next level, because also the self-energy has poles that suffer from the approximations, but the resulting spectra are overall more decent than those resulting from a direct expansion of the Green's function. One may note that other analogues of an auxiliary system are possible, in particular the cumulant representation 57,59 of the Green's function, where one approximates the cumulant C in G ≡ G 0 e C .
Comparison with DFT expansions in terms of the interaction 39 would be an interesting subject on its own. However, in the spirit of standard DFT approximations, and keeping in mind that the realistic possibility to tabulate a model system gives a huge advantage to such an approximation, in the following we will rather focus on expansions around an interacting model system with simplified potential.
To expand around an interacting model system, we have several possibilities, because the self-energy can be treated as functional of the potential, of the interacting or non-interacting Green's function, or even of the density. Most naturally, we would consider Σ xc as functional of G, which would yield Even in a model, even the first derivative will in general have to be approximated, for example, using the GW approximation 60 it consists of all diagrams in the model, except for the GW diagram that is partially evaluated in the real system. For more advanced self-energies in the derivative, calculations become cumbersome: although the derivative Ξ m is calculated in the model and might therefore be tabulated, to evaluate Ξ m (1, 2;3,4)G(3,4) would require more computational effort than the GW approximation. The usefulness of this expansion is therefore not established, although it might be worthwhile to explore it further.

Diagrammatic expansions for explicit density functionals
In spite of all possible shortcomings, expansions remain one of the most promising ways to improve upon a given approximation. After our brief discussion of expansions in terms of the Coulomb interaction or an interacting Green's function, we will dedicate the present section to expansions in terms of the density.

Motivation
We postulate that one very strong criterion for the choice of an approximation strategy is our ability to avoid redundancy: the idea is to profit from already existing calculations, which therefore have to be carried out only once and forever, serving numerous different applications. As discussed above, this is most easily achieved by using an interacting model system with a simple, but tunable potential. Such a model can be considered as the starting point around which we will expand. The most straightforward move would be to expand directly in the difference between the real and the model potentials, leading to corrections that are linear and higher order response. Alternatively, thanks to DFT we can expand in the difference of real and model densities. This option may have a major advantage: the density is the diagonal of the one-body Green's function, and as we have seen, G gives direct access to several interesting observables, and to others, though less immediately, through MBPT. With the tight relation between the Green's function and the density, there is hope that one can have error cancelling, and/or use exact constraints to improve results. Expanding in the density around a model system will yield observables or auxiliary potentials as explicit functionals of the the model density and of variations around it. These functionals will consist of universal building blocks stemming from the model, which are calculated once and forever, and simple integrals to be performed with the density variations. Here "universal" means that these building blocks do not depend on the real system of interest, within the class of real systems that may be described starting from a given model. Of course, the ambition is to find models that are general enough to serve for a huge class of real systems. The HEG plays a special role here, since it is probably the most widely used model in electronic structure calculations; it has also be shown that it can be considered to be the first step in an expansion in 61 . However, with the advent of computer power and storage capacities, it will be most interesting to examine also alternative models, which may have more degrees of freedom.

Expansion in the density around a model system
The strategy is to express an object O[n] in a functional Taylor expansion around O[n m ]. To be general enough, we suppose that O may be non-local in space and time, which includes the onebody Green's function. The expansion reads where the density difference ∆n is defined as We now insert the external potential v ext in a chain rule 62,63 . This is interesting for later discussions, and it corresponds to a choice that one would often make in practice. The chain rule yields: where is the inverse of the static response function of the model system. In the next terms, higher order response functions appear, schematically: where the second-order response function is Note that all response functions are static. O depends only on a time difference, and is fully described by a static ground state, because the system time-independent. Since there is a one-to-one relation between the interacting density, the external potential and the non-interacting density, the expansion can be made in either the interacting or the non-interacting density.

Expansion in the non-interacting density
Let us start with the expansion in the non-interacting density. We can use the above equations, interpreting the densities as noninteracting ones in the real and the model system, n → n 0 and n m → n m 0 , and all model response functions as non-interacting model response functions, χ m → χ m 0 . The latter can be expressed in terms of non-interacting Green's function where τ indicates time differences. The order N response function carries a prefactor N!. The expansion can be most conveniently depicted by diagrams, as indicated in Fig. 1. It contains the linear response term which we depict as a fat circle (∆n 0 ) attached to a wiggly line for (χ m 0 ) −1 (a "balloon"), and larger balloons consisting of a loop built with several Green's functions G m 0 (thin lines) that is also attached to a (χ m 0 ) −1 . Balloons can be attached to other balloons, as shown in Fig. 1 for the first three orders. The bar at the bottom represents O[n m ] or, when balloons are attached to it, its derivatives.
The diagrammatic rules to build the series are the following: • Build a small balloon L and large balloons ( • Attach small or large balloons to a G m 0 , or to O m .
• Large balloons are always dressed by at least two other balloons that can be large and/or small.
• Each large balloon contains a response function beyond linear response; it carries a factor i, and a prefactor N! corresponding to the order of the response function.
• To order O(N), one has a total of N small balloons, and one has to add an overall prefactor 1/N!. If this expansion in the non-interacting n 0 is performed for an interacting system, all the interaction effects are contained in O[n m ] and its derivatives with respect to the external potential. These derivatives are calculated once and forever in the model, as well as all response functions (which remain non-interacting), whereas all the information about the particular system is contained in the non-interacting density. This is a way to separate interaction effects from the specification of a particular system. As can be seen from Fig. 1, once the model is tabulated the main workload consists in evaluating the integrals that attach balloons to each other or to O m . The difficulty of the initial step -the evaluation of the model -depends on the object O, for which one has to be able to calculate derivatives with respect to the external potential in an interacting model system.

Expansion in the interacting density
To expand in the interacting density, one has to replace χ m 0 by χ m and n m 0 by n m . The resulting expansion is depicted in Fig. 2. Now fat circles stand for the interacting density difference and wiggly lines for (χ m ) −1 . Large balloons contain higher order interacting model response functions, which are depicted by a filled area, since they are no longer simple integrals of products of G m 0 nor of G m . The bottom bar represents O m and its derivatives, as before.
The expansion in the interacting density does not separate interaction and system effects: O m and its derivatives as well as the response functions contain only interaction effects, but the density difference contains both system and interaction effects, which makes its calculation not trivial. For the evaluation of some observables O it may be sufficient to use Kohn-Sham DFT with an approximate functional; if O is the KS v xc itself, the calculation may be done self-consistently.
In the chain rule, one can also use the full KS potential instead of the external one. All direct and inverse response functions are then non-interacting KS ones. When the model is chosen to be the HEG, the Hartree and exchange-correlation potentials do One may wonder why the expansion in n should be better than the expansion in n 0 , which is appealing because of the separation of system-and interaction-effects. While it is difficult to make a general point, some insight can be gained when we choose O to be the Green's function. Since G gives direct access to many interesting observables, it is also a key example.

Expansion of the one-body Green's function
When O in Eq. 7 or in its interacting counterpart is the one-body G, the derivatives of G make generalized response functions appear, with a non-locality as indicated by the pre-superscript with reference to the number of space arguments: The whole expansion is hence a chain of direct and inverse response functions. This allows one to make a quick consistency check: the density resulting from this expansion is the diagonal of the final Green's function. This transforms the generalized nonlocal into an ordinary response function; graphically, it consists in closing the bottom bar in Fig. 2 into a filled area. Taking into account the prefactors as explained above, together with the fact that χ m (χ m ) −1 = 1, one can see that the two terms in the third line cancel, as well as each term in the forth line with the term below it in the fifth line. The same holds for higher orders. The final result is the sum of the first two terms, which yield consistently the exact density. This consistency is an important feature of the expansion; it means for example that if one is interested in O the density matrix, its diagonal is automatically fixed to the exact result, to any order of the expansion. This is instead not true when one expands in the non-interacting n 0 , unless O is the noninteracting Green's function. It is a strong argument to expand in a density that is consistent with the object that is to be expanded. It also gives an argument in favour of expanding in the density instead of expanding in the external potential.
The result of these expansions are explicit functionals of the model density and of the density difference, but not of n m + ∆n. If n m is fixed and ∆n freely variable, this is equivalent to a density functional. If model density and density difference are set up as complementary components of the density, this is instead not true for a finite-order expansion.

Other expansion schemes
One can obtain the balloon expansion also by iterating the Dyson equation twice. Let us use for simplicity the expansion of the noninteracting Green's function, and expand the Dyson equation: where v ext is the difference between the real and the model potential. Taking the diagonal, we obtain and iteration leads to For the Green's function order by order in the potential, we have To infinite order, this is exactly the balloon expansion for the non-interacting Green's function, Fig. 1, with L the small balloons. One can now identify terms; for example, the first secondorder term corresponds to the first second-order diagram of Fig.  1. However, while the balloon expansion is order by order in the density, here one could stop the series at a given order in v ext , itself calculated to a given order in ∆n 0 . The resulting density is not necessarily consistent. For example, a first-order v ext,(1) used to build G 0 to second order misses the second of the second-order terms, which is necessary to cancel the spurious second-order contribution which introduces an error in the density that is of second order in ∆n. This is important, since one might be tempted to use the exact relation between external potential and noninteracting Green's function, G 0 = ((G m 0 ) −1 − v ext [n]) −1 , to build a density functional. This inconsistency is a general property of Dyson equations when their kernel is determined perturbatively; in some cases, e.g. the calculation of quasi-particle energies, solving the Dyson equation improves results, whereas in other cases, e.g. satellites, it deteriorates results 48,64,65 . This has to be investigated on a case-by-case basis.
There are many ways to iterate the Dyson equation explicitly or implicitly. For example, the exact KS potential for small systems has been obtained through numerical iteration [66][67][68][69][70][71] The same procedure can be used to write the KS potential as functional of the density: for example starting with a HEG v  KS ] −1 (r, r,t,t + ), and so on. However, the resulting functionals imply inversions leading to Green's functions, and they would therefore be computationally less efficient if used straightforwardly, calling for additional approximations for the evaluation of the Green's function. 72 5 Connector Theory approach The model system deserves particular attention in the present work. A good choice optimizes the starting point of an expansion and improves its convergence. If the model is flexible enough, it might even be possible to tune it such that it yields directly the desired result of the real system, without expanding beyond the zero-order. This is the hypothesis underlying connector theory 16 (CT), which we will briefly outline in the following, limiting ourselves to aspects that are important for the present work; a more complete description can be found in Ref. 16 .
The aim is to tune the model system such that O[Q] = O[Q m ], where Q can be the external potential, interacting or noninteracting density or Green's function, etc. If the model is flexible enough such that the equation can be satisfied in principle, one can try find the one or more Q m for which the equality holds. Of course, in order to find the exact Q m one would have to know the final solution. However, as explained in Ref. 16  Expansions are good candidates for such approximations. In particular, if Q depends on a set of coordinates x, Q(x), and if Q m is instead simply one number, the first order connector Q c reads where f (x ) ≡ (δ O/δ Q(x )) |Q=Q m . If this first order approximation to the connector is meaningful, f (x ) is the quantity that determines the range of x that is important to describe O. Therefore The workload to evaluate this first-order connector is determined by the calculation of the first order. However, using the connector can improve the result with respect to that of the firstorder expansion itself, as we will illustrate.
6 Spatially non-local observables: the density matrix To become concrete, we will study the one-body density matrix (DM), which is linked to the Green's function as n(r, r ) = −iG(r,t, r ,t + ). The DM is a fundamental ingredient for the understanding and calculations of many-electron systems [73][74][75] : for example, it gives direct access to the kinetic and exchange energies, and the occupation numbers, which are its eigenvalues, indicate the degree of correlation of a system. At the same time, it is non-local, in the sense that it depends on two space arguments: it is therefore not straightforward 19,28,38,[76][77][78] to propose an approximation such as the LDA, since one has to define nearsightedness in a more general way.

Illustration: the single-electron density matrix
In order to get more insight, it is useful to look at a situation where analytically exact results are available. The simplest case, independently of the external potential, is that of one single electron.

Exact relations
For one electron, the occupied orbital is the square root of the density, and the exact density matrix reads n(r 1 , r 2 ) = n(r 1 )n(r 2 ).
This DM depends only on the densities in r 1 and r 2 : this is a generalized nearsightedness for an object evaluated in points (r 1 , r 2 ). Starting from the one-electron Schrödinger equation for the occupied orbital, − ∇ 2 2 n(r) + v ext (r) n(r) = ε n(r).
one obtains the external potential, a part from the constant ε, v ext (r; [n]) = 1 n(r) the potential is not perfectly nearsighted in the density, but depends locally on the density and its gradients. For the single electron the DM is a frequency integral of the Green's function it depends on the external potential in a non-local way, and it would be difficult to use nearsightedness to build a potential functional. Let us now consider the exact expressions as benchmark for approximations based on expansions.

Expansions
To get some insight, let us first suppose that we use a Taylor expansion around a model system to approximate the DM.

First order expansion in the density. A first order expansion of (22) around a model density yields
The same expression can be obtained from the balloon expansion, Sec. 4, for the time-diagonal of the non-interacting Green's function. The static direct and inverse response functions of the model in terms of the model one-electron eigenvalues ε n and wavefunctions φ n are: where v (c) are the only occupied state (empty states). For simplicity, all wavefunctions are taken to be real. The inverse response function is * (28) and the derivative of the density matrix, For the first order correction we evaluate where the orthonormality of orbitals and the completeness relation were used, and integration with ∆n(r 4 ) yields because dr 4 ∆n(r 4 ) = 0. This equals therefore (26), the expansion of the exact result (22) to first order in ∆n.
One might think that the nearsightedness is a direct consequence of the short range of the response. However, let us examine 3 χ 0 . For the one-electron case, we have χ 0 (r 2 , r), (32) * As can be seen from (30), strictly speaking this inverse is missing a contribution, which does however not contribute to the first order result (31). In other words, it does not change (χ m 0 ) −1 ∆n as long as the number of electrons remains unchanged between the model and the real system. which immediately leads to the generalized nearsightedness, because χ 0 χ −1 0 is always short ranged. Note that the response function alone could even be very long ranged. Of course, if 3 χ 0 is not of form (32) but the response functions are (generalized) shortranged, the resulting DM is also generalized nearsighted. It remains to be seen what is dominating the nearsightedness range in systems of more than one electron. Note that here "nearsightedness" only applies to the difference ∆n, not to the density itself; this is perfectly suitable for our purpose of finding approximations given the results of the model system, but different from the more general question of the behaviour of the system. If the model is the HEG, then "nearsightedness" refers to the nearsightedness with respect to density variations, which is closer to the general concept.
Here we have examined only the linear response. From the direct expansion of (22) and from the balloon expansion we can, however, expect that the structure of the problem does not change at higher orders, and similar arguments should apply.

Choice of the zero order.
The result of Eq. 26 and higher orders depend on the choice of the model density n m (r), to be precise, on ∆n(r 1 )/n m (r 1 ) and ∆n(r 2 )/n m (r 2 ). These ratios must be small for the expansion to converge. This can be used as guideline to choose the zero-order of the expansion.
The simplest choice would be to start the expansion with one and the same HEG to approximate all elements of the DM, n m (r 1 ) = n m (r 2 ) = n m , for example, the average density of the system. However, whereas this might be good if the density is quasi-homogeneous, in strongly inhomogeneous systems one could have ∆n(r 1 )/n m ≥ 1, which would lead to divergence of the series. Instead of the average density, one could take the average between the highest and lowest occurring density, (n max +n min )/2. In that case ∆n(r 1 )/n m ≤ 1 everywhere. Still, the series will converge less well when the density is close to the maximum or minimum in points r 1 and/or r 2 , and extremely slowly when the minimum density is zero. So, such an expansion would work straightforwardly for systems where the density variation has a small amplitude compared to the average density.
The choice of the HEG as model system can be optimized further by allowing a different HEG for each pair (r 1 , r 2 ). The most natural choice would be a homogeneous density defined as n m r 1 r 2 (r) = (n(r 1 )+n(r 2 ))/2, independent of r but different for each pair (r 1 , r 2 ). In that case, ∆n(r 1 )/n m (r 1 ) = (n(r 1 ) − n(r 2 ))/(n(r 1 ) + n(r 2 )) ≤ 1, as desired, and similarly for r 2 . Again, however, one would expect bad convergence when one of the densities approaches zero. Note that taking n( r 1 +r 2 2 ) would be meaningless. To go further, let us make the hypothesis that even beyond a single electron the density matrix is generalized nearsighted, in the sense that for a matrix element at given (r 1 , r 2 ) only the density near r 1 and r 2 is important. Supposing, as in the LDA, that close to those points the density is slowly varying, one could build an inhomogeneous model, with the requirement that (n(r i ) − n m r 1 r 2 (r i ))/n m r 1 r 2 (r i ) 1 for i = 1, 2. The density n m r 1 r 2 (r) of this model would in general be different for every (r 1 , r 2 ). Although inhomogeneous, such a model should still be simple enough to be solved with advanced methods for a whole series of its parameters. An example could be crystals with only one fourier component. Such an enterprise would not have been thinkable when DFT was born, but is today within reach; we will come back to this in the outlook.

Higher orders.
If the model is homogeneous, for whatsoever n m the first order Eq. 26 yields n(r 1 , r 2 ) = n(r 1 ) + n(r 2 ) 2 .
For a general non-homogeneous model the result reads: The second order correction is of course small when ∆n is small, but it is also small when the ratio between the real and the model densities is similar in r 1 and r 2 , which is a less severe requirement that may be satisfied by a non-homogeneous model.

Connector approximations.
In the connector philosophy, we do not want to expand the result around some point, but optimize the zero-order, meaning that we find a model system that gives directly the correct result. As explained above, one way to find the parameters of this model is to use a low-order expansion. So the expansion is not used to produce directly the final result, but only to determine the parameters of the model, from which the result is taken. We will examine the connector approximation for the single-electron case, where we do not need it of course, but where we can hope to learn something about connector approximations based on expansions.
The DM for the single homogeneous electron is peculiar, because its density matrix is linear in the density: it is simply the density itself. Therefore, the first-order expansion is exact in the homogeneous model, so which is equal to the first-order approximation itself: since the approximation equals the exact solution in the model, there is no error cancelling, and nothing is gained by the connector. In order for the first-order connector to be useful in the single-electron example, the model system must therefore be inhomogeneous.
In the following, we will test different approximations, including various connector approximations, for single-electron systems. Instead, in Subsec. 6.2 we will use the exact single-electron results as approximation for the CT approach to the many-electron case.

Performance of the approximations in the singleelectron case.
Although the single-electron case can be solved exactly, testing approximations to it will give guidelines for the many-electron case. We therefore start with a numerical illustration for a single electron with density n(r) = A 1 cos(a 1 r) + A 2 cos(a 2 r) + A 3 cos(a 3 r) + B where a = ax in a cube of side length 2π.
To study the performance of different approximations, we show in Fig. 3 the relative error of the approximate density matrix with respect to the exact one, with r 1 and r 2 in directionx. To characterize the results in a single number, the captions also show the mean relative errors (MRE). The left column displays the result of first-and second-order expansions around different starting homogeneous densities. The same scale has been imposed to all results, which implies that in the bright yellow regions the error exceeds the maximum error set by the scale. The upper panel is the first order result. The error is largest in points (r 1 , r 2 ) for which the density is very small in one point and and large in the other, which may lead to ∆n(r i )/n m ≥ 1. As pointed out above, the firstorder result does not depend on the homogeneous starting point. The next two panels are second order results, starting from the mean density of the system, or from n m r 1 r 2 = 0.5(n(r 1 ) + n(r 2 )), respectively. In the first case, the error is larger than that of the first order result, with maxima where one of the two densities is small, i.e. in regions where already the first order is problematic. In the second case instead, where a different homogeneous system is chosen for each pair (r 1 , r 2 ), the result improves significantly.
The right column of Fig. 3 shows the result of the first order connector approximation. As explained above, we have to use an inhomogeneous model for a meaningful connector approximation based on a first-order expansion. We choose a model system with a density that is constrained to be much simpler than that of the real system: we allow only one Fourier component, n m (r) = A m cos(a m r) + B m , with parameters A m , a m and B m that can be varied to match the connector equality. To approximate this equality, as one would do in a real material, we use a first order expansion around A m = 0 and A i = 0, i.e. a homogeneous density for the real and the model DM: this is the same approximation as the one that gave the first panel in the left column, but now, it is used within CT.
The three model parameters A m , a m and B m cannot be uniquely defined by the connector condition . This require to impose additional constraints and solve for the remaining free model parameter(s), to finally obtain the set A c , a c , B c . Since the connector approach has never been studied for such a case, it is interesting to explore what is the most promising way to set the model parameters. Fig. 3 shows that the CT result depends strongly on the way the parameters are set: the best result is obtained when the average density of the model, B c r 1 r 2 , is used as free parameter to optimize each pair of points (r 1 , r 2 ), while the amplitude and periodicity of the dominant Fourier component equals that of the real system. As can be seen by comparing the left and right panels of Fig. 3a, using the same first-order expansion the CT result is clearly superior with respect to the direct approximation while, when the model results are tabulated, the workload to calculate a real system is the same in both cases. The CT result worsens when the amplitude of the first fourier component is set to an average value, but it is still better than the direct approximation. The worst results are obtained when the average density is kept fixed, and the connector is set by the amplitude of the oscillation, as shown in the last panel. In all cases, the worst results are obtained when one of the densities is small. Note that in the last panel the CT equation would even yield regions with negative densities; in those cases, the density matrix has been set to zero. The best result is the second-order expansion in panel (c) in the left column; however, it would require in practice a higher computational effort. When the average density B is much lower, e.g., B = 1, results become worse. Still, for n m = (n max + n min )/2 the expansion converges, though more slowly. Altogether, the results indicate that one may obtain practicable approximations using expansions, even around a homogeneous system, and that in particular CT combined with a low-order expansion is a promising direction.

The many-electron density matrix.
Real systems cannot be solved analytically in general, but we can try to use insight from the single-electron case to approximate the density matrix of real systems † . We will do this in the following for the non-interacting, or Kohn-Sham, case.

Nearsightedness of the many-electron density matrix.
The perfect generalized nearsightedness of the density matrix is destroyed as soon as one has a system with more than one orbital, even without interactions. To understand this better, we will look at the first-order expansion. In the case of one orbital, nearsightedness came through cancellations in 3 χ m 0 (χ m 0 ) −1 : the crucial point was the cancellation of an occupied orbital in the integral of (28) with (29), which allowed one to integrate and have a cancellation of the energy denominators transition by transition, and finally use the completeness relation yielding δ -functions. In the case of more orbitals, only when the different orbitals are localized in completely different regions of space and can therefore be treated separately, things will work in a similar way, and generalized nearsightedness will hold. Indeed, in that situation, in each region the density matrix is that of one electron, and off-diagonal elements connecting different regions vanish. Note again that general nearsightedness of the density matrix does not require χ m 0 (r 1 , r 2 ) to be ultra short ranged in |r 1 − r 2 |, as discussed earlier.

The many-electron density matrix of real systems: CT approximations
One might expect that generalized nearsightedness also holds in the many-electron system, as a first guess. To put things on a firmer ground, CT approximations might be helpful to improve the results when the nearsightedness breaks down, since they are † Also in the asymptotic region r 1 , r 2 → ∞ of a finite system the density matrix n(r 1 , r 2 ) behaves like 79 n(r 1 )n(r 2 ) . able to interpolate between the nearsighted and farsighted cases. We will concentrate on CT in the following. First, note that expanding is only one possibility to realise the connector strategy, which requires more generally to do the same approximation on the real and the model system. In order to illustrate the CT for real systems in the probably simplest possible case, and in order to draw maximum benefit from the analytic single-electron case, we will in the following use the single electron as approximation to build the connector.
We will study two real systems, bulk silicon 80 and solid 3 He 81 , on the Kohn-Sham level. Our model system will be homogeneous, characterized by its density. In the single electron approximation, the connector equality reads n c r 1 r 2 = n(r 1 )n(r 2 ).

1-15 | 11
The density matrix n(r 1 , r 2 ) in a pair of points is then obtained in the CT approximation by evaluating the density matrix of the HEG with homogeneous density n c r 1 r 2 . Using the HEG with geometric mean density is quite intuitive, and has been used e.g. as ingredient for an energy functional describing van der Waals interaction 82 .
Figs. 4 and 5 show, for helium and silicon, respectively, the Kohn-Sham density matrix n(r 1 , r 2 ), the direct single-electron approximation n(r 1 , r 2 ) ≈ n(r 1 )n(r 2 ), and the CT approximation n(r 1 , r 2 ) ≈ n(|r 1 − r 2 |; [n c r 1 r 2 ]). Figs. 4 and 5 also show results of the CT approximation with guesses for the connector extrapolated from the single electron. The relative error on the exchange energy calculated over one unit cell is given in in Table 1. Errors are now much larger than in the previous subsection, since the single-electron approximation is far from an infinite many-electron system, but trends are as before: CT significantly improves the results with respect to a direct use of the approximation. Again, errors in the exchange energy are smaller than the error in the density matrix itself. For silicon, which is closer to the HEG, i.e., closer to the model used by CT, even with this very simple approximation the error might even be acceptable for some applications. Also for helium the connector improves the results, but compared to silicon the direct approximation is better and the CT result is worse, because helium is closer to the single electron than to the HEG.

Comments on the interacting density matrix
Finding simple density functionals for the non-interacting density matrix is an important topic, since it could for example allow one to perform orbital-dependent Kohn-Sham calculations, or speed up calculations based on hybrid functionals. The practical part of the present work was focused on that topic, and to go beyond would exceed its frame. Therefore, we simply complement the discussion by a few remarks concerning the interacting case. As discussed in Sec. 4, one could expand the interacting density matrix with respect to the interacting or the non-interacting density. While with the second choice one remains, besides the derivatives of the density matrix itself, with non-interacting response functions, the former choice has the advantage of yielding the exact diagonal, the density, to all orders.
One may expect that the generalized nearsightedness is strongly weakened by interactions, since the full response χ feels the long-range Coulomb interaction, contrary to χ 0 . However, as explained earlier, nearsightedness is not based on the response function itself, but on the chain of response functions and their inverse, which gives some hope also in the interacting case. Moreover, using the CT one does not have to guess nearsightedness, but can let the connector set the pertinent range. It seems therefore worthwhile to explore this route further in the future.

Time non-local observables: spectral functions
With the density matrix half way between the density and the one-body Green's function, it is a unifying topic of DFT and GFFT. As we have seen, building density functionals for the DM seems to be within reach. It is, instead, much more difficult for the full Green's function, and in the following we will briefly outline why.
To understand the difficulty, it is again sufficient to look at the case of a single electron. Its Green's function reads This expression is not separable (and therefore, not nearsighted in the generalized sense) because of the sum over empty states. However, it is nearsighted in the generalized sense if one is only interested in electron removal spectra. In that case this suggests to express the non-local Green's function as functional of the local Green's function. This very simple starting point may indicate a direction towards approaches like DMFT, where functionals of the local Green's function are approximated using nearsightedness. Still, this is not an explicit density functional. If one examines the balloon expansion Sec. 4, the difficulty becomes clearer. With respect to the DM, the only change is in the outmost response functions, which carry now two different times at their endpoints.
In the one-electron case, and to first order in the density, let us look at the case t 1 > t 2 which corresponds to electron addition. The calculations proceed as before, but an additional phase factor e −iε c (t 1 −t 2 ) appears in (29). This prevents us from using the completeness relation; after a frequency fourier transform the non-nearsightedness of the resulting contributions is given by ω−ε c −iη ∆n(r 6 ), similarly for r 1 → r 2 . The correction with respect to G m 0 obtained in this way is a non-nearsighted renormalization of the weight of each pole. This expression is nearsighted only for frequencies very far from all electron addition energies ε c . In time space, this corresponds to a neglect of the phase factor, which is legitimate only for very short time differences, shorter than the inverse of the maximum addition energy. To see how the near-sighted result of the density matrix emerges, one can take the imaginary part of the correction in frequency space, and integrate over frequency.
Moreover, 3 χ m 0 has now an additional contribution, due to the fact that one can have t 1 > t 3 > t 2 , which was not possible for t 1 = t 2 . In this case, 3 χ m 0 consists of only conduction states c and one has to integrate , where states c are conduction states that appear in (χ m 0 ) −1 . Contrary to before, now the integral does not lead to δ cc , and there is no indication for nearsightedness. Another important point is that this contribution carries a factor (t 1 −t 2 )e −iε c (t 1 −t 2 ) , which leads to a double pole, 1 (ω−ε c −iη) 2 . Higher order corrections add higher order poles. This is a remnance of the usual expansion of a Dyson equation, where the kernel shifts the poles through a series of multiple poles. The convergence of such a series close to the poles is extremely slow, which indicates that a straightforward expansion of the Green's function in the density is not convenient to get spectral functions.
Therefore, we can conclude that calculating spectral functions, or more generally, parts of Green's functions that exhibit time differences, most likely requires functionals that have a significantly non-local dependence on the density. Moreover, the straightfor-ward balloon expansion in the density difference will have difficulties to shift the poles of the Green's function from the model ones to the poles of the Green's function of the real system. Indeed, usually the problem of shifting poles is overcome when one expands the inverse Green's function, which is equivalent in the present case to expanding the external potential and, in the interacting case, the difference between the self-energies of the real and the model system, with respect to the density.

Conclusions and Outlook
Building functionals of Green's functions or of the density instead of the full many-body wavefunction is an efficient way to calculate observables, if appropriate functionals can be found. Here we have examined two ingredients that have been used succesfully, namely, expansions and the use of model systems. Functionals of the one-body Green's function are most often known as expansions in the Coulomb interaction, and density functionals are often based on the use of the HEG as model system. These choices are, however, not a fatality. In particular, with today's computer and storage capacities it is interesting to explore whether the massive use of interacting model systems beyond the HEG could be a promising way to go, since this allows one in principle to calculate most of the interaction effects once and forever, and the calculations for a variety of real materials would consist essentially in putting together the Lego pieces obtained from the model. The resulting functionals can be functionals of the model density and density differences, or functionals of the model Green's function and its difference to the Green's function of the real system. This direction is particularly appealing to build density functionals, since relatively simple expressions can be obtained. In order to examine what could be promising, we have concentrated on the density matrix as object to be approximated. From our analysis and numerical results, we conclude that there is indeed hope to design practical approximations for the density matrix as explicit functional of the density. On one hand, expansions may converge fast for systems where the variations in the density are clearly smaller than the average density. The choice of the starting point of such an expansion is crucial. In the simplest case one expands around the HEG, and the best results are obtained when a different HEG is taken for each pair of points (r 1 , r 2 ) in n(r 1 , r 2 ), namely, a HEG with density n m r 1 r 2 = (n(r 1 ) + n(r 2 ))/2. With respect to the expansions, the use of connector theory improves the results while requiring the same computational effort, and appears the way to go. It also allows one to use different approximations, not based on expansions, which might be helpful for systems with large variations with respect to the average density.
Finally, neither expansions nor CT approximations are bound to the HEG as model system. As an outlook, we have examined single elements of the density matrix of helium within connector theory and the single electron approximation. In Sec. 6 the CT calculations were done using the HEG as model system, and we can directly compare the results to what is obtained when the model is instead a system with periodic density, described by one Fourier component in each direction. Such a model is described by a limited number of parameters, and could still be tabulated if it turns out to be useful. Using the single electron approximation 1-15 | 13 and imposing the CT equality requires n(r 1 )n(r 2 ) = n m (r 1 )n m (r 2 ), which can be fulfilled in many ways. The most intuitive is the generalized LDA, n(r i ) = n m (r i ) for i = 1, 2. We have examined a few elements of the density matrix where the real and model density in the two points are similar, and found that the inhomogeneous model improves the CT results by an order of magnitude. This is to be expected, since CT heavily relies on the model and becomes exact in the limit where the model equals the real system. However, it is encouraging to see that the improvement happens so fast.
For our discussions, we have chosen crude approximations, namely, low-order expansions or the single electron approximation. This has allowed us to highlight differences. One could do much better without losing efficiency; for example, for a periodic system the single electron approximation could be replaced by that of a periodic array of single electrons. We expect that in this way the errors on quantities such as the exchange energy could be brought into an acceptable range quite easily, and that largescale calculations could benefit from tabulated model results, in the same way as DFT in its beginnings has enormously profited from the QMC results in the HEG, through the LDA.

Conflicts of interest
There are no conflicts to declare.