# Dynamical mean-field theory for Rényi entanglement entropy and mutual
Information in Hubbard Model

Surajit Bera1, Arijit Haldar2,3 and Sumilan Banerjee1 1Centre for Condensed
Matter Theory, Department of Physics, Indian Institute of Science, Bangalore
560012, India
2 S. N. Bose National Centre for Basic Sciences JD Block, Sector-III, Salt
Lake City, Kolkata - 700 106, India
3Department of Physics, University of Toronto, 60 St. George Street, Toronto,
Ontario, M5S 1A7, Canada<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Quantum entanglement, lacking any classical counterpart, provides a
fundamental new route to characterize the quantum nature of many-body states.
In this work, we discuss an implementation of a new path integral method
[Phys. Rev. Res. 2, 033505 (2020)] for fermions to compute entanglement for
extended subsystems in the Hubbard model within dynamical mean field theory
(DMFT) in one and two dimensions. The new path integral formulation measures
entanglement by applying a “kick” to the underlying interacting fermions. We
show that the Rényi entanglement entropy can be extracted efficiently within
the DMFT framework by integrating over the strength of the kick term. Using
this method, we compute the second Rényi entropy as a function of subsystem
size for metallic and Mott insulating phases of the Hubbard model. We explore
the thermal entropy to entanglement crossover in the subsystem Rényi entropy
in the correlated metallic phase. We show that the subsystem-size scaling of
second Rényi entropy is well described by the crossover formula which
interpolates between the volume-law thermal Rényi entropy and the universal
boundary-law Rényi entanglement entropy with logarithmic violation, as
predicted by conformal field theory. We also study the mutual information
across the Mott metal-insulator transition.

##  I Introduction

Entanglement, arguably the strangest aspect of quantum mechanics, signifies
the existence of true non-local quantum correlations. As a result, it has
found enormous applications for characterizing quantum many-body states in
condensed matter and high-energy physics, and as a resource for quantum
computation Nielsen and Chuang (2000). In condensed matter systems,
entanglement can be used to distinguish various kinds of symmetry-broken and
topological states, gapped or gapless phases Laflorencie (2016) etc. For
instance, entanglement provides an unambiguous indicator of topological order
Levin and Wen (2006); Kitaev and Preskill (2006) in quantum ground states.
Entanglement has also emerged as an important measure for distinguishing high-
energy states as well as non-equilibrium dynamics. For example, entanglement
can be used to classify dynamical phases of isolated quantum systems as
ergodic or many-body localized Nandkishore and Huse (2015); Altman and Vosk
(2015); Abanin  _et al._ (2019).

Entanglement of a quantum system is quantified in terms of various measures,
e.g., von Neumann and Rényi entanglement entropies, mutual information and
entanglement negativity Laflorencie (2016); Calabrese and Cardy (2004); Vidal
and Werner (2002). These measures can be calculated by partitioning the
overall system into two subsystems and computing the reduced density matrix of
one of the subsystems by tracing over the other. To this end, the dependence
of entanglement entropy on the size and geometry of the subsystem under
various partitioning of the system are used to classify quantum many-body
states and their non-local entanglement properties. For example, ground-states
of gapped bosonic and fermionic systems in $d$ dimensions follow the so-called
‘area law’ or ‘boundary-law’ for entanglement entropy ($\sim L^{d-1}$) of a
subsystem with length $L$ Laflorencie (2016); Eisert  _et al._ (2010); Casini
and Huerta (2009); Calabrese and Cardy (2004). In contrast, critical states in
one dimension (1d) and fermionic systems with Fermi surface, i.e., standard
metals, in any dimension exhibit a logarithmic violation Laflorencie (2016);
Korepin (2004); Gioev and Klich (2006a); Li  _et al._ (2006); Swingle (2010,
2012a, 2012b); Swingle and Senthil (2013) of the area law, namely the
subsystem entanglement entropy scales as $L^{d-1}\ln{L}$. These
characterizations of the many-body ground states are mainly obtained through
powerful analytical results based on conformal field theory (CFT) methods
Calabrese and Cardy (2004, 2009) and related arguments Swingle (2010, 2012a,
2012b); Swingle and Senthil (2013), as well as numerical results for non-
interacting systems Li  _et al._ (2006); Casini and Huerta (2009). For the
latter, entanglement measures can be computed efficiently using the
correlation matrix of the subsystem Casini and Huerta (2009). However,
numerical computations of entanglement entropy is much more challenging for
interacting systems, typically limited to small systems accessible via exact
diagonalization (ED) or 1d systems through density matrix renormalization
group (DMRG) or heavily numerical and sophisticated quantum Monte Carlo (QMC)
techniques Hastings  _et al._ (2010); Humeniuk and Roscilde (2012); Grover
(2013); Assaad  _et al._ (2014); Broecker and Trebst (2014); Wang and Troyer
(2014); Assaad (2015); D’Emidio (2020).

The above numerical methods have provided many useful insights into
entanglement characteristics of interacting systems. However, there is a lack
of complementary quantum many-body methods, e.g., mean-field theories,
perturbation expansions, and other approximations, for computing entanglement
entropy of interacting systems, unlike those for usual thermodynamic,
spectroscopic and transport properties. The CFT techniques employ a replica
path integral approach Calabrese and Cardy (2004, 2009) where bosonic and
fermionic fields are defined on a non-trivial space-time manifold with
complicated boundary conditions. The latter are often hard to implement within
the standard quantum many-body methodology. To circumvent this difficulty, a
new path integral approach was first developed in ref.Chakraborty and
Sensarma, 2021a for bosons and was subsequently extended to fermionic systems
Haldar  _et al._ (2020); Moitra and Sensarma (2020); Chakraborty and Sensarma
(2021b). In particular, ref.Haldar  _et al._ , 2020 employed this method to
compute Rényi entanglement entropy of Fermi and non-Fermi liquid states of
strongly interacting fermions described by Sachdev-Ye-Kitaev (SYK) and related
models. The new method Haldar  _et al._ (2020) replaces the complicated
boundary conditions in the replica field theory for entanglement Calabrese and
Cardy (2004) by a fermionic self-energy that acts as a non-equilibrium kick.
Using this new path integral formalism, here we develop a dynamical mean-field
theory (DMFT) for Rènyi entanglement entropy in the paradigmatic Hubbard model
Hubbard (1963); Fazekas (1999); Auerbach (1994) of strongly correlated
electrons.

In the last three decades, single-site DMFT approximation and its cluster
extensions Georges  _et al._ (1996); Kotliar  _et al._ (2006), have gained
popularity as a very successful approach to describe Mott metal-insulator
transition and other associated electronic strong correlation phenomena, both
in and out-of-equilibrium Aoki  _et al._ (2014). The integration of DMFT with
first-principle electronic structure methods, like density functional theory
(DFT), has provided a viable route to compute and predict properties of
strongly correlated materials Kotliar  _et al._ (2006). In the DMFT
formulation, the strongly correlated lattice problem is reduced to a problem
of a single impurity or cluster of sites coupled to a self-consistent bath
Georges  _et al._ (1996). The original single-site implementation of DMFT
neglects spatial correlations but captures non-trivial local dynamical quantum
correlations and becomes exact in infinite dimension $d\to\infty$ Georges  _et
al._ (1996). The later cluster extensions of DMFT Hettler  _et al._ (2000);
Kotliar  _et al._ (2001); Maier  _et al._ (2005); Kotliar  _et al._ (2006),
along with state-of-the-art impurity solver like continuous-time quantum Monte
Carlo (CTQMC), incorporates back some of the spatial correlations, and can
even provide a good description of the properties of one-dimensional systems
Bolech  _et al._ (2003); Capone  _et al._ (2004); Kotliar  _et al._ (2006).

In this work, we compute entanglement properties of the correlated metallic
and insulating phases across the Mott metal-insulator transition in the
Hubbard model. We use the cavity method Georges  _et al._ (1996) for the
entanglement path integral of ref.Haldar  _et al._ , 2020 to derive single-
site DMFT self-consistency equations for obtaining the second Rényi entropy
$S_{A}^{(2)}$ of a contiguous subsystem $A$. We show that the Rényi entropy
can be extracted by integrating over the strength of a non-equilibrium ‘kick’
perturbation acting on the imaginary-time evolution. Remarkably, this only
requires the knowledge of onsite single-particle Green’s function for the
subsystem, even in the interacting system, albeit in the presence of the kick.
Due to the entanglement cut(s) and the non-equilibrium kick, both the lattice
and time translation symmetry are broken in the entanglement path integral.
Thus the single-site DMFT is implemented as an inhomogeneous non-equilibrium
DMFT. To this end, we develop an efficient recursive Green’s function method
to solve the DMFT self-consistency equations. Given the computational
complexity of the problem, we only consider the Hubbard model at half filling
and employ a simple DMFT impurity solver, namely the iterative perturbation
theory (IPT). The latter is known to work very well when compared to more
accurate exact diagonalization and QMC impurity solvers for the half-filled
Hubbard model in equilibrium Georges  _et al._ (1996). Our DMFT formulation is
general and can be extended in future to incorporate cluster generalizations
of DMFT Hettler  _et al._ (2000); Kotliar  _et al._ (2001); Maier  _et al._
(2005); Kotliar  _et al._ (2006) and more accurate impurity solvers, e.g.,
CTQMCGull  _et al._ (2011).

Using the inhomogeneous non-equilibrium single-site DMFT for subsystem Rényi
entropy, we compute $S_{A}^{(2)}$ in the Hubbard model as a function of
temperature $T$, interaction $U$, and linear size of the subsystem $N_{A}$ in
1d, and for 2d cylindrical geometry. In particular, we ask how the
entanglement properties of correlated metal described by completely local
self-energy approximation within single-site DMFT compare with those expected
from CFT Calabrese and Cardy (2004); Korepin (2004); Calabrese and Cardy
(2009) and related arguments Swingle (2010, 2012a, 2012b); Swingle and Senthil
(2013). At high temperature, subsystem Rényi entropy $S_{A}^{(2)}$ is
dominated by thermal entropy and, at low temperature, by entanglement Korepin
(2004); Calabrese and Cardy (2004, 2009); Ryu and Takayanagi (2006); Swingle
and Senthil (2013). We indeed find a crossover from thermal to entanglement
behavior in $S_{A}^{(2)}$. Specifically, we show that this crossover in DMFT
metallic state is well described by the known CFT crossover formula in 1d and
its higher dimensional generalization Korepin (2004); Calabrese and Cardy
(2004); Ryu and Takayanagi (2006). We also compare our results for
$S_{A}^{(2)}(N_{A},T)$ with that available from QMC Broecker and Trebst
(2014).

As a measure of entanglement at finite temperature Wolf  _et al._ (2008), we
extract Rényi mutual information between the subsystem $A$ and the rest of the
system from $S_{A}^{(2)}(N_{A})$. We find that the mutual information has a
hysteresis across the first-order Mott metal-insulator transition Georges  _et
al._ (1996) in the $U-T$ plane, culminating at the critical point where the
transition becomes second order. There have been previous studies Udagawa and
Motome (2015); Walsh  _et al._ (2019a, b, 2020) of von-Neumann entropy, mutual
information as well as entanglement spectrum of a single site, or a few sites
within a cluster, via cellular DMFT (CDMFT). Such _local_ entanglement
measures can be computed within the usual equilibrium DMFT formulation.
However, the full subsystem size dependence of the entanglement entropy and
mutual information cannot be obtained through such equilibrium DMFT. On the
contrary, the general method we develop here can be applied for extended
subsystems of arbitrary size and shape and requires an entirely different
implementation through the new path integral technique Haldar  _et al._ (2020)
and non-equilibrium kick term.

The rest of the paper is organized as follows. In Sec.II, we briefly revise
the general path integral formalism of ref.Haldar  _et al._ , 2020 for the
second Rényi entropy and discuss how the Rényi entropy can be extracted by
integrating over the strength of a non-equilibrium kick term. The DMFT
approximation for the entanglement path integral is discussed in Sec.III in
the context of half-filled Hubbard model. The numerical solution of the DMFT
self-consistency equations and some of the benchmarks through the comparisons
with the non-interacting limit and previous QMC simulations are discussed in
Sec.IV. We then discuss our main results for the subsystem size dependence of
Rényi entropy in 1d and 2d Hubbard model and the thermal entropy to
entanglement crossover in Secs.V,VI. In Sec.VII, we discuss the mutual
information across the Mott metal-insulator transition. We summarize our
results and discuss possible future scopes and extensions of our work in
Sec.VIII. The details of the numerical implementations of the DMFT equations,
benchmarks and analysis of the results are given in the Appendices A-G.

##  II The path integral for subsystem Rényi entropy

In this section, we briefly discuss the path integral formalism Haldar  _et
al._ (2020) for subsystem Rényi entropy of fermions in a thermal state. For
concreteness, we consider a system of spin-1/2 fermions on a lattice with $N$
sites in thermal state at a temperature $T$ described by a density matrix
$\rho=e^{-\beta\mathcal{H}}/Z$ with Hamiltonian $\mathcal{H}$. Here
$\beta=1/T$ ($k_{\mathrm{B}}=1$) and $Z$ is the partition function. We obtain
the reduced density matrix for subsystem $A$ with $N_{A}$ site by integrating
the degrees of freedom of the rest of the system $B$, i.e., $\rho_{A}={\rm
Tr_{B}}(\rho)$. Then, the $n$-th Rényi entropy for the subsystem $A$ is

$\displaystyle S^{(n)}_{A}=\frac{1}{1-n}\ln{\rm
Tr_{A}}\big{[}\rho^{n}_{A}\big{]}.$ (1)

In this work, we only consider the second Rényi entropy $S_{A}^{(2)}$ for
simplicity. The method can be generalized easily for higher-order Rényi
entropies. The path integral is constructed based on the following identity
Haldar  _et al._ (2020)

$\displaystyle e^{-S_{A}^{(2)}}\equiv\operatorname{Tr}_{A}[\rho_{A}^{2}]=\int
d^{2}\xi
f(\xi)\operatorname{Tr}[\rho_{A}D(\xi_{1})]\operatorname{Tr}[\rho_{A}D(\xi_{2})],$
(2)

where ‘$\mathrm{Tr}$’ denotes trace over the entire system, and
$D(\xi_{\alpha})=\exp{(\sum_{i\in
A,\sigma}c_{i\sigma}^{\dagger}\xi_{i\sigma\alpha})}\exp{(-\sum_{i\in
A,\sigma}\bar{\xi}_{i\sigma\alpha}c_{i\sigma}^{\dagger})}$ ($\alpha=1,2$) is
normal-ordered _fermionic_ displacement operator Cahill and Glauber (1969)
written in terms of fermionic creation and annihilation operators
$c_{i\sigma}^{\dagger},\ c_{i\sigma}$ at site $i$, in the $A$ subsystem, with
spin $\sigma=\uparrow,\downarrow$. Also, $\bar{\xi}_{i\sigma\alpha},\
\xi_{i\sigma\alpha}$ are static auxiliary Grassmann fields in $A$ with
$d^{2}\xi=\prod_{i\in
A,\sigma\alpha}d\bar{\xi}_{i\sigma\alpha}d\xi_{i\sigma\alpha}$, and

$\displaystyle f(\xi)=2^{N_{A}}e^{-\frac{1}{2}\sum_{i\in
A,\sigma}(\bar{\xi}_{i\sigma 1}\xi_{i\sigma 1}+\bar{\xi}_{i\sigma
2}\xi_{i\sigma 2}-\bar{\xi}_{i\sigma 1}\xi_{i\sigma 2}+\bar{\xi}_{i\sigma
2}\xi_{i\sigma 1})}$ (3)

is a Gaussian factor. The term $\mathrm{Tr}[\rho_{A}D(\xi)]$, called the
characteristic function, can be written in terms of a coherent-state path
integral Haldar  _et al._ (2020),

$\displaystyle e^{-S_{A}^{(2)}}=\frac{1}{Z^{2}}\int d^{2}\xi
f(\xi)\mathcal{D}(\bar{c},c)e^{-(\mathcal{S}+\mathcal{S}_{\xi})}$ (4)

where
$\mathcal{S}=\int_{0}^{\beta}[\sum_{i\sigma\alpha}\bar{c}_{i\sigma\alpha}(\tau)(\partial_{\tau}+\mu)c_{i\sigma\alpha}(\tau)+\mathcal{H}(\bar{c},c)]$,
and $\mathcal{S}_{\xi}=\int_{0}^{\beta}\sum_{i\in
A,\sigma}[\bar{c}_{i\sigma\alpha}(\tau)\delta(\tau-\tau_{0}^{+})\xi_{i\sigma\alpha}-\bar{\xi}_{i\sigma\alpha}\delta(\tau-\tau_{0})c_{i\sigma\alpha}(\tau)]$
is a source term which acts on the $A$ subsystem at imaginary time $\tau_{0}$
and originates from the displacement operator in Eq.(2). The imaginary time
$\tau_{0}$ is arbitrary and can be placed anywhere on the thermal cycle
$0\leq\tau<\beta$. The fermionic fields have the usual anti-periodic boundary
condition $c(\tau+\beta)=-c(\tau)$. For non-interacting systems, it is
straightforward to integrate out Chakraborty and Sensarma (2021a); Haldar  _et
al._ (2020); Moitra and Sensarma (2020) the fermionic fields $\bar{c},c$ and
the auxiliary fields $\bar{\xi},\xi$ to obtain the Rényi entropy. As discussed
in ref.Haldar  _et al._ , 2020, for interacting systems treated within some
non-perturbative approximations, like in large-$N$ models, it is advantageous
to first integrate out the Gaussian auxiliary fields in Eq.(2) and obtain

$\displaystyle
e^{-S_{A}^{(2)}}\equiv\frac{Z_{A}^{(2)}}{Z^{2}}=\frac{1}{Z^{2}}\int\mathcal{D}(\bar{c},c)e^{-(\mathcal{S}+\mathcal{S}_{\mathrm{kick}})},$
(5)

where

$\displaystyle\mathcal{S}_{\mathrm{kick}}=\sum_{i\in
A,\alpha\beta\sigma}\bar{c}_{i\sigma\alpha}(\tau)M_{\alpha\beta}\delta(\tau-\tau_{0}^{+})\delta(\tau^{\prime}-\tau_{0})c_{i\sigma\beta}(\tau^{\prime}),$
(6)

henceforth referred as the _kick_ term, corresponds to an effective time-
dependent self-energy for the fermions at $\tau_{0}$. The matrix

$\displaystyle M$ $\displaystyle=$ $\displaystyle\begin{bmatrix}1&1\\\
-1&1\end{bmatrix}$ (7)

couples the two replicas $\alpha=1,2$. In ref.Haldar  _et al._ , 2020, we have
used the path integral representation of Eq.(5) to evaluate $S_{A}^{(2)}$ of
the SYK model and its several extensions. Below we show that the same
representation can be utilized to formulate a DMFT for Rényi entanglement
entropy in the Hubbard model.

###  II.1 Subsystem Rényi entropy via integration of the kick term

Using Eq.(5), the 2nd Rényi entropy $S_{A}^{(2)}$ can be formally written as

$\displaystyle S_{A}^{(2)}=\beta(\Omega^{(2)}_{A}-2\Omega),$ (8)

where we define $\Omega_{A}^{(2)}\equiv-T\ln Z_{A}^{(2)}$ and $\Omega=-T\ln Z$
is the thermodynamic grand potential. However, direct computation of both
$\Omega_{A}^{(2)}$ and $\Omega$ for interacting systems is difficult in
general and typically requires thermodynamic or coupling constant integration
Haldar  _et al._ (2020); Fetter and Walecka (1971). Here we find a new way to
extract $S_{A}^{(2)}$ by using the kick term in Eq.(5). We consider the
following quantity

$\displaystyle
e^{-S_{A}^{(2)}(\lambda)}=\frac{Z_{A}^{(2)}(\lambda)}{Z^{2}}=\frac{1}{Z^{2}}\int\mathcal{D}(\bar{c},c)e^{-(\mathcal{S}+\lambda\mathcal{S}_{\mathrm{kick}})},$
(9)

which reduces to $S_{A}^{(2)}(\lambda=1)=S_{A}^{(2)}$, the second Rényi
entropy, for $\lambda=1$, and $S_{A}^{(2)}(\lambda=0)=0$. In the above, by
taking the derivative with respect to $\lambda$, we get

$\displaystyle\partial_{\lambda}S_{A}^{(2)}(\lambda)=\frac{\int\mathcal{D}(\bar{c},c)e^{-(\mathcal{S}+\lambda\mathcal{S}_{\rm
kick})}\mathcal{S}_{\rm
kick}}{\int\mathcal{D}(\bar{c},c)e^{-(\mathcal{S}+\lambda\mathcal{S}_{\rm
kick})}}=\langle\mathcal{S}_{\mathrm{kick}}\rangle_{Z_{A}^{(2)}(\lambda)},$
(10)

i.e., the expectation value of the kick term with respect to the effective
partition function $Z_{A}^{(2)}(\lambda)$. Integrating the above equation over
$\lambda$ from 0 to 1, we obtain an expression for the second Rényi entropy

$\displaystyle
S_{A}^{(2)}=\int_{0}^{1}d\lambda\langle\mathcal{S}_{\mathrm{kick}}\rangle_{Z_{A}^{(2)}(\lambda)}.$
(11)

The great advantage of the above expression is that the kick term is quadratic
in Grassmann variables $\bar{c},c$. In the above, we have assumed that no
phase transition occurs as we vary $\lambda$. Such transition as a function of
$\lambda$ in the entanglement action might be present and will be interesting
to study in future. Assuming that there are no transitions with $\lambda$, the
Rényi entropy can be obtained from Eq.(11) using

$\displaystyle\langle\mathcal{S}_{\rm
kick}\rangle_{Z_{A}^{(2)}(\lambda)}=\sum_{i\in
A,\alpha\beta\sigma}M_{\alpha\beta}G_{i\sigma\beta,i\sigma\alpha}(\tau_{0},\tau^{+}_{0}),$
(12)

where the imaginary-time local single-particle Green’s function

$\displaystyle
G_{i\sigma\alpha,i\sigma\beta}(\tau,\tau^{\prime})=-\langle\mathcal{T}_{\tau}c_{i\sigma\alpha}(\tau)\bar{c}_{i\sigma\beta}(\tau^{\prime})\rangle_{Z_{A}^{(2)}(\lambda)}.$
(13)

The Green’s function, however, needs to be evaluated in the presence of the
kick term with variable $\lambda$. In the next section, we show how the
Green’s function can be obtained through the DMFT approximation.

##  III Dynamical mean-field theory for the second Rényi entropy in the
Hubbard model

We consider the nearest-neighbor Hubbard model

$\displaystyle\mathcal{H}=\sum_{ij,\sigma}t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}-\mu\sum_{i}n_{i}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}$
(14)

where $t_{ij}$ is the hopping amplitude between lattice sites ($i=1,\dots,N$)
on 1d and 2d square lattice, $\mu$ is the chemical potential, and $U$ is the
onsite repulsive interaction strength between fermions with opposite spins
$\sigma=\uparrow,\downarrow$. Here,
$n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ and
$n_{i}=\sum_{\sigma}n_{i\sigma}$ are the electronic number operators. We write
down the entanglement action $\mathcal{S}_{\lambda}=\mathcal{S}+\lambda
S_{\mathrm{kick}}$ of Eq.(9) for the Hubbard model as

$\displaystyle\mathcal{S}_{\lambda}$ $\displaystyle=-\int_{0}^{\beta}d\tau
d\tau^{\prime}\sum_{ij,\sigma}\bar{c}_{i\sigma\alpha}(\tau)G^{-1}_{0,i\alpha,j\beta}(\tau,\tau^{\prime})c_{j\sigma\beta}(\tau^{\prime})$
$\displaystyle+\int_{0}^{\beta}d\tau
U\sum_{i\alpha}n_{i\uparrow\alpha}(\tau)n_{i\downarrow\alpha}(\tau),$ (15)

where $G^{-1}_{0,i\alpha,j\beta}(\tau,\tau^{\prime})$ is inverse non-
interacting lattice Green’s function in the presence of entanglement cut(s)
between $A$ and the rest of the systems, i.e.

$\displaystyle G^{-1}_{0,i\alpha,j\beta}(\tau,\tau^{\prime})$
$\displaystyle=-[(\partial_{\tau}-\mu)\delta_{ij}+t_{ij}]\delta(\tau-\tau^{\prime})\delta_{\alpha,\beta}$
$\displaystyle-\lambda\delta_{i\in
A}\delta_{ij}M_{\alpha,\beta}\delta(\tau-\tau_{0}^{+})\delta(\tau^{\prime}-\tau_{0}).$
(16)

As mentioned earlier in Sec.I, the self-energy kick, which only acts on $A$
($\delta_{i\in A}=1$ for $i\in A$ and zero otherwise), breaks both lattice and
time-translation symmetry in this formulation. As a result, we construct a
single-site inhomogeneous non-equilibrium DMFT. We use the cavity method
Georges  _et al._ (1996) to reduce the lattice problem into effective single-
site problems for each of the sites $i=1,\dots,N$, described by the generating
functions

$\displaystyle
Z^{(2)}_{\lambda,i}=\int\mathcal{D}(\bar{c},c)e^{-\mathcal{S}_{\lambda,i}}$
(17)

where the effective action $\mathcal{S}_{\lambda,i}$ is given by

$\displaystyle\mathcal{S}_{\lambda,i}$ $\displaystyle=-\int_{0}^{\beta}d\tau
d\tau^{\prime}\sum_{\sigma\alpha\beta}\bar{c}_{\sigma\alpha}(\tau)\mathcal{G}^{-1}_{i,\alpha\beta}(\tau,\tau^{\prime})c_{\sigma\beta}(\tau^{\prime})$
$\displaystyle+\int_{0}^{\beta}d\tau
U\sum_{\alpha}n_{\uparrow\alpha}(\tau)n_{\downarrow\alpha}(\tau).$ (18)

Here $\mathcal{G}_{i,\alpha\beta}(\tau,\tau^{\prime})$ is the dynamical
_Weiss_ field, such that

$\displaystyle\mathcal{G}^{-1}_{i}(\tau,\tau^{\prime})$
$\displaystyle=-(\partial_{\tau}-\mu)\delta(\tau-\tau^{\prime})\mathcal{I}-\Delta_{i}(\tau,\tau^{\prime})$
$\displaystyle-\lambda\delta_{i\in
A}M\delta(\tau-\tau_{0}^{+})\delta(\tau^{\prime}-\tau_{0}),$ (19)

is a $2\times 2$ matrix in the entanglement replica space, and $\mathcal{I}$
is the identity matrix in the same space. In the above, we have also assumed a
paramagnetic state. Of course, like in equilibrium DMFT Georges  _et al._
(1996), the formulation can be easily extended, to describe entanglement in
ordered states, such as the antiferromagnetic Néel state in the Hubbard model.
The matrix $\Delta_{i}(\tau,\tau^{\prime})$ in Eq.(III) is the hybridization
function which can be expressed in terms of the lattice Green’s function as
discussed below. The impurity Green’s function is related to the Wiess field
via the Dyson equation,

$\displaystyle
G^{-1}_{i}(\tau,\tau^{\prime})=\mathcal{G}^{-1}_{i}(\tau,\tau^{\prime})-\Sigma_{i}(\tau,\tau^{\prime}),$
(20)

where $\Sigma_{i}(\tau,\tau^{\prime})$ is the impurity self-energy. The
Green’s function can be obtained by solving the impurity problem using some
approximate or exact impurity solvers Georges  _et al._ (1996); Aoki  _et al._
(2014), e.g. CTQMC Gull  _et al._ (2011). In this work, for simplicity, and as
a first attempt to compute entanglement via DMFT within the new formalism
Haldar  _et al._ (2020), we use iterative perturbation theory (IPT) Georges
_et al._ (1996) to obtain the self-energy in Eq.(20).

We consider the particle-hole symmetric half-filling case with the chemical
potential $\mu=U/2$. At half-filling, IPT, which essentially retains the self-
energy up to second-order in $U$, is known to work very well Georges  _et al._
(1996) in equilibrium, especially in the metallic phase. As well known, in
this case, IPT coincides with the exact result for both $U\to 0$ and
$U\to\infty$, i.e., the atomic limit, and thus it interpolates well between
the two limits even at intermediate $U$. For the effective non-equilibrium
problem [Eq.(17)], we also use the IPT as an approximate impurity solver. The
IPT self-energy in our case is given by

$\displaystyle\Sigma_{i,\alpha\beta}(\tau,\tau^{\prime})$
$\displaystyle=UG_{{ii},{\alpha\beta}}(\tau,\tau^{+})\delta(\tau^{\prime}-\tau^{+})\delta_{\alpha\beta}$
$\displaystyle-U^{2}\tilde{\mathcal{G}}_{i,\beta\alpha}^{2}(\tau,\tau^{\prime})\tilde{\mathcal{G}}_{i,\alpha\beta}(\tau^{\prime},\tau).$
(21)

Here the first term is Hartree self-energy, and the second one is the second-
order self-energy obtained using Hartree corrected Green’s function

$\displaystyle\tilde{\mathcal{G}}^{-1}_{i,\alpha\beta}(\tau,\tau^{\prime})$
$\displaystyle=$
$\displaystyle\mathcal{G}^{-1}_{i,\alpha\beta}(\tau,\tau^{\prime})-UG_{{ii},{\alpha\beta}}(\tau,\tau^{+})\delta(\tau^{\prime}-\tau^{+})\delta_{\alpha\beta}.$
(22)

Within the single-site DMFT approximation, we assume the self-energy in the
lattice problem to be local and the same as the impurity self-energy. Thus the
lattice Green’s function $G_{i\alpha,j\beta}(\tau,\tau^{\prime})$ is obtained
from the lattice Dyson equation

$\displaystyle\int_{0}^{\beta}$ $\displaystyle
d\tau^{\prime\prime}\sum_{k\gamma}\big{[}G^{-1}_{0,i\alpha,k\gamma}(\tau,\tau^{\prime\prime})-\delta_{ik}\Sigma_{i,\alpha\gamma}(\tau,\tau^{\prime\prime})]G_{k\gamma,j\beta}(\tau^{\prime\prime},\tau^{\prime})$
$\displaystyle=\delta_{ij}\delta_{\alpha\beta}\delta(\tau-\tau^{\prime})$ (23)

The DMFT loop is closed by relating the lattice Green’s function with the
hybridization function $\Delta_{i}(\tau,\tau^{\prime})$. The latter can be
obtained via the cavity method (see Appendix-A.1.1) in terms of the cavity
Green’s function as

$\displaystyle\Delta_{i,\alpha\beta}(\tau,\tau^{\prime})=\sum_{jl}t_{ij}t_{il}G^{(i)}_{j\alpha,l\beta}(\tau,\tau^{\prime}),$
(24)

where the cavity Green’s function, obtained with the $i$-th site removed from
the original lattice, is related to full lattice Green’s function via

$\displaystyle G^{(i)}_{j\alpha,l\beta}(\tau,\tau^{\prime})$
$\displaystyle=G_{j\alpha,l\beta}(\tau,\tau^{\prime})$ $\displaystyle-\int
d\tau_{1}d\tau_{2}\sum_{\gamma\delta}$ $\displaystyle
G_{j\alpha,i\gamma}(\tau,\tau_{1})[G_{i\gamma,i\delta}(\tau_{1},\tau_{2})]^{-1}G_{i\delta,l\beta}(\tau_{2},\tau^{\prime}).$
(25)

The above closes the DMFT self-consistency loop. For our numerical
computations, we further make the large-connectivity Bethe lattice
approximation Georges  _et al._ (1996) for the Cavity Green’s function. As a
result, for the model [Eq.(14)] with only nearest-neighbor hopping, Eq.(24)
becomes

$\displaystyle\Delta_{i,\alpha\beta}(\tau,\tau^{\prime})=t^{2}\sum_{j}^{\prime}G_{j\alpha,j\beta}(\tau,\tau^{\prime}),$
(26)

where $\sum_{j}^{\prime}$ indicates that the summation is over only the
nearest-neighbors of $i$. This approximation makes the computation easier
keeping the essential features of the finite dimensionality through the
lattice Green’s function.

Computationally, the most expensive part of the DMFT loop here is the
inversion of Eq.(III) to obtain the lattice Green’s function
$\boldsymbol{\mathrm{G}}$, a matrix in indices
$(i\alpha\tau,j\beta\tau^{\prime})$. As discussed in the next section and in
Appendix A, we discretize the imaginary time and use a recursive Green’s
function method for large systems to obtain $\boldsymbol{\mathrm{G}}$. We also
benchmark our results by doing direct inverse in Eq.(III) for small systems.

##  IV Numerical solution of DMFT equations to obtain $S_{A}^{(2)}$

We solve the DMFT self-consistency Eqs.(III, 20, III, III, 24, 26) by
discretizing them in imaginary time $\tau$ with discretization step
$\delta\tau$ as discussed in Appendix A in detail. To evaluate $S_{A}^{(2)}$
from Eq.(11), we perform the DMFT for a range of values of $\lambda$ between 0
to 1 with step $\delta\lambda$ as discussed in Appendix B. After obtaining the
DMFT self-consistent Green’s function [Eq.(13)] for different $\lambda$, we
use Eqns.(12) and (11) to compute $S_{A}^{(2)}(\delta\tau)$ for a given
discretization $\delta\tau$. We compute $S_{A}^{(2)}(\delta\tau)$ for
different values of $\delta\tau$ and extrapolate to $\delta\tau\to 0$ limit to
finally obtain $S_{A}^{(2)}$. In Appendix D, we discuss the extrapolation of
$S_{A}^{(2)}(\delta\tau)$ to $\delta\tau\to 0$. As mentioned in the preceding
section, we employ a recursive Green’s function method to obtain the lattice
Green’s function from Eq.(III) (see Appendix A.1). By using the recursive
Green’s function method we can compute $S_{A}^{(2)}$ for reasonably large
systems, $N\leq 100$ in 1d, and $N\leq 20\times 20$ in 2d up to low
temperatures ($T\geq 0.05$, in units of nearest-neighbor hopping amplitude
$t$). The results reported in the main text are for periodic boundary
condition (PBC). We also discuss some results for open boundary condition
(OBC) in Appendix F.

Figure 1: The extrapolated $S_{A}^{(2)}(\delta\tau\to 0)$ (open circles) for
$U=0$ is compared with the $S_{A}^{(2)}$ from correlation matrix calculation
(closed circle+line, ‘corr’) for $N=10$ and three different temperatures
$T=0.2,0.5,1.0$.

###  IV.1 Comparison with the non-interacting limit and QMC

To benchmark the kick integration method of Eq.(11) and the extrapolation
$S_{A}^{(2)}(\delta\tau\to 0)$, we first compare the results for $S_{A}^{(2)}$
for the non-interacting case ($U=0$) with that calculated directly using the
correlation matrix $C_{ij}=\mathrm{Tr}[\rho c_{i}^{\dagger}c_{j}]$ for $i,j\in
A$. The correlation matrix can be easily evaluated using the single-particle
eigenenergies and eigenfunctions of the tight-binding model of Eq.(14) for
$U=0$. The second Rényi entropy is obtained from
$S_{A}^{(2)}=-\mathrm{Tr}\ln[(1-C)^{2}+C^{2}]$ Casini and Huerta (2009);
Haldar  _et al._ (2020). As shown in Fig.1, $S_{A}^{(2)}(\delta\tau\to 0)$ for
different temperatures matches very well with the corresponding $S_{A}^{(2)}$
from correlation matrix calculation for a non-interacting system of size
$N=10$.

Figure 2: DMFT results for $S_{A}^{(2)}$ as a function of subsystem size
$N_{A}$ in 1d Hubard model are compared with QMC data from ref.Broecker and
Trebst, 2014 at three temperatures $T=1.0,0.5,0.2$ for system size $N=32$ and
$U=2$.

We next compare our DMFT results for $S_{A}^{(2)}$ with QMC data taken from
ref.Broecker and Trebst, 2014 for $N=32$ and $U=2.0$. In Fig.(2), the
$S_{A}^{(2)}$ as a function of subsystem size $N_{A}$ for this two method is
shown for high ($T=1.0$) to intermediate temperatures ($T=0.5,0.2$). In high
$T$, the $S_{A}^{(2)}$ the DMFT and QMC results coincide. Even at intermediate
temperatures, the comparison is reasonable given the fact that 1d is the
worst-case scenario for a mean-field approach like single-site DMFT, and that
too, employing an approximate impurity solver like IPT. Nevertheless, it has
been shown Bolech  _et al._ (2003); Capone  _et al._ (2004); Kyung  _et al._
(2006) that cluster extension of DMFT can capture some of the subtle
Luttinger-liquid physics in 1d arising from long-distance correlations. Hence,
cluster extensions of our DMFT approach will be able to provide in the future
a good description of entanglement properties even in 1d.

Given the above benchmarks, in the next sections, we study the subsystem-size
dependence and the entropy to entanglement crossover of $S_{A}^{(2)}$, first
in 1d, and then for the 2d Hubbard model.

##  V $S_{A}^{(2)}$ in 1d Hubbard model

In 1d, single-site DMFT gives rise to a metal-insulator transition at finite
$U$ Bolech  _et al._ (2003); Kyung  _et al._ (2006) at half filling, unlike
the exact Bethe ansatz solution Lieb and Wu (1968). The latter leads to a
metallic state only at $U=0$ and gapped states for any $U>0$. The metallic
state in DMFT is a relic of the infinite dimension inherent in the local self-
energy approximation in single-site DMFT, even though some effects of finite
dimension are fed back through the lattice self-consistency. We first look
into $S_{A}^{(2)}(N_{A},T)$ of this mean-field metallic state in 1d since the
DMFT for entanglement formulated in Sec.III is numerically much easier to
implement in 1d than in higher dimensions.

Figure 3: Result for the second Rényi entropy $S_{A}^{(2)}$ in 1d Hubbard
model with periodic boundary condition. (a) $S_{A}^{(2)}$ as a function of
subsystem size $N_{A}$ is shown for $U=0,0.5,2$ and total system size $N=50$
at temperature $T=0.05$. $S_{A}^{(2)}(N_{A})$ for $U=0,2,3$ and $N=100$ at (b)
$T=0.1$, and (c) $T=0.2$. The $S_{A}^{(2)}$ for $U=0.0$ are calculated using
the correlation matrix and that for $U\neq 0$ using DMFT.

In our DMFT formalism [Sec.III], the subsystem Rényi entropy is obtained from
an imaginary-time path integral. Thus we perform the calculations at finite
temperature with a finite discretization $\delta\tau$. To obtain the ground-
state entanglement, we need to take the $T\to 0$ or $\beta\to\infty$ and
$\delta\tau\to 0$ limit. This is not straightforward since at finite
temperature $S_{A}^{(2)}$ contains both thermal and entanglement entropy
contributions, and we need to _disentangle_ these two contributions as the
$T\to 0$ limit is taken. As discussed below, we find that the thermal to
entanglement crossover in $S_{A}^{(2)}$ for the DMFT metallic state can be
described by the crossover function known from CFT Korepin (2004); Calabrese
and Cardy (2004, 2009); Ryu and Takayanagi (2006).

We show $S_{A}^{(2)}$ in Fig.3 as a function of subsystem size $N_{A}$ for 1d
Hubbard model with periodic boundary condition. In Fig.3(a), the result for
$S_{A}^{(2)}$ vs. $N_{A}$ is shown at low temperature $T=0.05$ for the total
system size $N=50$ and interaction strengths $U=0,0.5,2$. Figs.3(b, c) show
$S_{A}^{(2)}(N_{A})$ at relatively higher temperatures, $T=0.1,0.2$, for
$N=100$ and $U=0,2,3$. Here the $U=0$ results are computed using the
correlation matrix approach discussed in the preceding section, and
$S_{A}^{(2)}$ for non-zero $U$ is obtained through DMFT. At higher
temperatures, $S_{A}^{(2)}$ for $N_{A}\gg 1$ is seen to vary linearly with
$N_{A}$, i.e. a _volume law_ scaling. This indicates the dominance of thermal
entropy at higher temperatures. An _arc-like_ feature emerges at a lower
temperatures. This is the hallmark of entanglement contribution to subsystem
Rényi entropy. Thus, the change of linear to arc-like behavior results from an
entropy to entanglement crossover, as we discuss below.

Gapless systems in 1d, like critical bosonic or spin chains, and gapless
fermionic chains, exhibit the logarithmic violation of the area-law scaling of
entanglement. These systems are usually described by 1+1 D CFT characterized
by some central charge $c$ Cardy (1996); Calabrese and Cardy (2004). The Rényi
entropy $S^{(n)}_{A}$ at $T=0$ for a thermodynamically large system
($N\to\infty$) with one gapless mode is given by the CFT formula for $N_{A}\gg
1$ Holzhey  _et al._ (1994); Vidal  _et al._ (2003); Calabrese and Cardy
(2004)

$\displaystyle
S^{(n)}_{A}=\frac{1}{2}\bigg{(}1+\frac{1}{n}\bigg{)}\left(\frac{c}{6}\right)\log(N_{A})+b^{\prime}.$
(27)

The logarithmic term above is universal with the central charge $c$, and
$b^{\prime}$ is a subleading non-universal constant originating from high-
energy degrees of freedom. For finite $N$, and system with periodic boundary
condition Calabrese and Cardy (2004), the above formula is modified to

$\displaystyle
S^{(n)}_{A}=\frac{1}{2}\bigg{(}1+\frac{1}{n}\bigg{)}\left(\frac{c}{6}\right)\log\bigg{[}\frac{N}{\pi}\sin\big{(}\frac{\pi
N_{A}}{N}\big{)}\bigg{]}+b^{\prime},$ (28)

Similarly, one can obtain the Rényi entropy Korepin (2004); Calabrese and
Cardy (2004) at finite temperature and $N\to\infty$ as

$\displaystyle
S^{(n)}_{A}=\frac{1}{2}\big{(}1+\frac{1}{n}\big{)}\left(\frac{c}{6}\right)\log\bigg{[}\frac{v\beta}{\pi}\sinh\bigg{(}\frac{\pi
N_{A}}{v\beta}\bigg{)}\bigg{]}+b$ (29)

where $v$ is a velocity and $b$ is some non-universal constant. For non-
interacting spinless fermions, $v=v_{\mathrm{F}}$ is the Fermi velocity. In
this case, $S_{A}^{(n)}$ is obtained by adding the contributions of the two
gapless chiral modes, i.e., the left ($L$) and right ($R$) movers, at the two
Fermi points, each with the central charge $c_{L}=c_{R}=c=1$. Eq.(29) reduce
to Eq.(27) for $\beta\to\infty$, i.e., at zero temperature, to give us the
ground-state Rényi entanglement entropies. For $\beta\to 0$, Eq.(29)
reproduces thermal Rényi entropy, $S^{(n)}\simeq(1/2)(1+1/n)(\pi cT/6v)N_{A}$,
e.g. the thermal entropy $S(T)=(\pi cT/6v)N_{A}$ ($n=1$), scaling linearly
with subsystem size. Therefore, the above formula [Eq.(29)] can be viewed as a
crossover formulaSwingle and Senthil (2013) from thermal to entanglement
entropy. The same low-energy degrees of freedom give rise to the universal
part of entanglement and thermal entropy and thus leads to the smooth
crossover. The CFT formulas [Eqs.(27),(28),(29)] are also applicable for
gapless states of interacting fermions in 1d, i.e., for a Luttinger liquid. In
this case, the effect of interaction only enters in the crossover formula
[Eq.(29)] through the renormalized Fermi velocity $v$, whereas the central
charge remains unchanged.

As we discuss in the next section in more detail, for $d>1$, the Fermi liquid
state of interacting fermions at $T=0$ also obeys the universal logarithmic
scaling Swingle (2010) of Eq.(27), which is independent of any Fermi liquid
corrections or Landau parameters. The effect of interaction again only appears
Swingle (2010) in the thermal to entanglement crossover [Eq.(29)] through the
renormalized $v$. The single-site DMFT [Sec.III] with the IPT approximation
Kotliar  _et al._ (2006) is designed to give rise to a Fermi liquid metallic
state even in 1d. Thus we describe our DMFT results for $S_{A}^{(2)}(N_{A},T)$
(with the spin degeneracy taken into account in it) in the paramagnetic
metallic state of 1d Hubbard model using the CFT expression [Eq.(29)] with $c$
replaced by $2c$ for the two chiral modes. We note that Eq.(29) is valid for
thermodynamically large systems ($N\to\infty$). For finite $N$ and $T$, the
analytical expression for $S_{A}^{(n)}(N_{A},T,N)$ is not known Calabrese and
Cardy (2004); Ryu and Takayanagi (2006); Chen  _et al._ (2020) to the best of
our knowledge. Nevertheless, we use Eq.(29) to describe our DMFT data for
relatively large systems like $N=50,100$, as shown in Fig.3, assuming the $N$
to be large enough for neglecting finite $N$ corrections. Alternatively, we
can consider the above crossover function [Eq.(29)] as a fitting function. We
can fit our $S_{A}^{(2)}$ as a function of $N_{A}$ for a fixed $T$ using
Eq.(29) and three fitting parameters $c$, $v$ and $b$ (Appendix E). However,
to reduce the number of fitting parameters, we independently extract the ratio
$(c/v)$ by fitting the low-temperature specific heat (per site) $c_{V}$ from
equilibrium DMFT calculations [see Appendix E.1] with the CFT expression
$c_{V}=(\pi T/3)(c/v)$.

Figure 4:  (a) The second Rényi entropy $S_{A}^{(2)}$ (line + circle) as a
function subsystem size $N_{A}$ for interaction strength $U=2.0$ and system
size $N=50$ at $T=0.05,0.1,0.15$ is shown with the fitted CFT crossover
function (line) of Eq.(29). (b) The extracted central charge $c$ (dashed line
+ circle) as a function of $T$ is shown for different interaction strengths
$U=0,0.5,2$, and compared with the CFT value $c=1$ (solid line). (c) The
central charge $c$ as a function of $U$ is shown at $T=0.1$ for two system
sizes.

With the $(c/v)$ ratio fixed this way, we fit Eq.(29) to our data with two
parameters $c$ and $b$. As shown in Fig.4(a), $S_{A}^{(2)}$ follows the
crossover function quite well. The extracted central charge $c$ is shown as a
function $T$ for fixed $U$, and as a function of $U$ for a fixed $T$ in
Figs.4(b) and (c), respectively. The non-universal fitting parameter $b$ and
the extracted normalized Fermi velocity $v$ are shown in Appendix E.1. We find
that the extracted central charge is close to the free-fermion value $c=1$.
With decreasing temperature the extracted central charge approaches $c=1$ for
the non-interacting ($U=0$) and weakly interacting ($U=0.5$) systems
[Fig.4(b)], and seems to deviate slightly for relatively stronger interaction
$U=2$. However, the deviation might be an artifact of employing the
$N\to\infty$ formula [Eq.(29)] for finite $N$. In Fig.4(c), we see that $N=50$
and $N=100$ give very similar values of $c\simeq 1$ at $T=0.1$ as a function
of $U$, thus assuring convergence with $N$. In summary, we conclude that the
entanglement properties of DMFT Fermi liquid, captured through
$S_{A}^{(2)}(N_{A},T)$, and accessed within the local self-energy
approximation match quite well with that of CFT.

##  VI $S_{A}^{(2)}$ in 2d Hubbard model

Figure 5: The cylindrical subsytem $A$, with $N_{A}$ sites in $x$ direction
and $N_{y}$ sites in the $y$ direction, used for computing entanglement in 2d
Hubbard model, is shown. The periodic boundary condition is applied to the
whole system in both $x$ and $y$ directions.

In this section, we discuss the DMFT results for $S_{A}^{(2)}$ in 2d Hubbard
model for (Fermi liquid) metallic and Mott phases. We consider the system with
periodic boundary conditions in both $x$ and $y$ directions, i.e. a torus
geometry for the system. We subdivide the system along the $x$ axis, i.e., the
entanglement cut is parallel to the $y$ axis like in a cylindrical subsystem
geometry as shown in Fig.5. Due to the entanglement cut we break the
translational symmetry in the $x$-direction while retaining the translation
symmetry in the $y$ direction, along which periodic boundary condition is
applied. As discussed in the Appendix G, the periodic boundary condition along
the $y$-axis allows the wavevector $k_{y}$ to be a good quantum number, and
makes the inversion of the lattice Green’s function of Eq.III easier. In this
case, Eq.III can be decoupled for each $k_{y}$ mode. We discuss below our DMFT
results for $S_{A}^{(2)}$ and its dependence on the subsystem size $N_{A}$,
temperature, and interaction strength.

Most systems in higher dimension ($d>1$) follow boundary law scaling of
entanglement. Even the higher dimensional CFT gives strict boundary law
scaling. However, there are several important exceptions Swingle (2010, 2012a,
2012b); Swingle and Senthil (2013) with logarithmic violation of the boundary
law, e.g., free fermions, Fermi liquids, Weyl fermions in a magnetic field,
non-Fermi liquids with critical Fermi surface, and Bose metals. Here we focus
on Fermi liquid metallic state as captured within DMFT. The underlying reason
behind the violation of the area law is that these systems with Fermi surface
can be effectively described as a collection of patches on the Fermi surface
Swingle (2010); Ding  _et al._ (2012). Each of these Fermi surface patches
acts as a one-dimensional gapless chiral mode described by 1+1 D CFT. These
modes are chiral as they can only propagate with Fermi velocity radially
outward to Fermi surface at very low temperatures. Then, the scaling of
entanglement entropy with $N_{A}$ is simply the one-dimensional logarithmic
scaling multiplied by the number of gapless 1+1 D CFT modes Swingle (2010).
The counting of the number of these mode depends on both the geometry of the
Fermi surface and real space boundary Swingle (2010).

As discussed in the preceding section, for one dimension, we have both right
and left movers mode with central charge $c_{L}=c_{R}=c=1$ and the scaling of
Rényi entanglement entropy is given by Eq.27. For the chiral mode in $d>1$, we
have either $c_{L}=0,~{}c_{R}=c$ or $c_{R}=0,~{}c_{L}=c$, and hence the
contribution (per spin component) of each chiral mode to the Rényi entropy at
$T=0$ is still given by Eq.(27). The counting of the mode is obtained from the
Widom formula Gioev and Klich (2006b); Calabrese  _et al._ (2012); Swingle
(2010, 2012a, 2012b); Leschke  _et al._ (2014), originally developed in the
context of signal processing Widom (1982),

$\displaystyle N_{\rm modes}=\frac{1}{(2\pi)^{d-1}}\frac{1}{2}\int_{\partial
A_{x}}\int_{\partial
A_{k}}dA_{x}dA_{k}|\hat{\boldsymbol{n}}_{x}\cdot\hat{\boldsymbol{n}}_{k}|.$
(30)

The integrals are over the real-space boundary $\partial A_{x}$ of the
subsystem and the Fermi surface $\partial A_{k}$. $\hat{\boldsymbol{n}}_{x}$
and $\hat{\boldsymbol{n}}_{k}$ are the unit normals to the real-space boundary
and the Fermi surface, respectively. Here, the _flux factor_
$|\hat{\boldsymbol{n}}_{x}\cdot\hat{\boldsymbol{n}}_{k}|$ counts the fraction
of modes perpendicular to real-space boundary coming from a Fermi surface
patch at $\boldsymbol{k}$. The Widom formula has been verified numerically Li
_et al._ (2006) for free fermions in $d>1$. For Fermi liquids, where only
forward scattering is relevant, the Widom formula is expected to remain valid
Swingle (2010); Ding  _et al._ (2012) with the same $c$ in Eq.(27) modulo
possible modification of the Fermi surface geometry due to interaction if any.
Going beyond Fermi liquids, the Widom formula may get violated Shao  _et al._
(2015); Hu  _et al._ (2020) or modified Swingle and Senthil (2013), e.g., as
in the case of gapless states of composite fermions in the fractional quantum
Hall regime, quantum spin liquids, and non-Fermi liquids.

For the square lattice Hubbard model [Eq.(14)] that we consider here, the non-
interacting dispersion is $\varepsilon_{k}=-2t(\cos k_{x}+\cos k_{y})$. Thus
we can compute the $N_{\rm modes}$ for the cylindrical subsystem ($N_{A}\times
N_{y}$) as discussed in Appendix H. In this case, $N_{\rm modes}$ is given by
$2N_{y}$, where $N_{y}$ is the number of sites in the $y$ direction along the
entanglement cut. Therefore, taking the spin degeneracy into account, we
expect

$\displaystyle S_{A}^{(2)}/N_{y}=\frac{c}{2}\log(N_{A})+b^{\prime}$ (31)

at $T=0$. Moreover, like in the 1d case [Eq.(29)], we expect entropy to
entanglement crossover at finite temperature for a thermodynamically large
system to be given by

$\displaystyle
S_{A}^{(2)}/N_{y}=\frac{c}{2}\log\bigg{[}\frac{v\beta}{\pi}\sinh\big{(}\frac{\pi
N_{A}}{v\beta}\big{)}\bigg{]}+b$ (32)

where $b$ is again another non-universal constant. Hence the Rényi entropy per
unit length along $y$-direction i.e $S_{A}^{(2)}/N_{y}$ in this cylindrical
subsystem geometry almost have same form as the 1d crossover formula
[Eq.(29)]. As in the 1d case, we do not have any crossover formula that
interpolates between entanglement and thermal entropy for finite $N$ and
finite temperature. For the correlated metallic state obtained in our DMFT
calculations for the 2d Hubbard model, we verify the crossover formula
[Eq.(32)]. Again the effect of interaction only enters in the crossover
formula via the velocity $v$ for a Fermi liquid. We note that a universal
entropy-entanglement crossover formula may be valid more generally, even
beyond Fermi liquids. Similar crossover formulae, constrained by the
temperature dependence of thermal entropy, have been proposed Swingle and
Senthil (2013) to hold even for gapless fermionic systems not described by
Fermi liquid theory or devoid of quasiparticles, e.g., non-Fermi liquids with
critical Fermi surfaces.

Figure 6:  (a) $S_{A}^{(2)}/N_{y}$ for the metallic state of 2d Hubbard model
is shown for a $N=20\times 20$ lattice as a function of the subsystem length
along $x$-axis $N_{A}$ for temperatures $T=0.1,0.15,0.2$ at fixed interaction
strength $U=4$. The results (circle) fit well with the crossover function
(line) of Eq.(32). (b) $S_{A}^{(2)}/N_{y}$ (fill circle, open circle,
triangle) as a function of $N_{A}$ is shown for different interaction
strengths $U=0,4,6$ at fixed $T=0.1$. The corresponding fits (line) to the
crossover function are also shown.

We compute the second Rényi entropy $S_{A}^{(2)}$ as a function of subsystem
size $N_{A}N_{y}$ for a total system size $N_{x}N_{y}$. Here $N_{A}$ is the
number of lattice sites in the subsystem $A$ in the $x$ direction, and $N_{x}$
and $N_{y}$ are the total number of sites in the $x$ and $y$ directions,
respectively. We vary the subsystem size by varying $N_{A}$ while keeping
$N_{y}$ fixed. In Fig.6(a), $S_{A}^{(2)}/N_{y}$, computed from DMFT, is shown
as a function of $N_{A}$ at low temperatures, $T=0.1,0.15,0.2$, for $U=4$ and
system size $N=20\times 20$. As shown in the figure, we fit the DMFT results
with the crossover formula of Eq.(32). We find the formula to fit very well to
our result as shown in the [Fig.(6)(a)]. Similarly, in Fig.(6)(b), the
$S_{A}^{(2)}/N_{y}$ is shown for different interaction strengths $U=0,4,6$ at
fixed $T=0.1$ with the corresponding fits to the crossover formula [Eq.(32)]
for a $20\times 20$ system.

Figure 7:  (a) The extracted central charge $c$ as a function of $T$ is shown
for different interactions $U=0,2,4,6$ and system size $N=20\times 20$, and
compared with the CFT value $c=1$. (b) The variation of the $c$ as a function
of $U$ is shown at fixed $T=0.1$ for different system sizes $N$, as indicated
in the legends.

The $c$ extracted from the above fittings is shown in Fig.7. More details are
given in Appendix H. Fig.7(a) shows the extracted $c$ as a function of $T$ for
$U=0,2,4,6$, and compares with the expected CFT value $c=1$. We see that the
calculated $c\simeq 0.9<1$, and $c$ does not vary much with $T$. We find that
$c\simeq 0.9<1$ even for $U=0$, i.e., $c$ deviates from 1 by more or less the
same amount even for the non-interacting case for the system sizes accessed.
Thus, presumably, the deviation from the CFT value stems from the application
of the crossover formula [Eq.(32)] for the thermodynamic limit to the finite
systems. In Fig.7(b), where we plot $c$ as a function of $U$ for different $N$
and $T=0.1$, we observe that $c$ decreases slightly for larger $U$. $c$ tends
to increase very slowly with increasing $N$, implying that $c$ might approach
the expected value of $1$ for larger systems. However, in the absence of an
analytical crossover function for finite $N$ and $T$, and for the accessible
system sizes in our calculations, it is hard to extrapolate $c$ to the
$N\to\infty$. Nevertheless, we can conclude that modulo finite-size effects,
our DMFT results for $S_{A}^{(2)}(N_{A},T,N)$ for the metallic state of the 2d
Hubbard model at half filling are consistent with the Widom formula and the
expectations from CFT Swingle (2010).

We also compute the Rényi entropy deep inside the insulating phase at low
temperature, as shown in Fig.8 for $U=14$. We find that $S_{A}^{(2)}/N_{y}$ is
linear in $N_{A}$ with a slope $\gtrsim\ln 2$, with a very weak dependence on
$T$. This is expected for $T\ll U$, due to non-zero $T\to 0$ residual
entropyGeorges  _et al._ (1996), arising from spin degeneracy of the
paramagnetic Mott insulating solution in the single-site DMFT, .

Figure 8: The $S_{A}^{(2)}/N_{y}$ as a function of $N_{A}$ is shown for
several $T$’s at fixed $U=14$. This is computed from system size $20\times
20$.

##  VII Mutual information across Mott transition in 2d Hubbard model

Figure 9: (a) Double occupancy $D$, computed within single-site DMFT, as a
function of $U$ is shown at $T=0.14$ for 2d Hubbard model. The legend ‘met. to
ins.’ implies increasing $U$ from the metal to insulator and ‘ins. to met.’
decreasing $U$ the other way round. (b) The phase diagram for 2d half-filled
Hubbard model constructed based on the double occupancy $D$ is shown.
$(U_{c},T_{c})$ represents Mott critical point where the first-order line, or
the coexistence region, ends. In the co-existence region, both metal and
insulator solutions exist.

Here we discuss the second Rényi mutual information as an entanglement and
information theoretic measure at finite temperature in the temperature vs.
interaction ($T-U$) phase diagram of the Hubbard model. The second Rényi
mutual information,

$\displaystyle I(A,B)$ $\displaystyle=S^{(2)}_{A}+S^{(2)}_{B}-S^{(2)}_{A\cup
B}$ $\displaystyle=-(\ln Z_{A}^{(2)}+\ln Z_{B}^{(2)}-\ln Z_{A\cup
B}^{(2)}+2\ln Z)$ (33)

is obtained from the combination of the Rényi entropies of a subsystem $A$,
its complement $B$, and the whole system $A\cup B$. While the entanglement
entropy can characterize pure states, e.g., ground state and quantum phase
transitions between ground states at $T=0$, the mutual information is a better
information theoretic measure for finite-temperature phases and phase
transitions Melko  _et al._ (2010); Singh  _et al._ (2011); Iaconis  _et al._
(2013); Stéphan  _et al._ (2014). The mutual information is dominated by
entanglement contribution when classical correlations are short-ranged, e.g.,
at low temperature away from a finite-temperature critical point Wolf  _et
al._ (2008). Moreover, different parts of mutual information can exhibit
critical properties Melko  _et al._ (2010); Singh  _et al._ (2011); Iaconis
_et al._ (2013); Stéphan  _et al._ (2014) at temperatures related to the
finite-temperature critical point, e.g., at critical temperature $T_{c}$, and
at $2T_{c}$ due to critical behaviours of $Z_{A}^{(2)},~{}Z_{B}^{(2)}$ in
Eq.(33) from the edges and corners of the subsystem $A,~{}B$.

For a pure-state density matrix, $I(A,B)=2S_{A}^{(2)}$. The $S_{A}^{(2)}$ for
thermal density matrix at finite temperature contains both entanglement and
thermal entropy contributions. However the (Rényi) mutual information, by
construction, naturally excludes the volume-law thermal entropy of the
subsystem and its complement. Thus mutual information follows in general an
area-law scaling with subsystem size and captures both quantum (entanglement)
and classical correlations between the subsystems. The study of Mott
transition through the lens of mutual information is less explored in
literature. In refs.Walsh  _et al._ , 2019a, b, 2020, the authors have studied
the Mott transition in the 2d Hubbard model using equilibrium CDMFT through
the mutual information of a single site and the rest of the system. They
detect first-order phase transitions and supercritical regime for $T>T_{c}$.
The calculation of the singe-site mutual information only requires the
knowledge of occupation and double occupancy, that can be computed within
equilibrium DMFT. The subsystem size scaling of mutual information cannot be
captured within such equilibrium DMFT. Within our new path integral approach
we can easily study the subsystem size scaling of mutual information across
Mott transition, as we discuss below.

Figure 10: The bipartite mutual information for $N_{A}=N_{B}$ per site of the
subsystem along $y$ axis, i.e., $I(A,B)/N_{y}$ as a function of interaction
$U$ at $T=0.16<T_{c}$ is shown for a $N=16\times 16$ lattice. ‘Met. to ins’
refers to sweeping $U$ from the metal to insulator phase and ‘ins. to met’ for
the other way round.

To characterize the metal, insulator and the metal-insulator phase transition
at finite temperature, we first briefly discuss the $T-U$ phase diagram of the
half-filled 2d Hubbard model within large connectivity Bethe lattice
approximation in DMFTGeorges  _et al._ (1996). We draw the $T-U$ phase diagram
by monitoring the equal time correlation function $D=\langle
n_{\uparrow}n_{\downarrow}\rangle$ corresponding to the double occupancy of a
site within equilibrium DMFT (Appendix E.1) as a function of $U$ for different
temperatures. A representative plot for double occupancy $D$ vs. $U$ at
$T=0.14$ is shown in Fig.9(a). The hysteresis behaviour is due to the
coexistence of both metal and insulator solution across the first-order Mott
metal-to-insulator transition. The area of the hysteresis loop shrinks to zero
as the first order transition ends at finite temperature Mott critical point
($U_{c},T_{c}$). By monitoring $D(U)$, we obtain the $T-U$ phase diagram, as
shown in Fig.9(b). The critical point is at $(U_{c},T_{c})\approx(9.9,0.18)$.

In Fig.10, the mutual information per site along the $y$ durection for the
equal bi-partition $N_{A}=N_{B}$, i.e., $I(A,B)/N_{y}$ is shown at
$T=0.16<T_{c}$ for system size $N=16\times 16$. As mentioned earlier, for the
first-order Mott transitions, two co-existing solutions, metal and insulator,
appear in the phase diagram as indicated in Fig.9(b). In equilibrium DMFT, for
$T<T_{c}$, starting from the insulating solution at large $U$ and on
decreasing the interaction slowly, a sudden jump to the metallic solution
occurs at $U_{c1}(T)$, i.e., at the limit of metastability of the insulating
phase, as usual in a first-order transition. Similarly, sweeping $U$ from the
metallic side leads to a jump to the insulating solution at $U_{c2}(T)$. For
$T=0.16$, $U_{c1}$ and $U_{c2}$ computed from double occupancy within the
equilibrium calculation are shown in Fig.10 as dotted vertical lines. The DMFT
for Rényi entropy also leads to similar hysteresis behaviour in the mutual
information, as shown in Fig.10, where the calculation of $I(A,B)$ vs. $U$ is
done in steps of $\delta U=0.2$. Following the behaviour of double occupancy,
the bipartite mutual information $I(A,B)/N_{y}$ also jumps across
$U_{c1},U_{c2}$. Thus, the mutual information between two extended subsystems
can detect the first order nature of phase transitions, like the single-site
mutual information Walsh  _et al._ (2019a, b, 2020). Finite but weak
correlations, indicated by the non-zero mutual information, persist even in
the insulating phase for $U\gtrsim U_{c2}$, as can be seen in Fig.10. We
expect these correlations to approach zero for large interaction strengths
$U\gg W$, where $W$ is the non-interacting bandwidth.

The calculation of mutual information near the critical point $(U_{c},T_{c})$
through the extrapolation $\delta\tau\to 0$ becomes challenging due to
multiple solutions as well as very close numerical values of
$S^{(2)}(\delta\tau)$ for different $\delta\tau$’s. For this reason, we
present the data for a fixed $\delta\tau=0.029$, without any extrapolation, in
this section.

Figure 11: The mutual information $I(A,B)/N_{y}$ for $N=16\times 16$ as a
function of the length of subsystem $N_{A}$ along $x$-direction is shown here
across the Mott metal-insulator transition for different interaction strength
$U$ ($U=8.0,9.0,9.6,10.0,11.0$) and temperatures. $U$ gradually increases from
top to bottom, as indicated by the arrow, in each plot for a fixed
temperature.

The subsystem size scaling of mutual information is expected to be of the
formMelko  _et al._ (2010); Singh  _et al._ (2011)
$I(A,B)=I(N_{A},N_{y},\beta)\simeq
a(\beta,N_{A})N_{y}+d(\beta,N_{A})+\mathcal{O}(N^{-1}_{y})$, with coefficients
$a$ and $b$ weakly dependent on $N_{A}$. In our subsystem geometry, the
dominant contribution to $I(A,B)$ comes from the interface of the two
subsystem $A$ and $B$, and leads to the leading area law ($\propto N_{y}$) for
the mutual information with the coefficient $a(N_{A},\beta)$. The latter is
expected to approach a constant value with subsystem size $N_{A}$ for
sufficiently large $N_{A}$. The other term $d(N_{A},\beta)$ can appear due to
the corner contribution or, for a finite system, from the degeneracy of
thermodynamic state, arsing from symmetry breakingMelko  _et al._ (2010);
Singh  _et al._ (2011) or configurational entropy. In our subsystem geometry
[Fig.5], the corner contribution is absent. The finite-temperature Mott
transition is similar to a liquid-gas transition Kotliar  _et al._ (2000).
Thus, for the second Rényi mutual information, constant term $d(\beta,N_{A})$
can appear between $T_{c}<T<2T_{c}$ and $T<T_{c}$ from effective Ising-like
symmetry breaking in different parts of $I(A,B)$ along the first-order
transition line in the $T-U$ plane. We show $I(A,B)/N_{y}$ for $N=16\times 16$
lattice as a function of $N_{A}$ in Fig.11 for different $U$s across Mott
transition at several temperatures near Mott critical point. The arrows in the
plots indicate increasing values of interaction over the range $U=8-11$. In
Fig.11, we see that $I(A,B)/N_{y}$ becomes more or less independent of $N_{A}$
for $N_{A}\sim N/2$, except near the critical point $(U_{c},T_{c})$, where
more complex dependence on $N_{A}$ is seen.

Figure 12: The non-monotonic temperature dependence of $I(A,B)/N_{y}$ is shown
for $U=10.0\approx U_{c}$. This is computed from system size $16\times 16$.
Here $T$ is plotted in the logarithmic scale.

Finally, in the Fig.12, we show the mutual information $I(A,B)/N_{y}$ for the
bi-partition $N_{A}=N_{B}$ as a function of temperature (in logarithmic scale)
near $U\simeq U_{c}=10.0$. The $I(A,B)$ shows a non-monotonic behaviour with
$T$. In particular the mutual information appears to dip between $\sim T_{c}$
and $\sim 2T_{c}$. A substantial $I(A,B)$, indicating correlations, persists
even far above $T_{c}$ till $T\lesssim W$ in the supercritical regime Walsh
_et al._ (2019a, b, 2020). In future, it will be interesting to study the
system-size scaling of $I(A,B)$ with $N_{y}$ or $N$ to understand the nature
of correlations in the supercritical regime contributing to the mutual
information of an extended subsystem.

##  VIII Conclusions and discussion

To summarize, we have developed a DMFT formalism and its numerical
implementation for computing Rényi entanglement entropy and mutual information
in strongly correlated electronic systems described by the Hubbard model. We
show that the scaling of the Rényi entropy with subsystem size for an extended
subsystem can be used to characterize correlated Fermi liquid metallic state
in the half-filled Hubbard model. In particular, we show that the subsystem-
size scaling of Rényi entropy follows the entropy to entanglement crossover
formula expected from CFT and related arguments even in the presence of strong
electronic correlations captured by local self-energy approximation within the
single-site DMFT. We also show how the first-order Mott transition and Mott
critical point are manifested in the temperature, interaction and subsystem-
size dependence of Rényi mutual information.

Here, as a first attempt to implement the entanglement path integral formalism
of ref.Haldar  _et al._ , 2020 within DMFT, we use an approximate impurity
solver, namely the IPT Georges  _et al._ (1996). An immediate extension of our
work would be to employ the CTQMC impurity solverGull  _et al._ (2011) in the
entanglement DMFT framework. Our entanglement path integral formalism is
naturally suited for such a purpose and only requires the incorporation of the
local kick self-energy [Eq.(III)] in the impurity action for the CTQMC
solverGull  _et al._ (2011). Another interesting, albeit more challenging,
future direction would be to capture short-range correlations via cluster
extensionKotliar  _et al._ (2001, 2006) of the DMFT formalism. On a different
note, it will be interesting to explore the connections between real-space and
momentum-spaceFlynn  _et al._ (2023) entanglement in a Fermi liquid.

In this work, we have discussed the implementation of the entanglement path
integral for thermal density matrix suitable to approach entanglement of many-
body ground states of the Hubbard model in the $T\to 0$ limit. However, like
the usual non-equilibrium DMFTAoki  _et al._ (2014), the entanglement path
integral and the DMFT can be extended to non-equilibrium situations via
Schwinger-Keldysh formalism, as discussed in ref.Haldar  _et al._ , 2020;
Moitra and Sensarma, 2020. Moreover, the non-equilibrium formulation can be
further generalized to incorporate non-unitary dynamics, e.g., the dynamics
and the steady states of Hubbard model under repeated projective or weak
measurements, to study entanglement transitions similar to that seen in random
quantum circuitsFisher  _et al._ (2022). In recent years, DMFT, with its
integration with other first-principle electronic structure methods Kotliar
_et al._ (2006), has become one the most practical approaches to describe
realistic strongly correlated systems. Thus, our DMFT formulations, along with
its possible extensions discussed above, might lead to a viable route to
computing entanglement properties of strongly correlated materials in future.

###### Acknowledgements.

SB acknowledges support from SERB (ECR/2018/001742), DST, India and QuST, DST,
India.

##  Appendix A Numerical solution of the DMFT equations

The DMFT equations [Eqs.(III, 20, III, III, 24, 26)] for computing
$S_{A}^{(2)}$ are numerically much more challenging to solve compared to usual
equilibrium DMFT equations Georges  _et al._ (1996) since the Green’s function
$G_{i\sigma\alpha,j\sigma\beta}(\tau,\tau^{\prime})$ [Eq.(13)] is a matrix in
both space and time. The space translation symmetry is broken by the
entanglement cut, and the time translation symmetry is broken by choice of
time $\tau_{0}$ for inserting the auxiliary fields in the path integral
[Eq.(4)]. However, for the cylindrical subsystem geometry in Fig.5, considered
for the 2d Hubbard model, the system retains the translation symmetry parallel
to the ($y$) direction of the entanglement cut. In this case, one can use
Fourier transform along the $y$ direction, as we discuss later.

To solve the DMFT equations in imaginary time $\tau$ without time-translation
invariance, we discretize the DMFT equations in imaginary time. We divide the
time interval $[0,\beta)$ at inverse temperature $\beta=1/T$ into $N_{\tau}$
segments with discretization step $\delta\tau=\beta/N_{\tau}$. However, while
discretizing we have to ensure the appropriate anti-periodic boundary
conditions on the fermionic Green’s function, namely

$\displaystyle G(\tau+\beta,\tau^{\prime})$
$\displaystyle=-G(\tau,\tau^{\prime})$ (34a) $\displaystyle
G(\tau,\tau^{\prime}+\beta)$ $\displaystyle=-G(\tau,\tau^{\prime}),$ (34b)

which is equivalent to the anti-periodic boundary conditions on Grassmann
variables, $c_{N_{\tau}}=-c_{0}$ and $\bar{c}_{N_{\tau}}=-\bar{c}_{0}$, in the
time-discretized form. Here we have suppressed the space, spin and replica
indices for brevity. We use the indices $(n,m)$ running from $n,m=0$ to
$N_{\tau}$ for $(\tau,\tau^{\prime})$. We also write the following useful
discretization rules

$\displaystyle(\partial_{\tau}\bar{c}(\tau))c(\tau)$
$\displaystyle=\frac{(\bar{c}_{n+1}-\bar{c}_{n})c_{n}}{\delta\tau},$ (35a)
$\displaystyle\bar{c}(\tau)c(\tau)$ $\displaystyle=\bar{c}_{n+1}c_{n}.$ (35b)

The above rules arise since the creation operator $\bar{c}$ always appears
slightly later in time than the annihilation operator $c$ in the path
integral. Using the above rules, e.g., we can write Eq.(III) as

$\displaystyle S_{\lambda}$
$\displaystyle=\delta\tau^{2}\sum^{N_{\tau}}_{n,m=1}\sum_{ij,\sigma,\alpha\beta,nm}\bar{c}_{i\sigma\alpha
n}[-G^{-1}_{0,i\alpha n,j\beta m}]c_{j\sigma\beta m}$ (36)
$\displaystyle+U\delta\tau\sum^{N_{\tau}}_{n=1}\sum_{i\alpha}\bar{c}_{i\uparrow\alpha,n+1}c_{i\uparrow\alpha
n}\bar{c}_{i\downarrow\alpha,n+1}c_{i\downarrow\alpha n}.$

The inverse of lattice Green’s function appearing above is given by

$\displaystyle-G^{-1}_{0,i\alpha n,j\beta m}$ $\displaystyle=g^{-1}_{\alpha
m,\beta
n}\delta_{ij}+\frac{1}{\delta\tau}t_{ij}\zeta_{mn}\delta_{\alpha\beta}$
$\displaystyle+\lambda\delta_{i\in
A}\delta_{ij}M_{\alpha\beta}\frac{\delta_{m,p+1}\delta_{np}}{\delta\tau^{2}}$
(37)

where the index $p\in[0,N_{\tau})$ is arbitrary depending on $\tau_{0}$. Here

$\displaystyle g^{-1}_{\alpha m,\beta n}$
$\displaystyle=\frac{1}{\delta\tau^{2}}(\zeta_{mn}-\delta_{mn})\delta_{\alpha\beta}-\mu\frac{1}{\delta\tau}\zeta_{mn}\delta_{\alpha\beta}$
(38a) with $\displaystyle\zeta_{mn}$ $\displaystyle=\delta_{m,n+1}\hskip
11.38092ptn<N_{\tau}$ $\displaystyle=-1\hskip 11.38092ptn=N_{\tau}$ (38b)

Similarly, Eq.(III) becomes

$\displaystyle\mathcal{G}_{i\alpha m,\beta n}^{-1}$
$\displaystyle=g^{-1}_{\alpha m,\beta n}-\Delta_{i\alpha m,\beta
n}-\lambda\delta_{i\in
A}M_{\alpha\beta}\frac{\delta_{m,p+1}\delta_{np}}{\delta\tau^{2}}$ (39)

The IPT self-energy [Eq.(III)] is obtained as

$\displaystyle\Sigma_{i\alpha m,\beta n}$ $\displaystyle=UG_{i\alpha m,\beta
n}\frac{\delta_{m,n-1}}{\delta\tau}-U^{2}\mathcal{\tilde{G}}_{i\alpha m,\beta
n}^{2}\mathcal{\tilde{G}}_{i\beta m,\alpha n},$ (40a) where
$\displaystyle\mathcal{\tilde{G}}_{i\alpha m,\beta n}=\mathcal{G}_{i\alpha
m,\beta n}-UG_{ii,\alpha m,\beta n}\frac{\delta_{m,n-1}}{\delta\tau}.$ (40b)

The lattice Green’s function is obtained as

$\displaystyle\sum_{n_{1}}\sum_{k\gamma}G^{-1}_{i\alpha m,k\gamma
n_{1}}G_{k\gamma n_{1},j\beta n}$
$\displaystyle=\delta_{ij}\delta_{\alpha\beta}\frac{\delta_{mn}}{\delta\tau^{2}}$
(41a) where, within the local self-energy approximation, the Dyson equation is
$\displaystyle G^{-1}_{i\alpha m,j\beta n}=G^{-1}_{0,i\alpha m,j\beta
n}-\delta_{ij}\Sigma_{i\alpha m,\beta n}$ (41b)

Now one can obtain the lattice Green’s function from the above equation. The
lattice Green’s function determines the hybridization function through
Eqs.(24,III).

###  A.1 Recursive Green’s function method

The most computationally expensive part of the DMFT steps discussed above is
the inversion of $G^{-1}$, a matrix of dimension $\sim N_{x}N_{\tau}\times
N_{x}N_{\tau}$, to obtain $G$ via Eq.(41a) for a lattice with $N_{x}$ sites
along the ($x$) direction of partitioning. A direct inverse with the DMFT
self-consistency loop is only feasible for relatively small systems of size
$N_{x}\lesssim 12$ with $N_{\tau}\lesssim 1000$. For larger systems, we use a
recursive Green’s function method along the $x$ direction. The recursive
method can be implemented for the open boundary condition (OBC) as well as for
the periodic boundary condition (PBC). In dimension $d>1$ with translation
symmetry, we can make simplifications by using Fourier transform in the
transverse momenta (see the discussion later). To demarcate between directions
with translation symmetry and the direction of recursion $x$, we denote the
lattice sites by $\boldsymbol{i}=(i,\boldsymbol{i}_{\perp})$ below.

To set up the recursive method, we rewrite Eq.(41a) as a matrix equation

$\displaystyle\boldsymbol{\mathrm{G^{-1}}}\boldsymbol{\mathrm{G}}$
$\displaystyle=\boldsymbol{\mathrm{I}}$ (42a)
$\displaystyle\boldsymbol{\mathrm{I}}_{\boldsymbol{i}\alpha
m,\boldsymbol{j}\beta n}$
$\displaystyle=\delta_{\boldsymbol{i}\boldsymbol{j}}\delta_{mn}\delta_{\alpha\beta},$
(42b) where $\displaystyle(\boldsymbol{\mathrm{G}})_{\boldsymbol{i}\alpha
m,\boldsymbol{j}\beta n}$ $\displaystyle=\delta\tau G_{\boldsymbol{i}\alpha
m,\boldsymbol{j}\beta n},$ (42c)

and similarly for $\boldsymbol{\mathrm{G^{-1}}}$. We separate the system
(spatially) as a system “$S$” and the rest “$R$”, and write

$\displaystyle\bigg{[}\begin{pmatrix}(\boldsymbol{\mathrm{G}}_{R})^{-1}&-\boldsymbol{\mathrm{T}}_{RS}\\\
-\boldsymbol{\mathrm{T}}_{RS}&(\boldsymbol{\mathrm{G}}_{S})^{-1}\end{pmatrix}\bigg{]}\begin{pmatrix}\boldsymbol{\mathrm{G}}_{R}^{R+S)}&\boldsymbol{\mathrm{G}}_{RS}^{(R+S)}\\\
\boldsymbol{\mathrm{G}}_{SR}^{(R+S)}&\boldsymbol{\mathrm{G}}_{S}^{(R+S)}\end{pmatrix}$
$\displaystyle=\boldsymbol{\mathrm{I}}$ (43)

Here $\boldsymbol{\mathrm{G}}_{R}$ and $\boldsymbol{\mathrm{G}}_{S}$ are
Green’s functions of the system and the rest in the absence of any coupling
between them. $\boldsymbol{\mathrm{T}}$ connects the systems and the rest, and
$\boldsymbol{\mathrm{G}}^{(R+S)}$ is the full Green’s function of the combined
system. From the above, we get

$\displaystyle\boldsymbol{\mathrm{G}}_{R}^{(R+S)}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{R}+\boldsymbol{\mathrm{G}}_{R}\boldsymbol{\mathrm{T}}_{RS}\boldsymbol{\mathrm{G}}^{(R+S)}_{SR}$
(44a) $\displaystyle\boldsymbol{\mathrm{G}}^{(R+S)}_{RS}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{R}\boldsymbol{\mathrm{T}}_{RS}\boldsymbol{\mathrm{G}}^{(R+S)}_{S}$
(44b) $\displaystyle\boldsymbol{\mathrm{G}}^{(R+S)}_{SR}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{S}\boldsymbol{\mathrm{T}}^{\dagger}_{RS}\boldsymbol{\mathrm{G}}^{(R+S)}_{R}$
(44c) $\displaystyle\boldsymbol{\mathrm{G}}_{S}^{(R+S)}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{S}+\boldsymbol{\mathrm{G}}_{S}\boldsymbol{\mathrm{T}}^{\dagger}_{RS}\boldsymbol{\mathrm{G}}^{(R+S)}_{RS}$
(44d)

####  A.1.1 Derivation of Eq.(III) for the cavity Green’s function

For this case, we take $S$ as a single site $i$ and $R$ as rest of the system.
Using Eqs.(44a,44b), we obtain

$\displaystyle\boldsymbol{\mathrm{G}}_{R}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{R}^{(R+S)}-\boldsymbol{\mathrm{G}}_{RS}^{(R+S)}(\boldsymbol{\mathrm{G}}_{S}^{(R+S)})^{-1}\boldsymbol{\mathrm{G}}_{SR}^{(R+S)}$
(45)

The above leads to Eq.(III) when we identify
$\boldsymbol{\mathrm{G}}^{(i)}=\boldsymbol{\mathrm{G}}_{R}$, i.e., the cavity
Green’s function with $i$-th site removed, and
$\boldsymbol{\mathrm{G}}=\boldsymbol{\mathrm{G}}^{(R+S)}$, the Green’s
function of the whole lattice.

####  A.1.2 Recursive solution of Eq.(42)

We imagine successively building the system along the $x$ direction from the
left, starting from the first layer at $i=l=1$ and then adding successive
layers till $l=N_{x}$. Imagine that at the $l$-th step of recursion, we have
only left $l+1$ layers, and we separate the system into left $l$ layers
(“$L$”), i.e., $R$ of the preceding section, and add one layer (system “$S$”)
more. We denote the Green’s function of the left $l$ layers as
$\boldsymbol{\mathrm{G}}^{(l)}$ and that of $l+1$ layers as
$\boldsymbol{\mathrm{G}}^{(l+1)}$. From Eqs.(42)

$\displaystyle\bigg{[}\begin{pmatrix}(\boldsymbol{\mathrm{G}}^{(l)}_{L})^{-1}&-\boldsymbol{\mathrm{T}}_{LS}\\\
-\boldsymbol{\mathrm{T}}^{\dagger}_{LS}&(\boldsymbol{\mathrm{G}}_{S})^{-1}\end{pmatrix}\bigg{]}\begin{pmatrix}\boldsymbol{\mathrm{G}}_{L}^{(l+1)}&\boldsymbol{\mathrm{G}}_{LS}^{(l+1)}\\\
\boldsymbol{\mathrm{G}}_{SL}^{(l+1)}&\boldsymbol{\mathrm{G}}_{S}^{(l+1)}\end{pmatrix}$
$\displaystyle=\boldsymbol{\mathrm{I}}.$ (46)

where the coupling between $L$ and $S$ is given by

$\displaystyle\boldsymbol{\mathrm{T}}_{\boldsymbol{i},\boldsymbol{j}}$
$\displaystyle=\boldsymbol{\zeta}t_{\boldsymbol{i},\boldsymbol{j}}\hskip
14.22636pti\leq l,j=l+1$ $\displaystyle=0\hskip 14.22636pt\text{otherwise},$
(47a) with $\displaystyle\boldsymbol{\zeta}_{\alpha n,\beta m}$
$\displaystyle=\zeta_{nm}\delta_{\alpha\beta}.$ (47b) The inverse Green’s
functions of $L$ and $S$ in the absence of any coupling between them are
$\displaystyle\bigg{(}\boldsymbol{\mathrm{G}}^{(l)}_{L}\bigg{)}^{-1}_{\boldsymbol{i}\alpha
m,\boldsymbol{j}\beta n}$
$\displaystyle=\boldsymbol{\mathrm{G}}^{-1}_{\boldsymbol{i}\alpha
m,\boldsymbol{j}\beta n}\hskip 14.22636pti,j\leq l$ (47c)
$\displaystyle(\boldsymbol{\mathrm{G}}_{S})^{-1}_{\boldsymbol{i}\alpha
m,\boldsymbol{j}\beta n}$
$\displaystyle=\boldsymbol{\mathrm{G}}^{-1}_{\boldsymbol{i}\alpha
m,\boldsymbol{j}\beta n}\hskip 14.22636pti,j=l+1$ (47d)

Here $i\leq l$ indicates that the site $\boldsymbol{i}$ belongs to a layer
from $1$ to $l$. The full Green’s function that we eventually want to
calculate is $\boldsymbol{\mathrm{G}}=\boldsymbol{\mathrm{G}}^{(N_{x})}$. We
can rewrite Eq.(46) as

$\displaystyle\boldsymbol{\mathrm{G}}_{L}^{(l+1)}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{L}^{(l)}+\boldsymbol{\mathrm{G}}_{L}^{(l)}\boldsymbol{\mathrm{T}}_{LS}\boldsymbol{\mathrm{G}}^{(l+1)}_{SL}$
(48a) $\displaystyle\boldsymbol{\mathrm{G}}^{(l+1)}_{LS}$
$\displaystyle=\boldsymbol{\mathrm{G}}^{(l)}_{L}\boldsymbol{\mathrm{T}}_{LS}\boldsymbol{\mathrm{G}}^{(l+1)}_{S}$
(48b) $\displaystyle\boldsymbol{\mathrm{G}}^{l+1}_{SL}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{s}\boldsymbol{\mathrm{T}}^{\dagger}_{SL}\boldsymbol{\mathrm{G}}^{(l+1)}_{L}$
(48c) $\displaystyle\boldsymbol{\mathrm{G}}_{S}^{(l+1)}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{S}+\boldsymbol{\mathrm{G}}_{S}\boldsymbol{\mathrm{T}}^{\dagger}_{LS}\boldsymbol{\mathrm{G}}^{(l+1)}_{LS}$
(48d)

We now rewrite Eq.(48a) keeping only the index $i$ for the layers, where all
other indices $\alpha$, $m$ and $\boldsymbol{i}_{\perp}$ are implicit in the
matrices and contracted for matrix multiplications.

$\displaystyle\boldsymbol{\mathrm{G}}^{(l+1)}_{i,j}$
$\displaystyle=\boldsymbol{\mathrm{G}}^{(l)}_{i,j}+\sum_{i_{1}\leq
l}\boldsymbol{\mathrm{G}}^{(l)}_{i,i_{1}}\boldsymbol{\mathrm{T}}_{i_{1},l+1}\boldsymbol{\mathrm{G}}^{(l+1)}_{l+1,j}\hskip
11.38092pti,j\leq l$ (49)

From Eq.(48b), we get

$\displaystyle\boldsymbol{\mathrm{G}}^{(l+1)}_{i,l+1}$
$\displaystyle=\sum_{i_{1}\leq
l}\boldsymbol{\mathrm{G}}^{(l)}_{i,i_{1}}\boldsymbol{\mathrm{T}}_{i_{1},l+1}\boldsymbol{\mathrm{G}}^{(l+1)}_{l+1,l+1}\hskip
11.38092pti\leq l$ (50)

Since $G^{*}_{i\alpha m,j\beta n}=G_{j\beta n,i\alpha m}$, we can obtain from
the above

$\displaystyle\boldsymbol{\mathrm{G}}^{(l+1)}_{l+1,i}$
$\displaystyle=\sum_{i_{1}\leq
l}\boldsymbol{\mathrm{G}}^{(l+1)}_{l+1,l+1}\boldsymbol{\mathrm{T}}_{l+1,i_{1}}\boldsymbol{\mathrm{G}}^{(l)}_{i_{1},i}\hskip
11.38092pti\leq l$ (51)

Using the above in Eq.(49), we obtain for $i,j\leq l$

$\displaystyle\boldsymbol{\mathrm{G}}^{(l+1)}_{i,j}$
$\displaystyle=\boldsymbol{\mathrm{G}}^{(l)}_{i,j}+\sum_{i_{1},i_{2}\leq
l}\boldsymbol{\mathrm{G}}^{(l)}_{i,i_{1}}\boldsymbol{\mathrm{T}}_{i_{1},l+1}\boldsymbol{\mathrm{G}}^{(l+1)}_{l+1,l+1}\boldsymbol{\mathrm{T}}_{l+1,i_{2}}\boldsymbol{\mathrm{G}}^{(l)}_{i_{2},j}$
(52)

In the above equation, the only unknown quantity is
$\boldsymbol{\mathrm{G}}^{(l+1)}_{l+1,l+1}$. This can be obtained as follows.
From Eqs.(48b,48d), we get

$\displaystyle\boldsymbol{\mathrm{G}}^{(l+1)}_{S}$
$\displaystyle=\boldsymbol{\mathrm{G}}_{S}+\boldsymbol{\mathrm{G}}_{S}\boldsymbol{\mathrm{T}}^{\dagger}_{LS}\boldsymbol{\mathrm{G}}^{(l)}_{L}\boldsymbol{\mathrm{T}}_{LS}\boldsymbol{\mathrm{G}}^{(l+1)}_{S}$
(53)

The above can be written in the form of a Dyson equation,

$\displaystyle\big{(}\boldsymbol{\mathrm{G}}^{l+1}_{S}\big{)}^{-1}$
$\displaystyle=(\boldsymbol{\mathrm{G}}^{(l+1)}_{l+1,l+1}\big{)}^{-1}=\boldsymbol{\mathrm{G}}_{S}^{-1}-\boldsymbol{\Sigma}^{(l)}$
(54)

with the self energy

$\displaystyle\boldsymbol{\Sigma}^{(l)}$ $\displaystyle=\sum_{i_{1},i_{2}\leq
l}\boldsymbol{\mathrm{T}}_{l+1,i_{1}}\boldsymbol{\mathrm{G}}^{(l)}_{i_{1},i_{2}}\boldsymbol{\mathrm{T}}_{i_{2},l+1}$
(55)

Hence $\boldsymbol{\mathrm{G}}^{(l+1)}_{l+1,l+1}$ can be obtained using
Eq.(54). Thus from Eqs.(50, 51, 52, 54, 55), we can construct the complete
Green’s function matrix $\boldsymbol{\mathrm{G}}^{(l+1)}_{\boldsymbol{i}\alpha
m,\boldsymbol{j}\beta n}$ ($i,j\leq l+1$) of the system of $l+1$ layers from
that of $l$ layers. The process can be applied recursively, starting with
$l=1$ and continuing till $l=N_{x}-1$, which will yield us the Green’s
function of system size $N_{x}$, i.e., $\boldsymbol{\mathrm{G}}^{(N_{x})}$.
However, for the DMFT self-consistency Eqs.(III, 20, III, III, 24, 26), one
does not need the full Green’s function at each DMFT iteration, only certain
elements. In particular, if we only consider the nearest neighbor hopping, we
will need to keep track of the onsite, nearest, and next-nearest neighbor
Green’s functions to complete the DMFT loop using the cavity Eqs.(24,III). For
the Bethe lattice approximation, nearest-neighbour sites of $i$-th site are
disconnected, and from Eq.(24) we get

$\displaystyle\Delta_{i,\alpha\beta}(\tau,\tau^{\prime})=t^{2}\sum_{j}^{\prime}G^{(i)}_{j\alpha,j\beta}(\tau,\tau^{\prime}).$
(56)

Here $\sum^{\prime}_{j}$ indicates only summation over nearest neighbours of
$i$. Furthermore, in the limit of large connectivity Georges  _et al._ (1996),
$\boldsymbol{\mathrm{G}}^{(i)}=\boldsymbol{\mathrm{G}}$. Thus, in the large-
connectivity Bethe lattice approximation, we only need to compute onsite
elements of the lattice Green’s function during the DMFT self-consistency
loop.

For periodic boundary condition (PBC), we need to incorporate the hopping
matrix element between site $1$ and site $N_{x}$. We can implement this in the
recursive procedure by changing hopping coupling matrix in the last iteration
accordingly when we add the $l=N_{x}-1$ layer with single layer system $S$ to
form the required system size $l+1=N_{x}$. In particular, we can explicitly
write the hopping coupling matrix for nearest neighbour for PBC below

$\displaystyle\boldsymbol{\mathrm{T}}_{\boldsymbol{i},\boldsymbol{j}}$
$\displaystyle=\boldsymbol{\zeta}t_{\boldsymbol{i},\boldsymbol{j}}\hskip
14.22636pti=l;j=l+1\hskip 7.11317pt\text{if}\hskip 7.11317ptl\leq N_{x}-2$
$\displaystyle=\boldsymbol{\zeta}t_{\boldsymbol{i},\boldsymbol{j}}\hskip
14.22636pti=1,l;j=l+1\hskip 7.11317pt\text{if}\hskip 7.11317ptl=N_{x}-1$
$\displaystyle=0\hskip 14.22636pt\text{otherwise},$ (57)

Figure A1: The convergence of $S_{A}^{(2)}$ with respect to discretization
$\delta\lambda$ to evaluate the integral in Eq.(12) for $U=2.0$, $T=0.05$ and
system size $N=30$ in 1d Hubbard model.

##  Appendix B Calculation of $S_{A}^{(2)}$ from “kick term” integration
method

We numerically solve the self-consistent DMFT Eqs.(39, 40a, 41a,24) for
discretized values of $\lambda\in[0,1]$ with uniform step $\delta\lambda$ to
obtain $G_{i\sigma\alpha,j\sigma\beta}(\tau_{0},\tau^{+}_{0})$ [Eq.(13)],
where we choose $\tau_{0}=0$. We then compute $\langle S_{\rm
kick}\rangle_{Z_{A}^{(2)}(\lambda)}$ using Eq.(12). To obtain $S_{A}^{(2)}$,
we integrate $\langle S_{\rm kick}\rangle_{Z_{A}^{(2)}(\lambda)}$ over
$\lambda$ from 0 to 1 using numerical interpolation over the range of
$\lambda$.

Figure A2: The $S_{A}^{(2)}(\delta\tau)$ obtained using the direct and
recursive inversion in the non-interacting case is shown for $T=0.3,0.5$ and
system size $N=10$.

We have used $\delta\lambda=0.1$ for most of our calculations. We have varied
$\delta\lambda$ to check the convergence of $S_{A}^{(2)}$ with
$\delta\lambda$. We benchmark the $S_{A}^{(2)}$ computed this way in the non-
interacting case by comparing with $S_{A}^{(2)}$ obtained directly from the
correlation matrix calculations, as discussed in Sec.IV.1. For the interacting
case, we have numerically checked the convergence by taking different values
of $\delta\lambda$, and various numerical interpolation schemes. As an
example, in Fig.A1, we show the convergence of the $S_{A}^{(2)}$ with
different $\delta\lambda=0.02,0.04,0.1$, for $U=2.0$, $T=0.05$, and $N=30$ in
1d.

##  Appendix C Benchmark of the recursive Green’s function method

We can invert $\boldsymbol{\mathrm{G}}^{-1}$ in Eq.(41a) directly for small
systems ($N\lesssim 12$) to obtain the lattice Green’s function. We use the
recursive method for larger systems to invert $\boldsymbol{\mathrm{G}}^{-1}$
and obtain onsite lattice Green’s function within the DMFT loop. We benchmark
the recursive method by comparing computed $S_{A}^{(2)}$ with that obtained
from direct inversion for $N=10$ in 1d, as shown in Fig.A2 for the non-
interacting case.

Figure A3: The linear extrapolations of $S_{A}^{(2)}(\delta\tau)$ are shown
here for few subsystem sizes. The $S_{A}^{(2)}(\delta\tau)$ as a function of
$\delta\tau$ and the corresponding linear extrapolation (shown as ‘extrap.’ in
legend) are shown here in (a) for $U=0.5$ and subsystem sizes
$N_{A}=2,7,10,15$ from system size $N=30$, and in (b) for $U=2$ and subsystem
sizes $N_{A}=3,8,11,16$ from system $N=50$ for at $T=0.05$.

##  Appendix D The extrapolation of $S_{A}^{(2)}(\delta\tau)$ to
$\delta\tau\to 0$ limit

As discussed in appendix A, we solve the DMFT self-consistency equations by
discretizing them in imaginary time with discretization step
$\delta\tau=\beta/N_{\tau}$. To approach the continuum limit $\delta\tau\to
0$, we compute $S_{A}^{(2)}(\delta\tau)$ for a few values $\delta\tau$ and
then linearly extrapolate it to $\delta\tau\to 0$. We take $\delta\tau$ over
the range $0.015$ to $0.075$. In most of our calculations, we take four values
of $\delta\tau$ in the above range, particularly between $0.02$ to $0.04$, and
then do the linear extrapolation to $\delta\tau\to 0$. It becomes
progressively more challenging to take $\delta\tau$ in the above range for low
temperatures $T<0.05$ as the size ($2NN_{\tau}\times 2NN_{\tau}$) of the
Green’s function matrix becomes very large. Hence we restrict our DMFT
calculations up to $T=0.05$. We show $S_{A}^{(2)}(\delta\tau)$ as a function
of $\delta\tau$ with the linear extrapolation in Figs.A3(a,b) for a few
subsystem sizes for interactions $U=0.5,2$ and system size $N=30$ as an
illustration of the linear extrapolation.

Figure A4: The $S_{A}^{(2)}$ vs. $N_{A}$ in 1d Hubbard model for $U=2$ and
$N=50$ is shown for different temperatures with the corresponding fit to
Eq.(29) with $c,v,b$ as free parameters.

##  Appendix E  Entanglement to entropy crossover in 1d Hubbard model

Here we discuss the fitting of $S_{A}^{(2)}$ in 1d Hubbard model with the
crossover function of Eq.(29). We first fit $S_{A}^{(2)}(N_{A},T)$ with
Eq.(29) by varying all the parameters $c,v,b$, as shown in Fig.(A4). The
variations of the extracted fitting parameters $c,v,b$ with temperature are
shown in Fig.(A5). As evident from the figure, the expected CFT crossover
formula [Eq.(29)] describes $S_{A}^{(2)}$ for the DMFT metallic state of 1d
Hubbard model quite well. The extracted central charge slowly approaches the
CFT value $c=1$ with decreasing temperature [Fig.(A5)(a)] and converges well
with system size $N$ [Fig.(A5)(d)].

Figure A5: The parameters $c,v,b$ extracted from crossover function [Eq.(29)]
fitting, e.g., in Fig.A4, are shown here for 1d Hubbard model. (a) The central
charge $c$ as a function of $T$ for a few $U$ is shown and compared with the
CFT value $c=1$. (b) The renormalized Fermi velocity $v$ as a function of $T$
is shown for different $U$, and compared with the non-interacting value
$v_{F}$. (c) The non-universal parameter $b$ as a function of $T$. (d) System
size dependence of extracted $c$ vs. $U$ at $T=0.1$.

Here we have treated $c,v,b$ in Eq.(29) as fitting parameters to describe
$S_{A}^{(2)}(N_{A},T)$ in the 1d Hubbard model. In the next section, we first
fix the ratio $(c/v)$ using the specific heat calculated from equilibrium DMFT
and then fit our results for $S_{A}^{(2)}$ with Eq.(29) treating only $c$ and
$b$ as fitting parameters.

###  E.1 Calculation of the ratio $(c/v)$ from equilibrium DMFT

At low temperature ($T\to 0$), the specific heat $c_{V}$ can be obtained from
CFT Korepin (2004); Calabrese and Cardy (2004) as

$\displaystyle c_{V}=\frac{\pi T}{3}\left(\frac{c}{v}\right).$ (58)

We compute the specific heat $c_{V}$ for the paramagnetic metallic state in 1d
Hubbard model from equilibrium DMFT calculation, as described below.

In equilibrium, due to time translation symmetry,
$G_{ij}(\tau,\tau^{\prime})=G_{ij}(\tau-\tau^{\prime})$. Thus we can write the
DMFT self-consistency equations Georges  _et al._ (1996) as follows

$\displaystyle\mathcal{G}_{i}^{-1}(\tau-\tau^{\prime})$
$\displaystyle=-(\partial_{\tau}-\mu)\delta(\tau-\tau^{\prime})-\Delta_{i}(\tau-\tau^{\prime})$
(59) $\displaystyle G^{-1}_{i}(\tau)$
$\displaystyle=\mathcal{G}_{i}^{-1}(\tau)-\Sigma_{i}(\tau)$ (60)

where $\mathcal{G}_{i}(\tau)$ and $G_{i}(\tau)$ are the bare and full impurity
Green’s functions at site $i$. Furthermore, since the model [Eq.(14)] is space
translation invariant, all the sites are equivalent, unlike for the DMFT in
the presence of entanglement cut in Sec.III.

Figure A6: The specific heat ($c_{V}$) as a function of temperature $T$ for 1d
Hubbard model is shown for different interaction $U=0.5,2,3$. The low-
temperature linear fit to $c_{V}(T)\propto T$ is shown and the extracted
slopes are indicated in the legends.

The hybridization function is given by

$\displaystyle\Delta_{i}(\tau)=\sum_{jl}t_{ij}t_{il}G^{(i)}_{jl}(\tau)$ (61)

We use the large-connectivity Bethe lattice approximation for the cavity
Green’s function, i.e., $G^{(i)}_{jl}=G_{jl}$, and
$G^{(i)}_{jl}\simeq\delta_{jl}G_{jj}$. Therefore the hybridization function
for nearest-neighbor hopping becomes

$\displaystyle\Delta_{i}(\tau)=zt^{2}G(\tau)$ (62)

where $G(\tau)$ is the onsite lattice Green’s function and $z=2d$ is the
coordination number for $D$-dimensional hypercubic lattice. The lattice
Green’s function is obtained using the local self-energy approximation, i.e.,
the lattice self-energy $\Sigma_{ij}$ is replaced by the impurity self-energy,
$\Sigma_{ij}=\Sigma\delta_{ij}$, so that

$\displaystyle G(i\omega_{m})=\int
d\epsilon\frac{g(\epsilon)}{\imath\omega_{m}+\mu-\epsilon-\Sigma(\imath\omega_{m})}$
(63)

where $\omega_{m}=(2m+1)\pi/\beta$, with $m$ an integer, is the fermionic
Matsubara frequency, and $g(\epsilon)$ is the non-interacting density of
states per site. The self-energy within IPT approximation Georges  _et al._
(1996) is given by

$\displaystyle\Sigma(\tau)=Un-U^{2}\mathcal{\tilde{G}}^{2}(\tau)\mathcal{\tilde{G}}(-\tau)$
(64) $\displaystyle\mathcal{\tilde{G}}^{-1}(\tau)=\mathcal{G}^{-1}(\tau)-Un$
(65)

where $n$ is the occupation number of a site; $n=1/2$ at half filling. For
calculating thermodynamic properties such as specific heat, we solve the above
equilibrium DMFT self-consistency equation following standard procedure
Georges  _et al._ (1996) to obtain $G(\imath\omega_{m})$ and
$\Sigma(\imath\omega_{m})$.

Using these, we compute the internal energy $E$ Georges  _et al._ (1996) from

$\displaystyle\frac{E}{N}$
$\displaystyle=2T\sum_{m}\int^{\infty}_{-\infty}d\epsilon\frac{\epsilon
D(\epsilon)}{i\omega_{m}+\mu-\Sigma_{(}i\omega_{m})-\epsilon}$
$\displaystyle+T\sum_{m}\Sigma(i\omega_{m})G(i\omega_{m}).$ (66)

The specific heat per site is obtained by taking the numerical derivative of
internal energy, i.e., $c_{V}=(1/N)\big{(}\partial E/\partial T\big{)}$.

Figure A7: Temperature dependence of (a) $b$ and (b) $v$ extracted by fixing
$(c/v)$ ratio from specific heat and fitting CFT crossover formula [Eq.(29)]
to the DMFT results for $S_{A}^{(2)}(N_{A},T)$ in the 1d Hubbard model for
$U=0.0,0.5,2$ and system size $N=50$. The non-interacting Fermi velocity
$v_{F}$ is shown in plot (b) for comparision with renormalised velocity $v$.

We calculate the specific heat as a function of temperature for 1d and 2d
Hubbard models. We show $c_{V}$ vs. $T$ in 1d Hubbard model in Fig.A6. We
extract $(c/v)$ from the slope of the linear fit to $c_{V}(T)$ at low
temperature using Eq.(58), as shown in the figure. For the half-filled 2d
square-lattice Hubbard model with nearest-neighbor hopping, due to the van
Hove singularity of the non-interacting band at the Fermi energy $c_{V}$ has
$\ln T$ correction to Eq.(58), and the $(c/v)$ ratio cannot be estimated
reliably. Thus we only use the $(c/v)$ ratio in 1d to fit
$S_{A}^{(2)}(N_{A},T)$ with the crossover formula [Eq.(29)], as discussed in
Sec.V. The central charge $c$ extracted this way is shown in Fig.4. We show
the non-universal constant $b$ and the velocity $v$, obtained using the
$(c/v)$ ratio and $c$, in Figs.A7(a,b). As evident, $v$ extracted this way for
weak interaction matches quite well at low temperature with non-interacting
$v_{\mathrm{F}}$, unlike the $v$ extracted by fitting the CFT formula with
three parameters $c,~{}v,~{}b$ in Fig.A5(b). For the 2d Hubbard model, we
extract $c$ by fitting the crossover formula [Eq.(29)] to the computed
$S_{A}^{(2)}(N_{A},T)$ using $c,v,b$ as free parameters, as discussed in
Sec.VI.

##  Appendix F $S_{A}^{(2)}(N_{A},T)$ for open boundary condition (OBC) in 1d
Hubbard model

Figure A8: The $S_{A}^{(2)}$ as a function of $N_{A}$ is shown here for
different $T$ for open boundary condition (OBC) for $U=2.0$. The system size
is $N=30$ here.

We show the second Rényi entropy $S_{A}^{(2)}$ computed via DMFT for open
boundary condition (OBC) in the 1d Hubbard model for system size $N=30$ in
Fig.A8. We see quite large oscillations in $S_{A}^{(2)}$ between the odd and
even subsystem sizes at low temperatures for OBC. Such oscillations are
present for periodic boundary condition also, e.g., in Fig.3, but are much
weaker. These oscillations, with frequency $2k_{\mathrm{F}}$ determined by the
Fermi wave vector $k_{\mathrm{F}}$, are expectedCalabrese  _et al._ (2010);
Swingle  _et al._ (2013) due to the subleading corrections to the CFT result
[Eq.(32)], and appear to be enhanced in OBC compared to that in PBC.

##  Appendix G DMFT equations for computing second Rényi entropy in 2d

In two dimensions (2d), we take a cylindrical geometry for the subsystem $A$
to compute the second Rényi entropy, as shown in Fig.5. We take the
entanglement cut parallel to the $y$ axis, i.e., partition the system along
the $x$ direction. The most difficult part in solving the non-equilibrium DMFT
Eqs.(III, 20, III, III, 24, 26) is the inversion of the inverse lattice
Green’s function. We rewrite Eq.(III) below as

$\displaystyle\int_{0}^{\beta}$ $\displaystyle
d\tau^{\prime\prime}\sum_{\boldsymbol{r}^{\prime\prime}\gamma}\big{[}G^{-1}_{0,\boldsymbol{r}\alpha,\boldsymbol{r}^{\prime\prime}\gamma}(\tau,\tau^{\prime\prime})-\delta_{\boldsymbol{r},\boldsymbol{r}^{\prime\prime}}\Sigma_{\boldsymbol{r},\alpha\gamma}(\tau,\tau^{\prime\prime})]$
$\displaystyle
G_{\boldsymbol{r}^{\prime\prime}\gamma,\boldsymbol{r}^{\prime}\beta}(\tau^{\prime\prime},\tau^{\prime})=\delta_{\boldsymbol{r},\boldsymbol{r}^{\prime}}\delta_{\alpha\beta}\delta(\tau-\tau^{\prime})$
(67)

where $\boldsymbol{r}=(x,y)$ represents two dimensional co-ordinates. Due to
translation symmetry in the $y$ direction,
$G_{\boldsymbol{r}\alpha,\boldsymbol{r}^{\prime}\beta}(\tau,\tau^{\prime})=G_{x\alpha,x^{\prime}\beta,y-y^{\prime}}(\tau,\tau^{\prime})$.
As a result, Green’s function can be represented using Fourier transform along
the $y$ direction with momentum $k_{y}$,

$\displaystyle
G_{x\alpha,x^{\prime}\beta,y-y^{\prime}}(\tau,\tau^{\prime})=\frac{1}{N_{y}}\sum_{k_{y}}e^{-ik_{y}(y-y^{\prime})}G_{x\alpha,x^{\prime}\beta,k_{y}}(\tau,\tau^{\prime}).$
(68)

Thus, from Eq.(G), we can write

$\displaystyle\int_{0}^{\beta}$ $\displaystyle
d\tau^{\prime\prime}\sum_{x^{\prime\prime}\gamma}\big{[}G^{-1}_{0,x\alpha,x^{\prime\prime}\gamma,k_{y}}(\tau,\tau^{\prime\prime})-\delta_{xx^{\prime\prime}}\Sigma_{x,\alpha\gamma}(\tau,\tau^{\prime\prime})]$
$\displaystyle
G_{x^{\prime\prime}\gamma,x^{\prime}\beta,k_{y}}(\tau^{\prime\prime},\tau^{\prime})=\delta_{xx^{\prime}}\delta_{\alpha\beta}\delta(\tau-\tau^{\prime})$
(69)

for each $k_{y}$ mode, where

$\displaystyle G^{-1}_{0,x\alpha,x^{\prime}\gamma,k_{y}}(\tau,\tau^{\prime})=$
$\displaystyle\big{[}(-\partial_{\tau}+\mu-2t\cos{k_{y}})\delta_{xx^{\prime}}\big{.}$
$\displaystyle\big{.}-t_{xx^{\prime}}\big{]}\delta_{\alpha\beta}\delta(\tau-\tau^{\prime})$
(70)

where $t_{xx^{\prime}}$ is hopping amplitude along $x$ direction. For nearest-
neighbor hopping, we have $t_{x,x\pm 1}=t$, and $t_{xx^{\prime}}=0$ otherwise.
The hybridization function with large-connectivity Bethe lattice approximation
is given by

$\displaystyle\Delta_{x,\alpha\beta}$
$\displaystyle=t^{2}\bigg{[}G_{x-1,\alpha,x-1,\beta,y-y^{\prime}=0}+G_{x+1,\alpha,x+1,\beta,y-y^{\prime}=0}\bigg{.}$
$\displaystyle\bigg{.}+2G_{x\alpha,x\beta,y-y^{\prime}=0}\bigg{]}.$ (71)

In the above equation, we have omitted the time arguments
$(\tau,\tau^{\prime})$ for notational convenience. The Green’s function
$G_{x\alpha,x\beta,y-y^{\prime}=0}(\tau,\tau^{\prime})$ is obtained from
Eq.(68) as

$\displaystyle
G_{x\alpha,x\beta,y-y^{\prime}=0}(\tau,\tau^{\prime})=\frac{1}{N_{y}}\sum_{k_{y}}G_{x\alpha,x\beta,k_{y}}(\tau,\tau^{\prime})$
(72)

We use the recursive Green’s function method in the $x$ direction for
$G_{x\alpha,x\beta,k_{y}}(\tau,\tau^{\prime})$, as described in Appendix A, to
obtain the lattice Green’s function for each $k_{y}$ mode.

##  Appendix H Widom formula for $S^{(2)}_{A}(N_{A},T)$ in 2d

Figure A9:  A schematic of the Fermi sea at half-filling (in gray color) for
non-interacting 2d tight binding model. The real-space subsystem of
cylindrical geometry is shown in the top left corner. $\hat{n}_{k}$ is a unit
vector perpendicular to Fermi-surface and $\hat{n}_{x}$ is perpendicular to
real-space boundaries of the subsystem.

As we discussed in the main text, in the Widom formula, the effective number
of modes from Fermi surface assuming made of independent patches is given by
Eqn.30 and we rewrite it here

$\displaystyle N_{\rm modes}=\frac{1}{(2\pi)^{d-1}}\frac{1}{2}\int_{\partial
A_{x}}\int_{\partial
A_{k}}dA_{x}dA_{k}|\hat{\boldsymbol{n}}_{x}\cdot\hat{\boldsymbol{n}}_{k}|$
(73)

We use the 2d non-interacting dispersion $\varepsilon_{k}=-2t\cos k_{x}-2t\cos
k_{y}$ for nearest neighbour hopping. The Fermi-surface at half filling is
shown in Fig.A9. The unit normal to Fermi-surface
$\hat{\boldsymbol{n}}_{k}=\frac{1}{\sqrt{2}}(\pm\hat{x}\pm\hat{y})$. For
cylindrical geometry, we have two interface parallel to $y$-axis and hence the
unit normal to real space $\hat{\boldsymbol{n}}_{x}=\pm\hat{x}$. Hence, we get
$|\hat{\boldsymbol{n}}_{k}\cdot\hat{\boldsymbol{n}}_{x}|=\frac{1}{\sqrt{2}}$.
Therefore,

$\displaystyle N_{\rm modes}=\frac{1}{4\pi}\times
2N_{y}\times(4\sqrt{2}\pi)\times\frac{1}{\sqrt{2}}=2N_{y}$ (74)

where the contribution $2N_{y}$ comes from integration over real space
boundary and $4\sqrt{2}\pi$ comes from integration over the Fermi-surface.

We fit the numerically computed $S^{(2)}/N_{y}$ data within DMFT to the Widom
crossover formula Eq.(32) and the extracted central charge $c$ as a function
of $T$, $U$ and system size are shown in the main text. In the Fig.A10, the
temperature dependence of renormalized velocity $v$ and non-universal constant
$b$ extracted by fitting to Eq.(32) are shown for different interaction $U$.
This parameters ($v,b$) are shown for system $20\times 20$.

Figure A10: The temperature dependence of (a) $v$, and (b) $b$, extracted by
fitting Widom crossover formula of Eq.(32) to $S^{(2)}/N_{y}$ for different
$U$ computed from system size $20\times 20$.

## References

  * Nielsen and Chuang (2000) Michael A. Nielsen and Isaac L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2000).
  * Laflorencie (2016) Nicolas Laflorencie, “Quantum entanglement in condensed matter systems,” Physics Reports 646, 1–59 (2016).
  * Levin and Wen (2006) Michael Levin and Xiao-Gang Wen, “Detecting topological order in a ground state wave function,” Phys. Rev. Lett. 96, 110405 (2006).
  * Kitaev and Preskill (2006) Alexei Kitaev and John Preskill, “Topological entanglement entropy,” Phys. Rev. Lett. 96, 110404 (2006).
  * Nandkishore and Huse (2015) Rahul Nandkishore and David A. Huse, “Many-body localization and thermalization in quantum statistical mechanics,” Annual Review of Condensed Matter Physics 6, 15–38 (2015).
  * Altman and Vosk (2015) Ehud Altman and Ronen Vosk, “Universal dynamics and renormalization in many-body-localized systems,” Annual Review of Condensed Matter Physics 6, 383–409 (2015).
  * Abanin  _et al._ (2019) Dmitry A. Abanin, Ehud Altman, Immanuel Bloch, and Maksym Serbyn, “Colloquium: Many-body localization, thermalization, and entanglement,” Rev. Mod. Phys. 91, 021001 (2019).
  * Calabrese and Cardy (2004) Pasquale Calabrese and John Cardy, “Entanglement entropy and quantum field theory,” Journal of Statistical Mechanics: Theory and Experiment 2004, P06002 (2004).
  * Vidal and Werner (2002) G. Vidal and R. F. Werner, “Computable measure of entanglement,” Phys. Rev. A 65, 032314 (2002).
  * Eisert  _et al._ (2010) J. Eisert, M. Cramer, and M. B. Plenio, “Colloquium: Area laws for the entanglement entropy,” Rev. Mod. Phys. 82, 277–306 (2010).
  * Casini and Huerta (2009) H Casini and M Huerta, “Entanglement entropy in free quantum field theory,” Journal of Physics A: Mathematical and Theoretical 42, 504007 (2009).
  * Korepin (2004) V. E. Korepin, “Universality of entropy scaling in one dimensional gapless models,” Phys. Rev. Lett. 92, 096402 (2004).
  * Gioev and Klich (2006a) Dimitri Gioev and Israel Klich, “Entanglement entropy of fermions in any dimension and the widom conjecture,” Phys. Rev. Lett. 96, 100503 (2006a).
  * Li  _et al._ (2006) Weifei Li, Letian Ding, Rong Yu, Tommaso Roscilde, and Stephan Haas, “Scaling behavior of entanglement in two- and three-dimensional free-fermion systems,” Phys. Rev. B 74, 073103 (2006).
  * Swingle (2010) Brian Swingle, “Entanglement entropy and the fermi surface,” Phys. Rev. Lett. 105, 050502 (2010).
  * Swingle (2012a) Brian Swingle, “Conformal field theory approach to fermi liquids and other highly entangled states,” Phys. Rev. B 86, 035116 (2012a).
  * Swingle (2012b) Brian Swingle, “Rényi entropy, mutual information, and fluctuation properties of fermi liquids,” Phys. Rev. B 86, 045109 (2012b).
  * Swingle and Senthil (2013) Brian Swingle and T. Senthil, “Universal crossovers between entanglement entropy and thermal entropy,” Phys. Rev. B 87, 045123 (2013).
  * Calabrese and Cardy (2009) Pasquale Calabrese and John Cardy, “Entanglement entropy and conformal field theory,” Journal of Physics A Mathematical General 42, 504005 (2009), arXiv:0905.4013 [cond-mat.stat-mech] .
  * Hastings  _et al._ (2010) Matthew B. Hastings, Iván González, Ann B. Kallin, and Roger G. Melko, “Measuring renyi entanglement entropy in quantum monte carlo simulations,” Phys. Rev. Lett. 104, 157201 (2010).
  * Humeniuk and Roscilde (2012) Stephan Humeniuk and Tommaso Roscilde, “Quantum monte carlo calculation of entanglement rényi entropies for generic quantum systems,” Phys. Rev. B 86, 235116 (2012).
  * Grover (2013) Tarun Grover, “Entanglement of interacting fermions in quantum monte carlo calculations,” Phys. Rev. Lett. 111, 130402 (2013).
  * Assaad  _et al._ (2014) Fakher F. Assaad, Thomas C. Lang, and Francesco Parisen Toldin, “Entanglement spectra of interacting fermions in quantum monte carlo simulations,” Phys. Rev. B 89, 125121 (2014).
  * Broecker and Trebst (2014) Peter Broecker and Simon Trebst, “Rényi entropies of interacting fermions from determinantal quantum monte carlo simulations,” Journal of Statistical Mechanics: Theory and Experiment 2014, P08015 (2014).
  * Wang and Troyer (2014) Lei Wang and Matthias Troyer, “Renyi entanglement entropy of interacting fermions calculated using the continuous-time quantum monte carlo method,” Phys. Rev. Lett. 113, 110401 (2014).
  * Assaad (2015) Fakher F. Assaad, “Stable quantum monte carlo simulations for entanglement spectra of interacting fermions,” Phys. Rev. B 91, 125146 (2015).
  * D’Emidio (2020) Jonathan D’Emidio, “Entanglement entropy from nonequilibrium work,” Phys. Rev. Lett. 124, 110602 (2020).
  * Chakraborty and Sensarma (2021a) Ahana Chakraborty and Rajdeep Sensarma, “Nonequilibrium dynamics of renyi entropy for bosonic many-particle systems,” Phys. Rev. Lett. 127, 200603 (2021a).
  * Haldar  _et al._ (2020) Arijit Haldar, Surajit Bera, and Sumilan Banerjee, “Rényi entanglement entropy of fermi and non-fermi liquids: Sachdev-ye-kitaev model and dynamical mean field theories,” Phys. Rev. Res. 2, 033505 (2020).
  * Moitra and Sensarma (2020) Saranyo Moitra and Rajdeep Sensarma, “Entanglement entropy of fermions from wigner functions: Excited states and open quantum systems,” Phys. Rev. B 102, 184306 (2020).
  * Chakraborty and Sensarma (2021b) Ahana Chakraborty and Rajdeep Sensarma, “Renyi entropy of interacting thermal bosons in the large-$n$ approximation,” Phys. Rev. A 104, 032408 (2021b).
  * Hubbard (1963) J. Hubbard, “Electron Correlations in Narrow Energy Bands,” Proceedings of the Royal Society of London Series A 276, 238–257 (1963).
  * Fazekas (1999) Patrick Fazekas,  _Lecture Notes on Electron Correlation and Magnetism_ (WORLD SCIENTIFIC, 1999) https://www.worldscientific.com/doi/pdf/10.1142/2945 .
  * Auerbach (1994) Assa Auerbach,  _Interacting Electrons and Quantum Magnetism_ (Springer-Verlag, 1994).
  * Georges  _et al._ (1996) Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg, “Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions,” Rev. Mod. Phys. 68, 13–125 (1996).
  * Kotliar  _et al._ (2006) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, “Electronic structure calculations with dynamical mean-field theory,” Rev. Mod. Phys. 78, 865–951 (2006).
  * Aoki  _et al._ (2014) Hideo Aoki, Naoto Tsuji, Martin Eckstein, Marcus Kollar, Takashi Oka, and Philipp Werner, “Nonequilibrium dynamical mean-field theory and its applications,” Rev. Mod. Phys. 86, 779–837 (2014).
  * Hettler  _et al._ (2000) M. H. Hettler, M. Mukherjee, M. Jarrell, and H. R. Krishnamurthy, “Dynamical cluster approximation: Nonlocal dynamics of correlated electron systems,” Phys. Rev. B 61, 12739–12756 (2000).
  * Kotliar  _et al._ (2001) Gabriel Kotliar, Sergej Y. Savrasov, Gunnar Pálsson, and Giulio Biroli, “Cellular dynamical mean field approach to strongly correlated systems,” Phys. Rev. Lett. 87, 186401 (2001).
  * Maier  _et al._ (2005) Thomas Maier, Mark Jarrell, Thomas Pruschke, and Matthias H. Hettler, “Quantum cluster theories,” Rev. Mod. Phys. 77, 1027–1080 (2005).
  * Bolech  _et al._ (2003) C. J. Bolech, S. S. Kancharla, and G. Kotliar, “Cellular dynamical mean-field theory for the one-dimensional extended hubbard model,” Phys. Rev. B 67, 075110 (2003).
  * Capone  _et al._ (2004) M. Capone, M. Civelli, S. S. Kancharla, C. Castellani, and G. Kotliar, “Cluster-dynamical mean-field theory of the density-driven Mott transition in the one-dimensional Hubbard model,” Phys. Rev. B 69, 195105 (2004), arXiv:cond-mat/0401060 [cond-mat.str-el] .
  * Gull  _et al._ (2011) Emanuel Gull, Andrew J. Millis, Alexander I. Lichtenstein, Alexey N. Rubtsov, Matthias Troyer, and Philipp Werner, “Continuous-time monte carlo methods for quantum impurity models,” Rev. Mod. Phys. 83, 349–404 (2011).
  * Ryu and Takayanagi (2006) Shinsei Ryu and Tadashi Takayanagi, “Aspects of holographic entanglement entropy,” Journal of High Energy Physics 2006, 045 (2006), arXiv:hep-th/0605073 [hep-th] .
  * Wolf  _et al._ (2008) Michael M. Wolf, Frank Verstraete, Matthew B. Hastings, and J. Ignacio Cirac, “Area laws in quantum systems: Mutual information and correlations,” Phys. Rev. Lett. 100, 070502 (2008).
  * Udagawa and Motome (2015) Masafumi Udagawa and Yukitoshi Motome, “Entanglement spectrum in cluster dynamical mean-field theory,” Journal of Statistical Mechanics: Theory and Experiment 2015, P01016 (2015).
  * Walsh  _et al._ (2019a) C. Walsh, P. Sémon, D. Poulin, G. Sordi, and A.-M. S. Tremblay, “Local entanglement entropy and mutual information across the mott transition in the two-dimensional hubbard model,” Phys. Rev. Lett. 122, 067203 (2019a).
  * Walsh  _et al._ (2019b) C. Walsh, P. Sémon, D. Poulin, G. Sordi, and A.-M. S. Tremblay, “Thermodynamic and information-theoretic description of the mott transition in the two-dimensional hubbard model,” Phys. Rev. B 99, 075122 (2019b).
  * Walsh  _et al._ (2020) C. Walsh, P. Sémon, D. Poulin, G. Sordi, and A.-M. S. Tremblay, “Entanglement and classical correlations at the doping-driven mott transition in the two-dimensional hubbard model,” PRX Quantum 1, 020310 (2020).
  * Cahill and Glauber (1969) K. E. Cahill and R. J. Glauber, “Density operators and quasiprobability distributions,” Phys. Rev. 177, 1882–1902 (1969).
  * Fetter and Walecka (1971) A. L. Fetter and J. D. Walecka,  _Quantum Theory of Many-Particle Systems_ (McGraw-Hill, Boston, 1971).
  * Kyung  _et al._ (2006) B. Kyung, G. Kotliar, and A.-M. S. Tremblay, “Quantum monte carlo study of strongly correlated electrons: Cellular dynamical mean-field theory,” Phys. Rev. B 73, 205106 (2006).
  * Lieb and Wu (1968) Elliott H. Lieb and F. Y. Wu, “Absence of mott transition in an exact solution of the short-range, one-band model in one dimension,” Phys. Rev. Lett. 20, 1445–1448 (1968).
  * Cardy (1996) John L Cardy,  _Scaling and renormalization in statistical physics_, Cambridge lecture notes in physics (Cambridge Univ. Press, Cambridge, 1996).
  * Holzhey  _et al._ (1994) Christoph Holzhey, Finn Larsen, and Frank Wilczek, “Geometric and renormalized entropy in conformal field theory,” Nuclear Physics B 424, 443–467 (1994), arXiv:hep-th/9403108 [hep-th] .
  * Vidal  _et al._ (2003) G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, “Entanglement in quantum critical phenomena,” Phys. Rev. Lett. 90, 227902 (2003).
  * Chen  _et al._ (2020) Xiao Chen, Yaodong Li, Matthew P. A. Fisher, and Andrew Lucas, “Emergent conformal symmetry in nonunitary random dynamics of free fermions,” Phys. Rev. Res. 2, 033017 (2020).
  * Ding  _et al._ (2012) Wenxin Ding, Alexander Seidel, and Kun Yang, “Entanglement entropy of fermi liquids via multidimensional bosonization,” Phys. Rev. X 2, 011012 (2012).
  * Gioev and Klich (2006b) Dimitri Gioev and Israel Klich, “Entanglement entropy of fermions in any dimension and the widom conjecture,” Phys. Rev. Lett. 96, 100503 (2006b).
  * Calabrese  _et al._ (2012) P. Calabrese, M. Mintchev, and E. Vicari, “Entanglement entropies in free-fermion gases for arbitrary dimension,” Europhysics Letters 97, 20009 (2012).
  * Leschke  _et al._ (2014) Hajo Leschke, Alexander V. Sobolev, and Wolfgang Spitzer, “Scaling of rényi entanglement entropies of the free fermi-gas ground state: A rigorous proof,” Phys. Rev. Lett. 112, 160403 (2014).
  * Widom (1982) Harold Widom, “On a class of integral operators with discontinuous symbol,” in  _Toeplitz Centennial: Toeplitz Memorial Conference in Operator Theory, Dedicated to the 100th Anniversary of the Birth of Otto Toeplitz, Tel Aviv, May 11–15, 1981_, edited by I. Gohberg (Birkhäuser Basel, Basel, 1982) pp. 477–500.
  * Shao  _et al._ (2015) Junping Shao, Eun-Ah Kim, F. D. M. Haldane, and Edward H. Rezayi, “Entanglement entropy of the $\nu=1/2$ composite fermion non-fermi liquid state,” Phys. Rev. Lett. 114, 206402 (2015).
  * Hu  _et al._ (2020) Wen-Jun Hu, Yi Zhang, Andriy H. Nevidomskyy, Elbio Dagotto, Qimiao Si, and Hsin-Hua Lai, “Fractionalized excitations revealed by entanglement entropy,” Phys. Rev. Lett. 124, 237201 (2020).
  * Melko  _et al._ (2010) Roger G. Melko, Ann B. Kallin, and Matthew B. Hastings, “Finite-size scaling of mutual information in monte carlo simulations: Application to the spin-$\frac{1}{2}$ $xxz$ model,” Phys. Rev. B 82, 100409 (2010).
  * Singh  _et al._ (2011) Rajiv R. P. Singh, Matthew B. Hastings, Ann B. Kallin, and Roger G. Melko, “Finite-temperature critical behavior of mutual information,” Phys. Rev. Lett. 106, 135701 (2011).
  * Iaconis  _et al._ (2013) Jason Iaconis, Stephen Inglis, Ann B. Kallin, and Roger G. Melko, “Detecting classical phase transitions with renyi mutual information,” Phys. Rev. B 87, 195134 (2013).
  * Stéphan  _et al._ (2014) Jean-Marie Stéphan, Stephen Inglis, Paul Fendley, and Roger G. Melko, “Geometric mutual information at classical critical points,” Phys. Rev. Lett. 112, 127204 (2014).
  * Kotliar  _et al._ (2000) G. Kotliar, E. Lange, and M. J. Rozenberg, “Landau theory of the finite temperature mott transition,” Phys. Rev. Lett. 84, 5180–5183 (2000).
  * Flynn  _et al._ (2023) Michael O. Flynn, Long-Hin Tang, Anushya Chandran, and Chris R. Laumann, “Momentum space entanglement of interacting fermions,” Phys. Rev. B 107, L081109 (2023).
  * Fisher  _et al._ (2022) Matthew P. A. Fisher, Vedika Khemani, Adam Nahum, and Sagar Vijay, “Random Quantum Circuits,” arXiv e-prints , arXiv:2207.14280 (2022), arXiv:2207.14280 [quant-ph] .
  * Calabrese  _et al._ (2010) Pasquale Calabrese, Massimo Campostrini, Fabian Essler, and Bernard Nienhuis, “Parity effects in the scaling of block entanglement in gapless spin chains,” Phys. Rev. Lett. 104, 095701 (2010).
  * Swingle  _et al._ (2013) Brian Swingle, Jeremy McMinis, and Norm M. Tubman, “Oscillating terms in the renyi entropy of fermi gases and liquids,” Phys. Rev. B 87, 235112 (2013).