# A versatile unitary transformation framework for an optimal bath
construction in density-matrix based quantum embedding approaches

Quentin Marécat<EMAIL_ADDRESS>ICGM, Université de
Montpellier, CNRS, ENSCM, 34000 Montpellier (France)  Matthieu Saubanère
ICGM, Université de Montpellier, CNRS, ENSCM, 34000 Montpellier (France)

###### Abstract

The performance of embedding methods is directly tied to the quality of the
bath orbitals construction. In this paper, we develop a versatile framework,
enabling the investigation of the optimal construction of the orbitals of the
bath. As of today, in state-of-the-art embedding methods, the orbitals of the
bath are constructed by performing a Singular Value Decomposition (SVD) on the
impurity-environment part of the 1RDM, as originally presented in Density
Matrix Embedding Theory (DMET). Recently, the equivalence between the SVD
protocol and the use of unitary transformation, the so-called Block-
Householder transformation, has been established. We present a generalization
of the Block-Householder transformation by introducing additional flexible
parameters. The additional parameters are optimized such that the bath-
orbitals fulfill physically motivated constrains. The efficiency of the
approach is discussed and examplified in the context of the half-filled
Hubbard model in one-dimension.

##  I Introduction

Quantum chemistry is built on the foundation of solving the Schrödinger
equation, an endeavor that quickly becomes computationally prohibitive for
systems of practical interest. The challenge comes from the exponentially
scaling cost with system size. In that context, quantum embedding methods Sun
and Chan (2016); Wasserman and Pavanello (2020) have emerged as a powerful
tool, particularly for the study of strongly correlated systems Kotliar  _et
al._ (2001); Zheng (1993); Ma  _et al._ (2021), where traditional methods
often fall short Cohen  _et al._ (2008); Burke (2012). In short, an embedding
protocol consists in partitioning the original extended system into fragments.
Each fragment, often referred to as an impurity, is complemented with bath
orbitals to define an effective reduced system described by an effective
Hamiltonian. The reduced system, comprising the fragment+bath (denoted the
cluster), contains $N_{i}$ orbitals of the fragment, which interact with the
additional $N_{b}$ bath orbitals. Ideally, the fragment and bath reduced
system is entirely decoupled from the rest of the system, which is referred to
as the environment of the cluster. The cornerstone of these embedding methods
is the construction of bath orbitals associated with a fragment in order to
derive the effective Hamiltonian of the reduced system. This effective
Hamiltonian captures the physics of the larger original environment while
integrating out most of its degrees of freedom. Consequently, it provides a
highly efficient approach to treat localized electron correlation without the
need to explicitly consider the entire system, thus balancing computational
feasibility with accuracy. Different embedding strategies have been proposed
in the literature to construct such a reduced effective Hamiltonian.

Among these embedding strategies, some rely on the Green’s function formalism,
which focuses on single-particle excitations and can naturally incorporate
electronic correlation effects. In particular, within the dynamical mean-field
theory (DMFT) formalism Müller-Hartmann (1989); Georges  _et al._ (1996);
Georges and Kotliar (1992), the effective Hamiltonian is an effective Anderson
impurity model (AIM) Anderson (1961), and is derived self-consistently using
the Green’s function of the system. Note that various alternative embedding
approaches using the Green’s function have been developed to derive effective
AIM Potthoff (2003); Sarker (1988); Mazouin  _et al._ (2019) or other
effective Hamiltonians Nguyen Lan  _et al._ (2016); Lupo  _et al._ (2021).

In parallel, embedding approaches have been proposed within the one-body
Reduced Density-Matrix formalism (1RDM). The aim is to design a protocol to
construct the bath orbitals, and the corresponding effective Hamiltonian of
the cluster, that become functional of the 1RDM. Among the different
approaches recently proposed Lanatà (2023); Senjean (2019); Sekaran  _et al._
(2022); Mitra  _et al._ (2023), the pioneering work of Knizia et al. Knizia
and Chan (2012) proposes to define the effective Hamiltonians by means of the
Schmidt decomposition of a single Slater Determinant (SD) $|\Phi\rangle$ that
is univocally associated with an idempotent 1RDM $\gamma^{0}$, such that the
Schmidt decomposition of $|\Phi\rangle$ defines a projector via a Singular
Value Decomposition (SVD) of the fragment-environment part of $\gamma^{0}$
Wouters  _et al._ (2016); Yalouz  _et al._ (2022); Sekaran  _et al._ (2023).
Indeed, SVD is known to yield the most compact representation of all
informations contained within the fragment-environment 1RDM. This approach
imposes that the bath constructed through the SVD contains as many orbitals as
the impurity fragment, and also that only idempotent 1RDM can be used to
describe the extended system. Over the past decade, different variations and
improvements of the original DMET have been investigated and benchmarked Ayral
_et al._ (2017); Cancès  _et al._ (2023); Sun  _et al._ (2020); Ye  _et al._
(2018); Hermes and Gagliardi (2019). Recently, Sekaran et al. Sekaran  _et
al._ (2021); Yalouz  _et al._ (2022) proposed to use a specific unitary
transformation, defined as a functional of the 1RDM, to construct the bath
orbitals Sekaran  _et al._ (2021). The aforementioned unitary transformation
is known as the Block-Householder transformation. They demonstrated that the
SVD of the fragment-environment 1RDM is equivalent to the Block-Householder
transformation, even for a non-idempotent 1RDM Sekaran  _et al._ (2023).
Despite the appealing compactness of the resulting sub-space of the bath, it
should be noted that there is no inherent physical reason to believe that such
decomposition yields the ”best” effective Hamiltonian for reproducing the
interactions between the fragment and its environment.

In this contribution, we propose a versatile framework that generalizes the
block-Householder transformation and introduces additional degrees of freedom
for the construction of the orbitals of the bath and the derivation of the
effective Hamiltonian. As a result, additional constraints are required for
the bath orbitals such as maximally disentangling the cluster from the
environment, or matching density matrices. The effects of the additional
constraints to optimally construct the orbitals of the bath are benchmarked on
the well-known but non-trivial half-band filled one-dimensional Hubbard model
Hubbard (1126), following the divide and conquer algorithm proposed in a
previous work by the authors Marécat  _et al._ (2023).

##  II Theory

Let us consider the paramagnetic Hubbard ring given by

$\displaystyle\hat{H}=-t\sum_{<ij>\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma},+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}$
(1)

where $\hat{c}^{\dagger}_{i\sigma}$ ($\hat{c}_{i\sigma}$) corresponds to the
creation (annihilation) of an electron of spin $\sigma$ on the $i$-th orbital,
$\hat{n}_{i\sigma}$ is the counting operator equal to
$\hat{c}^{\dagger}_{i\sigma}c_{i\sigma}$. Indexes $<ij>$ refer at nearest
neighbor orbitals. $-t$ corresponds to the hopping integral, while $U$ stands
for the Coulomb integral.

###  II.1 Quantum bath from the Block-Householder transformation

In this section we recall the generic construction of the unitary Block-
Householder transformation $\mathbf{R}^{\sigma}$, following the work of
Sekaran et al. Sekaran  _et al._ (2021). The unitary transformation
$\mathbf{R}^{\sigma}$ is defined as a functional of the spin 1RDM
$\gamma^{\sigma}$. It performs the rotation of the mono-electronic basis set
such that the system is divided into a compact subset of $N_{i}+N_{b}$
orbitals, the cluster, interacting as few as possible with a large number of
$N_{e}$ orbitals. More precisely, the unitary transformation is design such
that (i) the fragment in the full system is the same as in the cluster

$\displaystyle\tilde{c}_{i\sigma}^{\dagger}=\hat{c}_{i\sigma}^{\dagger},\quad\tilde{c}_{i\sigma}=\hat{c}_{i\sigma}\quad\forall
i\in\rm fragment,$ (2)

where

$\displaystyle\tilde{c}_{i\sigma}^{\dagger}=\sum_{k}\mathbf{R}^{\sigma\dagger}_{ik}\hat{c}_{k\sigma}^{\dagger},$
$\displaystyle\quad\tilde{c}_{i\sigma}=\sum_{k}\mathbf{R}^{\sigma}_{ik}\hat{c}_{k\sigma},$
(3)

stand for the creation and annihilation operators of an electron in the $k$-th
orbital with spin $\sigma$ expressed in the $\mathbf{R}^{\sigma}$
representation, and (ii) the fragment is fully disconnected from the
environment at the one-body level,

$\displaystyle\sum_{k}\mathbf{R}^{\sigma}_{ik}\gamma^{\sigma}_{ik}=0\quad\forall
i\in\text{ fragment},\;\forall l\notin\rm cluster.$ (4)

Numerous unitary transformations satisfy conditions (2) and (4) resulting in
identical bath orbitals subspace Töws and Pastor (2011); Schade and Blöchl
(2018). Among all these transformations, we focus explicitly on the Block-
Householder transformation $\mathbf{R}^{\sigma}$ defined using an auxiliary
matrix $\mathbf{V}[\gamma^{\sigma}]\in\mathbb{R}_{NN_{i}}$, with

$\displaystyle\mathbf{R}^{\sigma}=Id.-2\mathbf{V}\left(\mathbf{V}^{T}\mathbf{V}\right)^{-1}\mathbf{V}^{T}.$
(5)

By construction, $\mathbf{R}^{\sigma}$ is a normal involution (i.e.
$\mathbf{R}^{\sigma^{-1}}$=$\mathbf{R}^{\sigma}$,
$\mathbf{R}^{\sigma^{T}}=\mathbf{R}^{\sigma}$ and
det$(\mathbf{R}^{\sigma})=(-1)^{N_{i}}$) with eigenvalues $\\{1,-1\\}$, and
dim(Ker($\mathbf{R}^{\sigma}+Id.$))$=r=N_{i}$. The matrix $\mathbf{V}$ is
given as following

$\mathbf{V}=\left(\begin{array}[]{l}\mathbf{0}^{r}_{N_{i}}\\\
\gamma^{r}_{N_{i}:2N_{i}}+\mathbf{X}^{r}_{N_{i}}\\\
\gamma^{r}_{2N_{i}:N}\end{array}\right).$ (6)

The superscript $\;{}^{r}\;$ (underscript $\;{}_{N_{i}}\;$) stands for th
number of columns (lines) of the associated matrix, and the notation
$\mathbf{A}^{c}_{i:j}$ corresponds to consider only the elements from $i$-th
to the $j$-th lines for first $c$-th columns of the matrix $\mathbf{A}$. For
example, the square matrix $\gamma^{r}_{N_{i}:2N_{i}}$ refers to the first
$r=N_{i}$-th columns of $\gamma^{\sigma}$ from the $N_{i}$-th element to the
$2N_{i}$-th. The null matrix $\mathbf{0}^{r}_{N_{i}}$ allows to fix the
identity in equation (2) on the subspace corresponding to the fragment, and
$\mathbf{X}^{r}_{N_{i}}$ is a square matrix to be determined in order to
satisfy condition (4). To that aim, Rotella et al. Rotella and Zambettakis
(1999) proposed a systematic construction with

$\displaystyle\mathbf{X}^{r}_{N_{i}}=\mathbf{P}^{T}\sqrt{\mathbf{D}}\mathbf{P}\gamma^{r}_{N_{i}},$
(7)

and $\sqrt{\mathbf{D}}={\rm diag}^{N}_{i=1}\\{\sqrt{d_{i}}\\}$, where the
nonnegative scalar $d_{i}$ and the orthogonal matrix $\mathbf{P}$ are defined
by

$\displaystyle
Id.+\left(\gamma^{r}_{2N_{i}:N}\gamma^{r^{-1}}_{N_{i}:2N_{i}}\right)^{T}\left(\gamma^{r}_{2N_{i}:N}\gamma^{r^{-1}}_{N_{i}:2N_{i}}\right)=\mathbf{P}^{T}\mathbf{D}\mathbf{P},$
(8)

and $\mathbf{D}$ ($\mathbf{P}$) is eigenvalues (eigenvectors) matrix of the
left sided matrix in equation (8), respectively. The construction of
$\mathbf{X}^{r}_{N_{i}}$ in equation (7) holds only if the square matrix
$\gamma^{r}_{N_{i}:2N_{i}}$ is invertible. Then we obtain

$\left(\mathbf{R}^{\sigma}\gamma^{\sigma}\right)^{r}=\left(\begin{array}[]{l}\gamma^{r}_{N_{i}}\\\
-\mathbf{X}^{r}_{N_{i}}\\\ \mathbf{0}^{r}_{N_{e}}\end{array}\right).$ (9)

Equation (9) is strictly equivalent as the conditions defined in equations (2)
and (4), i.e. the transformation preserves the orbitals of the fragment in the
cluster and disconnects the fragment from the environment at the one-body
level.

Recent study expounded on the mathematical equivalence between the bath
orbitals constructed using the Block-Householder transformation of the 1RDM
and the SVD of the fragment-environment 1RDM Sekaran  _et al._ (2023). As a
result of this equivalence, the Block-Householder transformation provides the
most compact representation, encapsulating information pertaining to the one-
body level interactions between the fragment and its environment. Once the
Householder transformation has been presented to construct the orbitals of the
bath, we briefly recall the procedure based on $\mathbf{R}^{\sigma}$ in order
to define an effective, yet approximated, embedded Hamiltonian. Details and
discussion on the straightforward generalization to the ab-initio Hamiltonian
can be found in Ref. Marécat  _et al._ (2023). Following equation (3), we
perform the transformation of the full Hamiltonian (1) to derive the
Hamiltonian $\tilde{H}$,

$\displaystyle\tilde{H}=\sum_{ij\sigma}\tilde{t}_{ij}\tilde{c}^{\dagger}_{i\sigma}\tilde{c}_{j\sigma}+\sum_{ijkl}\tilde{U}_{ijkl}\tilde{c}^{\dagger}_{i\uparrow}\tilde{c}_{j\uparrow}\tilde{c}^{\dagger}_{k\downarrow}\tilde{c}_{l\downarrow},$
(10)

where single- and two-bodies integrals in the block-Householder representation
$\tilde{t}_{ij}$ and $\tilde{U}_{ijkl}$, respectively, are calculated in the
same manner,

$\displaystyle\tilde{t}_{ij}$
$\displaystyle=\sum_{kl}t_{kl}\mathbf{R}^{\sigma\dagger}_{ik}\mathbf{R}^{\sigma}_{jl},$
(11) $\displaystyle\tilde{U}_{ijkl}$
$\displaystyle=U\sum_{m}\mathbf{R}^{\sigma\dagger}_{im}\mathbf{R}^{\sigma}_{jm}\mathbf{R}^{\bar{\sigma}\dagger}_{km}\mathbf{R}^{\bar{\sigma}}_{lm}.$
(12)

Considering that $\mathbf{R}^{\sigma}$ maximally uncouples the fragment from
the environment though the orbitals of the bath, an approximation for the
effective Hamiltonian on the cluster $\bar{H}^{c}$ is obtained by projecting
the Hamiltonian (10) onto the cluster orbitals. Part of the cluster-
environment two-body interactions are taken into account at a mean-field level
in $\bar{H}^{ce}_{\rm MF}$. Finally, an effective chemical potential $\mu_{\rm
emb}$ is added to preserve the total number of electrons in the fragment. All
together, the cluster effective Hamiltonian reads

$\bar{H}^{c}_{\rm eff}=\bar{H}^{c}+\bar{H}^{ce}_{\rm MF}+\mu_{\rm
emb}\sum_{i\in{\rm bath}}\tilde{c}_{i}^{\dagger}\tilde{c}_{i},$ (13)

and

$\displaystyle\bar{H}^{ce}_{\rm MF}[\gamma^{\sigma}]=$
$\displaystyle\sum_{i,j\in{\rm
cluster}\atop\sigma}\tilde{c}^{\dagger}_{i\sigma}\tilde{c}_{j\sigma}\sum_{(kl)}\tilde{U}_{ijkl}\sum_{mn}\mathbf{R}^{\bar{\sigma}\dagger}_{km}\mathbf{R}^{\bar{\sigma}}_{nl}\gamma^{\bar{\sigma}}_{kl}$
$\displaystyle+{\rm h.c.},$ (14)

where the notation $(kl)$ refers to pairs of orbitals in $\mathbf{R}^{\sigma}$
representation such that at least $k$ or $l$ belongs to the environment.
Likely influenced by the non-interacting character of the bath in DMFT, the
non-locals interaction integrals $\tilde{U}_{ijkl}$ that naturally arise in
the bath have been initially neglected in DMET, leading to a Anderson impurity
model in the cluster Knizia and Chan (2012). This approximation is called non-
interacting bath (NIB) in contrast to the interacting bath (IB) version of
DMET where $\tilde{U}_{ijkl}$ are explicitly taken into account Wouters  _et
al._ (2016). The similar distinction between NIB an IB can also be considered
using the Block-Householder transformation to construct the orbitals of the
bath Sekaran  _et al._ (2021); Marécat  _et al._ (2023).

Interestingly for non-interacting Hamiltonian ($U=0$), for which the
associated 1RDM is idempotent ($\gamma^{\sigma 2}=\gamma^{\sigma}$), the
transformation $\mathbf{R}^{\sigma}$ leads to a perfect decoupling of 1RDMs of
the clusters (Eq. (4) being fulfilled $\forall i\in$ the cluster instead of
the fragment), respectively. It corresponds to an exact factorization of the
underlying associated wave-function $|\Psi\rangle=|\phi\rangle$, where
$|\phi\rangle$ refers to a single-Slater determinant. The factorization of the
wave-function gives

$\mathbf{R}|\phi\rangle=\hat{\mathcal{A}}|\tilde{\phi}^{c}\rangle|\tilde{\phi}^{e}\rangle,$
(15)

where $|\tilde{\phi}^{c}\rangle$ ($|\tilde{\phi}^{e}\rangle$) is an anti-
symmetrized product of orbitals that belong solely to the cluster
(environment). Consequently, in this specific case, the projected cluster
Hamiltonian proposed in equation (13) can be used to extract exact local
properties of the impurity site, where the ground-state wave function
$|\bar{\phi}\rangle$ is equal to $|\tilde{\phi}^{c}\rangle$ Sekaran  _et al._
(2021).

In a more general context, particularly with interacting cases, the effective
cluster Hamiltonian, as defined in equation (13), can serve as a useful
approximation for local properties of the fragment. This Hamiltonian is de
facto a functional of the 1RDM via the definition of
$\mathbf{R}^{\sigma}[\gamma^{\sigma}]$. Several studies have developed self-
consistent schemes predicated on the local cluster 1RDM calculated with
equation (13), as seen in DMET Knizia and Chan (2012) and more recent divide
and conquer algorithms Marécat  _et al._ (2023). Regardless of the self-
consistent matching or conquer strategy employed, the efficiency of the method
hinges crucially on the ability of the effective cluster Hamiltonian (13) to
locally mimic the full Hamiltonian on the fragment. On this subject, there is
no unambiguously definitive approach to determine what exactly constitutes
”mimicry”. To address this challenge, we propose in the following sections to
generalize equation (6), introducing variational parameters to explore
different flavours of embedding.

###  II.2 Quantum bath from a versatile unitary transformation framework

Following the philosophy of exact diagonalization solver in DMFT Caffarel and
Krauth (0307), we would like to control the number of bath orbitals,
independently of the number of impurity orbitals in the fragment, in order to
systematically have a better description of the interactions of the fragment
with the environment. As shown schematically in Fig. 1, the Block-Householder
transformation leads the fragment (dark blue square) unchanged in the cluster
(gray square), and interacts solely with the bath (orange square), where the
number of bath orbitals is the same as in the fragment. In the following, we
give a general and flexible definition of the unitary matrix
$\mathbf{R^{\prime}}^{\sigma}[\gamma^{\sigma}]$ to obtain an optimized
cluster, composed of $N_{i}$ impurities coupled to $N_{b}\geq N_{i}$ bath
orbitals.

Figure 1: Schematic representation of the Block-Householder transformation and
the generalization, both functionals of the 1RDM. Small dark blue squares
correspond to the fragment 1RDM in the system, identical as in the cluster,
and orange squares correspond to bath 1RDM in the cluster. The gray squares
correspond to the cluster.

In what follows, spin index $\sigma$ is omitted for clarity. Similarly as the
Block-Householder transformation, we use the auxiliary matrix
$\mathbf{V^{\prime}}$ with

$\displaystyle\mathbf{R^{\prime}}^{\sigma}=Id.-2\mathbf{V^{\prime}}\left(\mathbf{V^{\prime}}^{T}\mathbf{V^{\prime}}\right)^{-1}\mathbf{V^{\prime}}^{T}.$
(16)

The auxiliary matrix $\mathbf{V^{\prime}}\in\mathbb{R}_{Nr}$, where $N_{i}\leq
r\leq N$ refers to the rank of the auxiliary matrix, is constructed as
following

$\displaystyle\mathbf{V^{\prime}}=\left(\begin{array}[]{l}\mathbf{0}^{r}_{N_{i}}\\\
\gamma^{r}_{N_{b}}-\mathbf{X^{\prime}}^{r}_{N_{b}}\\\
\gamma^{r}_{N_{e}}\end{array}\right),$ (20)

and indexes $N_{i}$, $N_{b}$, $N_{e}$ refer to the number of impurity orbitals
in the fragment, the number of bath orbitals and the number of orbitals in the
environment, respectively. $\gamma^{r}_{N_{b}}$ ($\gamma^{r}_{N_{e}}$)
corresponds to the $r$ first columns of the spin 1RDM $\gamma^{\sigma}$ and
the $N_{b}$ lines ($N_{e}$) after the $N_{i}$-th ($N_{i}+N_{b}$-th) first one.
At this point, condition (2) (preservation of the fragment in the cluster) is
already satisfied. The matrix $\mathbf{X^{\prime}}^{r}_{N_{b}}$ must be
determined in order to preserve condition (4), i.e. it must satisfy the
following equation,

$\displaystyle(\mathbf{X^{\prime}}^{r}_{N_{b}}-\gamma^{r}_{N_{b}})^{T}(\mathbf{X^{\prime}}^{N_{i}}_{N_{b}}+\gamma^{N_{i}}_{N_{b}})=\gamma^{r^{T}}_{N_{e}}\gamma^{N_{i}}_{N_{e}},$
(21)

where the superscript $N_{i}$ refers to the $N_{i}$ first column of the
matrix. In pratice, the nonlinear equation (21) can be solved numerically. In
the following, we discuss only the single impurity case $N_{i}=1$, where
equation (21) can be solved analytically. In this case, the matrix
$\mathbf{X^{\prime}}^{N_{i}}_{N_{b}}$ is reduced to a vector, and equation
(21) becomes

$\displaystyle(\mathbf{X^{\prime}}^{1}_{N_{b}}-\gamma^{1}_{N_{b}})^{T}(\mathbf{X^{\prime}}^{1}_{N_{b}}+\gamma^{1}_{N_{b}})=\gamma^{1^{T}}_{N_{e}}\gamma^{1}_{N_{e}}$
(22)
$\displaystyle(\mathbf{X^{\prime}}^{j}_{N_{b}}-\gamma^{j}_{N_{b}})^{T}(\mathbf{X^{\prime}}^{1}_{N_{b}}+\gamma^{1}_{N_{b}})=\gamma^{j^{T}}_{N_{e}}\gamma^{1}_{N_{e}},\,\forall\,1<j\leq
r.$ (23)

Equation (22) is nonlinear and constrains the norm of the vector
$\mathbf{X^{\prime}}^{1}_{N_{b}}$. More precisely, this equation is fulfilled
for any vectors $\mathbf{X^{\prime}}^{1}_{N_{b}}$ satisfying

$||\mathbf{X^{\prime}}^{1}_{N_{b}}||=||\gamma^{1}_{N_{b}}||+||\gamma^{1}_{N_{e}}||,$
(24)

where $||\mathbf{v}||=\mathbf{v}^{T}\mathbf{v}$ refers to the norm of the
vector $\mathbf{v}$. Equation (23) exists only for $r>N_{i}$ and is linear for
$\mathbf{X^{\prime}}^{j}_{N_{b}}$, corresponding to a scalar-product
preservation, and can be written as follows,

$\displaystyle\mathbf{X^{\prime}}^{j^{T}}_{N_{b}}\mathbf{X^{\prime}}^{1}_{N_{b}}=\gamma^{j^{T}}_{N_{b}}(\mathbf{X^{\prime}}^{1}_{N_{b}}+\gamma^{1}_{N_{b}})+||\gamma^{1}_{N_{e}}||.$
(25)

We introduce a spherical representation which allows to express the vectors
$\mathbf{X^{\prime}}^{j}_{N_{b}}$ with lengths $l^{j}$ and a complete set of
angles $\\{\theta\\}$. In this representation, length $l^{1}$ ($l^{j}$) is
used to fulfill norm preservation (24) (scalar product (25)) for
$\mathbf{X^{\prime}}^{1}_{N_{b}}$ ($\mathbf{X^{\prime}}^{j}_{N_{b}}$)
respectively. The set of angles $\\{\theta\\}$ are thus completely free and
can take any value between $[-\pi,\pi[$. Consequently, we have a complete set
of $\\{\theta\\}$ parameters to construct the corresponding set of different
auxiliary matrices defined in equation (20) satisfying the conditions (2) and
(4). The number of free parameters is equal to $r\times(N_{b}-N_{i})$, with
$r\geq N_{i}$. In Fig. 2, we illustrate all vectors that can be obtain for
$N_{i}=1$, $N_{b}=3$ and $r=2$. Dark blue vectors correspond to the special
Block-Householder solution, where $N_{b}=1$ corresponds to the second impurity
in fragment. The rank two Block-Householder transformation preserves the
identity over the second impurity. As a result, the projection of the first
(second) Block-Householder vector over the second impurity $i=2$ axis gives
$\gamma^{\sigma}_{12}$ ($\gamma^{\sigma}_{22}$) , as shown with the dashed
dark blue line in Fig. 2. Beyond the special Block-Householder transformation,
all vectors belonging to the orange sphere are solution of equation (24) such
as presented with $\mathbf{X^{\prime}}^{1}_{N_{b}}$ (orange vector) for
example. The second vector (light blue vector) norm depends on its scalar
product with $\mathbf{X^{\prime}}^{1}_{N_{b}}$. Thus the norm (and the
direction) is fixed using equation (25) and all other angles left are free.
More generally, when $r>1$, every vectors $\mathbf{X^{\prime}}^{j}_{N_{b}}$ ,
$1<j\leq r$ are independent from each other and are correlated to
$\mathbf{X^{\prime}}^{1}_{N_{b}}$ solely. In the case of $r=N_{i}$, only one
solution is available (with 0 degrees of freedom) and leads to the unique
Block-Householder transformation.

Note that by considering $r$, $N_{b}$ and $N^{\prime}_{b}>N_{b}$, the
auxiliary matrix $\mathbf{V^{\prime}}[\mathbf{X^{\prime}}^{r}_{N_{b}}]$ space
is not included into
$\mathbf{V^{\prime}}[\mathbf{X^{\prime}}^{r}_{N^{\prime}_{b}}]$, i.e. any rank
$r$ auxiliary matrices expressed using equation (20) with $N_{b}$ bath cannot
be expressed with an auxiliary matrix with a greater number of bath
$N^{\prime}_{b}$. Similarly, if we consider $N_{b}$, $r$ and $r^{\prime}>r$,
we get
$\mathbf{V^{\prime}}[\mathbf{X^{\prime}}^{r}_{N_{b}}]\notin\mathbf{V^{\prime}}[\mathbf{X^{\prime}}^{r^{\prime}}_{N_{b}}]$,
which means that different ranks $r$ lead to a a specific definition of the
cluster.

$N_{b}=2$$N_{b}=3$$N_{b}=1$$\mathbf{X^{\prime}}^{1}_{N_{b}}$$l^{1}$: Eq.
(24)$\mathbf{X^{\prime}}^{2}_{N_{b}}$$l^{2}$[$\mathbf{X^{\prime}}^{1}_{N_{b}}$]:
Eq. (25)$\theta^{1}_{2}$$\theta^{1}_{1}$

Figure 2: Schematic representation of vectors
$\mathbf{X^{\prime}}^{j}_{N_{b}}$ (light blue and orange vectors) in spherical
representation for $N_{i}=1$, $N_{b}=3$ and $r=2$. Dark blue vectors
correspond to exact Block-Householder vectors from a density- matrix of a non-
interacting 1D-Hubbard model. The length $l^{1}$ is fixed according to
equation (24), while lengths $l^{j}$, $r\geq j>1$ depends on the vector
$\mathbf{X^{\prime}}^{1}_{N_{b}}$ following the equation (25). All angles
parameters $\\{\theta^{j}_{i}\\}$, $1\leq i\leq r-N_{i}$ are free. For
example, all vectors belonging to the orange sphere are available as a choice
of $\mathbf{X^{\prime}}^{1}_{N_{b}}$. The number of free parameters is equal
to $r\times(N_{b}-N_{i})=4$ (orange and light blue dashed curves).

At this stage we have proposed a generic construction of unitary
transformations following equations (16) and (20) that generalize the block-
Householder construction. Indeed, given an arbitrary number of bath orbitals
$N_{b}$ ($N_{b}\geq N_{i}$) and rank $r$ of matrix $\mathbf{V^{\prime}}$
($N_{i}\leq r\leq N$), we show that we can construct many unitary
transformations that fulfill conditions (2) and (4) up to
$r\times(N_{b}-N_{i})$ free parameters. In the next section we propose to use
these additional parameters in order to add physically motivated criteria to
design bath orbitals.

###  II.3 Optimization of free parameters

The free parameters $\\{\theta\\}$ are variationaly optimized to adjust the
bath orbitals using physical insights, unlike the systematic Block-Householder
transformation. In what follows, the transformation
$\mathbf{R^{\prime}}^{\sigma}$ is a functional of the spin 1RDM, but also a
function of a full set of free parameters $\\{\theta\\}$ and is denoted as
$\mathbf{R^{\prime}}^{\sigma}[\gamma^{\sigma}](\\{\theta\\})$, or in a more
compact notation $\mathbf{R^{\prime}}^{\sigma}(\theta)$. The transformation of
$\gamma^{\sigma}$ with the transformation
$\mathbf{R^{\prime}}^{\sigma}(\theta)$ is denoted as
$\tilde{\gamma}^{\sigma}$.

According to the construction of $\mathbf{R^{\prime}}^{\sigma}(\theta)$
fulfilling constraint (4), the fragment is disconnected from the environment
at the one-body level. In the non-interacting case, we show analytically that
the Block-Householder ($r=N_{i}$) disconnects bath orbitals from the
environment. However, this is not the case for the correlated 1RDM. Therefore,
the set of variational parameters is used to minimize the value of the so-
called buffer-zone $\Delta$, which gives a quantitative insight into the
disentanglement of the environment cluster at the one body level and leads to
the saddle-point equation

$\displaystyle\Delta^{2}(\theta)=\sum_{b\in\rm bath}\sum_{e\in\rm
env}\tilde{\gamma}^{\sigma^{2}}_{ce},$ (26)
$\displaystyle\dfrac{\partial\Delta(\theta_{i})}{\partial\theta_{i}}=0\quad\forall\theta_{i}\in\\{\theta\\}.$
(27)

Following a similar philosophy, one could minimize the single-particle Von-
Neumann entropy of the truncated 1RDM of the cluster.

From the medium to the strong correlated regime, the buffer-zone $\Delta$ is
not able to give a quantitative value of the entanglement between the cluster
and the environment, where two-body interactions dominate at this regime. In
the following, we propose to minimize the square of the Hartree (i.e mean-
field) contribution energy between the cluster and the environment

$\displaystyle E^{2}_{\rm H}(\theta)=\sum_{c\in{\rm
cluster}\atop\sigma}\left[\sum_{(jkl)}\left(\tilde{U}(\theta)_{cjkl}\tilde{\gamma}^{\sigma}_{cj}\tilde{\gamma}^{\sigma}_{kl}\right)^{2}\right.$
(28) $\displaystyle\left.+\sum_{e\in{\rm
env}}\left(\tilde{t}_{ce}\tilde{\gamma}^{\sigma}_{ce}\right)^{2}\right],$ (29)
$\displaystyle\dfrac{\partial E^{2}_{\rm
H}(\theta_{i})}{\partial\theta_{i}}=0\quad\forall\theta_{i}\in\\{\theta\\},$
(30)

where $c$ ($e$) belongs to the cluster (environment), respectively, and the
notation $(jkl)$ refers to pairs of orbitals in
$\mathbf{R^{\prime}}^{\sigma}(\theta)$ representation such that at least $j$,
$k$ or $l$ belongs to the environment, and $\tilde{t}_{ce}$
($\tilde{U}(\theta)_{cjkl}$) defined in equation (11) (equation (12)),
respectively. In a practical way, the evaluation of $E^{2}_{\rm H}$ using
equation (28) is numerically more expensive than the evaluation of $\Delta$.

Finally, inspired by the DMET matching, we propose to design the
transformation in order to enforce the matching between density matrix
elements connected to the fragment in the transformation space, and in the
cluster

$\underset{\theta}{\rm
min}\sum_{ij}\left(\tilde{\gamma}^{\sigma}_{ij}-\bar{\gamma}^{\sigma
c}_{ij}\right)^{2}\quad\forall i\in{\rm fragment,}\forall j\in{\rm cluster},$
(31)

with $\bar{\gamma}^{\sigma c}$ refers as the ground-state spin $\sigma$ 1RDM
of the cluster obtained by solving $\bar{H}^{c}_{eff}$ defined in equation
(13).

As illustrated in Fig. 3, we test the cost functions proposed in equations
(26), (28) and (31) with respect to two free parameters $\\{\theta\\}$ for
$N_{i}=1$, $N_{b}=3$ and $r=2$. A non-interacting one body reduced density
matrix (idempotent) is used as a 1RDM test to obtain the transformation
$\mathbf{R^{\prime}}^{\sigma}(\theta)$. The first rank $r=1$ vector is fixed
using Block-Householder solution. In that case, two free parameters are
optimized, and cost functions can be represented in spherical coordinates. The
minimization of the buffer zone $\Delta$ (Fig. 3(a)) strictly cancels the
value of $\Delta$ and corresponds to the Block-Householder solution, as
discussed in Sec. II.1. This result is attributed to the fact that the trial
density matrix is idempotent. However, this particular value is enclaved
between regions of higher $\Delta$ values, which can make optimization
difficult depending on the starting point of the numerical minimization.
Subsequently, the minimization of the mean-field term between the cluster and
the environment (Fig. 3(b)) also yields the same solution in this case.
However, the landscape is distinct from the previously studied cost function.
Finally, the matching of the density matrices (Fig. 3(c)) presents a very
specific landscape. Indeed, such a landscape is numerically very challenging
to explore in order to obtain the global minimum. Contrary to the cost
functions studied previously, the Block-Householder solution (dark blue
triangle) is not the global minimum. It appears that there exists a continuous
set of minimums, which further complicates the exploration of the landscape.
Altogether, it highlights the non-trivial character of the resulting landscape
that might display many local minima and quasi-flat regions. It results that
the numerical optimization of the $\\{\theta\\}$ parameters might become
challenging for large amount of parameters.

(a)

(b)

(c)

Evaluation of the cost functions with respect to free parameters
$\\{\theta\\}$ for $N_{i}=1$, $N_{b}=3$ and $r=2$ using a non-interacting one
body reduced density matrix from a 1D Hubbard chain. The vector corresponding
to the rank $r=1$ is fixed and corresponds to the Block-Householder vector
(see Fig. 2), and the second is presented with the dark blue triangle. In this
case, two angles are available, and are represented using spherical
coordinates. In (a) we evaluate the buffer-zone $\Delta$ (see Eq. (26)). In
(b), we evaluate the Hartree contribution $E^{2}_{\rm H}$ (see Eq. (28)) for
$U/t$ = 8, and in (c) we evaluate the density matrix matching (see Eq. (31))
for $U/t$ = 8.

Figure 3:

##  III Results and Discussion

In this section, the various cost functions outlined in Section II.3 are
evaluated against the homogeneous paramagnetic Hubbard model at half-filling,
and compared to exact Bethe-Ansatz (BA) results Lieb and Wu (1968); Ogata and
Shiba (1990). The non-idempotent N-representable 1RDM space is spanned using
the self-consistent protocol presented in Marécat  _et al._ (2023). The
various results pertain to the case of a single impurity. The generalization
to multiple impurities, discussed in Section II.2, will not be covered in this
paper. We recall that the Block-Householder transformation corresponds to a
particular case, where the set of parameters $\\{\theta\\}$ is set according
to the equation (7). Additionally, the variational parameters are optimized
using the numpy python library Virtanen  _et al._ (2020), particularly with
the L-BFGS-B method, which is similar to the conjugate gradient optimization
method. Finally, as explained in Marécat  _et al._ (2023), a damping parameter
is added to avoid drastic changes of the 1RDM and convergence issues, meaning
that a fraction of the previous 1RDM obtained is kept in the new 1RDM. In the
results presented here we used a damping of $60\%$.

In Fig. 4, we show the relative error of the kinetic energy $\Delta
E_{k}=100\times\left(E^{\rm BA}_{K}-E_{K}\right)/E^{\rm BA}_{K}$, where BA
refers to the Bethe Ansatz solution, and the relative error of the double
occupation $\Delta d=100\times\left(d^{\rm BA}-d\right)/d^{\rm BA}$ as a
function of the relative correlation strength $U/(U+4t)$, where $4t$
corresponds to the non-interacting band width. This figure presents the effect
of the various cost functions outlined in Sec. II.3 (color-coded lines). In
this case, the single impurity $N_{i}=1$ is embedded with three IB orbitals
$N_{b}=3$ for a vector $\mathbf{V^{\prime}}$ of rank $r=2$. Thus, there is a
number of parameters ${\theta}$ equal to $2\times(3-1)=4$ to optimize. In the
non-interacting limit $U/t\rightarrow 0$, the kinetic energy and the double
occupation are correctly reproduced for all cost functions. In regards to the
atomic limit $U/t\rightarrow\infty$, this is also the case with an asymptotic
limit of $-8$ln$(2)t^{2}/U$ for the kinetic energy, and $4$ln$(2)t^{2}/U^{2}$
for the double occupation López-Sandoval and Pastor (2002). For intermediate
regimes $U/4t\simeq 1$, there is a strong competition between electronic
delocalization which increases kinetic energy, and the electron-electron
repulsion strength which penalises the number of double occupation. For state-
of-the-art embedding methods, such as DMET Knizia and Chan (2012) or the
projected site-occupation embedding theory (PSOET) Senjean (2019), describing
accuratly this regime is very challenging. For weakly to intermediate
correlated regimes ($U/(U+4t)<0.6$, $U/t<6$), the results from the
minimization of constrain (28) (dark blue line) are consistently better than
those obtained by minimizing the constrain (26) (orange line), which are in
turn better than the Block-Householder solutions (yellow line). However, the
numerical cost associated with the minimization of constrain (28) scale as
$\mathcal{O}\left(N_{c}N^{4}_{e}\right)$, which is significantly higher than
the cost associated with constrain (26) scaling as
$\mathcal{O}\left(N_{c}N^{2}_{e}\right)$, where $N_{e}$ ($N_{c}$) represents
the number of orbitals belonging to the environment of the cluster (cluster),
respectively. Concerning the one body reduced density matrix matching proposed
in equation (31) (blue line), the results obtained are similar to the
constrain (28) from low to middle correlated regime, and deviate for values of
$U/(U+4t)\simeq 0.6$ for kinetic energy and double occupancy. For strongly
correlated regimes ($U/(U+4t)\geq 0.6$), the interaction energy $E_{\rm
int}=U\times d$ dominates. Although the double occupancy is similar for the
Block-Householder method and the buffer-zone minimization at this regime, the
minimization of the mean-field term improves the results. The double occupancy
is badly described by the density matrix matching and exhibits nonphysical
numerical instabilities for values of $U/(U+4t)\simeq 1$.

Figure 4: Relative error for the kinetic energy $\Delta E_{k}$ (top panel) and
per site double occupation $\Delta d$ (bottom panel) in percent with respect
to correlation strength $U/(U+4t)$ for one impurity orbital in the fragment,
three bath orbitals and a rank two vector. Colored lines correspond to
different cost functions. Black dashed line corresponds to Bethe Ansatz.

In Fig. 5, we show the evolution of the different variational parameters
${\theta}$ for different cost functions as a function of relative correlation
strengh $U/(U+4t)$ for a two rank $r$ and a three bath orbitals calculation.
In the left panel, we show the two variational parameters associated with rank
$r=1$, i.e., the first column of the matrix $\mathbf{V^{\prime}}$ in equation
(20), while in the right panel, the parameters are associated with rank $r=2$.
Solid lines correspond to the first angles for all ranks, and dotted lines to
the second ones. We found that all parameters except for the first of rank two
(solid line in the right panel) do not change significantly with respect to
$U/t$. Interpreting this result is challenging. It might be due to specific
symmetries of the system, such as translation invariance and electron-hole
symmetry at half-filling. Regarding the first parameters of rank two, they all
follow the same trend, with the Hartree cost function slightly lower than the
others. In light of these results, it is conceivable to simplify the
variational optimization of the different parameters by considering only a
reduced set of parameters (in this case the first of the second rank) varying
significantly in the process. This simplification could greatly improve the
numerical optimization of the parameters, and therefore the numerical
efficiency of the method in general.

Figure 5: Optimized free angle parameters with respect to the relative
correlation strengh $U/(U+4t)$ for a rank $r=2$ and a three bath orbitals
calculation, for different cost functions (colored dashed lines). The left
panel corresponds to the first rank, while the right panel to the second.
Solid lines correspond to the first parameter of each rank, and dashed lines
to the second.

As originally presented in DMET Knizia and Chan (2012), the effective
Hamiltonian was an AIM, which implies that interaction terms within the bath
are neglected, a scenario referred to as the NIB approximation. In this
context, we present in Fig. 6 the relative error of the ground-state energy
per site $E_{gs}=E_{k}+Ud$, with respect to the correlation strength for a
system with a single impurity in the fragment, three bath orbitals, and a
rank-two vector. In the low correlation regime ($U/(U+4t)<0.3$), both NIB
(represented by dashed colored lines) and IB (represented by full colored
lines) yield similar results across all presented cost functions. However, for
the intermediate to strong correlation regime ($U/(U+4t)>0.6$), the NIB
approximation fails to provide an accurate description of the ground-state
energy of the system. It should be noted that the Block IB (full yellow line)
appears to offer the most accurate representation of the ground-state,
contradicting the observations made in Fig. 4. However, this apparent accuracy
is misleading and arises from a larger error compensation between the kinetic
and interaction energies. Additionally, we found the same error compensation
for the NIB approximation, leading to inadequate results for both kinetic and
interaction energies. It is also noteworthy to optimize variational parameters
to minimize the strength of Coulomb repulsion in the bath. This method could
potentially close the gap between finite interacting baths in DMET or DaC
approaches, and the infinite but non-interacting baths present in DMFT.

Figure 6: Relative error of the ground state energy per site $\Delta
E_{gs}=100\times(1-E_{gs}/E_{BA})$ with respect to correlation strength
$U/(U+4t)$ for one impurity and three bath orbitals and a rank two vector.
Colored lines correspond to different cost functions for the IB case, and
dashed lines for the NIB case. Black solid line corresponds to Bethe Ansatz.

We present in Fig. 7 relevant optimized parameter (the first parameters of
rank two, see results presented in Fig. 5) with respect to the relative
correlation strengh for different cost functions (colored lines). The two
scenarios, NIB (dashed lines) and the IB (full lines), are presented. For the
NIB approximation, optimized parameters follow the same trend as the IB
approximation, but are overestimated for all relative correlation strengh. As
explained in Fig. 6, both NIB and IB yield similar results across all
presented cost functions in the low correlated regime and equivalent
parameters for both non-interacting ($U/t\rightarrow 0$) and atomic
($U/t\rightarrow\infty$) limits. Interestingly, the interaction terms in the
bath (IB), that emerge naturally with an embedding scheme based on a unitary
transformation, are mandatory to obtain a correct description, in particular
at large $U/t$.

Figure 7: Optimized free angle parameters with respect to the relative
correlation strengh $U/(U+4t)$ for the first angle of the rank two vector of
the set of parameters $\\{\theta\\}$, for different cost functions (colored
lines). Solid lines correspond to the IB case, and dashed lines to the NIB
case.

In the following, we focus on the IB case and the minimization of the buffer
zone (see Eq. (26)) to explore the influence of the number of bath orbitals
and the rank of the matrix $\mathbf{V^{\prime}}$ used to define the unitary
transformation. In Fig. 8, we present the per site kinetic energy scaled with
the non-interacting kinetic energy per site $E^{0}_{K}=-4/\pi$ (upper panel)
and the per site double occupation (lower panel) as a function of relative
repulsion strength $U/(U+4t)$. Results correspond to the cases $N_{i}=1$,
$N_{b}=3$ and are given for different rank $r$ of the vector
$\mathbf{V^{\prime}}$, ranging from one up to three. For instance, the rank
two (blue line) corresponds to the results previously shown in Fig. 4 (orange
line). The number of variational parameters is equal to $2$ for $r=1$, equal
to $4$ for $r=2$, and equal to $6$ for $r=3$. For a rank $1$ vector, we are in
the special case where the rank equals the number of impurities $N_{i}$. In
this case, the construction of vector $\mathbf{V^{\prime}}$ needs to satisfy
only the norm preservation in equation (24). As demonstrated previously in
Section II.2, increasing the rank of vector $\mathbf{V^{\prime}}$ does not
systematically improve the solutions, as the accessible solution spaces are
disjoint. Indeed, we show that the rank $2$ results are the closest to the
exact results for both kinetic energy and double occupation, and this applies
to all correlation regimes $U/t$. For strongly correlated regimes, rank $3$
slightly improves the results of rank $1$, but it exhibits numerical
instabilities due to the optimization of a larger amount of variational
parameters. According to this figure, it is not necessary to increase the rank
of vector $\mathbf{V^{\prime}}$ in order to systematically improve the
results. Moreover, the best results are obtained for $r=2$, corresponding to
the number of singular values in DMET Knizia and Chan (2012), or the number of
columns of the vector $\mathbf{V^{\prime}}=\mathbf{V}$ for the Block-
Householder method presented here Marécat  _et al._ (2023).

Figure 8: Renormalized kinetic energy per site $E_{k}/E^{0}_{k}$ (top panel)
and per site double occupation $d$ (bottom panel) with respect to correlation
strength $U/(U+4t)$ for one impurity and three bath orbitals and the
minimization of the buffer zone as a cost function. Colored lines correspond
to the rank of the vector up to three. Black solid line correspond to Bethe
Ansatz.

In Fig. 9, we present the scaled kinetic energy per site (upper panel) and the
per site double occupation (lower panel) as a function of relative repulsion
strength $U/(U+4t)$. The results are shown for a rank $r=3$ vector
$\mathbf{V^{\prime}}$, for different number of orbitals in the bath, ranging
from one up to five (colored lines). For spin symmetry reasons, we only
consider cases where the number of orbitals in the cluster is even. The case
with a single orbital in the bath (yellow line) is very particular, as there
are no variational parameters to be optimized in this case. In the other
cases, the number of variational parameters corresponds to $6$ for $N_{b}=3$
orbitals, and $12$ for $N_{b}=5$ orbitals. Importantly, increasing the number
of orbitals in the bath systematically improves the kinetic energy and double
occupation for all correlation regimes $U/t$. However, for strongly correlated
regimes $U/t>>1$, we observe oscillations of the solutions for five bath
orbitals. In the latter case, there are a large number of variational
parameters, and their optimization is numerically challenging with the
proposed iterative process.

Figure 9: Renormalized kinetic energy per site $E_{k}/E^{0}_{k}$ (top panel)
and per site double occupation $d$ (bottom panel) with respect to correlation
strength $U/(U+4t)$ for a rank three vector and the minimization of the buffer
zone as a cost function. Colored lines correspond to the number of bath
orbitals up to five. Black solid line corresponds to Bethe Ansatz.

##  IV Conclusion

In this study, we have thoroughly investigated the performance of various 1RDM
based embedding methods to construct the orbitals of the bath. A particular
emphasis is placed on the characteristics of the resulting reduced and
effective Hamiltonian. Indeed this Hamiltonian is tasked with accurately
reproducing the interactions between the fragment of the system and its
environment within a downscaled cluster.

While DMET employs the SVD of the fragment-environment 1RDM to define the
effective Hamiltonian, we have demonstrated that the compact subspace is not
the optimal setting for deriving the effective Hamiltonian. By generalizing
the Block-Householder equations, we introduce a significant amount of
additional flexible parameters, notably by adding bath orbitals that are
nearly independent of the number of fragment orbitals or by exploring
different transformation domains via rank augmentation. To efficiently
leverage these additional degrees of freedom, we have proposed cost functions
that, in most cases, effectively disconnect the cluster, containing an integer
number of electrons, from the environment.

These cost functions were tested on the half-filled Hubbard model Hamiltonian,
with a single impurity orbital in the fragment for which the equations are
simplified. The results showed significant improvements over the Block-
Householder outcomes.

Nevertheless, these improvements imply numerical optimization, which often
proves challenging due to the complex landscape of cost functions. The
complexity of these landscapes likely contributes to the fact that we can
currently only achieve half-filled results. Moreover, this complexity
occasionally makes it difficult to obtain continuous solutions for all
relative correlation strengths, resulting in certain non-physical
instabilities. Therefore, we encourage further research into the development
of efficient cost functions that can derive an optimized effective Hamiltonian
to describe the fragment, offering a smoother landscape than its counterparts.
A compelling challenge for future research would be to test the method on
multiple-impurity fragments and, importantly, to derive linearized equations
to define the unitary transformation more effectively.

###### Acknowledgements.

The authors would like to thank the ANR (Grant No. ANR-19-CE29-0002 DESCARTES
project) for funding.

## References

  * Sun and Chan (2016) Q. Sun and G. K.-L. Chan, Acc. Chem. Res. 49, 2705 (2016).
  * Wasserman and Pavanello (2020) A. Wasserman and M. Pavanello, Int. J. Quantum Chem. 120, e26495 (2020).
  * Kotliar  _et al._ (2001) G. Kotliar, S. Y. Savrasov, G. Pálsson, and G. Biroli, Phys. Rev. Lett. 87, 186401 (2001).
  * Zheng (1993) H. Zheng, Phys. Rev. B 48, 14868 (1993).
  * Ma  _et al._ (2021) H. Ma, N. Sheng, M. Govoni, and G. Galli, J. Chem. Theory Comput. 17, 2116 (2021).
  * Cohen  _et al._ (2008) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008).
  * Burke (2012) K. Burke, J. Chem. Phys. 136, 150901 (2012).
  * Müller-Hartmann (1989) E. Müller-Hartmann, Zeitschrift für Physik B Condensed Matter 74, 507 (1989).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Georges and Kotliar (1992) A. Georges and G. Kotliar, Phys. Rev. B 45, 6479 (1992).
  * Anderson (1961) P. W. Anderson, Phys. Rev. 124, 41 (1961).
  * Potthoff (2003) M. Potthoff, Eur. Phys. J. B 32, 429 (2003).
  * Sarker (1988) S. Sarker, J. Phys.: Condens. Matter 21, L667 (1988).
  * Mazouin  _et al._ (2019) L. Mazouin, M. Saubanère, and E. Fromager, Phys. Rev. B 100, 195104 (2019).
  * Nguyen Lan  _et al._ (2016) T. Nguyen Lan, A. A. Kananenka, and D. Zgid, J. Chem. Theory Comput. 12, 4856 (2016).
  * Lupo  _et al._ (2021) C. Lupo, F. Jamet, T. Tse, I. Rungger, and C. Weber, Nat. Comput. Sci. 1, 410 (2021).
  * Lanatà (2023) N. Lanatà, arXiv:2305.11895 (2023), 10.48550/arXiv.2305.11895.
  * Senjean (2019) B. Senjean, Phys. Rev. B 100, 035136 (2019).
  * Sekaran  _et al._ (2022) S. Sekaran, M. Saubanère, and E. Fromager, Computation 10, 45 (2022).
  * Mitra  _et al._ (2023) A. Mitra, M. R. Hermes, and L. Gagliardi, J. Chem. Theory Comput. (2023), 10.1021/acs.jctc.3c00247.
  * Knizia and Chan (2012) G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
  * Wouters  _et al._ (2016) S. Wouters, C. A. Jiménez-Hoyos, Q. Sun, and G. K.-L. Chan, J. Chem. Theory Comput. 12, 2706 (2016).
  * Yalouz  _et al._ (2022) S. Yalouz, S. Sekaran, E. Fromager, and M. Saubanère, J. Chem. Phys. 157, 214112 (2022).
  * Sekaran  _et al._ (2023) S. Sekaran, O. Bindech, and E. Fromager, J. Chem. Phys. 159, 034107 (2023).
  * Ayral  _et al._ (2017) T. Ayral, T.-H. Lee, and G. Kotliar, Phys. Rev. B 96, 235139 (2017).
  * Cancès  _et al._ (2023) E. Cancès, F. Faulstich, A. Kirsch, E. Letournel, and A. Levitt, arXiv:2305.16472 (2023), 10.48550/arXiv.2305.16472.
  * Sun  _et al._ (2020) C. Sun, U. Ray, Z.-H. Cui, M. Stoudenmire, M. Ferrero, and G. K.-L. Chan, Phys. Rev. B 101, 075131 (2020).
  * Ye  _et al._ (2018) H.-Z. Ye, M. Welborn, N. D. Ricke, and T. Van Voorhis, J. Chem. Phys. 149, 194108 (2018), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5053992/16647937/194108_1_online.pdf .
  * Hermes and Gagliardi (2019) M. R. Hermes and L. Gagliardi, J. Chem. Theory Comput. 15, 972 (2019).
  * Sekaran  _et al._ (2021) S. Sekaran, M. Tsuchiizu, M. Saubanère, and E. Fromager, Phys. Rev. B 104, 035121 (2021).
  * Hubbard (1126) J. Hubbard, Proc. R. Soc. London, Ser. A 276, 238 (19631126).
  * Marécat  _et al._ (2023) Q. Marécat, B. Lasorne, E. Fromager, and M. Saubanère, arXiv:2306.07641 (2023), 10.48550/arXiv.2306.07641.
  * Töws and Pastor (2011) W. Töws and G. Pastor, Phys. Rev. B 83, 235101 (2011).
  * Schade and Blöchl (2018) R. Schade and P. E. Blöchl, Phys. Rev. B 97, 245131 (2018).
  * Rotella and Zambettakis (1999) F. Rotella and I. Zambettakis, Appl. Math. Lett. 12, 29 (1999).
  * Caffarel and Krauth (0307) M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 (19940307).
  * Lieb and Wu (1968) E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).
  * Ogata and Shiba (1990) M. Ogata and H. Shiba, Phys. Rev. B 41, 2326 (1990).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * López-Sandoval and Pastor (2002) R. López-Sandoval and G. M. Pastor, Phys. Rev. B 66, 155118 (2002).