{'fraction_non_alphanumeric': 0.03232528457190487, 'fraction_numerical': 0.010314917558802844, 'mean_word_length': 4.680970701390944, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-Consistent Density-Functional Embedding: a Systematic Approach for Density Functionals

A frequently found complaint about popular density-functional approximations is the lack of a systematic path how to improve results obtained with a selected functional. In the present work, we introduce a Self-Consistent Density-Functional Embedding technique, which allows to improve a chosen exchange-correlation functional in a systematic way. Our method adopts ideas from other embedding approaches and combines them with the scheme of density inversion methods, that are based on the rigorous relation of density and potential in density functional theory. With our approach, convergence to exact results can be achieved by increasing the size of the embedded system. While this implies to face exponential cost with the size of the embedded system, we find that already for small embedded fragments accurate results can be reached. We illustrate our approach for molecular bond stretching and demonstrate that it reproduces the known steps and peaks that are present in the exact exchange-correlation potential.


I. INTRODUCTION
Over the past decades, Density Functional Theory (DFT) has become a well established and successful method able to accurately describe molecular and condensed matter systems. One reason for its success can be attributed to its computational efficiency as all physical observables O[n] of interest are functionals of the ground-state density n(r) instead of the many-body wave function [1]. The most popular technique to find the density of the system accurately is the Kohn-Sham (KS) DFT, where the density of the full interacting system is computed via an auxiliary non-interacting system [2]. All interactions and correlations of the interacting system are mimicked by the so-called exchange-correlation potential v xc [n] which can be approximated in many different ways [2][3][4][5][6]. A remaining challenge is to find the functionals describing the other wanted observables O[n].
Despite its efficiency, there are still open issues in DFT which need to be addressed. One of the big issues is that, once a self-consistent density has been found and an approximate functional for a certain observable O [n] has been selected, a pathway to systematically improve the obtained results is missing.
Another issue with DFT is that, although significant progress in functional development over the years has been achieved, approximate DFT functionals usually still struggle to describe systems with strongly correlated electrons [7]. The dissociation limit of the H 2 molecule is a good example for a simple system that is not easy to describe with commonly used approximate DFT functionals. The functionals that are optimized to be able to mimic the dissociation of H 2 [8][9][10][11][12][13] struggle to describe other physical systems [14].
There are alternative methods that are able to describe strongly correlated electrons accurately. One big group are wave function methods, such as full configuration interaction (FCI) methods [15] and density-matrixrenormalization group (DMRG) [16]. These methods, although becoming more and more efficient, struggle with the problem of having to store the full (approximate) wave function and thus are only able to describe relatively small systems.
A pathway to use accurate (wave-function) methods on a larger scale is provided by embedding theories. The general idea behind embedding consists of dividing a system into one or more fragments of interest and an environment, which is then considered only indirectly. With this partition the need of performing an expensive calculation on the full system is circumvented. An established group of embedding theories are various density-functional embedding methods [17][18][19][20][21][22] that have been successfully applied to a large range of complex systems. They provide ways of calculating a system which is weakly bounded to an environment by representing the environment by an external field. Opposed to that, embedding methods such as dynamical-meanfield theory (DMFT) [23][24][25], density-matrix-embedding theory (DMET) [26][27][28], and density-embedding theory (DET) [29,30] consider correlations between system and environment more explicitly and, thus, are successful in describing systems with strongly correlated electrons. In the latter two methods, only a part of the system is described accurately while the rest of the system is described with a lower level calculation. Here, the challenge is the connection between the high-level and the low-level calculation. [26,[28][29][30][31][32] All mentioned embedding methods are tailored to describe the behavior of the fragments accurately. Opposed to that, we use in the present work the embedding idea to improve our large scale description of the full system by including insights from small fragments. To this end we introduce a feedback algorithm, which combines embedding approaches with density inversions based on the one-to-one correspondence of density and potential in exact DFT. This results in a self-consistent densityfunctional embedding (SDE) technique, which allows to explicitly construct a family of novel density functionals with increasing accuracy. In our method, the global Kohn-Sham system is self-consistently optimized by including findings from accurate wave functions for different fragments. Once the optimal Kohn-Sham system is obtained, we gain information about observables from those interacting fragment wave functions.
The paper is organized as follows. In section II, we introduce the proposed SDE method step by step. In section III, we present the Hamiltonian for two electrons in a heteroatomic model system in one and two dimensions, which we use to benchmark our approach. The results for the energy and the Kohn-Sham potential of the introduced systems are shown in section IV and our findings are summarized in section V.

II. THEORY
The fundamental idea of the SDE approach is to improve the KS description of a system by including results from wave function methods, thus being able to describe strongly correlated systems. By dividing the system into small parts, which we from now on call fragments or impurities, and solving each of these parts independently, the calculation becomes feasible even for big system sizes.
The SDE method is depicted schematically in Fig. 1. It consists of the following parts, to each of which we will assign a distinct subsection: 1. The full system is described in terms of its groundstate density n(r) by means of KS DFT. The foundations of KS DFT are given in subsection II A. Figure 1. General SDE algorithm: the full system is described by its electronic density n(r) on a DFT level (left hand side, in turquoise). This system is then divided into impurities. For each impurity (orange) the system is projected onto an embedded system (right hand side). The embedded system consists of the impurity, which remains unchanged by the projection and the part of the system that includes interaction and correlation with the impurity (dark violet on right hand side). The remaining part of the system has approximately no interaction and correlation with the impurity. Each of the embedded systems is then solved on a wave-function level, yielding an accurate density and exchange-correlation potential on the impurity. These accurate local properties are then used to improve the global description of the full system. The whole process is repeated until self-consistency in the density is reached.
5. Finally, accurately calculated properties of the fragments are used to improve the KS description of the full system. This procedure is done in a selfconsistent manner and we explain the algorithm in subsection II E.
As we divide our system into fragments in real space, we will, for the sake of convenience, consider only systems that are discretized on a real space lattice throughout the paper.

A. Density Functional Theory
The mathematical foundation of ground-state DFT is the Hohenberg-Kohn theorem [1]. This theorem states, that there is a one-to-one correspondence between the external potential v ext (r) of a given system and its groundstate wave function |Ψ , as well as there is a one-to-one correspondence between the ground-state wave function of this system and its ground-state density n(r). This implies that all observables O of a many-body system are uniquely determined by its ground-state density, although the functionals O[n] themselves are in general unknown. Specifically, when trying to find a functional for the ground-state energy E[n], the kinetic energy contribution T [n] of this functional is the most difficult part and all direct approximations are so far insufficient [35]. Although difficult to approximate directly in the interacting system of interest, there is a method to indirectly approximate this term through the KS construction [2]. This construction is based on the Hohenberg-Kohn theorem, from which follows that for any interacting system there exists one and only one auxiliary non-interacting system with local multiplicative potential, which has the same density, but different ground-state wave function and external potential.
Based on this connection, the ground-state density of an interacting system is found by solving the uncoupled one-body KS equations (atomic units are used throughout the paper) The ground state of the auxiliary KS system is given by a Slater determinant built from the orbitals ϕ i (r), from which the kinetic energy (of the non-interacting system) can be obtained in a straight-forward manner. The term v S (r) is the KS potential. It includes not only the external potential of the interacting system v ext (r) but also an additional term, the Hartree-exchange-correlation potential v Hxc [n](r) that accounts for the missing interactions and the kinetic correlations.
In the KS scheme, v Hxc [n](r) is the only term which has to be approximated. Hence, the two important goals in KS DFT are to find an accurate approximation for this mapping and then to describe the observables of interest, such as the energy.
The energy and other observables can, but don't have to be obtained from the Hartree-exchange-correlation potential. Within our approach, we will provide an algorithm to find the mapping v Hxc [n](r) but calculate observables in a different manner as will be explained in section II D.

B. Continuous partition
We continue by considering the problem of dividing the full problem into fragments. Generally, the fragments have to cover the full system and should be selected small enough to be calculated with required accuracy.
In embedding approaches like subsystem DFT [22] and also in the framework of partition DFT [33], the system is divided into non-overlapping fragments, which are weakly bounded to one another. In other words, the partition is dictated by density distribution and correlations within the system and cannot be chosen arbitrar- ily. Therefore, those approaches are applicable to a small range of (weakly correlated) systems.
In DMET [26][27][28] the system is also divided into nonoverlapping fragments, which are called impurities. The partition itself can be chosen arbitrarily, as particle transfer between impurity and the rest of the system is possible within this approach. The size of the fragments is dictated mostly by the correlation length in the system [26]. The amount of correlation, which can be described with the DMET method is limited by the size of the impurity. Thus, by increasing the impurity size a convergence towards the exact solution can be achieved, which makes the method systematically improvable. However, in DMET, dividing the system into non-overlapping fragments causes artificial discontinuities in local observables [31], which we would like to avoid.
In SDE, we include correlations as in DMET (see section II C) and we also introduce a partition that guarantees that all impurities connect smoothly to one another. We define a continuous partition, where the system is covered by overlapping impurities as is depicted in Figure  2. In practice, we sweep through the system by just going one site further for each impurity calculation. When computing local observables such as the density, we only take into account the site in the middle of each impurity. Real space lattices, for which we formulate the theory, have an intrinsic discontinuity (due to discretization of the real space) and our partition procedure is constructed such that the local observables remain as continuous as they can possibly be on a real space lattice. By selecting the lattice spacing approriately, the accuracy can then be improved to reach the desired accuracy. In practice, this has to be balanced with the computational cost as for any real-space method.

C. Projection onto the embedded system
Having decided on how to divide the system into impurities, we now treat each impurity separately and find an effective description for the corresponding embedded system (see Figure 1). We want the embedded system to be such that it describes the physics on the impurity as accurately as possible. As depicted in Figure 1, we have to project the full system onto an embedded system for each impurity.
Out of a manifold of possible projections [21,23,26] we adopt here the the projection used in DMET [26][27][28] as it provides an efficient way of including static correlations between impurity and the rest of the system, which we call bath from now on.
The DMET method can be understood as a complete active space (CAS) calculation under the assumption that the impurity basis functions are always in the active space. What then remains to be found are the orbitals that build up the remaining part of the active space.
We are solving a system with the general Hamiltonian H which includes one-particle terms, accounting for the kinetic energy and the external potential of the system, and two-particle terms describing the electron-electron interaction. In the following, we will denote the number of lattice sites of our real-space basis by N and the number of impurity lattice sites by N imp .
The ground state wave function ofĤ can be exactly split up as [36] Here, the |Ã a (which in Figure 1 are depicted in orange) are the many body states located on the impurity and the|B b (which in Figure 1 in the middle are depicted in turquoise) are defined on the bath.
To project the bath states onto a basis which is optimized with respect to the interaction of the bath with the impurity, we use the Schmidt decomposition [26]. The wave function is written in terms of many-body states |A i |B i , where Σ i is a diagonal matrix containing the singular values of the decomposition. The number of states needed to describe the wave function in this new basis is 4 2·Nimp opposed to the 4 N basis functions needed to set up the full Fock space of the problem.
In order to write the exact decomposition Eq. (3) explicitly, the wave function |Ψ has to be known which Figure 3. Visualization of the decomposition of the system into impurity and bath and the projection onto embedding and environment part. The dots depict the lattice sites, which correspond to our chosen initial basis set and the crosses the orbitals after projecting. The embedding and the environment part in the system are completely decoupled. In order to describe the physics of the impurity part of the system, the computation of the embedding part is sufficient.
makes it unfeasible for realistic systems. Thus, the projection has to be obtained from an approximated wave function. In practice, the bath orbitals |B i are obtained from a non-interacting system described by a Slater determinant |Φ . Similar to the interacting system, we split up the Slater determinant into impurity and bath like in Eq. (3). For any Slater determinant, it can be shown [28] that this equation can then further be decomposed as (see Figure Here, the embedding part and the environment part do not interact with each other. Hence, in order to describe the physics on the impurity, only the embedding part of the full system has to be considered. The two contributions in the embedding part, |A i and |B l , on the other hand are fully correlated, and also particle transfer between |A i and |B l is possible. From the decomposed Kohn-Sham Slater determinant, we define a projection Which, applied to the original HamiltonianĤ yields The HamiltonianĤ emb now approximately describes the physics on the impurity and the interaction of the bath with the impurity. Its ground state can be calculated accurately with wave-function based methods. Note that, while the projection P is approximated, as it is built from a Slater determinant, we apply it to the exact HamiltonianĤ including interactions. Therefore, also the embedding HamiltonianĤ emb includes interactions. The difference between the projector in SDE and DMET is the choice of the underlying mean-field system. While in the DMET scheme, any Slater determinant can be chosen, we use the KS Slater determinant to obtain the projection. This makes our calculation of the full system in principle exact, as the KS Slater determinant reproduces the density of the interacting system and, hence, is sufficient to obtain all other observables (knowing their functionals O[n]). We will also use the one-to-one relation between ground-state density and ground-state wave function to improve the projection, which we explain in the section II E.

D. Impurity calculation
Given a projection, we obtain the embedding Hamil-tonianĤ i emb for each impurity i as described in section II C and then diagonalize it to obtain the embedding wave function |Ψ i emb of this embedded system. In the present work, we use exact diagonalization (ED) to solve for the ground-state wave function of the embedded system. We emphasize that also other solvers, such as DMRG [16,37,38], coupled cluster [39][40][41], selective CI approaches [42], or Monte-Carlo methods [43,44] can be used for the impurity calculation.
The correlated embedding wave functions can then be used to calculate the energy of the full system E or any other non-local local observable. As described in reference [27],the energy of the full system E can be approximated as a sum of impurity energies, which are calculated by taking a partial trace of the corresponding embedding density matrixρ i emb . As in the SDE approach for each impurity i only one site α i is considered for obtaining properties of the full system (see section II B), we adopt the energy formula from reference [27] to The formula above can be applied to any other observable. Thus, we circumvent the usual problem in DFT of finding explicit functional dependence O[n] between an observable of interest O and the density n by simply using the embedding wave functions instead of the density.
Before moving on to improving the KS description of the full system, we have to add an additional constrain to the impurity calculations. As in DMET [28] or partition DFT [34], we have to make sure that, when patching the system back together, we retain the correct particle number N in the full system Following reference [28], we achieve this by adding and self-consistently optimizing a chemical potential µ to the Figure 4. Visualization of the SDE algorithm: First, an initial guess for the Kohn-Sham system is made from which the projection is built. Then, for each impurity the embedding Hamiltonian is calculated and the corresponding ground-state wave function and density are computed. The density of every embedded system is then inverted and yields an updated vHxc on each site independently. This potential is then used to update the Kohn-Sham system. This is repeated until selfconsistency. An additional self-consistency cycle is added in order to maintain the correct particle number.
embedding Hamiltonian of each impuritŷ wheren α denotes the density operator on site α and the index α runs over all impurity sites. The constant µ in Eq. (9) is added only to the impurity part of the embedding Hamiltonian in order to achieve a correct particle distribution between impurity and environment. In other words, the chemical potential is a Lagrange multiplier, which assures that the constraint in Eq. (8) is fulfilled.

E. Self-consistency
A remaining question is how to improve the KS description of the full system with correlated embedding wave functions. From the embedding wave functions, we calculate the density n i emb of the embedded system i. From the Hohenberg-Kohn theorem it follows that in real space to every interacting density n i emb there belongs a local KS potentialv i S, emb =v i ext, emb +v i Hxc, emb with an identical non-interacting density.
The Hartree-exchange-correlation potentialv i Hxc for -- Figure 5. Visualization of the 1D H2 molecule. The real space is discretized on a lattice with N sites. The two atoms are modeled through a symmetric double well potential v 1D ext .
each embedded system can be obtained either by analytical [45] or numerical inversion [46][47][48], or with any robust minimization as in DMET [28]. We represent the global density and potential on each site α i by their value obtained in the corresponding embedded system i This yields the new KS Hamiltonianĥ +v Hxc , which is then used to calculate a new set of projections P i . This is done until convergence (see algorithm in Figure 4).
Eventually, we obtain an accurate density and KS potential for the considered systems as well as for observables, which are calculated with embedding wave functions as described in Eq. (7).
The SDE algorithm can be improved systematically by increasing the impurity size and it converges to the exact solution.

III. H2 IN ONE AND TWO DIMENSIONS
The SDE approach so far is valid for all closed systems that can be represented by a time independent Schrödinger equation. In order to benchmark our method and to show its efficiency, we describe the two-electron bond stretching of a heteroatomic molecule in one and two dimensions (see Fig. 5).
We model this system with the following Hamiltonian [49] on a 1D/2D real-space lattice [50] whereĉ † i,σ andĉ i,σ are the usual creation and annihilation operators of an electron with spin σ on lattice site i and n i,σ =ĉ † i,σĉ i,σ is the corresponding density operator. In 2D the index i becomes a double index with The lattice spacing ∆x is determined by the box size L in direction x and the number of lattice points N and as external potential we employ a double well potential v ext . The first part of the Hamiltonian takes into account the kinetic energy of the molecule by means of a nextneighbors hopping term. The second term in Eq. (11) is the external potential which mimics the ions of the molecule and depends on the considered dimension. In the one dimensional case, the external potential on each point is given by . The numbers z 1 and z 2 determine the depth of each well respectively. In our case they take values between 0 and 2 and we will characterize the potential by their difference ∆z = z 1 − z 2 . In the two-dimensional case the external potential takes the form accounting for both, the charge distribution of the ions in x and y direction. The third term of the Hamiltonian takes into account the interaction of the electrons. We model the electronic interaction as well as the core potentials by the soft-Coulomb interaction, which avoids the singularity at zero distance. In order to do so, we include a softening parameter α = 1.
One reason for choosing a problem that only includes two electrons is that for this example we can analytically invert the density n of the problem to yield the Hartreeexchange-correlation potential v Hxc . As the ground state of a two-electron problem is always a singlet it is valid that Inserting this property into the Kohn-Sham equations Eq. (1) yields [45] v Hxc (r) = 1 2 The formula above is given in the real space domain but it can be applied to any quantum lattice system, as there is a one-to-one correspondence between density and potential for those systems [51]. The exact inversion formula can therefore be applied to every embedded system with two electrons, hence, to every embedded system resulting from our model.

A. Dissociation of the H2 molecule
Common DFT functionals like the local density approximation (LDA [2]), or generalized gradient approximations (GGA [3,4]) fail to describe the dissociation limit of the H 2 molecule. This failure is attributed to the so-called static correlation error, which is related to fractional spin states [14]. These states arise, when we consider the parts into which the closed-shell H 2 molecule dissociates. One half of the molecule is given by an H atom with half of electron with spin up and half of electron with spin down. For the exact functional this artificial atom has the same energy as the usual H atom with one electron with either spin up or spin down [52] so that the total energy of the dissociated H 2 is twice the energy of an H atom. Common approximate functionals, however, violate this condition and predict wrong energies for fractional spin states resulting in the wrong dissociation limit.
Although there are methods such as the strictlycorrelated electron functional [9], functionals based on the random phase approximation (RPA) [8,53] and on GW combined with RPA [54], or the exchangecorrelation potential by Baerends et. al. [55,56], which were designed to overcome these issues, modeling the bond stretching of H 2 remains a challenging test for any new functional.
In Fig. 6, we show how our SDE methods performs in this test case. We plot the ground-state energy of the Hamiltonian described in section III calculated with exact diagonalization, and compare with energies obtained with one dimensional LDA-DFT [57] and with SDE. Our method quantitatively agrees with the exact solution at a much lower cost than the exact diagonalization technique. While the costs grow exponentially with the number of orbitals N in the exact case (4 N ), in the SDE method, only the costs of impurity calculations are growing exponentially with the impurity size. The costs for the underlying mean field calculation are growing quadratically with the amount of orbitals, multiplied by the amount of impurities, which is also N , and the needed self-consistency iterations η yielding a scaling of 4 2·Nimp · N 3 · η. For comparison, the numerical costs of an LDA calculation grow quadratically with the amount of lattice sites times the self-consistency iterations (N 2 · η). We can therefore state that SDE is quantitatively accurate for the present case, but with a much lower cost than the exact method.

B. Peaks and steps in the Kohn-Sham potential
For the H 2 model, the Kohn-Sham system needs to describe the repulsion of the two electrons. As the system does not include an actual interaction term, this repulsion needs to be mimicked by the Kohn-Sham potential. As has been investigated in previous work [45,58,59], we expect to see a peak that prevents the two electrons from being at the same atom. In Fig. 7, we plot the density and the Kohn-Sham potential obtained with SDE for impurity sizes of 5 and 9 sites and compare them with the exact density and the exact Kohn-Sham potential. The density from the SDE calculations for both impurity sites agrees quantitatively with the exact density. We also see peaks at position x = 0 in the Kohn-Sham potential for both SDE calculations. This peak is overestimated for small impurity sizes (N imp = 5), but agrees quantitatively with the exact solution for bigger impurity sizes (N imp = 9).
Additionally, by considering an asymmetric external potential, we can study more general models such as the bond stretching of heteroatomic molecules. In Fig. 8, also accurately agreeing with the exact solution, we observe that the electron density is not symmetric anymore. In addition to the peak observed in the symmetric case in Fig. 7, in the asymmetric case we observe a step between the two wells in the Kohn-Sham potential that accounts for the asymmetric density distribution. and FCI (dashes black line). The exact and the SDE solution agree quantitatively. The SDE Kohn-Sham potential shows the expected peak in the center which mimics the electronelectron interaction. For small impurity sizes, this behavior is overestimated, but converges quickly to a quantitatively exact result for bigger impurity sizes. The following set of parameters has been used: N = 120, L = 20, d = 10

C. Convergence behavior
The assumption made in the presented SDE method is due to the projection P which is not defined with the exact wave function but with the Kohn-Sham Slater determinant. As we do not project the impurity, this assumption is only made for the bath, i.e. for the interactions of the rest of the system with the impurity. In contrast to conventional DFT functionals, we can therefore improve our SDE functional systematically by increasing the impurity size. Once the investigated property of the system M does not change any more (with respect to the required accuracy), we know that we have an accurate estimate for this property.
In Fig. 9 and 10, we see the deviation of our results from the exact solution for different properties of the system, integrated over the whole lattice: where ∆x is the lattice constant. In Fig. 9, we plot the deviation of the density ∆n and the Kohn-Sham potential ∆v S between the SDE calculation and the exact result. We consider two different core distances (d = 0 and d = 10), which correspond to weak and strong correlation between the electrons. In both cases and for both chosen properties, we observe a decrease in ∆M with increasing impurity size up to a quan- titative agreement of the two solutions, with ∆M ≤ within the wanted accuracy. Already for the smallest considered impurity size N imp = 3, the deviations are relatively small, that is of the order of the fourth digit for the density ∆n ≤ 10 −4 and of the order of the first digit for the Kohn-Sham potential ∆v S = 10 −1 . In Fig. 9, we show the deviation of the total energy E 0 of the SDE method from the exact calculation. Again, we consider one example with weakly static correlated electrons and one example with strongly static correlated electrons. For weakly correlated electrons, the difference in energy decreases and already for an impurity size of N imp = 7, the deviation from the exact solution is below chemical accuracy of 1.6 mhartree.
For strongly correlated electrons, we observe that the SDE energy becomes smaller than the exact energy for a range of impurities between N imp = 9 and N imp = 20. This is because the SDE method is not variational. While for each single considered impurity, we create a wave function (that is variational), we do not have the interacting wave function of the whole system and calculate the energy (and every other observable) of each impurity individually. The total energy is then, due to the continuous overlapping of the impurities, the sum of all site energies. This is why the estimate can also be lower than the real energy. Also for this observable though, already for small impurities our estimate is of order ∆E 0 ≤ 10 −5 which is far below chemical accuracy.
As discussed in section II B, we do not have a wave function for the whole system, but a wave function for  . Integrated deviation of the density (upper graph) and the Kohn-Sham potential (lower graph) of the SDE calculation from the exact solution for weakly static correlated (d = 0) and strongly static correlated electrons (d = 10). In both cases, we observe a decrease in the error between the two calculations. While in the weakly correlated case the error estimate is higher for small impurities and decreases faster, in the strongly correlated case already the calculations for small impurities are very good and decrease slower. Already for (Nimp = 3), the error is of the order of ∆n ≤ 10 −4 . Parameters for d = 0: N = 120, L = 10, ∆z = 0, α = 1; parameters for d = 10: N = 120, L = 20, ∆z = 0, α = 1 every impurity. As consequence, the total particle number is not necessarily correct. The employed optimization of the chemical potential leads to the correct number for N up to a desired accuracy (| N −N | < 10 −5 ). As the energy difference is of the same order of magnitude, we further rescale the energy with respect to the particle number to achieve a better convergence behavior, as can also be seen in Fig. 10. Nonetheless, the calculated energy can still be lower than the exact energy, meaning that we still observe the non-variational nature of our approximation.

D. Application to lattice systems in 2D
It is straightforward to apply our SDE method to higher-dimensional models. To demonstrate the universality of our approach, we here discuss the H 2 molecule and a model heteroatomic molecule in two dimensions.
In Fig. 11, we plot the density n, the Kohn-Sham rescaled SDE(5) Figure 10. Difference of the total energy between the SDE and the exact solution ∆E0 with and without rescaling with respect to the particle number. We consider two different core distances (d = 0, upper graph and d = 10, lower graph), which correspond to weak and strong correlation between the electrons. For the weakly static correlated system, already for Nimp = 9, the error between the two calculations is below our selected accuracy limit. For strongly static correlated electrons, d = 10, we observe that the energy estimate of the SDE calculations for Nimp ≥ 9 is too low compared to the exact solution. The deviation in energy is very low for small impurity sizes (∆E0 ≤ 10 −5 ). Parameters for d = 0: N = 120, L = 10, ∆z = 0, α = 1; parameters for d = 10: N = 120, L = 20, ∆z = 0, α = 1 potential v S , the external potential v ext , the Hartreeexchange-correlation potential v Hxc , and deviations from the exact solution ∆n and ∆v Hxc for the two-dimensional H 2 model.
We observe a homogeneous density distribution around the two core potentials that is consistent with the external potential. The Hartree-exchange-correlation potential which mimics the interactions of the electrons as well the kinetic correlations in the interacting case, shows a peak in the middle of the molecule. Our observations are consistent with the exact solution of this problem.
For a model heteroatomic molecule, we also plot the same properties as for H 2 . The density for the heteroatomic molecule in the two-dimensional case is asymmetrically distributed between the two cores, again consistent with the external potential. In the Hartreeexchange-correlation potential, additional to the peak accounting for the interaction of the electrons, we also observe a step that accounts for the asymmetric distribution of the density.  Figure 11. The H2 molecule in two dimensions. Plotted are the density n, the Hartree-exchange-correlation potential vHxc, as well as their difference from the exact reference ∆n and ∆vHxc, respectively, the Kohn-Sham potential vS, and the external potential vext with SDE(4 × 4). We observe a homogeneous density consistent with the external potential. vHxc shows the peak accounting for the interactions of the two electrons. We observe good agreement with the exact reference. The following set of parameters has been used: Nx = 40, Ny = 20, Lx = 20, Ly = 10, d = 10, ∆z = 0

V. CONCLUSION AND OUTLOOK
We present our self-consistent density-functional embedding (SDE) approach, which is a way to systematically improve functionals for DFT. In SDE the Kohn-Sham potential is found iteratively through an embedding scheme. The observables are calculated through the chosen embedding which avoids the needs of an explicit functionalÔ[n]. SDE yields accurate results for twoelectron systems in the one-and two dimensions for moderate impurity sizes. Not only we can reproduce the exact potential energy surfaces of this system, but also the peaks and steps in the Kohn-Sham potential predicted by the exact solution. Additionally, the SDE method is systematically improvable by increasing the size of the impurity and converges to the exact solution.
The proposed SDE method is general and hence valid for larger particle numbers and higher dimensions. A wide range of solvers based on DMRG [16,37,38], coupled cluster [39][40][41], selective CI [42], or quantum Monte-Carlo [43,44] can be included into SDE to calculate larger impurity sizes and particle numbers. Further, in order to treat larger particle numbers, the analytic inversion scheme in Eq. (16) has to be substituted by a numeric one, as e.g. proposed in [46][47][48], or simply be replaced by robust optimization schemes as in conventional DMET [28]. Plotted are the density n, the Hartree-exchange-correlation potential vHxc, as well as their difference from the exact reference ∆n and ∆vHxc, respectively, the Kohn-Sham potential vS, and the external potential vext with SDE(4 ×4). We observe an asymmetric density consistent with the external potential. vHxc again shows the peak accounting for the interactions of the two electrons. Additionally, a step accounting for the asymmetric distribution of the density can be observed. Again, we observe good agreement with the exact reference. The following set of parameters has been used: Nx = 40, Ny = 20, Lx = 20, Ly = 10, d = 10, ∆z = 0. 5 We expect to face one challenge with respect to the treatment of larger systems and that is the storage and projection of the electron-electron interaction termŴ of the Hamiltonian in Eq. (11), which numerically is stored in a large tensor of fourth order (and thus is also growing by fourth order with respect to the system size). In order to treat larger systems, we either have to find an efficient way of storing the interaction tensor of the original system and then project it to the embedded system or we have to employ the non-interacting bath picture from DMET [28], that circumvents the treatment of the interaction tensor for the full system altogether.
In this work we provide a first benchmark for a promising group of methods that combine functional methods with embedding schemes, yielding systematically improvable functionals that converge to the exact ones. Work to extend the method to larger systems is underway.