††thanks: These two authors contributed equally to this work††thanks: These
two authors contributed equally to this work

# Influence functional of many-body systems:
temporal entanglement and matrix-product state representation

Michael Sonner  Department of Theoretical Physics, University of Geneva, Quai
Ernest-Ansermet 30, 1205 Geneva, Switzerland  Alessio Lerose  Department of
Theoretical Physics, University of Geneva, Quai Ernest-Ansermet 30, 1205
Geneva, Switzerland  Dmitry A. Abanin  Department of Theoretical Physics,
University of Geneva, Quai Ernest-Ansermet 30, 1205 Geneva, Switzerland

(August 27, 2024; August 27, 2024)

###### Abstract

Feynman-Vernon influence functional (IF) was originally introduced to describe
the effect of a quantum environment on the dynamics of an open quantum system.
We apply the IF approach to describe quantum many-body dynamics in isolated
spin systems, viewing the system as an environment for its local subsystems.
While the IF can be computed exactly only in certain many-body models, it
generally satisfies a self-consistency equation, provided the system, or an
ensemble of systems, are translationally invariant. We view the IF as a
fictitious wavefunction in the temporal domain, and approximate it using
matrix-product states (MPS). This approach is efficient provided the temporal
entanglement of the IF is sufficiently low. We illustrate the broad
applicability of the IF approach by analyzing several models that exhibit a
range of dynamical behaviors, from thermalizing to many-body localized. In
particular, we study the non-equilibrium dynamics in the quantum Ising model
in both Floquet and Hamiltonian settings. We find that temporal entanglement
entropy may be significantly lower compared to the spatial entanglement and
analyze the IF in the continuous-time limit. We simulate the thermodynamic-
limit evolution of local observables in various regimes, including the
relaxation of impurities embedded in an infinite-temperature chain, and the
long-lived oscillatory dynamics of the magnetization associated with the
confinement of excitations. Furthermore, by incorporating disorder-averaging
into the formalism, we analyze discrete time-crystalline response using the IF
of a bond-disordered kicked Ising chain. In this case, we find that the
temporal entanglement entropy scales logarithmically with evolution time. The
IF approach offers a new lens on many-body non-equilibrium phenomena, both in
ergodic and non-ergodic regimes, connecting the theory of open quantum systems
theory to quantum statistical physics.

##  I Introduction

The problem of a quantum system interacting with an environment has played an
important role since the early days of quantum mechanics Breuer and
Petruccione (2007). Describing dynamics of such open quantum systems is
essential for understanding diverse phenomena including the process of quantum
measurement, transport in nanostructures, and thermalization Weiss (2012).

It is an extremely challenging task to describe sufficiently complex,
realistic environments exactly. Therefore much work focused on studying
simplified models of environments, and on developing approximations for open
system dynamics. A model that played a special role is that of an environment
made of harmonic oscillators Feynman and Vernon (1963) which do not interact
among themselves. Interestingly, the dynamics of a seemingly simple open
quantum system – a two-level system interacting with such a bosonic
environment (or a bath) – exhibits a rich variety of regimes, depending on the
spectral density of the oscillators Leggett  _et al._ (1987).

More recently, experimental advances have enabled realization of many-body
quantum systems which are well isolated from an external environment Bloch
_et al._ (2008); Blatt and Roos (2012). Such setups have brought into focus
the problem of highly non-equilibrium dynamics in closed many-body systems.
Various universality classes, distinguished by their drastically different
dynamical behavior, have been discovered.

The most common class is that of thermalizing, or ergodic dynamics: starting
from generic initial states, local subsystems reach an effectively thermal
state under unitary dynamics D’Alessio  _et al._ (2016). Remarkably,
ergodicity can be broken by disorder-induced many-body localization (MBL)
Nandkishore and Huse (2015); Abanin  _et al._ (2019); Alet and Laflorencie
(2018), as first foreseen by Anderson in the paper reporting discovery of
single-particle localization Anderson (1958). Persistence of MBL in
periodically driven (Floquet) systems Ponte  _et al._ (2015); Lazarides  _et
al._ (2015); Abanin  _et al._ (2016); Bordia  _et al._ (2017) has opened the
door to novel non-equilibrium phases such as discrete time crystals Khemani
_et al._ (2016); Else  _et al._ (2016) and anomalous Floquet insulators Nathan
_et al._ (2019). Much attention has also focused on investigating dynamics in
integrable one-dimensional systems, which at long times relax to a generalized
Gibbs ensemble Essler and Fagotti (2016), and on weak ergodicity breaking
mechanisms such as quantum scars Serbyn  _et al._ (2020).

At an intuitive level, the character of the system’s dynamics is determined by
its properties as a quantum environment. An ergodic system acts as a “good”
thermal bath for its subsystems, while MBL and other non-ergodic systems fail
to do so. This observation links the dynamics of a many-body system to that of
an open quantum system; however, the environment is now itself a complex
interacting system, rather than a set of non-interacting harmonic oscillators.

The ability of a many-body system to act as a thermal bath has been
characterized using a number of eigenstate and dynamical probes Abanin  _et
al._ (2019); Alet and Laflorencie (2018); D’Alessio  _et al._ (2016). The
eigenstate thermalization hypothesis (ETH) Deutsch (1991); Srednicki (1994),
as well as its breakdown in MBL systems provide a particularly useful tool.
One may also directly study the dynamics of physical observables following a
quantum quench using numerical techniques such as exact diagonalization and
tensor networks Verstraete  _et al._ (2008); Paeckel  _et al._ (2019), and
verify whether they settle to thermal values. The quantum quench setup is
routinely employed experimentally in different platforms, including cold
atoms, trapped ions, and superconducting qubits.

In this paper, building on very recent ideas Lerose  _et al._ (2021); Ye and
Kin-Lic Chan (2021) and earlier related work Bañuls  _et al._ (2009), we
employ the Feynman-Vernon influence functional to analyze the problem of
quantum dynamics and characterize the properties of a many-body system as an
environment. The influence functional (IF) for an open quantum system is
obtained by performing the Keldysh path integral over the environment degrees
of freedom $\\{q_{j}(\tau),\bar{q}_{j}(\tau)\\}$, subject to a time-dependent
trajectory of the system $Q(\tau),\bar{Q}(\tau)$ (here and below each
trajectory has a forward and a backward part, and coordinates without/with a
bar parametrize the former/latter). The idea of this approach is illustrated
in Fig. 1a: the effect of the environment is completely described by the IF
$\mathcal{I}[Q(\tau),\bar{Q}(\tau)]$, which weighs the trajectories of the
system. The explicit expression for the IF, available for certain simplified
models such as a bath of harmonic oscillators, has provided a starting point
for describing a two-level system subject to dissipation Leggett  _et al._
(1987).

---
Figure 1: a) Idea of the Feynman-Vernon influence functional for a quantum
system interacting with a bath, which consists of non-interacting harmonic
oscillators. Considering the time evolution in the Keldysh path integral
representation, the bath degrees of freedom
$\\{q_{j}(\tau),\bar{q}_{j}(\tau)\\}$ are integrated out. This yields a
contribution $\mathcal{I}[Q(\tau),\bar{Q}(\tau)]$ in the path integral over
the system coordinates only, which is the influence functional. b) A tensor
network describing evolution of a spin chain, in the discrete path integral
representation. Time runs upwards and the foreground (background) represents
the forward (backward) branch of time evolution. Contraction of the parts of
the network contained in the blue-shaded boxes defines the left and right
influence matrices – tensors acting on the space of trajectories of spin
$j=0$. For the case of Hamiltonian evolution, when the tensor network
represents a discretization of the continuous time evolution, the influence
matrices approximate the corresponding continuous-time influence functionals.
c) The influence functional of a homogeneous spin system is an eigenvector of
the dual transfer-matrix represented by the red-shaded box; equivalently, it
satisfies a linear self-consistency equation (7).

The approach proposed in Ref. Lerose  _et al._ (2021) is based on the
observation that the IF of a spatially homogeneous many-body system, while
difficult to compute directly, satisfies a self-consistency equation. In the
tensor-network language, this equation amounts to finding an eigenvector of
the dual transfer matrix arising in the path-integral description of the
system’s dynamics Bañuls  _et al._ (2009), as schematically illustrated in
Fig. 1b,c. This eigenvector can be efficiently sought in the form of a matrix-
product state (MPS), provided its temporal entanglement remains low. This is
analogous to how MPS-based algorithms such as the density matrix
renormalization group rely on the spatial entanglement of ground states being
low Verstraete  _et al._ (2008). However, while the latter is guaranteed by
rigorous results concerning area-law entanglement scaling in gapped systems,
relatively little is known about the scaling of temporal entanglement.

It was appreciated early on Müller-Hermes  _et al._ (2012) that for certain
dynamical problems the approaches based on transverse contractions of the
time-evolution tensor network be computationally more favorable compared to
the conventional methods such as the time-evolving block decimation (TEBD).
The recent work Lerose  _et al._ (2021), which considered periodically driven
(Floquet) systems, found that in a family of thermalizing models, including
so-called dual-unitary circuits Akila  _et al._ (2016); Bertini  _et al._
(2019); Piroli  _et al._ (2020), the discrete-time IF – or influence matrix
(IM) – can be found exactly and has vanishing temporal entanglement.
Furthermore, Ref. Klobas  _et al._ (2020) found a different exact MPS solution
for an integrable quantum cellular automaton. The IF approach to a Floquet-MBL
model was shown to be efficient as well Sonner  _et al._ (2020). In a
different direction, Refs. Strathearn  _et al._ (2018); Cygorek  _et al._
(2021); Ye and Kin-Lic Chan (2021) proposed a tensor-network-based compression
of the IF of baths of noninteracting particles.

In this paper, we extend the IF approach and demonstrate its utility for
several classes of non-equilibrium problems of current interest, both in
ergodic and non-ergodic regimes. We first employ it to describe thermalization
dynamics in Floquet and Hamitonian many-body systems. Focusing on spin-$1/2$
chains, we apply the IF formalism to compute the behavior of dynamical
correlators in an infinite-temperature ensemble. Furthermore, we analyze the
time-dependent observables following a quantum quench starting from simple
pure initial states. In particular, we simulate the long-lived oscillations of
local observables that arise due to confinement of excitations in the quantum
Ising chain in a tilted field Kormos  _et al._ (2017); James  _et al._ (2019);
Mazza  _et al._ (2019); Lerose  _et al._ (2020).

In a different direction, below we adapt the IF approach to analyze the
dynamics of non-equilibrium, non-ergodic phases. We focus on a one-dimensional
discrete time-crystal (DTC) enabled by disorder-induced MBL. Our approach
allows us to perform exact averaging over the disordered spin-spin couplings
(see below). We then consider the self-consistency equation for the
disordered-averaged IF, and solve it using an MPS representation.

In this paper we will be particularly interested in the scaling of temporal
entanglement (TE) entropy with the evolution time. We find that in all
problems considered, TE remains relatively low up to sizable evolution times,
which indicates efficiency of the MPS representation of the IF. Furthermore,
we investigate the behavior of TE in Hamiltonian systems as the continuum
limit is taken by decreasing the time step in the Trotterized evolution
operator. Surprisingly, TE entropy decreases to zero in the limit of vanishing
time steps; however, as we argue below, small singular values may not always
be truncated to faithfully approximate the IF of a Hamiltonian system.

The rest of the paper is organized as follows. In the subsequent Section II,
we introduce the IF approach for both Floquet and Hamiltonian dynamics. We
also discuss the MPS representation and the method we use to numerically
compute the IF. Subsequently, in Section III we analyze thermalization
dynamics in a (non-integrable) quantum Ising model, presenting results for
quenches from different initial ensembles/states. Both Floquet and Hamiltonian
dynamics will be considered. Further, in Section IV we turn to non-ergodic
dynamics and discuss how to incorporate disorder-averaging into the IF
formalism, considering discrete time crystal (DTC) as an application. We
conclude in Section V by discussing future directions and potential of the IF
approach to quantum many-body dynamics and open quantum systems.

##  II Influence functional for many-body dynamics

We start our analysis by formulating the influence functional approach to
describing dynamics of a many-body system. Throughout the paper, we focus on
one-dimensional lattice systems of spins with local Hilbert space dimension
$q$ subject to local two-body interactions. We note that the extensions to
bosons/fermions, longer-range $k$-body interactions, and to higher-dimensional
systems are straightforward. We consider models in which the time evolution
can be represented by a general “brick-wall” quantum circuit (Fig. 1b). This
class includes Floquet models, and trotterized Hamiltonian evolution. We note
that Ref. Lerose  _et al._ (2021) introduced the IF approach for a narrower
class of Floquet models.

The central idea of the approach is summarized in Fig. 1b,c. The evolution
operator over one period of the brick-wall circuit is given by:

$U=U_{e}U_{o},\;\;U_{o}=\prod_{j}U_{2j-1,2j},\;\;U_{e}=\prod_{j}U_{2j,2j+1},$
(1)

where the two operators correspond to odd and even layers of the quantum
circuit. The two-body evolution operators (unitary gates) in the above
expression define the model in the case of Floquet dynamics. For the case of
Hamiltonian evolution, generated by a two-body local Hamiltonian,

$H=\sum_{i}H_{i,i+1},$ (2)

the gates are obtained by trotterizing the time evolution, with a time step
$\epsilon\to 0$:

$U_{i,i+1}=e^{-iH_{i,i+1}\epsilon}.$ (3)

The physical evolution time $t$ in the Hamiltonian case is the product of the
time step $\epsilon$ and the number of time steps $T$ in the tensor network in
Fig. 1b, i.e., $t=T\epsilon$. We note that the brick-wall quantum circuit can
also be used to encode higher-order Trotter schemes; here for simplicity we
will consider only the lowest-order scheme.

The time-evolved density matrix of the system after $T$ steps of evolution is
related to the initial state density matrix $\rho^{0}$ as follows:

$\rho^{T}=U^{T}\rho^{0}{U^{\dagger}}^{T}.$ (4)

For simplicity, we will consider initial density matrices which can be
represented as a tensor product of individual spin density matrices, i.e.,

$\rho^{0}=\bigotimes_{j}\rho_{j}^{0}\,.$ (5)

Rewritten via a discrete path integral over intermediate configurations of the
system, the time evolution in Eq. (4) can be represented by a tensor network,
as illustrated in Fig. 1b. We denote the forward and backward trajectory of
spin $j$ by $\sigma_{j}^{\tau},\bar{\sigma}_{j}^{\tau}$,
$s_{j}^{\tau},\bar{s}_{j}^{\tau}$, as shown in the Figure. The gates acting in
the forward part of the network are $U_{j,j+1}$ and the ones in the backward
part of the network are $U^{*}_{j,j+1}$.

We will be interested in describing the evolution of a sufficiently small
subsystem, e.g. a single spin $j=0$. The dynamics can be characterized by
correlation functions of the form

$\displaystyle\langle\hat{O}_{0}(T)\hat{O}_{0}(0)\rangle,$ (6)

where $\hat{O}_{0}$ is an operator representing a local observable at site
$j=0$, see Fig. 1b. To calculate such correlators we begin by tracing out the
environment degrees of freedom, i.e., the other spins; the spins with $j>0$
can be viewed as a “right” environment, and the ones with $j<0$ as a “left”
environment. This procedure yields right and left influence matrices
$\mathcal{I}_{R(L)}$, which are formally introduced by contracting tensor
networks contained in the blue-shaded boxes in Fig. 1b, with a boundary
defined by spin $j=0$ trajectory; in this contraction, the summation over all
final states of the spins in the environment is carried out. We can now
calculate any correlator of the form Eq. (6) by contracting the IMs with a
tensor network for evolution of spin $j=0$, with operators $\hat{O}_{0}$
inserted at $\tau=0,\,T$. In the Hamiltonian case, the IMs approximate the
corresponding continuous-time IFs as $\epsilon\to 0$. (In the following we
will often use the names IM and IF interchangeably.)

In certain cases, which include dual-unitary circuits/perfect dephasers Akila
_et al._ (2016); Bertini  _et al._ (2019); Piroli  _et al._ (2020); Lerose
_et al._ (2021) and certain integrable systems Klobas  _et al._ (2020); Lerose
_et al._ (2021), the contraction of the tensor network that defines the IF can
be carried out analytically. Generally, the IF for an environment with two
additional spins can be obtained from the previous IF by applying one vertical
slice of the tensor network in Fig. 1c (red shaded box). These vertical slices
define the dual transfer matrices Bañuls  _et al._ (2009), which act on the
$q^{4T}$-dimensional space of single-spin trajectories. As discussed below,
for certain systems contracting the tensor network in the space direction
using dual transfer matrices proves advantageous compared to the standard
contraction in time direction.

In translationally invariant systems, the dual transfer matrix $\tilde{T}$ is
the same at every step. Notice that the strict light cone in our model means
that the IF that describes evolution over time $T$ can depend at most on the
motion of closest $2T$ spins. Hence the IF of a homogeneous system is a unique
eigenvector of $\tilde{T}$ with an eigenvalue $1$. The eigenvalue equation for
the right IF

$\tilde{\mathcal{T}}|\mathcal{I}\rangle=|\mathcal{I}\rangle,$ (7)

(where we omitted the subscript $R$ for simplicity) has a simple physical
interpretation: it can be viewed as a self-consistency equation for the IF. In
essence, a spin evolving under the influence of the right environment, exerts
the same influence on its left neighbor. The left influence functional clearly
satisfies an analogous equation with a slightly modified transfer matrix.

We note that this approach can be extended to describe ensembles of disordered
models with spatially uncorrelated, translationally-invariant distribution: in
that case, the transfer matrix in Eq. (7) is disorder-averaged, and its
eigenvector corresponds to an IF of a disorder-averaged environment. In Ref.
Sonner  _et al._ (2020), a particular Floquet-MBL model with on-site disorder
fields was considered (see also Refs. Sünderhauf  _et al._ (2018); Chan  _et
al._ (2020); Garratt and Chalker (2020)). Below in Section IV we will discuss
a generalization of the method to a model with disordered two-body
interactions, which exhibits discrete time crystal behavior.

A promising approach to solving the self-consistency equation (7) is based on
approximating the influence matrix by an MPS with a fixed bond dimension
Bañuls  _et al._ (2009); Müller-Hermes  _et al._ (2012); Lerose  _et al._
(2021), similarly to conventional techniques for two-dimensional systems
Haegeman and Verstraete (2017). This approach is expected to be efficient when
TE entropy is low. In that case, the IF formulation may allow one to analyze
time-dependent correlation functions which are out of reach of standard
methods such as TEBD, exact diagonalization, etc.

Interestingly, temporal entanglement can be low in different dynamical
regimes. A striking example is provided by perfect-dephaser Floquet circuits:
for an infinite-temperature initial density matrix and certain pure initial
states in MPS form, the IF is a product state with zero TE Lerose  _et al._
(2021); Piroli  _et al._ (2020); in this case, the environment effectively
dephases local spins at every step of Floquet evolution. TE remains low in the
vicinity of such solvable points, which ensures the efficiency of the MPS
representation. Furthermore, the TE was found to be low up to long times in
certain non-integrable Hamiltonian models Müller-Hermes  _et al._ (2012);
Hastings and Mahajan (2015), and in proximity to integrability Lerose  _et
al._ (2021). Ref. Sonner  _et al._ (2020) further demonstrated that the
approach remains efficient in the ergodicity-breaking Floquet-MBL phase,
thanks to slow, sub-linear scaling of TE of the disorder-averaged IF. Despite
these results, the full potential of the IF-based approach to describe quantum
dynamics remains to be investigated; the key open question concerns the
scaling of TE with the evolution time and the formulation of suitable
numerical algorithms to approximate it.

In what follows, we will analyze several Hamiltonian and Floquet models, both
in the ergodic and non-ergodic regimes. We will compute dynamical correlation
functions using the IF approach for both mixed and pure initial states of the
system. We will also investigate the scaling of TE, and its dependence on the
system’s parameters, and on the time evolution step $\epsilon$, for the case
of Hamiltonian dynamics.

##  III Thermalizing dynamics in Floquet and Hamiltonian systems

In this Section, we will study thermalization dynamics in non-disordered
Floquet and Hamiltonian models in different setups using the IF approach. We
will consider relaxation of observables in a chaotic system initialized in an
infinite-temperature ensemble, both for a spin that is part of the system, and
for an impurity embedded in a many-body environment. Further, quench dynamics
from a pure initial state in a model with confined excitations will be
studied. For all examples considered, we will investigate the behavior of TE
entropy with evolution time.

###  III.1 Floquet dynamics

Figure 2: Simulation of the kicked Ising model Eq.(8) with parameters
$J=0.8,g=0.7236,h=0.6472$ Kim  _et al._ (2014) for an infinite-temperature
initial state. Left: Autocorrelation function of the $\sigma^{z}$ operator as
a function of time. Exponential relaxation signals chaotic dynamics. The blue
crosses correspond to exact diagonalization results, which coincide with the
MPS calculations. The results are well converged up to time $T\approx 20$.
Right: Temporal entanglement entropy as a function of evolution time for
different bond dimensions. Dashed lines are maximal entanglement entropies
which are encountered during iteration process starting from open boundary IF.
Starting from the PD IF, the entropy of the final IF is also the maximal
entropy (solid lines). This illustrates the importance of the choice of the
initial IF in the context of the iterative MPS method.

We start by illustrating the use of the method for Floquet systems. As an
example of a thermalizing Floquet model, we choose the kicked Ising model for
spins-$1/2$, which has been extensively studied in the literature Kim  _et
al._ (2014); Akila  _et al._ (2016); Bertini  _et al._ (2019). The Floquet
operator of this model is given by:

$U=e^{-ig\sum_{j}\sigma_{j}^{x}}\,e^{-iJ\sum_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}-ih\sum_{j}\sigma_{j}^{z}},$
(8)

where $\sigma_{j}^{\alpha}$, $\alpha=x,y,z$ are standard spin-$1/2$ Pauli
operators. Depending on the values of parameters $J,g,h$, the model displays a
rich variety of dynamical behavior. For $h=0$, it can be mapped onto a free-
fermion model using a Jordan-Wigner transformation, and is therefore
integrable. In this case, it can be shown that the IF obeys an area-law
temporal entanglement, even in the limit of very long evolution times,
$T\to\infty$ Lerose  _et al._ (2021).

Furthermore, the model has special properties at the self-dual points
$|J|=|g|=\pi/4$ Akila  _et al._ (2016); Bertini  _et al._ (2019). In that
case, even though dynamics at $h\neq 0$ are chaotic, many quantities can be
computed analytically, including various dynamical correlation functions. At
these points, the IF has a perfect-dephaser form Lerose  _et al._ (2021). When
parameters of the model are detuned slightly away from these PD points, the
temporal entanglement remains low, showing an apparent volume-law scaling with
$T$, albeit with a small prefactor; thus, the IF can be well-approximated by a
MPS with a moderate bond dimension Lerose  _et al._ (2021).

Here, we focus on the parameter choice $J=0.8,g=0.7236,h=0.6472$, considered
in Ref. Kim  _et al._ (2014). It was found that such a model displays chaotic
dynamics, and, moreover, the deviations from the eigenstate thermalization
hypothesis quickly become negligible already at moderate finite sizes. We
computed the IF for the infinite-temperature initial density matrix,

$\rho_{j}^{0}=\frac{1}{2}\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)$ (9)

by iteratively applying the dual transfer matrix to an initial “boundary” IF.
As discussed below, certain choices of the boundary IF can significantly
improve the performance of the MPS approach.

To perform the tensor network contraction, we set up a code using the tenpy
library Hauschild and Pollmann (2018). The diagonal form of interactions in
the $z$-basis allows us to identify the input and output legs at each time
step (denoted by $s$ and $\sigma$ in Fig. 1). Thus, compared to the general
circuits, the number of degrees of freedom in the dual transfer matrix and in
the IM is reduced by a factor of $2$ Lerose  _et al._ (2021); moreover, the
maximum velocity of propagation is $1$ rather than $2$. The absence of
correlations in the initial state of the chain combined with the strict light
cone ensures that the IF reaches its self-consistent thermodynamic limit form
after at most $T$ iterations of the dual transfer matrix $\tilde{T}$. The
latter can be straightforwardly represented as an MPO of bond dimension $4$
for the model (8). We select a maximum bond dimension $\chi$ and iteratively
apply this MPO to a boundary vector expressed in an MPS form, truncating the
result to $\chi$ using SVD compression when necessary.

The choice of the boundary IF does not affect the final IF obtained after $T$
iterations, provided the calculations are carried out exactly. However,
different choices generally lead to different intermediate IFs encountered
during iterations. Physically, the $\ell$-th intermediate IF describes the
effect of an environment that consists of $\ell$ spins, which evolve under
unitary evolution set by the brick-wall quantum circuit, and are subject to an
external bath at the boundary; the properties of this bath are determined by
the chosen initial IF. In order for MPS-based methods to be efficient, it is
important that the TE of the intermediate IF remains low enough, such that
truncation does not incur a significant error. However, this is not always the
case. Interestingly, seemingly the most natural choice of the boundary IF that
corresponds to open boundary conditions,
$\mathcal{I}[\sigma_{\tau},\bar{\sigma}_{\tau}]=1$, turns out not to be the
most efficient one: in this case TE grows as a function of $\ell$ first before
decreasing again to the final value. We found that the maximal TE entropy of
intermediate IFs can be very high; its growth as a function of evolution time
is shown by the dashed lines in Fig. 2b. We note that a similar observation
was reported in Ref. Ye and Kin-Lic Chan (2021) for different models.

However, the problem of the “entanglement barrier” described above can be
circumvented by choosing a different boundary IF. For the thermalizing Floquet
model, we chose the perfect dephaser IF
$I[\sigma_{\tau},\bar{\sigma}_{\tau}]=\prod\delta_{\sigma_{\tau},\bar{\sigma}_{\tau}}$
as a boundary condition. Physically, this corresponds to a boundary where the
spin is measured (in the $\sigma^{z}$-basis) every period. We confirmed that
this choice indeed leads to a monotonic growth of entanglement as a function
of iteration number $\ell$, and therefore the “entanglement barrier” can be
avoided.

Using the MPS representation of the IF we computed the dynamical correlation
function:

$\mathcal{C}_{zz}(t)=\Big{\langle}\hat{\sigma}^{z}_{0}(t)\hat{\sigma}^{z}_{0}(0)\Big{\rangle}\,.$
(10)

The resulting exponential decay of this correlator, illustrated in Fig. 2a, is
consistent with the scrambling of local quantum information expected in this
chaotic system.

Fig. 2b shows the scaling of TE entropy (for a bipartition at $\tau=T/2$) with
the evolution time computed with the MPS approach for three choices of the
maximum bond dimension. The convergence of TE for different bond dimensions up
to $T\approx 20$ indicates that the MPS form provides a good approximation to
the exact IF. For times up to $T=7$ this was verified by comparing the results
to exact ones found with the standard sparse linear algebra library ARPACK
(blue crosses in Fig. 2). TE entropy follows an apparent volume-law scaling in
the accessible time interval. We note that at present it is unclear whether
this scaling persists as $T\to\infty$.

###  III.2 Hamiltonian dynamics

Figure 3: Left: Autocorrelation $\mathcal{C}_{zz}$ of Hamiltonian model (11)
for parameters $J=1.0,g=\sqrt{2},h=0.681$ at infinite temperature with
different impurities. The parameters were chosen to be incomensurate to avoid
any resonances which may affect dynamical properties. The result for a
homogeneous chain is plotted by the solid line. We can observe the
perturbative (Born) limit by tuning the coupling of the impurity spin to the
spin chain (dotted lines) and the memoryless (Markovian) limit by tuning the
driving of the impurity (dashed lines) as well as the combined Born-Markov
limit (dash-dotted line). The calculation is fully converged for bond
dimensions $\chi=64,128$ and time steps $\epsilon=0.01,0.02,0.04$. Right: TE
entropy scales to zero as a function of the time step $\epsilon$, in a manner
that is approximately consistent with a power-law dependence on $\epsilon$ in
the parameter range considered.

We can use a related approach to study Hamiltonian dynamics, by viewing Eq.
(8) as a (second-order) trotterization of the continuous time unitary flow
$U=e^{-iHt}$ generated by the Hamiltonian

$H=-J\sum_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}-h\sum_{j}\sigma_{j}^{z}-g\sum_{j}\sigma_{j}^{x}.$
(11)

This Hamiltonian describes a quantum Ising chain in a tilted magnetic field.
The identification is realized by rescaling the parameters in the Hamiltonian
Eq. (11) with the time step $0<\epsilon\ll 1$ in order to obtain the
dimensionless parameters $J,g,h\mapsto J\epsilon,g\epsilon,h\epsilon$ in Eq.
(8); in this case, the physical time is $t=T\epsilon$. Higher-order Trotter
schemes can be implemented analogously.

The model in Eq. (11), despite its apparent simplicity, exhibits a rich
variety of phenomena, and therefore it provides a fruitful playground in
quantum many-body physics, both in and out-of-equilibrium. Similarly to its
Floquet counterpart, it is integrable when the magnetic field is purely
transverse ($h=0$), and exhibits a paradigmatic quantum phase transition
between a paramagnetic $|g|>J$ and a ferromagnetic $|g|<J$ phase,
corresponding to a spontaneous breaking of the $\mathbb{Z}_{2}$ symmetry. The
free fermionic quasiparticles of the model change their nature across the
transition, from non-topological, spin-wave-like to topological, domain-wall-
like excitations interpolating between the two degenerate vacua Sachdev
(2011). In the latter case, the longitudinal field $h$ explicitly breaks the
symmetry, lifting the degeneracy of the two magnetized ground states and
inducing a first-order transition. Accordingly, a confining potential
$V(r)\propto hr$ arises between a pair of domain walls, which bind together
into “mesonic” composite particles McCoy and Wu (1978). The model with
nonvanishing fields $g,h\neq 0$ is non-integrable and its dynamics are
believed to be generally chaotic Kim  _et al._ (2014), although anomalously
slow or even suppressed relaxation after a quench has been reported even far
away from solvable limits Bañuls  _et al._ (2011); Lin and Motrunich (2017);
Wurtz and Polkovnikov (2020) and related to confinement of excitations Kormos
_et al._ (2017); James  _et al._ (2019); Mazza  _et al._ (2019); Lerose  _et
al._ (2020).

First, we study the local relaxation of a slower impurity spin at position
$j=0$ embedded in the quantum Ising chain at infinite temperature. This
impurity is characterized by having smaller values of the fields $g_{0}=\alpha
g$, $h_{0}=\alpha h$ and of the couplings $J_{-1,0}=J_{0,1}=\beta J$ to the
neighboring spins. Note that the IF only depends on $\beta$ and not on
$\alpha$. Tuning $\beta$ and $\alpha$ to be much smaller than $1$, the
impurity dynamics approaches a typical open quantum system (OQS) setup, i.e.,
that of a two-level system weakly coupled to a much faster environment, which
allows for the celebrated Born (for $\beta\ll 1$) and Markov (for $\alpha\ll
1$) approximations, respectively. Putting $\beta$ and $\alpha$ to $1$,
instead, we retrieve the standard homogeneous quantum many-body system setup,
where the distinction between the relevant subsystem and its bath is purely
conceptual. Fig. 3 illustrates representative results of these computations.
The left panel shows the infinite-temperature time-dependent autocorrelator
$\mathcal{C}_{zz}(t)$ [cf. Eq. (10)] converged with respect to both Trotter
step and bond dimension, for decreasing $\alpha$ and $\beta$ from the
homogeneous case to an OQS limit. The right panel shows the scaling of TE
entropy $S$ converged with respect to the bond dimension vs the physical
evolution time $t$ (i.e., for $T=t/\epsilon$ steps), for a range of decreasing
Trotter steps $\epsilon$. Interestingly, for fixed $t$, the value of $S$
scales to zero for $\epsilon\to 0$ with an apparent power law with an exponent
$\simeq 1.8$. We will further comment on this observation below.

Figure 4: Comparison of TEBD and IM methods for the Hamiltonian model (11) at
$J=1.0,g=0.25,h=0.4$ undergoing a quench from the polarized state (cf Refs.
Bañuls  _et al._ (2011); Kormos  _et al._ (2017)). In both cases, bond
dimension $\chi=64$ was used. Left: While both curves are converged with
respect to time step and bond dimension, the TEBD results still show small
finite size effects. Right: The real space entanglement entropy of the TEBD
MPS is significantly higher than the TE entropy for all three time steps,
indicating the lower computational cost of the IM approach.

Second, we study a prototypical quantum quench setup, where the (now impurity-
free) system is initialized in a fully polarized state in the positive $z$
direction. We focus on the ferromagnetic phase, setting $g=0.25J$ and
$h=0.4J$. The evolution of the ferromagnetic order parameter
$\langle\sigma^{z}_{0}(t)\rangle$ in Fig. 4 displays clear, long-lived
quasiperiodic oscillations, consistent with the discrete tower of mesonic
particles excited by the sudden quench Kormos  _et al._ (2017). The right
panel shows that the same scaling of the maximum TE entropy with the Trotter
step is found even for quench dynamics, cf. Fig. 3. In this case, furthermore,
it makes sense to compare spatial entanglement entropy to temporal
entanglement entropy. As the plot clearly shows, the former – obtained by TEBD
converged with system size and bond dimension – significantly exceeds the
latter, while still remaining relatively low as a consequence of confinement
Kormos  _et al._ (2017); Lerose  _et al._ (2020).

###  III.3 Continuous-time limit

The simulations presented in this Section were obtained with a second-order
Trotter scheme with time steps $\epsilon=0.01,0.02,0.04$. The Trotter error
was found to be negligible in all cases for $\epsilon\leq 0.04$: in the
examples shown in Figs. 3, 4, time-dependent observables were perfectly
converged with respect to decreasing $\epsilon$. Thus, effectively, the
continuous time limit is well captured by the discretized evolution.

On the methodogical side, it is interesting to discuss the non-trivial scaling
of temporal entanglement when $\epsilon$ is decreased. Surprisingly, as shown
in the right panel of Fig. 3, the maximum temporal entanglement entropy
$S_{t}$ at fixed physical evolution time $t$ scales to zero for $\epsilon\to
0$. This behavior was found in all our simulations of Hamiltonian dynamics. To
gain an insight into the origin of this behavior, below we analyze several
solvable limits.

First, we study the continuous-time limit of the IM in the limiting case
$g=0$, where quantum fluctuations are completely suppressed and all product
states in the $\sigma^{z}$-basis are eigenstates of the Floquet
operator/Hamiltonian for all $J$, $h$. Let us first consider the Floquet
problem of Eq. (8), with an infinite-temperature initial ensemble. The
analytical form of the IM can be found exactly:

$I[\sigma_{\tau},\bar{\sigma}_{\tau}]=\cos\Big{[}J\sum_{\tau}(\sigma_{\tau}-\bar{\sigma}_{\tau})\Big{]}.$
(12)

It is straightforward to verify that this expression satisfies the self-
consistency equation for $g=0$. Physically, this expression can be thought of
as the influence functional due to an interaction with a constant classical
magnetic field along $\hat{z}$-axis, which takes values $+J$ or $-J$ with
probability $1/2$. To compute the von Neumann entropy of $I$, we need to
normalize it as a pure wavefunction in the folded Hilbert space, i.e.,

$\ket{\Psi}=\frac{1}{\sqrt{2[1+\cos^{2T}(2J)]}}\bigg{(}\bigotimes_{\tau}\frac{1}{2}\begin{pmatrix}1\\\
1\end{pmatrix}_{cl,\tau}\otimes\begin{pmatrix}e^{2iJ}\\\
e^{-2iJ}\end{pmatrix}_{q,\tau}\;+\;\bigotimes_{\tau}\frac{1}{2}\begin{pmatrix}1\\\
1\end{pmatrix}_{cl,\tau}\otimes\begin{pmatrix}e^{-2iJ}\\\
e^{+2iJ}\end{pmatrix}_{q,\tau}\bigg{)}$ (13)

where we have factored the four-dimensional folded Hilbert space at time
$\tau$ into the tensor product of a “classical”
$\mathrm{Span}(\ket{\uparrow\uparrow},\ket{\downarrow\downarrow})$ and a
“quantum” $\mathrm{Span}(\ket{\uparrow\downarrow},\ket{\uparrow\downarrow})$
two-dimensional subsystems. The chosen prefactor ensures the normalization
condition, $\braket{\Psi}{\Psi}=1$.

Our aim is to compute the bipartite von Neumann entanglement entropy of this
normalized state wavefunction with respect to a bipartition $(1,M),(M+1,T)$.
First, we notice that the classical sector does not contribute to
entanglement, since the corresponding degrees of freedom are in a product
state. Second, restricted to the quantum sector, $\ket{\Psi}$ can be thought
of as an equal-weight quantum superposition of two pure spin-coherent states
of a spin of magnitude $T/2$, fully polarized along the directions
$\cos(2J)\hat{x}\pm\sin(2J)\hat{y}$, respectively. These two spin-coherent
states can be pictured as Gaussian wavepackets on the surface of the Bloch
sphere, centered around those two directions, with transverse fluctuations of
width equal to $1/\sqrt{T}$ (relative to the unit radius of the Bloch sphere).
Since the two polarization directions form an angle $4J$, it is easy to see
that the overlap of the two spin-coherent states decreases exponentially with
$J^{2}T$. For $J^{2}T\ll 1$, the two states have large overlap, so
$\ket{\Psi}$ is close to a product state, giving a vanishing von Neumann
entanglement entropy. In contrast, for $J^{2}T\gg 1$, the two states are
effectively orthogonal, such that $\ket{\Psi}$ becomes a GHZ state with von
Neumann entanglement entropy $\log 2$. The exact calculation gives the binary
entropy

$S_{M,T}(J)=-P\log P-(1-P)\log(1-P),$ (14)

with

$P=1-\frac{\cos^{2M}(2J)+\cos^{2(T-M)}(2J)}{1+\cos^{2T}(2J)}\,,$ (15)

which exhibits the expected features. Setting $x=J^{2}T$, $0\leq f=M/T\leq 1$,
this quantity becomes a smooth function of $x$ and $f$ in the limit
$T\to\infty$.

Figure 5: Behavior of the single-particle entanglement spectrum of the IM of
the transverse-field quantum Ising chain. The largest eigenvalues $0\leq
p_{n}\leq 1/2$ of the correlation matrix of a half-time interval originating
from the IM wavefunction $\ket{\Psi}$ (see the text) are reported, rescaled by
$\epsilon^{2}$, as a function of $\epsilon$, for three values of $g/J$
reported above the plots. In all cases, the results indicate the scaling
behavior $p_{n}\sim\epsilon^{2}P_{n}$, where $P_{n}$ are
$\epsilon$-independent quantities.

We are now in a position to understand the continuous-time limit in this
solvable case. Setting $J=\epsilon$, $T=t/\epsilon$, we find that the variable
$J^{2}T=\epsilon t$ flows to zero for any fixed value of the physical
evolution time $t$. Since $P\sim\epsilon^{2}$ as $\epsilon\to 0$, the von
Neumann entropy of any extensive bipartition converges to zero in the
continuous-time limit. It is instructive to note that in this case the
entanglement spectrum is composed of only two values, $P$ and $1-P$. When
$\epsilon$ is sufficiently small, one might be tempted to truncate the
corresponding singular value (e.g., when performing a standard SVD
compression). This seemingly harmless approximation, however, would produce
incorrect results. Indeed, after such a truncation, the IM becomes a product
state relative to the considered bipartition, say at the bond $(M,M+1)$:
physically, this corresponds to “refreshing the environment”, i.e., tracing it
out after $M$ time steps and replacing it with another infinite-temperature
environment before time step $M+1$. Performing a similar truncation in other
bonds, we end up with an IM corresponding to a “noisy” classical magnetic
field along $\hat{z}$ switching between $\pm J$, whereas the actual exact IM
in Eq. (12) represents the influence of a constant field. We thus conclude
that in the continuous-time limit, truncation of small Schmidt values might be
dangerous, and may not lead to an accurate representation of system’s
dynamics. Furthermore, counter-intuitively, the mere value of temporal entropy
is not representative of the extent of temporal correlations and of the
computational resources needed to parametrize the IM.

We have further confirmed this picture by an analytical calculation of the IF
in the integrable limit $h=0$ Lerose  _et al._ (2021). Making use of the
quadratic fermionic representation of the model, $S_{t}$ is reduced to the
entropy of a Bardeen-Cooper-Schrieffer-like wavefunction on the Keldysh
contour, due to its Gaussian form. In this case, as shown in Fig. 5, the
scaling of $S_{t}$ with $\epsilon$ stems from an $\epsilon^{2}$ global scaling
of the single-particle entanglement spectrum, $p_{n}\sim\epsilon^{2}P_{n}$,
which gives

$S_{\tau,t}\equiv-\sum_{n}p_{n}\log p_{n}+(1-p_{n})\log(1-p_{n})\\\
\sim\epsilon^{2}[\log(1/\epsilon^{2})+\hat{S}],$ (16)

where $\hat{S}=-\sum_{n}P_{n}\log P_{n}+(1-P_{n})\log(1-P_{n})$ is a finite,
$\epsilon$-independent quantity. Once again, our tentative conclusion that all
the singular values could be safely truncated as $\epsilon\to 0$ would lead to
an incorrect description of dynamics, as the resulting product-state form of
the IM would be incompatible with the existence of finite temporal
correlations in the continuous-time limit.

The behavior of the IM at $\epsilon\to 0$ in the solvable cases discussed
above indicates that caution must be exercised when applying conventional
compression techniques to influence functionals. Thus, in the future it would
be beneficial to develop other numerical schemes that take into account the
structural constraints of IMs in the continuous-time limit.

##  IV Non-ergodic dynamics: discrete time crystals

Figure 6: Depiction of the tensor network corresponding to the discrete path
integral of the Floquet model with fully random $J$. Averaging over $J$ can be
implemented by an operator $\delta_{\sum s\sigma\sum{\bar{s}\bar{\sigma}}}$.
This operator can be implemented as an MPO by introducing a virtual index
which at each spin corresponds to the difference in partial sums up to this
spin. The condition can be enforced by only allowing partial sum values that
are consistent with it. The maximal bond dimension of the resulting operator
is $t$. Figure 7: Discrete time crystal simulated using influence functional
approach. We use a Floquet model in Eq. (17) with parameters
$g=\frac{\pi}{2}-0.1,h=0.3$, and $J$ fully random. The initial state
corresponds to an infinite temperature ensemble. Left: $C_{zz}$ correlations
show persistent oscillations. For larger time, the amplitude of oscillations
decreases for bond dimension $\chi=32$, which is attributed to the effects of
the compression. Right: TE grows logarithmically with evolution time. Two
curves for different bond dimensions coincide.

The influence functional approach allows to treat the effects of spatially
uncorrelated randomness via exact ensemble averaging, denoted by
$\mathbb{E}\big{(}\cdot\big{)}$. As in the Schwinger-Keldysh path integral
formalism Kamenev (2011), ensemble-averaged observables can be found from the
ensemble-averaged IF, without the need of introducing replicas. Provided
randomness in different spatial points is uncorrelated (or short-range
correlated), and provided its distribution is spatially homogeneous, the
ensemble-averaged IF satisfies a self-consistency equation defined by the
transfer matrix $\mathbb{E}\big{(}\tilde{T}\big{)}$. This applies to various
kinds of randomness. Above we have already introduced ensemble-averaging over
statistical mixtures of initial states (in particular, infinite-temperature
density matrices). Random spatiotemporal noise can also be easily
incorporated; in Ref. Lerose  _et al._ (2021), we showed that circuits with
fully random Haar-distributed kicks and fully random diagonal interactions
behave as perfect dephasers (almost surely in the ensemble distribution in the
limit of large local Hilbert space). Furthermore, in Ref. Sonner  _et al._
(2020), the formalism was applied to systems with quenched disorder. In this
case, the random variables are constant in time, and exact averaging over them
introduces non-local in time effective interactions in the problem.

One of the most appealing aspects of the IF approach is that the self-
consistent IF fully encodes the ability of a quantum many-body system to
provide an efficient thermal bath. In particular, it should be possible to
classify universality classes of quantum dynamics based on the structure of
their IF, e.g., to distinguish between ergodicity vs non-ergodicity. It is
known that strongly disordered quantum lattice systems can break ergodicity
via the mechanism of many-body localization (MBL) Abanin  _et al._ (2019).
Reference Sonner  _et al._ (2020) began to explore the MBL transition using
the IF approach. The model considered in Ref. Sonner  _et al._ (2020) is the
kicked Ising chain of Eq. (8) with random longitudinal fields $h_{j}$
uniformly distributed in $[0,2\pi)$. This model exhibits a Floquet-MBL
transition upon varying the relative strength of the transverse field and/or
of interactions Zhang  _et al._ (2016a); Sonner  _et al._ (2020).

Discrete time crystals (DTCs) represent a notable regime of non-equilibrium
quantum dynamics. DTC behavior in periodically driven quantum lattice systems
is associated with a robust spatiotemporal order that spontaneously breaks the
discrete time-translation symmetry. A DTC phase requires a coherent periodic
drive and an underlying physical mechanism protecting the system from
indefinite heating, such as MBL or prethermalization Abanin  _et al._ (2017);
Else  _et al._ (2017). DTC response is robust against arbitrary (weak)
perturbations of the system interactions and of the driving protocol Khemani
_et al._ (2016); Else  _et al._ (2016). Below we show that the IM approach
provides a suitable tool for this phenomenon as well. Indeed, the IM formalism
naturally incorporates the necessary ingredients to observe time-crystalline
behavior, namely discrete time-translation symmetry, thermodynamic limit, and
disorder averaging.

DTC response can be observed, in particular, in a kicked Ising model with
random interactions and single-spin kicks close to perfect $\pi$-rotations
around the $\hat{x}$ axis Khemani  _et al._ (2016); Else  _et al._ (2016). We
consider the model in Eq. (8), introducing bond-dependent couplings $J_{j}$
uniformly distributed in $[0,2\pi]$, and setting $g=\pi/2-\varepsilon$:

$U=e^{-i(\pi/2-\varepsilon)\sum_{j}\sigma_{j}^{x}}\,e^{-\sum_{j}J_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}-h\sum_{j}\sigma_{j}^{z}}.$
(17)

For $\varepsilon=0$, the kick commutes with the interaction term: every
initial product state in the $z$-basis undergoes perfect flips of all the
spins each period, and thus returns to the initial configuration every other
driving period, exhibiting a trivial DTC dynamics. When $\varepsilon$ is non-
zero, but sufficiently small, the subharmonic DTC response of the system
survives many-body quantum fluctuations. The robustness of the DTC is made
possible by MBL, which prevents the domain-wall excitations generated by the
periodic quenches from spreading across the chain and melting the
spatiotemporally ordered pattern, as it would happen in the absence of
disorder.

To analyze DTC behavior using the IF formalism, we first note that averaging
over $J$ in the dual transfer matrix can be done exactly:

$\int\frac{\mathrm{d}J}{2\pi}\prod_{\tau}\exp\left[iJ(\sigma^{\tau}s^{\tau}-\bar{\sigma}^{\tau}\bar{s}^{\tau})\right]=\delta_{\sum_{\tau}{\sigma^{\tau}s^{\tau}},\sum_{\tau}{\bar{\sigma}^{\tau}\bar{s}^{\tau}}}\,.$
(18)

The resulting operator can be interpreted as an adjacency matrix, which allows
a transition between trajectories $\sigma_{\tau},\bar{\sigma}_{\tau}$ and
$s_{\tau},\bar{s}_{\tau}$ as long as the total number of spins flipped on the
forward trajectory is equal to the total number of spins flipped on the
backward trajectory. We can express this operator as a MPO as depicted in Fig.
6: Each tensor in the MPO is given by
$\delta_{B+\sigma^{\tau}s^{\tau}-\bar{\sigma}^{\tau}\bar{s}^{\tau},B^{\prime}}$
where the virtual indices $B,B^{\prime}$ carry the running sums of flipped
spins on the forward minus those on the backward trajectory, such that the
averaging condition (18) is globally enforced. This procedure yields an MPO
with bond dimension $t$. Since bond dimension of the operator is relatively
large, we cannot apply the MPO to the MPS and compress the result. Instead,
the compression needs to happen during the application of the operator
Stoudenmire and White (2010).

To detect subharmonic response, we computed the (disorder-averaged) IM using
an MPO-based algorithm, and used it to calculate the dynamical correlation
function $\mathcal{C}_{zz}(T)$. The model (17), similar to closely related
models considered in Refs. Khemani  _et al._ (2016); Else  _et al._ (2016),
exhibits a DTC phase at sufficiently small detuning from a perfect
$\pi$-rotation, $\varepsilon$. This is evident in Fig. 7, which illustrates
$\mathcal{C}_{zz}(T)$ for $\varepsilon=0.1,h=0.3$. The correlator exhibits
persistent oscillations with an amplitude that stays sizeable at the times
simulated. Note the slight decay at times $T>20$ for the smaller bond
dimension; we attribute it to the errors introduced by compression to an MPS
form. Indeed, increasing the bond dimension $\chi$ restores the stable
oscillatory behavior for longer times. We have further observed that for
$\varepsilon\gtrsim 0.3$ the behavior changes qualitatively, and the
correlation function decays at short, $\chi$-independent time scales,
suggesting that the system enters a thermalizing phase.

To assess the efficiency of the method, we have also analyzed the scaling of
TE entropy with $T$ (right panel of Fig. 7, converged with respect to $\chi$),
finding it to be approximately logarithmic. To understand the origin of this
behavior, let us consider the trivial DTC point $\varepsilon=0$, i.e.,
$g=\pi/2$. In this case quantum fluctuations are completely suppressed, and,
analogously to Eq. (11), the exact IM for spin $j$ reads

$I[\sigma_{\tau},\bar{\sigma}_{\tau}]=\cos\big{[}J_{j}\sum_{\tau}(-)^{\tau}(\sigma_{\tau}-\bar{\sigma}_{\tau})\big{]}.$
(19)

Taking the exact average over the random coupling $J_{j}$, we obtain

$\mathbb{E}\big{(}I[\sigma_{\tau},\bar{\sigma}_{\tau}]\big{)}=\delta_{\sum_{\tau}(-)^{\tau}\sigma_{\tau}\,,\,\sum_{\tau}(-)^{\tau}\bar{\sigma}_{\tau}}.$
(20)

Viewed as a wavefunction in the folded multi-time Hilbert space, similarly to
what done in Eq. (12), we find that

$\mathbb{E}\big{(}I[\sigma_{\tau},\bar{\sigma}_{\tau}]\big{)}\mapsto\Psi=\bigg{(}\bigotimes_{\tau}\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
1\end{pmatrix}_{cl,\tau}\bigg{)}\otimes\bigg{\lvert}\frac{T}{2},0\bigg{\rangle}_{q}$
(21)

where $\ket{\frac{T}{2},0}_{q}$ represents the Dicke state of $T$ spins-$1/2$
with staggered collective spin magnitude $S=\frac{T}{2}$ and magnetization
$S_{z}=0$. The entanglement entropy of this wavefunction scales as
$\frac{1}{2}\log T$. This form of temporal entanglement stems from the
effective long-range interactions in time arising from exact disorder
averaging. (We observe that a similar conclusion applies to the case $g=0$,
where the collective spin magnetization is uniform rather than staggered.)

We finally note that Ref. Sonner  _et al._ (2020), which considered Floquet-
MBL phase in a random-$h$ kicked Ising model, found that the scaling of TE
entropy was strongly sublinear in the MBL phase, qualitatively similar to what
is reported above. However, we note that there a variational DMRG approach to
computing the IM was used, while here we employed iterative algorithm. In the
future, it would be interesting to compare the performance of the two
approaches to solving self-consistency equation for different models.

##  V Summary and outlook

In this paper, we characterized dynamics of a many-body system via its
influence functional – an object commonly used for describing dynamics of open
quantum systems. IF provides a tool for quantifying the ability (or lack
thereof) of a many-body systems to serve as a thermal bath for its parts
Lerose  _et al._ (2021); Sonner  _et al._ (2020). One of our goals was to
illustrate the versatility of this new approach, by applying it to describe
dynamics in very different regimes and setups. In particular, we studied
thermalization in chaotic Floquet and Hamiltonian models, as well as its
absence, accompanying the DTC response, in disordered systems.

A central role in our analysis was played by the temporal entanglement of the
IF “wavefunction”, which is an analogue of the conventional real-space
entanglement of many-body wavefunctions. Similarly to the case of ground
states in many-body systems, temporal entanglement serves as a measure of the
computational complexity of parametrizing the IF by a MPS. As discussed above,
in several physical situations of interest, temporal entanglement remains low,
providing new corners where many-body dynamics can be efficiently simulated
(and, in some cases even solvable). Furthermore, we expect that other measures
of temporal entanglement of IF could provide a witness of universal aspects of
dynamics, for example, distinguishing between ergodic and non-ergodic
dynamics.

Looking forward, the IF formalism appears suitable for studying entanglement
transitions in hybrid quantum circuits composed of unitary gates and
projective measurements or local dissipation Skinner  _et al._ (2019); Li  _et
al._ (2019): indeed, both randomness and dissipation can be naturally
incorporated in this approach. We expect that the character of the asymptotic
state of such systems is encoded in the structure of the self-consistent,
ensemble-averaged influence matrix. On a more practical level, the development
of efficient numerical schemes to approximate the IF of generic quantum many-
body systems, while challenging, is a promising direction for future research,
which may provide new insights into several long-standing problems in quantum
dynamics, including the nature of the MBL-thermal transition Potter  _et al._
(2015); Vosk  _et al._ (2015); Zhang  _et al._ (2016b); Dumitrescu  _et al._
(2017); Thiery  _et al._ (2018); Goremykina  _et al._ (2019).

Finally, we hope that the approach outlined above will stimulate new fruitful
connections between theory of open quantum systems and quantum many-body
dynamics Strathearn  _et al._ (2018); Cygorek  _et al._ (2021); Ye and Kin-Lic
Chan (2021). For example, one can envision that standard concepts of the
former – such as memory kernels Makarov and Makri (1994) and non-Markovianity
measures Breuer  _et al._ (2016) – may be turned into new tools for the
latter. Progress in this direction may enable classification of universal
basins of attraction in quantum many-body dynamics. In addition, ideas and
tools from quantum many-body physics, including tensor-networks, could widen
the scope of open quantum systems theory, for example, allowing one to
describe the effects of truly interacting, quantum chaotic environments.

##  VI Acknowledgments

This work was supported by the Swiss National Science Foundation and by the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No. 864597). We thank
Lorenzo Piroli and Soonwon Choi for useful discussions. Computations were
performed at the University of Geneva on the “Baobab” and “Yggdrasil” HPC
clusters.

## References

  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, 2007).
  * Weiss (2012) Ulrich Weiss,  _Quantum Dissipative Systems_, 4th ed. (WORLD SCIENTIFIC, 2012) https://www.worldscientific.com/doi/pdf/10.1142/8334 .
  * Feynman and Vernon (1963) R.P Feynman and F.L Vernon, “The theory of a general quantum system interacting with a linear dissipative system,” Annals of Physics 24, 118 – 173 (1963).
  * Leggett  _et al._ (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, Matthew P. A. Fisher, Anupam Garg, and W. Zwerger, “Dynamics of the dissipative two-state system,” Rev. Mod. Phys. 59, 1–85 (1987).
  * Bloch  _et al._ (2008) Immanuel Bloch, Jean Dalibard, and Wilhelm Zwerger, “Many-body physics with ultracold gases,” Rev. Mod. Phys. 80, 885–964 (2008).
  * Blatt and Roos (2012) R. Blatt and C. F. Roos, “Quantum simulations with trapped ions,” Nature Physics 8, 277–284 (2012).
  * D’Alessio  _et al._ (2016) Luca D’Alessio, Yariv Kafri, Anatoli Polkovnikov, and Marcos Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Advances in Physics 65, 239–362 (2016).
  * 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).
  * 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).
  * Alet and Laflorencie (2018) Fabien Alet and Nicolas Laflorencie, “Many-body localization: An introduction and selected topics,” Comptes Rendus Physique 19, 498 – 525 (2018).
  * Anderson (1958) P. W. Anderson, “Absence of diffusion in certain random lattices,” Phys. Rev. 109, 1492–1505 (1958).
  * Ponte  _et al._ (2015) Pedro Ponte, Z. Papić, François Huveneers, and Dmitry A. Abanin, “Many-body localization in periodically driven systems,” Phys. Rev. Lett. 114, 140401 (2015).
  * Lazarides  _et al._ (2015) Achilleas Lazarides, Arnab Das, and Roderich Moessner, “Fate of many-body localization under periodic driving,” Phys. Rev. Lett. 115, 030402 (2015).
  * Abanin  _et al._ (2016) Dmitry A. Abanin, Wojciech De Roeck, and François Huveneers, “Theory of many-body localization in periodically driven systems,” Annals of Physics 372, 1 – 11 (2016).
  * Bordia  _et al._ (2017) Pranjal Bordia, Henrik Lüschen, Ulrich Schneider, Michael Knap, and Immanuel Bloch, “Periodically driving a many-body localized quantum system,” Nature Physics 13, 460–464 (2017).
  * Khemani  _et al._ (2016) Vedika Khemani, Achilleas Lazarides, Roderich Moessner, and S. L. Sondhi, “Phase structure of driven quantum systems,” Phys. Rev. Lett. 116, 250401 (2016).
  * Else  _et al._ (2016) Dominic V. Else, Bela Bauer, and Chetan Nayak, “Floquet time crystals,” Phys. Rev. Lett. 117, 090402 (2016).
  * Nathan  _et al._ (2019) Frederik Nathan, Dmitry Abanin, Erez Berg, Netanel H. Lindner, and Mark S. Rudner, “Anomalous floquet insulators,” Phys. Rev. B 99, 195133 (2019).
  * Essler and Fagotti (2016) Fabian H L Essler and Maurizio Fagotti, “Quench dynamics and relaxation in isolated integrable quantum spin chains,” Journal of Statistical Mechanics: Theory and Experiment 2016, 064002 (2016).
  * Serbyn  _et al._ (2020) Maksym Serbyn, Dmitry A. Abanin, and Zlatko Papić, “Quantum Many-Body Scars and Weak Breaking of Ergodicity,” arXiv e-prints , arXiv:2011.09486 (2020), arXiv:2011.09486 [quant-ph] .
  * Deutsch (1991) J. M. Deutsch, “Quantum statistical mechanics in a closed system,” Phys. Rev. A 43, 2046–2049 (1991).
  * Srednicki (1994) Mark Srednicki, “Chaos and quantum thermalization,” Phys. Rev. E 50, 888–901 (1994).
  * Verstraete  _et al._ (2008) F. Verstraete, V. Murg, and J.I. Cirac, “Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems,” Advances in Physics 57, 143–224 (2008).
  * Paeckel  _et al._ (2019) Sebastian Paeckel, Thomas Köhler, Andreas Swoboda, Salvatore R. Manmana, Ulrich Schollwöck, and Claudius Hubig, “Time-evolution methods for matrix-product states,” Annals of Physics 411, 167998 (2019).
  * Lerose  _et al._ (2021) Alessio Lerose, Michael Sonner, and Dmitry A. Abanin, “Influence matrix approach to many-body floquet dynamics,” (2021), arXiv:2009.10105 [cond-mat.str-el] .
  * Ye and Kin-Lic Chan (2021) Erika Ye and Garnet Kin-Lic Chan, “Constructing Tensor Network Influence Functionals for General Quantum Dynamics,” arXiv e-prints , arXiv:2101.05466 (2021), arXiv:2101.05466 [quant-ph] .
  * Bañuls  _et al._ (2009) M. C. Bañuls, M. B. Hastings, F. Verstraete, and J. I. Cirac, “Matrix product states for dynamical simulation of infinite chains,” Phys. Rev. Lett. 102, 240603 (2009).
  * Müller-Hermes  _et al._ (2012) Alexander Müller-Hermes, J Ignacio Cirac, and Mari Carmen Banuls, “Tensor network techniques for the computation of dynamical observables in one-dimensional quantum spin systems,” New Journal of Physics 14, 075003 (2012).
  * Akila  _et al._ (2016) M Akila, D Waltner, B Gutkin, and T Guhr, “Particle-time duality in the kicked ising spin chain,” Journal of Physics A: Mathematical and Theoretical 49, 375101 (2016).
  * Bertini  _et al._ (2019) Bruno Bertini, Pavel Kos, and Toma ž Prosen, “Exact correlation functions for dual-unitary lattice models in $1+1$ dimensions,” Phys. Rev. Lett. 123, 210601 (2019).
  * Piroli  _et al._ (2020) Lorenzo Piroli, Bruno Bertini, J. Ignacio Cirac, and Toma ž Prosen, “Exact dynamics in dual-unitary quantum circuits,” Phys. Rev. B 101, 094304 (2020).
  * Klobas  _et al._ (2020) Katja Klobas, Bruno Bertini, and Lorenzo Piroli, “Exact thermalization dynamics in the “Rule 54” Quantum Cellular Automaton,” arXiv e-prints , arXiv:2012.12256 (2020), arXiv:2012.12256 [cond-mat.stat-mech] .
  * Sonner  _et al._ (2020) Michael Sonner, Alessio Lerose, and Dmitry A. Abanin, “Characterizing many-body localization via exact disorder-averaged quantum noise,” arXiv e-prints , arXiv:2012.00777 (2020), arXiv:2012.00777 [cond-mat.dis-nn] .
  * Strathearn  _et al._ (2018) A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, “Efficient non-Markovian quantum dynamics using time-evolving matrix product operators,” Nature Communications 9, 3322 (2018).
  * Cygorek  _et al._ (2021) Moritz Cygorek, Michael Cosacchi, Alexei Vagov, Vollrath Martin Axt, Brendon W. Lovett, Jonathan Keeling, and Erik M. Gauger, “Numerically-exact simulations of arbitrary open quantum systems using automated compression of environments,” arXiv e-prints , arXiv:2101.01653 (2021), arXiv:2101.01653 [quant-ph] .
  * Kormos  _et al._ (2017) M. Kormos, M. Collura, G. Takács, and P. Calabrese, “Real time confinement following a quantum quench to a non-integrable model,” Nature Physics 13, 246–249 (2017).
  * James  _et al._ (2019) Andrew J. A. James, Robert M. Konik, and Neil J. Robinson, “Nonthermal states arising from confinement in one and two dimensions,” Phys. Rev. Lett. 122, 130603 (2019).
  * Mazza  _et al._ (2019) Paolo Pietro Mazza, Gabriele Perfetto, Alessio Lerose, Mario Collura, and Andrea Gambassi, “Suppression of transport in nondisordered quantum spin chains due to confined excitations,” Phys. Rev. B 99, 180302 (2019).
  * Lerose  _et al._ (2020) Alessio Lerose, Federica M. Surace, Paolo P. Mazza, Gabriele Perfetto, Mario Collura, and Andrea Gambassi, “Quasilocalized dynamics from confinement of quantum excitations,” Phys. Rev. B 102, 041118 (2020).
  * Lerose  _et al._ (2021) Alessio Lerose, Michael Sonner, and Dmitry A. Abanin, “to appear,” arXiv e-prints (2021).
  * Sünderhauf  _et al._ (2018) Christoph Sünderhauf, David Pérez-García, David A. Huse, Norbert Schuch, and J. Ignacio Cirac, “Localization with random time-periodic quantum circuits,” Phys. Rev. B 98, 134204 (2018).
  * Chan  _et al._ (2020) Amos Chan, Andrea De Luca, and J. T. Chalker, “Spectral lyapunov exponents in chaotic and localized many-body quantum systems,” (2020), arXiv:2012.05295 [cond-mat.stat-mech] .
  * Garratt and Chalker (2020) S. J. Garratt and J. T. Chalker, “Many-body delocalisation as symmetry breaking,” (2020), arXiv:2012.11580 [cond-mat.stat-mech] .
  * Haegeman and Verstraete (2017) Jutho Haegeman and Frank Verstraete, “Diagonalizing transfer matrices and matrix product operators: A medley of exact and computational methods,” Annual Review of Condensed Matter Physics 8, 355–406 (2017), https://doi.org/10.1146/annurev-conmatphys-031016-025507 .
  * Hastings and Mahajan (2015) M. B. Hastings and R. Mahajan, “Connecting entanglement in time and space: Improving the folding algorithm,” Phys. Rev. A 91, 032306 (2015).
  * Kim  _et al._ (2014) Hyungwon Kim, Tatsuhiko N. Ikeda, and David A. Huse, “Testing whether all eigenstates obey the eigenstate thermalization hypothesis,” Phys. Rev. E 90, 052105 (2014).
  * Hauschild and Pollmann (2018) Johannes Hauschild and Frank Pollmann, “Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy),” SciPost Phys. Lect. Notes 5 (2018), 10.21468/SciPostPhysLectNotes.5, code available from https://github.com/tenpy/tenpy, arXiv:1805.00055 .
  * Sachdev (2011) S. Sachdev,  _Quantum Phase Transitions_ (Cambridge University Press, 2011).
  * McCoy and Wu (1978) Barry M. McCoy and Tai Tsun Wu, “Two-dimensional ising field theory in a magnetic field: Breakup of the cut in the two-point function,” Phys. Rev. D 18, 1259–1267 (1978).
  * Bañuls  _et al._ (2011) M. C. Bañuls, J. I. Cirac, and M. B. Hastings, “Strong and weak thermalization of infinite nonintegrable quantum systems,” Phys. Rev. Lett. 106, 050405 (2011).
  * Lin and Motrunich (2017) Cheng-Ju Lin and Olexei I. Motrunich, “Quasiparticle explanation of the weak-thermalization regime under quench in a nonintegrable quantum spin chain,” Phys. Rev. A 95, 023621 (2017).
  * Wurtz and Polkovnikov (2020) Jonathan Wurtz and Anatoli Polkovnikov, “Emergent conservation laws and nonthermal states in the mixed-field ising model,” Phys. Rev. B 101, 195138 (2020).
  * Kamenev (2011) A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, 2011).
  * Zhang  _et al._ (2016a) Liangsheng Zhang, Vedika Khemani, and David A. Huse, “A floquet model for the many-body localization transition,” Phys. Rev. B 94, 224202 (2016a).
  * Sonner  _et al._ (2020) Michael Sonner, Maksym Serbyn, Zlatko Papic, and Dmitry A. Abanin, “Thouless energy in a floquet model of many-body localization,” (2020), to appear .
  * Abanin  _et al._ (2017) Dmitry A. Abanin, Wojciech De Roeck, Wen Wei Ho, and Fran çois Huveneers, “Effective hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems,” Phys. Rev. B 95, 014112 (2017).
  * Else  _et al._ (2017) Dominic V. Else, Bela Bauer, and Chetan Nayak, “Prethermal phases of matter protected by time-translation symmetry,” Phys. Rev. X 7, 011026 (2017).
  * Stoudenmire and White (2010) E M Stoudenmire and Steven R White, “Minimally entangled typical thermal state algorithms,” New Journal of Physics 12, 055026 (2010).
  * Skinner  _et al._ (2019) Brian Skinner, Jonathan Ruhman, and Adam Nahum, “Measurement-induced phase transitions in the dynamics of entanglement,” Phys. Rev. X 9, 031009 (2019).
  * Li  _et al._ (2019) Yaodong Li, Xiao Chen, and Matthew P. A. Fisher, “Measurement-driven entanglement transition in hybrid quantum circuits,” Phys. Rev. B 100, 134306 (2019).
  * Potter  _et al._ (2015) Andrew C. Potter, Romain Vasseur, and S. A. Parameswaran, “Universal properties of many-body delocalization transitions,” Phys. Rev. X 5, 031033 (2015).
  * Vosk  _et al._ (2015) Ronen Vosk, David A. Huse, and Ehud Altman, “Theory of the many-body localization transition in one-dimensional systems,” Phys. Rev. X 5, 031032 (2015).
  * Zhang  _et al._ (2016b) Liangsheng Zhang, Bo Zhao, Trithep Devakul, and David A. Huse, “Many-body localization phase transition: A simplified strong-randomness approximate renormalization group,” Phys. Rev. B 93, 224201 (2016b).
  * Dumitrescu  _et al._ (2017) Philipp T. Dumitrescu, Romain Vasseur, and Andrew C. Potter, “Scaling theory of entanglement at the many-body localization transition,” Phys. Rev. Lett. 119, 110604 (2017).
  * Thiery  _et al._ (2018) Thimothée Thiery, Fran çois Huveneers, Markus Müller, and Wojciech De Roeck, “Many-body delocalization as a quantum avalanche,” Phys. Rev. Lett. 121, 140601 (2018).
  * Goremykina  _et al._ (2019) Anna Goremykina, Romain Vasseur, and Maksym Serbyn, “Analytically solvable renormalization group for the many-body localization transition,” Phys. Rev. Lett. 122, 040601 (2019).
  * Makarov and Makri (1994) Dmitrii E. Makarov and Nancy Makri, “Path integrals for dissipative systems by tensor multiplication. Condensed phase quantum dynamics for arbitrarily long time,” Chemical Physics Letters 221, 482–491 (1994).
  * Breuer  _et al._ (2016) Heinz-Peter Breuer, Elsi-Mari Laine, Jyrki Piilo, and Bassano Vacchini, “Colloquium: Non-markovian dynamics in open quantum systems,” Rev. Mod. Phys. 88, 021002 (2016).


supplemental_material

# Energies and spectra of solids from the algorithmic inversion of localized
$GW$

Tommaso Chiarotti<EMAIL_ADDRESS>Theory and Simulations of
Materials (THEOS) and National Centre for Computational Design and Discovery
of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, 1015
Lausanne, Switzerland  Andrea Ferretti  Centro S3, CNR–Istituto Nanoscienze,
41125 Modena, Italy  Nicola Marzari  Theory and Simulations of Materials
(THEOS) and National Centre for Computational Design and Discovery of Novel
Materials (MARVEL), École Polytechnique Fédérale de Lausanne, 1015 Lausanne,
Switzerland

###### Abstract

Energy functionals of the Green’s function are designed to provide spectral
and thermodynamic properties of materials. In this work, following a
localized-GW perspective, we design an approximation to the $\Phi$ term in the
Klein functional that generalizes the energy functional of DFT+U to host a
dynamical screened-potential $U(\omega)$. The resulting effective Hamiltonians
feature dynamical potentials, and, in order to solve the problem, we extend to
crystals the algorithmic inversion method – a numerical approach aimed at
solving the resulting Dyson equation. To demonstrate the method, we compute
the spectrum and energetics of SrVO3, finding results in excellent agreement
with experiments and state-of-the-art computational methods at a negligible
computational cost.

###### pacs:

Introduction. The computation of accurate and predictive materials properties
has been and remains a formidable task marzari_electronic-structure_2021.
Though density-functional theory (DFT) has provided a huge step forward in the
prediction of a number of ground-state properties
hohenberg_inhomogeneous_1964; kohn_self-consistent_1965, addressing
spectroscopic quantities remains challenging borghi_koopmans-compliant_2014;
ferretti_bridging_2014; perdew_understanding_2017; nguyen_koopmans-
compliant_2018; de_gennaro_blochs_2022; colonna_koopmans_2022. Green’s
function methods have been used to overcome these limitations, including many-
body perturbation theory (MPBT) approaches such as $GW$ hedin_effects_1970;
reining_gw_2018, dynamical mean-field theory georges_dynamical_1996;
kotliar_electronic_2006 (DMFT), spectral functional theories
savrasov_spectral_2004; vanzini_connector_2022 (SFT), or electron-boson
interaction schemes Caruso2018PRB; Zhou2020PNAS. In all of these methods,
dynamical, frequency-dependent self-energies or potentials appear.

The Luttinger-Ward (LW) and Klein functionals luttinger_ground-state_1960;
baym_conservation_1961; baym_self-consistent_1962; martin_interacting_2016;
stefanucci_nonequilibrium_2013 are energy functionals of the Green’s function
that are variational and result in conserving dynamical potentials at self-
consistency. Though explicitly known diagrammatically, the exact $\Phi$ term
of the functional is computationally inaccessible and needs to be approximated
almbladh_variational_1999; dahlen_variational_2004. The choice of the
approximation for $\Phi$ determines the physics accessible to the functional,
ranging from describing long-range plasmonic effects, as in $GW$
hedin_effects_1970; reining_gw_2018, to strongly-correlated local interactions
as in DMFT georges_dynamical_1996; kotliar_electronic_2006. The stationarity
condition of the functional yields the Dyson equation, involving the dynamical
self-energy (as derivative of the $\Phi$ term wrt $G$) and the interacting
propagator martin_interacting_2016; stefanucci_nonequilibrium_2013. The
functional and its derivative determine the thermodynamics and the spectrum of
the material, also allowing, at self-consistency, for the use of the Hellman-
Feynman theorem to compute ground-state quantities such as forces and phonons
kotliar_electronic_2006.

Mainly due to the presence of dynamical quantities, the applications of a
scheme similar to DFT, for which both spectroscopic and thermodynamics are
computed at high numerical accuracy, are limited chiarotti_unified_2022. In
the context of $GW$, such calculations were done for model systems, such as
the homogeneous electron gas holm_fully_1998, or Hubbard chains
honet_exact_2022, (AF: I’ll fix refs) and later extended to solids using an
imaginary axis formalism kutepov_ground-state_2009; kutepov_electronic_2012;
Grumet2018PRB; yeh_fully_2022. The latter approach has been proven very
effective for the prediction of ground-state properties, but with a limited
accuracy on the spectrum due to the analytic continuation to the real axis
kutepov_one-electron_2017. In the DMFT framework, an efficient implementation
of a LW functional has been worked out, allowing for the calculation of some
ground-state quantities for crystals, together with the possibility of
accessing the spectrum haule_free_2015; kocer_efficient_2020. This comes
though at the price of numerical limitations related to the cost of the
solvers kotliar_electronic_2006, and of extra conceptual difficulties such as
the definition of the impurity model or the extension of standard electronic
structure constructions (eg Born-effective charges kocer_efficient_2020 or
topological invariants XXX). (AF: forse ci possiamo fermare all’impurity
model…)

In the Letter, within the Klein functional formalism, we propose an exchange-
correlation $\Phi$ term that generalizes the DFT+U energy functional of
Dudarev et al. dudarev_electron-energy-loss_1998 to host a dynamical screened
potential $U(\omega)$ instead of a static $U$. In doing so, the derivative of
the functional yields the localized-$GW$ self-energy of Ref. sakuma_ab_2014,
together with a proper choice for the double-counting term. In order to
address the resulting electronic structure problem – involving dynamical
potentials – we analyze the Dyson equation as a non-linear eigenvalue problem,
and solve it via an exact mapping to a non-interacting system. We do so by
generalizing the so called “algorithmic inversion method on sum-over-poles”
(AIM-SOP) chiarotti_unified_2022 to the case of realistic systems. In
practice, the single-particle Green’s function, along with the Dyson orbitals
and excitation energies of the system, is found and used to calculate accurate
spectral and thermodynamic properties of the material. The procedure avoids,
e.g., problems with the quasi-particle solution of the Dyson equation, or
numerical difficulties when integrating the Green’s function to obtain
thermodynamic quantities chiarotti_unified_2022. We apply the framework to
SrVO3 finding results in excellent agreement with experiments or state-of-the-
art methods for the one-particle spectra, bulk-modulus, and $\Gamma$ phonons,
at a negligible computational cost.

Algorithmic inversion method. First, we introduce a general framework to solve
electronic structure problems involving dynamical potentials (here represented
by $\Sigma(\omega)$). This task is equivalent to find the Green’s function
solving the Dyson equation,
$\hat{G}(\omega)=\left[\omega-\hat{h}_{0}-\hat{\Sigma}(\omega)\right]^{-1},$
which can be mapped math-paper to the non-linear eigenvalue problem

yields a non-linear eigenvalue problem when trying to find the sum-over-poles
form chiarotti_unified_2022 of the Green’s function
$\hat{G}(\omega)=\sum_{l}\frac{\hat{A}_{l}}{\omega-z_{l}}$. The solution gives
access to the Dyson orbitals $\ket{\psi_{l}}$, excitation energies
$\epsilon_{l}$, and particle lifetimes $\delta_{l}$—related to the amplitude,
energy, and decay of the quasi-particles—of the system. Indeed, by solving the
right (left) non-linear eigenvalue equation
$\left[\hat{h}_{0}-\hat{\Sigma}(\omega)\right]\ket{\psi_{l}}=\omega\ket{\psi_{l}}$,
it is possible to see that

$\sum_{l}\frac{\ket{\psi_{l}^{R}}\bra{\psi_{l}^{L}}}{\omega-
z_{l}}=\left[\omega-\hat{h}_{0}-\hat{\Sigma}(\omega)\right]^{-1},$ (1)

where $\ket{\psi_{l}^{R}}$ ($\bra{\psi_{l}^{L}}$) is the right (left)
eigenvector with eigenvalue $z_{l}=\epsilon_{l}+i\delta_{l}$ of the non-linear
eigenvalue equation. This can be seen, e.g., by Taylor expanding around each
$z_{i}$ the operator $\omega-\hat{h}_{0}-\hat{\Sigma}(\omega)=\left(\omega-
z_{i}\right)\ket{\psi_{i}^{R}}\bra{\psi_{i}^{L}}+\sum_{j}f^{i}_{j}(\omega)\ket{\psi_{j}^{R}}\bra{\psi_{j}^{L}}$,
with $f^{i}_{i}(\omega)\propto o\left(\omega-z_{i}\right)$, thus showing that
the operator in the rhs of Eq. (1) has one-order poles in $z_{i}$.
Furthermore, taking into account the decay of the self-energy farid_gw_1988
and that $\hat{1}=\sum_{l}\ket{\psi_{l}^{R}}\bra{\psi_{l}^{L}}$, the lhs and
rhs of Eq. (1) have both $\omega\to\infty$ limit equal to $\hat{1}$. Thus, the
equality in (1) is proven (on any basis) since two meromorphic functions with
same poles and $\omega\to\infty$ limit must be the same thomas_calculus_1988.
(TC: Paragraph to be rewritten also in view of the Supplemental Material)

Within Green’s function frameworks, a natural way of expressing time-ordered
propagators, e.g., the dynamical part of the self-energy—the static
contribution is thought to be included in $\hat{h}_{0}$—, is provided by the
truncated Lehmann representation engel_calculation_1991 or the so called “sum-
over-poles” form chiarotti_unified_2022:

$\hat{\Sigma}(\omega)=\sum_{l=1}^{N}\frac{\hat{\Gamma}_{l}}{\omega-\Omega_{l}},$
(2)

where the amplitudes $\hat{\Gamma}_{l}$ are operators acting on the one-
particle Hilbert space defined by $\hat{h}_{0}$, and
$\Omega_{l}=\omega_{l}+i\delta_{l}$ are scalars respecting the time-ordering
of the operator, $\delta_{l}\gtrless 0$ when $\Omega_{l}\gtrless\mu$ with
$\mu$ the chemical potential of the system. Introducing a vectorial
representation, the self-energy can be written as:

$\begin{aligned}
\hat{\Sigma}(\omega)=\mathbf{\sqrt{\hat{\Gamma}^{\dagger}}}\left(\omega-\mathbf{\Omega}\right)^{-1}\mathbf{\sqrt{\hat{\Gamma}}},\end{aligned},$
(3)

where $\mathbf{\Omega}$ is a diagonal matrix with
$\Omega_{1}\hat{1},\dots,\Omega_{N}\hat{1}$ on the diagonal, and
$\mathbf{\sqrt{\hat{\Gamma}}}=\left(\sqrt{\hat{\Gamma}_{1}},\dots,\sqrt{\hat{\Gamma}_{N}}\right)$.
The Green’s function

$\hat{G}(\omega)=\left(\omega-\hat{h}_{0}-\mathbf{\sqrt{\hat{\Gamma}^{\dagger}}}\left(\omega-\mathbf{\Omega}\right)^{-1}\mathbf{\sqrt{\hat{\Gamma}}}\right)^{-1}$
(4)

corresponds to embed two non interacting systems with hamiltonians
$\hat{h}_{0}$ and $\mathbf{\Omega}$, and a right coupling between the two
$\mathbf{\sqrt{\hat{\Gamma}}}$. By going backwards in the embedding procedure,
the total Green’s function is
$\mathbf{\hat{G}}_{\text{tot}}=(\omega-\mathbf{\hat{H}}_{\text{AIM}})^{-1}$
with

$\mathbf{\hat{H}}_{\text{AIM}}=\begin{pmatrix}\hat{h}_{0}&\sqrt{\hat{\Gamma}_{1}}&\dots&\sqrt{\hat{\Gamma}_{N}}\\\
\sqrt{\hat{\Gamma}_{1}^{\dagger}}&\Omega_{1}\hat{1}&0&0\\\
\vdots&0&\ddots&0\\\
\sqrt{\hat{\Gamma}_{N}^{\dagger}}&0&\dots&\Omega_{N}\hat{1}\end{pmatrix},$ (5)

and in particular
$\hat{G}=\hat{P}_{0}(\omega-\mathbf{\hat{H}}_{\text{AIM}})^{-1}\hat{P}_{0}$,
with $\hat{P}_{0}$ projecting onto the $\hat{h}_{0}$ Hilbert space.
Consequently, the diagonalization of $\mathbf{\hat{H}}_{\text{AIM}}$ yields
the sum over poles for the Green’s function:

$\hat{G}(\omega)=\sum_{i=l}^{N+1}\frac{\ket{\psi^{R}_{l}}\bra{\psi^{L}_{l}}}{\omega-
z_{l}},$ (6)

where $\ket{\psi^{R}_{l}}$ ($\bra{\psi^{L}_{l}}$) are the projection of the
right (left) eigenvectors of $\mathbf{\hat{H}}_{\text{AIM}}$ (by
$\hat{P}_{0}$), and $z_{l}$ the eigenvalues. Hence, starting from poles and
residues of the self-energy this procedure—so called algorithmic inversion
method—grants the knowledge of the excitation energies, lifetimes, and
amplitudes (Dyson orbitals) of the system, i.e., the real and imaginary part
of the poles and their residues. Furthermore, the knowledge of the Green’s
function as sum over poles allows for the analytic calculations of its moments
and thus of the integrated thermodynamic quantities chiarotti_unified_2022.
The approach is the generalization to the non-homogeneous case of the one in
Ref. chiarotti_unified_2022, in which the homogeneous electron gas was
studied. A similar construction was used to efficiently invert the Dyson
equation in savrasov_many-body_2006, and used, e.g., in
kotliar_electronic_2006; wang_pole_2012; budich_fluctuation-induced_2012.
Other similar approaches are bintrim_full-frequency_2021; cho_simplified_2022.

Even though the algorithmic inversion method on sum over poles gives an
effective way of solving the Dyson equation mapping it into a linear
diagonalization problem, it requires the knowledge of the frequency-dependent
potential as a sum over poles. For this purpose, we generalize the DFT+U
energy functional of Dudarev et al. dudarev_electron-energy-loss_1998 to host
a frequency-dependent screened potential $U(\omega)$. We write a
localized-$GW$ Klein energy functional
$E_{\text{LocGW}}[\hat{G}]=E_{xc}[\rho]+E_{H}[\rho]+\Phi_{\text{LocGW}}[\mathbf{G}]-\Tr\left[\int
d\omega\ \frac{e^{i\omega 0^{+}}}{2\pi
i}\hat{G}_{0}^{-1}(\omega)\hat{G}(\omega)\right]-\Tr\left[\int d\omega\
\frac{e^{i\omega 0^{+}}}{2\pi i}\log{\left(\hat{G}(\omega)\right)}\right]$,
where $\hat{G}_{0}^{-1}=\omega-\hat{h}_{0}$, $\rho$ is the density of
$\hat{G}$, and $\mathbf{G}$ is a projection of the Green’s function of the
system $\hat{G}$ onto a localized “Hubbard” manifold. The generalized Hubbard
energy functional reads:

$\displaystyle\Phi_{\text{LocGW}}[\mathbf{G}]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int\frac{d\omega}{2\pi
i}\frac{d\omega^{\prime}}{2\pi i}\ e^{i\omega
0^{+}}e^{i\omega^{\prime}0^{+}}U(\omega^{\prime})\times$ $\displaystyle\times$
$\displaystyle\Tr\\{\mathbf{G}(\omega+\omega^{\prime})\left[\delta(\omega-c)\mathbf{1}-\mathbf{G}(\omega)\right]\Big{\\}}$

with $U(\omega)$ the screened potential $\hat{W}(\omega)$ of the system
projected onto the localized space and averaged over, and $c$ a constant that
helps to adjust the double counting. The Hubbard generalized correction
$\Phi_{\text{LocGW}}[\mathbf{G}]$ yields the self-energy

$\begin{aligned} \mathbf{\Sigma_{\text{LocGW}}}(\omega)=2\pi
i\frac{\delta\Phi_{\text{LocGW}}[\mathbf{G}]}{\delta\mathbf{G}}=\\\
\frac{1}{2}U(c)\mathbf{1}-\int d\omega^{\prime}\
\frac{e^{i\omega^{\prime}0^{+}}}{2\pi
i}U(\omega^{\prime})\mathbf{G}(\omega+\omega^{\prime})\end{aligned},$ (8)

which is composed by a double-counting term (first addend) and the
localized-$GW$ self-energy sakuma_ab_2014 (second addend). Here, we fix the
double counting parameter constant $c\to\infty$, $U(c)=U_{\infty}$, in
agreement with Ref. vanzini_towards_2021. Importantly, we note that
$\Phi_{\text{LocGW}}$ reduces to the Hubbard energy term of Ref.
dudarev_electron-energy-loss_1998 when a constant screened local interaction
equal to $U$ is considered. In addition, $\Phi_{\text{LocGW}}$ reduces to the
$GW_{0}$ functional (plus a double-counting term) in the limit of the Hubbard
manifold to be the entire crystal. Here, we suppose to have only one Hubbard
manifold on site $I$, the generalization to include multiples sites amounts to
sum over the sites $\Phi^{I}_{\text{LocGW}}$, each having their own projection
of $\mathbf{G}_{I}$ and a local-screened-average interaction $U_{I}(\omega)$.

Being it a Klein functional, with the screened interaction fixed to
$U(\omega)$, $E_{\text{LocGW}}[\hat{G}]$ can be made stationary by solving
self-consistently the Dyson equation (LABEL:DysonEq) with a static-term
$\hat{h}_{0}=\hat{h}_{KS}(\rho)$ and a self-energy
$\hat{\Sigma}_{\text{LocGW}}(\omega)=\sum_{m,m^{\prime}}\ket{\phi_{m}}\matrixelement{\phi_{m}}{\mathbf{\Sigma}_{\text{LocGW}}(\omega)}{\phi_{m^{\prime}}}\bra{\phi_{m^{\prime}}}$,
where $\phi_{m}$ span the localized Hubbard space. Here, the use of the
algorithmic inversion method is crucial and provides a “closed” formalism with
respect to the SOP form for the propagators. Starting with an initial Green’s
function on SOP, e.g., $\hat{G}=\hat{G}_{KS}$, and a SOP for the screened
interaction $U(\omega)$, the self-energy $\mathbf{\Sigma_{\text{LocGW}}}$ of
Eq. (8) and consequently $\hat{\Sigma}(\omega)$ can be written in
SOP—convolutions between two SOPs is a SOP chiarotti_unified_2022 and
projections are just linear operations on the residues. With a self-energy on
SOP, the algorithmic inversion can be used to find the SOP for the Green’s
function. The cycle can be repeated until self-consistency. Furthermore, the
SOP form for the Green’s function naturally allows for the accurate evaluation
of the generalized Hubbard energy of Eq. (Energies and spectra of solids from
the algorithmic inversion of localized $GW$), the chemical potential, and the
thermodynamics in general chiarotti_unified_2022.

As a case-study for the approach, we apply the formalism to the paradigmatic
example of SrVO3. Localized-$GW$ has already been proven to give qualitative
good spectral results sakuma_ab_2014. Here, we provide not only the spectrum
of the material but also integrated quantities, such as the total energy and
its derivatives (lattice parameters, bulk modulus, phonons), making use of the
functional $E_{\text{LocGW}}[\hat{G}]$. Similarly to choosing the screening
parameter $U$ in DFT+U, the choice of $U(\omega)$, accounting for the
dynamical screening of the electrons in the Hubbard localized manifold, is
needed. In accordance with the localized-$GW$ nature of the functional
miyake_effects_2013 and with vanzini_towards_2021, we consider $U(\omega)$ to
be the average over the Hubbard space of the screened interaction $W$ in the
random-phase approximation. Further details are given in the Supplemental
Material supp-info.

Though designed for self-consistency, here we limit the approach to a one-shot
calculation—i.e, a single-step in the stationarization of
$E_{\text{LocGW}}$—using as starting propagator the self-consistent Kohn-Sham
Green’s function computed using standard semi-local (PBEsol) DFT
perdew_generalized_1996. We need to evaluate the energy functional at
$E_{\text{LocGW}}[\hat{G}=\hat{G}_{\text{KS}}]$, and the same for the self-
energy, $\hat{\Sigma}_{\text{LocGW}}[\hat{G}=\hat{G}_{\text{KS}}]$. Exploiting
the algorithmic inversion we find the non-self-consistent Green’s function. At
this level, the localized-$GW$ Klein energy functional simplifies to
$E_{\text{LocGW}}[\hat{G}_{\text{KS}}]=E_{\text{DFT}}[\rho]+\Phi_{\text{LocGW}}[\mathbf{G}_{\text{KS}}]$,
simply correcting the DFT energy with a generalized Hubbard energy term. The
numerical details of the calculations are given in the Supplemental Material
supp-info.

Figure 1: Spectral function of SrVO3 (color plot) from this work compared to PBE (white-solid line) and PBE+U (white-dashed line). Only t2g bands are displayed. The chemical potential is shifted to $0$ in all the three cases. The black line passing from $0$ shows the analytic divergence of the spectral function in a metal at the chemical potential. The color map is logarithmic. | Occ. bw | Z | LS | US
---|---|---|---|---
exp. yoshida_direct_2005 | 0.7 | 0.55 | -1.5 |
exp. takizawa_coherent_2009 | 0.44 | 0.5 | -1.5 |
PBE | 1.0 | 1 |  |
PBE + U | 0.92 | 0.92 |  |
GW+DMFT tomczak_asymmetry_2014 | 0.5 | 0.5 | -1.6 | 2
GW+DMFT sakuma_ab_2014 | 0.6 | 0.6 | -1.5 | 2.5
GW+EDMFT boehnke_when_2016 |  |  | -1.7 | 2.8
This work | 0.6 | 0.6 | -1.7/-3.0 | 2.2/3.9
Table 1: Occupied bandwidth (bw), renormalization factor ($Z$), energy of
lower (LS) and upper satellites (US), and total energy $E_{\text{tot}}$ for
SrVO3 from experiments and different theoretical frameworks. Our results are
presented in the last line. The full spectral function of this work is
displayed in Fig 1. All energies (bw, LS, US) are in eV.

In Fig. 1 we show the spectral function resulting from the calculation.
Contrary to what happens in PBEsol DFT+U (PBE+U) that only shifts the chemical
potential, the spectral function coming from $\Sigma_{\text{LocGW}}$ shows a
reduction of the bandwidth along with a renormalization of the full width of
the t2g bands. In Table 1 we report the band width (bw), the renormalization
factor ($Z$), and the positions of the lower/upper satellites (LS/US). While
PBEsol DFT (PBE) and PBE+U give an overestimation of the bandwidth and $Z$,
the results from the generalized Hubbard approach are basically reproducing
the localized-$GW$ data of Ref. miyake_effects_2013 and are in general good
agreement with experiments tomczak_asymmetry_2014; yoshida_direct_2005;
takizawa_coherent_2009 and DMFT sakuma_ab_2014; boehnke_when_2016.

Figure 2: Equation of state for SrVO3 calculated using PBE (blue), PBE+U
(orange), and the generalized Hubbard approach (green). Data (crosses) are
fitted using the Birch-Murnaghan curve (in solid). The values and errors
obtained using the fit are displayed in the legend. For reference, in dashed
black we mark the experimental volume, and report the bulk modulus in the
legend. Experimental values are taken from maekawa_physical_2006.

In Fig. 2 we compare the different methods to calculate the equation of state
for SrVO3. In the legend we report the estimated values for the equilibrium
volume and the bulk modulus, using a birch-Murnagan third-order function for
the fitting birch_finite_1947. It can be observed that the generalized Hubbard
method does give a better prediction for the bulk-modulus, while
overcorrecting and worsening the value for the equilibrium lattice parameter.
This worsening can be due to the lack of self consistency in our calculations.
In the Supplemental Material supp-info we interpret in terms of charge
transfer the lowering of the bulk modulus.

(THz) |  PBE |  LDA+EDMFT kocer_efficient_2020 |  This work
---|---|---|---
$\omega_{0}$ | -0.01 | 0 | 0.1
$\omega_{1}$ | 4.86 | 4.3 | 3.9
$\omega_{2}$ | 10.3 | 9.0 | 9.4
$\omega_{3}$ | 11.0 | 11.2 | 12.1
$\omega_{4}$ | 17.3 | 18.9 | 19.4
Table 2: $\Gamma$ phonons of SrVO3 with different approaches. Last column are
the calculation with the generalized Hubbard functional introduced in this
work. The data from Ref kocer_efficient_2020 have been obtained digitalizing a
plot.

Using the relaxed structures, and taking advantage of the cheap computational
cost, in Table 2 we calculate the $\Gamma$ phonons for SrVO3 in the
generalized Hubbard framework. Since Hellmann-Feynman does not apply due to
the non-self-consistent nature of the calculations, we evaluate the forces for
the different displacements of the atoms with finite differences. Similarly to
more expensive state-of-the-art LDA+EDMFT kocer_efficient_2020, the
generalized Hubbard approach clearly improves the PBE frequencies.

In summary, in the Letter we provide a Green’s function framework to study
correlated materials, giving results similar to state-of-the-art computational
approaches, like DMFT, at a negligible computational cost. Designing a novel
energy Klein functional, and generalizing the “algorithmic inversion method on
sum over poles” of chiarotti_unified_2022 to crystals, we give an accurate,
closed, and computationally cheap way of calculating accurate spectroscopic
quantities, like particle energies and lifetimes, and thermodynamic quantities
like the bulk-modulus or the phonons. The existence of a Klein functional
guarantees the possibility to perform self-consistent calculation, and thus
the use Hellmann-Feynman theorem to calculate forces in a Green’s function
formalism.

##  I Acknowledgments

We acknowledge stimulating discussions with Mario Caserta and Marco Vanzini.
This work was supported by the Swiss National Science Foundation (SNSF)
through grant No. 200021-179138 (T.C.) and its National Centre of Competence
in Research MARVEL on ‘Computational Design and Discovery of Novel Materials’
(N.M.), and from the EU Commission for the MaX Centre of Excellence on
‘Materials Design at the eXascale’ under grant No. 824143 (N.M., A.F.).