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

Following the usual construction of quantum physics, we will start with the single-particle space of N sites, which we denote by

$$\begin{equation}{h}_{1}\cong {\mathbb{C}}^{N}\end{equation} \tag{ 1 }$$

with the usual inner product and ≅ meaning isometrically isomorphic. A Hamiltonian ${\hat{h}}^{(1)}$ on this space can be represented in the standard (site) basis |i⟩ as a Hermitian N × N matrix ${h}^{(1)}(i,j)=\langle i\vert {\hat{h}}^{(1)}j\rangle$. With the N eigenfunctions of this matrix ⟨i|ϕ_μ ⟩ = ϕ_μ (i) and their eigenenergies _μ the Hamiltonian can be equivalently represented as

$$\begin{equation}{h}^{(1)}(i,j)=\sum\limits _{\mu =1}^{N}{{\epsilon}}_{\mu }{\phi }_{\mu }(i){\phi }_{\mu }^{{\ast}}(j).\end{equation} \tag{ 2 }$$

While we will give several explicit examples for spinless fermions in the appendix (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

$$\begin{equation}{\mathcal{H}}_{1}={h}_{1}\otimes {\mathbb{C}}^{2}\cong {\mathbb{C}}^{2N}.\end{equation} \tag{ 3 }$$

Here the standard (site-spin) basis is denoted as |z⟩ ≡ |iσ⟩ and a Hamiltonian ${\hat{H}}^{(1)}$ can be represented as a 2N × 2N Hermitian matrix that reads in eigenrepresentation

$$\begin{equation}{H}^{(1)}(z,{z}^{\prime })=\sum\limits _{\mu =1}^{2N}{{\epsilon}}_{\mu }{\phi }_{\mu }(z){\phi }_{\mu }^{{\ast}}({z}^{\prime }).\end{equation} \tag{ 4 }$$

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 ${\mathcal{H}}_{M}={\mathcal{H}}_{1}\otimes \cdots \otimes {\mathcal{H}}_{1}$, which has dimensions (2N)^M and standard basis states of the form $\vert {z}_{1}\dots {z}_{M}\left.\right)=\vert {z}_{1}\rangle \otimes \cdots \otimes \vert {z}_{M}\rangle$. We want to emphasize that we denote the distinguishable-particle (non-symmetrized) basis with $\vert \cdot \left.\right)$ while we later denote the indistinguishable-particle (anti-symmetrized) basis with |⋅⟩. In this space the non-interacting M-particle Hamiltonian is defined as ${\hat{H}}^{(M)}={\hat{H}}^{(1)}\otimes {\hat{\mathbb{1}}}^{(1)}\otimes \cdots \otimes {\hat{\mathbb{1}}}^{(1)}+\cdots +{\hat{\mathbb{1}}}^{(1)}\otimes \cdots \otimes {\hat{\mathbb{1}}}^{(1)}\otimes {\hat{H}}^{(1)}$, where ${\hat{\mathbb{1}}}^{(1)}$ is the identity of ${\mathcal{H}}_{1}$. If we denote $\vert {\phi }_{{\mu }_{1}}\rangle \otimes \cdots \otimes \vert {\phi }_{{\mu }_{M}}\rangle =\vert {\mu }_{1}\dots {\mu }_{M}\left.\right)$ it can be expressed as ${\hat{H}}^{(M)}={\sum }_{{\mu }_{1}\dots {\mu }_{M}=1}^{2N}({{\epsilon}}_{{\mu }_{1}}+\cdots +{{\epsilon}}_{{\mu }_{M}})\vert {\mu }_{1}\dots {\mu }_{M}\left.\right)\left(\right.{\mu }_{1}\dots {\mu }_{M}\vert$, which with the expression in the standard basis $({z}_{1}\dots {z}_{M}\vert {\mu }_{1}\dots {\mu }_{M})={\phi }_{{\mu }_{1}}({z}_{1})\dots {\phi }_{{\mu }_{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 ${\mathcal{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

$$\begin{equation}\vert {z}_{1}\dots {z}_{M}\rangle =\frac{1}{\sqrt{M!}}\sum\limits _{\mathfrak{p}}\hspace{2pt}\sigma (\mathfrak{p})\vert \mathfrak{p}({z}_{1})\rangle \otimes \cdots \otimes \vert \mathfrak{p}({z}_{M})\rangle ,\end{equation} \tag{ 5 }$$

where the sum goes over all permutations $\mathfrak{p}$ of the M indices and $\sigma (\mathfrak{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 ${\hat{H}}^{(M)}$ is denoted as |μ₁...μ_M ⟩. The number of such states is $\left(\genfrac{}{}{0pt}{}{2N}{M}\right)$. If we now look for the eigenstate of ${\hat{H}}^{(M)}$, however, restricted on this fermionic subspace, we will find all Slater determinants of the non-interacting Hamiltonian, i.e.

$$\begin{align}~{{\Phi}}({z}_{1}\dots {z}_{M})& =\left({z}_{1}\dots {z}_{M}\vert {k}_{1}\dots {k}_{M}\rangle \right.\\ & =\frac{1}{\sqrt{M!}}\sum\limits _{\mathfrak{p}}\hspace{2pt}\sigma (\mathfrak{p}){\phi }_{\mathfrak{p}({\mu }_{1})}({z}_{1})\dots {\phi }_{\mathfrak{p}({\mu }_{M})}({z}_{M}).\end{align} \tag{ 6 }$$

Instead of working in the higher-dimensional space ${\mathcal{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

$$\begin{equation}{\mathcal{H}}_{M}^{\mathrm{F}}={\mathcal{H}}_{1}\wedge \cdots \wedge {\mathcal{H}}_{1}\cong {\mathbb{C}}^{\left(\genfrac{}{}{0pt}{}{2N}{M}\right)},\end{equation} \tag{ 7 }$$

which is just the span of all the anti-symmetrized states. The Hamiltonian in this space can then be represented by

$$\begin{equation}{\hat{H}}_{\mathrm{F}}^{(M)}=\sum\limits _{{\mu }_{1}=1}^{2N}\cdots \sum\limits _{{\mu }_{M}{ >}{\mu }_{M-1}}^{2N} \left({{\epsilon}}_{{\mu }_{1}}+\cdots +{{\epsilon}}_{{\mu }_{M}}\right)\vert {\mu }_{1}\dots {\mu }_{M}\rangle \langle {\mu }_{1}\dots {\mu }_{M}\vert .\end{equation} \tag{ 8 }$$

That is, in accordance to the smaller dimension the sums with respect to eigenstates are nested, i.e. μ₁ <...< μ_M . Furthermore, with respect to the anti-symmetrized spin-basis states |z₁...z_M ⟩ the Slater determinants are now

$$\begin{align}{\Phi}({z}_{1}\dots {z}_{M})& =\langle {z}_{1}\dots {z}_{M}\vert {\mu }_{1}\dots {\mu }_{M}\rangle \\ & =\sum\limits _{\mathfrak{p}}\hspace{2pt}\sigma (\mathfrak{p}){\phi }_{\mathfrak{p}({\mu }_{1})}({z}_{1})\dots {\phi }_{\mathfrak{p}({\mu }_{M})}({z}_{M}).\end{align} \tag{ 9 }$$

Since in the following we will work (almost) exclusively with the anti-symmetrized spaces, our Slater determinants will not have the factor $1/\sqrt{M!}$.

Let us next go one step further and relax the fixed number of particles restriction. To this end we construct the Fock space

$$\begin{equation}\mathcal{F}=\underset{M=0}{\overset{2N}{\oplus }}{\mathcal{H}}_{M}^{\mathrm{F}}\cong {\mathbb{C}}^{{2}^{2N}},\end{equation} \tag{ 10 }$$

where the Fock-space dimension is determined by the binomial equality ${\sum }_{M=0}^{2N}\left(\genfrac{}{}{0pt}{}{2N}{M}\right)={2}^{2N}$. In an overloading of symbols we also denote |z₁...z_M ⟩ ≡ |∅⟩₀⊕...|z₁...z_M ⟩_M ⋯ ⊕ |∅⟩_2N , where |∅⟩ is the null vector in the respective spaces and accordingly also $\vert {\Phi}\rangle \in \mathcal{F}$. The non-interacting Hamiltonian can be defined straightforwardly by $\hat{H}={\bigoplus}_{M=0}^{2N}{\hat{H}}_{\mathrm{F}}^{(M)}$. Yet instead of this expression we would like to use creation ${\hat{c}}_{z}^{{\dagger}}$ and annihilation operators ${\hat{c}}_{z}$, which obey the usual anti-commutation relations $\left\{{\hat{c}}_{z},{\hat{c}}_{{z}^{\prime }}^{{\dagger}}\right\}={\delta }_{z{z}^{\prime }}$ such that $\vert {z}_{1}\dots {z}_{M}\rangle ={\hat{c}}_{{z}_{M}}^{{\dagger}}\dots {\hat{c}}_{{z}_{1}}^{{\dagger}}\vert 0\rangle$, where $\vert 0\rangle \in {\mathcal{H}}_{0}^{\mathrm{F}}$ is the vacuum state. With these we can then define the creation and annihilation operators for the single-particle eigenstates

$$\begin{equation}{\hat{\phi }}_{\mu }^{{\dagger}}=\sum\limits _{z=1}^{2N}{\phi }_{\mu }(z){\hat{c}}_{z}^{{\dagger}},\end{equation} \tag{ 11 }$$

and accordingly for ${\hat{\phi }}_{\mu }$, which allows us to express

$$\begin{equation}{\hat{H}}_{s}=\sum\limits _{\mu =1}^{2N}{{\epsilon}}_{\mu }{\hat{\phi }}_{\mu }^{{\dagger}}{\hat{\phi }}_{\mu }=\sum\limits _{z,{z}^{\prime }=1}^{2N}{H}^{(1)}(z,{z}^{\prime }){\hat{c}}_{z}^{{\dagger}}{\hat{c}}_{{z}^{\prime }}.\end{equation} \tag{ 12 }$$

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 $\mathcal{F}$ by re-labeling as follows (see appendix A.1.1 for an explicit example):

$$\begin{align}\begin{aligned}\vert {F}_{1}\rangle & =\vert 0\rangle r\\ \vert {F}_{2}\rangle & ={\hat{c}}_{1{\uparrow}}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{3}\rangle & ={\hat{c}}_{1{\downarrow}}^{{\dagger}}\vert 0\rangle \\ \dots \\ \vert {F}_{{2}^{2N}}\rangle & ={\hat{c}}_{1{\uparrow}}^{{\dagger}}\dots {\hat{c}}_{N{\downarrow}}^{{\dagger}}\vert 0\rangle .\end{aligned}\end{align} \tag{ 13 }$$

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 equation (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 appendix A.1.1). Similar problems arise with a different construction for the Fock space, which uses local Fock spaces ${\mathcal{F}}_{i}\cong {\mathbb{C}}^{4}$, i.e. ${\mathcal{F}}_{i}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert 0{\rangle }_{i},\vert {\uparrow}{\rangle }_{i},\vert {\downarrow}{\rangle }_{i},\vert {\uparrow}{\downarrow}{\rangle }_{i}\right\}$, such that ${\mathcal{F}}^{\prime }={\bigotimes}_{i=1}^{N}{\mathcal{F}}_{i}\cong \mathcal{F}$. This allows to use a local site-spin basis $\vert {\nu }_{1}\rangle \otimes \cdots \otimes \vert {\nu }_{N}\rangle \in {\mathcal{F}}^{\prime }$. Yet again, this basis is not explicitly anti-symmetrized ⁴ . This can also be seen by the local creation ${\hat{a}}_{i\sigma }^{{\dagger}}$ and annihilation ${\hat{a}}_{i\sigma }$ operators, which locally anti-commute, i.e. $\left\{{\hat{a}}_{i\sigma }^{{\dagger}},{\hat{a}}_{i{\sigma }^{\prime }}\right\}={\delta }_{\sigma ,{\sigma }^{\prime }}$, yet when extended to all of ${\mathcal{F}}^{\prime }$ for i ≠ i' actually commute, i.e. $[{a}_{i\sigma }^{{\dagger}},{\hat{a}}_{{i}^{\prime }{\sigma }^{\prime }}]=0$. As a result, the Hamiltonian of equation (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 ${\hat{c}}_{i\sigma }=\mathrm{exp}(\mathrm{i}\pi {\sum }_{{\sigma }^{\prime }}{\sum }_{{k}^{\prime }{< }i}{\hat{a}}_{{k}^{\prime }{\sigma }^{\prime }}^{{\dagger}}{\hat{a}}_{{k}^{\prime }{\sigma }^{\prime }})\enspace {\hat{a}}_{i\sigma }$ 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.

Let us next consider the form of the Hamiltonian of equation (12) restricted to a subspace of $\mathcal{F}$. We will not consider just any subspace but we choose a different single-particle basis with creation operators ${\hat{\varphi }}_{~{k}}^{{\dagger}}$ and an M − 2n state $\vert ~{K}\rangle$ such that we have

$$\begin{equation}\mathcal{E}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert ~{K}\rangle ,{\hat{\varphi }}_{1}^{{\dagger}}\vert ~{K}\rangle ,\dots ,{\hat{\varphi }}_{4n}^{{\dagger}}\dots {\hat{\varphi }}_{1}^{{\dagger}}\vert ~{K}\rangle \right\}\cong {\mathbb{C}}^{{2}^{4n}}.\end{equation} \tag{ 14 }$$

Here we have chosen all $~{\mu }\in \left\{1,\dots ,4n\right\}$ such that ${\hat{\varphi }}_{~{\mu }}\vert ~{K}\rangle =0$, and for the explicit example in the appendix the number of basis functions 4n are 2n without spin. The subspace $\mathcal{E}$ is then its own Fock space of lower dimension with the new vacuum state $\vert ~{0}\rangle =\vert ~{K}\rangle$. To determine the Hamiltonian on this subspace we can define a projector onto $\mathcal{E}$ which we denote by ${P}_{\mathcal{E}}$ and then find ${\hat{H}}_{s}^{\prime }={P}_{\mathcal{E}}{\hat{H}}_{s}{P}_{\mathcal{E}}$. We can either do so by labeling the states similarly to equation (13) by $\left\{\vert {~{F}}_{1}\rangle ,\dots ,\vert {~{F}}_{{2}^{4n}}\rangle \right\}$ and have a representation in an ordered basis (see appendix A.3) or we use the representation in terms of the anti-symmetrized Fock-state basis $\vert {~{\mu }}_{1}\dots {~{\mu }}_{l}~{K}\rangle$. In the latter case, using that we only have contributions for equal number of particles and at most one $~{\mu }\ne {~{\mu }}^{\prime }$, we find with ${H}^{\prime }(~{\mu },{~{\mu }}^{\prime })={\sum }_{{z}_{1},{z}_{2}}{H}^{(1)}({z}_{1},{z}_{2}){\varphi }_{~{\mu }}^{{\ast}}({z}_{1}){\varphi }_{{~{\mu }}^{\prime }}({z}_{2})$ and ${\Delta}{\epsilon}=\langle ~{K}\vert {\hat{H}}_{s}~{K}\rangle$

$$\begin{equation}{\hat{H}}_{s}^{\prime }=\sum\limits _{~{\mu },{~{\mu }}^{\prime }=1}^{4n}{H}^{\prime }(~{\mu },{~{\mu }}^{\prime }){\hat{\varphi }}_{~{\mu }}^{{\dagger}}{\hat{\varphi }}_{{~{\mu }}^{\prime }}+\frac{{\Delta}{\epsilon}}{2n}\hat{~{N}}.\end{equation} \tag{ 15 }$$

Here ${\Delta}{\epsilon}\hat{~{N}}$, with $\hat{~{N}}={\sum }_{~{\mu }}\hspace{2pt}{\hat{\varphi }}_{~{\mu }}^{{\dagger}}{\hat{\varphi }}_{~{\mu }}$ the particle number operator in $\mathcal{E}$, acts as a chemical potential and takes into account the energy due to $\vert ~{K}\rangle$. Alternatively, we could have just used the identity operator on $\mathcal{E}$ and just added ${\Delta}{\epsilon}{\hat{\mathbb{1}}}^{\mathcal{E}}$. If we go beyond non-interacting Hamiltonians we usually add a two-particle interaction term of the form $\hat{W}={\sum }_{{z}_{1},{z}_{2},{z}_{2},{z}_{1}}{W}^{(2)}({z}_{1},{z}_{2},{z}_{3},{z}_{4}){\hat{c}}_{{z}_{1}}^{{\dagger}}{\hat{c}}_{{z}_{2}}^{{\dagger}}{\hat{c}}_{{z}_{2}}{\hat{c}}_{{z}_{1}}$. We first represent the interaction term in creation and annihilation operators that contain the above ${\hat{\varphi }}_{1}^{{\dagger}}$ to ${\hat{\varphi }}_{4n}^{{\dagger}}$, which leads to

$$\begin{align}& {W}^{(2)}(\mu ,\nu ,\xi ,o)\\ & \hspace{25.0pt}=\sum\limits _{{z}_{1},{z}_{2}=1}^{2N}{\varphi }_{\mu }^{{\ast}}({z}_{1}){\varphi }_{\nu }^{{\ast}}({z}_{2}){W}^{(2)}({z}_{1},{z}_{2},{z}_{2},{z}_{1}){\varphi }_{\xi }({z}_{2}){\varphi }_{o}({z}_{1}).\end{align} \tag{ 16 }$$

Here μ, ν, ξ, o go from 1 to 2N. The first 4n correspond to the ones used in $\mathcal{E}$ and the ones from (4n + 1) to (2n + 1 + M) build up $\vert ~{K}\rangle$. Next we rearrange the resulting $\hat{W}$ that acts on all of $\mathcal{F}$ in sums that go from 1 to 4n and sums that go from 4n + 1 to 2N. Since we have a fixed $\vert ~{K}\rangle$ 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 $\mathcal{E}$ thus becomes

$$\begin{align}{\hat{W}}^{\prime }& =\sum\limits _{~{\mu },~{\nu },~{\xi },~{o}=1}^{4n}{W}^{(2)}(~{\mu },~{\nu },~{\xi },~{o}){\hat{\varphi }}_{~{\mu }}^{{\dagger}}{\hat{\varphi }}_{~{\nu }}^{{\dagger}}{\hat{\varphi }}_{~{\xi }}{\hat{\varphi }}_{~{o}}+\frac{\langle ~{K}\vert \hat{W}~{K}\rangle }{2n}\hat{~{N}}\\ & \quad +\sum\limits _{~{\mu },~{\xi }=1}^{4n}\left[\sum\limits _{\nu =4n+1}^{2n+1+M}\left({W}^{(2)}(~{\mu },\nu ,\nu ,~{\xi })-{W}^{(2)}(~{\mu },\nu ,~{\xi },\nu )\right.\left.+{W}^{(2)}(\nu ,~{\mu },~{\xi },\nu )-{W}^{(2)}(\nu ,~{\mu },\nu ,~{\xi })\right)\right]{\hat{\varphi }}_{~{\mu }}^{{\dagger}}{\hat{\varphi }}_{~{\xi }}.\end{align} \tag{ 17 }$$

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. equations (16) and (17) of reference [37].

Let us finally comment also on some general properties of the 1RDM that will become important. For any density matrix (mixed state) $\hat{\rho }={\sum }_{l}\hspace{2pt}{w}_{l}\vert {{\Psi}}_{l}\rangle \langle {{\Psi}}_{l}\vert$ with ∑_l w_l = 1 and $\vert {{\Psi}}_{l}\rangle \in \mathcal{F}$, the 1RDM is given by $\gamma ({z}_{1},{z}_{2})=\mathrm{Tr}(\hat{\rho }{\hat{c}}_{{z}_{1}}^{{\dagger}}{\hat{c}}_{{z}_{2}})={\sum }_{\mu =1}^{2N}{n}_{\mu }{\psi }_{\mu }^{{\ast}}({z}_{1}){\psi }_{\mu }({z}_{2})$, 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={\sum }_{z=1}^{2N}\gamma (z,z)$ of the density matrix. Of specific interest are here pure states in the M-particle sector of $\mathcal{F}$, where one can distinguish between interacting M-particle states |Ψ⟩ with usually 0 < n_μ ⩽ 1 and non-interacting (Slater determinant) wave functions |Φ⟩ with n₁ = ⋯ = n_M = 1 and the rest zero. This implies that the natural orbitals are equivalent to the orbitals of the Slater determinant, e.g. $\langle {\mu }_{1}\dots {\mu }_{M}\vert {\hat{c}}_{{z}_{1}}^{{\dagger}}{\hat{c}}_{{z}_{2}}{\mu }_{1}\dots {\mu }_{M}\rangle ={\sum }_{i=1}^{M}{\phi }_{{\mu }_{i}}^{{\ast}}({z}_{1}){\phi }_{{\mu }_{i}}({z}_{2})$. Additionally, it also implies that a 1RDM of an interacting system cannot be reproduced by a single Slater determinant ⁵ .

The original (undivided) problem is usually to solve an M-particle problem on ${\mathcal{H}}_{M}^{\mathrm{F}}$ with a general (usually interacting) Hamiltonian ${\hat{H}}_{\mathrm{F}}^{M}$. For the DMET embedding procedure it then becomes necessary to lift this problem into Fock space. That is, we consider a Hermitian Hamiltonian of the form

$$\begin{align}\hat{H}& =\sum\limits _{{z}_{1},{z}_{2}}\hspace{2pt}{H}^{(1)}({z}_{1},{z}_{2}){\hat{c}}_{{z}_{1}}^{{\dagger}}{\hat{c}}_{{z}_{2}}\\ & \quad +\sum\limits _{{z}_{1},{z}_{2}}{W}^{(2)}({z}_{1},{z}_{2},{z}_{2},{z}_{1}){\hat{c}}_{{z}_{1}}^{{\dagger}}{\hat{c}}_{{z}_{2}}^{{\dagger}}{\hat{c}}_{{z}_{2}}{\hat{c}}_{{z}_{1}}.\end{align} \tag{ 18 }$$

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 $\left(\genfrac{}{}{0pt}{}{2N}{M}\right){\times}\left(\genfrac{}{}{0pt}{}{2N}{M}\right)$ 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.

$$\begin{equation}\vert {\Psi}\rangle =\sum\limits _{i=1}^{{2}^{2N}}{{\Psi}}_{i}\vert {F}_{i}\rangle .\end{equation} \tag{ 19 }$$

Since each $\vert {F}_{i}\rangle =\vert {F}_{i}^{A}\rangle \otimes \vert {F}_{i}^{B}\rangle$, where $\vert {F}_{i}^{A}\rangle \in {\mathcal{F}}_{A}\cong {\mathbb{C}}^{{2}^{2n}}$ and $\vert {F}_{i}^{B}\rangle \in {\mathcal{F}}_{B}\cong {\mathbb{C}}^{{2}^{2(N-n)}}$ belong to the impurity A and the environment B, respectively, we can re-express the ground state in a new basis

$$\begin{equation}\vert {\Psi}\rangle =\sum\limits _{i=1}^{{2}^{2n}}\sum\limits _{j=1}^{{2}^{2(N-n)}}{{\Psi}}_{ij}\vert {F}_{i}^{A}\rangle \otimes \vert {F}_{j}^{B}\rangle .\end{equation} \tag{ 20 }$$

Of course, since it is an M-particle problem most contributions in the full Fock space are zero (see appendix A.2.1 for an explicit example). The expansion coefficients Ψ_ij are then called the connection matrix between $\vert {F}_{i}^{A}\rangle$ and $\vert {F}_{j}^{B}\rangle$. Alternatively we could also use, e.g. the local basis |ν₁⟩ ⊗ ⋯ ⊗ |ν_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 equation (20) in a diagonal form (see appendix A.2.1 for an explicit example). This procedure can be done efficiently with a singular value decomposition (SVD) [12] of Ψ_ij . Assuming without loss of generality n ⩽ (N − n), this leads to

$$\begin{equation}{{\Psi}}_{ij}=\sum\limits _{\alpha =1}^{{2}^{2n}}\sum\limits _{\beta =1}^{{2}^{2(N-n)}}{U}_{i\alpha }{{\Lambda}}_{\alpha \beta }{V}_{\beta j}^{{\dagger}}.\end{equation} \tag{ 21 }$$

Here, U_iα and ${V}_{\beta j}^{{\dagger}}$ are matrix elements of unitary matrices $U\in {\mathbb{C}}^{{2}^{2n}}{\times}{\mathbb{C}}^{{2}^{2n}}$ and $V\in {\mathbb{C}}^{{2}^{2(N-n)}}{\times}{\mathbb{C}}^{{2}^{2(N-n)}}$, and Λ_αβ is a rectangular diagonal (2²ⁿ × 2^2(N−n))-dimensional matrix with 2n real values λ_α on its diagonal. Plugging equation (21) into equation (20) then yields

$$\begin{align}\vert {\Psi}\rangle & =\sum\limits _{i=1}^{{2}^{2n}}\sum\limits _{j=1}^{{2}^{2(N-n)}}\sum\limits _{\alpha =1}^{{2}^{2n}}{U}_{i\alpha }{\lambda }_{\alpha }{V}_{\alpha j}^{{\dagger}}\vert {F}_{i}^{A}\rangle \otimes \vert {F}_{j}^{B}\rangle ,\\ & =\sum\limits _{\alpha =1}^{{2}^{2n}}{\lambda }_{\alpha }\underset{=\vert {A}_{\alpha }\rangle }{\underbrace{\sum\limits _{i=1}^{{2}^{2n}}{U}_{i\alpha }\vert {F}_{i}^{A}\rangle }}\otimes \underset{=\vert {B}_{\alpha }\rangle }{\underbrace{\sum\limits _{j=1}^{{2}^{2(N-n)}}{V}_{\alpha j}^{{\dagger}}\vert {F}_{j}^{B}\rangle }},\\ & =\sum\limits _{\alpha =1}^{{2}^{2n}}{\lambda }_{\alpha }\vert {A}_{\alpha }\rangle \otimes \vert {B}_{\alpha }\rangle .\end{align} \tag{ 22 }$$

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 appendix A.2.2 for an explicit example). The new states |B_α ⟩ (which are now only 2²ⁿ 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 $\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert {A}_{\alpha }\rangle \otimes \vert {B}_{\beta }\rangle \vert \alpha ,\beta \in \left\{1,\dots ,{2}^{2n}\right\}\right\}\in \mathcal{F}$, which by construction contains the M-particle ground state of interest, and define a corresponding projector

$$\begin{equation}\hat{P}=\sum\limits _{\alpha ,\beta =1}^{{2}^{2n}}\vert {A}_{\alpha }\rangle \otimes \vert {B}_{\beta }\rangle \langle {A}_{\alpha }\vert \otimes \langle {B}_{\beta }\vert \end{equation} \tag{ 23 }$$

we can define a 2⁴ⁿ × 2⁴ⁿ embedded Hamiltonian by

$$\begin{equation}{\hat{H}}^{\prime }=\hat{P}\hat{H}\hat{P}.\end{equation} \tag{ 24 }$$

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 appendix A.3.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}\ll \left(\genfrac{}{}{0pt}{}{2N}{M}\right)$.

In the case that the Hamiltonian is non-interacting, i.e. takes the form of equation (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

$$\begin{equation}\vert {\Phi}\rangle =\left(\sum\limits _{z=1}^{2N}\underset{={C}_{z1}}{\underbrace{{\phi }_{1}(z)}}{\hat{c}}_{z}^{{\dagger}}\right)\dots \left(\sum\limits _{z=1}^{2N}\underset{={C}_{zM}}{\underbrace{{\phi }_{M}(z)}}{\hat{c}}_{z}^{{\dagger}}\right)\vert 0\rangle .\end{equation} \tag{ 25 }$$

If we use for the (2N × M)-dimensional matrix C_zμ the division , and employ an SVD of the impurity submatrix ${C}_{z\mu }^{A}={\sum }_{x=1}^{2n}{\sum }_{l=1}^{M}{U}_{zx}^{A}{{\Lambda}}_{xl}{V}_{l\mu }^{A{\dagger}}$, where ${U}_{A}\in {\mathbb{C}}^{2n}{\times}{\mathbb{C}}^{2n}$ and ${V}^{A}\in {\mathbb{C}}^{M}{\times}{\mathbb{C}}^{M}$, we find

Here, due to Λ being a rectangular diagonal 2n × M matrix (assuming that 2n ⩽ M) with 2n entries λ_x on the diagonal, we find $U\cdot \lambda \in {\mathbb{C}}^{2n}{\times}{\mathbb{C}}^{2n}$ (see appendix A.2.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 orbitals ${~{C}}_{z\mu }$ 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 each ${~{C}}_{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 ${\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{A}={({\sum }_{z=1}^{2n}\vert {~{\varphi }}_{\mu }(z){\vert }^{2})}^{1/2}$ as well as on B via ${\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{B}={({\sum }_{z=2n+1}^{2N}\vert {~{\varphi }}_{\mu }(z){\vert }^{2})}^{1/2}$ this leads to

$$\begin{align}\begin{aligned}{\varphi }_{\mu }^{A}(z)& =\frac{1}{{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{A}}{~{\varphi }}_{\mu }(z){{\Theta}}_{+}(2n-z),\\ {\varphi }_{\mu }^{B}(z)& =\frac{1}{{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{B}}{~{\varphi }}_{\mu }(z){{\Theta}}_{+}(z-2n-1),\end{aligned}\end{align} \tag{ 26 }$$

where Θ₊(z) is the Heaviside step function which is 1 for z ⩾ 0 and zero else. For μ > 2n we have ${~{\varphi }}_{\mu }(z)={\varphi }_{\mu }^{B}(z)$ by construction. Defining with these states

$$\begin{align}\begin{aligned}{\hat{\varphi }}_{\mu }^{A,{\dagger}}=\sum\limits _{z=1}^{2N}{\varphi }_{\mu }^{A}(z){\hat{c}}^{{\dagger}}(z)\\ {\hat{\varphi }}_{\mu }^{B,{\dagger}}=\sum\limits _{z=1}^{2N}{\varphi }_{\mu }^{B}(z){\hat{c}}^{{\dagger}}(z)\end{aligned}\end{align} \tag{ 27 }$$

and accordingly ${\hat{\varphi }}_{\mu }^{{\dagger}}$ for μ > 2n we can express the non-interacting ground-state wave function as

$$\begin{equation}\vert {\Phi}\rangle =\left({\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\hat{\varphi }}_{1}^{A,{\dagger}}+{\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}{\hat{\varphi }}_{1}^{B,{\dagger}}\right)\dots \left({\Vert}{~{\varphi }}_{2n}{{\Vert}}_{A}{\hat{\varphi }}_{2n}^{A,{\dagger}}+{\Vert}{~{\varphi }}_{2n}{{\Vert}}_{B}{\hat{\varphi }}_{2n}^{B,{\dagger}}\right){\hat{\varphi }}_{2n+1}^{{\dagger}}\dots {\hat{\varphi }}_{M}^{{\dagger}}\vert 0\rangle .\end{equation} \tag{ 28 }$$

This leads to 2²ⁿ terms which are equivalent to the ones from equation (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 equation (23) (see appendix A.2.2 for an explicit example). Instead we use that the above defined creation operators of equation (27) span a subspace of the form of equation (14) and the projection thus leads to a Hamiltonian of the form of equation (15). Since the new non-interacting Hamiltonian can be determined as a 4n × 4n matrix of the form ${H}_{s}^{\prime }={C}_{\text{CAS}}^{{\dagger}}{H}_{s}{C}_{\text{CAS}}$ with the complete active space (CAS) matrix

we see that ${H}_{s}^{\prime }({z}_{1},{z}_{2})={H}_{s}({z}_{1},{z}_{2})$ for z₁ and z₂ restricted to z ∈ A, i.e. on the impurity the Hamiltonian has the same form (see also appendix A.3.2 for an explicit example). Restricting now to the 2n particle subspace in the Fock space $\mathcal{E}$ gives back the ground-state wave function of the original problem, provided we also know the form of $\vert ~{0}\rangle \equiv \vert ~{K}\rangle$ 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 1RDM on the impurity A since $\vert ~{K}\rangle$ has zero contribution on A. We furthermore see that this procedure, in contrast to the one of equation (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

$$\begin{align}\gamma ({z}_{1},{z}_{2})& =\sum\limits _{k=1}^{M}{C}_{{z}_{1}k}^{{\ast}}{C}_{k{z}_{2}}=\sum\limits _{\mu =1}^{M}{~{C}}_{{z}_{1}\mu }^{{\ast}}{~{C}}_{\mu {z}_{2}}\\ & =\sum\limits _{\mu =1}^{2n}{~{\varphi }}_{\mu }^{{\ast}}({z}_{1}){~{\varphi }}_{\mu }({z}_{2})+\sum\limits _{\mu =2n+1}^{M}{~{\varphi }}_{\mu }^{{\ast}}({z}_{1}){~{\varphi }}_{\mu }({z}_{2}).\end{align} \tag{ 30 }$$

Using that in the sub-matrix γ_env(z₁, z₂) of γ(z₁, z₂), with z₁ and z₂ in B ≡ {2n + 1, ..., 2N}, only the ${~{\varphi }}_{\mu }(z)$ and ${\varphi }_{\mu }^{B}(z)$ from equation (26) contribute, we find that

$$\begin{align}{\gamma }_{\text{env}}({z}_{1},{z}_{2})& =\sum\limits _{\mu =1}^{2n}{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{B}^{2}{\varphi }_{\mu }^{B,{\ast}}({z}_{1}){\varphi }_{\mu }^{B}({z}_{2})\\ & \quad +\sum\limits _{\mu =2n+1}^{M}{\varphi }_{\mu }^{B,{\ast}}({z}_{1}){\varphi }_{\mu }^{B}({z}_{2}).\end{align} \tag{ 31 }$$

Thus diagonalizing γ_env(z₁, z₂) and only keeping those eigenfunctions ${\varphi }_{\mu }^{B}(z)$ that have eigenvalues (natural occupation numbers) $0{< }{n}_{\mu }^{B}={\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{B}^{2}{< }1$ gives us directly the non-trivial entries of the matrix C_CAS. See appendix A.3.2 for an example of this construction. For later use we define here also the impurity 1RDM γ_imp(z₁, z₂), which is the sub-matrix of γ(z₁, z₂) with z₁ and z₂ restricted to A ≡ {1, ..., 2n}. Furthermore, we define

$$\begin{align}{\gamma }_{\text{emb}}({z}_{1},{z}_{2})& =\underset{={\gamma }_{\text{imp}}({z}_{1},{z}_{2})}{\underbrace{\sum\limits _{\mu =1}^{2n}{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{A}^{2}{\varphi }_{\mu }^{A,{\ast}}({z}_{1}){\varphi }^{A}({z}_{2})}}\\ & \quad +\sum\limits _{\mu =1}^{2n}{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{B}^{2}{\varphi }_{\mu }^{B,{\ast}}({z}_{1}){\varphi }^{B}({z}_{2})\end{align} \tag{ 32 }$$

the embedded 1RDM, which can also be found by calculating the 2n-particle ground state of the embedded Hamiltonian ${\hat{H}}_{s}^{\prime }$ and excluding the orbitals of $\vert ~{K}\rangle =\vert ~{0}\rangle$ (also called unentangled occupied/core orbitals).

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 equation (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 $\left(\genfrac{}{}{0pt}{}{2N}{M}\right)$. The goal is now to find an approximate projection $\hat{P}$. If we just use any approximate version of the form of equation (23) we work in a sub-space of the full Fock space with the dimension 2⁴ⁿ . 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 appendix A.3 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 $\mathcal{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 equation (18) with W⁽²⁾ ≡ 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 $\left(\genfrac{}{}{0pt}{}{4n}{2n}\right)$, where we have assumed above that 2n ⩽ M holds. However, if 2n > M the dimension becomes $\left(\genfrac{}{}{0pt}{}{4n}{M}\right)$ (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 section 4.3).

A standard self-consistency condition is then

$$\begin{equation}{\gamma }_{\text{imp}}^{s}({z}_{1},{z}_{2})={\gamma }_{\text{imp}}^{\prime }({z}_{1},{z}_{2}),\end{equation} \tag{ 33 }$$

where ${\gamma }_{\text{imp}}^{s}({z}_{1},{z}_{2})$ is the 1RDM on the impurity of the auxiliary non-interacting system that provides the approximate mean-field projector ${\hat{P}}_{s}$, and γ_imp'(z₁, z₂) is the 1RDM on the impurity of the projected interacting problem with Hamiltonian ${\hat{P}}_{s}\hat{H}{\hat{P}}_{s}$ 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₁, z₂) and thus also the approximate interacting wave function is close to the exact γ_imp(z₁, z₂) and the exact interacting wave function |Ψ⟩.

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 appendix B.2 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 equation (18) without interaction (although there are different choices). The resulting ${\hat{P}}_{s}^{(0)}$ is then used to solve ${\hat{P}}_{s}^{(0)}\hat{H}{\hat{P}}_{s}^{(0)}$, from which we can determine ${\gamma }_{\text{imp}}^{(0)}({z}_{1},{z}_{2})$. In a next step a non-interacting system is constructed such that it reproduces the interacting 1RDM submatrix ${\gamma }_{\text{imp}}^{(0)}({z}_{1},{z}_{2})$ on the impurity A. First we diagonalize on A

$$\begin{equation}{\gamma }_{\text{imp}}^{(0)}({z}_{1},{z}_{2})=\sum\limits _{\mu =1}^{2n}{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{A}^{2}{\varphi }_{\mu }^{A,{\ast}}({z}_{1}){\varphi }_{\mu }^{A}({z}_{2}),\end{equation} \tag{ 34 }$$

where we have denoted the corresponding natural occupation numbers and natural orbitals in accordance to equation (26). Now we only need to reverse the steps that led to equation (26). Firstly we choose 2n arbitrary states ${\varphi }_{\mu }^{B}(z)$ that are orthonormalized on B. Since B has a size of (2N − 2n) we have as many choices. With ${\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{B}^{2}=1-{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{A}^{2}$ we then define for μ ∈ {1, ..., 2n} states

$$\begin{equation}{~{\varphi }}_{\mu }(z)={\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{A}{\varphi }_{\mu }^{A}(z)+{(1-{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{A}^{2})}^{1/2}{\varphi }_{\mu }^{B}(z)\end{equation} \tag{ 35 }$$

(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 ${\varphi }_{\mu }^{B}(z)$ (of course orthogonal to the previous 2n) and define for μ ∈ {2n + 1, ..., M} states

$$\begin{equation}{~{\varphi }}_{\mu }(z)\equiv {\varphi }_{\mu }^{B}(z).\end{equation} \tag{ 36 }$$

We have thus constructed M orthonormal single-particle states ${~{\varphi }}_{\mu }(z)$ with μ ∈ {1, ..., M} on ${\mathcal{H}}_{1}$. Since ${\mathcal{H}}_{1}$ has a dimension of 2N, we are left with (2N − M) further orthonormal states that we again order arbitrarily and denote by ${~{\varphi }}_{\mu }(z)$ for μ ∈ {M + 1, ..., 2N}. As a final step we choose arbitrary energies ${~{{\epsilon}}}_{\mu }\in \mathbb{R}$ such that

$$\begin{equation}{~{{\epsilon}}}_{1}{\leqslant}\cdots {\leqslant}{~{{\epsilon}}}_{M}{< }{~{{\epsilon}}}_{M+1}{\leqslant}\cdots {\leqslant}{~{{\epsilon}}}_{2N}.\end{equation} \tag{ 37 }$$

With these ingredients we find a single-particle Hamiltonian

$$\begin{equation}{~{H}}^{(1)}(z,{z}^{\prime })=\sum\limits _{\mu =1}^{2N}{~{{\epsilon}}}_{\mu }{~{\varphi }}_{\mu }^{{\ast}}({z}^{\prime }){~{\varphi }}_{\mu }(z)\end{equation} \tag{ 38 }$$

and a corresponding Fock-space Hamiltonian with ${\hat{~{\varphi }}}_{\mu }^{{\dagger}}={\sum }_{z=1}^{2N}{~{\varphi }}_{\mu }(z){\hat{c}}_{z}^{{\dagger}}$ that has as its M-particle ground state $\vert ~{{\Phi}}\rangle ={\hat{~{\varphi }}}_{M}^{{\dagger}}\dots {\hat{~{\varphi }}}_{1}^{{\dagger}}\vert 0\rangle$. And by construction ${\gamma }_{s}({z}_{1},{z}_{2})=\langle ~{{\Phi}}\vert {\hat{c}}_{{z}_{1}}^{{\dagger}}{\hat{c}}_{{z}_{2}}~{{\Phi}}\rangle \equiv {\gamma }_{\text{imp}}^{(0)}({z}_{1},{z}_{2})$ if restricted to z₁ and z₂ ∈ A.

Let us note that we have just shown that there are infinitely many ${~{H}}^{(1)}(z,{z}^{\prime })$ that reproduce a given impurity 1RDM. Except of ${\varphi }_{\mu }^{A}(z)$ every other part of our construction is completely arbitrary. Yet different choices generate different projections ${\hat{P}}_{s}^{(1)}$ and corresponding subspaces ${\mathcal{E}}^{(1)}$. And if we now proceed with our iteration, each of this projector will lead to a different ${\hat{P}}_{s}^{(1)}\hat{H}{\hat{P}}_{s}^{(1)}$ and consequently different |Ψ⁽¹⁾⟩ as well as ${\gamma }_{\text{imp}}^{(1)}({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 [8, 33]. 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 equation (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 section 2.3), prohibits usually the use of an auxiliary non-interacting system in 1RDM functional theory [8, 33]. 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 [11]. This will be discussed briefly also in section 6.

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 $\mathcal{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 $\vert ~{K}\rangle$. 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 equation (21), 2n ⩽ N (such that the impurity is smaller or equal to the rest of the system). Only upon increasing the dimension of the CAS 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 section 3.2 is not possible anymore. Until now we have assumed that 2n ⩽ M while for 2n = M all orbitals contribute to the CAS and $\vert ~{0}\rangle \equiv \vert 0\rangle$. 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 equation (23). For a single Slater determinant we know from equation (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⁴ⁿ -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.

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 equation (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 equation (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 equation (4) and the auxiliary Hamiltonian is given by equation (38). But instead of enforcing that $\vert {\Phi}\rangle ={\hat{\phi }}_{M}^{{\dagger}}\dots {\hat{\phi }}_{1}^{{\dagger}}\vert 0\rangle$ and $\vert ~{{\Phi}}\rangle ={\hat{~{\varphi }}}_{M}^{{\dagger}}\dots {\hat{~{\varphi }}}_{1}^{{\dagger}}\vert 0\rangle$ share the same impurity 1RDM, we choose that $\vert ~{{\Phi}}\rangle$ reproduces the impurity 1RDM of $\vert {{\Phi}}^{\prime }\rangle ={\hat{\phi }}_{M+1}^{{\dagger}}\dots {\hat{\phi }}_{2}^{{\dagger}}\vert 0\rangle$. That is, it is not the ground state of the Hamiltonian of equation (4) but an excited state. Furthermore, in the construction that leads to the auxiliary Hamiltonian of equation (38) we choose all ${\varphi }_{\mu }^{B}(z)$ such that

$$\begin{equation}{\phi }_{1}(z)\hspace{2pt}\perp \hspace{2pt}\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{{\varphi }_{1}^{A}(z),\dots ,{\varphi }_{M}^{B}(z)\right\}.\end{equation} \tag{ 39 }$$

If N is large enough, i.e. 2N > 2n + M, this is always possible. The approximate projector ${\hat{P}}_{s}$ and its subspace $\mathcal{E}$ then exclude the actual ground state |Φ⟩ of the M-particle sector of the Hamiltonian of equation (4) and a minimization leads to |Φ'⟩ and the corresponding projection ${\hat{P}}_{s}^{\prime }$. This implies that ${\hat{P}}_{s}^{\prime }\hat{H}{\hat{P}}_{s}^{\prime }$ and ${\hat{P}}_{s}\hat{~{H}}{\hat{P}}_{s}$ share the same impurity 1RDMs and the self-consistency condition of equation (33) is fulfilled. And instead of |Φ⟩ we find |Φ'⟩ at the fixed point (see also appendix B.2 for an explicit example). Realizing that we can easily construct a fixed-point solution that is even further away from |Φ⟩ by choosing the ${\varphi }_{\mu }^{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 section 5. Alternatively, to overcome these ambiguities, the self-consistency condition is adapted or a global iteration is employed instead. We discuss these two options first.

The crucial problem of the self-consistency condition of equation (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

$$\begin{equation}\mathrm{min}\enspace {\Vert}{\gamma }_{\text{emb}}^{s}-{\gamma }_{\text{emb}}^{\prime }{{\Vert}}_{2},\end{equation} \tag{ 40 }$$

where ${\gamma }_{\text{emb}}^{s}({z}_{1},{z}_{2})$ is the 1RDM of the auxiliary non-interacting system that provides the approximate mean-field projector ${\hat{P}}_{s}$, and γ_emb'(z₁, z₂) is the 1RDM of the projected interacting problem with Hamiltonian ${\hat{P}}_{s}\hat{H}{\hat{P}}_{s}$ in the respective M-particle sector (see also appendix B.2 for an explicit example). If the full projector is used then z₁ and z₂ 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/core orbitals) φ_μ (z) for μ ∈ {2n + 1, ..., M} are discarded. In this case z₁ and z₂ correspond to the original lattice sites, only for z₁ and z₂ in A (see for an example the embedded 1RDM in CAS representation in equation (B.30) and in spatial representation in equation (B.31) and then compare with the original equation (B.2), which is identical with the one of the full projection) ⁶ .

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 ${\varphi }_{\mu }^{B}(z)$ are determined by the self-consistency condition of equation (40). The exact solution of the minimum condition of equation (40) is always zero. However, this leads to the impractical case of a highly degenerate non-interacting system ⁷ . Restricting instead to only allow for a single Slater determinant for the case of 2n ⩽ M to construct ${\gamma }_{\text{emb}}^{s}({z}_{1},{z}_{2})$ (in which case the minimum of equation (40) is non-zero in general [2, 21, 38]) 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 (${\gamma }_{\text{emb}}^{s}({z}_{1},{z}_{2})$ = ${\gamma }_{\text{emb}}^{\prime }({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.

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 [37].

Firstly, even though we can find for every A_x potentially many auxiliary non-interacting Hamiltonian ${~{H}}_{x}^{(1)}(z,{z}^{\prime })$ 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 section 2.3, that interacting and non-interacting Hamiltonians cannot share the same 1RDM. Instead, similar to section 4.5, 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 section 4.4) 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.

A.1.1. The basic spaces, Hamiltonians and ground-state representations

Following the construction discussed in section 2.1 we first set-up the single-particle Hamiltonian. In our case of spinless particles on a five-site lattice the single-particle space is ${\mathcal{H}}^{1}\cong {h}_{1}\cong {\mathbb{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

$$\begin{equation}{\boldsymbol{H}}^{(1)}\cong {\boldsymbol{h}}^{(1)}=\left(\begin{matrix}\hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ A.1 }$$

Diagonalizing the matrix H⁽¹⁾(i, j) we find five five-dimensional orthonormal eigenvectors ϕ_μ (i), which allows to express ${\hat{H}}^{(1)}(i,j)={\sum }_{\mu =1}^{5}\hspace{2pt}{{\epsilon}}_{\mu }{\phi }_{\mu }(i){\phi }_{\mu }^{{\ast}}(j)$. Furthermore, they allow us to build the anti-symmetrized M-particle space ${\mathcal{H}}_{M}^{\mathrm{F}}$ of dimension $\left(\genfrac{}{}{0pt}{}{5}{M}\right)$ by constructing all possible M-particle Slater determinants as well as to setup the (non-interacting) 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 $\vert i,j,k\left.\right)=\vert i\rangle \otimes \vert j\rangle \otimes \vert k\rangle$ becomes

$$\begin{align}~{{\Phi}}(i,j,k)& =\left(\right.i,j,k\vert \mu ,\nu ,\xi \rangle \\ & =\frac{1}{\sqrt{3!}}\enspace \mathrm{det}\enspace \left\vert \begin{matrix}\hfill {\phi }_{\mu }(i)\hfill & \hfill {\phi }_{\mu }(j)\hfill & \hfill {\phi }_{\mu }(k)\hfill \\ \hfill {\phi }_{\nu }(i)\hfill & \hfill {\phi }_{\nu }(j)\hfill & \hfill {\phi }_{\nu }(k)\hfill \\ \hfill {\phi }_{\xi }(i)\hfill & \hfill {\phi }_{\xi }(j)\hfill & \hfill {\phi }_{\xi }(k)\hfill \end{matrix}\right\vert .\end{align} \tag{ A.2 }$$

This is only non-zero if μ ≠ ν ≠ ξ, which means we have $\left(\genfrac{}{}{0pt}{}{5}{3}\right)=10$ such wave functions. Alternatively, we can construct the three-particle Slater determinants in the non-symmetrized basis $\vert i,j,k\left.\right)$ of the anti-symmetrized site basis |i', j', k'⟩ as

$$\begin{equation}(i,j,k\vert {i}^{\prime },{j}^{\prime },{k}^{\prime }\rangle =\frac{1}{\sqrt{3!}}\enspace \mathrm{det}\enspace \left\vert \begin{matrix}\hfill {\delta }_{i{i}^{\prime }}\hfill & \hfill {\delta }_{j{i}^{\prime }}\hfill & \hfill {\delta }_{k{i}^{\prime }}\hfill \\ \hfill {\delta }_{i{j}^{\prime }}\hfill & \hfill {\delta }_{j{j}^{\prime }}\hfill & \hfill {\delta }_{k{j}^{\prime }}\hfill \\ \hfill {\delta }_{i{k}^{\prime }}\hfill & \hfill {\delta }_{j{k}^{\prime }}\hfill & \hfill {\delta }_{k{k}^{\prime }}\hfill \end{matrix}\right\vert .\end{equation} \tag{ A.3 }$$

We therefore find the above Slater determinant in the anti-symmetrized basis with i < j < k as

$$\begin{align}{\Phi}(i,j,k)& =\langle i,j,k\vert \mu ,\nu ,\xi \rangle \\ & =\mathrm{det}\enspace \left\vert \begin{matrix}\hfill {\phi }_{\mu }(i)\hfill & \hfill {\phi }_{\mu }(j)\hfill & \hfill {\phi }_{\mu }(k)\hfill \\ \hfill {\phi }_{\nu }(i)\hfill & \hfill {\phi }_{\nu }(j)\hfill & \hfill {\phi }_{\nu }(k)\hfill \\ \hfill {\phi }_{\xi }(i)\hfill & \hfill {\phi }_{\xi }(j)\hfill & \hfill {\phi }_{\xi }(k)\hfill \end{matrix}\right\vert .\end{align} \tag{ A.4 }$$

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⁵. Defining the anti-commuting creation and annihilation operators $\left\{{\hat{c}}_{i},{\hat{c}}_{j}^{{\dagger}}\right\}={\delta }_{ij}$, we can create upon acting on the vacuum state |0⟩ an orthonormal basis of 2⁵ states. For instance, we have

$$\begin{equation}{\hat{c}}_{k}^{{\dagger}}{\hat{c}}_{j}^{{\dagger}}{\hat{c}}_{i}^{{\dagger}}\vert 0\rangle =\vert \varnothing {\rangle }_{0}\oplus \vert \varnothing {\rangle }_{1}\oplus \vert \varnothing {\rangle }_{2}\oplus \vert i,j,k{\rangle }_{3}\oplus \vert \varnothing {\rangle }_{4}\oplus \vert \varnothing {\rangle }_{5}\equiv \vert i,j,k\rangle ,\end{equation} \tag{ A.5 }$$

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 equation (A.3) as well as to the three-particle Fock state of equation (A.5). Dimensionally these two states are different since they belong to different spaces. While $\vert i,j,k\rangle \in {\mathcal{H}}_{3}^{\mathrm{F}}$ is a $\left(\genfrac{}{}{0pt}{}{5}{3}\right)$-dimensional vector with a single non-zero entry, $\vert i,j,k\rangle \in \mathcal{F}$ is a 2⁵-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 three-particle Hamiltonian reads ${\sum }_{\mu =1}^{5}{\sum }_{\nu { >}\mu }^{5}{\sum }_{\xi { >}\nu }^{5}({{\epsilon}}_{\mu }+{{\epsilon}}_{\nu }+{{\epsilon}}_{\xi })\vert \mu ,\nu ,\xi \rangle \langle \mu ,\nu ,\xi \vert$. Expressing the eigenstates in the anti-symmetrized site basis |i, ..., k⟩ we find the Fock-space Hamiltonian in the concise form of

$$\begin{equation}\hat{H}=-\sum\limits _{\langle i,j\rangle }({\hat{c}}_{i}^{{\dagger}}{\hat{c}}_{j}+{\hat{c}}_{j}^{{\dagger}}{\hat{c}}_{i}),\end{equation} \tag{ A.6 }$$

where ⟨i, j⟩ indicates summation only over next neighbors. If we then consider the three-particle subspace, the minimal-energy solution is simply

$$\begin{equation}\vert {\Phi}\rangle =\prod\limits _{\mu =1}^{3}{\hat{\phi }}_{\mu }^{{\dagger}}\vert 0\rangle ={\hat{\phi }}_{1}^{{\dagger}}{\hat{\phi }}_{2}^{{\dagger}}{\hat{\phi }}_{3}^{{\dagger}}\vert 0\rangle ,\end{equation} \tag{ A.7 }$$

where every orbital creation operator is defined by

$$\begin{equation}{\hat{\phi }}_{\mu }^{{\dagger}}=\sum\limits _{k=1}^{5}\underset{={C}_{k\mu }}{\underbrace{{\phi }_{\mu }(k)}}{\hat{c}}_{k}^{{\dagger}}.\end{equation} \tag{ A.8 }$$

More compactly this reads as

$$\begin{equation}\vert {\Phi}\rangle =\prod\limits _{\mu =1}^{3}\sum\limits _{k=1}^{5}{C}_{k\mu }{\hat{c}}_{k}^{{\dagger}}\vert 0\rangle ,\end{equation} \tag{ A.9 }$$

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

$$\begin{equation}\mathbf{C}=\left(\begin{matrix}\hfill \frac{1}{\sqrt{12}}\hfill & \hfill -0.5\hfill & \hfill \frac{1}{\sqrt{3}}\hfill \\ \hfill 0.5\hfill & \hfill -0.5\hfill & \hfill 0\hfill \\ \hfill \frac{1}{\sqrt{3}}\hfill & \hfill 0\hfill & \hfill -\frac{1}{\sqrt{3}}\hfill \\ \hfill 0.5\hfill & \hfill 0.5\hfill & \hfill 0\hfill \\ \hfill \frac{1}{\sqrt{12}}\hfill & \hfill 0.5\hfill & \hfill \frac{1}{\sqrt{3}}\hfill \end{matrix}\right).\end{equation} \tag{ A.10 }$$

The three lowest eigenvalues of (A.1) are ($-\sqrt{3}$, −1, 0), so the ground state energy of the three particle system is $E=-\sqrt{3}-1$. We note that |Φ⟩ as defined in equation (A.9) is defined only in Fock space. Only upon projecting onto the three-particle subspace we find a wave function of the form of equation (A.4). 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 equation (13). For instance we can define

$$\begin{align}\begin{aligned}\vert {F}_{1}\rangle & =\vert 0\rangle \\ \vert {F}_{2}\rangle & ={\hat{c}}_{1}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{3}\rangle & ={\hat{c}}_{2}^{{\dagger}}\vert 0\rangle \\ \dots \\ \vert {F}_{i}\rangle & ={\hat{c}}_{l}^{{\dagger}}{\hat{c}}_{j}^{{\dagger}}{\hat{c}}_{k}^{{\dagger}}\dots \vert 0\rangle \\ \dots \\ \vert {F}_{{2}^{5}}\rangle & ={\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{2}^{{\dagger}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle .\end{aligned}\end{align} \tag{ A.11 }$$

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 equation (19) as

$$\begin{equation}\vert {\Phi}\rangle =\sum\limits _{i=1}^{{2}^{5}}{{\Phi}}_{i}\vert {F}_{i}\rangle =\sum\limits _{i=1}^{10}{{\Phi}}_{i}^{\prime }\vert {F}_{i}^{\prime }\rangle ,\end{equation} \tag{ A.12 }$$

where we have used the prime to denote only the three-particle basis functions $\vert {F}_{i}^{\prime }\rangle$ in Fock space:

$$\begin{align}\begin{aligned}\vert {F}_{1}^{\prime }\rangle & ={\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{2}^{{\dagger}}{\hat{c}}_{3}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{2}^{\prime }\rangle & ={\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{2}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{3}^{\prime }\rangle & ={\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{2}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{4}^{\prime }\rangle & ={\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{5}^{\prime }\rangle & ={\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{6}^{\prime }\rangle & ={\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{7}^{\prime }\rangle & ={\hat{c}}_{2}^{{\dagger}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{8}^{\prime }\rangle & ={\hat{c}}_{2}^{{\dagger}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{9}^{\prime }\rangle & ={\hat{c}}_{2}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ \vert {F}_{10}^{\prime }\rangle & ={\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle .\end{aligned}\end{align} \tag{ A.13 }$$

In this basis we can find a different expression for the Fock-space Hamiltonian of equation (A.6), which is implicitly restricted to the three-particle subspace. Diagonalizing the resulting 10 × 10 matrix with matrix elements $\langle {F}_{i}^{\prime }\vert \hat{H}\vert {F}_{j}^{\prime }\rangle$ leads to the following expansion coefficients ${{\Phi}}_{i}^{\prime }$ that appear in equation (A.12):

$$\begin{align}\begin{aligned}{{\Phi}}_{1}^{\prime }& ={{\Phi}}_{10}^{\prime }=0.105\enspace 66\\ {{\Phi}}_{2}^{\prime }& ={{\Phi}}_{3}^{\prime }={{\Phi}}_{6}^{\prime }={{\Phi}}_{7}^{\prime }={{\Phi}}_{9}^{\prime }=0.288\enspace 67\\ {{\Phi}}_{4}^{\prime }& ={{\Phi}}_{8}^{\prime }=0.394\enspace 34\\ {{\Phi}}_{5}^{\prime }& =0.5.\end{aligned}\end{align} \tag{ A.14 }$$

To see that this agrees with the definition of equation (A.9), we compare with C_kμ . To do so we carry out the sum and the product appearing in equation (A.9) and the wave function is then expressed in the Fock space basis (A.13). We find that

$$\begin{align}{{\Phi}}_{i}^{\prime }=\mathrm{det}\enspace \left\vert \begin{matrix}\hfill {C}_{j,1}\hfill & \hfill {C}_{j,2}\hfill & \hfill {C}_{j,3}\hfill \\ \hfill {C}_{k,1}\hfill & \hfill {C}_{k,2}\hfill & \hfill {C}_{k,3}\hfill \\ \hfill {C}_{l,1}\hfill & \hfill {C}_{l,2}\hfill & \hfill {C}_{l,3}\hfill \end{matrix}\right\vert \end{align} \tag{ A.15 }$$

which is the coefficient associated to a basis function $\vert {F}_{i}^{\prime }\rangle$, with the sites j, k, l occupied. For instance, for $\vert {{\Phi}}_{2}^{\prime }\rangle$ we have j = 1, k = 2, l = 4. Carrying out this procedure for all the terms appearing in equation (A.9) one can verify that it is the same wave function as in equation (A.12). 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. equation (A.15). Moreover, diagonalizing (A.6) in the Fock space basis gives again the same ground-state energy E as the sum of the three lowest orbital energies of (A.1), i.e. $E=-\sqrt{3}-1$.

In the previous part of the appendix we have highlighted the connections between the single-particle, the multi-particle 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 $\mathcal{E}$.

A.2.1. 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 an 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 equation (20) in a way that it will comprise of basis functions $\vert {F}_{i}^{A}\rangle$ that belong only to the impurity and $\vert {F}_{i}^{B}\rangle$ that will belong only to the environment. The number of linearly independent $\vert {F}_{i}^{A}\rangle$ that we get is 2², while for $\vert {F}_{i}^{B}\rangle$ is 2³. For this we define two new vacua |0⟩_A and |0⟩_B and define ${\hat{c}}_{1}^{{\dagger}}$ and ${\hat{c}}_{2}^{{\dagger}}$ only on |0⟩_A (which leads to a Fock space ${\mathcal{F}}_{A}$) and accordingly for sites 3, 4 and 5 that constitute B and ${\mathcal{F}}_{B}$. By dimensional correspondence we find $\mathcal{F}\cong {\mathcal{F}}_{A}\otimes {\mathcal{F}}_{B}$. The Fock states of the respective Fock spaces are

$$\begin{align}\begin{aligned}\vert {F}_{1}^{A}\rangle & =\vert 0{\rangle }_{A}\\ \vert {F}_{2}^{A}\rangle & ={\hat{c}}_{2}^{{\dagger}}\vert 0{\rangle }_{A}\\ \vert {F}_{3}^{A}\rangle & ={\hat{c}}_{1}^{{\dagger}}\vert 0{\rangle }_{A}\\ \vert {F}_{4}^{A}\rangle & ={\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{2}^{{\dagger}}\vert 0{\rangle }_{A}\end{aligned}\end{align} \tag{ A.16 }$$

and

$$\begin{align}\begin{aligned}\vert {F}_{1}^{B}\rangle & =\vert 0{\rangle }_{B}\\ \vert {F}_{2}^{B}\rangle & ={\hat{c}}_{3}^{{\dagger}}\vert 0{\rangle }_{B}\\ \vert {F}_{3}^{B}\rangle & ={\hat{c}}_{4}^{{\dagger}}\vert 0{\rangle }_{B}\\ \vert {F}_{4}^{B}\rangle & ={\hat{c}}_{5}^{{\dagger}}\vert 0{\rangle }_{B}\\ \vert {F}_{5}^{B}\rangle & ={\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}\vert 0{\rangle }_{B}\\ \vert {F}_{6}^{B}\rangle & ={\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0{\rangle }_{B}\\ \vert {F}_{7}^{B}\rangle & ={\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0{\rangle }_{B}\\ \vert {F}_{8}^{B}\rangle & ={\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0{\rangle }_{B}.\end{aligned}\end{align} \tag{ A.17 }$$

Similar to the local creation and annihilation operators (see section 2.1), 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 $\mathcal{F}$

$$\begin{equation}\vert {\Phi}\rangle =\sum\limits _{i=1}^{{2}^{2}}\sum\limits _{j=1}^{{2}^{3}}{{\Phi}}_{i,j}\vert {F}_{i}^{A}\rangle \otimes \vert {F}_{j}^{B}\rangle \end{equation} \tag{ A.18 }$$

the coefficients need to take care of the proper symmetry. The induced basis of $\mathcal{F}$ then corresponds to the previously introduced basis of equation (A.11). This means that the coefficients Φ_i,j are identical with the coefficients Φ_i that appear in equation (A.12). This means that there will be only 10 non-zero entries of Φ_i,j . Recasting now equation (A.18) in terms of only its non-zero entries we find that

$$\begin{align}\vert {\Phi}\rangle & ={{\Phi}}_{10}^{\prime }{\overbrace{{{\Phi}}_{1,8}}}\vert {F}_{10}^{\prime }\rangle {\overbrace{\vert {F}_{1}^{A}\rangle \otimes \vert {F}_{8}^{B}\rangle }}+{{\Phi}}_{7}^{\prime }{\overbrace{{{\Phi}}_{2,5}}}\vert {F}_{7}^{\prime }\rangle {\overbrace{\vert {F}_{2}^{A}\rangle \otimes \vert {F}_{5}^{B}\rangle }}\\ & \quad +{{\Phi}}_{8}^{\prime }{\overbrace{{{\Phi}}_{2,6}}}\vert {F}_{8}^{\prime }\rangle {\overbrace{\vert {F}_{2}^{A}\rangle \otimes \vert {F}_{6}^{B}\rangle }}+{{\Phi}}_{9}^{\prime }{\overbrace{{{\Phi}}_{2,7}}}\vert {F}_{9}^{\prime }\rangle {\overbrace{\vert {F}_{2}^{A}\rangle \otimes \vert {F}_{7}^{B}\rangle }}\\ & \quad +{{\Phi}}_{4}^{\prime }{\overbrace{{{\Phi}}_{3,5}}}\vert {F}_{4}^{\prime }\rangle {\overbrace{\vert {F}_{3}^{A}\rangle \otimes \vert {F}_{5}^{B}\rangle }}+{{\Phi}}_{5}^{\prime }{\overbrace{{{\Phi}}_{3,6}}}\vert {F}_{5}^{\prime }\rangle {\overbrace{\vert {F}_{3}^{A}\rangle \otimes \vert {F}_{6}^{B}\rangle }}\\ & \quad +{{\Phi}}_{6}^{\prime }{\overbrace{{{\Phi}}_{3,7}}}\vert {F}_{6}^{\prime }\rangle {\overbrace{\vert {A}_{3}\rangle \otimes \vert {F}_{7}^{B}\rangle }}+{{\Phi}}_{1}^{\prime }{\overbrace{{{\Phi}}_{4,2}}}\vert {F}_{1}^{\prime }\rangle {\overbrace{\vert {F}_{4}^{A}\rangle \otimes \vert {F}_{2}^{B}\rangle }}\\ & \quad +{{\Phi}}_{2}^{\prime }{\overbrace{{{\Phi}}_{4,3}}}\vert {F}_{2}^{\prime }\rangle {\overbrace{\vert {A}_{4}\rangle \otimes \vert {F}_{3}^{B}\rangle }}+{{\Phi}}_{3}^{\prime }{\overbrace{{{\Phi}}_{4,4}}}\vert {F}_{3}^{\prime }\rangle {\overbrace{\vert {F}_{4}^{A}\rangle \otimes \vert {F}_{4}^{B}\rangle }},\end{align} \tag{ A.19 }$$

where the matrix elements ${{\Phi}}_{i}^{\prime }$ are given by equation (A.15). Having written the wave function of our example in the form of equation (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)

$$\begin{equation}\mathbf{\Phi }=\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {{\Phi}}_{10}^{\prime }\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {{\Phi}}_{7}^{\prime }\hfill & \hfill {{\Phi}}_{8}^{\prime }\hfill & \hfill {{\Phi}}_{9}^{\prime }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {{\Phi}}_{4}^{\prime }\hfill & \hfill {{\Phi}}_{5}^{\prime }\hfill & \hfill {{\Phi}}_{6}^{\prime }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {{\Phi}}_{1}^{\prime }\hfill & \hfill {{\Phi}}_{2}^{\prime }\hfill & \hfill {{\Phi}}_{3}^{\prime }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ A.20 }$$

The connecting matrix for our example reads

$$\begin{equation}\mathbf{\Phi }=\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.105\enspace 66\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.288\enspace 68\hfill & \hfill 0.394\enspace 34\hfill & \hfill 0.288\enspace 68\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.394\enspace 34\hfill & \hfill 0.5\hfill & \hfill 0.288\enspace 68\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0.105\enspace 66\hfill & \hfill 0.288\enspace 68\hfill & \hfill 0.288\enspace 68\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ A.21 }$$

Rotating with U and V^†, which are defined in equation (22) and the discussion after it, we obtain the following states on the impurity

$$\begin{equation}\vert {A}_{1}\rangle =-0.629\enspace 78\vert {F}_{2}^{A}\rangle -0.776\enspace 78\vert {F}_{3}^{A}\rangle \end{equation} \tag{ A.22 }$$

$$\begin{equation}\vert {A}_{2}\rangle =-\vert {F}_{4}^{A}\rangle \end{equation} \tag{ A.23 }$$

$$\begin{equation}\vert {A}_{3}\rangle =-\vert {F}_{1}^{A}\rangle \end{equation} \tag{ A.24 }$$

$$\begin{equation}\vert {A}_{4}\rangle =0.776\enspace 78\vert {F}_{2}^{A}\rangle -0.629\enspace 78\vert {F}_{3}^{A}\rangle \end{equation} \tag{ A.25 }$$

and on the bath:

$$\begin{align}\vert {B}_{1}\rangle =& -0.542\enspace 83\vert {F}_{5}^{B}\rangle -0.708\enspace 12\vert {F}_{6}^{B}\rangle \\ & -0.451\enspace 557\vert {F}_{7}^{B}\rangle \end{align} \tag{ A.26 }$$

$$\begin{align}\vert {B}_{2}\rangle =& -0.250\enspace 56\vert {F}_{2}^{B}\rangle -0.6846\vert {F}_{3}^{B}\rangle \\ & -0.6846\vert {F}_{4}^{B}\rangle \end{align} \tag{ A.27 }$$

$$\begin{equation}\vert {B}_{3}\rangle =-\vert {F}_{8}^{B}\rangle \end{equation} \tag{ A.28 }$$

$$\begin{align}\vert {B}_{4}\rangle =& -0.486\enspace 54\vert {F}_{5}^{B}\rangle -0.173\enspace 10\vert {F}_{6}^{B}\rangle \\ & +0.856\enspace 34\vert {F}_{7}^{B}\rangle \end{align} \tag{ A.29 }$$

$$\begin{align}\vert {B}_{5}\rangle =& -0.436\enspace 46\vert {F}_{2}^{B}\rangle -0.325\enspace 17\vert {F}_{3}^{B}\rangle \\ & +0.484\enspace 93\vert {F}_{4}^{B}\rangle +0.468\enspace 61\vert {F}_{5}^{B}\rangle -0.468\enspace 61\vert {F}_{6}^{B}\rangle \\ & +0.171\enspace 523\vert {F}_{7}^{B}\rangle \end{align} \tag{ A.30 }$$

$$\begin{align}\vert {B}_{6}\rangle =& -0.646\enspace 38\vert {F}_{2}^{B}\rangle -0.088\enspace 58\vert {F}_{3}^{B}\rangle \\ & +0.325\enspace 17\vert {B}_{4}\rangle -0.468\enspace 61\vert {F}_{5}^{B}\rangle +0.468\enspace 61\vert {F}_{6}^{B}\rangle \\ & -0.171\enspace 523\vert {F}_{7}^{B}\rangle \end{align} \tag{ A.31 }$$

$$\begin{align}\vert {B}_{7}\rangle =& -0.573\enspace 51\vert {F}_{2}^{B}\rangle +0.646\enspace 38\vert {F}_{3}^{B}\rangle \\ & -0.436\enspace 46\vert {F}_{4}^{B}\rangle +0.171\enspace 52\vert {F}_{5}^{B}\rangle -0.171\enspace 52\vert {F}_{6}^{B}\rangle \\ & +0.062\enspace 78\vert {F}_{7}^{B}\rangle \end{align} \tag{ A.32 }$$

$$\begin{equation}\vert {B}_{8}\rangle =\vert {F}_{1}^{B}\rangle \end{equation} \tag{ A.33 }$$

and a new connecting matrix which is diagonal

$$\begin{equation}{\Lambda}=\left(\begin{matrix}0.899\enspace 19& 0.& 0.& 0.& 0.& 0.& 0.& 0.\\ 0.& 0.4217& 0.& 0.& 0.& 0.& 0.& 0.\\ 0.& 0.& 0.105\enspace 66& 0.& 0.& 0.& 0.& 0.\\ 0.& 0.& 0.& 0.049\enspace 55& 0.& 0.& 0.& 0.\end{matrix}\right).\end{equation} \tag{ A.34 }$$

The wave function after the SVD reads

$$\begin{align}\vert {\Phi}\rangle =& \sum\limits _{i=1}^{4}{\lambda }_{i}\vert {A}_{i}\rangle \otimes \vert {B}_{i}\rangle \end{align} \tag{ A.35 }$$

$$\begin{align}& =0.899\enspace 19\cdot \vert {A}_{1}\rangle {\overbrace{(-0.629\enspace 78\vert {F}_{2}^{A}\rangle -0.776\enspace 78\vert {F}_{3}^{A}\rangle )}}\\ & \quad \otimes \vert {B}_{1}\rangle {\overbrace{(-0.542\enspace 83\vert {F}_{5}^{B}\rangle -0.708\enspace 12\vert {F}_{6}^{B}\rangle -0.451\enspace 56\vert {F}_{7}^{B}\rangle )}}\\ & \quad +0.4217\cdot \vert {A}_{2}\rangle {\overbrace{(-\vert {F}_{4}^{A}\rangle )}}\\ & \quad \otimes \vert {B}_{2}\rangle {\overbrace{(-0.250\enspace 56\vert {F}_{2}^{B}\rangle -0.684\enspace 55\vert {F}_{3}^{B}\rangle -0.684\enspace 55\vert {F}_{4}^{B}\rangle )}}\\ & \quad +0.105\enspace 66\cdot \vert {A}_{3}\rangle {\overbrace{(-\vert {F}_{1}^{A}\rangle )}}\otimes \vert {B}_{3}\rangle {\overbrace{(-\vert {F}_{8}^{B}\rangle )}}\\ & \quad +0.049\enspace 55\cdot \vert {A}_{4}\rangle {\overbrace{(0.776\enspace 78\vert {F}_{2}^{A}\rangle -0.629\enspace 778\vert {F}_{3}^{A}\rangle )}}\\ & \quad \otimes \vert {B}_{4}\rangle {\overbrace{(-0.486\enspace 54\vert {F}_{5}^{B}\rangle -0.173\enspace 10\vert {F}_{6}^{B}\rangle +0.856\enspace 34\vert {F}_{7}^{B}\rangle )}}.\end{align} \tag{ A.36 }$$

If we compare to equation (A.19), we see that this is of course the same wave function.

A.2.2. Projection via singular-value decomposition in the impurity submatrix

Here we perform the projection via an SVD for the non interacting system on the orbital coefficient matrix of the creation operators. Due to our choice of A we have

$$\begin{equation}{\mathbf{C}}^{A}=\left(\begin{matrix}\hfill \frac{1}{\sqrt{12}}\hfill & \hfill -0.5\hfill & \hfill \frac{1}{\sqrt{3}}\hfill \\ \hfill 0.5\hfill & \hfill -0.5\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ A.37 }$$

which corresponds to the first two lines of the matrix defined in equation (A.10). Performing an SVD on C^A we get the factorization for its matrix elements

$$\begin{equation}{C}_{k\mu }^{A}={C}_{k,\in \mu }(k\in A)=\sum\limits _{k=1}^{2}\sum\limits _{\nu =1}^{3}{U}_{\text{k}\;\text{i}}{{\Lambda}}_{i\nu }{V}_{\nu \mu }^{{\dagger}},\end{equation} \tag{ A.38 }$$

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:

$$\begin{equation}\mathbf{U}=\left(\begin{matrix}\hfill -0.776\enspace 78\hfill & \hfill -0.629\enspace 78\hfill \\ \hfill -0.629\enspace 78\hfill & \hfill 0.776\enspace 78\hfill \end{matrix}\right)\end{equation} \tag{ A.39 }$$

$$\begin{equation}\mathbf{\Lambda }=\left(\begin{matrix}\hfill 0.993\enspace 17\hfill & \hfill 0.\hfill & \hfill 0.\hfill \\ \hfill 0.\hfill & \hfill 0.424\enspace 60\hfill & \hfill 0.\hfill \end{matrix}\right)\end{equation} \tag{ A.40 }$$

$$\begin{equation}\mathbf{V}=\left(\begin{matrix}\hfill -0.542\enspace 83\hfill & \hfill 0.486\enspace 54\hfill & \hfill -0.684\enspace 55\hfill \\ \hfill 0.708\enspace 12\hfill & \hfill -0.173\enspace 10\hfill & \hfill -0.684\enspace 55\hfill \\ \hfill -0.451\enspace 56\hfill & \hfill -0.856\enspace 34\hfill & \hfill -0.250\enspace 56\hfill \end{matrix}\right).\end{equation} \tag{ A.41 }$$

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

$$\begin{align}~{\mathbf{C}}& =\mathbf{C}\cdot \mathbf{V}\\ & =\left(\begin{matrix}\hfill \frac{1}{\sqrt{12}}\hfill & \hfill -0.5\hfill & \hfill \frac{1}{\sqrt{3}}\hfill \\ \hfill 0.5\hfill & \hfill -0.5\hfill & \hfill 0\hfill \\ \hfill \frac{1}{\sqrt{3}}\hfill & \hfill 0\hfill & \hfill -\frac{1}{\sqrt{3}}\hfill \\ \hfill 0.5\hfill & \hfill 0.5\hfill & \hfill 0\hfill \\ \hfill \frac{1}{\sqrt{12}}\hfill & \hfill 0.5\hfill & \hfill \frac{1}{\sqrt{3}}.\hfill \end{matrix}\right)\cdot \left(\begin{matrix}\hfill -0.542\enspace 83\hfill & \hfill 0.486\enspace 54\hfill & \hfill -0.684\enspace 55\hfill \\ \hfill 0.708\enspace 12\hfill & \hfill -0.173\enspace 10\hfill & \hfill -0.684\enspace 55\hfill \\ \hfill -0.451\enspace 56\hfill & \hfill -0.856\enspace 34\hfill & \hfill -0.250\enspace 56\hfill \end{matrix}\right)\\ & =\left(\begin{matrix}\hfill -0.771\enspace 47\hfill & \hfill -0.267\enspace 41\hfill & \hfill 0\hfill \\ \hfill -0.625\enspace 48\hfill & \hfill 0.329\enspace 82\hfill & \hfill 0\hfill \\ \hfill -0.052\enspace 70\hfill & \hfill 0.775\enspace 31\hfill & \hfill -0.250\enspace 56\hfill \\ \hfill 0.082\enspace 64\hfill & \hfill 0.156\enspace 72\hfill & \hfill -0.684\enspace 55\hfill \\ \hfill -0.063\enspace 35\hfill & \hfill -0.440\enspace 50\hfill & \hfill -0.684\enspace 55\hfill \end{matrix}\right).\end{align} \tag{ A.42 }$$

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 equation (A.9) is

$$\begin{align}\vert {\Phi}\rangle & =\prod\limits _{\mu =1}^{3}\sum\limits _{k=1}^{5}{~{C}}_{k\mu }{\hat{c}}_{k}^{{\dagger}}\vert 0\rangle \\ & =\left(\prod\limits _{\mu =1}^{2}\sum\limits _{k=1}^{5}{~{C}}_{k\mu }{\hat{c}}_{k}^{{\dagger}}\right)\underset{={\hat{~{\varphi }}}_{3}^{{\dagger}}}{\underbrace{\left(\sum\limits _{k=3}^{5}{~{C}}_{k3}{\hat{c}}_{k}^{{\dagger}}\right)}}\vert 0\rangle \\ & ={\hat{~{\varphi }}}_{1}^{{\dagger}}{\hat{~{\varphi }}}_{2}^{{\dagger}}{\hat{~{\varphi }}}_{3}^{{\dagger}}\vert 0\rangle .\end{align} \tag{ A.43 }$$

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 equation (27) by ${\hat{~{\varphi }}}_{3}^{{\dagger}}={\hat{\varphi }}_{3}^{B,{\dagger}}$.

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

$$\begin{equation}{\hat{\varphi }}_{\mu }^{B,{\dagger}}=\sum\limits _{k=3}^{5}{\hat{c}}_{k}^{{\dagger}}\frac{{~{C}}_{k\mu }}{\sqrt{{\sum }_{l=3}^{5}\vert {~{C}}_{l\mu }{\vert }^{2}}}=\sum\limits _{k=3}^{5}{\hat{c}}_{k}^{{\dagger}}\underset{={\varphi }_{\mu }^{B}(k)}{\underbrace{\frac{{~{C}}_{k\mu }}{{\Vert}{~{\varphi }}_{\mu }{{\Vert}}_{B}}}}.\end{equation} \tag{ A.44 }$$

The normalization factors that appear in the denominator of equation (A.44) read

$$\begin{align}{\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}& =\sqrt{\sum\limits _{l=3}^{5}\vert {~{C}}_{l1}{\vert }^{2}}=0.116\enspace 70\\ {\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}& =\sqrt{\sum\limits _{l=3}^{5}\vert {~{C}}_{l2}{\vert }^{2}}=0.905\enspace 38\end{align} \tag{ A.45 }$$

such that

$$\begin{equation}{\hat{\varphi }}_{1}^{B,{\dagger}}=-0.451\enspace 56{\hat{c}}_{3}^{{\dagger}}+0.708\enspace 12{\hat{c}}_{4}^{{\dagger}}-0.542\enspace 83{\hat{c}}_{5}^{{\dagger}}\end{equation} \tag{ A.46 }$$

$$\begin{equation}{\hat{\varphi }}_{2}^{B,{\dagger}}=0.856\enspace 34{\hat{c}}_{3}^{{\dagger}}+0.173\enspace 10{\hat{c}}_{4}^{{\dagger}}-0.486\enspace 54{\hat{c}}_{5}^{{\dagger}}.\end{equation} \tag{ A.47 }$$

In a similar manner (see equation (27)) we can define the renormalization factors ${\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}=0.993\enspace 17$, ${\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}=0.424\enspace 60$ and the 2 impurity orbitals as

$$\begin{equation}{\hat{\varphi }}_{1}^{A,{\dagger}}=-0.776\enspace 78{\hat{c}}_{1}^{{\dagger}}-0.629\enspace 78{\hat{c}}_{2}^{{\dagger}}\end{equation} \tag{ A.48 }$$

$$\begin{align}\\ {\hat{\varphi }}_{2}^{A,{\dagger}}& =-0.629\enspace 78{\hat{c}}_{1}^{{\dagger}}+0.776\enspace 77{\hat{c}}_{2}^{{\dagger}}.\end{align} \tag{ A.49 }$$

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

$$\begin{align}\vert {\Phi}\rangle & ={\hat{\stackrel{~}{\varphi }}}_{1}^{{\dagger}}{\overbrace{\left({\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\hat{\varphi }}_{1}^{A,{\dagger}}+{\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}{\hat{\varphi }}_{1}^{B,{\dagger}}\right)}}\\ & \quad \cdot {\hat{\stackrel{~}{\varphi }}}_{2}^{{\dagger}}{\overbrace{\left({\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}{\hat{\varphi }}_{2}^{A,{\dagger}}+{\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}{\hat{\varphi }}_{2}^{B,{\dagger}}\right)}}{\hat{\varphi }}^{B,{\dagger}}\vert 0\rangle \\ & ={\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \\ & \quad +{\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \\ & \quad +{\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \\ & \quad +{\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle .\end{align} \tag{ A.50 }$$

Since we have now four terms in accordance to equation (A.35), we can individually compare. Firstly we find that

$$\begin{equation}{\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}=0.421\enspace 70={\lambda }_{2}\end{equation} \tag{ A.51 }$$

and the corresponding vectors can be associated as (note the anti-symmetrization)

$$\begin{align}{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}\vert 0\rangle & =\underset{=-1}{\underbrace{({\varphi }_{1}^{A}(1){\varphi }_{2}^{A}(2)-{\varphi }_{2}^{A}(1){\varphi }_{1}^{A}(2))}}{\hat{c}}_{1}^{{\dagger}}{\hat{c}}_{2}^{{\dagger}}\vert 0\rangle \\ & \equiv \vert {A}_{2}\rangle ,\end{align} \tag{ A.52 }$$

and

$$\begin{align}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle & =\left(-0.250\enspace 56{\hat{c}}_{3}^{{\dagger}}-0.684\enspace 55{\hat{c}}_{4}^{{\dagger}}-0.684\enspace 55{\hat{c}}_{5}^{{\dagger}}\right)\vert 0\rangle \\ & \equiv \vert {B}_{2}\rangle .\end{align} \tag{ A.53 }$$

We note here again that the left-hand sides of the above equivalence relations are vectors that are defined on a Fock space of the full lattice, while the right-hand sides are defined on Fock spaces of sub-lattices. Consequently they are not the same vectors since they are defined on dimensionally different spaces, yet they describe the same physical states. This allows us to associate

$$\begin{equation}{\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \equiv {\lambda }_{2}\vert {A}_{2}\rangle \otimes \vert {B}_{2}\rangle \end{equation} \tag{ A.54 }$$

while dimensionally (and also physically) ${\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}\vert 0\rangle \otimes {\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle$ would not make sense. This highlights again the subtleties that arise when mixing different representations of the same physical state. We can then proceed by associating the other states in a similar manner. Since ${\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}=0.899\enspace 19={\lambda }_{1}$,

$$\begin{align}{\hat{\varphi }}_{1}^{A,{\dagger}}\vert 0\rangle & =\left(-0.776\enspace 77{\hat{c}}_{1}^{{\dagger}}-0.629\enspace 78{\hat{c}}_{2}^{{\dagger}}\right)\vert 0\rangle \\ & \equiv \vert {A}_{1}\rangle \end{align} \tag{ A.55 }$$

and

$$\begin{align*}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle & =\underset{=-0.542\enspace 83}{\underbrace{\left({\varphi }_{2}^{B}(3){\varphi }_{3}^{B}(4)-{\varphi }_{2}^{B}(4){\varphi }_{3}^{B}(3)\right)}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}\vert 0\rangle \\ & \quad +\underset{=-0.708\enspace 12}{\underbrace{\left({\varphi }_{2}^{B}(3){\varphi }_{3}^{B}(5)-{\varphi }_{2}^{B}(5){\varphi }_{3}^{B}(3)\right)}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ & \quad +\underset{=-0.451\enspace 56}{\underbrace{\left({\varphi }_{2}^{B}(4){\varphi }_{3}^{B}(5)-{\varphi }_{2}^{B}(5){\varphi }_{3}^{B}(4)\right)}}{\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ & \equiv \vert {B}_{1}\rangle .\end{align*}$$

we have ${\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \equiv {\lambda }_{1}\vert {A}_{1}\rangle \otimes \vert {B}_{1}\rangle$. Further, since ${\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}=0.049\enspace 55={\lambda }_{4}$,

$$\begin{equation}{\hat{\varphi }}_{2}^{A,{\dagger}}\vert 0\rangle =\left(0.776\enspace 78{\hat{c}}_{2}^{{\dagger}}-0.629\enspace 78{\hat{c}}_{1}^{{\dagger}}\right)\vert 0\rangle \equiv \vert {A}_{4}\rangle \end{equation} \tag{ A.56 }$$

and

$$\begin{align*}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle & =\underset{=0.486\enspace 54}{\underbrace{\left({\varphi }_{1}^{B}(3){\varphi }_{3}^{B}(4)-{\varphi }_{1}^{B}(4){\varphi }_{3}^{B}(3)\right)}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}\vert 0\rangle \\ & \quad +\underset{=0.173\enspace 09}{\underbrace{\left({\varphi }_{1}^{B}(3){\varphi }_{3}^{B}(5)-{\varphi }_{1}^{B}(5){\varphi }_{3}^{B}(3)\right)}}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ & \quad +\underset{=-0.856\enspace 34}{\underbrace{\left({\varphi }_{1}^{B}(4){\varphi }_{3}^{B}(5)-{\varphi }_{1}^{B}(5){\varphi }_{3}^{B}(4)\right)}}{\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ & \equiv -\vert {B}_{4}\rangle ,\end{align*}$$

we have ${\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \equiv {\lambda }_{4}\vert {A}_{4}\rangle \otimes \vert {B}_{4}\rangle$. Finally, since ${\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}=0.105\enspace 66={\lambda }_{3}$,

$$\begin{equation}\vert 0\rangle \equiv -\vert {A}_{3}\rangle \end{equation} \tag{ A.57 }$$

and

$$\begin{align*}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle & =\underset{=1}{\underbrace{\mathrm{det}\enspace \left\vert \begin{matrix}\hfill {\varphi }_{1}^{B}(3)\hfill & \hfill {\varphi }_{1}^{B}(4)\hfill & \hfill {\varphi }_{1}^{B}(5)\hfill \\ \hfill {\varphi }_{2}^{B}(3)\hfill & \hfill {\varphi }_{2}^{B}(4)\hfill & \hfill {\varphi }_{2}^{B}(5)\hfill \\ \hfill {\varphi }_{3}^{B}(3)\hfill & \hfill {\varphi }_{3}^{B}(4)\hfill & \hfill {\varphi }_{3}^{B}(5)\hfill \end{matrix}\right\vert }}{\hat{c}}_{3}^{{\dagger}}{\hat{c}}_{4}^{{\dagger}}{\hat{c}}_{5}^{{\dagger}}\vert 0\rangle \\ & \equiv \vert {B}_{3}\rangle ,\end{align*}$$

we have ${\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}{\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \equiv {\lambda }_{3}\vert {A}_{3}\rangle \otimes \vert {B}_{3}\rangle$. Thus we have explicitly verified that performing the SVD in the orbital coefficient matrix with a subsequent division in impurity and environment orbitals is equivalent to performing the SVD in the connecting matrix in the Fock-space representation.

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 appendix A.4), the latter approach is the one that is employed in practice and is only exact for non-interacting systems.

A.3.1. Exact embedded Hamiltonian via the Fock-space projection

The Fock-space projection according to equation (23) reads in our case

$$\begin{equation}\hat{P}=\sum\limits _{\alpha =1}^{4}\sum\limits _{\beta =1}^{4}\vert {A}_{\alpha }\rangle \otimes \vert {B}_{\beta }\rangle \langle {A}_{\alpha }\vert \otimes \langle {B}_{\beta }\vert ,\end{equation} \tag{ A.58 }$$

with |A_α ⟩ given by equations (A.22)–(A.25) and |B_β ⟩ given by equations (A.26)–(A.29). 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 associate

$$\begin{align}\hat{P}& \equiv {\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{B}{\hat{\varphi }}_{1}^{A}\\ & \quad +{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{1}^{A}\\ & \quad +{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{B}{\hat{\varphi }}_{1}^{B}{\hat{\varphi }}_{1}^{A}\\ & \quad +{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{1}^{B}{\hat{\varphi }}_{1}^{A}\\ & \quad +{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{B}{\hat{\varphi }}_{2}^{A}{\hat{\varphi }}_{1}^{A}\\ & \quad +{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{A}{\hat{\varphi }}_{1}^{A}\\ & \quad +{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{B}{\hat{\varphi }}_{1}^{B}{\hat{\varphi }}_{2}^{A}{\hat{\varphi }}_{1}^{A}\\ & \quad +{\hat{\varphi }}_{1}^{A,{\dagger}}{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{1}^{B}{\hat{\varphi }}_{2}^{A}{\hat{\varphi }}_{1}^{A}\\ & \quad +{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}\\ & \quad +{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{B}\\ & \quad +{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{B}{\hat{\varphi }}_{1}^{B}\\ & \quad +{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{1}^{B}\\ & \quad +{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{A}\\ & \quad +{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{B}{\hat{\varphi }}_{2}^{A}\\ & \quad +{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{2}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{2}^{B}{\hat{\varphi }}_{1}^{B}{\hat{\varphi }}_{2}^{A}\\ & \quad +{\hat{\varphi }}_{2}^{A,{\dagger}}{\hat{\varphi }}_{1}^{B,{\dagger}}{\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle \langle 0\vert {\hat{\varphi }}_{3}^{B}{\hat{\varphi }}_{1}^{B}{\hat{\varphi }}_{2}^{A}.\end{align} \tag{ A.59 }$$

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 $\langle {\varphi }_{1}^{A}{\varphi }_{2}^{B}{\varphi }_{3}^{B}\vert \hat{H}{\varphi }_{1}^{A}{\varphi }_{3}^{B}\rangle$ 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 becomes

$$\begin{align}{\hat{P}}^{\prime }& =\vert {A}_{1}\rangle \otimes \vert {B}_{1}\rangle \langle {B}_{1}\vert \otimes \langle {A}_{1}\vert \\ & \quad +\vert {A}_{2}\rangle \otimes \vert {B}_{2}\rangle \langle {B}_{2}\vert \otimes \langle {A}_{2}\vert \\ & \quad +\vert {A}_{3}\rangle \otimes \vert {B}_{3}\rangle \langle {B}_{3}\vert \otimes \langle {A}_{3}\vert \\ & \quad +\vert {A}_{4}\rangle \otimes \vert {B}_{4}\rangle \langle {B}_{4}\vert \otimes \langle {A}_{4}\vert \\ & \quad +\vert {A}_{1}\rangle \otimes \vert {B}_{4}\rangle \langle {B}_{4}\vert \otimes \langle {A}_{1}\vert \\ & \quad +\vert {A}_{4}\rangle \otimes \vert {B}_{1}\rangle \langle {B}_{1}\vert \otimes \langle {A}_{4}\vert .\end{align} \tag{ A.60 }$$

This yields a 6 × 6 Hamiltonian matrix ${\hat{P}}^{\prime }\hat{H}{\hat{P}}^{\prime }$ (see the Jupyter notebook for the explicit matrix). Diagonalizing provides its lowest eigenvalue as $E=-\sqrt{3}-1$ with the eigenstate Φ = −0.899 19|A₁⟩ ⊗ |B₁⟩ − 0.421 70|A₂⟩ ⊗ |B₂⟩ − 0.105 66|A₃⟩ ⊗ |B₃⟩ − 0.049 55|A₄⟩ ⊗ |B₄⟩. This agrees with equation (A.35). 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.

A.3.2. 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 section 3.2. 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. equation (A.1) 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

$$\begin{equation*}\gamma (i,j)=\langle {\Phi}\vert {\hat{c}}_{j}^{{\dagger}}{\hat{c}}_{i}\vert {\Phi}\rangle =\sum\limits _{\mu =1}^{3}{C}_{j\mu }^{\mathrm{T}}{C}_{i\mu }.\end{equation*}$$

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}

$$\begin{equation}{\gamma }_{\text{env}}(i,j)=\sum\limits _{\mu =1}^{3}{C}_{j\mu }^{\mathrm{T}}{C}_{i\mu }\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\enspace i,j\in B.\end{equation} \tag{ A.61 }$$

Diagonalizing this 3 × 3 submatrix gives

$$\begin{align*}{n}_{1}& =0.013\enspace 62\equiv {\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}^{2}\\ {\varphi }_{1}^{B}& =-0.451\enspace 56\vert 3\rangle +0.708\enspace 12\vert 4\rangle -0.542\enspace 83\vert 5\rangle \\ {n}_{2}& =0.819\enspace 71\equiv {\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}^{2}\\ {\varphi }_{2}^{B}& =0.856\enspace 34\vert 3\rangle +0.173\enspace 10\vert 4\rangle -0.486\enspace 54\vert 5\rangle \\ {n}_{3}& =1\equiv {\Vert}{~{\varphi }}_{3}{{\Vert}}_{B}^{2}\\ {\varphi }_{3}^{B}& =0.250\enspace 56\vert 3\rangle +0.684\enspace 55\vert 4\rangle +0.684\enspace 55\vert 5\rangle ,\end{align*}$$

where we have introduced the notation ${\hat{c}}_{i}^{{\dagger}}\vert 0\rangle =\vert i\rangle$. The orbital ${\varphi }_{3}^{B}$ with occupation number 1 is called in the DMET literature unentangled occupied environmental orbital. It agrees with ${~{\varphi }}_{3}\equiv {\varphi }_{3}^{B}$ from appendix A.2.2. The two orbitals that have eigenvalues (occupations) between 0 and 1 are called the bath orbitals and agree with the corresponding ones from appendix A.2.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 ${\varphi }_{1}^{A}$ and ${\varphi }_{2}^{A}$ from appendix A.2.2, we can equivalently use

$$\begin{align*}{\varphi }_{1}^{A}& =\vert 1\rangle ,\\ {\varphi }_{2}^{A}& =\vert 2\rangle .\end{align*}$$

Discarding the unentangled occupied environmental orbital ${\varphi }_{3}^{B}$ that constitutes the vacuum state $\vert ~{0}\rangle ={\hat{\varphi }}_{3}^{B,{\dagger}}\vert 0\rangle$ of the Fock space $\mathcal{E}$ (see also section 2.2), we are left with the 4 orbitals of the CAS, i.e. ${\varphi }_{1}^{\text{CAS}}={\varphi }_{1}^{A}$, ${\varphi }_{2}^{\text{CAS}}={\varphi }_{2}^{A}$, ${\varphi }_{3}^{\text{CAS}}={\varphi }_{1}^{B}$ and ${\varphi }_{4}^{\text{CAS}}={\varphi }_{2}^{B}$. The corresponding 5 × 4 CAS matrix ${C}_{k\mu }^{\text{CAS}}\equiv {\varphi }_{\mu }^{\text{CAS}}(k)$ takes the form of equation (29). With this the embedded single-particle Hamiltonian becomes

$$\begin{align}{\mathbf{h}}_{\mathbf{s}}^{\prime }& ={[{\mathbf{C}}^{\text{CAS}}]}^{\mathrm{T}}{\mathbf{h}}_{\mathbf{s}}{\mathbf{C}}^{\text{CAS}}\\ & =\left(\begin{matrix}\hfill 0.0\hfill & \hfill -1.0\hfill & \hfill 0.0\hfill & \hfill 0.0\hfill \\ \hfill -1.0\hfill & \hfill 0.0\hfill & \hfill 0.451\enspace 56\hfill & \hfill -0.856\enspace 338\hfill \\ \hfill 0.0\hfill & \hfill 0.451\enspace 56\hfill & \hfill 1.408\enspace 29\hfill & \hfill -0.089\enspace 73\hfill \\ \hfill 0.0\hfill & \hfill -0.856\enspace 34\hfill & \hfill -0.089\enspace 73\hfill & \hfill -0.128\enspace 02\hfill \end{matrix}\right).\end{align} \tag{ A.62 }$$

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

$$\begin{align*}{{\epsilon}}_{1}^{\prime }& =-1.372\enspace 56\\ {\varphi }_{1}^{\prime }& =-0.513\enspace 87{\varphi }_{1}^{\text{CAS}}-0.705\enspace 32{\varphi }_{2}^{\text{CAS}}\\ & \quad +0.099\enspace 10{\varphi }_{3}^{\text{CAS}}-0.478\enspace 17{\varphi }_{4}^{\text{CAS}}\end{align*}$$

$$\begin{align*}{{\epsilon}}_{2}^{\prime }& =-0.079\enspace 22\\ {\varphi }_{2}^{\prime }& =0.634\enspace 51{\varphi }_{1}^{\text{CAS}}+0.050\enspace 27{\varphi }_{2}^{\text{CAS}}\\ & \quad -0.061\enspace 64{\varphi }_{3}^{\text{CAS}}-0.768\enspace 81{\varphi }_{4}^{\text{CAS}}\end{align*}$$

$$\begin{align*}{{\epsilon}}_{3}^{\prime }& =1.0\\ {\varphi }_{3}^{\prime }& =-0.5{\varphi }_{1}^{\text{CAS}}+0.5{\varphi }_{2}^{\text{CAS}}\\ & \quad -0.625\enspace 47{\varphi }_{3}^{\text{CAS}}-0.329\enspace 82{\varphi }_{4}^{\text{CAS}}\end{align*}$$

$$\begin{align*}{{\epsilon}}_{4}^{\prime }& =1.732\enspace 05\\ {\varphi }_{4}^{\prime }& =-0.288\enspace 68{\varphi }_{1}^{\text{CAS}}+0.5{\varphi }_{2}^{\text{CAS}}\\ & \quad +0.771\enspace 47{\varphi }_{3}^{\text{CAS}}-0.267\enspace 41{\varphi }_{4}^{\text{CAS}}.\end{align*}$$

Lifting the single-particle Hamiltonian to the Fock-space $\mathcal{E}$ we have (see also equation (15))

$$\begin{equation}{\hat{H}}^{\prime }=\sum\limits _{~{k}=1}^{4}\sum\limits _{{~{k}}^{\prime }=1}^{4}{h}_{s}^{\prime }(~{k},{~{k}}^{\prime }){\hat{\varphi }}_{~{k}}^{\text{CAS},{\dagger}}{\hat{\varphi }}_{{~{k}}^{\prime }}^{\text{CAS}}+\frac{{\Delta}{\epsilon}}{2}\sum\limits _{~{k}=1}^{4}{\hat{\varphi }}_{~{k}}^{\text{CAS},{\dagger}}{\hat{\varphi }}_{~{k}}^{\text{CAS}}\end{equation} \tag{ A.63 }$$

$$\begin{equation}=\sum\limits _{\mu =1}^{4}{{\epsilon}}_{\mu }^{\prime }{\hat{\varphi }}_{\mu }^{\prime {\dagger}}{\hat{\varphi }}_{\mu }^{\prime }+\frac{{\Delta}{\epsilon}}{2}\sum\limits _{\mu =1}^{4}{\hat{\varphi }}_{\mu }^{\prime {\dagger}}{\hat{\varphi }}_{\mu }^{\prime }.\end{equation} \tag{ A.64 }$$

Because we have discarded the unentangled occupied environmental orbital the sought-after ground state is given by the lowest two-particle eigenstate of ${\hat{H}}^{\prime }$ which leads to E' = ₁' + ₂' = −1.451 79 and $\vert {{\Phi}}^{\prime }\rangle ={\hat{\varphi }}_{1}^{\prime {\dagger}}{\hat{\varphi }}_{2}^{\prime {\dagger}}\vert ~{0}\rangle$. Because in our case ${\Delta}{\epsilon}=\langle ~{0}\vert \hat{H}~{0}\rangle =E-{E}^{\prime }=-1.280\enspace 26$. For the orbitals we find

$$\begin{align*}{\varphi }_{1}^{\prime }& =0.513\enspace 87\vert 1\rangle +0.705\enspace 32\vert 2\rangle \\ & \quad +0.454\enspace 22\vert 3\rangle +0.012\enspace 60\vert 4\rangle -0.178\enspace 85\vert 5\rangle ,\end{align*}$$

$$\begin{align*}{\varphi }_{2}^{\prime }& =0.634\enspace 51\vert 1\rangle +0.050\enspace 27\vert 2\rangle \\ & \quad -0.630\enspace 53\vert 3\rangle -0.176\enspace 73\vert 4\rangle +0.407\enspace 52\vert 5\rangle .\end{align*}$$

If we again disregard the unentangled occupied environmental orbital then the resulting Slater determinant is

$$\begin{equation}{~{{\Phi}}}^{\prime }(k,l)=\frac{1}{\sqrt{2}}\left({\varphi }_{1}^{\prime }(k){\varphi }_{2}^{\prime }(l)-{\varphi }_{1}^{\prime }(l){\varphi }_{2}^{\prime }(k)\right).\end{equation} \tag{ A.65 }$$

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 ${~{{\Phi}}}^{\prime }(1,2)=-0.298\enspace 187$ we can compare to, e.g. ${\sum }_{k=3}^{5}~{{\Phi}}(1,2,k)=0.683\enspace 012$ or some arbitrary combination with k of equation (A.2). However, if we also use that we know the discarded orbital, i.e. the form of the vacuum state $\vert ~{0}\rangle$, we find instead

$$\begin{align}~{{\Phi}}(i,j,k)& =\left(i,j,k\vert 1,2,3\rangle \right.\\ & =\frac{1}{\sqrt{3!}}\enspace \mathrm{det}\enspace \left\vert \begin{matrix}\hfill {\varphi }_{1}^{\prime }(i)\hfill & \hfill {\varphi }_{1}^{\prime }(j)\hfill & \hfill {\varphi }_{1}^{\prime }(k)\hfill \\ \hfill {\varphi }_{2}^{\prime }(i)\hfill & \hfill {\varphi }_{2}^{\prime }(j)\hfill & \hfill {\varphi }_{2}^{\prime }(k)\hfill \\ \hfill {\varphi }_{3}^{\prime }(i)\hfill & \hfill {\varphi }_{3}^{\prime }(j)\hfill & \hfill {\varphi }_{3}^{\prime }(k)\hfill \end{matrix}\right\vert \end{align} \tag{ A.66 }$$

and have recovered the full wave function.

Let us next see explicitly, why for an interacting system only the (less convenient) projection in Fock space works. The simplest wave function that exhibits the main feature of an interacting many-body wave function (multi-reference character) is the linear combination of two Slater determinants. Here we will use two Slater determinants build from orbitals of the non-interacting Hamiltonian of equation (A.1),

$$\begin{equation}\vert {{\Phi}}_{1}\rangle ={\hat{\phi }}_{1}^{{\dagger}}{\hat{\phi }}_{2}^{{\dagger}}{\hat{\phi }}_{3}^{{\dagger}}\vert \mathrm{v}\mathrm{a}\mathrm{c}\rangle ,\end{equation} \tag{ A.67 }$$

$$\begin{equation}\vert {{\Phi}}_{2}\rangle ={\hat{\phi }}_{1}^{{\dagger}}{\hat{\phi }}_{4}^{{\dagger}}{\hat{\phi }}_{5}^{{\dagger}}\vert \mathrm{v}\mathrm{a}\mathrm{c}\rangle .\end{equation} \tag{ A.68 }$$

Our model interacting wave function is then

$$\begin{equation}\vert {\Psi}\rangle ={\nu }_{1}\vert {{\Phi}}_{1}\rangle +{\nu }_{2}\vert {{\Phi}}_{2}\rangle \end{equation} \tag{ A.69 }$$

with (real) ${\nu }_{1}^{2}+{\nu }_{2}^{2}=1$. We can define, similar to the coefficient matrix C in equation (A.10), the coefficient matrix of the multi-determinant wave function as

$$\begin{equation}\mathbf{D}=\left(\begin{matrix}\hfill \frac{1}{\sqrt{12}}\hfill & \hfill -0.5\hfill & \hfill \frac{1}{\sqrt{3}}\hfill & \hfill -0.5\hfill & \hfill \frac{1}{\sqrt{12}}\hfill \\ \hfill 0.5\hfill & \hfill -0.5\hfill & \hfill 0\hfill & \hfill 0.5\hfill & \hfill -0.5\hfill \\ \hfill \frac{1}{\sqrt{3}}\hfill & \hfill 0\hfill & \hfill -\frac{1}{\sqrt{3}}\hfill & \hfill 0\hfill & \hfill \frac{1}{\sqrt{3}}\hfill \\ \hfill 0.5\hfill & \hfill 0.5\hfill & \hfill 0\hfill & \hfill -0.5\hfill & \hfill -0.5\hfill \\ \hfill \frac{1}{\sqrt{12}}\hfill & \hfill 0.5\hfill & \hfill \frac{1}{\sqrt{3}}\hfill & \hfill 0.5\hfill & \hfill \frac{1}{\sqrt{12}}\hfill \end{matrix}\right),\end{equation} \tag{ A.70 }$$

where the last two columns correspond to ϕ₄ and ϕ₅. With this we can determine the 1RDM of the interacting wave function

$$\begin{align}\gamma (i,j)& =\langle {\Psi}\vert {\hat{c}}_{j}^{{\dagger}}{\hat{c}}_{i}\vert {\Psi}\rangle \\ & ={\nu }_{1}^{2}\langle {{\Phi}}_{1}\vert {\hat{c}}_{j}^{{\dagger}}{\hat{c}}_{i}\vert {{\Phi}}_{1}\rangle +{\nu }_{2}^{2}\langle {{\Phi}}_{2}\vert {\hat{c}}_{j}^{{\dagger}}{\hat{c}}_{i}\vert {{\Phi}}_{2}\rangle \\ & ={\nu }_{1}^{2}\sum\limits _{\mu =1}^{3}{D}_{j,\mu }^{\mathrm{T}}{D}_{i,\mu }+{\nu }_{2}^{2}({D}_{j2}^{\mathrm{T}}{D}_{i2}+{D}_{j4}^{\mathrm{T}}{D}_{i4}+{D}_{j5}^{\mathrm{T}}{D}_{i5}).\end{align} \tag{ A.71 }$$

If we then fix the missing values, e.g. ν₁ = 0.8 and ν₂ = 0.6, we can calculate numerically the 1RDM and determine its environmental submatrix γ^env(i, j) which would correspond to i, j ∈ {3, 4, 5}. The eigenvalues and eigenvectors of γ^env(i, j) are

$$\begin{align*}{n}_{1}& =0.365\enspace 30\equiv {\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}^{2}\\ {\varphi }_{1}^{B}& =-0.451\enspace 53\vert 3\rangle +0.678\enspace 86\vert 4\rangle -0.579\enspace 02\vert 5\rangle \end{align*}$$

$$\begin{align*}{n}_{2}& =0.591\enspace 50\equiv {\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}^{2}\\ {\varphi }_{2}^{B}& =0.889\enspace 05\vert 3\rangle +0.287\enspace 34\vert 4\rangle -0.356\enspace 41\vert 5\rangle \end{align*}$$

$$\begin{align*}{n}_{3}& =0.816\enspace 53\equiv {\Vert}{~{\varphi }}_{3}{{\Vert}}_{B}^{2}\\ {\varphi }_{3}^{B}& =0.075\enspace 58\vert 3\rangle +0.675\enspace 70\vert 4\rangle +0.733\enspace 29\vert 5\rangle .\end{align*}$$

While before we did go on by discarding the orbital with n = 1, here we do not find such an unentangled occupied environmental orbital. Thus the procedure that uses the impurity submatrix does in general not work for interacting systems.

Here we show explicitly non-uniqueness of the approximate projection by constructing two non-interacting system that have the same impurity 1RDM ${\gamma }_{\text{imp}}^{s}(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:

$$\begin{equation}{\boldsymbol{h}}^{s}=\left(\begin{matrix}\hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ B.1 }$$

Diagonalizing it leads to six eigenstates {ϕ₁, ..., ϕ₆}, and choosing the three lowest eigenstates leads {ϕ₁, ϕ₂, ϕ₃} with ground-state energy E = −4. The 1RDM

$$\begin{equation}{\boldsymbol{\gamma }}^{s}=\left(\begin{matrix}\hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0\hfill & \hfill -0.166\enspace 66\hfill & \hfill 0\hfill & \hfill 0.333\enspace 33\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0\hfill & \hfill -0.166\enspace 66\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0\hfill & \hfill -0.166\enspace 66\hfill \\ \hfill -0.166\enspace 66\hfill & \hfill 0\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -0.166\enspace 666\hfill & \hfill 0\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill & \hfill 0.333\enspace 333\hfill \\ \hfill 0.333\enspace 333\hfill & \hfill 0\hfill & \hfill -0.166\enspace 666\hfill & \hfill 0\hfill & \hfill 0.333\enspace 333\hfill & \hfill 0.5\hfill \end{matrix}\right).\end{equation} \tag{ B.2 }$$

The resulting impurity 1RDM on A ≡ {1, 2} is then

$$\begin{equation}{\boldsymbol{\gamma }}_{\text{imp}}^{s}=\left(\begin{matrix}\hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill \end{matrix}\right).\end{equation} \tag{ B.3 }$$

To then construct a different system with the same impurity 1RDM as its three-particle ground state we first diagonalize ${\boldsymbol{\gamma }}_{\text{imp}}^{s}$ of equation (B.3) and get the eigenvalues and the eigenvectors of the impurity 1RDM as

$$\begin{align}\begin{aligned}{n}_{1}^{\text{imp}}& =0.166\enspace 66\equiv {\Vert}{~{\varphi }}_{1}{{\Vert}}_{A}^{2}\\ {\varphi }_{1}^{A}& =0.707\enspace 11\vert 1\rangle -0.707\enspace 11\vert 2\rangle .\end{aligned}\end{align} \tag{ B.4 }$$

$$\begin{align}\begin{aligned}{n}_{2}^{\text{imp}}& =0.833\enspace 33\equiv {\Vert}{~{\varphi }}_{2}{{\Vert}}_{A}^{2}\\ {\varphi }_{2}^{A}& =0.707\enspace 11\vert 1\rangle +0.707\enspace 11\vert 2\rangle \end{aligned}.\end{align} \tag{ B.5 }$$

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 equation (B.2) and use two of them to build our CAS space, i.e.

$$\begin{equation}{\boldsymbol{\gamma }}_{\text{env}}^{s}=\left(\begin{matrix}\hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0\hfill & \hfill -0.166\enspace 66\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill & \hfill 0.333\enspace 333\hfill \\ \hfill -0.166\enspace 666\hfill & \hfill 0\hfill & \hfill 0.333\enspace 333\hfill & \hfill 0.5\hfill \end{matrix}\right).\end{equation} \tag{ B.6 }$$

This leads to

$$\begin{equation}{n}_{1}=0\equiv {\Vert}{~{\varphi }}_{1}{{\Vert}}_{B}^{2}\end{equation} \tag{ B.7 }$$

$$\begin{equation}{\varphi }_{1}^{B}=-0.316\enspace 23\vert 3\rangle +0.632\enspace 46\vert 4\rangle -0.632\enspace 46\vert 5\rangle +0.316\enspace 23\vert 6\rangle \end{equation} \tag{ B.8 }$$

$$\begin{align}\begin{aligned}{n}_{2}& =0.166\enspace 66\equiv {\Vert}{~{\varphi }}_{2}{{\Vert}}_{B}^{2}\\ {\varphi }_{2}^{B}& =0.632\enspace 46\vert 3\rangle -0.316\enspace 23\vert 4\rangle -0.316\enspace 23\vert 5\rangle +0.632\enspace 46\vert 6\rangle \end{aligned}\end{align} \tag{ B.9 }$$

$$\begin{equation}{n}_{3}=0.833\enspace 33\equiv {\Vert}{~{\varphi }}_{3}{{\Vert}}_{B}^{2}\end{equation} \tag{ B.10 }$$

$$\begin{equation}{\varphi }_{3}^{B}=-0.632\enspace 46\vert 3\rangle -0.316\enspace 23\vert 4\rangle +0.316\enspace 23\vert 5\rangle +0.632\enspace 46\vert 6\rangle \end{equation} \tag{ B.11 }$$

$$\begin{align}\begin{aligned}{n}_{4}& =1.0\equiv {\Vert}{~{\varphi }}_{4}{{\Vert}}_{B}^{2}\\ {\varphi }_{4}^{B}& =0.316\enspace 23\vert 3\rangle +0.632\enspace 46\vert 4\rangle +0.632\enspace 46\vert 5\rangle +0.316\enspace 23\vert 6\rangle .\end{aligned}\end{align} \tag{ B.12 }$$

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

$$\begin{equation}{\varphi }_{1}^{\text{CAS}}=\vert 1\rangle ,\end{equation} \tag{ B.13 }$$

$$\begin{equation}{\varphi }_{2}^{\text{CAS}}=\vert 2\rangle .\end{equation} \tag{ B.14 }$$

Alternatively, we could have used the two orbitals ϕ^A of the impurity submatrix as we have discussed in the previous part of the appendix. 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

$$\begin{equation}{\varphi }_{3}^{\text{CAS}}={\varphi }_{2}^{B},\end{equation} \tag{ B.15 }$$

$$\begin{equation}{\varphi }_{4}^{\text{CAS}}={\varphi }_{3}^{B}.\end{equation} \tag{ B.16 }$$

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 ${\varphi }_{3}^{B}$ and ${\varphi }_{4}^{B}$ and define

$$\begin{equation}{~{\varphi }}_{1}^{\prime }=\sqrt{{n}_{1}^{\text{imp}}}{\varphi }_{1}^{A}+\underset{=\sqrt{0.833\enspace 33}}{\underbrace{\sqrt{1-{n}_{1}^{\text{imp}}}}}{\varphi }_{1}^{B}\end{equation} \tag{ B.17 }$$

$$\begin{equation}{~{\varphi }}_{2}^{\prime }=\sqrt{{n}_{2}^{\text{imp}}}{\varphi }_{2}^{A}+\underset{=\sqrt{0.166\enspace 66}}{\underbrace{\sqrt{1-{n}_{2}^{\text{imp}}}}}{\varphi }_{4}^{B}\end{equation} \tag{ B.18 }$$

as well as

$$\begin{equation}{~{\varphi }}_{3}^{\prime }={\varphi }_{2}^{B}\end{equation} \tag{ B.19 }$$

$$\begin{equation}{~{\varphi }}_{4}^{\prime }={\varphi }_{3}^{B}.\end{equation} \tag{ B.20 }$$

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 ${{\epsilon}}_{1}^{\prime }{\leqslant}{{\epsilon}}_{2}^{\prime }{\leqslant}{{\epsilon}}_{3}^{\prime }{< }{{\epsilon}}_{4}^{\prime }{\leqslant}{{\epsilon}}_{5}^{\prime }={{\epsilon}}_{6}^{\prime }=0$. For definiteness, we choose ₁' = −4, ${{\epsilon}}_{2}^{\prime }=-3$, ${{\epsilon}}_{3}^{\prime }=-2$, ${{\epsilon}}_{4}^{\prime }=-1$ and ${{\epsilon}}_{5}^{\prime }={{\epsilon}}_{6}^{\prime }=0$. With this the new Hamiltonian is

$$\begin{align}{\boldsymbol{h}}_{\text{new}}^{s}& =\sum\limits _{\mu =1}^{4}{{\epsilon}}_{\mu }^{\prime }{~{\varphi }}_{\mu }^{\prime }(i){~{\varphi }}_{\mu (j)}^{\prime {\ast}}\\ & =\left(\begin{matrix}\hfill -1.583\enspace 33\hfill & \hfill -0.916\enspace 67\hfill & \hfill -0.583\enspace 33\hfill & \hfill 0.166\enspace 67\hfill & \hfill -1.166\enspace 67\hfill & \hfill 0.083\enspace 33\hfill \\ \hfill -0.916\enspace 67\hfill & \hfill -1.583\enspace 33\hfill & \hfill 0.083\enspace 33\hfill & \hfill -1.166\enspace 67\hfill & \hfill 0.166\enspace 67\hfill & \hfill -0.583\enspace 33\hfill \\ \hfill -0.583\enspace 33\hfill & \hfill 0.083\enspace 33\hfill & \hfill -1.583\enspace 33\hfill & \hfill 0.766\enspace 67\hfill & \hfill -0.166\enspace 67\hfill & \hfill -0.116\enspace 67\hfill \\ \hfill 0.166\enspace 67\hfill & \hfill -1.166\enspace 67\hfill & \hfill 0.766\enspace 67\hfill & \hfill -1.833\enspace 33\hfill & \hfill 1.033\enspace 33\hfill & \hfill -0.166\enspace 67\hfill \\ \hfill -1.166\enspace 67\hfill & \hfill 0.166\enspace 67\hfill & \hfill -0.166\enspace 67\hfill & \hfill 1.033\enspace 33\hfill & \hfill -1.833\enspace 33\hfill & \hfill 0.766\enspace 67\hfill \\ \hfill 0.083\enspace 33\hfill & \hfill -0.583\enspace 33\hfill & \hfill -0.116\enspace 67\hfill & \hfill -0.166\enspace 67\hfill & \hfill 0.766\enspace 67\hfill & \hfill -1.583\enspace 33\hfill \end{matrix}\right).\end{align} \tag{ B.21 }$$

If we diagonalize the Hamiltonian and take the lowest three eigenvectors we find the three-particle ground-state 1RDM

$$\begin{equation}{\boldsymbol{\gamma }}_{\text{new}}^{s}=\left(\begin{matrix}\hfill 0.500\enspace 00\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.166\enspace 67\hfill & \hfill 0.000\enspace 00\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.000\enspace 00\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.500\enspace 00\hfill & \hfill 0.000\enspace 00\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.000\enspace 00\hfill & \hfill 0.166\enspace 67\hfill \\ \hfill 0.166\enspace 67\hfill & \hfill 0.000\enspace 00\hfill & \hfill 0.500\enspace 00\hfill & \hfill -0.333\enspace 33\hfill & \hfill 0.000\enspace 00\hfill & \hfill 0.333\enspace 33\hfill \\ \hfill 0.000\enspace 00\hfill & \hfill 0.333\enspace 33\hfill & \hfill -0.333\enspace 33\hfill & \hfill 0.500\enspace 00\hfill & \hfill -0.166\enspace 67\hfill & \hfill 0.000\enspace 00\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.000\enspace 00\hfill & \hfill 0.000\enspace 00\hfill & \hfill -0.166\enspace 67\hfill & \hfill 0.500\enspace 00\hfill & \hfill -0.333\enspace 33\hfill \\ \hfill 0.000\enspace 00\hfill & \hfill 0.166\enspace 67\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.000\enspace 00\hfill & \hfill -0.333\enspace 33\hfill & \hfill 0.500\enspace 00\hfill \end{matrix}\right)\end{equation} \tag{ B.22 }$$

and the corresponding ground-state energy is ${E}^{\prime }={{\epsilon}}_{1}^{\prime }+{{\epsilon}}_{2}^{\prime }+{{\epsilon}}_{3}^{\prime }=-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)

$$\begin{equation}{\varphi }_{3}^{\prime \mathrm{C}\mathrm{A}\mathrm{S}}={\varphi }_{4}^{B},\end{equation} \tag{ B.23 }$$

$$\begin{equation}{\varphi }_{4}^{\prime \mathrm{C}\mathrm{A}\mathrm{S}}={\varphi }_{1}^{B},\end{equation} \tag{ B.24 }$$

in contrast to the ones of equations (B.15) and (B.16). However, ${\varphi }_{4}^{B}$ and ${\varphi }_{1}^{B}$ were orbitals that did not belong to the original CAS. That means that also the projection constructed from the Hamiltonian ${\boldsymbol{h}}_{\text{new}}^{s}$ will look very different from the one of h ^s of equation (B.1). 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.

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

$$\begin{align}{\hat{\boldsymbol{h}}}^{\text{tar}}& \equiv \sum\limits _{\mu =1}^{6}{{\epsilon}}_{\mu }^{\text{tar}}{\phi }_{\mu }^{\mathrm{t}\mathrm{a}\mathrm{r}}(i){\phi }_{\mu }^{\text{tar}}(j)\\ & \equiv \left(\begin{matrix}0.916\enspace 67& -0.583\enspace 333& 0.166\enspace 667& -0.083\enspace 333& -0.583\enspace 333& -0.833\enspace 333\\ -0.583\enspace 333& 0.916\enspace 667& -0.833\enspace 333& -0.583\enspace 333& -0.083\enspace 333& 0.166\enspace 667\\ 0.166\enspace 667& -0.833\enspace 333& -0.233\enspace 333& -0.033\enspace 333& -0.633\enspace 333& 0.566\enspace 667\\ -0.083\enspace 333& -0.583\enspace 333& -0.033\enspace 333& -0.683\enspace 333& 1.016\enspace 667& -0.633\enspace 333\\ -0.583\enspace 333& -0.083\enspace 333& -0.633\enspace 333& 1.016\enspace 667& -0.683\enspace 333& -0.033\enspace 333\\ -0.833\enspace 333& 0.166\enspace 667& 0.566\enspace 667& -0.633\enspace 333& -0.033\enspace 333& -0.233\enspace 333\end{matrix}\right),\end{align} \tag{ B.25 }$$

where we have defined the orthogonal set of eigenfunctions as

$$\begin{align}\left\{{\phi }_{1}^{\text{tar}}\right.& =\left.{\varphi }_{1}^{B},{\phi }_{2}^{\text{tar}}={\phi }_{1},{\phi }_{3}^{\text{tar}}={\phi }_{2},{\phi }_{4}^{\text{tar}}={\phi }_{3},{\phi }_{5}^{\text{tar}}=\frac{1}{\vert {\phi }_{4}-\langle {\phi }_{1}^{\text{tar}}\vert {\phi }_{4}\rangle {\phi }_{1}^{\text{tar}}\vert }\right.\\ & \quad {\times}\left.({\phi }_{4}-\langle {\phi }_{1}^{\text{tar}}\vert {\phi }_{4}\rangle {\phi }_{1}^{\text{tar}}),\right.\\ {\phi }_{6}^{\text{tar}}& =\left.\frac{1}{\vert {\phi }_{5}-\langle {\phi }_{1}^{\text{tar}}\vert {\phi }_{5}\rangle {\phi }_{1}^{\text{tar}}\vert }({\phi }_{5}-\langle {\phi }_{1}^{\text{tar}}\vert {\phi }_{5}\rangle {\phi }_{1}^{\text{tar}})+\frac{1}{\vert {\phi }_{5}-\langle {\phi }_{5}^{\text{tar}}\vert {\phi }_{5}\rangle {\phi }_{5}^{\text{tar}}\vert }\right.\\ & \quad {\times}\left.({\phi }_{5}-\langle {\phi }_{5}^{\text{tar}}\vert {\phi }_{5}\rangle {\phi }_{5}^{\text{tar}})\right\}\end{align} \tag{ B.26 }$$

and {ϕ₁, ..., ϕ₅} are the five lowest eigenstates of equation (B.1). Further we have chosen

$$\begin{equation*}{{\epsilon}}_{1}^{\text{tar}}=-2,{{\epsilon}}_{2}^{\text{tar}}=-1,{{\epsilon}}_{3}^{\text{tar}}=-0.5,{{\epsilon}}_{4}^{\text{tar}}=0.5,{{\epsilon}}_{5}^{\text{tar}}=1,{{\epsilon}}_{6}^{\text{tar}}=2.\end{equation*}$$

For this Hamiltonian the three-particle ground-state energy is E^tar = −3.5 and the corresponding 1RDM is

$$\begin{equation}{\boldsymbol{\gamma }}^{\text{tar}}=\left(\begin{matrix}\hfill 0.497\enspace 96\hfill & \hfill 0.354\enspace 83\hfill & \hfill 0.023\enspace 54\hfill & \hfill -0.164\enspace 63\hfill & \hfill -0.021\enspace 50\hfill & \hfill 0.309\enspace 80\hfill \\ \hfill 0.354\enspace 83\hfill & \hfill 0.273\enspace 54\hfill & \hfill 0.085\enspace 37\hfill & \hfill -0.021\enspace 50\hfill & \hfill 0.059\enspace 80\hfill & \hfill 0.247\enspace 96\hfill \\ \hfill 0.023\enspace 54\hfill & \hfill 0.085\enspace 37\hfill & \hfill 0.328\enspace 50\hfill & \hfill 0.109\enspace 80\hfill & \hfill 0.447\enspace 96\hfill & \hfill 0.004\enspace 83\hfill \\ \hfill -0.164\enspace 63\hfill & \hfill -0.021\enspace 50\hfill & \hfill 0.109\enspace 80\hfill & \hfill 0.897\enspace 96\hfill & \hfill -0.045\enspace 17\hfill & \hfill 0.223\enspace 54\hfill \\ \hfill -0.021\enspace 50\hfill & \hfill 0.059\enspace 80\hfill & \hfill 0.447\enspace 96\hfill & \hfill -0.045\enspace 17\hfill & \hfill 0.673\enspace 54\hfill & \hfill -0.114\enspace 63\hfill \\ \hfill 0.309\enspace 80\hfill & \hfill 0.247\enspace 96\hfill & \hfill 0.004\enspace 83\hfill & \hfill 0.223\enspace 54\hfill & \hfill -0.114\enspace 63\hfill & \hfill 0.3285\hfill \end{matrix}\right)\end{equation} \tag{ B.27 }$$

with the impurity 1RDM as

$$\begin{equation}{\boldsymbol{\gamma }}_{\text{imp}}^{\text{tar}}=\left(\begin{matrix}\hfill 0.497\enspace 96\hfill & \hfill 0.354\enspace 83\hfill \\ \hfill 0.354\enspace 83\hfill & \hfill 0.273\enspace 54\hfill \end{matrix}\right).\end{equation} \tag{ B.28 }$$

Next we assume that a DMET iteration step led us to an auxiliary Hamiltonian of the form of equation (B.1). So we follow the DMET procedure and determine the CAS of this auxiliary Hamiltonian (see equations (B.13) to (B.16)) and define the embedded Hamiltonian

$$\begin{align}{\mathbf{h}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}& ={[{\mathbf{C}}^{\text{CAS}}]}^{\mathrm{T}}{\mathbf{h}}^{\text{tar}}{\mathbf{C}}^{\text{CAS}}\\ & \quad {\times}\left(\begin{matrix}\hfill 0.000\enspace 00\hfill & \hfill -1.0000\hfill & \hfill -0.632\enspace 46\hfill & \hfill 0.632\enspace 46\hfill \\ \hfill -1.000\enspace 00\hfill & \hfill 0.0000\hfill & \hfill -0.632\enspace 46\hfill & \hfill 0.632\enspace 46\hfill \\ \hfill -0.632\enspace 46\hfill & \hfill -0.632\enspace 46\hfill & \hfill 0.600\enspace 00\hfill & \hfill 0.000\enspace 00\hfill \\ \hfill 0.632\enspace 46\hfill & \hfill -0.632\enspace 46\hfill & \hfill 0.000\enspace 00\hfill & \hfill -0.600\enspace 00\hfill \end{matrix}\right),\end{align} \tag{ B.29 }$$

with C^CAS the 6 × 4 matrix constructed from these orbitals. Diagonalizing this Hamiltonian and keeping only the two lowest (embedded) orbitals we obtain an embedded 1RDM of the target Hamiltonian in the CAS basis as

$$\begin{align}{\boldsymbol{\gamma }}_{\text{CAS}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}& \equiv \sum\limits _{k=1}^{2}{\varphi }_{k}^{\text{emb}}(\mu ){\varphi }_{k}^{\text{emb}}(\nu )\\ & \equiv \left(\begin{matrix}\hfill 0.500\enspace 00\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.263\enspace 52\hfill & \hfill 0.263\enspace 52\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.500\enspace 00\hfill & \hfill 0.263\enspace 52\hfill & \hfill 0.263\enspace 52\hfill \\ \hfill 0.263\enspace 52\hfill & \hfill 0.263\enspace 52\hfill & \hfill 0.166\enspace 67\hfill & \hfill 0.000\enspace 00\hfill \\ \hfill -0.263\enspace 52\hfill & \hfill 0.263\enspace 52\hfill & \hfill 0.000\enspace 00\hfill & \hfill 0.833\enspace 33\hfill \end{matrix}\right),\end{align} \tag{ B.30 }$$

where ${\varphi }_{\mu }^{\text{emb}}$ are the two lowest eigenstates of equation (B.29). Transforming the 1RDM into the site basis by

$$\begin{equation}{\gamma }_{\text{emb}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}(i,j)=\sum\limits _{\mu ,\nu =1}^{4}{\gamma }_{\text{CAS}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}(\mu ,\nu ){\varphi }_{\mu }^{\text{CAS}}(i){\varphi }_{\nu }^{\text{CAS}}(j)\end{equation} \tag{ B.31 }$$

leads to

$$\begin{equation}{\boldsymbol{\gamma }}_{\text{emb}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}=\left(\begin{matrix}\hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.000\enspace 00\hfill & \hfill -0.166\enspace 67\hfill & \hfill 0.\hfill & \hfill 0.333\enspace 33\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.000\enspace 00\hfill & \hfill -0.166\enspace 67\hfill & \hfill 0.000\enspace 00\hfill \\ \hfill 0.000\enspace 00\hfill & \hfill 0.333\enspace 33\hfill & \hfill 0.400\enspace 00\hfill & \hfill 0.133\enspace 33\hfill & \hfill -0.200\enspace 00\hfill & \hfill -0.266\enspace 67\hfill \\ \hfill -0.166\enspace 67\hfill & \hfill 0.000\enspace 00\hfill & \hfill 0.133\enspace 33\hfill & \hfill 0.1\hfill & \hfill -0.066\enspace 67\hfill & \hfill -0.2\hfill \\ \hfill 0.\hfill & \hfill -0.166\enspace 67\hfill & \hfill -0.2\hfill & \hfill -0.066\enspace 67\hfill & \hfill 0.1\hfill & \hfill 0.133\enspace 33\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.000\enspace 00\hfill & \hfill -0.266\enspace 67\hfill & \hfill -0.200\enspace 00\hfill & \hfill 0.133\enspace 33\hfill & \hfill 0.400\enspace 00\hfill \end{matrix}\right).\end{equation} \tag{ B.32 }$$

We notice that the 1RDM ${\boldsymbol{\gamma }}_{\text{emb}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}$ constructed from the embedded Hamiltonian of equation (B.29) does not agree with the target 1RDM γ ^tar even on the impurity A. That is, the approximate impurity 1RDM is

$$\begin{equation}{\boldsymbol{\gamma }}_{\text{imp}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}=\left(\begin{matrix}\hfill 0.5\hfill & \hfill 0.333\enspace 33\hfill \\ \hfill 0.333\enspace 33\hfill & \hfill 0.5\hfill \end{matrix}\right),\end{equation} \tag{ B.33 }$$

while the 'exact' impurity 1RDM is given in equation (B.28). Yet it does agree with the 1RDM γ _imp of the auxiliary Hamiltonian of equation (B.1) on the impurity. So we have attained the convergence criterion

$$\begin{equation}{\boldsymbol{\gamma }}_{\text{imp}}^{s}={\boldsymbol{\gamma }}_{\text{imp}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}},\end{equation} \tag{ B.34 }$$

and thus our DMET iteration is finished. Besides that we find completely wrong 1RDMs, also the energy estimate is not necessarily good. To demonstrate this we are going to use the following formula to calculate first the energy of the fragment A:

$$\begin{equation}{{\epsilon}}_{\mathrm{f}}^{\text{exact}}=\sum\limits _{i=1,2,j=1-6}{\boldsymbol{h}}_{i,j}^{\text{tar}}{\gamma }^{\text{tar}}(j,i)=-0.55242,\end{equation} \tag{ B.35 }$$

where the expression for the fragment energy is taken from [37] (equation (25)). However, because in practice we do not have the correct 1RDM that corresponds to this Hamiltonian available we need to calculate the fragment energy using the embedded 1RDM:

$$\begin{equation}{{\epsilon}}_{\mathrm{f}}^{\text{emb}}=\sum\limits _{i=1,2,j=1-6}{\boldsymbol{h}}_{i,j}^{\text{tar}}{\gamma }_{\text{emb}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}(j,i)=-0.438\enspace 16.\end{equation} \tag{ B.36 }$$

Following the same procedure after adding to embedded 1RDM the environment orbital that we had originally discarded (so as to have three particles)

$$\begin{equation}{\gamma }_{\text{emb}}^{\mathrm{t}\mathrm{a}\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{t}}(j,i)={\boldsymbol{\gamma }}_{\text{emb}}^{\prime \mathrm{t}\mathrm{a}\mathrm{r}}(j,i)+{\phi }_{4}^{{\ast}B}(j){\phi }_{4}^{B}(i)\end{equation} \tag{ B.37 }$$

we obtain the same wrong fragment energy as in (B.36).

The reason why we can construct such a 'bad' fixed point is that we can instead of the ground state of a target Hamiltonian (in our case the three lowest orbitals of the Hamiltonian of equation (B.25)) end up in an excited state (in our case a Slater determinant that excludes the lowest-energy orbitals ${\phi }_{1}^{\text{tar}}$). 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 equations (B.13) to (B.16) is orthonormal to ${\phi }_{1}^{\text{tar}}\equiv {\varphi }_{3}^{B}$). 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 ${\phi }_{1}^{\text{tar}}\equiv {\varphi }_{3}^{B}$) 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.