# Householder-transformed density matrix functional embedding theory

Sajanthan Sekaran, Masahisa Tsuchiizu, M. Saubanère, Emmanuel Fromager

► To cite this version:

Sajanthan Sekaran, Masahisa Tsuchiizu, M. Saubanère, Emmanuel Fromager. Householder-transformed density matrix functional embedding theory. Physical Review B, 2021, 104 (3), 10.1103/PhysRevB.104.035121. [hal-03366562](https://hal.science/hal-03366562)

HAL Id: hal-03366562

https://hal.science/hal-03366562

Submitted on 5 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Householder transformed density matrix functional embedding theory

Sajjanthan Sekaran,¹ Masahisa Tsuchiizu,² Matthieu Saubanère,³ and Emmanuel Fromager¹

¹Laboratoire de Chimie Quantique, Institut de Chimie,
CNRS/Université de Strasbourg, 4 rue Blaise Pascal, 67000 Strasbourg, France

²Department of Physics, Nara Women's University, Nara 630-8506, Japan

³ICGM, Université de Montpellier, CNRS, ENSCM, Montpellier, France

Quantum embedding based on the (one-electron reduced) density matrix is revisited by means of the unitary Householder transformation. While being exact and equivalent to (but formally simpler than) density matrix embedding theory (DMET) in the non-interacting case, the resulting Householder transformed density matrix functional embedding theory (Ht-DMFET) preserves, by construction, the single-particle character of the bath when electron correlation is introduced. In Ht-DMFET, the projected “impurity+bath” cluster’s Hamiltonian (from which approximate local properties of the interacting lattice can be extracted) becomes an explicit functional of the density matrix. In the spirit of single-impurity DMET, we consider in this work a closed (two-electron) cluster constructed from the full-size non-interacting density matrix. When the (Householder transformed) interaction on the bath site is taken into account, per-site energies obtained for the half-filled one-dimensional Hubbard lattice match almost perfectly the exact Bethe Ansatz results in all correlation regimes. In the strongly correlated regime, the results deteriorate away from half-filling. This can be related to the electron number fluctuations in the (two-site) cluster which are not described neither in Ht-DMFET nor in regular DMET. As expected, the per-site energies dramatically improve when increasing the number of embedded impurities. Formal connections with density/density matrix functional theories have been briefly discussed and should be explored further. Work is currently in progress in this direction.

## I. INTRODUCTION

Quantum embedding has emerged over the last two decades as a viable strategy for modelling strong electron correlation in large molecules and extended systems [1]. The purpose of an embedding procedure is to replace the original full-size problem, for which an accurate solution to the Schrödinger equation is out of reach, by one or several simpler problems that preserve only a fragment of the original system. The fragment, which can be a single atomic site (often referred to as *impurity*) in a lattice, is embedded into a formal bath that is supposed to mimic the effects of the impurity’s environment. The mathematical construction of the bath depends on the choice of basic variable in the embedding procedure. Obviously, this choice is not unique, which explains the diversity of embedding schemes in the literature [1]. In the well-established dynamical mean-field theory (DMFT) [2–6], the so-called local Green function, which is evaluated on the impurity, is the quantity of interest. In this case, the non-interacting sites of the Anderson model (on which the Green function is mapped) represent the bath. In the simplified two-site formulation of DMFT [7], the latter reduces to a single site. Note that the fragmentation of a system, which is central in embedding calculations, allows for the combination of different electronic structure methods. We can refer, for example, to density functional theory (DFT) or the Green-function-based GW method on top of DMFT [8–14], but also self-energy embedding theory (SEET) [15–18] or the dynamical configuration interaction method [19]. Various combinations of DFT with many-body wave functions have also been explored

for model [20–23] and *ab initio* quantum chemical Hamiltonians [24–29], with the purpose of improving DFT in the description of strong electron correlation.

In recent years, *density matrix embedding theory* (DMET) [30–34] has attracted an increasing attention [35–39]. Its drastic simplification of the bath (when compared to DMFT) explains this success. Moreover, unlike DMFT, DMET is originally a frequency-independent (and therefore formally simpler) theory. Later on it has been extended to the description of dynamical properties [40] and non-equilibrium dynamics [41]. As shown in Refs. 39, 42–45, formal connections can be established between the two theories. Since, in DMET, the number of bath sites equals (at most) the number of impurity sites, the Schrödinger equation can be solved accurately (if not exactly) for the reduced-in-size “impurity+bath” system. Despite its name, DMET falls in many ways into the category of *wave function*-based methods. It is clear in its exact formulation, where the many-body bath states arise from the Schmidt decomposition of the full-system many-body wave function [30]. We note in passing that embedding theories similar in spirit to DMET but based on the exact factorization of the (electronic) many-body wave function have been proposed very recently [46, 47]. In standard implementations of DMET, the Schmidt decomposition is applied to a single Slater determinant Φ. In this case, the bath simply consists of effective sites (or orbitals) that can be determined numerically from the overlap matrix between the fragment and the occupied orbitals in Φ [48]. In other words, the quantum partitioning in standard DMET is done at the
single-particle level [49]. The (one-electron reduced) density matrix comes into play in the optimization of $\Phi$, through mapping constraints. As pointed out in previous works [36, 39, 42], representability issues may arise in this context, thus leading to difficulties in the numerical robustness of the method. Note that relaxed constraints have been used, like in density embedding theory (DET) [35], where only site occupations (i.e., diagonal elements of the density matrix) are mapped. A connection with DFT can be established in this case [50].

The performance of DMET can in principle be improved systematically by incorporating correlation into the bath [30, 36, 51, 52]. The latter is then described by many-body states. Nevertheless, the numerical efficiency and formal beauty of standard (mean-field-based) DMET calculations definitely rely on the one-electron character of the bath. We may actually wonder if DMET can be made formally exact and systematically improvable while preserving a single-particle quantum partitioning picture. The recently proposed energy-weighted DMET (EwDMET) [39, 44, 45] is an elegant way to incorporate quantum fluctuations through a (truncated) description of the effective dynamics. We want to follow a different path where, ideally, the fully static character of DMET would be preserved. The present work is a first step toward this goal. More precisely, we rewrite the embedding as a functional of the density matrix, thus bypassing the Schmidt decomposition of the reference (correlated or not) full-system wave function. This is achieved via a unitary Householder transformation [53]. In this approach, that we refer to as Householder transformed density-matrix functional embedding theory (Ht-DMFET), the bath orbital becomes a simple and analytical functional of the (possibly correlated) density matrix, thus greatly simplifying its construction, even in the commonly used non-interacting (or mean-field) case. In this work, we describe in detail the embedding of a single impurity in the one-dimensional (1D) Hubbard model. Note that we also implemented a multiple-impurity version of Ht-DMFET based on a block Householder transformation [54]. Its detailed description will be given in a separate paper.

The present paper is organized as follows. After introducing the Householder transformation in Sec. II A, we show how it can be applied to the density matrix in order to derive an alternative embedding procedure (Sec. II B). A comparison with the Schmidt decomposition is made in Sec. II C. The introduction into the theory of a correlation potential and related embedding constraints is also discussed (in Sec. II D). Exact and approximate formulations of Ht-DMFET are detailed in non-interacting (Sec. III A) and interacting (Sec. III B) cases, respectively. The present (single-shot) implementation is summarized in Sec. III C. Its connection to DMET is also discussed. Results obtained for the 1D Hubbard model are presented in Sec. IV. Conclusions and perspectives

are finally given in Sec. V.

## II. THEORY

For simplicity, embedding strategies will be discussed in the following for the 1D Hubbard Hamiltonian,

$$ \hat{H} = -t \sum_{\sigma=\uparrow,\downarrow} \sum_{i=0}^{L-1} (\hat{c}_{i\sigma}^{\dagger} \hat{c}_{(i+1)\sigma} + \text{H.c.}) + U \sum_{i=0}^{L-1} \hat{n}_{i\uparrow} \hat{n}_{i\downarrow}, \quad (1) $$

where $\hat{n}_{i\sigma} = \hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}$, $\hat{c}_{L\sigma}^{\dagger} \equiv \hat{c}_{0\sigma}^{\dagger} (-\hat{c}_{0\sigma}^{\dagger})$ when periodic (antiperiodic) boundary conditions are used, $t$ is the hopping parameter, and $U$ is the on-site repulsion parameter. A single impurity site (labelled as $i=0$) will be embedded. Extensions to more general (higher-dimension or quantum chemical) Hamiltonians with multiple impurities are left for future work.

### A. Householder transformation

Let us consider two different column vectors $X$ and $Y$ with the same norm

$$ |X| = \sqrt{X^\dagger X} = \sqrt{Y^\dagger Y} = |Y|. \quad (2) $$

What we will refer to as Householder transformation [53] in the following simply corresponds to the reflection that transforms $X$ into $Y$. As shown in Appendix A and sketched in Fig. 1, its matrix representation $P$, which fulfills

$$ Y = PX, \quad (3) $$

can be written explicitly as follows,

$$ P = I - 2vv^{\dagger}, \quad (4) $$

where $I$ is the identity matrix and

$$ v = \frac{X - Y}{|X - Y|} \quad (5) $$

is the normalized Householder vector. Note that $P$ is both unitary and hermitian. Note that real algebra will be used throughout the paper. As shown in Fig. 1, the Householder vector is, by construction, orthogonal to the reflection's hyperplane.

Let us now connect the Householder transformation to density matrix embedding. Following the basic idea of DMET [30], we select an atomic site (say $i=0$), referred to as impurity, and reduce drastically the size of its environment in the lattice, in order to describe electron correlation. For that purpose, we use the density matrix elements that connect the impurity to the other lattice sites. All these elements will be collected into the above-mentioned column vector $X$, i.e.,

$$ X^{\dagger} = [\gamma_{00}, \gamma_{10}, \dots, \gamma_{i0}, \dots], \quad (6) $$
where, as we focus here on singlet ground states, the (spin) density matrix elements are denoted without spin indices:

$$ \gamma_{ij} = \langle \hat{c}_{i\uparrow}^{\dagger} \hat{c}_{j\uparrow} \rangle = \langle \hat{c}_{i\downarrow}^{\dagger} \hat{c}_{j\downarrow} \rangle. \quad (7) $$

In the reduced-in-size system, which we refer to as Householder cluster, we would like to preserve the (spin up or down) occupation $\gamma_{00}$ of the impurity site. The simplest cluster we can think of is a dimer consisting of the impurity and a bath site that needs to be determined. Ideally, we would like the cluster to be, within the density matrix, disconnected from its environment (whether this can actually be achieved or not will be discussed later on). On that basis, we define our Householder transformation-based embedding as follows:

$$ \mathbf{X}^{\dagger} \xrightarrow{\mathbf{P}} \mathbf{Y}^{\dagger} = [\gamma_{00}, \xi, 0, 0, \dots, 0, \dots], \quad (8) $$

where all but the first two rows of $\mathbf{Y}$ are set to zero and the constraint of Eq. (2) reads

$$ \xi^2 = \sum_{j>0} \gamma_{j0}^2. \quad (9) $$

As a result,

$$ \begin{aligned} |\mathbf{X} - \mathbf{Y}|^2 &= (\gamma_{10} - \xi)^2 + \sum_{j>1} \gamma_{j0}^2 \\ &= (\gamma_{10} - \xi)^2 + \xi^2 - \gamma_{10}^2, \\ &= 2\xi(\xi - \gamma_{10}). \end{aligned} \quad (10) $$

Obviously, $\xi \neq 0$, otherwise the impurity would be disconnected from the rest of the lattice (and there would be no need for an embedding). In the special case where $\gamma_{j0} \stackrel{j \ge 1}{=} 0$, which means that the impurity's nearest-neighboring site $i=1$ is the bath, the transformation should still be defined, meaning that $|\mathbf{X} - \mathbf{Y}|$ should never vanish. The latter requirement fixes the sign of $\xi$, thus leading to the final expressions [see Eqs. (4), (5), (6), and (8)]:

$$ P_{ij} = \delta_{ij} - 2v_i v_j, \quad (11) $$

where

$$ \begin{aligned} v_0 &= 0, \\ v_1 &= \frac{\gamma_{10} - \xi}{\sqrt{2\xi(\xi - \gamma_{10})}}, \\ v_j &= \frac{\gamma_{j0}}{\sqrt{2\xi(\xi - \gamma_{10})}}, \end{aligned} \quad (12) $$

and

$$ \xi = -\operatorname{sgn}(\gamma_{10}) \sqrt{\sum_{j>0} \gamma_{j0}^2}. \quad (13) $$

<figure><img src="image_2.png" /><figcaption>FIG. 1: Geometrical interpretation of the Householder transformation.</figcaption></figure>

## B. Householder transformed density matrix

As sketched in Fig. 2, the Householder transformed density matrix can now be evaluated as follows,

$$ \tilde{\gamma} = \mathbf{P}^{\dagger}\gamma\mathbf{P} = \mathbf{P}\gamma\mathbf{P}, \quad (14) $$

or, equivalently,

$$ \tilde{\gamma}_{ij} = \langle \hat{d}_{i\sigma}^{\dagger} \hat{d}_{j\sigma} \rangle, \quad (15) $$

where, according to Eqs. (11), (12), and (13), the Householder transformed creation (annihilation) operators are functionals of the density matrix elements $\gamma_{ij}$:

$$ \begin{aligned} \hat{d}_{i\sigma}^{\dagger} &:= \sum_{j} P_{ij} \hat{c}_{j\sigma}^{\dagger} \\ &= \hat{c}_{i\sigma}^{\dagger} - 2v_i \sum_{j>0} v_j \hat{c}_{j\sigma}^{\dagger}. \end{aligned} \quad (16) $$

We stress that the transformation leaves the impurity unchanged,

$$ \hat{d}_{\text{imp}}^{\dagger} \equiv \hat{d}_{0\sigma}^{\dagger} = \hat{c}_{0\sigma}^{\dagger}, \quad (17) $$

while the bath orbital is constructed explicitly from the environment (in the lattice) of the impurity as follows,

$$ \hat{d}_{\text{bath}}^{\dagger} \equiv \hat{d}_{1\sigma}^{\dagger} = (1 - 2v_1^2) \hat{c}_{1\sigma}^{\dagger} - 2v_1 \sum_{j>1} v_j \hat{c}_{j\sigma}^{\dagger}. \quad (18) $$

As shown in Fig. 3, for small Hubbard rings, the nearest neighbors of the impurity contribute the most to the bath. Nevertheless, the delocalization of the latter over the lattice can be quite substantial, in particular in the weakly correlated regime. The impact of correlation on the bath varies with the lattice filling. For example, in the quarter-filled 12-site ring, the deviation from the non-interacting bath orbital remains relatively small
<figure><img src="image_3.png" /><figcaption>FIG. 2: Schematics of the Householder transformation applied to the density matrix.</figcaption></figure>

when entering the strongly correlated regime (see the top panel of Fig. 3). This is an important observation which makes the use of a mean-field bath in conventional DMET calculations [48] relevant. As expected from Eq. (12) and Appendix B [see Eq. (B.16)], in the half-filled case, only the lattice sites with odd indices contribute to the bath (see the bottom panel of Fig. 3). Interestingly, in this case, the bath delocalization reduces as correlation increases.

Returning to the general theory, we note that the inverse transformation (from the Householder representation to the lattice one) simply reads

$$ \sum_i P_{ki} \hat{d}_{i\sigma}^{\dagger} = \sum_{ij} P_{ki} P_{ij} \hat{c}_{j\sigma}^{\dagger} = \sum_{j} \delta_{kj} \hat{c}_{j\sigma}^{\dagger} = \hat{c}_{k\sigma}^{\dagger}, \quad (19) $$

or, equivalently,

$$ \hat{c}_{k\sigma}^{\dagger} = \hat{d}_{k\sigma}^{\dagger} - 2v_k \sum_{i>0} v_i \hat{d}_{i\sigma}^{\dagger}. \quad (20) $$

By construction, the first column (row) of the Householder transformed density matrix equals zero outside the cluster [see Fig. 2]. Indeed, according to Eqs. (8) and

$$ \begin{aligned} \tilde{\gamma}_{j0} &= \sum_{kl} P_{jk} \gamma_{kl} P_{l0} \\ &= \sum_{k} P_{jk} \gamma_{k0} \\ &= [\mathbf{P} \mathbf{X}]_j \\ &= Y_j \stackrel{j \ge 2}{=} 0. \end{aligned} \quad (21) $$

Moreover, the occupation of the impurity is invariant under the Householder transformation:

$$ \tilde{\gamma}_{00} = \sum_{kl} P_{0k} \gamma_{kl} P_{l0} = \gamma_{00}. \quad (22) $$

In the general (interacting) case, the exact Householder cluster is not disconnected from its environment. As shown in Figs. 4 and 5 for small Hubbard rings, the Householder transformed (ground-state) density matrix has a nonzero buffer sector $\{\tilde{\gamma}_{j1}\}_{j \ge 2}$, which is a signature of the cluster's entanglement with its environment. The existence of such a buffer can be related to the deviation of the density matrix from idempotency. Indeed, according to Eq. (21),

$$ \begin{aligned} [\tilde{\gamma}^2]_{j0} &= \sum_k \tilde{\gamma}_{jk} \tilde{\gamma}_{k0} = \tilde{\gamma}_{j0} \tilde{\gamma}_{00} + \tilde{\gamma}_{j1} \tilde{\gamma}_{10} \\ \stackrel{j \ge 2}{=} \tilde{\gamma}_{j1} \tilde{\gamma}_{10}, \end{aligned} \quad (23) $$
<figure><img src="image_4.png" /><figcaption>FIG. 3: Expansion of the exact (ground-state) bath orbital on the lattice [see Eq. (18)] for quarter-filled 12-site (top panel) and half-filled 10-site (bottom panel) Hubbard rings in various correlation regimes.</figcaption></figure>

thus leading to [see Eq.(8)]

$$ \tilde{\gamma}_{j1} \stackrel{j \ge 2}{\equiv} \frac{[\tilde{\gamma}^2]_{j0}}{\xi} \qquad (24) $$

or, equivalently,

$$ \tilde{\gamma}_{j1} \stackrel{j \ge 2}{\equiv} \frac{[\tilde{\gamma}^2 - \tilde{\gamma}]_{j0}}{\xi}. \qquad (25) $$

Another important observation [see Fig. 4] is that the Householder cluster is, in general, an open subsystem. This can also be related to the deviation from idempotency. Indeed, by considering the particular case $j = 1$ on the first line of Eq. (23), it comes

$$ \frac{\mathcal{N}^c}{2} := \tilde{\gamma}_{00} + \tilde{\gamma}_{11} = \frac{[\tilde{\gamma}^2]_{10}}{[\tilde{\gamma}]_{10}} = \frac{[\tilde{\gamma}^2]_{10}}{\xi}, \qquad (26) $$

where $\mathcal{N}^c$ is the total (spin-summed) number of electrons in the cluster. Interestingly, when the interacting lattice is half-filled, the fluctuations in the number of electrons

<figure><img src="image_5.png" /><figcaption>FIG. 4: (Top) Exact (ground-state) Householder transformed density matrix elements in the cluster sector and root mean square (RMS) of the elements in the buffer sector, both plotted as functions of $U/t$ for a quarter-filled 12-site Hubbard ring. (Bottom) Individual elements in the buffer sector for various $U/t$ values.</figcaption></figure>

within the Householder cluster vanish and the latter contains exactly two electrons for all $U/t$ values [see Fig. 5]. Moreover, in the buffer sector of the density matrix, elements with odd row (column) indices are zero [see the bottom panel of Fig. 5]. Nevertheless, even in this particular case, the cluster remains connected to its Householder environment as long as the lattice is interacting ($U/t \neq 0$). As proved analytically in Appendix B, these properties originate from the hole-particle symmetry of the Hubbard Hamiltonian.

### C. Comparison with the Schmidt decomposition

Let us consider the determinant expansion of any (correlated or not) wave function $\Psi$ in the Householder rep-
<figure><img src="image_6.png" /><figcaption>FIG. 5: Same as Fig. 4 for a half-filled 10-site Hubbard ring.</figcaption></figure>

resentation,

$$|\Psi\rangle = \sum_{\overline{\mathcal{I}\mathcal{B}}} \sum_{\mathcal{E} \in \overline{\mathcal{I}\mathcal{B}}} \lambda_{\mathcal{I}\mathcal{B}\mathcal{E}} |\mathcal{I}\mathcal{B}\mathcal{E}\rangle, \quad (27)$$

where $|\mathcal{I}\mathcal{B}\mathcal{E}\rangle = |\mathcal{I}\rangle|\mathcal{B}\rangle|\mathcal{E}\rangle \equiv \hat{\mathcal{I}}^{\dagger}\hat{\mathcal{B}}^{\dagger}\hat{\mathcal{E}}^{\dagger}|\text{vac}\rangle$ consists of impurity, bath, and Householder environment states, i.e., in a second-quantized language,

$$ \begin{aligned} \hat{\mathcal{I}}^{\dagger} &= (\hat{d}_{0\uparrow}^{\dagger})^{n_{0\uparrow}^{\mathcal{I}}} (\hat{d}_{0\downarrow}^{\dagger})^{n_{0\downarrow}^{\mathcal{I}}}, \\ \hat{\mathcal{B}}^{\dagger} &= (\hat{d}_{1\uparrow}^{\dagger})^{n_{1\uparrow}^{\mathcal{B}}} (\hat{d}_{1\downarrow}^{\dagger})^{n_{1\downarrow}^{\mathcal{B}}}, \\ \hat{\mathcal{E}}^{\dagger} &= (\hat{d}_{2\uparrow}^{\dagger})^{n_{2\uparrow}^{\mathcal{E}}} (\hat{d}_{2\downarrow}^{\dagger})^{n_{2\downarrow}^{\mathcal{E}}} \dots (\hat{d}_{(L-1)\uparrow}^{\dagger})^{n_{(L-1)\uparrow}^{\mathcal{E}}} (\hat{d}_{(L-1)\downarrow}^{\dagger})^{n_{(L-1)\downarrow}^{\mathcal{E}}}, \end{aligned} \quad (28) $$

$n_{i\sigma}^{\mathcal{I},\mathcal{B},\mathcal{E}} \in \{0,1\}$ being the occupation numbers. The notation “$\mathcal{E} \in \overline{\mathcal{I}\mathcal{B}}$” means that the summation over the environment states is restricted to those preserving the total number of electrons in the lattice. In order to evaluate the density matrix elements $\tilde{\gamma}_{j0} = \langle \Psi | \hat{d}_{j\sigma}^{\dagger} \hat{d}_{0\sigma} | \Psi \rangle$, we need to look at the following overlaps:

$$\langle \mathcal{I}' \mathcal{B}' \mathcal{E}' | \hat{d}_{j\sigma}^{\dagger} \hat{d}_{0\sigma} | \mathcal{I}\mathcal{B}\mathcal{E} \rangle \xrightarrow{j \ge 2} \delta_{\mathcal{B}'\mathcal{B}} \langle \mathcal{I}' | \hat{d}_{0\sigma} | \mathcal{I} \rangle \langle \mathcal{E}' | \hat{d}_{j\sigma}^{\dagger} | \mathcal{E} \rangle. \quad (29)$$

Those will contribute if the wave function exhibits, in its determinant expansion, an excitation from the impurity to the environment. By imposing the constraint in Eq. (21) through the Householder transformation, we aim at removing such excitations, thus forcing the impurity to exchange electrons only with the bath spin-orbital $d_{1\sigma}$. The corresponding wave function can be expanded as follows,

$$|\Psi\rangle = \sum_{\overline{\mathcal{B}}} |\mathcal{I}^{\mathcal{B}}\mathcal{B}\rangle \sum_{\mathcal{E}^{\mathcal{B}}} \lambda_{\mathcal{I}^{\mathcal{B}}\mathcal{B}\mathcal{E}^{\mathcal{B}}} |\mathcal{E}^{\mathcal{B}}\rangle, \quad (30)$$

where the occupation of the bath fixes the occupation of the impurity (hence the notation $\mathcal{I}^{\mathcal{B}}$) and the number of electrons in the environment. Eq. (30) can be seen as a partial Schmidt decomposition [30]. The latter would be complete if, for each bath state $\mathcal{B}$, the environment were described by a single quantum state, i.e., $\sum_{\mathcal{E}^{\mathcal{B}}} \lambda_{\mathcal{I}^{\mathcal{B}}\mathcal{B}\mathcal{E}^{\mathcal{B}}} |\mathcal{E}^{\mathcal{B}}\rangle \rightarrow \lambda_{\mathcal{B}} |\Psi_{\mathcal{E}}^{\mathcal{B}}\rangle$. Note that, in the present formalism, the bath refers to a one-electron quantum state, even for a correlated wave function. In exact DMET, the bath states are many-body quantum states [30].

As mentioned previously, in the general case, the buffer sector $\{\tilde{\gamma}_{j1}\}_{j \ge 2}$ is nonzero. Since

$$\langle \mathcal{I}' \mathcal{B}' \mathcal{E}' | \hat{d}_{j\sigma}^{\dagger} \hat{d}_{1\sigma} | \mathcal{I}\mathcal{B}\mathcal{E} \rangle \xrightarrow{j \ge 2} \delta_{\mathcal{I}'\mathcal{I}} \langle \mathcal{B}' | \hat{d}_{1\sigma} | \mathcal{B} \rangle \langle \mathcal{E}' | \hat{d}_{j\sigma}^{\dagger} | \mathcal{E} \rangle, \quad (31)$$

it means that the bath can in principle exchange electrons with the Householder environment. In the particular case where the lattice wave function $\Psi$ can be rewritten in a single-Slater-determinant form $\Phi$, which is standard in practical DMET calculations [48], the density matrix becomes idempotent and the buffer sector vanishes [see Eq. (25)]. Therefore, in this case, the cluster and its environment become disentangled, as illustrated in Fig. 2. Moreover, the cluster contains exactly two electrons [see Eq. (26)]. As a result, in the determinant expansion of Eq. (30), each contributing cluster determinant $|\mathcal{I}^{\mathcal{B}}\mathcal{B}\rangle$ can be determined solely from the occupation of the impurity (because the latter gives the occupation of the bath). By diagonalizing the (idempotent) Householder environment block of the density matrix, we obtain a single Slater determinant $\Phi_{\mathcal{E}}$ from which the full-system wave function can be determined as follows,

$$ \begin{aligned} |\Phi\rangle &= \left( \sum_{\mathcal{I}} \lambda_{\mathcal{I}} |\mathcal{I}\mathcal{B}^{\mathcal{I}}\rangle \right) |\Phi_{\mathcal{E}}\rangle \\ &= \sum_{\mathcal{I}} \lambda_{\mathcal{I}} |\mathcal{I}\rangle |\mathcal{B}^{\mathcal{I}}\Phi_{\mathcal{E}}\rangle, \end{aligned} \quad (32) $$

which is formally identical to a Schmidt decomposition where, in the DMET terminology, the many-body bath states are $|\mathcal{B}^{\mathcal{I}}\Phi_{\mathcal{E}}\rangle$. The occupied orbitals in $\Phi_{\mathcal{E}}$ are usually referred to as core embedding orbitals [41, 48]. Note that, if we diagonalize the Householder (two-electron and idempotent) cluster block of the density matrix,

$$\tilde{\gamma}^C = \begin{bmatrix} \gamma_{00} & \xi \\ \xi & 1 - \gamma_{00} \end{bmatrix}, \quad (33)$$
where the idempotency constraint reads

$$ \xi^2 = \gamma_{00} (1 - \gamma_{00}), \quad (34) $$

we will obtain two spin-orbitals. The first one,

$$ |\mu_{\sigma}\rangle = \frac{\xi}{\sqrt{1 - \gamma_{00}}} |c_{0\sigma}\rangle + \sqrt{1 - \gamma_{00}} |d_{1\sigma}\rangle, \quad (35) $$

will be occupied, while its orthogonal counterpart (within the cluster) will be unoccupied. In this particular case, the cluster wave function can be written in a single-Slater-determinant form,

$$ \sum_{\mathcal{I}} \lambda_{\mathcal{I}} |\mathcal{I}B^{\mathcal{I}}\rangle \equiv \prod_{\sigma=\uparrow,\downarrow} \hat{\mu}_{\sigma}^{\dagger} |\text{vac}\rangle = |\Phi^{\mathcal{C}}\rangle, \quad (36) $$

exactly like in standard (approximate) DMET [see Eq. (9) in Ref. 48]. In order to make the connection with DMET even more explicit, we can rewrite Eq. (35) as follows,

$$ \begin{aligned} |\mu_{\sigma}\rangle &= \langle d_{0\sigma}|\mu_{\sigma}\rangle |d_{0\sigma}\rangle + \langle d_{1\sigma}|\mu_{\sigma}\rangle |d_{1\sigma}\rangle \\ &= \langle c_{0\sigma}|\mu_{\sigma}\rangle |c_{0\sigma}\rangle + \langle d_{1\sigma}|\mu_{\sigma}\rangle |d_{1\sigma}\rangle, \end{aligned} \quad (37) $$

thus leading to

$$ \begin{aligned} |d_{1\sigma}\rangle &= \frac{1}{\sqrt{1 - \langle c_{0\sigma}|\mu_{\sigma}\rangle^2}} (|\mu_{\sigma}\rangle - \langle c_{0\sigma}|\mu_{\sigma}\rangle |c_{0\sigma}\rangle) \\ &= \frac{1}{\sqrt{1 - \langle c_{0\sigma}|\mu_{\sigma}\rangle^2}} \sum_{j \ge 1} \langle c_{j\sigma}|\mu_{\sigma}\rangle |c_{j\sigma}\rangle. \end{aligned} \quad (38) $$

If we now introduce, like in DMET, the overlap matrix between all the occupied spin-orbitals $\nu_{\sigma}$ in $\Phi$ ($\mu_{\sigma}$ being one of them),

$$ S_{\nu\nu'} = \langle \nu_{\sigma} | c_{0\sigma} \rangle \langle c_{0\sigma} | \nu_{\sigma'} \rangle, \quad (39) $$

it becomes clear that the Householder transformation is equivalent to (but simpler than) the Schmidt decomposition in the particular case of single determinant wave functions. Indeed, by construction, $\mu_{\sigma}$ is the only occupied spin-orbital that actually overlaps with the impurity, thus leading to the simplified expression

$$ S_{\nu\nu'} = V_{\nu}\lambda^2 V_{\nu'}, \quad (40) $$

where $V_{\nu} = \delta_{\nu\mu}$ and $\lambda = \langle c_{0\sigma}|\mu_{\sigma}\rangle$. Thus, we recover from Eq. (38) the DMET expression for the bath [see Eq. (11) in Ref. 48],

$$ |d_{1\sigma}\rangle = \sum_{j \ge 1} \sum_{\nu}^{\text{occ.}} \frac{\langle c_{j\sigma}|\nu_{\sigma}\rangle V_{\nu}}{\sqrt{1 - \lambda^2}} |c_{j\sigma}\rangle, \quad (41) $$

where, according to Eq. (35),

$$ \lambda^2 = \frac{\xi^2}{1 - \gamma_{00}} = \gamma_{00}. \quad (42) $$

In summary, when the lattice is described with a single Slater determinant, the Householder transformation and the Schmidt decomposition will lead to the same block-diagonalized density matrix, each block corresponding

to the Householder cluster and its environment, respectively. However, the Householder transformation is much simpler in this case since the block-diagonalized structure is recovered automatically (via an analytical unitary transformation) from the density matrix of the lattice.

Obviously, the simplification in Eq. (36) does not hold anymore when electron correlation is introduced within the cluster. Nevertheless, in this case, the (disentangled) wave function structure shown in Eq. (32) will be preserved and denoted later in Sec. III B as follows:

$$ |\Phi\rangle = |\Phi^C\rangle |\Phi_{\mathcal{E}}\rangle \rightarrow |\Psi^C\rangle |\Phi_{\mathcal{E}}\rangle. \quad (43) $$

Note that, as shown in Appendix C, changes in the Householder cluster (and/or environment) sector of the density matrix have no direct impact on the direction of the Householder vector and therefore no impact on the bath. In other words, even though correlation within the cluster and/or the environment modifies the density matrix in the lattice representation, the correlated and original (idempotent) density matrices will still share the same Householder transformation [defined according to Eqs. (11) and (12)] where the cluster remains disconnected from its environment.

### D. Embedding constraints, correlation potential, and connection to DFT

In practical (multiple-impurity [39]) DMET calculations, the bath is indirectly affected by the correlation within the cluster through ad hoc density matrix mapping constraints which can modify, via a (possibly nonlocal) correlation potential operator acting on the full lattice [48, 55], the reference Slater determinant $\Phi$ of Eq. (32). In the present formalism, it would induce a change in the Householder transformation and, consequently, it would make the embedding procedure self-consistent.

As stressed in previous works [36, 42], such a mapping cannot be exact, simply because, unlike a correlated density matrix, the reference (non-interacting) one is idempotent. From a formal point of view, rewriting the embedding procedure in the language of DFT, like in DET [35] or in the self-consistent density-functional embedding (SDE) approach [50], is more appealing because the (weaker) density-only constraint can be fulfilled, in principle exactly. More precisely, if we combine Kohn-Sham (KS) DFT [56] with a quasi-degenerate perturbation theory-like formalism [57, 58], the true interacting ground-state wave function $\Psi_0$ can be described by its projection onto a so-called model Hilbert space $\mathcal{M}$. The latter would be constructed from the (non-interacting KS) Householder cluster many-body states and the core embedding KS orbitals. While the (ground-state) KS determinant $|\Phi_{\text{KS}}\rangle \equiv |\Phi_{\text{KS}}^{\mathcal{C}}\rangle |\Phi_{\mathcal{E}}\rangle$ belongs to $\mathcal{M}$ [see Eqs. (32)]
and (36)], it will of course not be the case for the true solution $\Psi_0$. Its projection $|\Psi^C\rangle|\Phi_\epsilon\rangle$ onto $\mathcal{M}$ can be determined, in principle exactly, via the diagonalization of an effective Hamiltonian [57, 58]. We can then restore $\Psi_0$ by applying a wave operator $\hat{\Omega}$, i.e.,

$$|\Psi_0\rangle = \hat{\Omega} [ |\Psi^C\rangle |\Phi_\epsilon\rangle ], \quad (44)$$

thus making the embedding formally exact. In practice, one may opt for a perturbation expansion of the (local in this case) correlation potential which can be determined order by order from the density constraint $n_{\Psi_0} = n_{\Phi_{KS}}$, by analogy with Görling-Levy perturbation theory [59–61]. Work is currently in progress in this direction.

In the present work, where the embedding is applied to the uniform 1D Hubbard lattice, the density profile is trivially obtained from the fixed number of electrons $N$ in the lattice. The local correlation potential is an arbitrary constant and, therefore, what is referred to as single-shot embedding in the DMET terminology [48] is sufficient. In other words, as long as the number of electrons $N$ is fixed and we aim at mapping only the correct filling $n = N/L$ onto the embedded impurity, there is no correlation potential to optimize.

### III. SINGLE-SHOT DENSITY MATRIX FUNCTIONAL EMBEDDING

In the following, an exact formulation of Ht-DMFET is derived for a non-interacting lattice (Sec. III A). On that basis, we construct an approximate Ht-DMFET for the interacting lattice (Sec. III B). Finally, connections with DMET and its extensions are made in Sec. III C.

#### A. Exact embedding in the non-interacting case

In the non-interacting (NI) case, which is equivalent to the mean-field approximation in the uniform Hubbard model, the ground-state per-site energy reads [see Eqs. (1) and (7)]

$$e^{\text{NI}}(n) = -4t\gamma_{10}, \quad (45)$$

where we fix the number $N$ of electrons in the lattice and therefore the uniform filling $n = N/L = 2\gamma_{ii}$. Evaluating $\gamma_{10}$ requires solving the NI problem for the full system, which is computationally affordable for a large number of sites (we considered $L = 400$ in this work). For that purpose, we have to minimize the total NI energy

$$E^{\text{NI}} = 2 \sum_{i,j=0}^{L-1} h_{ij} \gamma_{ij}, \quad (46)$$

where

$$h_{ij} \stackrel{0<i<L-1}{=} -t[\delta_{j(i+1)} + \delta_{j(i-1)}], \quad (47)$$

and $h_{0(L-1)} \equiv \pm t$ (the sign depends on the boundary conditions). In practice, we simply need to diagonalize the hopping matrix $\mathbf{h} \equiv \{h_{ij}\}$ and construct the density matrix in the lattice representation from the occupied (orthonormal) eigen-spin-orbitals $|\kappa_\sigma\rangle = \sum_i C_{i\kappa} |c_{i\sigma}\rangle$ as follows,

$$\gamma_{ij} = \sum_{\kappa}^{\text{occ.}} C_{i\kappa} C_{j\kappa}. \quad (48)$$

We propose in this section to reformulate the NI problem into an embedded one. As such, it is in principle useless. However, when it comes to introduce electron correlation, which is of course our ultimate goal, this reformulation will provide a starting point for the embedding. It will also suggest how the latter can be systematically improved, as discussed further in Sec. III B.

Let us rewrite Eq. (46) in the Householder representation,

$$E^{\text{NI}} = 2 \sum_{kl} \tilde{h}_{kl} \tilde{\gamma}_{kl}, \quad (49)$$

where

$$\begin{aligned} \tilde{h}_{kl} &= \sum_{ij} P_{ki} h_{ij} P_{jl} \\ &= h_{kl} - 2 \sum_{i} v_i (v_k h_{il} + v_l h_{ik}) + 4 v_k v_l \sum_{ij} v_i h_{ij} v_j. \end{aligned} \quad (50)$$

Since, in this representation, the (idempotent) density matrix can be split into cluster and environment parts (see Fig. 2), the same applies to the NI energy:

$$E^{\text{NI}} = E_{C}^{\text{NI}} + E_{\mathcal{E}}^{\text{NI}}. \quad (51)$$

The cluster energy reads

$$E_{C}^{\text{NI}} = 2 \sum_{k,l=0}^{1} \tilde{h}_{kl} \tilde{\gamma}_{kl}, \quad (52)$$

or, equivalently,

$$E_{C}^{\text{NI}} = \langle \Phi^{C} | \hat{h}^{C} | \Phi^{C} \rangle, \quad (53)$$

where $\Phi^C$ is the (single determinant) two-electron cluster wave function introduced in Eq. (36). The non-interacting Householder cluster and the Hubbard dimer [62, 63] have formally identical Hamiltonians,

$$\begin{aligned} \hat{h}^C &= \sum_{\sigma} \sum_{i,j=0}^{1} \tilde{h}_{ij} \hat{d}_{i\sigma}^{\dagger} \hat{d}_{j\sigma} \\ &= \sum_{\sigma} [-\tilde{t} (\hat{d}_{0\sigma}^{\dagger} \hat{d}_{1\sigma} + \text{H.c.}) + \tilde{\epsilon}_1 \hat{d}_{1\sigma}^{\dagger} \hat{d}_{1\sigma}], \end{aligned} \quad (54)$$

where, according to Eq. (50),

$$\tilde{t} = -\tilde{h}_{01} = t + 2v_1 \sum_i v_i h_{i0} \quad (55)$$
and

$$ \tilde{\epsilon}_1 = \tilde{h}_{11} = -4v_1 \sum_i v_i \left( h_{i1} - v_1 \sum_j v_j h_{ij} \right). \quad (56) $$

On the other hand, the energy of the environment is an explicit functional of the environment's density matrix $\tilde{\gamma}^{\mathcal{E}} = \{\tilde{\gamma}_{ij}\}_{i>1, j>1}$:

$$ E_{\mathcal{E}}^{\text{NI}} = 2 \sum_{k,l=2}^{L-1} \tilde{h}_{kl} \tilde{\gamma}_{kl} \equiv 2 \text{Tr} [\tilde{h}^{\mathcal{E}} \tilde{\gamma}^{\mathcal{E}}], \quad (57) $$

where Tr denotes the trace. The total ground-state NI energy can be reached variationally, in principle exactly, as follows,

$$ E_0^{\text{NI}} = \min_{\mathbf{v}} \{E_C^{\text{NI}}[\mathbf{v}] + E_{\mathcal{E}}^{\text{NI}}[\mathbf{v}]\}, \quad (58) $$

where

$$ E_C^{\text{NI}}[\mathbf{v}] = \min_{\Phi^C} \langle \Phi^C | \hat{h}^C[\mathbf{v}] | \Phi^C \rangle \quad (59) $$

and

$$ E_{\mathcal{E}}^{\text{NI}}[\mathbf{v}] = 2 \min_{\tilde{\gamma}^{\mathcal{E}}} \text{Tr} [\tilde{h}^{\mathcal{E}}[\mathbf{v}] \tilde{\gamma}^{\mathcal{E}}]. \quad (60) $$

Dependencies in the Householder vector $\mathbf{v}$ have been introduced, for clarity. In order to evaluate the per-site energy [see Eq. (45)], we can switch to the Householder representation,

$$ \begin{aligned} \gamma_{10} &= [\mathbf{P}\tilde{\gamma}\mathbf{P}]_{10} = \sum_{ij} P_{1i}\tilde{\gamma}_{ij}P_{j0} \\ &= \sum_{i} P_{1i}\tilde{\gamma}_{i0} = \sum_{i=0}^{1} P_{1i}\tilde{\gamma}_{i0} \\ &= \sum_{i,j=0}^{1} P_{1i}\tilde{\gamma}_{ij}P_{j0} \\ &= \sum_{ij} P_{1i} \langle \Phi^C | \hat{d}_{i\sigma}^{\dagger} \hat{d}_{j\sigma} | \Phi^C \rangle P_{j0} \\ &= \langle \hat{c}_{1\sigma}^{\dagger} \hat{c}_{0\sigma} \rangle_{\Phi^C}, \end{aligned} \quad (61) $$

thus leading to the final expression

$$ e^{\text{NI}}(n) = -4t \langle \hat{c}_{1\sigma}^{\dagger} \hat{c}_{0\sigma} \rangle_{\Phi^c}. \quad (62) $$

According to Eq. (62), the per-site NI energy can be evaluated directly from the cluster, which is obviously a huge simplification of the full-size problem. The exact NI per-site energy is recovered when the minimizing Householder vector $\mathbf{v}$ in Eq. (58) is employed, thus providing the optimal bath orbital. As readily seen from the latter equation, the Householder vector connects the cluster to its environment energy wise. In the present formalism, the optimal cluster (or, equivalently, the optimal bath) minimizes the sum of the cluster and environment energies. At the NI level, Eqs. (58)–(60) are much more complicated than Eq. (48) implementation wise, especially because the Hamiltonian of the environment should in

principle be diagonalized for each trial Householder vector. However, once the full-size NI problem is solved [with Eq. (48)], Eqs. (58)–(60) can be used for describing two-electron interactions. In the simplest embedding scheme, which is described in the present work, electron correlation is introduced within the cluster while freezing the Householder vector to its NI value. The embedding might then be systematically improved by (i) updating the Householder vector variationally, (ii) describing correlations between the cluster and the environment and, ultimately, (iii) describing correlation within the environment. Such refinements are left for future work.

## B. Correlation on top of the non-interacting Householder embedding

By analogy with DMET [30], we now introduce the exact two-electron repulsion operator within the Householder cluster. For simplicity, we keep on using the Householder vector $\mathbf{v}$ evaluated from the NI density matrix of Eq. (48). First we need to rewrite the on-site repulsion operator in the Householder representation [see Eq. (19)],

$$ U \sum_{i=0}^{L-1} \hat{n}_{i\uparrow} \hat{n}_{i\downarrow} = \sum_{jklm} \tilde{U}_{jklm} \hat{d}_{j\uparrow}^{\dagger} \hat{d}_{k\uparrow} \hat{d}_{l\downarrow}^{\dagger} \hat{d}_{m\downarrow}, \quad (63) $$

where

$$ \tilde{U}_{jklm} = U \sum_{i=0}^{L-1} P_{ij} P_{ik} P_{il} P_{im}. \quad (64) $$

After projecting onto the cluster, we obtain the following expression for the interacting cluster Hamiltonian:

$$ \begin{aligned} \tilde{\mathcal{H}}^C &= \hat{h}^C + \sum_{j,k,l,m=0}^{1} \tilde{U}_{jklm} \hat{d}_{j\uparrow}^{\dagger} \hat{d}_{k\uparrow} \hat{d}_{l\downarrow}^{\dagger} \hat{d}_{m\downarrow} \\ &= \hat{h}^C + U \hat{d}_{0\uparrow}^{\dagger} \hat{d}_{0\uparrow} \hat{d}_{0\downarrow}^{\dagger} \hat{d}_{0\downarrow} + \tilde{U}_1 \hat{d}_{1\uparrow}^{\dagger} \hat{d}_{1\uparrow} \hat{d}_{1\downarrow}^{\dagger} \hat{d}_{1\downarrow}, \end{aligned} \quad (65) $$

where, according to Eq. (64), the effective interaction in the bath is determined from (the NI approximation to) the full lattice as follows,

$$ \begin{aligned} \tilde{U}_1 &= U \sum_{i=1}^{L-1} P_{i1}^4 \\ &= U \left( 1 - 8v_1^2 + 24v_1^4 - 32v_1^6 + 16v_1^4 \sum_{i=1}^{L-1} v_i^4 \right). \end{aligned} \quad (66) $$

In the DMET terminology [35, 48], using the complete embedding Hamiltonian of Eq. (65) is referred to as interacting bath (IB) embedding. In this case, the interaction $\tilde{U}_1$ on the bath site is taken into account. We should stress that, in the present work, like in regular DMET, the bath orbital (which is described by the operators $\hat{d}_{1\sigma}^{\dagger}$ and $\hat{d}_{1\sigma}$) is determined from the density matrix of the full non-interacting lattice, according to Eq. (18). Therefore, the acronym IB should not be
confused with the level of calculation of the bath orbital. In the non-interacting bath (NIB) formulation [35, 48], which is commonly used in DMET and is also tested in Sec. IV, the interaction $\tilde{U}_1$ on the bath site is simply neglected. Note that, in the context of DMET, correlating the bath would require the computation of a correlated many-body wave function for the full system [52] or the introduction of frequency dependencies into the theory [39]. In the present density matrix functional formulation, correlation might be introduced into the bath simply by employing in the Householder transformation [see Eqs. (11)–(13)] a correlated full-system-size density matrix. This can be achieved at a reasonable computational cost [64] with low-level natural orbital functionals (NOFs). Their (sometimes poor [65, 66]) description of strongly correlated energies may then be improved via the embedding, thus avoiding the use of more sophisticated NOFs which can be more difficult to converge. In this context, we may actually proceed with successive Householder transformations since the Householder “impurity+bath” cluster will in principle not be disconnected anymore from its environment, as discussed in Sec. II B. Applying a second Householder transformation to the “bath+buffer+environment” block of the Householder transformed density matrix (see Fig. 2) would generate a second bath orbital. The interacting lattice problem can then be projected onto the enlarged “impurity+two bath orbitals” cluster. Applying further Householder transformations would generate more bath orbitals and, ultimately, make the embedding exact (because equivalent to the full lattice diagonalization problem). Interestingly, in EwDMET, the enlarged number of bath orbitals is determined from the order to which fragment spectral moments should be reproduced. A formal connection between Ht-DMFET and EwDMET might be established at this level. We leave the development of such an extension for future work. Note that, within a hybrid NOF/Householder scheme, we will still be able to choose between IB and NIB formulations. In the rest of the paper, we employ an uncorrelated bath (which is identical to the non-interacting bath for uniform systems).

Let us return to the embedding Hamiltonian in Eq. (65). As the impurity occupation may deviate from the lattice filling $n$ when solving the interacting problem within the cluster, we introduce and adjust a chemical potential $\tilde{\mu}^{\text{imp}}$ on the impurity site, in complete analogy with DMET [48], such that the (two-electron) cluster wave function

$$ \Psi^C = \arg \min_{\Psi} \langle \Psi | \hat{H}^C - \tilde{\mu}^{\text{imp}} \sum_{\sigma} \hat{d}_{0\sigma}^{\dagger} \hat{d}_{0\sigma} | \Psi \rangle \quad (67) $$

reproduces the desired occupation $n$, i.e.

$$ \langle \Psi^C | \sum_{\sigma} \hat{d}_{0\sigma}^{\dagger} \hat{d}_{0\sigma} | \Psi^C \rangle \stackrel{!}{=} n. \quad (68) $$

Once the constraint in Eq. (68) is fulfilled, we obtain an

approximate correlated expression for the per-site energy,

$$ e(n) \approx -4t \langle \hat{c}_{1\sigma}^{\dagger} \hat{c}_{0\sigma} \rangle_{\Psi c} + U \langle \hat{n}_{0\uparrow} \hat{n}_{0\downarrow} \rangle_{\Psi c}, \quad (69) $$

where the density matrix element can be evaluated as follows [see Eq. (20)],

$$ \begin{aligned} \langle \hat{c}_{1\sigma}^{\dagger} \hat{c}_{0\sigma} \rangle_{\Psi c} &= (1 - 2v_1^2) \langle \hat{d}_{1\sigma}^{\dagger} \hat{d}_{0\sigma} \rangle_{\Psi c} \\ &= (1 - 2v_1^2) \tilde{\gamma}_{10}^{\Psi^c}, \end{aligned} \quad (70) $$

and the impurity double occupation simply reads $\langle \hat{n}_{0\uparrow} \hat{n}_{0\downarrow} \rangle_{\Psi c} = \langle \hat{d}_{0\uparrow}^{\dagger} \hat{d}_{0\downarrow} \hat{d}_{0\downarrow}^{\dagger} \hat{d}_{0\uparrow} \rangle_{\Psi c} = d_0^{\Psi^c}$. Note that the Hamiltonian in Eq. (67) describes an asymmetric Hubbard dimer whose two-electron singlet ground-state energy can be expressed analytically in terms of $\tilde{t}$, $U$, $\tilde{U}_1$, $\tilde{\epsilon}_1$, and $\tilde{\mu}^{\text{imp}}$. For that purpose, we simply need to shift the local one-electron potential (on both impurity and bath sites) by the “constant” $(\tilde{\mu}^{\text{imp}} - \tilde{\epsilon}_1)/2$ and to symmetrize the dimer interaction wise [63],

$$ \begin{aligned} \hat{H}^C - \tilde{\mu}^{\text{imp}} \sum_{\sigma} \hat{d}_{0\sigma}^{\dagger} \hat{d}_{0\sigma} &\rightarrow \sum_{\sigma} -\tilde{t} (\hat{d}_{0\sigma}^{\dagger} \hat{d}_{1\sigma} + \text{H.c.}) \\ &+ U_{\text{eff}} \sum_{i=0}^{1} \hat{d}_{i\uparrow}^{\dagger} \hat{d}_{i\downarrow} \hat{d}_{i\downarrow}^{\dagger} \hat{d}_{i\downarrow} + \frac{\Delta v_{\text{eff}}}{2} \sum_{\sigma} (\hat{d}_{1\sigma}^{\dagger} \hat{d}_{1\sigma} - \hat{d}_{0\sigma}^{\dagger} \hat{d}_{0\sigma}), \end{aligned} \quad (71) $$

where the effective interaction and potential read

$$ U_{\text{eff}} = \frac{U + \tilde{U}_1}{2} \quad (72) $$

and

$$ \Delta v_{\text{eff}} = \tilde{\mu}^{\text{imp}} + \tilde{\epsilon}_1 + \frac{\tilde{U}_1 - U}{2}, \quad (73) $$

respectively. Note that the Hamiltonians in the left- and right-hand side of Eq. (71) share the same ground-state wave function $\Psi^C$ [63]. The corresponding ground-state energies are immediately deduced from the analytical expression given in Ref. 62. The chemical potential $\tilde{\mu}^{\text{imp}}$ that fulfills Eq. (68) or, equivalently, the effective potential $\Delta v_{\text{eff}}$, can be determined straightforwardly, for example, by Lieb maximization [63]. Both $\tilde{\gamma}_{10}^{\Psi^c}$ and $d_0^{\Psi^c}$ can then be evaluated analytically via the Hellmann-Feynman theorem (see Appendix A in Ref. 38). In the DMET terminology [48], Eqs. (67)–(69) describe a single-shot embedding.

### C. Summary of the present implementation and connection to DMET

We have described the single-shot embedding of a single impurity site in the particular case of a 1D Hubbard lattice. The full procedure can be summarized as follows. First we solve the non-interacting problem. In the present work we diagonalize the bare hopping matrix [see Eq. (47)] and construct the (idempotent) ground-state density matrix of the full lattice for a given and fixed number $N = nL$ of electrons. The latter density matrix gives immediately access to the bath orbital thanks
to the Householder transformation [see the expression in Eq. (18); see also Eqs. (12) and (13)]. Then we project the original interacting lattice Hamiltonian of Eq. (1) onto the “impurity+bath” many-body subspace, which gives the cluster Hamiltonian expression of Eq. (65). At this level, we can decide to keep the interaction in the bath (IB formulation) or to remove it (NIB formulation). Finally, a chemical potential is introduced and adjusted on the embedded impurity to ensure that its occupation matches the lattice filling n [see Eqs. (67) and (68)]. A correlated per-site energy can then be evaluated from the (interacting) cluster many-body wave function [see Eq. (69)]. For analysis purposes, NIB results obtained with two or three impurities are presented in Sec. IV. In the latter case, we used a block Householder transformation [54] where the column vector $X-Y$ in Eq. (5) is replaced by a $L \times N_{imp}$ matrix which is denoted $V$ in Ref. 54, $N_{imp}$ being the number of impurities. The Householder transformation matrix is then modified as follows,

$$P = I - 2vv^{\dagger} \rightarrow I - 2V [V^{\dagger}V]^{-1} V^{\dagger}. \quad (74)$$

A detailed derivation of multiple-impurity Ht-DMFET will be presented in a separate work.

Let us stress that, at a given level of approximation (we performed a single-shot embedding for simplicity but stronger mapping constraints could of course be employed), the IB formulations of Ht-DMFET and standard DMET are formally equivalent. Indeed, for a non-interacting (or mean-field) lattice, the bath orbitals constructed from the Householder transformation and the Schmidt decomposition are identical, as shown in Sec. IIC. As a result, when projecting the lattice interactions onto the “impurity+bath” cluster, both approaches will lead to the same embedding Hamiltonian [the one in Eq. (65)]. More generally, as long as the bath is uncorrelated, which means that we describe the full lattice with an idempotent (non-interacting or mean-field) density matrix, Ht-DMFET becomes a simpler to implement (but equivalent) version of DMET. Note also that the Householder transformation can in principle be substituted for the Schmidt decomposition when constructing hole and particle bath states in EwDMET [see Eqs. (25) and (26) in Ref. 45]. In this case, we would need to construct density matrices for cationic and anionic systems [see Eq. (10) in Ref. 44].

Finally, if we just see in the Householder transformation a simplification of the Schmidt decomposition in the particular case of non-interacting systems, it becomes clear that, by complete analogy with DMET, Ht-DMFET can be extended to more general Hamiltonians like, for example, quantum chemical ones written in a localized orbital basis. In this case, several fragments (each consisting of multiple localized impurity orbitals) would be employed. More precisely, once a mean-field density matrix has been computed for the full molecule,

it can be Householder transformed with respect to a given fragment, thus providing a set of bath orbitals from which correlated local properties (on the fragment) can be evaluated. Each fragment will have its own Householder transformation but all will apply to the same (full-system) density matrix. Obviously, the block version of the Householder transformation [54] would be employed in this case. Moreover, as already mentioned in Sec. III B, the Householder transformation offers the possibility to correlate the bath through the density matrix, a feature that DMET does not have and from which we may benefit. Work is currently in progress in this direction.

## IV. RESULTS AND DISCUSSION

In this section, Ht-DMFET is applied to a large ($L = 400$) uniform Hubbard ring. The hopping parameter is set to $t = 1$. In order to remove pathological degeneracies from the N-electron full-size non-interacting calculation, periodic/antiperiodic boundary conditions are used when $N/2$ is odd/even. Comparison is made with the Bethe Ansatz (BA) results which are exact in the thermodynamic limit [67]. We note that the single-impurity DMET results presented in Ref. 30 were obtained without the interaction in the bath. The interacting-bath results presented in the following will be of particular interest in this respect. For analysis purposes, multiple-impurity Ht-DMFET calculations have been performed at the NIB level [see Eq. (74)]. In the latter case, the correlated embedded impurity problem has been solved numerically at the density matrix renormalization group (DMRG) level with the Block code [68–72].

Let us first discuss the single-impurity results. We plot in Fig. 6 the deviation in occupation of the (two-electron) Householder cluster’s impurity site from the filling $n$ of the lattice when no additional chemical potential $\tilde{\mu}^{imp}$ is introduced. As expected from Sec. III A, the correct occupation $n$ is recovered in the non-interacting $U = 0$ case (not shown). When $U > 0$, we observe, in the NIB formulation, a systematic depletion of the impurity site as $U$ increases. For given $U$ and $n$ values, adding the ($U$- and $n$-dependent) interaction $\tilde{U}_1$ [see Eq. (66)] on the bath site has the opposite effect. The fact that the interactions on the impurity and the bath do not compensate each other occupation wise can be understood as follows. Let us consider, for example, the half-filled case ($n = 1$) for which $\tilde{\varepsilon}_1 = 0$. When $\tilde{\mu}^{imp} = 0$, in the NIB case, the cluster is equivalent to a symmetric-in-interaction Hubbard dimer with effective interaction $U_{eff} = U/2$ and potential $\Delta v_{eff} = -U/2$ [see Eqs. (72) and (73)], thus leading to $\Delta v_{eff}/U_{eff} \stackrel{\text{NIB}}{=} -1$, which obviously will not give the symmetric density profile expected for $n = 1$. As shown in Fig. 1 of Ref. 73, the occupation of the impurity will decrease with $U_{eff}$ (and therefore $U$) and tend to 0.5 as $U/t \rightarrow +\infty$, which is in agreement with Fig. 6. In the IB
<figure><img src="image_7.png" /><figcaption>FIG. 6: Householder cluster's (single) impurity site occupation plotted as a function of the lattice filling $n$ in various correlation regimes for $\tilde{\mu}^{\text{imp}} = 0$ [see Eq. (67)]. Both interacting (solid lines) and non-interacting (dotted lines) bath cases are shown for analysis purposes. The reference black straight line corresponds to the desired situation (which is ultimately reached by adjusting $\tilde{\mu}^{\text{imp}}$) where the embedded impurity occupation matches the lattice filling.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 8: Ht-DMFET (single) impurity site double occupation $\langle \hat{n}_{0\uparrow} \hat{n}_{0\downarrow} \rangle_{\psi c}$ plotted as a function of the interaction strength at half-filling ($n = 1$). Comparison is made with the exact Bethe Ansatz (BA) result. Non-interacting bath (NIB) results are shown for analysis purposes.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 7: (Single-impurity) Ht-DMFET per-site energy plotted as a function of the interaction strength at half-filling. Results obtained for $\tilde{\mu}^{\text{imp}} = 0$ (dotted lines) and/or without interaction in the bath (blue lines) are shown for analysis purposes. Comparison is made with the exact Bethe Ansatz (BA) result.</figcaption></figure>

case (see Sec. III B), the situation is quite different since

$$ \frac{\Delta v_{\text{eff}}}{U_{\text{eff}}} \stackrel{\text{IB}}{=} \frac{\tilde{U}_1 - U}{\tilde{U}_1 + U} = \frac{\sum_{i=1}^{L-1} P_{i1}^4 - 1}{\sum_{i=1}^{L-1} P_{i1}^4 + 1}, \quad (75) $$

thus leading to $|\Delta v_{\text{eff}}/U_{\text{eff}}| < 1$. In this respect, we are closer to a symmetric dimer, which is an improvement. Nevertheless, since $\Delta v_{\text{eff}} \neq 0$, the correct occupation ($n = 1$) will not be recovered, in general. By computing $\tilde{U}_1$ numerically, we obtained $\tilde{U}_1 \approx U/3$ at half-filling, which gives $\Delta v_{\text{eff}}/U_{\text{eff}} \stackrel{\text{IB}}{=} -1/2$. In this case, as $U$ (and therefore $U_{\text{eff}}$) starts deviating from 0, there will still be an electron depletion of the impurity site [see the top right panel of Fig. 1 in Ref. 73]. Interestingly, for larger $U$ values, the impurity occupation will start increasing and will approach (from below) the correct 1.0 value, as shown in the bottom panels of Fig. 1 in Ref. 73. This is again in perfect agreement with the IB results of Fig. 6. In summary, when no chemical potential is introduced on the impurity, the IB and NIB impurity occupation profiles differ substantially. This difference is driven by the effective ratio $\Delta v_{\text{eff}}/U_{\text{eff}}$. Despite this difference, the IB scheme cannot reproduce the correct occupation and, as a result, introducing a chemical potential $\tilde{\mu}^{\text{imp}}$ remains necessary, like in DMET [48].

Its importance in the calculation of per-site energies is illustrated in Fig. 7 in the particular case of a half-filled lattice. Once a proper $\tilde{\mu}^{\text{imp}}$ value is employed (which will be the case in the rest of the discussion), thus ensuring that the filling and the impurity occupation match, the error becomes substantial in the strongly correlated regime only if the interaction in the bath is neglected. We note that, in the latter case, we reproduce the single-impurity DMET results of Ref. 30, as expected from the analysis in Sec. II C. The agreement with the BA results is almost perfect in all correlation regimes
once the interaction in the bath is restored. This success
may be related to the fact that, like in our approximate
Ht-DMFET scheme, the true (correlated) Householder
cluster contains exactly two electrons in the half-filled
case, as a consequence of the hole-particle symmetry (see
Appendix B). Nevertheless, even though the interaction
in the bath improves on the impurity double occupation,
the error remains substantial when electron correlation
is strong, as shown in Fig. 8. In order to further reduce
the error, more impurities should be introduced into
the cluster [30]. Therefore, the success of the present
(single-impurity) Ht-DMFET at half-filling relies also
on error cancellations in the evaluation of the (total)
per-site energy.

Away from half-filling, the performance of (single-
impurity) Ht-DMFET deteriorates as U/t increases, as
shown in Fig. 9, probably because fluctuations in the
number of electrons within our (“single impurity+single
bath”) cluster are not allowed in our approximate
embedding. As discussed in Sec. II B, away from
half-filling, the cluster becomes an open subsystem
as soon as U/t deviates from zero. Surprizingly, in
this density regime, per-site energies are in better
agreement with the BA values when the interaction
in the bath is neglected. Again, in the latter case, we
recover the single-impurity DMET results of Ref. 30.
As expected [30, 35] and shown in the bottom panel of
Fig. 9, the results dramatically improve when a larger
fragment (consisting of two or three impurities) is em-
bedded, even at the simplest NIB level of approximation.

Finally, we investigate in Fig. 10 the density-driven
Mott-Hubbard transition via the evaluation of the
density-functional μ(n) = ∂e(n)/∂n chemical potential
from the Ht-DMFET energy expression of Eq. (69). As
expected from Ref. 30, at the single-impurity level, there
is no gap opening when the interaction in the bath is ne-
glected. Restoring the interaction in the bath has actu-
ally no impact on the transition. In the light of Sec. IIB,
we can reasonably assume that Ht-DMFET fails in this
case because it relies on a closed two-electron “single im-
purity+single bath” cluster. Already at the NIB level of
approximation, the embedding of a larger fragment (con-
sisting of two or three impurities) substantially improves
the results. Nevertheless, even in this case, the gap re-
mains closed, which is in perfect agreement with the DET
results of Ref. 35. As we perform single-shot embeddings
(where we only require the embedded impurity to repro-
duce the correct filling n), we expect from Ref. 35 the
transition to be better described at the multiple-impurity
level when the interactions in the bath are taken into
account. It would also be interesting to see how Ht-
DMFET performs when a correlated (through the den-
sity matrix) bath is employed. This is left for future
work.

<figure><img src="image_10.png" /><figcaption>FIG. 9: Ht-DMFET per-site energies plotted as a function of the lattice filling $n$ for various correlation regimes. Results obtained with a single impurity are shown as (colored) solid lines. The blue color corresponds to the non-interacting bath (NIB) case. In the strongly correlated $U/t = 8$ regime (bottom panel), NIB results obtained with two ($N_{imp} = 2$) and three ($N_{imp} = 3$) impurities are also shown (as points), for analysis purposes (see Sec. III C for further details). Comparison is made with the exact Bethe Ansatz (BA) results (black solid lines). In the weakly $U/t = 1$ correlated case (top panel), exact and approximate results are almost indistinguishable.</figcaption></figure>

V. CONCLUSIONS AND PERSPECTIVES

Similar in spirit to DMET, a (static and zero-
temperature) single-impurity Householder transformed
density matrix functional embedding theory (Ht-
DMFET) has been derived. The theory has been applied
to the 1D Hubbard model. In the non-interacting case,
the formal reduction of the full lattice to a two-electron
dimer is exact. Thanks to the Householder transforma-
tion, the bath site can be determined (analytically) from
the density matrix of the (full) lattice. Alternatively,
one may determine, in principle exactly, the House-
holder vector **v** (which defines the transformation) by
minimizing the sum of the **v**-dependent Householder
cluster and environment energies. While the two-site
“impurity+bath” cluster problem is trivially solved, the
ground-state energy of the cluster’s environment must
<figure><img src="image_11.png" /><figcaption>FIG. 10: Lattice filling plotted, *via* the relation $\mu \equiv \mu(n) = \partial e(n)/\partial n$, as a function of the (lattice) chemical potential $\mu$ at the Ht-DMFET level of calculation for various correlation regimes. (Single-impurity) non-interacting bath (NIB) results are shown as solid blue lines. In the strongly correlated $U/t = 8$ case, NIB results obtained with $N_{imp} = 2$ and $N_{imp} = 3$ impurities are also shown (as points), for analysis purposes. Comparison is made with the exact Bethe Ansatz (BA) results.</figcaption></figure>

be evaluated for each trial vector $v$. Even though such a strategy is uselessly complicated in practice, because the non-interacting full-size problem can be solved directly, it is enlightening in many ways. First, it clearly shows that the optimal cluster cannot be determined without learning from its environment (and therefore from the full lattice). The two subsystems “communicate” through the Householder vector. Secondly, the resulting variational character of the bath might be exploited when electron correlation is introduced into the theory. One could also extract an approximate Householder vector from many-body perturbation theory. This is left for future work. Finally, the single occupied spin-orbital that overlaps with the impurity is, from the beginning, automatically separated from the core embedding orbitals which, in the terminology of Ht-DMFET, belong to the Householder environment. As a result, for non-interacting systems, Ht-DMFET is equivalent to DMET. It is however simpler since the embedding is constructed analytically from the lattice representation of the density matrix.

Starting from the exact non-interacting and closed two-electron Householder cluster, correlation can be introduced straightforwardly by Householder transforming (and projecting onto the cluster) the on-site two-electron repulsion operator, which is defined in the lattice representation. At half-filling ($n = 1$), the

resulting (approximate) Ht-DMFET per-site energies are in almost perfect agreement with the Bethe Ansatz (BA) results in all correlation regimes provided that (i) a chemical potential is introduced on the impurity site, like in DMET, thus ensuring that the correct filling is reproduced, and (ii) the interaction in the bath is taken into account. The good performance of Ht-DMFET in this case can be partly related to the fact that, at half-filling, the true (correlated) Householder cluster contains exactly two electrons, as a consequence of the hole-particle symmetry. Away from half-filling, the deviation from the BA results becomes substantial in the strongly correlated regime. The results dramatically improve when a larger fragment (consisting of two or three impurities) is embedded. In the latter case, a block Householder transformation was employed [54]. The failure of the single-impurity embedding away from half-filling may originate from the fact that, unlike in our approximate Ht-DMFET, the number of electrons in the exact (correlated) Householder cluster can be fractional. This is probably the reason why, like single-impurity DMET with a non-interacting bath [30], (single-impurity) Ht-DMFET fails in describing the density-driven Mott-Hubbard transition, whether the bath is interacting or not.

In the light of recent advances in DMET and related approaches, several extensions of the present work can already be foreseen. A multiple-impurity implementation of Ht-DMFET that is applicable to more general Hamiltonians (like the quantum chemical ones) is the most urgent. Preserving a Householder cluster that is disconnected from its environment is convenient in practice but only relevant in the non-interacting case. Applying the Householder transformation to the Kohn-Sham density matrix would make the approach formally exact. Introducing site (or orbital) occupation mapping constraints in a self-consistency loop, like in DET [35] or SDE [50], would make complete sense in this context. Work is currently in progress in these directions.

## ACKNOWLEDGMENTS

The authors would like to thank Laurent Mazouin, Bruno Senjean, and Mauricio Rodriguez-Mayorga for fruitful discussions. They also thank LabEx CSC (ANR-10-LABX-0026-CSC) and ANR (ANR-19-CE29-0002 DESCARTES project) for funding.
# Appendix A: Derivation of the Householder transformation

We start from the decomposition

$$Y = (I - v\bar{v}^\dagger)Y + v\bar{v}^\dagger Y. \quad (A.1)$$

Using real algebra and the constraint $|X| = |Y|$ gives, according to Eq. (5),

$$ 
\begin{aligned}
v^\dagger Y &= \frac{X^\dagger Y - |Y|^2}{|X - Y|} \\
&= \frac{Y^\dagger X - |X|^2}{|X - Y|} \quad (A.2) \\
&= \frac{(Y^\dagger - X^\dagger) X}{|X - Y|} \\
&= -v^\dagger X,
\end{aligned}
$$

and

$$(I - v\bar{v}^\dagger) v = v - (v^\dagger v) v = 0, \quad (A.3)$$

or, equivalently,

$$(I - v\bar{v}^\dagger) X = (I - v\bar{v}^\dagger) Y. \quad (A.4)$$

Combining Eqs. (A.1), (A.2), and (A.4) leads to

$$Y = (I - v\bar{v}^\dagger) X - v\bar{v}^\dagger X \equiv PX, \quad (A.5)$$

where $P = I - 2v\bar{v}^\dagger$.

# Appendix B: Exact properties of the Householder cluster for a half-filled 1D interacting lattice

Let us evaluate the matrix element

$$[\tilde{\gamma}^2]_{j0} = [P\gamma^2 P]_{j0} = \sum_{klm} P_{jk}\gamma_{kl}\gamma_{lm}P_{m0}, \quad (B.1)$$

or, equivalently,

$$ 
\begin{aligned}
[\tilde{\gamma}^2]_{j0} &= \sum_{kl} P_{jk} \gamma_{kl} \gamma_{l0} \\
&\stackrel{j \ge 0}{=} \sum_{l} P_{j1} \gamma_{1l} \gamma_{l0} + \sum_{k>1} \sum_{l} P_{jk} \gamma_{kl} \gamma_{l0} \\
&\stackrel{j \ge 0}{=} P_{j1} \gamma_{10} (\gamma_{00} + \gamma_{11}) + P_{j1} \sum_{l>1} \gamma_{1l} \gamma_{l0} \quad (B.2) \\
&\quad + \gamma_{00} \sum_{k>1} P_{jk} \gamma_{k0} + \gamma_{10} \sum_{k>1} P_{jk} \gamma_{k1} \\
&\quad + \sum_{k>1} \sum_{l>1} P_{jk} \gamma_{kl} \gamma_{l0}.
\end{aligned}
$$

Since $\gamma_{ii} = n_\sigma$ is the uniform spin occupation in the system and $n = 2n_\sigma$ is the lattice filling, it comes

$$ 
\begin{aligned}
[\tilde{\gamma}^2]_{j0} &\stackrel{j \ge 0}{=} n P_{j1} \gamma_{10} + n \sum_{k>1} P_{jk} \gamma_{k0} \\
&\quad + P_{j1} \sum_{l>1} \gamma_{1l} \gamma_{l0} + \gamma_{10} \sum_{k>1} \delta_{jk} \gamma_{k1} \\
&\quad - \frac{2v_j \gamma_{10}}{\sqrt{2\xi(\xi - \gamma_{10})}} \sum_{k>1} \gamma_{k0} \gamma_{k1} \quad (B.3) \\
&\quad + \sum_{k>1} \sum_{1<l \neq k} \delta_{jk} \gamma_{kl} \gamma_{l0} \\
&\quad - \frac{2v_j}{\sqrt{2\xi(\xi - \gamma_{10})}} \sum_{k>1} \sum_{1<l \neq k} \gamma_{k0} \gamma_{kl} \gamma_{l0}.
\end{aligned}
$$

Moreover,

$$ 
\begin{aligned}
P_{j1}\gamma_{10} + \sum_{k>1} P_{jk}\gamma_{k0} &= \delta_{j1}\gamma_{10} - 2v_j \frac{(\gamma_{10} - \xi)\gamma_{10}}{\sqrt{2\xi(\xi - \gamma_{10})}} \\
&\quad - 2v_j \sum_{k>1} \frac{\gamma_{k0}^2}{\sqrt{2\xi(\xi - \gamma_{10})}} \quad (B.4) \\
&\quad + \sum_{k>1} \delta_{jk}\gamma_{k0},
\end{aligned}
$$

where, according to Eq. (9),

$$ 
\begin{aligned}
& \frac{(\gamma_{10} - \xi)\gamma_{10}}{\sqrt{2\xi(\xi - \gamma_{10})}} + \sum_{k>1} \frac{\gamma_{k0}^2}{\sqrt{2\xi(\xi - \gamma_{10})}} \\
&= \frac{1}{\sqrt{2\xi(\xi - \gamma_{10})}} (\xi^2 - \xi\gamma_{10}) \quad (B.5) \\
&= \frac{1}{2\sqrt{2\xi(\xi - \gamma_{10})}},
\end{aligned}
$$

thus leading to

$$ 
\begin{aligned}
P_{j1}\gamma_{10} + \sum_{k>1} P_{jk}\gamma_{k0} &= \delta_{j1}\gamma_{10} - v_j \sqrt{2\xi(\xi - \gamma_{10})} \\
&\quad + \sum_{k>1} \delta_{jk}\gamma_{k0},
\end{aligned}
\quad (B.6)$$

and, consequently, to the final expression

$$ 
\begin{aligned}
[\tilde{\gamma}^2]_{j0} &\stackrel{j \ge 0}{=} n \left( \gamma_{10} \delta_{j1} - \sqrt{2\xi(\xi - \gamma_{10})} v_j \right) \\
&\quad + \sum_{k>1} \delta_{jk} \left( n \gamma_{k0} + \gamma_{10} \gamma_{k1} + \sum_{1<l \neq k} \gamma_{kl} \gamma_{l0} \right) \\
&\quad + \left[ P_{j1} - \frac{2v_j \gamma_{10}}{\sqrt{2\xi(\xi - \gamma_{10})}} \right] \sum_{k>1} \gamma_{k0} \gamma_{k1} \\
&\quad - \frac{4v_j}{\sqrt{2\xi(\xi - \gamma_{10})}} \sum_{1<l<k} \gamma_{k0} \gamma_{kl} \gamma_{l0}.
\end{aligned}
\quad (B.7)$$

As shown in the following, the last two terms on the right-hand side of Eq. (B.7) vanish at half-filling (i.e., when $n = 1$). To prove this, let us consider the modified lattice Hamiltonian

$$ \hat{H} \rightarrow \hat{H}_{kl}(\eta) = \hat{H} + \eta \sum_{\sigma} (\hat{c}_{k\sigma}^{\dagger} \hat{c}_{l\sigma} + \hat{c}_{l\sigma}^{\dagger} \hat{c}_{k\sigma}), \quad (B.8) $$
where $0 \le k < l < L$ are fixed site labels. According to the Hellmann-Feynman theorem, the exact N-electron density matrix element for the Hubbard Hamiltonian $\hat{H}$ can be determined as follows:

$$ \gamma_{kl}^{N} := \langle \hat{c}_{k\sigma}^{\dagger} \hat{c}_{l\sigma} \rangle_{\hat{H}}^{N} = \frac{1}{4} \left. \frac{dE_{kl}^{N}(\eta)}{d\eta} \right|_{\eta=0}, \quad (B.9) $$

where $E_{kl}^N(\eta)$ is the N-electron energy for the Hamiltonian $\hat{H}_{kl}(\eta)$. If we now apply the following hole-particle transformation [21],

$$ \hat{c}_{i\sigma} \rightarrow \hat{b}_{i\sigma} = (-1)^i \hat{c}_{i\sigma}^{\dagger}, \quad (B.10) $$

the latter Hamiltonian can be rewritten as follows,

$$ \begin{aligned} \hat{H}_{kl}(\eta) = & -t \sum_{\sigma} \sum_{i=0}^{L-1} (\hat{b}_{i\sigma}^{\dagger} \hat{b}_{(i+1)\sigma} + \text{H.c.}) \\ & + U \sum_{i=0}^{L-1} \hat{b}_{i\uparrow}^{\dagger} \hat{b}_{i\uparrow} \hat{b}_{i\downarrow}^{\dagger} \hat{b}_{i\downarrow} \\ & + U \left( L - \sum_{\sigma} \sum_{i=0}^{L-1} \hat{b}_{i\sigma}^{\dagger} \hat{b}_{i\sigma} \right) \\ & + (-1)^{l-k-1} \eta \sum_{\sigma} (\hat{b}_{k\sigma}^{\dagger} \hat{b}_{l\sigma} + \hat{b}_{l\sigma}^{\dagger} \hat{b}_{k\sigma}). \end{aligned} \quad (B.11) $$

Note that $\hat{b}_{L\sigma} = (-1)^L \hat{c}_{L\sigma}^\dagger = \pm(-1)^L \hat{c}_{0\sigma}^\dagger = \pm(-1)^L \hat{b}_{0\sigma}$, which means that the same (periodic or antiperiodic) boundary conditions can be used in the particle [Eq. (B.8)] or hole [Eq. (B.11)] Hamiltonians for a total even number L of sites. In this case, we can conclude from Eq. (B.11) that

$$ E_{kl}^{2L-N}(\eta) = E_{kl}^{N}((-1)^{l-k-1}\eta) + U(L-N), \quad (B.12) $$

thus leading to [see Eq. (B.9)]

$$ \gamma_{kl}^{2L-N} = (-1)^{l-k-1} \gamma_{kl}^{N}. \quad (B.13) $$

Therefore, in the half-filled $n = N/L = 1$ case, we have

$$ \gamma_{kl}^{N=L} \times (1 + (-1)^{l-k}) = 0, \quad (B.14) $$

that we simply denote

$$ \gamma_{kl} \times (1 + (-1)^{l-k}) \stackrel{n=1}{=} 0. \quad (B.15) $$

In conclusion, at half-filling, a density matrix element equals zero if the two indices differ by an even number:

$$ \gamma_{kl} \stackrel{n=1}{=} 0 \quad \text{if} \quad l-k = 2p > 0. \quad (B.16) $$

If we now return to Eq. (B.7), we immediately see that, for $k > 1$,

$$ \gamma_{k0} \gamma_{k1} \stackrel{n=1}{=} 0 \quad (B.17) $$

and

$$ \gamma_{k0} \gamma_{kl} \gamma_{l0} \stackrel{n=1}{=} 0, \quad (B.18) $$

so that

$$ \begin{aligned} [\tilde{\gamma}^2]_{j0} \stackrel{j>0, n=1}{=} & (\gamma_{10} \delta_{j1} - \sqrt{2\xi(\xi - \gamma_{10})} v_j) \\ & + \sum_{k>1} \delta_{jk} \left( \gamma_{k0} + \gamma_{10}\gamma_{k1} + \sum_{1<l \neq k} \gamma_{kl}\gamma_{l0} \right). \end{aligned} \quad (B.19) $$

We can now extract interesting information from our simplified Eq. (B.19). For example, for $j = 1$, we obtain, according to Eq. (12),

$$ [\tilde{\gamma}^2]_{10} \stackrel{n=1}{=} \gamma_{10} - \sqrt{2\xi(\xi - \gamma_{10})} v_1 = \xi, \quad (B.20) $$

which, according to Eq. (26), leads to $\mathcal{N}^C = 2$. Moreover,

$$ \begin{aligned} [\tilde{\gamma}^2]_{(2p+1)0} \stackrel{p>0, n=1}{=} & -\gamma_{(2p+1)0} + \gamma_{(2p+1)0} \\ & + \gamma_{10}\gamma_{(2p+1)1} + \sum_{1<l \neq 2p+1} \gamma_{(2p+1)l}\gamma_{l0}, \end{aligned} \quad (B.21) $$

thus leading to [see Eqs. (24) and (B.16)]

$$ \tilde{\gamma}_{(2p+1)1} \stackrel{p>0, n=1}{=} 0. \quad (B.22) $$

On the other hand, the Householder transformed density matrix elements in column “1” are a priori nonzero for even row indices,

$$ \begin{aligned} \xi \tilde{\gamma}_{(2p)1} \stackrel{p>0}{=} [\tilde{\gamma}^2]_{(2p)0} \stackrel{p>0, n=1}{=} & \gamma_{10} \gamma_{(2p)1} + \sum_{1<l \neq 2p} \gamma_{(2p)l} \gamma_{l0} \\ & \neq 0, \end{aligned} \quad (B.23) $$

which means that the cluster is not disconnected from its environment, even though it contains exactly two electrons.

## Appendix C: “Direct” correlation effects on the Householder transformation

Let us consider the idempotent non-interacting density matrix $\gamma$ in the lattice representation, $P = P[\gamma]$ [see Eqs. (11)–(13)] the corresponding Householder transformation matrix, and $\tilde{\gamma} = P\gamma P$. We want to see how the Householder vector is directly affected (i.e., through modifications of the density matrix in the lattice representation) by changes

$$ \tilde{\gamma} \rightarrow \tilde{\gamma} + \Delta \quad (C.1) $$

that may occur in the Householder representation of the density matrix, for example, when electron correlation is introduced. For that purpose, we need to return to the lattice representation,

$$ \begin{aligned} \gamma \rightarrow \bar{\gamma} &= P(\tilde{\gamma} + \Delta)P \\ &= \gamma + P\Delta P, \end{aligned} \quad (C.2) $$

thus leading to

$$ \gamma_{i0} \rightarrow \bar{\gamma}_{i0} = \gamma_{i0} + \sum_{j} P_{ij} \Delta_{j0} \quad (C.3) $$

$$ = \gamma_{i0} + \Delta_{i0} - 2v_i \sum_{j>0} v_j \Delta_{j0}, \quad (C.4) $$
or, in a more compact form,

$$X \to \bar{X} = X + P\Delta_0, \qquad (C.5)$$

where $\Delta_0^\dagger = [\Delta_{00}, \Delta_{10}, \dots, \Delta_{j0}, \dots]$. As readily seen from Eq. (C.4), correlating within the Householder environment (which would imply $\Delta_{j0} = 0, \forall j$) has not impact on $X$ and, therefore, no impact on the Householder transformation. However, changes occur on all rows of $X$ when the density matrix is modified within the cluster, i.e., when $\Delta_{j0} = \sum_{i=0}^1 \delta_{ij}\Delta_{i0}$. As we will see, even though $X$ and $Y = PX$ change in this case, $X-Y$ will not change direction.

To prove the above statement, we need to evaluate the change $Y^\dagger \to \bar{Y}^\dagger = [\bar{\gamma}_{00}, \bar{\xi}, 0, \dots, 0]$ in $Y$. Since, by construction,

$$|\bar{Y}|^2 = |\bar{X}|^2 = |P\bar{X}|^2 = |Y + \Delta_0|^2, \qquad (C.6)$$

it comes from Eq. (C.4) that

$$\begin{aligned} \bar{\xi}^2 &= |\bar{Y}|^2 - \bar{\gamma}_{00}^2 \\ &= \sum_i (Y_i + \Delta_{i0})^2 - (\gamma_{00} + \Delta_{00})^2 \\ &= (\xi + \Delta_{10})^2 + \sum_{i>1} \Delta_{i0}^2, \end{aligned} \qquad (C.7)$$

or, equivalently,

$$\bar{\xi} = (\xi + \Delta_{10}) \sqrt{1 + \frac{\sum_{i>1} \Delta_{i0}^2}{(\xi + \Delta_{10})^2}}, \qquad (C.8)$$

thus leading to

$$\bar{\xi} - \xi = \xi \left( \sqrt{1 + \frac{\sum_{i>1} \Delta_{i0}^2}{(\xi + \Delta_{10})^2} - 1} + \Delta_{10} \sqrt{1 + \frac{\sum_{i>1} \Delta_{i0}^2}{(\xi + \Delta_{10})^2}} \right). \qquad (C.9)$$

Note that the sign of $\bar{\xi}$ has been chosen such that $\bar{\xi} = \xi$ when $\Delta_0 = 0$, as expected. Since

$$\begin{aligned} \bar{X} - \bar{Y} &= \bar{X} - Y - (\bar{Y} - Y) \\ &= (X - Y) + P\Delta_0 - \Delta_{00}e_0 \qquad (C.10) \\ &\quad + (\xi - \bar{\xi})e_1, \end{aligned}$$

where $e_i$ are the unit basis vectors ($e_0^\dagger = [1, 0, 0, \dots, 0]$, $e_1^\dagger = [0, 1, 0, \dots, 0]$), we can now determine the direction of the updated Householder vector:

$$\bar{v} \sim \bar{X} - \bar{Y} = (|X - Y| - 2v^\dagger \Delta_0) v + \sum_{i>0} [\Delta_{i0} + \delta_{i1} (\xi - \bar{\xi})] e_i. \qquad (C.11)$$

As readily seen from Eqs. (C.9) and (C.11), when correlation is treated within the Householder cluster only (i.e., when $\Delta_{i0} \stackrel{i \to 1}{=} 0$), then $\bar{\xi} - \xi = \Delta_{10}$ and

$$\bar{v} \sim (|X - Y| - 2v_1 \Delta_{10}) v \sim v, \qquad (C.12)$$

which, after normalization, leads to $\bar{v} = v$ and $\bar{P} = P$.

<footnote n="1">1 A. Wasserman and M. Pavanello, Int. J. Quantum Chem. **120**, e26495 (2020).</footnote>

<footnote n="2">2 A. Georges and G. Kotliar, Phys. Rev. B **45**, 6479 (1992).</footnote>

<footnote n="3">3 A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. **68**, 13 (1996).</footnote>

<footnote n="4">4 G. Kotliar and D. Vollhardt, Phys. Today **57**, 53 (2004).</footnote>

<footnote n="5">5 K. Held, Adv. Phys. **56**, 829 (2007).</footnote>

<footnote n="6">6 D. Zgid and G. K.-L. Chan, J. Chem. Phys. **134**, 094115 (2011).</footnote>

<footnote n="7">7 M. Potthoff, Phys. Rev. B **64**, 165114 (2001).</footnote>

<footnote n="8">8 G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. **78**, 865 (2006).</footnote>

<footnote n="9">9 P. Sun and G. Kotliar, Phys. Rev. B **66**, 085120 (2002).</footnote>

<footnote n="10">10 S. Biermann, F. Aryasetiawan, and A. Georges, Phys. Rev. Lett. **90**, 086402 (2003).</footnote>

<footnote n="11">11 K. Haule, Phys. Rev. Lett. **115**, 196403 (2015).</footnote>

<footnote n="12">12 L. Boehnke, F. Nilsson, F. Aryasetiawan, and P. Werner, Phys. Rev. B **94**, 201106 (2016).</footnote>

<footnote n="13">13 P. Werner and M. Casula, J. Phys. Condens. Matter **28**, 383001 (2016).</footnote>

<footnote n="14">14 F. Nilsson, L. Boehnke, P. Werner, and F. Aryasetiawan, Phys. Rev. Materials **1**, 043803 (2017).</footnote>

<footnote n="15">15 A. A. Kananenka, E. Gull, and D. Zgid, Phys. Rev. B **91**, 121111 (2015).</footnote>

<footnote n="16">16 T. N. Lan, A. A. Kananenka, and D. Zgid, J. Chem. Phys. **143**, 241102 (2015).</footnote>

<footnote n="17">17 D. Zgid and E. Gull, New J. Phys. **19**, 023047 (2017).</footnote>

<footnote n="18">18 T. N. Lan and D. Zgid, J. Phys. Chem. Lett. **8**, 2200 (2017).</footnote>

<footnote n="19">19 M. Dvorak and P. Rinke, Phys. Rev. B **99**, 115134 (2019).</footnote>

<footnote n="20">20 E. Fromager, Mol. Phys. **113**, 419 (2015).</footnote>

<footnote n="21">21 B. Senjean, N. Nakatani, M. Tsuchiizu, and E. Fromager, Phys. Rev. B **97**, 235105 (2018).</footnote>

<footnote n="22">22 B. Senjean, N. Nakatani, M. Tsuchiizu, and E. Fromager, Theor. Chem. Acc. **137**, 169 (2018).</footnote>

<footnote n="23">23 R. Requist and E. K. U. Gross, Phys. Rev. B **99**, 125114 (2019).</footnote>

<footnote n="24">24 J. Toulouse, F. Colonna, and A. Savin, Phys. Rev. A **70**, 062505 (2004).</footnote>

<footnote n="25">25 A. Ferté, E. Giner, and J. Toulouse, J. Chem. Phys. **150**, 084103 (2019).</footnote>

<footnote n="26">26 L. Gagliardi, D. G. Truhlar, G. Li Manni, R. K. Carlson, C. E. Hoyer, and J. L. Bao, Acc. Chem. Res. **50**, 66 (2016).</footnote>

<footnote n="27">27 S. Ghosh, P. Verma, C. J. Cramer, L. Gagliardi, and D. G. Truhlar, Chem. Rev. **118**, 7249 (2018).</footnote>

<footnote n="28">28 T. A. Wesolowski, S. Shedge, and X. Zhou, Chem. Rev. **115**, 5891 (2015).</footnote>

<footnote n="29">29 M. A. Mosquera, L. O. Jones, C. H. Borca, M. A. Ratner, and G. C. Schatz, J. Phys. Chem. A **123**, 4785 (2019).</footnote>
30 G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012).

31 G. Knizia and G. K.-L. Chan, J. Chem. Theory Comput. 9, 1428 (2013).

32 Q. Sun and G. K.-L. Chan, Acc. Chem. Res. 49, 2705 (2016).

33 X. Wu, Z.-H. Cui, Y. Tong, M. Lindsey, G. K.-L. Chan, and L. Lin, J. Chem. Phys. 151, 064108 (2019).

34 Z.-H. Cui, T. Zhu, and G. K.-L. Chan, J. Chem. Theory Comput. 16, 119 (2020).

35 I. W. Bulik, G. E. Scuseria, and J. Dukelsky, Phys. Rev. B 89, 035140 (2014).

36 T. Tsuchimochi, M. Welborn, and T. Van Voorhis, J. Chem. Phys. 143, 024107 (2015).

37 M. Welborn, T. Tsuchimochi, and T. Van Voorhis, J. Chem. Phys. 145, 074102 (2016).

38 B. Senjean, Phys. Rev. B 100, 035136 (2019).

39 P. V. Sriluckshmy, M. Nusspickel, E. Fertitta, and G. H. Booth, Phys. Rev. B 103, 085131 (2021).

40 G. H. Booth and G. K.-L. Chan, Phys. Rev. B 91, 155107 (2015).

41 J. S. Kretchmer and G. K.-L. Chan, J. Chem. Phys. 148, 054108 (2018).

42 T. Ayral, T.-H. Lee, and G. Kotliar, Phys. Rev. B 96, 235139 (2017).

43 T.-H. Lee, T. Ayral, Y.-X. Yao, N. Lanata, and G. Kotliar, Phys. Rev. B 99, 115129 (2019).

44 E. Fertitta and G. H. Booth, Phys. Rev. B 98, 235132 (2018).

45 E. Fertitta and G. H. Booth, J. Chem. Phys. 151, 014115 (2019).

46 L. Lacombe and N. T. Maitra, Phys. Rev. Lett. 124, 206401 (2020).

47 R. Requist and E. K. U. Gross, "Fock space embedding theory for strongly correlated topological phases," (2019), arXiv:1909.07933 [cond-mat.str-el].

48 S. Wouters, C. A. Jiménez-Hoyos, Q. Sun, and G. K.-L. Chan, J. Chem. Theory Comput. 12, 2706 (2016).

49 A. H. Mühlbach and M. Reiher, J. Chem. Phys. 149, 184104 (2018).

50 U. Mordovina, T. E. Reinhard, I. Theophilou, H. Appel, and A. Rubio, J. Chem. Theory Comput. 15, 5209 (2019).

51 Z. Fan and Q.-L. Jie, Phys. Rev. B 91, 195118 (2015).

52 M. R. Hermes and L. Gagliardi, J. Chem. Theory Comput. 15, 972 (2019).

53 A. S. Householder, J. ACM 5, 339 (1958).

54 F. Rotella and I. Zambettakis, Appl. Math. Lett. 12, 29 (1999).

55 X. Wu, M. Lindsey, T. Zhou, Y. Tong, and L. Lin, Phys. Rev. B 102, 085123 (2020).

56 W. Kohn and L. Sham, Phys. Rev. 140, A1133 (1965).

57 I. Lindgren and J. Morrison, *Atomic many-body theory* (Springer-Verlag Berlin Heidelberg, 1986).

58 C. Brouder, G. H. Duchamp, F. Patras, and G. Z. Tóth, Int. J. Quantum Chem. 112, 2256 (2012).

59 A. Görling and M. Levy, Phys. Rev. A 50, 196 (1994).

60 A. Görling and M. Levy, Int. J. Quantum Chem. 56, 93 (1995).

61 S. Ivanov and M. Levy, J. Chem. Phys. 116, 6924 (2002), https://doi.org/10.1063/1.1453952.

62 D. J. Carrascal, J. Ferrer, J. C. Smith, and K. Burke, J. Phys. Condens. Matter 27, 393001 (2015).

63 B. Senjean, M. Tsuchiizu, V. Robert, and E. Fromager, Mol. Phys. 115, 48 (2017).

64 K. J. H. Giesbertz, Phys. Chem. Chem. Phys. 18, 21024 (2016).

65 I. Mitxelena, M. Piris, and M. Rodriguez-Mayorga, J. Phys. Condens. Matter 29, 425602 (2017).

66 I. Mitxelena, M. Piris, and M. Rodríguez-Mayorga, J. Phys. Condens. Matter 30, 089501 (2018).

67 E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).

68 G. K.-L. Chan and M. Head-Gordon, J. Chem. Phys. 116, 4462 (2002).

69 G. K.-L. Chan, J. Chem. Phys. 120, 3172 (2004).

70 D. Ghosh, J. Hachmann, T. Yanai, and G. K.-L. Chan, J. Chem. Phys. 128, 144117 (2008).

71 S. Sharma and G. K.-L. Chan, J. Chem. Phys. 136, 124121 (2012).

72 S. Guo, M. A. Watson, W. Hu, Q. Sun, and G. K.-L. Chan, J. Chem. Theory Comput. 12, 1583 (2016).

73 K. Deur, L. Mazouin, and E. Fromager, Phys. Rev. B 95, 035120 (2017).