# Disentangling Interacting Systems with Fermionic Gaussian Circuits:
Application to the Single Impurity Anderson Model

Ang-Kun Wu  Department of Physics and Astronomy, Center for Materials Theory,
Rutgers University, Piscataway, New Jersey 08854, USA  Center for
Computational Quantum Physics, Flatiron Institute, New York, New York 10010,
USA  Matthew T. Fishman  Center for Computational Quantum Physics, Flatiron
Institute, New York, New York 10010, USA  J. H. Pixley  Department of Physics
and Astronomy, Center for Materials Theory, Rutgers University, Piscataway,
New Jersey 08854, USA  Center for Computational Quantum Physics, Flatiron
Institute, New York, New York 10010, USA  E. M. Stoudenmire  Center for
Computational Quantum Physics, Flatiron Institute, New York, New York 10010,
USA

###### Abstract

Tensor network representations of quantum many-body states provide powerful
tools for strongly correlated systems, tailored to capture local correlations
such as ground states exhibiting entanglement area laws. When applying tensor
network states to interacting fermionic systems, a proper choice of basis or
orbitals can greatly reduce the bond dimension of tensors and speed up
calculations. We introduce such a change of basis with unitary gates obtained
via compressing fermionic Gaussian states into quantum circuits corresponding
to various tensor networks. These circuits can reduce the ground state
entanglement entropy and improve the performance of algorithms such as the
density matrix renormalization group. We study the 1D single impurity Anderson
model to show the power of the method in improving computational efficiency
and interpreting impurity physics. Furthermore, fermionic Gaussian circuits
also show potential for suppressing entanglement during the time evolution of
a low-lying excited state that is used to compute the impurity Green’s
function. Lastly, we consider Gaussian multi-scale entanglement
renormalization ansatz (GMERA) circuits which compress fermionic Gaussian
states hierarchically. The emergent coarse-grained physical models from these
GMERA circuits are studied in terms of their entanglement properties and
suitability for performing time evolution.

##  I Introduction

Tensor network methods, originating from the density matrix renormalization
group (DMRG), provide a powerful set of tools for studying many-body quantum
systems White (1992, 1993); Schollwöck (2005); Hallberg (2006); Schollwöck
(2011). The organizing principle of these methods is entanglement. States with
low entanglement can be approximately factored into a contracted network of
small tensors, leading to fast and precise algorithms. The size of the indices
connecting the tensors, known as the bond dimension, is tied to the amount of
entanglement and directly controls the cost of tensor network calculations.

Entanglement is exquisitely sensitive to the basis in which a state is
represented. Non-interacting fermion systems exhibit some of the most highly
entangled ground states Evenbly and Vidal (2011); Barthel  _et al._ (2006);
Wolf (2006); Vidmar  _et al._ (2017); Bhattacharjee  _et al._ (2021); Bianchi
_et al._ (2021); Kiefer-Emmanouilidis  _et al._ (2020) yet become completely
unentangled when transformed into the momentum basis (assuming periodic
boundary conditions). More generally, low-energy states of local Hamiltonians
typically obey the _area law_ (up to logarithmic corrections), meaning the
entanglement of subsystem of volume $l^{d}$ in dimension $d$ scales as its
area $\sim l^{d-1}$. Though the area law explains why a wide class of low-
energy states have limited entanglement in real space, it does not imply that
real space is the basis with the _lowest_ entanglement.

The insight that another basis besides real space could yield more efficient
tensor network calculations has been put into practice for some time, a
notable example being the work by Krumnow, Eisert, and Legeza on computing
basis transformations Krumnow  _et al._ (2016, 2019, 2021). Their approach
optimizes local two-site unitary gates so as to minimize the entanglement
across a given bipartition of the system, often resulting in significantly
lower entanglement throughout the system. It has been successfully applied to
ground states of interacting one- and two-dimensional systems and to dynamical
systems.

In a separate advance, Fishman and White Fishman and White (2015) showed how
the problem of transforming Gaussian fermionic states into tensor networks
could be solved by finding quantum circuits—fermionic Gaussian circuits—which
provide a basis where the Gaussian state becomes unentangled. Their approach
has many appealing aspects, such as a fast, direct method for computing
circuits that preserve the locality of the Hamiltonian. These many-body
circuits are equivalent to canonical transformations of the single-particle
basis and map straightforwardly into known classes of tensor networks, like
matrix product states (MPS) White (1992, 1993); Rommer and Östlund (1997),
tree tensor networks (TTN) Shi  _et al._ (2006), and multi-scale entanglement
renormalization ansatz (MERA) networks Vidal (2007, 2008); Evenbly and Vidal
(2010).

The relevance of such Gaussian circuits for interacting, non-Gaussian systems
has been unclear. In this work, we synthesize the above ideas, utilizing basis
transformations obtained from non-interacting systems to reduce the
entanglement of interacting systems. We find that Gaussian circuit
transformations can significantly speed up tensor network calculations and
yield physical insights. To demonstrate our approach, we consider the single
impurity Anderson model in one dimension Anderson (1961), whose physical
properties and entanglement structure are known and provide a rich testbed for
our methods.

After demonstrating that Gaussian circuits can uncover a basis in which the
ground state has low entanglement, we show how certain physical properties
also become more transparent in the new basis. The same approach can also
reduce the entanglement of a system undergoing quench dynamics.

Beyond simpler Gaussian circuit transformations based on matrix product
states, we also consider multi-scale circuits known as Gaussian MERA (GMERA).
In addition to reducing entanglement, GMERA circuits provide a physical
interpretation of the transformed degrees of freedom. The basis introduced by
a conventional GMERA is analogous to a “star geometry” bath used in impurity
solvers Bauernfeind  _et al._ (2019); Wolf  _et al._ (2014) while a modified
boundary GMERA Evenbly  _et al._ (2010); Evenbly and Vidal (2014, 2015) is
analogous to a Wilson “chain geometry” bath Wilson (1975); Bulla  _et al._
(2008). The GMERA approach lets us identify an important subset of sites in
the transformed basis, and we show that reducing the bath to consist of these
sites only is sufficient for obtaining dynamical properties.

We begin by reviewing the technology of fermionic Gaussian circuits in Sec.
II. After introducing the impurity models in Sec. III, we explore how Gaussian
MPS (GMPS) circuits can drastically reduce their entanglement in Sec. IV,
resulting in much faster DMRG solutions without ever passing through a more
highly entangled basis. In Sec. V we show how the Gaussian circuit
transformations reveal Kondo screening physics and in Sec. VI we demonstrate
that the disentangling can succeed in systems undergoing dynamics, too.
Finally, Sec. VII explores two types of GMERA circuits, the physical
interpretation of the non-interacting bath following each GMERA
transformation, and the possibility of simulating the impurity model in a
smaller coarse-grained space.

##  II Review: Fermionic Gaussian Tensor Network Circuits

Before turning to interacting systems, we briefly review the technology of
compressing fermionic Gaussian states, following Ref. Fishman and White
(2015), through which _fermionic Gaussian circuits_ Fishman and White (2015);
Surace and Tagliacozzo (2022) (also known as _matchgate circuits_ Valiant
(2002); Jozsa and Miyake (2008)) are defined. Different layouts and patterns
of these circuits produce various Gaussian tensor networks, such as the GMPS
and the GMERA circuits.

Given a non-interacting system with a Hamiltonian

$\displaystyle H_{0}=\sum_{ij}h_{ij}c^{\dagger}_{i}c_{j},$ (1)

and $c_{i}^{\dagger}$ the spinless fermionic creation operator at site $i$,
one can straightforwardly construct its ground state by diagonalizing the
matrix $h=UDU^{\dagger}$ and filling the orbitals (columns of $U$) having the
lowest energies (eigenvalues of $h$). The ground state can be written as

$\displaystyle\ket{\psi_{0}}=\prod_{k=1}^{N_{F}}c_{k}^{\dagger}\ket{\Omega}$
(2)

where $N_{F}$ is the number of fermions and $\ket{\Omega}$ is the vacuum.

The construction of fermionic Gaussian circuits is based on the correlation
matrix

$\displaystyle\Lambda_{ij}=\bra{\psi_{0}}c_{i}^{\dagger}c_{j}\ket{\psi_{0}}=\sum_{k=1}^{N_{F}}U^{*}_{ik}U_{jk},$
(3)

where the goal is to approximately diagonalize $\Lambda$ with a circuit that
also preserves locality. In a basis where $\Lambda$ is diagonal, the state
will have zero entanglement Vidal  _et al._ (2003); Fishman and White (2015).

An important fact about the correlation matrix is that, for a pure and number
conserving fermionic state, its eigenvalues are all exactly $0$ or $1$,
corresponding to empty or filled orbitals. For 1D states obeying the
entanglement area law, such as ground states of gapped and local Hamiltonians,
there exists a contiguous subblock $\mathcal{B}$ of size $B$ that is
independent of system size $N$ in which one of the eigenvalues will be close
to 0 or 1. For ground states of gapless 1D systems, the block size needed to
find an eigenvalue close to 0 or 1 only grows modestly with system size as
$\log{N}$ Fishman and White (2015).

Here we give a brief summary of the deterministic algorithm for diagonalizing
a correlation matrix with local rotations. Local rotations of the form:

$\displaystyle R(\theta)=\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix}$ (4)

are used to transform $2\times 2$ blocks of $\mathcal{B}$, much like in the
$QR$ factorization of a matrix based on Givens rotations. The rotations
$R(\theta)$ are chosen to isolate an _inactive orbital_ or eigenvector of
$\mathcal{B}$ with eigenvalue close to $0$ or $1$ in the upper-left corner of
the block $\mathcal{B}$ of $\Lambda$. The rotation angle $\theta$ is
determined numerically from an eigenvector, for example if the eigenvector is
$v=(v_{1},v_{2})^{T}$, the angle $\theta=-\arctan\frac{v2}{v1}$ rotates the
vector to $v^{\prime}=(1,0)^{T}.$ Higher dimensional eigenvectors can be
rotated iteratively Fishman and White (2015). In the resulting basis, the
first site of the transformed block will correspond to the inactive orbital
which will be approximately disentangled from the rest of the system.

For each subblock $\mathcal{B}$, one can obtain a minimal set of local $B-1$
rotation gates which act together as $R_{\mathcal{B}}$ that isolate an
inactive orbital of $\mathcal{B}$ through $B-1$ rotations. After discarding
the row and column corresponding to the disentangled (inactive) site, one can
repeat the process on a new block, expanding the block as needed to find the
next inactive orbital by including adjacent rows and columns of the
correlation matrix. Once all blocks have had their most inactive orbital
disentangled, the entire correlation matrix will be approximately diagonalized
by a sequence of gates

$R\equiv R_{\mathcal{B}_{1}}R_{\mathcal{B}_{2}}\cdots R_{\mathcal{B}_{N-1}},$
(5)

as described in detail in Fishman and White (2015); Surace and Tagliacozzo
(2022). At the single-particle level, the matrix $R$ defines a transformation
of the $c_{i}$ operators into new operators

$f_{r}=\sum_{i}R_{ri}c_{i},$ (6)

which annihilate fermions in the transformed basis indexed by $r$. These
orbitals generally differ from the single-particle eigenstates, in that they
are local and do not diagonalize the quadratic Hamiltonian.

Once the set of single-particle transformations $R$ have been found, these
transformations can be “lifted” or promoted into many-body quantum gates
forming a quantum circuit. We show an example of an (inverse) GMPS circuit
acting on a many-body MPS in Fig. 1. Specifically, an individual matrix
$R(\theta)$ acting on the $(j,j+1)$ block of the correlation matrix can be
reinterpreted as a unitary gate $\hat{R}_{j,j+1}(\theta)$ acting on sites
$j,j+1$ and having the following form

$\hat{R}(\theta)=\begin{pmatrix}1&0&0&0\\\ 0&\cos\theta&-\sin\theta&0\\\
0&\sin\theta&\cos\theta&0\\\ 0&0&0&1\end{pmatrix}$ (7)

The $4$ by $4$ matrix above acts in the basis
$\\{\ket{\Omega},c^{\dagger}_{j}\ket{\Omega},c^{\dagger}_{j+1}\ket{\Omega},c^{\dagger}_{j}c^{\dagger}_{j+1}\ket{\Omega}\\}$
of the $j$ and $j+1$ sites. These $\hat{R}_{j,j+1}(\theta)$ are the blue
tensors shown in the circuit at the top of Fig. 1. The red ellipse in the
figure highlights the tensors acting within the first block
$\mathcal{B}_{1}=\\{1,2,3,4\\}$, assuming a block size of $B=4$.

Figure 1:  Demonstration of a fermionic Gaussian circuit computed from a
fermionic Gaussian state being applied to an MPS. The MPS ground state in the
original basis is composed of smaller tensors (orange squares) with physical
dimensions (black lines) and connected by bonds (red edges). After applying
the conjugated and reversed fermionic Gaussian gates
$\hat{R}(\theta_{i})^{\dagger}$ (blue rectangles) with fixed block size $B=4$,
the MPS is transformed into a new basis (green squares), hopefully with a
smaller bond dimensions $\chi$. The $B-1$ gates acting on a single block of
size $B=4$, denoted in the text by $\hat{R}_{\mathcal{B}}$, are marked by the
red ellipse. In the case when the original MPS is a fermionic Gaussian state
and the fermionic Gaussian circuit is constructed from its correlation matrix,
the disentangled state will have a bond dimension of $\chi=1$, approximately
disentangled up to controllable errors in the circuit construction.

Fermionic Gaussian circuits obtained as described above are called Gaussian
matrix product state (GMPS) circuits. A slight modification of the procedure
leads to Gaussian MERA (GMERA) circuits discussed in Sec. VII. Once the
circuits are obtained, they can be used to efficiently transform and
disentangle many-body states in the MPS representation as shown in Fig. 1.

The techniques above can be straightforwardly extended to apply to fermions
with spin or to systems not preserving fermion number Thoenniss  _et al._
(2022), though we will continue to use the spinless fermion formalism for
reasons explained in the next section.

##  III Impurity Model System

To explore the ability of Gaussian circuits to disentangle interacting
systems, we will study the 1D single impurity Anderson model (SIAM) Anderson
(1961). This is an ideal test system since the interaction is local and act
only on the impurity site. The properties of this model are well understood,
allowing us to physically interpret properties revealed by the
transformations.

The Hamiltonian for the SIAM is

$H_{\text{A}}=H_{\mathrm{imp}}+H_{\mathrm{hyb}}+H_{\mathrm{bath}}.$ (8)

with the impurity terms

$H_{\mathrm{imp}}=\epsilon_{d}(\hat{n}_{d\uparrow}+\hat{n}_{d\downarrow})+U\hat{n}_{d\uparrow}\hat{n}_{d\downarrow},$
(9)

describing a degenerate impurity orbital occupied by
$\hat{n}_{d\sigma}=d_{\sigma}^{\dagger}d_{\sigma}$ electrons having spin
$\sigma$, energy $\epsilon_{d}$, and an on-site Coulomb repulsion $U$. The
second term hybridizes the impurity with the first bath electron (denoted as
site $2$)

$H_{\text{hyb}}=-V\sum_{\sigma}(d_{\sigma}^{\dagger}c_{2\sigma}+\text{H.c.}),$
(10)

where $V$ is the hybridization strength. The bath electrons are taken to be
the 1D tight-binding model

$H_{\mathrm{bath}}=-t\sum_{r=2,\sigma}^{N-1}(c_{r\sigma}^{\dagger}c_{(r+1)\sigma}+\text{H.c.}),$
(11)

with hopping strength $t$ and total number of sites $N$, including the
impurity. For all cases, we keep the particle-hole symmetry for the impurity
$U=-2\epsilon_{d}$.

Figure 2:  Demonstration of the split-site geometry of the impurity model and
the corresponding GMPS circuit. In the bottom, the spilt-site geometry of the
impurity model $H_{\mathrm{A}}$ with $N=8$ is shown with all up spin degrees
freedom in the left-hand side while all down spins in the right with Eq. 12.
The corresponding GMPS circuit with block size $B=4$ can be obtained by
symmetrizing the circuit from the up spin part (non-interacting limit). To
maintain the locality of the interaction $U$, rotation gates at the impurity
sites are avoided. Blocks of the circuit start from the two ends of the chain,
each of which consist of $B-1$ unitary gates (blue rectangles).

Computationally, when using the spinful $S=1/2$ local Hilbert space
$\ket{0}$,$\ket{\uparrow}$,$\ket{\downarrow}$,$\ket{\uparrow\downarrow}$ for
each site in the original geometry of the model, the DMRG algorithm requires
high bond dimensions $\chi$ (thousands for system sizes of hundreds) to obtain
good precision. Without changing the physics, one can rewrite the model by
separating the spin up and down operators

$\begin{split}a_{-j}\equiv c_{j\uparrow},\quad a_{j}\equiv
c_{j\downarrow},\end{split}$ (12)

where $j=1,2,\cdots,N$ and $j=1$ denotes the impurity operator $a_{-1}\equiv
d_{\uparrow},a_{1}\equiv d_{\downarrow}$. Such a _split-site representation_
leads to more favorable entanglement and MPS bond dimensions. This “unfolded”
representation has been used often for studying impurity systems with MPS
techniques Saberi  _et al._ (2008); Ganahl  _et al._ (2015); Rams and Zwolak
(2020); Kohn and Santoro (2021). The geometry of this split-site
representation and the corresponding GMPS circuit are shown in Fig. 2. The
detailed model in terms of operators $a_{j}$ can be found in Appendix A.

Throughout, we will be especially interested in studying the von Neumann
entanglement entropy, defined as

$S=-\mathrm{Tr}\rho_{A}\ln\rho_{A}$ (13)

where $\rho_{A}=\mathrm{Tr}_{B}\big{[}\ket{\Psi}\bra{\Psi}\big{]}$ is the
reduced density matrix of system $A$ and $\ket{\Psi}$ is the quantum many-body
state of interest.

The MPS formalism allows one to obtain the entropy $S$ efficiently when $A$ is
either a contiguous region of sites of the form $1,2,\ldots,L_{A}$ or when $A$
is a small region of neighboring sites. Computing the entanglement entropy of
what would be a simple bipartitioning of the sites in real space (before
splitting up and down spin degrees of freedom) involves a more complicated ABA
partitioning in the split-site representation. We describe how to compute the
entanglement entropy with MPS in the split-site representation in Appendix B.

##  IV Disentangling Impurity Models with Gaussian MPS Circuits

aleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
bleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]

Figure 3:  Entanglement entropy and MPS bond dimension before and after the
fermionic Gaussian circuit transformation. Properties in the original basis
(dashed curves) are compared with the same properties in the transformed basis
(solid curves) when using the GMPS circuit for systems $N=128,V=t$. Systems
with small interaction $U/t=1$ are shown in blue and strong interaction
$U/t=10$ in red. The upper panel (a) shows the entanglement entropy $S$ as a
function of bipartite (down spins) site distance $j$ away from the impurity
site. Inset: entanglement entropy close to the impurity. The lower panel (b)
shows MPS bond dimensions $\chi(j)$ at each bond. The legend shows the maximum
MPS bond dimension $\chi$ (when using a fixed truncation error cutoff in DMRG
of $10^{-8}$).

We now explore the power of Gaussian MPS circuits for disentangling ground
states. After disentangling an interacting ground state MPS to demonstrate the
idea, we use the circuit to change the basis of the Hamiltonian _before_
computing the ground state to accelerate tensor network calculations.

###  IV.1 Disentangling the ground state

Using DMRG, we obtain the ground state of the single-impurity Anderson model
(SIAM) as an MPS. To disentangle this ground state MPS, we compute a GMPS
circuit from the $U=0$ limit of the SIAM, apply this circuit to the ground
state MPS, and study its entanglement in the transformed basis.

Figure 2 shows this process in more detail, for the split-site geometry we
used. The GMPS transformation starts from the ends of the finite leads and
proceeds toward the impurity. Any gates touching the impurity are omitted.
Because of the left-right (up and down spin) symmetry of the model we use the
same gates from both ends.

Figure 3(a) shows the entanglement entropy before and after the transformation
for both $U/t=1$ and $U/t=10$, showing only the down-spin (right-hand) side of
the system. The entanglement starts out rather high in the original basis,
with an “arc” shape characteristic of gapless systems. After the GMPS
transformation, the entanglement is reduced nearly to zero over most of the
system except for a region near the impurity. Remarkably, the same circuit
based on the $U=0$ limit works well in both cases. The MPS bond dimension
$\chi$ is also significantly reduced by the GMPS transformation (Fig. 3(b)).

Certain physical properties become apparent after the transformation, a
phenomenon we study in more detail in Section V below. In the SIAM, the
impurity is known to be screened by the bath electrons over a length scale
$\xi_{K}$, the Kondo screening length, and should be primarily entangled with
sites within this distance. The Kondo screening length is known to become
larger as $U/t$ increases. Correspondingly, one can see in Fig. 3 that the
region of finite entanglement near the impurity is in fact larger in the
$U/t=10$ system than in the $U/t=1$ system.

###  IV.2 Transforming the Hamiltonian

The speed of a DMRG calculation depends on the entanglement of the state being
calculated. But the approach of the previous section would not yield faster
ground state calculations since the state was computed before the
transformation was ever used. In this section, we exploit the dual nature of
fermionic Gaussian circuits as both many-body quantum circuits and single-
particle basis transformations Fishman and White (2015) to transform the
Hamiltonian before starting any DMRG calculations. It is not obvious this will
work, as the more complicated form of the transformed Hamiltonian could offset
gains from the less entangled state. But the tradeoff turns out to be very
favorable.

To transform the Hamiltonian, we interpret the GMPS circuit as a “matrix
circuit” of sparse Givens rotation matrices (Eq. 4) multiplied together to
form a unitary single-particle transformation $R$. Using the definition
$f_{r}=\sum_{j}R_{rj}c_{j}$ of the fermion operators in the transformed basis,
the quadratic part of the Hamiltonian transforms as

$\displaystyle\sum_{jj^{\prime}}t_{jj^{\prime}}c^{\dagger}_{j}c_{j^{\prime}}$
$\displaystyle=\sum_{rr^{\prime}}\Big{(}\sum_{jj^{\prime}}R_{rj}t_{jj^{\prime}}R^{*}_{r^{\prime}j^{\prime}}\Big{)}f^{\dagger}_{r}f_{r^{\prime}}$
$\displaystyle=\sum_{rr^{\prime}}\tilde{t}_{rr^{\prime}}f^{\dagger}_{r}f_{r^{\prime}}$
(14)

As in the previous section, any matrices (or “matrix gates”) in $R$ that act
on the impurity are dropped from the circuit. This choice keeps the impurity
part of the Hamiltonian $H_{\text{imp}}$ precisely the same after the change
of basis and keeps the interaction $U$ term local.

Running DMRG on the transformed Hamiltonian turns out to be much faster since
the ground state has very low entanglement—see Fig. 4 for timings. The median
sweep time is between 2.5 and 12.5 times _faster_ than in the original basis,
with greater speedups for the larger the system size.

Because the Hamiltonian is less local in the transformed basis, the bond
dimension of its matrix product operator (MPO) representation will necessarily
become larger than in the original basis. To compute the Hamiltonian MPO in
the transformed basis, we use a compression algorithm based on Ref. Chan  _et
al._ , 2016 that is publicly available through the ITensor software Fishman
_et al._ (2022a). The transformed MPO bond dimension grows logarithmically
with system size (lower inset of Fig. 4), reaching $\chi=26$ for the largest
system size of $N=512$. This is in contrast to a bond dimension of $\chi=4$ in
the original basis. An MPO bond dimension of $\chi=26$ is not very large in
practice, and DMRG calculations scale only quadratically with the MPO bond
dimension, versus as the cube of the MPS bond dimension, explaining why the
“disentangled” calculations in the transformed basis can be so much faster
despite a more complicated Hamiltonian.

We conclude this section by mentioning an interesting negative result. We
attempted to reduce the MPO bond dimension further by discarding any quadratic
Hamiltonian terms with an absolute value below a small threshold
$\tilde{t}_{ij}<\epsilon$. But this led to _larger_ MPO bond dimensions
compared to keeping all quadratic terms—see Appendix D for details. In
hindsight, this effect is probably due to the transformed $\tilde{t}_{ij}$
decaying smoothly, whereas truncating them would introduce an artificially
sharp step. Therefore a recommended practice is to keep all Hamiltonian terms
in the transformed basis when compressing the Hamiltonian into an MPO.

Figure 4:  Median DMRG sweep time in the original versus GMPS transformed
basis as a function of system size $N$. Each DMRG calculation involved 20
sweeps and a truncation cutoff of $10^{-6}$; the system used was
$U=t,V=t,\epsilon_{d}=-U/2.$ Top left inset: the speedup, defined as original
sweep time divided by sweep time in the transformed basis, as a function of
system size. Bottom right inset: the maximum bond dimensions $\chi$ of MPO as
a function of system sizes $N$, which asymptotically scale as $\log(N)$.

##  V Kondo screening physics in the transformed basis

Gaussian circuit transformations remove entanglement by localizing “inactive”
single-particle states. What physical properties become more transparent or
easier to obtain in the transformed basis? We will demonstrate that universal
Kondo physics emerges from entanglement in the transformed basis and show that
this basis is practically useful for computing other physical properties to
high accuracy.

Kondo physics emerges from the SIAM in the low-energy limit, where impurity
occupancy becomes pinned to $n_{d}=1$ and a Schrieffer-Wolff transformation
maps the SIAM to an effective Kondo model in which the impurity couples to the
bath only through a spin interaction (in the limit of particle hole symmetry).
The strength of this spin interaction is called the Kondo (exchange) coupling
and is given by

$\displaystyle J_{K}=8V^{2}D/U$ (15)

with $D$ the half bandwidth for the particle-hole symmetric case. A
renormalization group (RG) analysis of the Kondo model (in the metallic and
thermodynamic limit) shows the system always flows toward a strong-coupling
fixed point where the impurity is screened by the bath electrons. The
universal nature of the strong coupling Kondo screened fixed point implies
that correlation functions, as well as the entanglement entropy, are scaling
functions in terms of the Kondo screening length $\xi_{K}$ given by

$\displaystyle\xi_{K}\sim\frac{1}{\sqrt{UV^{2}}}\exp\Big{(}\frac{\pi
U}{8V^{2}}\Big{)}\ .$ (16)

up to third order in $J_{K}$ within the RG Hewson (1993); Boutin and Bauer
(2021); Coleman (2015). In real space, this length scale can be associated
with the size of the Kondo screening cloud.

aleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
bleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
cleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]

Figure 5:  Impurity entanglement entropy and Kondo screening length from
systems $N=400,V=t,\epsilon_{d}=-U/2$ with various $U\in[1,12]$. (a)
$\xi_{K}$, defined as the length where the entanglement entropy drops below
the threshold $S_{\mathrm{imp}}(\xi_{K})=0.1$. (b) The collapse of
entanglement entropy in the mapped GMPS basis, $\tilde{S}(\tilde{r})$, as a
function of scaled distance $\tilde{r}/\xi_{K}$ away from the impurity site,
where $\xi_{K}$ are obtained from (a). Inset: corresponding results for
impurity entanglement entropy $S_{\mathrm{imp}}(r)$ (with subtraction of bath
electron contribution) in the original basis. Note that in the original basis,
we take middle locations across the even-odd result average. (c) Orbitals from
the GMPS transformed basis with $N=100,V=t$. $r$ denotes the distance from the
impurity $r=1$ in real space. Inset: the peaks location $r_{p}$ in real space
versus orbital index $\tilde{r}$ in the GMPS transformed basis.

###  V.1 Kondo Screening Length From Entanglement

A direct measure of the Kondo screening cloud is to study the behavior of the
entanglement entropy at a distance $r$ away from the impurity Wagner  _et al._
(2018); Sørensen  _et al._ (2007a, b); Affleck  _et al._ (2009); Boutin and
Bauer (2021). It has been shown that the impurity entanglement entropy
exhibits universal scaling behavior in terms of $r/\xi_{K}$. The impurity
entanglement entropy in real space is defined as

$\begin{split}S_{\mathrm{imp}}(r,U)&=S(r,U)-S^{(0)}(r),\\\
S(r,U)&=-\mathrm{Tr}[\rho_{r}\ln\rho_{r}],\end{split}$ (17)

where $S^{(0)}(r)$ is entanglement entropy only from the bath ($U=0,V=0$),
which can be obtained from the correlation matrix (Appendix A of Fishman and
White (2015)), and $\rho_{r}$ is the reduced density matrix for the subsystem
of the first $r$ sites around the impurity, computed from the interacting
state at finite $U$. Computing the entanglement entropy in the split-site
representation is explained in more detail in the Appendix B. The impurity
entanglement entropy oscillates between even and odd sites Boutin and Bauer
(2021); Sørensen and Affleck (1996); Wagner  _et al._ (2018), so in what
follows we will always average over even and odd sites numerically. There
exists a scaling function such that

$S_{\mathrm{imp}}(r,U)=g\Big{(}\frac{r}{\xi_{K}(U)}\Big{)},$ (18)

where $g(x)$ is a universal scaling function.

The goal of the bath subtraction in Eq. (17) above is to reveal the location
of the screening cloud, where deep inside the screening cloud
$r\ll\xi_{K}(U)$, $S_{\mathrm{imp}}=\ln 2$, and outside of it
$r\gg\xi_{K}(U)$, $S_{\mathrm{imp}}\rightarrow 0$. In the GMPS transformed
basis, the decay of the entanglement away from the impurity is naturally
achieved and a subtraction of the non-interacting bath entanglement is no
longer needed (see Fig. 3). However, it is not obvious that the entanglement
in the transformed basis will exhibit universal scaling behavior similar to
Eq. 18 since the distance $\tilde{r}$ from the impurity in the transformed
basis is not the same as the real-space distance $r$.

Nevertheless, Fig. 5 shows we can obtain universal information about the Kondo
screening physics directly from the entanglement in the transformed basis. For
the results in panel (a), $\xi_{K}$ is determined in the GMPS transformed
basis by finding the value of $\tilde{r}$ such that
$S(\tilde{r}=\xi_{K})=0.1$, where $0.1$ is an arbitrary small threshold. We
also perform an identical procedure in the original real space basis but using
the quantity $S_{\text{imp}}$ Eq. (17). The values of $\xi_{K}$ determined
this way scale with $U$ as $\sqrt{U}\xi_{K}\sim\exp(CU)$ over a wide range of
$U$ with some constant $C$ (different for each basis) which fits qualitatively
well with the expression in Eq. (16).

In panel (b), we scale the $x$ axis by the Kondo length $\xi_{K}$ and plot the
transformed-basis entanglement against this rescaled axis. There is a clear
scaling collapse for a range of $U$, even better than for the real-space
$S_{\mathrm{imp}}$ results shown in the inset of Fig. 5(b), and without the
difficulty of subtraction and even-odd averaging necessary to analyze the
entanglement collapse in real space Wagner  _et al._ (2018). To understand
better what the $\tilde{r}$ distance means in the transformed basis, Fig. 5(c)
also shows the orbitals $R_{\tilde{r}r}$ in real space $r$ corresponding to a
specific site $\tilde{r}$ in the transformed basis.

###  V.2 Universal Collapse of Spin Correlators

Universal aspects of Kondo physics can also be realized through the
correlations between the impurity spin and the non-interacting bath, given by

$C(r)=\langle{\bf S}_{d}\cdot{\bf S}(r)\rangle=\langle\Psi_{0}|{\bf
S}_{d}\cdot{\bf S}(r)|\Psi_{0}\rangle,$ (19)

where ${\bf
S}_{d}=d^{{\dagger}}_{\alpha}(\bm{\sigma}_{\alpha\beta}/2)d_{\beta}$ and
$\mathbf{S}(r)=\frac{1}{2}\sum_{ss^{\prime}}c^{\dagger}_{r\alpha}\bm{\sigma}_{\alpha\beta}c_{r\beta}$
is the impurity and conduction electron spin density, respectively, and
$\bm{\sigma}$ denotes a vector of Pauli matrices. The spin correlations
exhibit a universal scaling collapse of the form Borda (2007); Holzner  _et
al._ (2009); Barzykin and Affleck (1998)

$C(r,U)=\frac{1}{\xi_{K}}f\left(\frac{r}{\xi_{K}}\right),$ (20)

where $f(x)$ is a universal scaling function.

In Fig. 6(a) we show the spin correlation function $C(r,U=10t)$ computed in
real-space for various DMRG truncation error cutoffs (smaller cutoff means
higher accuracy) or in the transformed basis using a moderate cutoff. We show
results for even sites only, to avoid oscillations. To obtain the transformed-
basis results, we first converge the DMRG calculation, apply the GMPS circuit
to map back to the original basis, then compute measurements. From the figure
we can see that the transformed-basis approach produces results matching the
highest-accuracy calculation in the original basis with a much more moderate
computational effort (smaller bond dimension).

The advantage of performing computations using the transformed basis becomes
even more noticeable when we consider a range of interactions $U$ in Fig. 6(b)
and plot the spin correlations with the axes rescaled by $\xi_{K}$ so as to
obtain a universal scaling collapse, using a value of $\xi_{K}$ determined
analytically from Eq. 16. The main figure shows the excellent collapse
obtained for $1\leq U/t\leq 12$ when using a DMRG cutoff of $10^{-6}$ and the
transformed-basis approach. The inset shows a similar calculation but
performed entirely in the original basis, exhibiting noticeable deviations
from the universal form at larger values of $r$.

aleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
bleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]

Figure 6:  Even site spin-spin correlator $C^{e}(r)$ (Eq. 19) and the scaling
collapse for $N=400,V=t,\epsilon_{d}=-U/2$ with various $U\in[1,12]$.
Empirically, the screening length is found to be $\xi_{K}=e^{0.6U}/\sqrt{U}$.
(a) The correlation function $C(r)$ with $U=10t$ from the DMRG in the GMPS
transformed basis (sold circles) and in the original basis with various DMRG
calculation cutoffs (dashed curves). (b) The correlation function $C(r)$ first
obtained from the DMRG algorithm in the GMPS transformed basis and then
transformed back to the original basis with the DMRG calculation
cutoff$=10^{-6}$. Inset: results in the original basis with the same axes,
parameters and DMRG calculation cutoff.

##  VI Time evolved states transformed by fermionic Gaussian circuits

An interesting question is whether we can use the GMPS circuit we computed
from the non-interacting ground state to help reduce the entanglement of
systems undergoing dynamics.

To understand the dynamics of the impurity, we consider time evolution of
state with a local quench at the impurity site at time $t=0$

$\ket{\Phi(t)}\equiv U(t)\ket{\Phi(0)},\quad\ket{\Phi(0)}\equiv
d_{\downarrow}\ket{\Psi_{0}}$ (21)

where $\ket{\Psi_{0}}$ is the many-body ground state and $U(t)$ is the time-
evolution operator. We focus on the spin-down impurity as the spin-up impurity
results are the same. The dynamics of such a locally quenched state can be
used to compute the impurity Green’s function.

As an initial test of the idea of disentangling dynamics using Gaussian
circuits, we consider the non-interacting, $U=0$ case only in this work. This
choice allows us to obtain precise results to arbitrary times. In contrast,
studying the interacting dynamics in the transformed basis would require using
the time-dependent variational principle (TDVP) algorithm Haegeman  _et al._
(2011, 2016), which can experience convergence issues when the MPS bond
dimension becomes small, such as after a disentangling transformation, and
requires additional subspace expansion techniques. Therefore we leave the
interacting case to a future work. Because entanglement behaves very similarly
between interacting and non-interacting systems, we believe the results below
are typical and do not depend on non-interacting physics in any special way.

aleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
bleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
cleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]

Figure 7:  The entanglement entropy of the quenched non-interacting system.
(a) The entanglement entropy as a function of bipartite sites for several time
snapshot for system size $N=128,V=0.2t$. Solid curves are the entanglement
after applying the GMPS circuit (disentangled) while dashed curves are the
original basis results. (b) The maximum entanglement entropy (across sites) as
a function of time with $N=128.$ Solid curves are results in the GMPS basis
while dashed curves are the original. (c) The time averaged maximum
entanglement entropy (over various $t$ range with $100$ data points) as a
function of system size $N$ with $V=0.2t$.

Since the fermionic Gaussian state is fully characterized by the correlation
matrix, the time evolution operator for the non-interacting system,
$U(t)=e^{-iht}$ , can be considered just a change of single-particle basis on
the correlation matrix. Here $h$ is the matrix of the single-particle
Hamiltonian coefficients from Eq. 1. Thus, at any time $t$, the state can be
fully represented by another correlation matrix

$[\Lambda_{q}(t)]_{ij}=[U(t)\Lambda_{q}(0)U^{\dagger}(t)]_{ij},$ (22)

where
$[\Lambda_{q}(0)]_{ij}=\langle\Phi_{0}|c_{i}^{\dagger}c_{j}|\Phi_{0}\rangle$
and
$[\Lambda_{q}(t)]_{ij}=\langle\Phi(t)|c_{i}^{\dagger}c_{j}|\Phi(t)\rangle$.
Any correlation matrix of a fermionic Gaussian state with low enough
entanglement can also be efficiently compressed into an MPS with a fermionic
Gaussian circuit. The entanglement entropy for the quenched state is denoted
as

$\begin{split}S(r,t)=-\mathrm{Tr}[\rho_{r}(t)\ln\rho_{r}(t)],\end{split}$ (23)

where $\rho_{r}(t)=\mathrm{Tr}_{\bar{r}}\rho,\rho=\ket{\Phi(t)}\bra{\Phi(t)}$
and $\mathrm{Tr}_{\bar{r}}$ denotes trace over complement sites from $r+1$ to
$N$. Note that the entanglement entropy can be computed directly from the
correlation matrix in a straightforward way (Appendix A of Fishman and White
(2015)).

Fig. 7 shows the behavior of the entanglement entropy for the time-evolved
quenched state $\ket{\Phi(t)}$. Panel (a) shows four snapshots of the
entanglement structure. Note that at time $t=0$, the quenched state
$\ket{\Phi(0)}$ has almost no entanglement after applying the GMPS circuit
from the ground state. As the state evolves, the transformed-basis
entanglement grows outward from the impurity and spreads across further sites.
Interestingly, the region where the transformed-basis entanglement is non-zero
corresponds to where the original-basis entanglement deviates from its time
$t=0$ value.

In Fig. 7(b) we focus on the maximum value of entanglement across all sites
and plot its time dependence. In the original basis, the entanglement starts
at a relatively large value around $1.1$ and grows to $1.5$, where it
saturates. (The oscillations are finite-size effects and occur at multiples of
the system size.) In contrast, in the transformed basis the maximum
entanglement starts out very small at early times then becomes nearly constant
($0.69$) at late times, at a value much lower than the original-basis
entanglement. The disentangling power of the circuit is even more prominent
when we look at the system size scaling of the maximum entanglement in panel
(c), which is averaged over time to smooth the finite size oscillations. The
GMPS circuit not only decreases the overall entanglement but also keeps the
growth of entanglement with system size under control. Our results demonstrate
that fermionic Gaussian circuits could be a powerful tool to compute the
unitary time evolution in interacting systems, though we leave this for future
work.

##  VII Entanglement Renormalization of the Bath Electrons

Figure 8:  Demonstration of the GMERA circuit. (a) the GMERA circuit obtained
by diagonalizing odd sites in each iteration with system size $N=16$ and block
size $B=4$. The impurity acts as a local spin in the left. The GMERA circuit
generates effective hybridizations from the impurity to exponential sites,
$r=2^{k}+1,k=0,1,\cdots$, marked by the red dots and red squares. (b) The
transformed hopping Hamiltonian ($U=0$) under the GMERA circuit from the bath
electrons with $N=80,V=t$, where the effective hybridizations (color bar) are
prominent. (c) The part of mutual information matrix $I(r_{1},r_{2})$ of the
ground state with $U=t$ for $N=80,V=t$. The peaks (color bar) in the mutual
information perfectly fit with exponential sites appeared in (a) and (b).

In the sections above, we considered fermionic Gaussian circuits with the
structure of the GMPS, which worked very well in disentangling the interacting
problem with local orbitals. However, the large degeneracy of the correlation
matrix eigenvalues gives us freedom to choose other diagonalization
procedures. In this section, we introduce other types of fermionic Gaussian
circuits with structures inspired by the multi-scale entanglement
renormalization ansatz (MERA) and discuss how the transformed Hamiltonians
have interesting emergent properties in the context of finding ground state
energies and performing time evolution.

###  VII.1 GMERA and Emergent Models

By introducing non-local unitary gates with a hierarchical structure, a class
of fermionic Gaussian circuits implementing the MERA tensor network structure
were proposed in Ref. Evenbly and Vidal, 2010 and a simplified construction
was introduced in Ref. Fishman and White, 2015. This fermionic Gaussian
circuit is called the Gaussian MERA (GMERA) circuit, and its construction can
be conveniently realized as a simple modification of the GMPS circuit
construction where the gates are reorganized into a hierarchical structure.
(The GMERA circuits are closely related to “deep MERA” circuits Kim and
Swingle (2017); Sewell  _et al._ (2022) but are composed of Gaussian
matchgates Eq. (7).) In contrast to the GMPS circuit which constructs a single
layer of gates in one long pass over the system, the GMERA circuit makes a
logarithmic number of passes over the system as shown in Fig. 8(a). In this
way, higher layers of the circuit introduce longer-range basis transformations
in real space, capturing longer-range correlations.

Here we are primarily interested in how the entanglement-based renormalization
group carried out by the GMERA coarse grains the impurity problem and gives
rise to an emergent effective Hamiltonian with fewer degrees of freedom, in
one case resembling a Wilson chain from the numerical renormalization group
approach. Similar to the GMPS construction, we only apply the GMERA circuit to
the bath electrons as shown in Fig. 8(a) in order to avoid spreading the
interaction terms over multiple sites in the transformed basis.

In each layer/iteration, a GMERA circuit starts from the bath sites closest to
the impurity, whose real space location scales exponentially as
$r=2,3,5,\cdots$ (or $r=2^{k}+1$ with $k=0,1,2,\cdots$). The entanglement
renormalization generates effective _bath hybridizations_ directly from the
impurity to these logarithmic number of sites, reflected by the first row and
column in the GMERA-transformed bath in Fig. 8(b).

The most significant hopping terms in the transformed bath form a pattern
resembling the so-called “star geometry” form of an impurity problem Wolf  _et
al._ (2014); Bauernfeind  _et al._ (2019)

$\begin{split}H^{\mathrm{star}}_{\mathrm{eff}}\equiv
H_{\mathrm{imp}}+\sum_{k\in\mathcal{A},\sigma}\bigg{[}\tilde{V}_{k}(d^{\dagger}_{\sigma}c_{k\sigma}+\mathrm{H.c.})+\tilde{\varepsilon}_{k}c_{k\sigma}^{\dagger}c_{k\sigma}\bigg{]},\end{split}$
(24)

where there are no hoppings between bath sites and only hybridizations $V_{k}$
to the impurity site. Here $k$ denotes the GMERA orbitals that run over the
set $\mathcal{A}$ of coarse-grained active orbitals that are strongly
entangled with the impurity. The matrix elements
$\tilde{\varepsilon}_{k},\tilde{V}_{k}$ are the transformed bath couplings to
the impurity, that is, the diagonal and the first row of Fig. 8(b).

Since the entanglement renormalization procedure does not project out diagonal
components of the Hamiltonian, the inactive (diagonal) set of orbitals,
$\mathcal{I}$, needs to be included explicitly as

$\begin{split}H^{\mathrm{GMERA}}_{\mathrm{eff}}=H^{\mathrm{star}}_{\mathrm{eff}}+\sum_{k\in\mathcal{I},\sigma}\tilde{\varepsilon}_{k}c_{k\sigma}^{\dagger}c_{k\sigma},\end{split}$
(25)

for the ground state energy to be preserved.

The decoupling of the complementary sites from the impurity can be further
verified by computing the mutual information matrix, which is defined as

$\begin{split}I(r_{1},r_{2})&=S(\rho_{\\{r_{1}\\}})+S(\rho_{\\{r_{2}\\}})-S(\rho_{\\{r_{1},r_{2}\\}}),\end{split}$
(26)

where $\rho_{\\{\cdots\\}}$ denotes the reduced density matrix for the sites
specified inside the curly brackets. In Fig. 8(c), numerical results including
a non-zero interaction $U$ indicate that the mutual information is strongly
peaked at these exponential sites while the majority of the remaining sites
share negligible mutual information with the impurity 111As a technical but
important detail, in order to completely diagonalize the correlation matrix
for a finite-sized system, the right boundary of the chain is treated by the
GMPS procedure, which creates large non-zero mutual information at the tail as
well as some finite values at limited intermediate sites..

Figure 9:  Demonstration of a boundary GMERA circuit. (a) the boundary GMERA
circuit obtained by diagonalizing even sites in each iteration with system
size $N=16$ and block size $B=4$. The impurity acts as a local spin in the
left. The boundary GMERA circuit generates long-range hoppings among
exponential sites, $r=2^{k},k=0,1,\cdots$, marked by the red dots and red
squares. (b) The transformed bath Hamiltonian under the boundary GMERA circuit
from bath electrons with $N=80,V=t$, where the effective hoppings (color bar)
are prominent. (c) The part of mutual information matrix $I(r_{1},r_{2})$ of
the ground state with $U=t$ for $N=80,V=t$. The peaks in the mutual
information (color bar) perfectly fit with exponential sites appeared in (a)
and (b).

An interesting variation of the GMERA circuit is to omit the part of each
layer closest to the impurity, as shown in Fig. 9(a). We call this pattern the
boundary GMERA (bGMERA) circuit. This is reminiscent of the boundary MERA
introduced in Ref. Evenbly  _et al._ , 2010; Evenbly and Vidal, 2014, 2015, a
modification of the traditional MERA tensor networks in the presence of
boundaries and impurities. Instead of fully diagonalizing the correlation
matrix, this procedure intentionally keeps the first site in each iteration
untreated. These untreated sites also scale exponentially in distance in the
transformed basis, $r=2^{k},$ with $k=1,2,\cdots$ and create long-range
hopping in the transformed bath Hamiltonian as shown in Fig. 9(b). The pattern
from the transformed Hamiltonian provides a coarse-grained version of the
chain-like geometry of the original model (Eq. 8,
$H^{\mathrm{chain}}=H_{\mathrm{A}}$)

$\begin{split}H^{\mathrm{chain}}_{\mathrm{eff}}&=H_{\mathrm{imp}}+\sum_{\sigma}\bigg{[}\tilde{V}(d^{\dagger}_{\sigma}c_{2\sigma}+\mathrm{H.c.})+\tilde{\varepsilon}_{2}c_{2\sigma}^{\dagger}c_{2\sigma}\bigg{]}\\\
&+\sum_{k\in\mathcal{A},\sigma}\bigg{[}\tilde{t}_{k}(c^{\dagger}_{k\sigma}c_{k+1\sigma}+\mathrm{H.c.})+\tilde{\varepsilon}_{k+1}c_{k+1\sigma}^{\dagger}c_{k+1\sigma}\bigg{]},\end{split}$
(27)

where $k$ is the boundary GMERA orbitals and the set $\mathcal{A}$ is the set
of coarse-grained active orbitals from boundary GMERA. The transformed
hybridization $\tilde{V}$ and matrix elements
$\tilde{\varepsilon}_{k},\tilde{t}_{k}$ can be read out from the transformed
quadratic part of the Hamiltonian, i.e., the diagonal and the prominent
hopping at sites $2^{k}$ of Fig. 9(b). Again, for the energy to match the
original, untransformed system, the inactive complementary set of orbitals,
$\mathcal{I}$ have a purely diagonal Hamiltonian which we can add to form the
total Hamiltonian as

$\begin{split}&H^{\mathrm{bGMERA}}_{\mathrm{eff}}=H^{\mathrm{chain}}_{\mathrm{eff}}+\sum_{k\in\mathcal{I},\sigma}\tilde{\varepsilon}_{k}c_{k\sigma}^{\dagger}c_{k\sigma}\
.\end{split}$ (28)

The mutual information, shown in Fig. 9(c), elucidates the pattern of
correlations between the active and inactive orbitals and the impurity. This
is consistent with previous results on boundary MERA where a Wilson-chain-like
structure Okunishi and Nishino (2010) naturally arises from the MERA tensor
network in the presence of impurities Evenbly  _et al._ (2010); Evenbly and
Vidal (2014, 2015). Here we also find the emergent long-range hopping terms
(significant off-diagonal entries in Fig. 9(b)) also form a chain of hoppings
between the active sites, with values decaying with distance from the
impurity. We plan to investigate the connection to the Wilson chain
construction in more depth in a future work.

The correctness of the above two emergent reduced effective Hamiltonians for
approximating the ground-state energy is surprisingly good, as shown in Fig.
10. For a system of size $N=129$, a reduced system of $N_{\mathrm{bGMERA}}=14$
sites can already accurately capture the ground state energy for a wide range
of interactions $U/t\in[1,10]$ with an error less than $0.1\%$ 222Since we
focus on the interaction at the impurity site $1$ and the impurity mutual
information (first row of the matrix) contains all necessary entanglement
information, emergent models are taken from the impurity mutual information
$I(1,r)$ with cutoffs, for both ground state energies and dynamics.. Note that
such great agreement of ground state energies is not obvious purely from the
transformed bath Hamiltonian since some neglected terms of the Hamiltonian
have coefficients with absolute values greater than $0.1$. Thus, entanglement
renormalization of the GMERA circuits yield non-trivial emergent models for
the interacting problem.

Figure 10:  The ground state energy as a function of interaction strength $U$
with $N=129,V=t,\epsilon_{d}=-U/2$. The GMERA energies are obtained by finding
the ground state in the emergent
Hamiltonian,$H^{\mathrm{chain}}_{\mathrm{eff}},H^{\mathrm{star}}_{\mathrm{eff}}$
, plus single-particle energies smaller than the chemical potential, $\mu=0$.
The relevant sites are $N_{\mathrm{GMERA}}=22,N_{\mathrm{bGMERA}}=14$,
respectively, which are selected from the mutual information with cutoff
$I(1,r)\geq 10^{-3}$. Inset: the relative error of the emergent Hamiltonian
benchmarked by the ground state in the original basis.

###  VII.2 Time Evolution under the GMERA Gates

Figure 11:  The dynamics of the impurity Green’s function (Eq. 29) for
emergent Hamiltonians with original system $N=256,V=t$. The blue solid curve
represents the Green’s function with the full size regardless of basis. The
different markers are emergent systems from GMERA (empty markers) and boundary
GMERA (bGMERA, solid markers) transformed basis, which is taken from the
mutual information of the quenched state $\ket{\Phi(t=10)}$. For each case,
the two different emergent system sizes (that are labeled in the legend) are
obtained at two mutual information cutoffs ($10^{-4}$ and $10^{-5}$). Inset:
the deviation of the emergent models away from the full system Green’s
function.

Another application of the GMERA circuits is to time evolve the quenched state
$\ket{\Phi(t)}$ from Eq. 21 and compare the impurity Green’s function in the
coarse-grained models

$\begin{split}G_{\downarrow}(t)&\equiv\langle
d_{\downarrow}^{\dagger}(t)d_{\downarrow}(0)\rangle\\\
&=e^{iE_{0}t/\hbar}\bra{\Psi_{0}}d_{\downarrow}^{\dagger}U(t)d_{\downarrow}\ket{\Psi_{0}}\\\
&=e^{iE_{0}t/\hbar}\langle\Phi(0)|\Phi(t)\rangle,\end{split}$ (29)

where $E_{0}$ is the energy of the ground state. This is an important quantity
both for understanding impurity physics and as an input to the dynamical mean-
field theory (DMFT) method Georges  _et al._ (1996). We still focus on time
evolving non-interacting impurity models as in Section VI. Previously, the
emergent active sites were identified by analyzing the mutual information of
ground states. Though the ground state shares very similar mutual information
peaks with the impurity quenched state at most times 333Note that the mutual
information $I(1,r)=0$ for all $r$ immediately after the quench, but then
becomes quickly stabilized with fixed peak locations except at certain
isolated times., we find it is more accurate to obtain emergent sites from the
quenched state mutual information.

The dynamics of the impurity Green’s function, Eq. 29, is shown in Fig. 11.
Given the same system size $N$, the absolute values of $|G_{\downarrow}(t)|$
are invariant under arbitrary unitary change of basis. In other words,
$|G_{\downarrow}(t)|$ stays the same for the GMPS, the GMERA and the boundary
GMERA circuits without coarse-graining. However, if we choose to only work
with the active sites $\mathcal{A}$ (from Eqs. 24 and 27) based on the mutual
information of $\ket{\Phi(t)}$ at a typical time, e.g. $t=10$, and time evolve
the impurity under the resulting emergent Hamiltonian, the impurity Green’s
function can be well approximated with an exponentially small system size
($N_{\mathrm{GMERA}}=45,N_{\mathrm{bGMERA}}=28$ as compared to the full size
$N=256$) for a wide range of time until finite size effects start to dominate.
Unlike calculating the ground state energies, the inactive sites (sites with
occupancies very close to $0$ or $1$) only provide a phase change in their
eigenstates, which cancel out for the expectation value of the Green’s
function and thus are irrelevant 444However, they are relevant when
determining particle numbers in coarse-grained models. Given the total
particle $N-1$ ($-1$ for impurity annihilation), inactive sites give close to
$0$ or $1$ diagonal matrix elements in the diagonalized correlation matrix,
whose summation $N_{\mathcal{I}}$ gives the corresponding particle number in
the coarse-grained model by $N_{\mathcal{A}}=N-N_{\mathcal{I}}-1$.. The fact
that the Green’s function can be obtained through a logarithmically small
system makes GMERA and boundary GMERA circuits a powerful tools for dynamics,
though their efficacy is yet to be seen for studying the dynamics of
interacting systems.

##  VIII Conclusion and Discussion

Fermionic Gaussian circuits obtained from non-interacting states can
drastically reduce the entanglement of the interacting single-impurity
Anderson model. Circuits known as Gaussian matrix product states (GMPS) are
able to efficiently disentangle many-body ground states across a wide range of
interactions realizing different fixed-point physics of the impurity problem.
Furthermore, transforming the Hamiltonian before performing a DMRG calculation
can yield much faster DMRG runtimes and more accurate properties for a given
truncation threshold.

We studied the Kondo screening cloud under the GMPS circuit transformation,
finding the entanglement remaining after a GMPS transformation harbors
universal information about Kondo screening physics. We also studied the
potential for GMPS transformations to reduce the entanglement of systems
undergoing dynamics, where the entanglement growth of a quenched, non-
interacting system was found to hardly increase with system size after the
transformation, in stark contrast to the behavior in the original real space
basis.

We then tested a variation of fermionic Gaussian circuits inspired by the
multi-scale entanglement renormalization ansatz (MERA) tensor network, called
the fermionic Gaussian MERA (GMERA). Such a renormalization introduces a
hierarchical structure which induces effective long-range interactions. In
this renormalized basis an effective reduced system is identified which is
comprised of a logarithmic number of sites that are highly entangled to the
impurity (or each other, depending on the construction), from which a smaller
effective Hamiltonian emerges and the rest of the sites can be treated as
decoupled free electron modes. Such a bold approximation works exceptionally
well in finding the ground state energies for a wide range of interactions and
can also efficiently simulate the impurity Green’s function dynamics, over a
time range subjected to finite size effect. The essential functionalities for
conveniently working with fermionic Gaussian circuits are now available in the
Julia ITensor ecosystem
555https://github.com/ITensor/ITensors.jl/tree/main/ITensorGaussianMPS.

There are many more directions to explore in applying fermionic Gaussian
circuit transformations to interacting systems. A direct extension of the
current work is to study the effects of using the GMPS circuit as a basis
transformation when time evolving interacting impurity problems with modern
time evolution methods able to handle the complicated transformed Hamiltonian.
Based on the results from non-interacting systems, we believe GMPS circuits
could also keep entanglement of the interacting problem under control. When
using the GMERA circuit, it will be interesting to see how well the emergent
sites from non-interacting quenched states still work in time evolution when
including interactions. In a separate but related vein, recent work has
demonstrated how to efficiently convert the so-called influence matrix in time
evolution of interacting impurity models Thoenniss  _et al._ (2022), which
shows a potential usage of fermionic Gaussian circuits in non-equilibrium
quantum dynamics under non-interacting fermionic reservoirs.

Some more challenging directions involve uniformly interacting systems, such
as in the one-dimensional Hubbard model, where the locality of the interaction
could be crucial. In general, we believe the application of fermionic Gaussian
circuits could push tensor network methods and DMRG toward handling much more
challenging and highly-entangled problems with a proper choice of the
fermionic Gaussian circuit transformation, potentially even in higher
dimensions, i.e., general tree tensor networks, MERA circuits in 2D systems
and Gaussian projected entangled pair states for 2D fermionic problems.

## Acknowledgements

We thank Xiaodong Cao, Chia-Min Chung, Jan von Delft, Örs Legeza, Olivier
Parcollet, and Andreas Weichselbaum for constructive discussions. All
calculations were performed with the ITensor software Fishman  _et al._
(2022a, b) in the Julia language Bezanson  _et al._ (2017). A.-K. W. and
J.H.P. are partially supported by National Science Foundation (NSF) CAREER
Grant No. DMR-1941569 and the Alfred P. Sloan Foundation through a Sloan
Research Fellowship.

##  Appendix A Split-site Representation

In Sec. III we introduced the split spin representation of the Anderson
impurity model, here in this Appendix we describe this in more detail. Without
changing the physics, the impurity model can be rewritten by separating the
spin up and down operators $a_{-j}\equiv c_{j\uparrow},a_{j}\equiv
c_{j\downarrow}$, which has been used often for studying impurity systems with
MPS techniques Saberi  _et al._ (2008); Ganahl  _et al._ (2015); Rams and
Zwolak (2020); Kohn and Santoro (2021).

Specifically, the model can be rewritten in terms of spinless fermion
operators by putting all up spins on the left side of the impurity and all
down spins on the right side. The impurity becomes two spinless fermion sites
in the center of the system. Quantitatively:

$\begin{split}H_{\mathrm{A}}&=H_{\uparrow}+H_{\mathrm{imp}}+H_{\downarrow},\end{split}$
(30)

with

$\begin{split}&H_{\uparrow}=-t\sum_{j=-2}^{-(N-1)}(a_{j}^{\dagger}a_{j-1}+\text{H.c.})-V(a_{-1}^{\dagger}a_{-2}+\text{H.c.})\\\
&H_{\mathrm{imp}}=\epsilon_{d}(\hat{n}_{-1}+\hat{n}_{1})+U\hat{n}_{-1}\hat{n}_{1}\\\
&H_{\downarrow}=-V(a_{1}^{\dagger}a_{2}+\text{H.c.})-t\sum_{j=2}^{N-1}(a_{j}^{\dagger}a_{j+1}+\text{H.c.}).\end{split}$
(31)

Here $a_{j}$ is a spinless fermonic operator. $j<0$ represents original site
$r=-j$ for up spins while $j>0$ is site $r=j$ for down spins, and $-1,1$ sites
are the impurity site for up and down impurity spins (index $0$ is skipped).
$\hat{n}_{j}=a^{\dagger}_{j}a_{j}$ is the occupancy operator. For all cases,
we keep the particle-hole symmetry for the impurity $U=-2\epsilon_{d}$.

As alluded to earlier, the primary advantage of this model geometry is that
the ground state bond dimension $\chi$ is decreased by approximately square
root of the original bond dimension (which becomes exact in the non-
interacting limit $U=0$) for the same MPS accuracy. The trade-off is that
computing operators such as the spin-spin correlators
$C(r)=\langle\mathbf{S}(1)\cdot\mathbf{S}(r)\rangle$ and the bipartite
impurity entanglement entropy $S(r)$ in real space become slightly more
complicated to compute.

##  Appendix B Entanglement entropy of a middle segment in an MPS (ABA
partition)

To obtain the entanglement entropy of a cloud of radius $r$ around the
impurity in real space, we need to divide the split-site geometry of the 1D
chain into three partitions, where the entanglement entropy of interest is
between the middle partition of sites $j\in[-r,r]$ and the rest of the system.
The orthogonality or the canonical form of the MPS is crucial in obtaining the
entanglement of this ABA partition. To begin with, it is helpful to revisit
how the bipartite entanglement is computed. When we decompose a 1D chain into
blocks of $A$ and $B$, where block $A$ goes from $1$ to $l$ while $B$ goes
from $l+1$ to $N$, we can order the MPS in the following way

$\begin{split}\ket{\Psi_{0}}&=\sum_{\sigma_{1},\cdots,\sigma_{2N}}c_{\sigma_{1},\cdots,\sigma_{2N}}\ket{\sigma_{1},\cdots,\sigma_{2N}}\\\
&=\sum_{\sigma_{1},\cdots,\sigma_{2N}}M^{\sigma_{1}}M^{\sigma_{2}}\cdots
M^{\sigma_{2N}}\ket{\sigma_{1},\cdots,\sigma_{2N}}\\\
&=\sum_{a_{l}}\ket{a_{l}}_{A}\ket{a_{l}}_{B},\end{split}$ (32)

where

$\begin{split}\ket{a_{l}}_{A}&=\sum_{\sigma_{1},\cdots,\sigma_{l}}(M^{\sigma_{1}}\cdots
M^{\sigma_{l}})_{1,a_{l}}\ket{\sigma_{1},\cdots,\sigma_{l}},\\\
\ket{a_{l}}_{B}&=\sum_{\sigma_{l+1},\cdots,\sigma_{N}}(M^{\sigma_{l+1}}\cdots
M^{\sigma_{N}})_{a_{l},1}\ket{\sigma_{l+1},\cdots,\sigma_{N}}.\end{split}$
(33)

The boundary tensors $M^{\sigma_{1}},M^{\sigma_{N}}$ give $1$ in the
subscript. Note that this decomposition
$\ket{\Psi_{0}}=\sum_{a_{l}}\ket{a_{l}}_{A}\ket{a_{l}}_{B}$ is not a Schmidt
decomposition, because, in the left-canonical form of $M^{\sigma}$, the states
$\ket{a_{l}}_{A}$ form an orthonormal set for system $A$ while states
$\ket{a_{l}}_{B}$ do not. Thus, to perform the Schmidt decomposition, we apply
the singluar value decomposition (SVD) to the index/bond between
$M^{\sigma_{l}}$ and $M^{\sigma_{l+1}}$, where the diagonal $S$ matrix gives
the singular values.

Figure 12:  Entanglement entropy in different geometries of the single
impurity Anderson model with $N=64,V=1.0t,U=8t,\epsilon_{d}=-U/2$. Solid blue
curve represents bipartite entropy in the electron name space of the
Hamiltonian where the impurity locates in the left. Red dashed curve is the
entanglement entropy of a middle segment of the split-site MPS with a radius
$r$ around the impurity sites. Green curve is the double of the right half of
the bipartite entropy of the split-site MPS, demonstrating that it is
quantitatively similar to the electron entanglement entropy sufficiently far
from the impurity.

For the ABA partition, where the middle part runs from $l_{1}+1$ to $l_{2}$,
we want to use a mixed-canonical MPS automatically facilitated by ITensor,
which allows us to represent the wavefunction amplitudes in Eq. (32) as

$\begin{split}c_{\sigma_{1},\cdots,\sigma_{N}}&=\sum_{\sigma_{1},\cdots,\sigma_{N}}L^{\sigma_{1}}\cdots
L^{\sigma_{l_{1}}}SR^{\sigma_{l_{1}+1}}\dots R^{\sigma_{N}},\\\ \end{split}$
(34)

where $L^{\sigma},R^{\sigma}$ represent left and right normalized,
$\sum_{\sigma}(L^{\sigma})^{\dagger}L^{\sigma}=\sum_{\sigma}R^{\sigma}(R^{\sigma})^{\dagger}=\mathbb{I}$
respectively with the orthogonality center at site $l_{1}+1$. The
decomposition we look for is

$\begin{split}\ket{\Psi_{0}}&=\sum_{b}\ket{a}_{A}\ket{b}_{B}\ket{c}_{C},\\\
\ket{a}_{A}&=\sum_{\sigma_{1},\cdots,\sigma_{l_{1}}}(L^{\sigma_{1}}\cdots
L^{\sigma_{l_{1}}})_{1,a}\ket{\sigma_{1},\cdots,\sigma_{l_{1}}},\\\
\ket{b}_{B}&=\sum_{\sigma_{l_{1}+1},\cdots,\sigma_{l_{2}}}(SR^{\sigma_{l_{1}+1}}\cdots
R^{\sigma_{l_{2}}})_{a,c}\ket{\sigma_{l_{1}+1},\cdots,\sigma_{l_{2}}},\\\
\ket{c}_{C}&=\sum_{\sigma_{l_{2}+1},\cdots,\sigma_{N}}(R^{\sigma_{l_{2}+1}}\cdots
R^{\sigma_{N}})_{c,1}\ket{\sigma_{l_{2}+1},\cdots,\sigma_{N}},\end{split}$
(35)

where for both subsystem $A$ and $C$, the states are orthogonal and states
$\ket{b}_{B}$ will have a size of $a\times c$. However, if we directly extract
the middle segment, we will still work on a huge tensor with physical space
$2^{l_{2}-l_{1}}$ for the SVD while most information within the middle segment
is redundant. Therefore, we can instead contract these physical indexes
$l_{1}+1,\cdots,l_{2}$ but only leave the $2$ bond indexes by multiplying with
complex conjugate of $\ket{b}_{B}$ as the compressed density matrix. In the
end, the SVD will work on a reduced tensor, consisting of $4$ indices
$a,a^{\prime},c,c^{\prime}$, by 2 indices $a,c$. Graphic instructions for this
calculation appear in the Appendix A of Ref. Rogerson  _et al._ , 2022.

When comparing the split-site representation of the Hamiltonian with the
original electron representation as in Fig. 12, it is clear that the bipartite
entanglement entropy in the electron space is the same as the middle segment
entropy of the split-site MPS. Sufficiently far away from the impurity, the
entanglement entropy is well described by twice that of the split-site fermion
MPS.

##  Appendix C Variations of the GMPS circuits

Figure 13:  Disentangling effect of different GMPS circuits for system
$N=128,V=t$ for the small repulsion $U=t$ (upper panel) and the large
repulsion $U=10t$ (lower panel) cases. The solid curves are the same results
as the main text disentangled by the GMPS circuit from the resonant level
model with impurity gates dropped in the end. The orange dashed curves are
results under the whole GMPS circuit without dropping gates. The green dotted
curves are results under the GMPS circuit obtained only from bath electrons.

In the main text, the GMPS circuits are obtained from the $U=0$ limit of the
impurity model (known as the resonant level model, RLM) with impurity gates
dropped so as to maintain locality of the Coulomb repulsion $U$. Such a
treatment benefits the computation without changing the qualitative
disentangling behavior. As shown in Fig. 13, the disentangled results are very
similar between the circuits with and without gates involving the impurity.
Though, the full circuit gives less entanglement at the impurity site in the
small $U$ limit, the entanglements away from the impurity are basically the
same with the circuit that avoids the impurity, which becomes more prominent
in large $U$ limit.

The GMPS circuit obtained from the bath electrons is also considered in Fig.
13. Compared to the resonant level model, disentangling results of the circuit
from the bath has higher entanglement over a large range of sites around the
impurity. Qualitatively, we attribute this to the electron MPS having to
represent the impurity bound state as a superposition over a large number of
bath sites, which is already taken into account exactly in the resonant level
model MPS.

##  Appendix D Hamiltonian MPO with hopping cutoff

Figure 14:  The maximum bond dimensions $\chi$ as a function of physical sizes
$N$ with various hopping cutoffs with $V=t,U=t$. Given a finite cutoff, the
maximum MPO bond $\chi$ dramatically increases at a large system size. Inset:
same results in log-$x$ scale.

In Sec. IV.2 of the main text, the GMPS circuits are used to transform the
Hamiltonian, which generate quadratic number of terms. To constrain the number
of non-zero hopping terms $t_{kl}$ in the transformed non-interacting
Hamiltonian,
$\tilde{H}_{\uparrow}=R^{\dagger}H_{\uparrow}R=\sum_{kl}t_{kl}f_{k}^{\dagger}f_{l}$,
we tried imposing a hopping size cutoff, whose effects on the resulting
compressed MPO are studied here. As shown in Fig. 14, the maximum MPO bond
dimensions does not change with system sizes in the original basis, while it
increases logarithmically in the system size $N$ asymptotically when we use
all non-zero terms in the transformed basis. Perhaps counterintuitively,
attempting to truncate small hopping terms according to a fixed cutoff leads
to MPOs with larger bond dimensions. Thus, small hopping terms have non-
trivial contribution in compressing the Hamiltonian into an MPO. This effect
can possibly be understood as the hoppings having a smooth decaying behavior,
so that truncating them introduces artificial “steps” that can make the MPO
compression algorithm less effective. Therefore in practice we keep all
transformed Hamiltonian terms.

## References

  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * White (1993) S. R. White, Phys. Rev. B 48, 10345 (1993).
  * Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
  * Hallberg (2006) K. A. Hallberg, Adv, Phys. 55, 477 (2006).
  * Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
  * Evenbly and Vidal (2011) G. Evenbly and G. Vidal, J. Stat. Phys. 145, 891 (2011).
  * Barthel  _et al._ (2006) T. Barthel, M.-C. Chung, and U. Schollwöck, Phys. Rev. A 74, 022329 (2006).
  * Wolf (2006) M. M. Wolf, Phys. Rev. Lett. 96, 010404 (2006).
  * Vidmar  _et al._ (2017) L. Vidmar, L. Hackl, E. Bianchi, and M. Rigol, Phys. Rev. Lett. 119, 020601 (2017).
  * Bhattacharjee  _et al._ (2021) B. Bhattacharjee, P. Nandy, and T. Pathak, Phys. Rev. B 104, 214306 (2021).
  * Bianchi  _et al._ (2021) E. Bianchi, L. Hackl, and M. Kieburg, Phys. Rev. B 103, L241118 (2021).
  * Kiefer-Emmanouilidis  _et al._ (2020) M. Kiefer-Emmanouilidis, R. Unanyan, J. Sirker, and M. Fleischhauer, SciPost Phys. 8, 083 (2020).
  * Krumnow  _et al._ (2016) C. Krumnow, L. Veis, O. Legeza, and J. Eisert, Phys. Rev. Lett. 117, 210402 (2016).
  * Krumnow  _et al._ (2019) C. Krumnow, J. Eisert, and Ö. Legeza, arXiv:1904.11999 (2019).
  * Krumnow  _et al._ (2021) C. Krumnow, L. Veis, J. Eisert, and O. Legeza, Phys. Rev. B 104, 075137 (2021).
  * Fishman and White (2015) M. T. Fishman and S. R. White, Phys. Rev. B 92, 075132 (2015).
  * Rommer and Östlund (1997) S. Rommer and S. Östlund, Phys. Rev. B 55, 2164 (1997).
  * Shi  _et al._ (2006) Y.-Y. Shi, L.-M. Duan, and G. Vidal, Phys. Rev. A 74, 022320 (2006).
  * Vidal (2007) G. Vidal, Phys. Rev. Lett. 99, 220405 (2007).
  * Vidal (2008) G. Vidal, Phys. Rev. Lett. 101, 110501 (2008).
  * Evenbly and Vidal (2010) G. Evenbly and G. Vidal, Phys. Rev. B 81, 235102 (2010).
  * Anderson (1961) P. W. Anderson, Phys. Rev. 124, 41 (1961).
  * Bauernfeind  _et al._ (2019) D. Bauernfeind, M. Aichhorn, and H. G. Evertz, arXiv:1906.09077 (2019).
  * Wolf  _et al._ (2014) F. A. Wolf, I. P. McCulloch, and U. Schollwöck, Phys. Rev. B 90, 235131 (2014).
  * Evenbly  _et al._ (2010) G. Evenbly, R. N. C. Pfeifer, V. Picó, S. Iblisdir, L. Tagliacozzo, I. P. McCulloch, and G. Vidal, Phys. Rev. B 82, 161107 (2010).
  * Evenbly and Vidal (2014) G. Evenbly and G. Vidal, J. Stat. Phys. 157, 931 (2014).
  * Evenbly and Vidal (2015) G. Evenbly and G. Vidal, Phys. Rev. B 91, 205119 (2015).
  * Wilson (1975) K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
  * Bulla  _et al._ (2008) R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
  * Surace and Tagliacozzo (2022) J. Surace and L. Tagliacozzo, SciPost Phys. Lect. Notes 54 (2022).
  * Valiant (2002) L. G. Valiant, SIAM J. Comput. 31, 1229 (2002).
  * Jozsa and Miyake (2008) R. Jozsa and A. Miyake, Proc. R. Soc. A 464, 3089 (2008).
  * Vidal  _et al._ (2003) G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. 90, 227902 (2003).
  * Thoenniss  _et al._ (2022) J. Thoenniss, A. Lerose, and D. A. Abanin, arXiv:2205.04995 (2022).
  * Saberi  _et al._ (2008) H. Saberi, A. Weichselbaum, and J. von Delft, Phys. Rev. B 78, 035124 (2008).
  * Ganahl  _et al._ (2015) M. Ganahl, M. Aichhorn, H. G. Evertz, P. Thunström, K. Held, and F. Verstraete, Phys. Rev. B 92, 155132 (2015).
  * Rams and Zwolak (2020) M. M. Rams and M. Zwolak, Phys. Rev. Lett. 124, 137701 (2020).
  * Kohn and Santoro (2021) L. Kohn and G. E. Santoro, Phys. Rev. B 104, 014303 (2021).
  * Chan  _et al._ (2016) G. K.-L. Chan, A. Keselman, N. Nakatani, Z. Li, and S. R. White, J. Chem. Phys. 145, 014102 (2016).
  * Fishman  _et al._ (2022a) M. Fishman, S. R. White, and E. M. Stoudenmire, SciPost Phys. Codebases , 4 (2022a).
  * Hewson (1993) A. C. Hewson,  _The Kondo Problem to Heavy Fermions_ (Cambridge University Press, Cambridge, UK, 1993).
  * Boutin and Bauer (2021) S. Boutin and B. Bauer, Phys. Rev. Research 3, 033188 (2021).
  * Coleman (2015) P. Coleman,  _Introduction to many-body physics_ (Cambridge University Press, 2015).
  * Wagner  _et al._ (2018) C. Wagner, T. Chowdhury, J. H. Pixley, and K. Ingersent, Phys. Rev. Lett. 121, 147602 (2018).
  * Sørensen  _et al._ (2007a) E. S. Sørensen, M.-S. Chang, N. Laflorencie, and I. Affleck, J. Stat. Mech.: Theory Exp. 2007, L01001 (2007a).
  * Sørensen  _et al._ (2007b) E. S. Sørensen, M.-S. Chang, N. Laflorencie, and I. Affleck, J. Stat. Mech.: Theory Exp. 2007, P08003 (2007b).
  * Affleck  _et al._ (2009) I. Affleck, N. Laflorencie, and E. S. Sørensen, J. Phys. A: Math Theor. 42, 504009 (2009).
  * Sørensen and Affleck (1996) E. S. Sørensen and I. Affleck, Phys. Rev. B 53, 9153 (1996).
  * Borda (2007) L. Borda, Phys. Rev. B 75, 041307 (2007).
  * Holzner  _et al._ (2009) A. Holzner, I. P. McCulloch, U. Schollwöck, J. von Delft, and F. Heidrich-Meisner, Phys. Rev. B 80, 205114 (2009).
  * Barzykin and Affleck (1998) V. Barzykin and I. Affleck, Phys. Rev. B 57, 432 (1998).
  * Haegeman  _et al._ (2011) J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
  * Kim and Swingle (2017) I. H. Kim and B. Swingle, arXiv:1711.07500 (2017).
  * Sewell  _et al._ (2022) T. J. Sewell, C. D. White, and B. Swingle, arXiv:2210.16419 (2022).
  * Note (1) As a technical but important detail, in order to completely diagonalize the correlation matrix for a finite-sized system, the right boundary of the chain is treated by the GMPS procedure, which creates large non-zero mutual information at the tail as well as some finite values at limited intermediate sites.
  * Okunishi and Nishino (2010) K. Okunishi and T. Nishino, Phys. Rev. B 82, 144409 (2010).
  * Note (2) Since we focus on the interaction at the impurity site $1$ and the impurity mutual information (first row of the matrix) contains all necessary entanglement information, emergent models are taken from the impurity mutual information $I(1,r)$ with cutoffs, for both ground state energies and dynamics.
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Note (3) Note that the mutual information $I(1,r)=0$ for all $r$ immediately after the quench, but then becomes quickly stabilized with fixed peak locations except at certain isolated times.
  * Note (4) However, they are relevant when determining particle numbers in coarse-grained models. Given the total particle $N-1$ ($-1$ for impurity annihilation), inactive sites give close to $0$ or $1$ diagonal matrix elements in the diagonalized correlation matrix, whose summation $N_{\mathcal{I}}$ gives the corresponding particle number in the coarse-grained model by $N_{\mathcal{A}}=N-N_{\mathcal{I}}-1$.
  * Note (5) https://github.com/ITensor/ITensors.jl/tree/main/ITensorGaussianMPS.
  * Fishman  _et al._ (2022b) M. Fishman, S. R. White, and E. M. Stoudenmire, SciPost Phys. Codebases , 4 (2022b).
  * Bezanson  _et al._ (2017) J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM Rev. 59, 65 (2017).
  * Rogerson  _et al._ (2022) D. Rogerson, F. Pollmann, and A. Roy, J. High Energ. Phys. 2022, 165 (2022).