Approximations based on density-matrix embedding theory for density-functional theories

Recently a novel approach to find approximate exchange-correlation functionals in density-functional theory (DFT) was presented (U. Mordovina et. al., JCTC 15, 5209 (2019)), which relies on approximations to the interacting wave function using density-matrix embedding theory (DMET). This approximate interacting wave function is constructed by using a projection determined by an iterative procedure that makes parts of the reduced density matrix of an auxiliary system the same as the approximate interacting density matrix. If only the diagonal of both systems are connected this leads to an approximation of the interacting-to-non-interacting mapping of the Kohn-Sham approach to DFT. Yet other choices are possible and allow to connect DMET with other DFTs such as kinetic-energy DFT or reduced density-matrix functional theory. In this work we give a detailed review of the basics of the DMET procedure from a DFT perspective and show how both approaches can be used to supplement each other. We do so explicitly for the case of a one-dimensional lattice system, as this is the simplest setting where we can apply DMET and the one that was originally presented. Among others we highlight how the mappings of DFTs can be used to identify uniquely defined auxiliary systems and auxiliary projections in DMET and how to construct approximations for different DFTs using DMET inspired projections. Such alternative approximation strategies become especially important for DFTs that are based on non-linearly coupled observables such as kinetic-energy DFT, where the Kohn-Sham fields are no longer simply obtainable by functional differentiation of an energy expression, or for reduced density-matrix functional theories, where a straightforward Kohn-Sham construction is not feasible.

Recently a novel approach to find approximate exchange-correlation functionals in densityfunctional theory was presented (U. Mordovina et. al., Journal of Chemical Theory and Computation 15, 5209 (2019)), which relies on approximations to the interacting wave function using density-matrix embedding theory (DMET). This approximate interacting wave function is constructed by using a projection determined by an iterative procedure that makes parts of the reduced density matrix of an auxiliary system the same as the approximate interacting density matrix. If only the diagonal of both systems are connected this leads to an approximation of the interactingto-non-interacting mapping of the Kohn-Sham approach to density-functional theory. Yet other choices are possible and allow to connect DMET with other density-functional theories such as kinetic-energy density functional theory or reduced density-matrix functional theory. In this work we give a detailed review of the basics of the DMET procedure from a density-functional perspective and show how both approaches can be used to supplement each other. We do so explicitly for the case of a one-dimensional lattice system, as this is the simplest setting where we can apply DMET and the one that was originally presented. Among others we highlight how the mappings of density-functional theories can be used to identify uniquely defined auxiliary systems and auxiliary projections in DMET and how to construct approximations for different density-functional theories using DMET inspired projections. Such alternative approximation strategies become especially important for density-functional theories that are based on non-linearly coupled observables such as kinetic-energy density-functional theory, where the Kohn-Sham fields are no longer simply obtainable by functional differentiation of an energy expression, or for reduced density-matrix functional theories, where a straightforward Kohn-Sham construction is not feasible.

I. INTRODUCTION
To ease access to readers unfamiliar with DMET we provide an extensive supplement where the many different concepts are explained with simple examples.

II. THEORETICAL SETTING
Let us, for simplicity and definiteness, choose in our considerations a finite, one-dimensional lattice system. Since we will be changing Hilbert spaces a lot in the following, let us introduce all of these spaces and how they are connected. At the same time we will also define the Hamiltonians and discuss their representations in the different spaces. Finally we will briefly discuss projections of Hamiltonians onto subspaces and some properties of the 1RDM.
A. From single-particle space to the fermionic Fock space Following the usual construction of quantum physics, we will start with the single-particle space of N sites, which we denote by with the usual inner product and ∼ = meaning isometrically isomorphic. A Hamiltonianĥ (1) on this space can be represented in the standard (site) basis |i as a hermitean N × N matrix h (1) (i, j) = i|ĥ (1) j . With the N eigenfunctions of this matrix i|φ µ = φ µ (i) and their eigenenergies ǫ µ the Hamiltonian can be equivalently represented as While we will give several explicit examples for spinless fermions in the supplement (to keep the dimensions small), in general we will consider spin 1/2 particles. All the results in the following will not depend on whether we include spin or not. The only difference lies in the dimensionalities of the objects that we consider. Since we will keep the spin dimension (a factor 2) explicit, it is usually easy to infer the spinless dimensions (else we state it explicitly). The single-particle space including spin we denote by Here the standard (site-spin) basis is denoted as |z ≡ |iσ and a HamiltonianĤ (1) can be represented as a 2N × 2N hermitean matrix that reads in eigenrepresentation So far no statistics of the particles have entered our construction. Now for the M -particle space the fermionic nature of our electrons will become important. It is common practice to construct the M -particle space in two consecutive steps. First we define the space of distinguishable particles as H M = H 1 ⊗· · ·⊗H 1 , which has dimensions (2N ) M and standard basis states of the form |z 1 ...z M ) = |z 1 ⊗· · ·⊗|z M . We want to emphasize that we denote the distinguishable-particle (non-symmetrized) basis with |·) while we later denote the indistinguishable-particle (anti-symmetrized) basis with |· . In this space the non-interacting M -particle Hamiltonian is defined asĤ (M) =Ĥ (1) ⊗1 (1) ⊗ · · · ⊗1 (1) + · · · + 1 (1) ⊗ · · · ⊗1 (1) ⊗Ĥ (1) , where1 (1) is the identity of H 1 . If we denote |φ µ1 ⊗ · · · ⊗ |φ µM = |µ 1 ...µ M ) it can be expressed asĤ (M) = 2N µ1...µM =1 (ǫ µ1 + · · · + ǫ µM )|µ 1 ...µ M )(µ 1 ...µ M |, which with the expression in the standard basis (z 1 ...z M |µ 1 ...µ M ) = φ µ1 (z 1 )...φ µM (z M ) leads to the eigenrepresentation in the spin-site basis. At this point one could wonder why we did introduce a space of distinguishable particles, when we anyway want to work with electrons? As we will show below, in quantum physics we often work explicitly in H M but restrict then the allowed states to the indistinguishable ones. Nevertheless, we can equivalently work with the Hilbert space of indistinguishable fermions, as we will also show below. Both approaches look formally similar but have some important differences, that we need to highlight for completeness and to avoid subtle errors. The first approach is straightforward. We make all anti-symmetric products for the standard basis where the sum goes over all permutations p of the M indices and σ(p) denotes whether the permutation is even (+) or odd (-). In a similar manner we can do that for any other basis, e.g., the eigenbasis of the non-interacting Hamiltonian H (M) is denoted as |µ 1 ...µ M . The number of such states is 2N M . If we now look for the eigenstate ofĤ (M) , however, restricted on this fermionic subspace, we will find all Slater determinants of the non-interacting Hamiltonian, i.e. Φ(z 1 ...z M ) = (z 1 ...z M |k 1 ...k M (6) Instead of working in the higher-dimensional space H M and then restricting the allowed states, it is also possible to work directly in the properly anti-symmetrized (fermionic) M -particle Hilbert space which is just the span of all the anti-symmetrized states. The Hamiltonian in this space can then be represented bŷ . . .
Since in the following we will work (almost) exclusively with the anti-symmetrized spaces, our Slater determinants will not have the factor 1/ √ M !. Let us next go one step further and relax the fixed number of particles restriction. To this end we construct the Fock space where the Fock-space dimension is determined by the binomial equality 2N M=0 2N M = 2 2N . In an overloading of symbols we also denote |z 1 ...z M ≡ |∅ 0 ⊕ . . . |z 1 ...z M M · · ·⊕ |∅ 2N , where |∅ is the null vector in the respective spaces and accordingly also |Φ ∈ F . The non-interacting Hamiltonian can be defined straightforwardly byĤ = 2N

M=0Ĥ
(M) F . Yet instead of this expression we would like to use creationĉ † z and annihilation operatorsĉ z , which obey the usual anti-commutation relations {ĉ z ,ĉ † z ′ } = δ zz ′ such that |z 1 ...z M =ĉ † zM ...ĉ † z1 |0 , where |0 ∈ H F 0 is the vacuum state. With these we can then define the creation and annihilation operators for the single-particle eigenstateŝ and accordingly forφ µ , which allows us to expresŝ Here the subindex s indicates in analogy to Kohn-Sham theory a non-interacting Hamiltonian. We will later see how to introduce interactions, which is the reason why a direct solution for even just the ground state becomes in practice unfeasible and we need to resort to approximations. Further, for later reference we want to introduce a basis for the Fock space F by re-labeling as follows (see Supp. SI A for an explicit example): While we did nothing intricate, this basis makes the anti-symmetry of the space implicit due to the fixed ordering of the creation operators. This implies that the Hamiltonian of Eq. (12) expressed in this basis will look quite different and the anti-symmetry of the fermionic wave functions will be carried over to the expansion coefficients (see Supp. SI A). Similar problems arise with a different construction for the Fock space, which uses local Fock spaces This allows to use a local site-spin basis |ν 1 ⊗ · · · ⊗ |ν N ∈ F ′ . Yet again, this basis is not explicitly anti-symmetrized 1 . This can also be seen by the local creationâ † iσ and annihilationâ iσ operators, which locally anti-commute, i.e., {â † iσ ,â iσ ′ } = δ σ,σ ′ , yet when extended to all of F ′ for i = i ′ actually commute, i.e., [a † iσ ,â i ′ σ ′ ] = 0. As a result, the Hamiltonian of Eq. (12) does not take the same form in terms of the local creation and annihilation operators except for special Hamiltonians like next-neighbor hopping (Hubbard) Hamiltonians. The connection follows the Jordan-Wigner transformationĉ iσ = exp(iπ σ ′ k ′ <iâ † k ′ σ ′âk ′ σ ′ )â iσ and accordingly for the creation operator. Furthermore it implies that for fermionic wave functions the expansion coefficients in this basis need to carry the missing anti-symmetry. Such an issue will appear later in our considerations when we want to express a fermionic wave function as an impurity and environment tensor product.

B. Hamiltonian restricted on Fock subspace
Let us next consider the form of the Hamiltonian of Eq. (12) restricted to a subspace of F . We will not consider just any subspace but we choose a different single-particle basis with creation operatorsφ † k and an M − 2n state |K such that we have Here we have chosen allμ ∈ {1, ..., 4n} such thatφμ|K = 0, and for the explicit exmaple in the supplement the number of basis functions 4n are 2n without spin. The subspace E is then its own Fock space of lower dimension with the new vacuum state |0 = |K . To determine the Hamiltonian on this subspace we can define a projector onto E which we denote by P E and then findĤ ′ s = P EĤs P E . We can either do so by labeling the states similarly to Eq. (13) by {|F 1 , . . . , |F 2 4n } and have a representation in an ordered basis (see Supp. SI C) or we use the representation in terms of the anti-symmetrized Fock-state basis |μ 1 . . .μ lK . In the latter case, using that we only have contributions for equal number of particles and at most oneμ =μ ′ , we find with H ′ (μ,μ ′ ) = z1,z2 H (1) (z 1 , z 2 )ϕ * µ (z 1 )ϕμ′ (z 2 ) and ∆ǫ = K |Ĥ sK Here ∆ǫN , withN = μφ † µφμ the particle number operator in E, acts as a chemical potential and takes into account the energy due to |K . Alternatively, we could have just used the identity operator on E and just added ∆ǫ1 E . If we go beyond non-interacting Hamiltonians we usually add a two-particle interaction term of the form . We first represent the interaction term in creation and annihilation operators that contain the aboveφ † 1 toφ † 4n , which leads to Here µ, ν, ξ, o go from 1 to 2N . The first 4n correspond to the ones used in E and the ones from (4n+1) to (2n+1+M ) build up |K . Next we rearrange the resultingŴ that acts on all of F in sums that go from 1 to 4n and sums that go from 4n + 1 to 2N . Since we have a fixed |K in all our states, the terms that have one index up to 4n and the other three are in (4n + 1) to (2N ) (and vice versa) are zero. The projection on E thus becomeŝ Let us note here that the terms of the projected interaction that we find here do not agree with the ones presented in, e.g., Eqs. (16) and (17) of Ref. [10].
C. Properties of the one-body reduced density matrix Let us finally comment also on some general properties of the 1RDM that will become important. For any density matrix (mixed state)ρ = l w l |Ψ l Ψ l | with l w l = 1 and |Ψ l ∈ F , the 1RDM is given by γ , where the latter expression is its diagonal representation in terms of the natural occupation numbers 0 ≤ n µ ≤ 1 and natural orbitals ψ µ (z). The diagonal provides the particle number N = 2N z=1 γ(z, z) of the density matrix. Of specific interest are here pure states in the M -particle sector of F , where one can distinguish between interacting M -particle states |Ψ with usually 0 < n µ ≤ 1 and non-interacting (Slater determinant) wave functions |Φ with n 1 = · · · = n M = 1 and the rest zero. This implies that the natural orbitals are equivalent to the orbitals of the Slater determinant, e.g., . Additionally, it also implies that a 1RDM of an interacting system cannot be reproduced by a single Slater determinant. 2

III. EXACT EMBEDDINGS VIA PROJECTIONS
The basic idea of DMET is that we divide the system into a part that we treat in detail -called the impurity -and a part that while coupled to the impurity is not treated in detail -called the environment. This division of the system into impurity and environment and the subsequent reformulation of the problem based on this division is called an embedding. While in practice the impurity is changed consecutively and the calculation is repeated such that we have treated all parts of the system in detail, the basic ingredient is the treatment of a single such impurity. In this section, where we discuss how this can be done exactly, we focus on the specific impurity A which is chosen to consist of the sites i ∈ {1, . . . , n} and the rest we denote by B. Thus the corresponding spin-site impurity is A = {1, . . . , 2n} and accordingly for B such that H 1 ∼ = A ⊕ B.

A. General embedding projections
The original (undivided) problem is usually to solve a M -particle problem on H F M with a general (usually interacting) HamiltonianĤ M F . For the DMET embedding procedure it then becomes necessary to lift this problem into Fock space. That is, we consider a hermitean Hamiltonian of the form We would then like to solve for the ground-state |Ψ in the M -particle sector. Without further simplifications this amounts to a diagonalization of a 2N M × 2N M dimensional matrix, which already for small systems becomes impossible to perform numerically exactly. We would like to reduce this prohibitively large dimensionality. To do so we assume we would know |Ψ and in a first step make the problem even more intractable by representing it in some Fock-space basis, e.g. Since belong to the impurity A and the environment B, respectively, we can re-express the ground state in a new basis Of course, since it is a M -particle problem most contributions in the full Fock space are zero (see Supp. SI B 1 for an explicit example). The expansion coefficients Ψ ij are then called the connection matrix between |F A i and |F B j . Alternatively we could also use, e.g., the local basis |ν 1 ⊗ · · · ⊗ |ν N to find such a basis for A and B, respectively.
We can then in a next step just keep those contributions that are non-zero, re-order and bring Eq. (20) in a diagonal form (see Supp. SI B 1 for an explicit example). This procedure can be done efficiently with a singular value decomposition (SVD) [22] of Ψ ij . Assuming without loss of generality n ≤ (N − n), this leads to Here, U iα and V † αj are matrix elements of unitary matrices U ∈ C 2 2n × C 2 2n and V ∈ C 2 2(N −n) × C 2 2(N −n) , and Λ αβ is a rectangular diagonal (2 2n × 2 2(N −n) )-dimensional matrix with 2n real values λ α on its diagonal. Plugging Eq. (21) into Eq. (20) then yields We have thus decomposed the ground-state wave function into the sum of tensor products of two different sets of wave functions |A α and |B α . The states |A α are defined exclusively on the impurity, while the states |B α are only defined on the environment (see Supp. SI B 2 for an explicit example). The new states |B α (which are now only 2 2n as opposed to 2 2(N −n) in (20)) are then the only ones still considered of the environment B and constitute what is called a bath for the impurity A. This construction of the impurity plus the bath is referred to in the DMET literature as the embedded system. If we next define a subspace of this embedded system span{|A α ⊗ |B β |α, β ∈ {1, . . . , 2 2n }} ∈ F , which by construction contains the M -particle ground state of interest, and define a corresponding projector we can define a 2 4n × 2 4n embedded Hamiltonian byĤ If we now restrict to only the M -particle sector and minimize the energy therein we get back the original wave function by construction. We note, however, that it is not a priori clear how many M -particle wave functions are in this subspace and it might be non-trivial to sort these wave functions (see also Supp. SI C 1 for an explicit example). Moreover, since the basis is not properly anti-symmetrized, only properly anti-symmetrized coefficients are allowed in the ensuing minimization. All of this implies that even with the exact states |A α and |B α this problem might be as hard to solve in practice as the original one if n is not chosen small enough, i.e., 2 4n ≪ 2N M .

B. Embedding projections from non-interacting systems
In the case that the Hamiltonian is non-interacting, i.e., takes the form of Eq. (12), and is non-degenerate we can express the embedded Hamiltonian in a more compact and simple form. This is due to the fact that also the ground state takes the much simpler form If we use for the (2N × M )-dimensional matrix C zµ the division and employ a SVD of the impurity submatrix Here, due to Λ being a rectangular diagonal 2n × M matrix (assuming that 2n ≤ M ) with 2n entries λ x on the diagonal, we find U · λ ∈ C 2n × C 2n (see Supp. SI B 2 for an explicit example). Note that we have overloaded the notation again by choosing the same notation for the matrices in the SVD as before in the Fock-space case. The differences (dimensions) should be obvious from the context. The rotation of orbitals that we performed implies that for µ ∈ {1, ..., 2n} the corresponding orbitalsC zµ have non-zero entries on A and B, while for µ ∈ {2n+1, ..., M } they only have non-zero entries on B. Based on these new orbitals we can introduce new creation and annihilation operators. Instead of defining such Fock-space operators for eachC zk , which would amount to M creation and annihilation operators, we define 4n + (M − 2n) by further dividing the first 2n orbitals into 2n that have non-zero values only on A and 2n that have non-zero values only on B. With the norm on A defined as φ µ A = ( 2n z=1 |φ µ (z)| 2 ) 1/2 as well as on B via φ µ B = ( 2N z=2n+1 |φ µ (z)| 2 ) 1/2 this leads to where Θ + (z) is the Heaviside step function which is 1 for z ≥ 0 and zero else. For µ > 2n we haveφ µ (z) = ϕ B µ (z) by construction. Defining with these statesφ and accordinglyφ † µ for µ > 2n we can express the non-interacting ground-state wave function as This leads to 2 2n terms which are equivalent to the ones from Eq. (22) for a non-interacting wave function. So we could now re-arrange the sum, express the different states as |A α ⊗ |B α and identify the corresponding λ α , which allows us to define a projection of the form of Eq. (23) we see that H ′ s (z 1 , z 2 ) = H s (z 1 , z 2 ) for z 1 and z 2 restricted to z ∈ A, i.e., on the impurity the Hamiltonian has the same form (see also Supp. SI C 2 for an explicit example). Restricting now to the 2n particle subspace in the Fock space E gives back the ground-state wave function of the original problem, provided we also know the form of |0 ≡ |K in terms of the original Fock space. If we do not know the form of this new vacuum state in terms of the original basis then we at least still get back the wave function on the impurity A since |K has zero contribution on A. We furthermore see that this procedure, in contrast to the one of Eq. (23), can only work in general for non-interacting problems. The reason being that an interacting wave function consists of (usually) all possible Slater determinants that we can construct and hence we cannot discard any of the original 2N orbitals and corresponding creation operators a priori .
Before we move on, let us highlight that there is a very elegant way to obtain the CAS and the corresponding matrix C CAS . If we use the previous SVD for C zµ , the 1RDM of the system can be brought into the form Using that in the sub-matrix γ env (z 1 , z 2 ) of γ(z 1 , z 2 ), with z 1 and z 2 in B ≡ {2n + 1, . . . , 2N }, only theφ µ (z) and ϕ B µ (z) from Eq. (26) contribute, we find that Thus diagonalizing γ env (z 1 , z 2 ) and only keeping those eigenfunctions ϕ B µ (z) that have eigenvalues (natural occupation numbers) 0 < n B µ = φ µ 2 B < 1 gives us directly the non-trivial entries of the matrix C CAS . See Supp. SI C 2 for an example of this construction. For later use we define here also the impurity 1RDM γ imp (z 1 , z 2 ), which is the sub-matrix of γ(z 1 , z 2 ) with z 1 and z 2 restricted to A ≡ {1, . . . , 2n}. Furthermore, we define the embedded 1RDM, which can also be found by calculating the 2n-particle ground state of the embedded Hamilto-nianĤ ′ s and excluding the orbitals of |K = |0 (also called unentangled occupied/core orbitals).

IV. MEAN-FIELD EMBEDDINGS, SELF-CONSISTENCY AND THE NON-INTERACTING v-REPRESENTABILITY ISSUE
So far we have only given some basic constituents that are part of the DMET procedure. Let us in the following connect them and discuss in more detail the fundamental algorithm of DMET. While there are many flavours available, we want to stick to the essentials and consider the standard choices where 1RDMs are matched in specific ways.
To this end will focus on matching 1RDMs locally (on each impurity).
A. Mean-field embedding via impurity one-body reduced density matrix As said before, we divide our problem in an impurity A and an environment B. To find the exact projector to perform the embedding onto A we would first need to solve the original interacting problem of the form of Eq. (18). This is of course not practical because the DMET procedure was developed to avoid exactly this unfeasible numerical task. Hence in the following we want to reduce the dimension of our problem which is 2N M . The goal is now to find an approximate projectionP . If we just use any approximate version of the form of Eq. (23) we work in a sub-space of the full Fock space with the dimension 2 4n . Already at this point we highlight that the moment we assume the size of A to be half the system, i.e., 2n = N , nothing is discarded and the projector becomes the identity, i.e., we are back in needing to solve the original problem. To find the approximate ground state (due to the approximate projection) we then need to restrict to those states that provide exactly M particles. To identify these states can be cumbersome (see also Supp. SI C for an explicit example) and hence it is desirable to have an ordering by particle number a priori. The non-interacting projections provide such an ordering, since they give rise to a new Fock space E and purpose-built Slater determinants. Hence, in practice a non-interacting projector is used. But instead of just, e.g., the projection from the Eq. (18) with W (2) ≡ 0, a self-consistency condition is enforced. Which condition and how it is enforced then connects DMET to different density-functional theories. With a non-interacting projector we therefore have the dimension 4n 2n , where we have assumed above that 2n ≤ M holds. However, if 2n > M the dimension becomes 4n M (which as one could verify gives back the original problem in the limit that 2n = N ). It is important to note here that an adaptation of how the approximate projection is determined in general would be needed for 2n > M (see discussion in Sec. IV C). A standard self-consistency condition is then where γ s imp (z 1 , z 2 ) is the 1RDM on the impurity of the auxiliary non-interacting system that provides the approximate mean-field projectorP s , and γ ′ imp (z 1 , z 2 ) is the 1RDM on the impurity of the projected interacting problem with HamiltonianP sĤPs in the respective M -particle sector. We note, however, that unless the impurity is half of the system size 2n = N (where one solves practically the original problem) it is not guaranteed that γ ′ imp (z 1 , z 2 ) and thus also the approximate interacting wave function is close to the exact γ imp (z 1 , z 2 ) and the exact interacting wave function |Ψ .
B. Non-interacting v-representability: ambiguities in the mean-field projection In order to attain self-consistency we need to define the mean-field Hamiltonian which gives the approximate projection iteratively. As we will show by the following construction there are ambiguities in this procedure as there are infinitely many non-interacting Hamiltonians that reproduce a given impurity 1RDM (see also Supp. SII B for an explicit example). We then discuss the connection of this result to the problem of the non-interacting v-representability of 1RDMs.
As an initial guess we can, e.g., solve Eq. (18) without interaction (although there are different choices). The resultingP (0) s is then used to solveP . In a next step a noninteracting system is constructed such that it reproduces the interacting 1RDM submatrix γ where we have denoted the corresponding natural occupation numbers and natural orbitals in accordance to Eq. (26). Now we only need to reverse the steps that led to Eq. (26). Firstly we choose 2n arbitrary states ϕ B µ (z) that are orthonormalized on B. Since B has a size of (2N − 2n) we have as many choices. With φ µ (where which state on A goes together with which state on B is again completely arbitrary). Since we have assumed 2n < M we have to choose (M − 2n) further arbitrary orthonormal orbitals ϕ B µ (z) (of course orthogonal to the previous 2n) and define for µ ∈ {2n + 1, . . . , M } states We have thus constructed M orthonormal single-particle statesφ µ (z) with µ ∈ {1, . . . , M } on H 1 . Since H 1 has a dimension of 2N , we are left with (2N − M ) further orthonormal states that we again order arbitrarily and denote byφ µ (z) for µ ∈ {M + 1, . . . , 2N }. As a final step we choose arbitrary energiesǫ µ ∈ R such that With these ingredients we find a single-particle Hamiltoniañ and a corresponding Fock-space Hamiltonian withφ † µ = 2N z=1φ µ (z)ĉ † z that has as its M -particle ground state imp (z 1 , z 2 ) if restricted to z 1 and z 2 ∈ A. Let us note that we have just shown that there are infinitely manyH (1) (z, z ′ ) that reproduce a given impurity 1RDM. Except of ϕ A µ (z) every other part of our construction is completely arbitrary. Yet different choices generate different projectionsP (1) s and corresponding subspaces E (1) . And if we now proceed with our iteration, each of this projector will lead to a differentP imp (z 1 , z 2 ). This is one reason why in practice the iteration might not converge. Such an ambiguity with respect to the non-interacting Hamiltonians is well known in reduced density-matrix functional theories [23,24]. It is called the non-interacting v-representability problem. It states that a non-interacting 1RDM can be generated by the ground state of many different non-interacting Hamiltonians that differ with regard to their non-local potentials v. It stems from the fact that for a non-degenerate non-interacting 1RDM only the first M orbitals are occupied. If we, however, consider a single-particle space of dimension 2N > M , the rest of the orbitals are not determined and we can thus have many Hamiltonians (see Eq. 4) that have the same non-interacting wave function as ground state. This, together with the fact that a non-degenerate non-interacting Hamiltonian cannot reproduce the 1RDM of an interacting system (see Sec. II C), prohibits usually the use of an auxiliary non-interacting system in 1RDM functional theory [23,24]. Instead one has to enforce representability conditions of the 1RDMs, which except for ensembles increase exponentially with the dimension of the single-particle space and the number of particles [25]. This will be discussed briefly also in Sec. VI.

C. Extension to the exact embedding projection
With regard to the accuracy of projecting the interacting problem with a non-interacting projector we want to highlight one specific detail. Since we solve for the ground state in the subspace E, we explicitly restrict the CAS in the M -particle sector to Slater determinants that all share the same (M − 2n) occupied orbitals. These "frozen" orbitals form |K . We expect that the thus constructed approximate interacting ground state is not very accurate if 2n is small compared to M . It is expected that for a more accurate approximation to the interacting ground state one needs to be close to 2n = M .
Of course, even in the case that 2n = M there is no guarantee that the resulting interacting ground-state wave function is well approximated. As discussed above Eq. (21), 2n ≤ N (such that the impurity is smaller or equal to the rest of the system). Only upon increasing the dimension of the complete active space to 2n = N (which corresponds to impurity being half the system size) one can guarantee to obtain the exact result. For this, however, one needs to adapt the DMET procedure in general and the projection using the CAS as described in Sec. III B is not possible anymore. Until now we have assumed that 2n ≤ M while for 2n = M all orbitals contribute to the complete active space and |0 ≡ |0 . This implies that without modifications the above procedure only works for M ≥ N , where the half-filling case 2n = M = N is still captured. Yet for M < N (which is the usual situation in quantum chemistry, since we usually approximate an infinite-dimensional problem N → ∞ by some finite value for N ) and 2n > M , we can no longer use the above introduced procedure, since we can at most define 2M orthonormal orbitals by dividing the full lattice into A and B. Hence, for 2n > M we cannot even resolve the identity on A in this way. In order to allow for an in principle exact limit of the DMET procedure with a mean-field projection for 2n > M we need to change the construction. The simplest way is to go back to the general form of the projection defined via Eq. (23). For a single Slater determinant we know from Eq. (28) that the rank of the connection matrix is at most 2M , i.e., only 2M of all the λ α are non-zero. Hence only a part of the projection onto a 2 4n -dimensional subspace of the Fock space is determined by Φ ij and the rest is arbitrary. This is why, if we want to control the rest of these dimensions by some self-consistency condition we need to work with multi-determinant mean-field wave functions. And this can only happen if the auxiliary non-interacting system is degenerate. Such a system can of course be engineered, yet becomes rather impractical and again leads to ambiguities. On the one hand, there are many multi-determinant wave functions that lead to the same impurity 1RDM as it is also the case with single determinant wave functions. By choosing the auxiliary non-interacting system to have a degenerate ground-state manifold that contains all of the necessary determinants, these wave functions can be turned into a ground state. Also, each multi-determinant wave function will lead to a different approximate projection. On the other hand, even then the rank of the connection matrix is not necessarily 2n. So there might be no clear advantage to enforce this self-consistency condition when approaching the exact projection for 2n = N .
D. Non-interacting v-representability: ambiguities in the fixed points Let us next consider the influence of the non-interacting v-representability problem on the fixed points. To do so, we employ the self-consistency condition of Eq. (33) for the special case where we apply the DMET procedure to a non-interacting reference system. While in practice not relevant, since one always solves a non-interacting system numerically exactly, it highlights potential pitfalls that arise due to the non-interacting v-representability issue. We will highlight in the following that we can find a fixed point that is an excited state of the target Hamiltonian. Still we see that the self-consistency condition of Eq. (33) is fulfilled, i.e. we have an auxiliary Hamiltonian which shares the same impurity 1RDM.
Assume that the target Hamiltonian has the form of Eq. (4) and the auxiliary Hamiltonian is given by Eq. (38).
. That is, it is not the ground state of the Hamiltonian of Eq. (4) but an excited state. Furthermore, in the construction that leads to the auxiliary Hamiltonian of Eq. (38) we choose all ϕ B µ (z) such that If N is large enough, i.e., 2N > 2n + M , this is always possible. The approximate projectorP s and its subspace E then exclude the actual ground state |Φ of the M -particle sector of the Hamiltonian of Eq. (4) and a mininmization leads to |Φ ′ and the corresponding projectionP ′ s . This implies thatP ′ sĤP ′ s andP sĤPs share the same impurity 1RDMs and the self-consistency condition of Eq. (33) is fulfilled. And instead of |Φ we find |Φ ′ at the fixed point (see also Supp. SII B for an explicit example). Realizing that we can easily construct a fixed-point solution that is even further away from |Φ by choosing the ϕ B µ (z) such that, e.g., all states φ µ (z) of |Φ do not appear in |Φ ′ (provided 2N > 2n + 2M ), the self-consistency condition does not automatically imply accuracy. We therefore do not only find multiple fixed points but also the fixed points can be far away from the exact result |Φ .
While the example is rather academic, it nicely illustrates a potential pitfall that the non-interacting v-representability poses also in the context of the DMET procedure. Here the results of density-functional theories and their mapping theorems can be potentially helpful. We will discuss this point in more detail in Sec. V. Alternatively, to overcome these ambiguities, the self-consistency condition is adapted or a global iteration is employed instead. We discuss these two options first.
E. Mean-field embedding via embedded one-body reduced density matrix The crucial problem of the self-consistency condition of Eq. (33) is that it has no unique solution due to the non-interacting v-representability problem. There are many non-interacting systems that produce a given impurity 1RDM. So it seems desirable to avoid this ambiguity. One way that is motivated by the numerical instability of the above procedure is to use the (in practice) more stable condition where γ s emb (z 1 , z 2 ) is the 1RDM of the auxiliary non-interacting system that provides the approximate mean-field projectorP s , and γ ′ emb (z 1 , z 2 ) is the 1RDM of the projected interacting problem with HamiltonianP sĤPs in the respective M -particle sector (see also Supp. SII B for an explicit example). If the full projector is used then z 1 and z 2 are defined on all of 2N . If instead, as is common practice, we build the projection using the CAS space, some of the bath orbitals (unentangled occupied/cor orbitals) ϕ µ (z) for µ ∈ {2n + 1, . . . , M } are discarded. In this case z 1 and z 2 correspond to the original lattice sites, only for z 1 and z 2 in A (see for an example the embedded 1RDM in CAS representation in Eq. (S101) and in spatial representation in Eq. (S102) and then compare with the original Eq. (S73), which is identical with the one of the full projection). 3 Which ever way we choose to determine the projection, we first note that we have to slightly modify our DMET procedure, since now also the ϕ B µ (z) are determined by the self-consistency condition of Eq. (40). The exact solution of the minimum condition of Eq. (40) is always zero. However, this leads to the impractical case of a highly degenerate non-interacting system. 4 . Restricting instead to only allow for a single Slater determinant for the case of 2n ≤ M to construct γ s emb (z 1 , z 2 ) (in which case the minimum of Eq. (40) is non-zero in general [5,7,26]) will again lead to a large ambiguity. To see this we again consider the case of a non-interacting reference system. If we choose, following the above considerations, an excited state of the reference system and construct an auxiliary Hamiltonian that has a ground-state wave function with a CAS that excludes orbitals appearing in the ground state of the reference system, we have found the minimum (γ s emb (z 1 , z 2 ) = γ ′ emb (z 1 , z 2 )). Yet this is again an undesirable fixed point. Thus this simple adaptation of the self-consistency condition is not yet enough to avoid potential problems of the non-interacting v-representability for 1RDMs.

F. Local vs global iterations
So far we have considered the situation of one impurity and investigated the ensuing self-consistency. While this is in principle enough, in practice several impurities A x with x ∈ {1, . . . , I} that together constitute the full lattice are used. This leads to yet a further large number of possible constructions and iteration procedures with different convergence criteria. It is then usually assumed that iterating locally until convergence and then step successively through all the impurities leads to the same result as when performing the iterations for all the impurities simultaneously [10].
Firstly, even though we can find for every A x potentially many auxiliary non-interacting HamiltonianH (1) x (z, z ′ ) that have the same 1RDM (from a non-degenerate ground state) on A x as the projected interacting problem, there is no procedure that somehow connects all of these auxiliary Hamiltonians and enforces that the interacting and non-interacting projected 1RDMs agree on the full lattice (for a non-degenerate ground state). The reason being, as discussed in Sec. II C, that interacting and non-interacting Hamiltonians cannot share the same 1RDM. Instead, similar to Sec. IV E, one can try to minimize the difference between the 1RDMs globally. This leads to a completely degenerate auxiliary system and in general there is no dimensional reduction. If we further enforce that we only allow for a single Slater determinant we will again find many fixed points. The reasoning is similar to the previous section. We can consider the case of two non-interacting systems on the full lattice, and can construct projectors that single out some excited state of the target system, and then build (following roughly the construction in Sec. IV D) an auxiliary system that has this state as its ground state (and generates the chosen projection). This underlines that all ambiguities due to the non-interacting v-representability that we encountered locally are also present globally. There are two main reasons for the discussed ambiguities. First, if we allow for a general non-local Hamiltonian of the form of Eq. (12), different such non-interacting Hamiltonians can have the same ground-state 1RDM. Second, unless we assume total degeneracy (which is rather impractical), a non-interacting system cannot reproduce the full 1RDM of an interacting system. These non-interacting v-representability issues are also the reason why there is no Kohn-Sham construction for 1RDM functional theory. A possible way to avoid theses ambiguities is to use the mapping theorems of density-functional theories that indicate that certain observables are representable in an interacting and a non-interacting system uniquely. For instance, instead of working with the 1RDM, we can consider only its diagonal, i.e., the (one-body spin) density. And following the usual mapping theorems we need to do this globally. In this case we can rely on the Hohenberg-Kohn mapping theorems that guarantee that there is only one auxiliary system that generates a specific density. And based on this uniqueness we have a unique auxiliary noninteracting system associated to any interacting one, at least for the global system. While this does not imply that the auxiliary projection is more accurate (for this we would need to consider the norm difference between the exact projection and the auxiliary one), we avoid the above ambiguities and can use this as a unique starting point for refinements.
The trick is thus to restrict to the density n(z) = γ(z, z) as well as the form of possible auxiliary systems. So far the auxiliary system allowed for any non-local single-particle Hamiltonian H (1) (z, z ′ ), which introduced the above discussed ambiguities. Yet to have the lattice analogue of the Hohenberg-Kohn mapping theorem we need to restrict to where we fix the hopping/kinetic term T (1) (z, z ′ ) to the one of the interacting reference system and we only allow to change the (spin-dependent) single-particle potential v(z). We note that the case of finding a projector based on the density together with the restriction to only local potentials is therefore not just a special case of the usual DMET procedure. Firstly, the basic local impurity construction of Sec. IV B is no longer possible. This is because the local potential cannot change the non-local hopping term and hence the density on the impurity depends also on (at least) the bath. So we can only follow the construction presented in Sec. IV E or directly enforce the same density globally, similar to Sec. IV F. Secondly, we avoid the major drawback of having a completely degenerate auxiliary system and do not need to enforce to only allow a single Slater determinant in the minimization. Further, the simple examples for multiple fixed points are ruled out. We (fortunately) lack the flexibility of the non-local auxiliary Hamiltonians. While the restriction to only the density n(z) has been used and discussed in the DMET literature [7], the ongoing discussion highlights that this case is special. The relation between DMET and density-based embedding theory is similar as the relation between 1RDM functional theory and density-functional theory. They are closely connected, yet call for quite different practical procedures and approximations. The use of auxiliary non-interacting systems in 1RDM functional theory is usually avoided, while in density-functional theory it is very natural and unambiguous. Similarly, the use of a non-interacting auxiliary system for the density-based procedure seems perfectly suited, while a procedure based on the 1RDM can lead to ambiguities as highlighted above. Indeed, borrowing from density-functional theory on a lattice, we know we can uniquely identify an auxiliary non-interacting systemĤ s [n] and its corresponding Kohn-Sham ground state Φ[n] (with the exact non-interacting projector P s [n]) from which we can (in principle) uniquely construct the exact interacting ground state Ψ[n] and consequently the exact projectorP [n]. And this holds irrespective of the size of the impurity. So, while in the general DMET procedure only increasing the impurity size can improve the reliability and accuracy, in the density-based embedding theory one can make the procedure exact for any impurity size. And similarly to the usual Kohn-Sham approach we can find the exact projection witĥ indication how to proceed towards interacting projections. Also, while using the general interacting projector of Eq. (23) leads to the aforementioned practical issues (symmetrization and unknown number of M -particle states), the non-interacting projection and its associated subspace E can be more practically adapted. For instance, one could aim at approximating the correlated M -particle states |A α ⊗ |B α directly from E. In this way one has direct control over symmetry and the number of particles. Besides the standard density-based functional theories there are also extensions that consider in addition to the density more complex objects, such as the current density or the kinetic-energy density. These objects are all related to parts of the full 1RDM and highlight that besides its diagonal one can potentially influence further parts of the 1RDM in an interacting as well as a non-interacting system. This then leads to new auxiliary non-interacting systems, whose auxiliary projections are potentially a better first guess to the exact projection than just connecting the density. The quantity we look at here specifically is the kinetic-energy density (for a definition of a Hubbard-type of Hamiltonian see Ref. [27] and in a continuum setting Ref. [28] (Chapter 8)). The kinetic-energy density for a Hamiltonian of the type of Eq. (18) would be where we used the decomposotion of Eq. (41). While this quantity is closely related to the 1RDM, we note that there are two main differences: (i) the T (1) (z 1 , z 2 ) (which for a Hubbard-type of Hamiltonian amounts to next-neighbour hopping term) is included and (ii) z 1 and z 2 do not take all the possible values but only the ones that appear in T (1) (z 1 , z 2 ) (for example in the standard Hubbard it will be only the next neighbors that appear). Then one only allows specific non-local potentials (of the same freedom as the interacting ones) by introducing a mean-field Hamiltonian of the type: Thus the target of the DMET procedure could be adapted so that the auxiliary system is constructed in such a way as to reproduce the density n and the kinetic-energy density K of the interacting system. The advantage of the kinetic-energy density with respect to the 1RDM is that it does not suffer from the idempotency issue, i.e. in general a non-interacting system can share the same ground-state kinetic-energy density as an interacting one [27]. However, the second question to make such a procedure well-defined is, whether the mapping between density and kinetic-energy density and local as well as non-local potential is one-to-one, i.e.

T
(1) The complication in showing that there is such a mapping lies in the fact that we consider a quantity K(z, z ′ ) that now includes the external control field T ke (z, z ′ ) as well as internal quantities, e.g., in the usual Hubbard case the first off-diagonal of the 1RDM. We therefore no longer have a simple linear structure as in density-functional theory, where external control field v(z) and the internal control objective n(z) are separate entities and are connected via a Legendre-Fenchel transformation [29,30]. This also makes the construction of approximations much more complicated. And it is for such problems, where density-functional methods can benefit strongly from the DMET procedure as we will discuss in Sec. VI. Although there is no general answer to the question whether the mapping of Eq. (45) exists, recent numerical considerations indicate that this might be the case under certain conditions [27]. Hence in analogy to the density-functional based approach, one could apply a kinetic-energy-density based approach where the exact projector isP It seems reasonable to assume that, since now the interacting and the non-interacting systems share more properties, also the zeroth order approximation to the mapping, i.e.,P Hxc [K, n] ≡ 0, is more accurate than the one from the density-functional based approach. Following this logic one can try to identify further potential mappings between the interacting and the auxiliary non-interacting system that allow to make both systems more and more alike. For instance, by including a Peierls phase in the hopping, corresponding to an external magnetic field, also the link current becomes potentially controllable [31]. Finally, there is yet a different direct way to overcome the ambiguities associated with the 1RDMs. If instead of zero temperature and definite number of particles one considers a (grand-)canonical setting, the inclusion of the entropy in the (grand-)canonical potential allows to reproduce any interacting (grand-)canonical ensemble by a unique non-interacting one [24]. The expressions for the non-interacting auxiliary system in the case of the grand-canonical situation are even analytical (see Ref. [24] in Sec. 2.4). This is in accordance to recent extensions of DMET to the (grand-)canonical setting [18].

VI. USING DENSITY-MATRIX EMBEDDING THEORY IN DENSITY-FUNCTIONAL THEORIES: A NOVEL APPROXIMATION SCHEME
Up until now we have focused on the DMET procedure and how we can understand certain subtleties connected to the non-interacting v-representability from a density-functional perspective. Having realized how the mappings of density-functional approaches appear in DMET provides us with a very interesting possibility. We can use the DMET methodology to directly approximate the interacting-to-non-interacting mapping that is the basis of the Kohn-Sham approach. Instead of indirectly connecting the interacting reference system with the auxiliary non-interacting system via the energy, we instead can directly connect a non-interacting wave function to an approximate interacting wave function that have the same target observable, e.g., in density-functional theory the density.
This idea has been realized in Ref. [5] for the standard case of density-functional theory. From the v-representability of the density in both (interacting and auxiliary non-interacting) systems we have v(z) We can of course also express this with the help of the exact embedded HamiltoniansĤ ′ [v] andĤ ′ s [v s ], respectively. This mapping directly defines the exact Hxc potential of density-functional theory by If we now approximate the embedded Hamiltonians via a self-consistent mean-field projection that makes n ′ (z) = n s (z) on the whole lattice we find the approximate mapping v(z) where the first part v(z) → n ′ (z) and the interacting wave function Ψ ′ is now only an approximation to v(z) → n(z) and Ψ. This means that we have now a new interacting mapping v → n ′ that we connect with the non-interacting system and we thus have the approximate Hxc potential Here we have indicated by v ′ that we now have a different interacting mapping (since we use the exact mean-field projection the non-interacting mapping v s is still exact) and with n ′ that we have in general also a different density at self consistency when compared to the exact density n. However, as has been demonstrated in Ref. [5], by increasing the size of the impurities A x the difference in density n − n ′ → 0. This implies that we can consistently increase the accuracy of the approximate v ′ Hxc even for strongly correlated problems. And we have access to an approximate interacting wave functions which allows to approximate many non-trivial observables that are hard to access in normal density-functional theory [32,33].
The above described procedure is a novel alternative to the usual way of obtaining density functionals. The common approach is to approximate the energy expression E[n] − E s [n] and then obtain the corresponding Hxc potential via functional derivative with respect to n(z). However, for the case of certain more complex functional variables like the kinetic-energy density K(z, z ′ ), the usual approach via the energy is no longer viable [27]. In this case the above procedure becomes instrumental to go beyond the few simple approximations known. Hence by using the approximate mappings induced by the approximate projection we find This allows to determine approximately the self-consistent effective hopping term (effective local mass) T (1) + T xc as well as the effective local potential v + v Hxc in the corresponding generalized Kohn-Sham equations with K(z, z ′ ) = M µ=1 (T (1) (z, z ′ ) + T (1) xc ([K, n]; z, z ′ ))ϕ * µ (z)ϕ µ (z ′ ) + c.c. and n(z) = M µ=1 ϕ * µ (z)ϕ µ (z). Finally, let us discuss how DMET can be used in density-matrix functional theories to find new approximation schemes. In all the above cases we do not only have access to the reduced variable under investigation, i.e., the density or the kinetic-energy density, but more importantly to an approximate interacting wave function Ψ ′ . With this we also have access to approximate interacting 1RDMs and two-body reduced density matrices (2RDMs). In this regard the DMET procedure provides a direct approximation to parts of the 2RDM and the corresponding interaction energy as well as to parts of the 1RDM and the corresponding kinetic energies. This is not as trivial as it initially sounds. In 1RDM and 2RDM functional theories it is exceedingly hard to guarantee that a trial density matrix, which is used to minimize the energy functional, corresponds to a physical interacting wave function. This crucial point has several layers of complexity attached. Firstly, while there are simple necessary and sufficient conditions known for a 1RDM to be representable by an ensemble of wave functions (ensemble N -representability) [34], for the 2RDM these conditions are infeasible in practice [35] and one hence uses rather only a subset [36]. To restrict the search space to only pure states, also for the 1RDM the conditions for pure-state N -representability become infeasible in practice [25]. Finally, if we only want to consider states that are due to the solution of an interacting Schrödinger-typ equation (v-representability) then no specific conditions are known in general. Yet using the DMET procedure we can use directly the approximate interacting 1RDMs and 2RDMs associated with Ψ ′ to minimize the energy as a functional of the respective reduced density matrices. Such an approach would suggest to adopt the usual DMET update procedure and not necessarily use a self-consistency condition. Specifically, after having made an initial guess for the auxiliary system and the corresponding auxiliary projectionP (0) s we obtain approximate 1RDMs and 2RDMs. For instance in the case of a Hubbard system with next-neighbor interaction and Hubbard on-site interaction already small impurities A x are enough to have access to all the necessary parts of the 1RDM and the 2RDM to calculate the energy. By, for instance, a downhill-simplex method [37] one then finds a modified 2RDM and the corresponding 1RDM that has lower energy. Since we have just done so by hand we are not guaranteed that it really corresponds to a wave function. If we construct a non-interacting system that shares some of the properties of this modified reduced density matrices we can use the resulting projectionP (1) s to find new physical reduced density matrices, which potentially have lower energy than the previous physical ones. In this way we can perform a minimization over v-representable 1RDMs and corresponding 2RDMs.

VII. CONCLUSION AND OUTLOOK
In this work we have highlighted how density-matrix embedding theory (DMET) and different density-functional theories can be used to supplement each other. For the simplest setting of one-dimensional finite lattices we have given a detailed review of the basics of DMET, which allowed us to directly connect this method with different density-functional-type theories. Certain ambiguities that appear in the DMET procedure could be traced back to well-known issues such as the non-interacting v-representability issue for one-body reduced density-matrix functional theory. This suggested to overcome these problems by employing appropriate mappings of density-functional theories, which guarantee unique auxiliary systems. On the other hand we could show that DMET can be used to approximate the interacting-to-non-interacting mapping fundamental to the Kohn-Sham construction directly, which provides an approximate interacting wave function from which advanced functional observables can be determined. Furthermore, the DMET procedure suggests itself as a new way to devise approximations in reduced density-matrix functional theories.
While our results are geared towards a specific setting and stay on a rather abstract level, we think that they show the potential in combining both approaches to the many-electron problem. The on-the-fly-construction of approximate interacting wave functions provides a novel paradigm in density(-matrix)-functional approximations. While in densityfunctional theories it is usually an energy expression that is approximated in terms of the functional variable or Kohn-Sham orbitals, we have seen here that DMET allows to approximate directly the interacting-to-non-interacting mapping. Considering the long and arduous history of devising more accurate density-functional approximations that also work for strongly-correlated systems this approach is promising. The approximate interacting wave functions and their reduced density matrices could also overcome in certain situations the drawback of density-matrix functional theories to enforce numerically expensive representability conditions. The main problem in both cases is of course how to treat more realistic many-electron problems in three spatial dimensions. But with the advances in the DMET procedure together with novel inversion schemes for the non-interacting mapping [38][39][40][41] it seems worthwhile to Supplemental material for: Approximations based on density-matrix embedding theory for density-functional theories

SHORT GUIDE TO THIS SUPPLEMENTAL MATERIAL
In this supplement we want to accompany the general discussion of the main text with simple, yet pedagogical examples. While the main text stays on an abstract level, we find it helpful to follow the discussion to a large part with explicit examples. This allows to focus on the essentials of the different ingredients of the DMET procedure, which for the simple systems presented in this supplement boil down to elementary matrix manipulations of relatively small matrices. Further, it allows to highlight further subtle issues, such as the proper anti-symmetrization of the physical wave function in different basis representations or that one cannot approximate the wave function of the original problem without the discarded core orbitals even on the impurity, by explicit calculations. For further convenience all of the presented results can be re-calculated with a publicly available code that can be found on GitHub.
In the following we will consider spinless fermions, i.e., the dimension of the different objects discussed here and in the main text differ by a factor of 2 in various places. In the main text we have always kept this factor explicit. This allows to directly compare with the abstract (spin-dependent) objects in the main text. In Supp. SI we start with a five-site example, give a simple non-interacting Hamiltonian and determine the three-particle ground state in different basis representations. We then present the different ways to perform the projections of the exact wave function (Supp. SI B) and of the Hamiltonian (Supp. SI C) to calculate the embedded system. We finally exemplify why for interacting systems only the projection in Fock space is applicable straightforwardly (Supp. SI D).
In Supp. SII we then consider a six-site example and then highlight first that we can find infinitely many different non-interacting Hamiltonians for a given impurity 1RDM. In Supp. SII B we then demonstrate that we can construct arbitrary fixed points of the DMET procedure if only the 1RDM on the impurity are matched.

SI. EXEMPLIFICATION OF THE DIFFERENT PROJECTIONS, SUBSPACES AND PROJECTED HAMILTONIANS
A. The basic spaces, Hamiltonians and ground-state representations Following the construction discussed in Sec. II A we first set-up the single-particle Hamiltonian. In our case of spinless particles on a five-site lattice the single-particle space is H 1 ∼ = h 1 ∼ = C 5 . The single-particle Hamiltonian includes in our case only next-neighbour hopping terms (with zero boundary conditions) and is expressed in the standard sites basis |i with i ∈ {1, . . . , 5} as (S1) Diagonalizing the matrix H (1) (i, j) we find five five-dimensional orthonormal eigenvectors φ µ (i), which allows to expressĤ (1) . Furthermore, they allow us to build the anti-symmetrized M -particle space H F M of dimension 5 M by constructing all possible M -particle Slater determinants as well as to setup the (noninteracting) M -particle Hamiltonians in this space. To be specific, for the three-particle case a Slater determinant in the non-symmetrized three-particle site basis |i, j, k) = |i ⊗ |j ⊗ |k becomes Φ(i, j, k) = (i, j, k|µ, ν, ξ This is only non-zero if µ = ν = ξ, which means we have 5 3 = 10 such wave functions. Alternatively, we can construct the three-particle Slater determinants in the non-symmetrized basis |i, j, k) of the anti-symmetrized site basis |i ′ , j ′ , k ′ as We therefore find the above Slater determinant in the anti-symmetrized basis with i < j < k as Let us next introduce the Fock space of the spinless five-site problem. If we sum over all possible Slater determinants from M = 0 to M = 5, the dimension of the resulting Fock space is 2 5 . Defining the anti-commuting creation and annihilation operators {ĉ i ,ĉ † j } = δ ij , we can create upon acting on the vacuum state |0 an orthonormal basis of 2 5 states. For instance, we haveĉ where we indicate the null vector in the respective subspace by |∅ M and we have overloaded the symbol |i, j, k as referring to the three-particle state of Eq. (S3) as well as to the three-particle Fock state of Eq. (S5). Dimensionally these two states are different since they belong to different spaces. While |i, j, k ∈ H F 3 is a 5 3 -dimensional vector with a single non-zero entry, |i, j, k ∈ F is a 2 5 -dimensional vector with a single non-zero entry.
Next we construct the many particle Hamiltonian by summing over all M -particle Hamiltonians, e.g., the threeparticle Hamiltonian reads 5 µ=1 5 ν>µ 5 ξ>ν (ǫ µ + ǫ ν + ǫ ξ )|µ, ν, ξ µ, ν, ξ|. Expressing the eigenstates in the antisymmetrized site basis |i, . . . , k we find the Fock-space Hamiltonian of Eq. (S3) in the concise form of where i, j indicates summation only over next neighbours. If we then consider the three-particle subspace, the minimal-energy solution is simply where every orbital creation operator is defined bŷ More compactly this reads as where C kµ are the overlap elements between the two different bases k|µ ≡ φ µ (k). In other words, C kµ gives the value the orbital µ has on site k, e.g., here we find The three lowest eigenvalues of (S1) are (− √ 3, -1,0), so the ground state energy of the three particle system is E = − √ 3 − 1. We note that |Φ as defined in Eq. (S9) is defined only in Fock space. Only upon projecting onto the three-particle subspace do we find a wave function of the form of Eq. (S4). Yet we in the following denote both wave functions with the same symbol |Φ since they correspond to the same physical object just represented in different spaces. Where the difference matters we will comment on it.
The issue of different spaces becomes even more clear once we choose to label the above Fock-state basis functions in a specific order similarly to Eq. (13). For instance we can define (S11) However, the basis functions that will have non-zero contributions to |Φ will be only 10, as it is in the three-particle subspace of the whole Fock space. Then the wave function |Φ in Fock space can be written as a linear combination of the basis functions |F i similarly to Eq. (19) as where we have used the prime to denote only the three-particle basis functions |F ′ i in Fock space: (S13) In this basis we can find a different expression for the Fock-space Hamiltonian of Eq. (S6), which is implicitly restricted to the three-particle subspace. Diagonalizing the resulting 10x10 matrix with matrix elements F ′ i |Ĥ|F ′ j leads to the following expansion coefficients Φ ′ i that appear in Eq. (S12): To see that this agrees with the definition of Eq. (S9), we compare with C kµ . To do so we carry out the sum and the product appearing in Eq. (S9) and the wave function is then expressed in the Fock space basis (S13). We find that (S15) which is the coefficient associated to a basis function |F ′ i , with the sites j, k, l occupied. For instance, for |Φ ′ 2 we have j=1, k=2, l=4. Carrying out this procedure for all the terms appearing in Eq. (S9) one can verify that it is the same wave function as in Eq. (S12). Here we point out that as long as we work with the creation and annihilation operators, which take into account the anti-symmetry by construction, we do not need to anti-symmetrize the coefficients C kµ . Once we have fixed a basis, e.g., |F i , the coefficients need to be anti-symmetrized, e.g., Eq. (S15). Moreover, diagonalizing (S6) in the Fock space basis gives again the same ground-state energy E as the sum of the three lowest orbital energies of (S1), i.e. E = − √ 3 − 1.

B. The different projections of the exact wave function
In the previous part of the supplement we have highlighted the connections between the single-particle, the multiparticle and the Fock-space perspectives. Now we will employ the above different representations of the same physical object, i.e., the three-particle wave function, to perform the exact projection onto an impurity subspace. To do so we choose the impurity to be A = {1, 2}, i.e., the first two sites. We then want to find a representation of the wave function such that only two orbitals have a contribution on A. This representation will then be used to define the exact projected Hamiltonian on a smaller Fock space E.

Projection via singular-value decomposition in Fock-space basis
First we show the projection of the exact wave function performed in a Fock-space basis. We will do so via a SVD in the connecting matrix that involves the Fock-space basis functions on the impurity and the corresponding ones on the environment. We can recast our wave function similarly to Eq. (20) in a way that it will comprise of basis functions |F A i that belong only to the impurity and |F B i that will belong only to the environment. The number of linearly independent |F A i that we get is 2 2 , while for |F B i is 2 3 . For this we define two new vacua |0 A and |0 B and defineĉ † 1 andĉ † 2 only on |0 A (which leads to a Fock space F A ) and accordingly for sites 3, 4 and 5 that constitute B and F B . By dimensional correspondence we find F ∼ = F A ⊗ F B . The Fock states of the respective Fock spaces are and Similar to the local creation and annhihilation operators (see Sec. II A), while the operators in each class A and B anti-commute, operators of different classes commute. So there is no automatic anti-symmetry of combined wave functions, i.e., when representing the three-particle wave function of F the coefficients need to take care of the proper symmetry. The induced basis of F then corresponds to the previously introduced basis of Eq. (S11). This means that the coefficients Φ i,j are identical with the coefficients Φ i that appear in Eq. (S12). This means that there will be only 10 non-zero entries of Φ i,j . Recasting now Eq. (S18) in terms of only its non-zero entries we find that where the matrix elements Φ ′ i are given by Eq. (S15). Having written the wave function of our example in the form of Eq. (20), we proceed in performing the SVD on the connecting matrix Φ with entries Φ i,j defined just above (all the other entries that this matrix has are zero) The connecting matrix for our example reads Rotating with U and V † , which are defined in Eq. (22) and the discussion after it, we obtain the following states on the impurity and on the bath: The wave function after the SVD reads If we compare to Eq. (S19), we see that this is of course the same wave function.

Projection via singular-value decomposition in the impurity submatrix
Here we perform the projection via a SVD for the non interacting system on the orbital coefficient matrix of the creation operators. Due to our choice of A we have which corresponds to the first two lines of the matrix defined in Eq. (S10). Performing a SVD on C A we get the factorization for its matrix elements where U (size: 2 × 2) and V (size: 3 × 3) are both orthonormal matrices and Λ is a 2 × 3 matrix with only 2 entries non-zero on the diagonal. These matrices for our example read: The matrix V is now the sought-after rotation matrix, which rotates the orbitals into a new basis. In this new basis, only the first two orbitals have overlap with the first two impurity sites. As V is a unitary matrix its determinant is one, which results in leaving the Slater determinant of the original wave function unchanged and the new "rotated" orbitals still orthonormal. Thus the impurity-projected representation of Eq. (S9) is We see in this basis that the third orbital is zero on A, i.e. the third orbital belongs purely to the environment (or it is an occupied entangled environment orbital as it is called in the DMET literature) and we hence denote it in accordance to Eq. (27) The first two orbitals have still contributions on the impurity and the environment. To construct the missing further two bath orbitals (we should have 3) we remove the part that belongs to the impurity and renormalize the resulting vectors and creation operators asφ .

(S44)
The normalization factors that appear in the denominator of Eq. (S44) read In a similar manner (see Eq. (27)) we can define the renormalization factors φ 1 A = 0.99317, φ 2 A = 0.42460 and the 2 impurity orbitals asφ If we now express Φ in these new orbitals (that are no longer normalized on the full lattice A + B but only on the respective sub-lattices) we find with the corresponding normalization coefficients Since we have now four terms in accordance to Eq. (S35), we can individually compare. Firstly we find that φ 1 A φ 2 A = 0.42170 = λ 2 (S51) and the corresponding vectors can be associated as (note the anti-symmetrization) (S53) C. The different constructions of the exact embedded system Next we construct the exact embedded system. We do so first by using the Fock-space projection and then we use the single-particle projection. While the first is the only possibility for doing the exact projection for an interacting system (see also Supp. SI D), the latter approach is the one that is employed in practice and is only exact for non-interacting systems.

Exact embedded Hamiltonian via the Fock-space projection
The Fock-space projection according to Eq. (23) reads in our casê with |A α given by Eqs. (S22)-(S25) and |B β given by Eqs. (S26)-(S29). It is instructive, however, to see how this projection looks like using the impurity plus bath orbitals (and also the environment one). Since we have already associated these states with each other in the previous section, we readily can associatê In this projection appear different terms that project to subspaces with a different number of particles. Since the Hamiltonian is particle-number conserving, contributions like ϕ A 1 ϕ B 2 ϕ B 3 |Ĥϕ A 1 ϕ B 3 are zero, yet we still have non-zero contributions within different particle-number subspaces. If we do not restrict at this point to only the three-particle subspace, a minimization of the projected Hamiltonian will lead to a different ground state with different number of particles. In this approach we therefore have to restrict by hand to those states such that the correct projector becomesP This yields a 6 × 6 Hamiltonian matrixP ′ĤP ′ (see the jupyter notebook for the explicit matrix). Diagonalizing provides its lowest eigenvalue as E = − √ 3 − 1 with the eigenstate Φ = −0.89919|A 1 ⊗ |B 1 − 0.42170|A 2 ⊗ |B 2 − 0.10566|A 3 ⊗ |B 3 − 0.04955|A 4 ⊗ |B 4 . This agrees with Eq. (S35). Alternatively we could have also restricted the projection further to only the states |A α ⊗ |B α with α ∈ {1, 2, 3, 4} (leading to a 4 × 4 Hamiltonian) and would have found the same ground state.

Exact embedded Hamiltonian via the single-particle projection
Since the above projection needs to be further restricted by hand to only the right particle sector, it is advantageous if one can directly construct the correct projector without further filtering. This can be done for non-interacting systems on the single-particle level as discussed in Sec. III B. Since one uses in practice always a non-interacting projection the following is the standard way in DMET to construct the embedded system.
We start from the non-interacting Hamiltonian of the full system, i.e., Eq. (S1) and construct the non-interacting 1RDM in the site basis from the three lowest energy orbitals (which are the ones that form the Slater determinant Φ). This will be a 5 × 5 (the number of sites) matrix As this is a non-interacting 1RDM its eigenvalues are 1 or 0. Next we consider a submatrix of this 1RDM which consists only of the sites that belong to the bath B. For our example this is the matrix defined above with i, j ∈ B ≡ {3, 4, 5} Diagonalizing this 3 × 3 submatrix gives where we have introduced the notationĉ † i |0 = |i . The orbital ϕ B 3 with occupation number 1 is called in the DMET literature unentangled occupied environmental orbital. It agrees withφ 3 ≡ ϕ B 3 from Supp. SI B 2. The two orbitals that have eigenvalues (occupations) between 0 and 1 are called the bath orbitals and agree with the corresponding ones from Sec. SI B 2.
Since they have zero contribution on the impurity A ≡ {1, 2} we need two further states (the size of the impurity) that are non-zero only on the impurity to express a 5 × 5 matrix. While we could use ϕ A 1 and ϕ A 2 from Sec. SI B 2, we can equivalently use Discarding the unentangled occupied environmental orbital ϕ B 3 that constitutes the vacuum state |0 =φ B, † 3 |0 of the Fock space E (see also Sec. II B), we are left with the 4 orbitals of the complete active space (CAS), i.e., ϕ CAS (k) takes the form of Eq. (29). With this the embedded single-particle Hamiltonian becomes While here we do not gain much in dimensionality, in the case that the impurity is much smaller than the original lattice, this leads indeed to a large reduction. The eigenvalues and eigenvectors of this single-particle embedded Hamiltonian are Lifting the single-particle Hamiltonian to the Fock-space E we have (see also Eq. (15)) Because we have discarded the unentangled occupied environmental orbital the sought-after ground state is given by the lowest two-particle eigenstate ofĤ ′ which leads to E ′ = ǫ ′ 1 + ǫ ′ 2 = −1.45179 and |Φ ′ =φ . Because in our case ∆ǫ = 0 |Ĥ0 = E − E ′ = −1.28026. For the orbitals we find If we again disregard the unentangled occupied environmental orbital then the resulting Slater determinant is While it gives the right impurity 1RDM, it does not give the wave function even on the impurity. Because the only non-trivial term isΦ ′ (1, 2) = −0.298187 we can compare to, e.g., 5 k=3Φ (1, 2, k) = 0.683012 or some arbitrary combination with k of Eq. (S2). However, if we also use that we know the discarded orbital, i.e., the form of the vacuum state |0 , we find insteadΦ (i, j, k) = (i, j, k|1, 2, 3 and have recovered the full wave function.
A. Non-uniqueness of the mean-field projection Here we show explicitly non-uniqueness of the approximate projection by constructing two non-interacting system that have the same impurity 1RDM γ s imp (i, j) but different orbitals and thus projections. We first make a random choice for a non-interacting system. Let us consider a translation invariant Hubbard system, i.e., we have only next-neighbor hopping with periodic boundary conditions: Diagonalizing it leads to six eigenstates {φ 1 , ..., φ 6 }, and choosing the three lowest eigenstates leads {φ 1 , φ 2 , φ 3 } with ground-state energy E = −4. The 1RDM To then construct a different system with the same impurity 1RDM as its three-particle ground state we first diagonalize γ s imp of Eq. (S74) and get the eigenvalues and the eigenvectors of the impurity 1RDM as Since B is four-dimensional we have four basis functions. We can choose problem adopted ones by just diagonalizing the environment 1RDM of the original γ s (i, j) of Eq. (S73) and use two of them to build our CAS space, i.e., We can construct the CAS orbitals that would be used in the auxiliary projection of the target Hamiltonian (where for the purpose of the example is not interacting but in a real application one would be interested in interacting Hamiltonians). The first two CAS orbitals can be always chosen as Alternatively, we could have used the two orbitals φ A of the impurity submatrix as we have discussed in the previous part of the supplement. Because the fourth eigenvector of the environment submatrix is discarded in the usual approximate projection (unentangled occupied/core orbital) and the first orbital is perpendicular to the subspace of the three lowest orbitals, we build the other CAS (environmental) orbitals from the remaining orbitals as While we do not need this CAS orbitals in this section, they will become important in the next. Further, they will show that we get a very different projection when we compare to the CAS from the different Hamiltonian that we construct next.
If we now take ϕ B 3 and ϕ B 4 and defineφ as well asφ Since we have now four orthogonal vectors we would still need to choose two orthonormal ones to fill up all of the six dimensional space. However, since we only want to construct a Hamiltonian that has the same impurity 1RDM in the ground-state three-particle sector we leave them undefined but instead choose a set of random numbers ǫ ′ 1 ≤ ǫ ′ 2 ≤ ǫ ′ 3 < ǫ ′ 4 ≤ ǫ ′ 5 = ǫ ′ 6 = 0. For definiteness, we choose ǫ ′ 1 = −4, ǫ ′ 2 = −3, ǫ ′ 3 = −2, ǫ ′ 4 = −1 and ǫ ′ 5 = ǫ ′ 6 = 0. With this the new Hamiltonian is and the corresponding ground-state energy is E ′ = ǫ ′ 1 + ǫ ′ 2 + ǫ ′ 3 = −9. By construction the 1RDM agrees on the impurity but the rest is different. Also, the CAS orbitals that are used in the projection will be different. In our case they become (besides the first two that are always the same) in contrast to the ones of Eqs. (S86) and (S87). However, ϕ B 4 and ϕ B 1 were orbitals that did not belong to the original CAS. That means that also the projection constructed from the hamiltonian h s new will look very different from the one of h s of Eq. (S72). Thus in this example we have highlighted that by the requirement that the impurity 1RDM is the same there is the possibility to construct completely different projections even in a very simple setting where also the target system is non-interacting and we just consider only a few sites.
B. Non-uniqueness of DMET fixed point Next we are going to demonstrate that besides the projection also the fixed point is arbitrary and that it can be arbitrarily far away from the "exact result". In our case the "exact result" is the three-particle ground state of the following "target" Hamiltonian h tar ≡ our case a Slater determinant that excludes the lowest-energy orbitals φ tar 1 ). We can engineer that by defining an auxiliary system that has the same impurity 1RDM as the excited state and a CAS that excludes the ground state of the system (in our case the CAS of Eqs. (S84) to (S87) is orthonormal to φ tar 1 ≡ ϕ B 3 ). We therefore see that by changing the eigenenergies of our auxiliary Hamiltonian in an almost arbitrary fashion as well as by choosing different eigenstates (yet still the CAS needs to be orthonormal to the lowest orbital φ tar 1 ≡ ϕ B 3 ) we can find even find many auxiliary systems that lead to this "bad" fixed point. Moreover, we can of course generate other "bad" fixed points by changing the excited state we target and construct the corresponding auxiliary systems.