###### Abstract

In this work, we derive a multi-fragment real-time extension of projected
density matrix embedding theory (pDMET) designed to treat non-equilibrium
electron dynamics in strongly correlated systems. As in the previously
developed static pDMET, real time pDMET partitions the total system into many
fragments; the coupling between each fragment and the rest of the system is
treated through a compact representation of the environment in terms of a
quantum bath. Real-time pDMET involves simultaneously propagating the
wavefunctions for each separate fragment-bath embedding system along with an
auxiliary mean-field wavefunction of the total system. The equations of motion
are derived by _(i)_ projecting the time-dependent Schrödinger equation in the
fragment and bath space associated with each separate fragment and by _(ii)_
enforcing the pDMET matching conditions between the global 1-particle reduced
density matrix (1-RDM) obtained from the fragment calculations and the mean-
field 1-RDM at all points in time. The accuracy of the method is benchmarked
through comparisons to time-dependent density-matrix renormalization group
(TD-DMRG) and time-dependent Hartree-Fock (TDHF) theory; the methods were
applied to a one- and two-dimensional single-impurity Anderson model and
multi-impurity Anderson models with ordered and disordered distributions of
the impurities. The results demonstrate a large improvement over TDHF and
rapid convergence to the exact dynamics with an increase in fragment size. Our
results demonstrate that real-time pDMET is a promising and flexible method
that balances accuracy and efficiency to simulate non-equilibrium electron
dynamics in heterogeneous systems of large size.

##  I Introduction

Recent developments in time-resolved spectroscopy have dramatically expanded
the opportunity for real-time observation of correlated electron dynamics in
complex molecular environments at the attosecond scalePeng  _et al._ (2019);
Nisoli  _et al._ (2017); Varillas  _et al._ (2022); Gallmann  _et al._ (2012).
Phenomena of interest include charge transfer following photoexcitation in
biologicalMånsson  _et al._ (2021); de la Cruz Valbuena  _et al._ (2019) and
photovoltaic systemsKrueger  _et al._ (2020); Fuemmeler  _et al._ (2016);
Ponseca  _et al._ (2017), electron and spin transport in molecular
junctionsCui  _et al._ (2019); Joachim  _et al._ (2000); Nitzan and Ratner
(2003); Zhang  _et al._ (2021); Dayen  _et al._ (2020) or quantum dots Singh
_et al._ (2019); Bande (2013), Auger processes in small molecules and
surfaces,Jahnke  _et al._ (2020); Demekhin  _et al._ (2009); Ouchi  _et al._
(2011), among others. New experimental techniques and a wide range of
applications make development of real-time non-equilibrium electronic
structure methods relevant, while bringing the challenge of accurate treatment
of electron correlation and large system size to the forefront.

Some of the methods that are currently used to address this question include
techniques based on non-equilibrium Green’s functionsAoki  _et al._ (2014);
Härtle  _et al._ (2008); Rammer (2018); Freericks  _et al._ (2006), semi-
classical approximationsLi  _et al._ (2013); Swenson  _et al._ (2011, 2012),
real-time Monte Carlo methodsCohen  _et al._ (2014a, b); Gull  _et al._
(2011); Schiró (2010); Werner  _et al._ (2009); Cohen  _et al._ (2015),
numerical path-integral techniquesMühlbacher and Rabani (2008); Weiss  _et
al._ (2009); Segal  _et al._ (2010) and wave function propagation-based
methodsLi  _et al._ (2020a); Goings  _et al._ (2018); Schlünzen  _et al._
(2019); Church and Rubenstein (2021); Thiel (2014); Chatterjee and Makri
(2019). While methods based on frequency domain response theory showed great
success in the calculation of response properties of complex systemsHelgaker
_et al._ (2012); Norman and Dreuw (2018); V. Toukach and P. Ananikov (2013),
real-time approaches show a great advantage in the description of nonlinear
phenomena and provide versatility in the calculated observablesGoings  _et
al._ (2018).

Methods like time dependent density matrix renormalization group (TD-
DMRG)Baiardi and Reiher (2019); Frahm and Pfannkuche (2019); Baiardi (2021);
Li  _et al._ (2020b), configuration interaction (TD-CI)Sonk  _et al._ (2011)
and other complete active space based time-dependent methods Durden and Levine
(2022); Peng  _et al._ (2018); Sato and Ishikawa (2013a); Miranda  _et al._
(2011); Kochetov and Bokarev (2022) achieve high accuracy and have been highly
applicable in the spectroscopic studies of small molecular systems. However,
their practical application is limited by high computational cost. Real-time
time-dependent density functional theory can be a powerful tool that decreases
computational cost associated with the large-scale systems, but the inherent
difficulty of reproducing some of the fundamental effects associated with
approximations of the exchange-correlation potential leaves room for further
improvements in real-time methodsLi  _et al._ (2020a); Isborn  _et al._
(2007).

In this work we address the challenge associated with treatment of a large
system at a high level of accuracy through the use of quantum embedding.
Embedding techniques are based on the highly accurate treatment of a small
region of interest (subsystem) and an approximate, inexpensive calculation of
the surrounding environment. In particular, we extend a specific form of
density matrix embedding theory (DMET)Knizia and Chan (2012, 2013); Wouters
_et al._ (2016a), projected DMET (pDMET)Wu  _et al._ (2019), for the treatment
of non-equilibrium electron dynamics.

Both forms of DMET are based on forming a compact, but approximate description
of the environment in terms of a set of embedding orbitals derived from the
Schmidt decomposition of a mean-field wavefunction for the total system. This
can be performed either for a single subsystem or in a multi-fragment
formulation that involves multiple subsystems each coupled to their own set of
embedding orbitals.Wouters  _et al._ (2016a); Wu  _et al._ (2019) The
difference between conventional DMET and pDMET is the treatment of self-
consistency between the mean-field calculation and the correlated calculations
of the subsystem and embedding orbitals. In the original formulation of DMET,
this self-consistency is achieved by adding the correlation potential, a
single-particle operator, to the Fock operator of the total system; the terms
of the correlation potential are varied to minimize the difference between the
one-particle reduced density matrices (1-RDM) of the mean-field and correlated
calculations.Wouters  _et al._ (2016a) In comparison, the self-consistency
condition is achieved in pDMET by projecting a global 1-RDM, constructed by
averaging over the correlated 1-RDMs of all fragments, onto the set of
idempotent matrices. While the original DMET matching condition provided a
numerical relationship between the mean-field and correlated calculations,
pDMET introduces an analytical correspondence between the global and mean-
field 1-RDMs. DMET has been successfully applied to a variety of strongly
correlated model systemsZheng and Chan (2016); Cui  _et al._ (2020); Faulstich
_et al._ (2022); Reinhard  _et al._ (2019), extended materials Pham  _et al._
(2020); Cui  _et al._ (2020) and spin systemsFan and Jie (2015); Gunst  _et
al._ (2017). It has also been extended to excited state calculations Tran  _et
al._ (2019); Mitra  _et al._ (2021); Qiao and Jie (2021), bond breaking
eventsPham  _et al._ (2018), and electron-phononSandhoefer and Chan (2016) and
long range interactionsRicke  _et al._ (2017).

In the context of dynamics, a real-time expression of DMET has been previously
formulated and was shown to have excellent results in the context of the
single impurity Anderson modelKretchmer and Chan (2018a). The equations of
motion for the embedding wave-function were derived using the time-dependent
variational principle (TDVP). However, the variational nature of the equations
of motion does not include a time-dependent condition that ensures consistency
between fragments, which limits the method to a single-fragment formulation.

In principle, a correlation potential can be introduced in a time-dependent
fashion to enforce self-consistency between the mean-field and correlated
1-RDMs throughout time. However, it was found that $v$-representability issues
almost always occurred after sufficiently long time-propagation, such that the
dynamics became completely independent of the correlation potentialKretchmer
and Chan (2018a).

Therefore, in this work, we take advantage of the analytic matching condition
in pDMET to derive a real-time formulation of pDMET which eliminates the need
for the TDVP and allows for a multi-fragment implementation. This is achieved
by simultaneously propagating the correlated wavefunctions associated with
each fragment and an auxiliary dynamics corresponding to a mean-field dynamics
for the total system, while simultaneously enforcing the pDMET matching
conditions at each point in time. The derived equations of motion provide a
unique real-time method for the simulation of non-equilibrium electron
dynamics that balances accuracy and computational efficiency in extended
systems. Each fragment is small in size, such that the the accurate dynamics
are only solved in a compact space. Furthermore, the dynamics within each
fragment can be treated simultaneously such that the method is parallelizable
to allow for efficient scaling with total system size.

To benchmark the real-time pDMET method, we apply it to the single- and multi-
impurity Anderson model in one and two-dimensions, where the impurities are
distributed both in ordered and disordered fashions. We observe a marked
improvement over time-dependent Hartree-Fock (TDHF) theory and excellent
agreement with the numerically exact TD-DMRG, while only explicitly
correlating a small fraction of the total number of orbitals. This agreement
occurs even in the strongly correlated regime for disordered systems, for
which a single-impurity formulation of the method is poorly suited.

##  II Method

In this section, we begin by first summarizing the static pDMET equations
necessary to provide a foundation for the subsequent presentation of the real-
time formulation; we direct readers to previous work for the full details of
the static pDMET method.Wu  _et al._ (2019) The following notation will be
utilized throughout all sections: the site basis, which is comprised of the
orthonormal real-space localized single-particle basis, will be indexed as
$p,q,r,s$; the embedding basis, which is comprised of the correlated set of
subsystem and associated embedding orbitals, are indexed by $a,b,c$ and $d$;
and the natural orbital basis, which is comprised of the eigenvectors of the
correlated global 1-RDM formed from the 1-RDMs of each fragment, are indexed
by $\mu$ and $\nu$.

###  II.1 Static pDMET

Let us consider a system spanned by an orthonormal single-particle basis,
where $N$ is the number of basis functions and $N_{e}$ is the total number of
electrons in the system. These basis function will be further referred to as
_sites_. The general Hamiltonian of the system is:

$\hat{H}=\sum_{pq}^{N}\sum_{\sigma}h_{pq}\hat{a}^{\dagger}_{p\sigma}\hat{a}_{q\sigma}+\frac{1}{2}\sum_{pqrs}V_{pqrs}\sum_{\sigma\tau}\hat{a}^{\dagger}_{p\sigma}\hat{a}_{r\tau}^{{\dagger}}\hat{a}_{s\tau}\hat{a}_{q\sigma},$
(1)

where $h_{pq}$ and $V_{pqrs}$ are one- and two - electron integrals and
$\hat{a}^{\dagger}_{p\sigma}$ and $\hat{a}_{p\sigma}$ are creation and
annihilation operators associated with site $p$ and spin $\sigma$,
respectively.

In the multi-fragment formulations of DMET or pDMET, the total system is
partitioned into $N_{frag}$ fragments that each consist of $N_{sub}$ subsystem
orbitals. Each fragment is treated as an independent embedding problem where
the orbitals that span the fragment are defined as a subsystem and the rest of
the system is its environment, pictorially illustrated in Figure 1. For each
fragment, the global wavefunction can be compactly written through its Schmidt
decomposition, such that

$|\Psi\rangle=\sum_{\alpha=1}^{N_{A}}\lambda_{\alpha}|A_{\alpha}\rangle|B_{\alpha}\rangle$
(2)

where $|A_{\alpha}\rangle$ are multi-electron states describing the subsystem
A and $|B_{\alpha}\rangle$ are a special set of multi-electron states
describing the environment. Note that the number of states used to describe
the environment is given by $N_{A}$: the same number of multi-electron states
used to describe the subsystem. If the states $|B_{\alpha}\rangle$ are
obtained from the Schmidt decomposition of the exact global wavefunction, they
exactly capture the entanglement with the environment.

In practical DMET calculations, the states $|B_{\alpha}\rangle$ are obtained
from the Schmidt decomposition of a mean-field wavefunction for the total
system, $|\Psi^{mf}\rangle$. The global wavefunction is thus given by a single
Slater determinant and can be described in terms of its single-particle
density matrix,

${\rho}^{mf}_{pq}=\langle\Psi^{mf}|\sum_{\sigma}\hat{a}^{\dagger}_{q\sigma}\hat{a}_{p\sigma}|\Psi^{mf}\rangle.$
(3)

In this case, the states $|B_{\alpha}\rangle$ are comprised of determinants
formed from a set of _single-particle_ embedding orbitals. These embedding
orbitals can be obtained in many mathematically equivalent ways.Wouters  _et
al._ (2016b) Here, we obtain the embedding orbitals associated with fragment
$A$ by diagonalizing the portion of $\rho^{mf}$ that corresponds solely to the
environment of fragment $A$, such that

$\sum_{q\in
env}^{N-N_{sub}}\rho^{mf}_{pq}(R^{emb}_{qa})^{A}=\lambda^{A}_{a}(R^{emb}_{pa})^{A},$
(4)

where both $p$ and $q$ are restricted to the environment of fragment $A$; $p$
and $q$ run over all indices except those that correspond to the subsystem
orbitals of fragment $A$. The $N-N_{sub}\times N-N_{sub}$ matrix
$(R^{emb})^{A}$ corresponds to the matrix of expansion coefficients that
define the embedding orbitals in terms of a linear combination of the site
orbitals associated with the environment of fragment $A$.

The embedding orbitals, $(R^{emb})^{A}$, can be classified into three types
based on their corresponding eigenvalues in Eq. (4). Bath orbitals have
eigenvalues between zero and two and are fully correlated with the subsystem
orbitals. Core orbitals have an eigenvalue equal to two and are fully occupied
in the states $|B_{\alpha}\rangle$. And virtual orbitals have an eigvenvalue
equal to zero and are fully unoccupied in the states $|B_{\alpha}\rangle$.
This implies that the full embedding wavefunction associated with each
fragment is a complete active space (CAS)-like wavefunction composed of the
subsystem and corresponding embedding orbitals.Wouters  _et al._ (2016b)

In this work we use the following notation for the CAS-like embedding
wavefunction associated with fragment $A$:

$|\Psi^{A}\rangle=\sum_{M}C^{A}_{M}|M^{A}\rangle,$ (5)

where $|M\rangle$ denotes a determinant in the embedding space and $C_{M}$ the
corresponding expansion coefficient.

In the static algorithms, the CAS-like embedding wavefunction for each
fragment can be calculated by finding the ground-state of an embedding
Hamiltonian, $H_{emb}^{A}$.Wouters  _et al._ (2016b) The embedding Hamiltonian
is obtained by projecting the original fully interacting Hamiltonian, Eq. (1),
into the CAS space. This projection can be performed by performing a change of
single-particle basis using the following rotation matrix, while including a
contribution from the core orbitals,Wouters  _et al._ (2016b)

$\displaystyle R^{A}=\begin{bmatrix}\mathds{1}_{N_{sub}\times N_{sub}}&0\\\
0&(R^{emb})^{A}\end{bmatrix}$ (6)

Here, the identity matrix denotes that the sub-system orbitals are the same in
the original site basis and the embedding basis; $(R^{emb})^{A}$ is the matrix
of embedding orbitals and is defined in Eq. (4).

Figure 1: Schematic representation of multi-fragmented static DMET. The system
is divided into $N_{frag}$ fragments, each containing $N_{sub}$ subsystem
orbitals (pictured in orange). The environment of each fragment is described
in terms of $N_{sub}$ bath orbitals (pictured in blue) that are correlated
with the subsystem orbitals. These bath orbitals are obtained from an
idempotent density matrix, $\rho^{mf}$, of the total system. An independent
correlated calculation is performed in the subsystem-bath space of each
fragment to form a correlated density matrix, $\rho_{emb}$. The mapping
$\mathcal{F}$ takes in the set of correlated $\rho_{emb}$ and generates a new
idempotent $\rho^{mf}$. The form of the mapping $\mathcal{F}$ is the core
difference between pDMET and previous formulations of the method.

The final step in the DMET algorithm is introducing a self-consistency between
the correlated fragment calculations and the original mean-field calculation
used to generate the embedding orbitals. For the specific case of pDMET this
is accomplished by first constructing a correlated 1-RDM in the global domain,
$\rho^{glob}$, using the 1-RDMs from all separate fragment
calculations.Wouters  _et al._ (2016b); Wu  _et al._ (2019)

$\displaystyle{\rho}^{glob}_{pq}={}$
$\displaystyle\sum_{ab}\frac{1}{2}\left[R_{pa}^{f(p)}{(\rho}_{ab}^{emb})^{f(p)}R^{*f(p)}_{qb}\right.$
(7)
$\displaystyle+\left.R_{pa}^{f(q)}({\rho}^{emb}_{ab})^{f(q)}R^{*f(q)}_{qb}\right]$

The index $f(p)$ corresponds to the fragment that includes site $p$ and
$({\rho}^{emb})^{f(p)}$ corresponds to the 1-RDM obtained from the embedding
wavefunction for fragment $f(p)$,

$\displaystyle({\rho}^{emb}_{ab})^{f(p)}$ $\displaystyle=$
$\displaystyle\langle\Psi^{f(p)}|\sum_{\sigma}\hat{a}^{\dagger}_{b\sigma}\hat{a}_{a\sigma}|\Psi^{f(p)}\rangle.$
(8)

The global 1-RDM, Eq. (7), is projected onto the space of idempotent matrices
with $N_{e}$ electrons to form a new guess for the mean-field wavefunction of
the total system.Wu  _et al._ (2019) In the spin restricted formulation, which
is utilized in this work, the corresponding idempotent matrix is formed from
the $N_{occ}=N_{e}/2$ eigenvectors of $\rho^{glob}$ with the highest
occupation. These eigenvectors of $\rho^{glob}$ are referred to as natural
orbitals, represented by the matrix $U$, and indexed by $\mu$ and $\nu$, such
that

$\sum_{q}^{N}\rho^{glob}_{pq}U_{q\mu}=\lambda_{\mu}U_{q\mu}$ (9)

and

$\rho^{mf}_{pq}=2\sum_{\mu}^{N_{occ}}U_{p\mu}U_{q\mu}^{*}$ (10)

Given the newly formed $\rho^{mf}$, the pDMET algorithm is iterated until
$\rho^{mf}$ does not change between iterations.

The alternative self-consistency procedure used in pDMET compared to
conventional DMET results in a significant reduction of the computational cost
for a large non-periodic lattice and provides a more robust analytical
relationship between the embedding and mean field representations of the
problem. The latter point is crucial for the derivation of a real-time multi-
fragmented DMET algorithm.

###  II.2 Real-time pDMET

We now describe the central methodological contribution of the work, a multi-
fragment real-time extension of pDMET. This is achieved through the
simultaneous propagation of _(i)_ all the CAS-like wavefunctions for each
fragment along with _(ii)_ an additional mean-field dynamics for the total
system. The equations of motions are derived to maintain the inherit
relationship between the mean-field and global description of the system. This
is achieved by preserving the pDMET self-consistencies that connect the mean-
field and embedding representations of the system at each point in time. The
full derivations of the working dynamical equations defined in this section
can be found in the Supplemental Information (SI).

We begin by first deriving the formal time-dependence of the embedding
wavefunction of fragment $A$ by projecting the time-dependent Schrödinger
equation into the embedding space,

$\displaystyle i|\dot{\Psi}^{A}\rangle$ $\displaystyle=$
$\displaystyle\hat{H}|\Psi^{A}\rangle$ (11) $\displaystyle
iP^{A}|\dot{\Psi}^{A}\rangle$ $\displaystyle=$ $\displaystyle
P^{A}\hat{H}P^{A}|\Psi^{A}\rangle$ (12) $\displaystyle
i\frac{\partial}{\partial
t}\left(P^{A}|\Psi^{A}\rangle\right)-i\dot{P}^{A}|\Psi\rangle$
$\displaystyle=$ $\displaystyle\hat{H}^{emb}|\Psi^{A}\rangle$ (13)
$\displaystyle i|\dot{\Psi}^{A}\rangle=\hat{H}^{emb}|\Psi^{A}\rangle$
$\displaystyle+$ $\displaystyle i\dot{P}^{A}|\Psi^{A}\rangle$ (14)

where $P^{A}$ is the formal projector into the embedding space of fragment
$A$, $H^{emb}=P^{A}HP^{A}$ is the Hamiltonian projected into the embedding
space, and we have used that $P^{A}|\Psi^{A}\rangle=|\Psi^{A}\rangle$Wouters
_et al._ (2016b). Eq. (14) indicates that the time-dependence of the embedding
wavefunction has two components, the first associated with the embedding
Hamiltonian and the second associated with the time-dependence of the quantum
bath given by $\dot{P}^{A}$.

For the case of a mean-field global wavefunction, $|\Psi^{A}\rangle$ is given
by Eq. (5), and the projector is $P=\sum_{M}|M^{A}\rangle\langle M^{A}|$.
Therefore, the time-dependence of $|\Psi^{A}\rangle$ is described in terms of
the time-dependence of the expansion coefficients, $C^{A}$, and determinants,
$|M^{A}\rangle$,

$\displaystyle i|\dot{\Psi}^{A}\rangle$ $\displaystyle=$
$\displaystyle\sum_{M}i\dot{C}^{A}_{M}|M^{A}\rangle+C^{A}_{M}i|\dot{M}^{A}\rangle$
(15) $\displaystyle=$
$\displaystyle\sum_{M}i\dot{C}^{A}_{M}|M^{A}\rangle+C^{A}_{M}\hat{X}^{A}|M^{A}\rangle,$
(16)

and the time-dependence of the projector is given by

$\displaystyle i\dot{P^{A}}$ $\displaystyle=$
$\displaystyle\hat{X}^{A}P^{A}-P^{A}\hat{X}^{A}.$ (17)

Here, we have introduced the single-particle operator $\hat{X}$ that governs
the time-dependence of the embedding orbitals that comprise the determinants,
such that

$i\dot{R}^{A}_{pa}=\sum_{q}X^{A}_{pq}R^{A}_{qa}.$ (18)

This operator is analogous to the operator employed in time-dependent CASSCF
(TD-CASSCF) that propagates the orbitals of the CASSCF-wavefunction.Sato and
Ishikawa (2013b)

Equation (16) indicates that the dynamics of the embedding wavefunctions are
given by the equations of motion for the expansion coefficients $C^{A}$ and
the matrix elements of the operator $\hat{X}^{A}$. The equations of motion of
the coefficients $C^{A}$ are derived by combining Eqs. (14), (16), (17) and
multiplying by $\langle N^{A}|$ yielding

$i\dot{C}^{A}_{N}=\sum_{M}(H^{emb}_{NM}-X^{A}_{NM})C^{A}_{M}.$ (19)

This equation effectively corresponds to the time-dependent full CI equations
of motion in the compact embedding space. In principle, further computational
efficiency could be obtained by restricting the number of allowed excitations
in $\hat{H}^{emb}$ and $\hat{X}^{A}$, albeit at decreased accuracy.

We now turn our attention to the derivation of the matrix elements of the
operator $\hat{X}^{A}$. First, in order to maintain the original definition of
the fragments $A$, the subsystem orbitals are defined as time-independent.
Therefore, $X^{A}_{pq}=0$ if orbitals $p\textrm{ or }q\in A$. Non-subsystem
matrix elements of the operator $\hat{X}^{A}$ are derived by taking the time-
derivative of the pDMET condition in Eq. (4), which yields

$X_{ba}^{A}=\sum_{pq}\frac{R_{pb}^{A*}i\dot{\rho}^{mf}_{pq}R_{qa}^{A}}{\lambda_{a}^{A}-\lambda_{b}^{A}}$
(20)

This matrix can be interpreted as a rotation of all embedding pDMET orbitals
in time, which allows for a few simplifications. Analogous to TD-CASSCF, any
intraspace rotations among embedding orbitals of the same type are
redundant.Sato and Ishikawa (2013b) We can therefore choose the elements of
$\hat{X}^{A}$ for these rotations. We make the common choice to set
$X^{A}_{ab}=0$ if $a$ and $b$ both correspond to either core orbitals or
virtual orbitals. However, unlike TD-CASSCF, the bath-bath rotations are not
arbitrary, and are given by Eq. (20), because we have the additional
constraint that the embedding orbitals must be eigenvectors of the mean-field
1-RDM as given in Eq. (4); this is not an issue for the core and virtual space
as those orbitals are degenerate. Additionally, though, we are able to set
$X^{A}_{ab}=0$ for $a=b$ even for the bath-orbitals as the diagonal elements
of $\hat{X}^{A}$ only affect the norm of the orbitals.

The time dependence of $\rho^{mf}$ referenced in equation (20) is derived by
enforcing the pDMET matching condition between the natural orbitals of the
$\rho^{glob}$ and $\rho^{mf}$, Eq. (10). To derive the time-dependence of
$\rho^{mf}$ we first introduce a single-particle operator, $\hat{G}$ that
governs the time-dependence of the natural orbitals:

$i\dot{U}_{p\mu}=\sum_{q}G_{pq}U_{q\mu}.$ (21)

The time-dependence of $\rho^{mf}$ can be obtained by taking the time-
derivative of Eq. (10) and using Eq. (21), yielding

$i\dot{\rho}^{mf}_{pq}=\sum_{r}[G_{pr}\rho^{mf}_{rq}-\rho^{mf}_{pr}G_{rp}].$
(22)

The matrix-elements of $\hat{G}$ are derived from the relationship between the
global 1-RDM and the natural orbitals, Eq. (9). Differentiating Eq. (9) with
respect to time and solving for the matrix-elements of $\hat{G}$ yields

$G_{\mu\nu}=\sum_{pq}\frac{U_{p\mu}^{*}i\dot{\rho}^{glob}_{pq}U_{q\nu}}{\lambda_{\nu}-\lambda_{\mu}}$
(23)

where $G_{\mu\nu}$ is defined in the natural orbital basis. The matrix
elements in the site-basis utilized in Eq. (22) can be obtained by a single-
particle rotation using the natural orbitals, $U$, as the rotation matrix. The
matrix elements of $\hat{G}$ take on a similar form as $\hat{X}$ due to the
analogous structure of Eqs. (4) and (9).

Analogous to the definition of $X^{A}_{ba}$, $G_{\mu\nu}=0$ for
$\lambda_{\mu}=\lambda_{\nu}$. Additionally, a regularization parameter
$\varepsilon_{G}$ is introduced such that $G_{\mu\nu}=0$ if
$|\lambda_{\mu}-\lambda_{\nu}|<\varepsilon_{G}$ to avoid numerical
instabilities in the dynamics. Similar parameters have been utilized in the
context of TD-CASSCF and MCTDH.Kretchmer and Chan (2018a, b); Manthe  _et al._
(1992); Meyer  _et al._ (1990)

The mean-field dynamics retain an implicit dependence on the time-evolution of
$\rho^{glob}$ through the form of $\hat{G}$, Eq. (23). The dynamics of
$\rho^{glob}$ are governed by the correlated dynamics of each fragment through
differentiation of Eq. (7),

$\displaystyle i\dot{\rho}^{glob}_{pq}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{abc}\left[{R}_{pc}^{f(p)}X_{ca}^{f(p)}(\rho_{ab}^{emb})^{f(p)}R^{*f(p)}_{qb}-{R}_{pa}^{f(p)}(\rho_{ab}^{emb})^{f(p)}X_{bc}^{f(p)}{R}^{*f(p)}_{qc}\right.$
(24)
$\displaystyle+\left.{R}_{pc}^{f(q)}X_{ca}^{f(q)}(\rho_{ab}^{emb})^{f(q)}R^{*f(q)}_{qb}-{R}_{pa}^{f(q)}(\rho_{ab}^{emb})^{f(q)}X_{bc}^{q}{R}^{*f(q)}_{qc}\right]$
$\displaystyle+\frac{1}{2}\sum_{ab}\left[R_{pa}^{f(p)}i(\dot{\rho}_{ab}^{emb})^{f(p)}R^{*f(p)}_{qb}+R_{pa}^{f(q)}i(\dot{\rho}_{ab}^{emb})^{f(q)}R^{*f(q)}_{qb}\right].$

Here we have used Eq. (18).

The time-dependence of $(\rho^{emb})^{A}$ is obtained by _(i)_ differentiating
Eq. (8) and _(ii)_ simplifying the expression based on the form of the 1-RDM
and 2-RDM for the CAS-like fragment wavefunction, analogously to the
derivation of the TD-CASSCF and static CASSCF theories.Sato and Ishikawa
(2013b); Helgaker  _et al._ (2014) The final expression for the time-
dependence of $(\rho^{emb})^{A}$ is given by

$\displaystyle i(\dot{\rho}^{emb}_{ab})^{A}$ $\displaystyle=$
$\displaystyle(F_{ba}-F_{ab}^{*})^{A}$ (25) $\displaystyle+$
$\displaystyle\sum_{c}\left((\rho_{ac}^{emb})^{A}X_{cb}^{A}-X_{ac}^{A}(\rho^{emb}_{cb})^{A}\right)$

where the matrix $F$ is the generalized Fock matrix.Helgaker  _et al._ (2014)

Plugging Eq. (25) into Eq. (24) and simplifying yields the final expression
for $\dot{\rho}^{glob}$,

$\displaystyle i\dot{\rho}^{glob}_{pq}$ $\displaystyle=$
$\displaystyle\sum_{ab}\frac{1}{2}\left[R_{pa}^{f(p)}(F_{ba}-F_{ab}^{*})^{f(p)}R^{f(p)}_{qb}\right.$
(26)
$\displaystyle+\left.R_{pa}^{f(q)}(F_{ba}-F_{ab}^{*})^{f(q)}R^{f(q)}_{qb}\right].$

Interestingly, the dependence on $R$ in Eqs. (24) and (25) cancel out such
that $\dot{\rho}^{glob}$ is independent of $\hat{X}$. This is an important
result because it means that the real-time pDMET equations of motion do not
involve an inversion process, which can be numerically unstable and expensive.

This methodology presents a robust set of equations of motion that describe
the coupled time evolution of the system at the mean-field and correlated
levels. In principle, it is only necessary to directly propagate the CAS
coefficients, $C_{M}^{A}$, for all fragments and the mean-field density
matrix, $\rho^{mf}$, using Eqs. (19) and (22); all other time-dependent
quantities can be calculated once these are known. However, to increase the
numerical stability of the method, we additionally directly propagate the
embedding orbitals, $R^{A}$, for all fragments, the natural orbitals, $U$, and
the global density matrix $\rho^{glob}$ using Eqs. (18), (21), and (26),
respectively. Full algorithmic details are provided in the SI. The working
equations of motion have been implemented using in-house code that has been
interfaced with the PySCF packageSun  _et al._ (2020). Publicly available
version of the code has been made available rtd .

##  III Results

The performance of real time pDMET is examined by simulating electron dynamics
in a set of 1-dimensional (1D) multi-impurity Anderson models (MIAMs) and
single-impurity Anderson models (SIAMs). The first set of calculations involve
suddenly changing the strength of the electron-electron interaction terms,
while a second set of calculations additionally include a laser pulse. We
conclude with a calculation on a 2-dimensional (2D) SIAM.

The SIAM models a system that consists of a single impurity with a local
Coulomb interaction ($U$) and a nearest-neighbor hopping between all sites
($t$).

The MIAM extends this model to include multiple impurity sites, $N_{imp}$,
that each has a local Coulomb interaction. The Hamiltonian that encompasses
both the SIAM and MIAM is given by

$H=-t\sum_{<ij>\sigma}(a_{i\sigma}^{\dagger}a_{j\sigma}+h.c.)+U\sum_{i\in
imp}n_{i\uparrow}n_{i\downarrow}$ (27)

The Coulomb interaction is described in terms of the occupation number
operator $n_{i,\sigma}=a_{i\sigma}^{\dagger}a_{i\sigma}$ on impurity site $i$.
The sum $\sum\limits_{<ij>}$ implies a sum over only nearest-neighbors. These
models provide a useful and challenging benchmark platform for evaluating the
feasibility and effectiveness of the real-time pDMET methodTitvinidze  _et
al._ (2015); Eickhoff and Anders (2020).

To generate non-equilibrium dynamics, we initialize the system in the non-
interacting ground state of the SIAM/MIAM ($U=0$, $t=1$) by performing a
static pDMET calculation using exact diagonalization to solve the embedding
problem and subsequently propagate the dynamics using a Hamiltonian with a
non-zero value of the interaction parameter ($U^{\prime}\neq U$). Initializing
the system at $U=0$ provides a useful benchmark for real-time pDMET because
pDMET is exact in the non-interacting limit. To examine the accuracy of the
dynamics, we benchmark the performance of real-time pDMET against TD-DMRG and
TDHF theory. A discussion of the numerical implementation and full
computational details are provided in the SI.

Figure 2: The time dependent occupancy on the single impurity, $n_{d}(t)$, for
the SIAM for which the strength of the Coulomb interaction is switched from
$U=0$ to (a) $U^{\prime}=2$, (b) $U^{\prime}=3$ and (c) $U^{\prime}=4$. Real-
time pDMET results for $N_{sub}=2$ (orange), $N_{sub}=4$ (blue), $N_{sub}=5$
(pink) are compared to TD-DMRG (black) and TDHF (green). The total size of the
system is $N=60$ with the single impurity located at $d=N/2$.

Figure 2 presents results for a SIAM consisting of $N=60$ sites in which the
initial Coulomb interaction is changed from $U=0$ to $U^{\prime}=$ (a) 2, (b)
3, and (c) 4. Firstly, Figure 2 shows that TDHF is able to reasonably capture
the dynamics when correlation effects are weak ($U^{\prime}=2$), but rapidly
declines in accuracy as the importance of correlation effects increases with
higher $U^{\prime}$. In contrast, real-time pDMET always converges to the TD-
DMRG results with increasing fragment size regardless of the magnitude of
$U^{\prime}$. While $N_{frag}=2$ is insufficient at higher $U^{\prime}$,
strongly correlated dynamics are fully reproduced with a fragment size as
small as $N_{frag}=5$. This implies that a reduced active space of 5 subsystem
and 5 bath orbitals for each fragment is sufficient to capture the correlated
electron dynamics of this $N=60$ system.

Figure 3: The time dependent occupancy on one of the 6 impurities, $n_{d}(t)$,
for the MIAM with (a) evenly distributed impurities and $d=65$, (b) unevenly
distributed impurities and $d=67$, and (c) unevenly distributed impurities and
$d=45$. Real-time pDMET results for $N_{sub}=5$ (pink) are compared to TD-DMRG
(black) and TDHF (green). In all panels the strength of the Coulomb
interaction is switched from $U=0$ to $U^{\prime}=3$ and the size of the
system is $N=120$.

Figure 3 presents results for the more challenging case of the MIAM in which
the Coulomb interaction is changed from $U=0$ to $U^{\prime}=3$. Figure 3(a)
corresponds to a system in which the impurities are evenly distributed
throughout the MIAM, while (b) and (c) present results for the time-dependent
occupancy on different impurity sites in a MIAM with unevenly distributed
impurities (Figure 4).

Analogous to the SIAM, TDHF fails to reproduce the dynamics in these strongly
correlated systems, while real-time pDMET shows excellent agreement with the
TD-DMRG results. The similar accuracy in parts (b) and (c) indicate that the
multi-fragment formulation allows for an even treatment of the entire system.
The results for the unevenly distributed MIAM are especially encouraging as
they indicate the ability for real-time pDMET to treat complex dynamics even
in disordered systems.

Figure 4: Pictorial representation of the multi-impurity Anderson model with
an (a) even and (b) uneven distribution of impurities (dark grey). Dynamics
are initialized by suddenly changing the local Coulomb interaction from $U=0$
to a finite $U^{\prime}$ only on the impurity sites; the interaction term
remains $U=0$ on all other sites (white).

We now turn our attention to results that additionally include an ultra-fast
laser pulse to examine the behavior of rt-pDMET under a more complex non-
equilibrium setting. Figure 5 presents the time-dependent occupancy for the
analogous $N=60$ SIAM system presented in Figure 2(b), but now including a
laser pulse that is restricted to interact with only the central 20 sites of
the system. The laser induces a time-dependence to the hopping terms of the
Hamiltonian associated with these sites as described in detail in the SI. In
addition to the laser pulse, we still include a sudden change from $U=0$ to
$U^{\prime}=3$ at $t=0$, such that the dynamics are initialized from an exact
solution.

Figure 5(a) shows the occupancy on the impurity site and is comparable to
Figure 2(b); the dynamics on this site are dominated by the sudden change in
$U$ versus the interaction with the laser pulse. In comparison, Figure 5(b)
shows the occupancy on a site removed from the central impurity but within the
region interacting with the laser pulse. Here, real-time pDMET (pink)
accurately captures the distinct change in occupancy of the laser-exposed site
at the pumping time at $t=5$, which can not be described by TDHF (green) even
on a qualitative level. Lastly, Figure 5(c) presents the occupancy on a site
outside of the interaction region with the laser pulse. Real-time pDMET shows
excellent agreement with TD-DMRG in capturing the delayed dynamics that don’t
appear until $t\sim 7.5$ associated with electron flow from the central
impurity site and the region interacting with the laser.

Figure 5:  The time-dependent occupancy on various sites for a SIAM exposed to
a laser pulse and a sudden change in the Coulomb interaction. (a) Site $d=30$,
which corresponds to the impurity site containing the on-site Coulomb
interaction. (b) Site $d=21$ located within the region exposed to the laser.
(c) Site $d=15$ positioned away from both the laser region and the impurity
site. Real-time pDMET results for $N_{sub}=5$ (pink) are compared to TD-DMRG
(black) and TDHF (green). Figure 6:  The time dependent occupancy on sites
within 2D SIAM of total size $N=64$, where electron dynamics are driven by
Coulomb interaction change from $U=0$ to $U^{\prime}=3$. (a) Site $d=(3,4)$,
which corresponds to the impurity site containing the on-site Coulomb
interaction. TDHF (green) results are compared to the real-time pDMET results
for three fragment sizes: $N_{sub}=1\times 1$ (yellow), $N_{sub}=1\times 2$
(orange), $N_{sub}=1\times 4$ (blue). (b) Two-dimensional electron density
distribution across the full system at four different points of time during
the simulation.

We conclude by illustrating the flexibility of real-time pDMET towards more
realistic systems beyond what can be easily treated by TD-DMRG. Specifically,
we examine the dynamics in a 2D SIAM initialized by a sudden switch in the
value of $U$ from 0 to 3 on one of the central sites. TD-DMRG is poorly suited
to treat multi-dimensional systems due to the structure of the matrix-product
state _Ansätze_ for the wavefunction.Baiardi and Reiher (2020) Figure 6(a)
shows the time-dependent occupation on the impurity that contains the on-site
Coulomb interaction. Real-time pDMET demonstrates a rapid convergence of the
dynamics with an increase in fragment size. The oscillatory behaviour observed
in the $N_{sub}=1\times 1$ results (yellow) is eliminated by a minimum
increase in fragment size to $N_{sub}=1\times 2$ (orange). Sequential increase
of the fragment size to $N_{sub}=1\times 4$ (blue) demonstrates convergence of
the dynamics that is distinct from the TDHF results (green).

Figure 6(b) shows the electron density in the full system at four
representative times throughout the simulation. The rapid change in the
occupation of the site with the Coulomb interaction is observed immediately
after the change in the Hamiltonian. Sequential redistribution of electron
density is observed from the upper left corner to the bottom right corner of
the system as a function of time associated with a response to the sudden
change in $U$ and reflections off of the boundaries of the system. This
electron flow is what leads to the recurrence peak around $t=8$ in Figure
6(a).

##  IV Conclusion

In this work, we present the development of a real-time extension of projected
density matrix embedding theory (real-time pDMET). The analytic matching
condition used to obtain self-consistency in static pDMET compared to the
correlation potential utilized in conventional DMET removes previous
theoretical limitations associated with real-time extensions of DMET. This
allows for a multi-fragmented formulation of the dynamics, analogous to static
pDMET and DMET.

The real-time pDMET method was benchmarked for non-equilibrium electron
dynamics in both single-impurity and multiple-impurity 1- and 2-dimensional
Anderson models (SIAM and MIAM). Dynamics were initialized through a sudden
change in the electron interaction term and through an external laser pulse.
In all cases, real-time pDMET converged rapidly to numerically exact TD-DMRG
results with increasing fragment size. Most notably real-time pDMET was able
to accurately capture the dynamics across all regions of space in a disordered
MIAM and exhibited rapid convergence of the dynamics in a 2D SIAM, for which
TD-DMRG is poorly suited.

Overall, the multi-fragmented formulation of real-time pDMET provides an
efficient and flexible framework to investigate time-dependent properties of
strongly correlated systems. The embedding framework is not limited by the
dimensionality of the system and allows for the future incorporation of other
wavefunction _Ansätze_ , such as coupled-cluster or matrix-product state
wavefunctions, as has been done in the context of static calculations. Future
development will involve an unrestricted formulation of the method, which will
allow for higher accuracy and explicit treatment of spin dynamics.

##  V Supplementary Material

The supplementary material contains a full derivation of the working dynamical
equations for the real-time pDMET method along with a discussion of the
numerical implementation and full computational details. Publicly available
version of the real-time pDMET code used to generate the results in this
manuscripts has been made availablertd .

###### Acknowledgements.

This work was partially supported by a Georgia Tech Seed grant. Computational
resources were supported in part through research cyberinfrastructure
resources and services provided by the Partnership for an Advanced Computing
Environment (PACE) at the Georgia Institute of Technology, Atlanta, Georgia,
USA.

Supporting Information

A multi-fragment real-time extension of projected density matrix embedding
theory: Non-equilibrium electron dynamics in extended systems

D. Yehorova, Joshua S. Kretchmer

## Detailed derivations of the equations of motion

Here we provide complete derivations for the real-time pDMET equations of
motion defined in the main text. The notation used throughout the following
derivations is consistent with the notation used in the methods section of
main text. The site basis is indexed by $p,q,r,s$; the embedding basis is
indexed by $a,b,c$ and $d$; and the natural orbital basis is indexed by $\mu$
and $\nu$.

Two-electron integrals are defined using chemist notation,

$V_{pqrs}=\int\int\phi_{p}^{*}(r_{1})\phi_{r}^{*}(r_{2})\frac{1}{r_{12}}\phi_{q}(r_{1})\phi_{s}(r_{2})dr_{1}dr_{2}.$
(28)

The 1- and 2-particle reduced density matrices are defined as:

$\displaystyle{\rho}_{pq}$ $\displaystyle=$
$\displaystyle\langle\Psi|\hat{a}^{\dagger}_{q\sigma}\hat{a}_{p\sigma}|\Psi\rangle$
(29) $\displaystyle\Gamma_{pqrs}$ $\displaystyle=$
$\displaystyle\sum_{\sigma\tau}\langle\Psi|a_{p\sigma}^{{\dagger}}a_{r\tau}^{{\dagger}}a_{s\tau}a_{q\sigma}|\Psi\rangle,$
(30)

respectively. Both the embedding and natural orbitals define a rotation
between the site and the corresponding basis.

$\displaystyle|a\rangle=\sum_{p}R_{pa}|p\rangle,$ $\displaystyle
R_{pa}=\langle p|a\rangle$ (31)
$\displaystyle|\mu\rangle=\sum_{p}U_{p\mu}|p\rangle,$ $\displaystyle
U_{p\mu}=\langle p|\mu\rangle$ (32)

###  V.1 Derivation of $X$ and $G$ matrices

The operators $\hat{X}$ and $\hat{G}$ describe time evolution of the embedding
and natural orbitals respectively. The following discussion focuses on the
derivation of the matrix elements of $X$, Eq. (20) in the main text; the
derivation of the matrix elements of $G$, Eq. (23) in the main text, follows
completely analogous steps.

The derivation of $X$ is based on two previously described conditions: the
pDMET condition that defines the relationship between $\rho^{mf}$ and $R$ and
the initial definition of $X$ as a propagator of the embedding orbitals.

$\displaystyle\sum_{q\in env}^{N-N_{A}}\rho^{mf}_{pq}R^{A}_{qa}$
$\displaystyle=$ $\displaystyle\lambda^{A}_{a}R^{A}_{pa}$ (33) $\displaystyle
i\dot{R}_{pa}^{A}$ $\displaystyle=$
$\displaystyle\sum_{q}X_{pq}^{A}R_{qa}^{A}$ (34)

Each fragment, $A$, has a corresponding set of embedding orbitals, which
implies that there is a unique $\hat{X}^{A}$ for each fragment. However, the
explicit reference to fragment $A$ will be suppressed in the following
derivation for convenience.

The derivation of the non-subsystem elements of $X$ begins by taking the time-
derivative of Eq. (33):

$\sum_{q\in
env}\left(\dot{\rho}^{mf}_{pq}R_{qa}+\rho^{mf}_{pq}\dot{R}_{qa}\right)=\dot{\lambda}_{a}R_{pa}+\lambda_{a}\dot{R}_{pa}$
(35)

The right side of the equation can be rewritten by explicitly defining
$\dot{\lambda}_{a}$ using Eq. (33):

$\displaystyle\lambda_{a}$ $\displaystyle=$ $\displaystyle\sum_{qr\in
env}R_{qa}^{*}\rho^{mf}_{qr}R_{ra}$ (36) $\displaystyle\dot{\lambda}_{a}$
$\displaystyle=$ $\displaystyle\sum_{qr\in
env}\left[\dot{R}_{qa}^{*}\rho^{mf}_{qr}R_{ra}+R_{qa}^{*}\rho^{mf}_{qr}\dot{R}_{ra}+R_{qa}^{*}\dot{\rho}^{mf}_{qr}R_{ra}\right]$
(37)

Using the definition of the matrix $X$ stated in Eq. (34), Eq. (37), and the
pDMET condition Eq. (33) yields

$i\dot{\lambda}_{a}=\sum_{q\in
env}\sum_{s}\left(-R^{*}_{sa}X_{sq}\lambda_{a}R_{qa}\right)+\sum_{r\in
env}\sum_{s}\left(\lambda_{a}R_{ra}^{*}X_{rs}R_{sa}\right)+\sum_{qr\in
env}\left(R_{qa}^{*}i\dot{\rho}_{qr}^{mf}R_{ra}\right)\\\ $ (38)

This expression can be further simplified by realizing that the subsystem
orbitals are time-independent; this condition implies that $X_{sq}=0$ if
either $s$ or $r$ is a subsystem orbital of the fragment. Therefore, both sums
can be restricted to only the environment orbitals. This simplification allows
for the rearrangement of indices such that the first two terms cancel,
reducing the definition of $i\dot{\lambda}_{a}$ in Eq. (38) to only the last
term.

The matrix $X$ is similarly introduced into the left side of Eq. (35) by using
Eqs. (33) and (37), such that Eq. (35) becomes

$\sum_{q\in env}\dot{\rho}^{mf}_{pq}R_{qa}+\sum_{q\in
env}\sum_{r}\rho^{mf}X_{qr}R_{ra}=\sum_{qr\in
env}\left(R_{qa}^{*}i\dot{\rho}_{qr}^{mf}R_{ra}\right)R_{pa}+\lambda_{a}\sum_{r}X_{pr}R_{ra}$
(39)

Multiplying on the left by $\langle b|$ using the definition (31), and then
further simplifying using the complex conjugate of Eq. 33, the fact that
$\rho^{mf}$ is Hermitian, $\lambda_{a}$ is real, and the unitarity of the
matrix $R$ yields

$\displaystyle\sum_{pq\in env}R_{pb}^{*}i\dot{\rho}^{mf}_{pq}R_{qa}+\sum_{q\in
env}\sum_{r}\lambda_{b}R_{pb}^{*}X_{qr}R_{ra}=$ (40) $\displaystyle\sum_{qr\in
env}\left(R_{qa}^{*}i\dot{\rho}_{qr}^{mf}R_{ra}\right)\delta_{ba}+\sum_{p\in
env}\sum_{r}\lambda_{a}R_{pb}^{*}X_{pr}R_{ra}$

All summations can now be restricted to the environment space by an analogous
argument used in Eq. (38), such that

$\displaystyle\sum_{pq\in
env}R_{pb}^{*}i\dot{\rho}^{mf}_{pq}R_{qa}+\sum_{qr\in
env}\lambda_{b}R_{pb}^{*}X_{qr}R_{ra}=$ $\displaystyle\sum_{qr\in
env}\left(R_{qa}^{*}i\dot{\rho}_{qr}^{mf}R_{ra}\right)\delta_{ba}+\sum_{pr\in
env}\lambda_{a}R_{pb}^{*}X_{pr}R_{ra}$ (41)

This equation can be simplified using that the matrix $X$ in the embedding
space is defined as

$\displaystyle X_{ba}$ $\displaystyle=$ $\displaystyle\sum_{pr\in
env}R_{pb}^{*}X_{pr}R_{ra}$ (42)

and that due to the redundancy of intraspace rotations, the nonzero terms of
$X$ correspond to $a\neq b$. Therefore, the final expression for the nonzero
elements correspond to

$X_{ba}=\sum_{pq\in
env}\frac{R_{pb}^{*}i\dot{\rho}^{mf}_{pq}R_{qa}}{\lambda_{a}-\lambda_{b}},$
(43)

which corresponds to Eq. (20) in the main text.

The derivation of the operator governing the time-dependence for the natural
orbital is completely analogous where equations (33) and (34) are replaced by
the corresponding relationships in natural orbital space and all operations
are carried out in the global domain.

$\displaystyle\sum_{q}\rho^{glob}_{pq}U_{q\mu}$ $\displaystyle=$
$\displaystyle\lambda_{\mu}U_{p\mu}$ (44) $\displaystyle i\dot{U}_{p\mu}$
$\displaystyle=$ $\displaystyle\sum_{q}G_{pq}U_{q\mu}$ (45)

The final expression of the matrix elements of the operator $G$ is:

$G_{\mu\nu}=\sum_{pq}\frac{U_{p\mu}^{*}i\dot{\rho}^{glob}_{pq}U_{q\nu}}{\lambda_{\nu}-\lambda_{\mu}},$
(46)

which corresponds to Eq. (23) in the main text.

###  V.2 Derivation of $\dot{\rho}^{mf}$

The time-dependence of $\rho^{mf}$ is directly obtained from the dynamics of
the natural orbitals $U$. This is observed after taking the derivative of the
pDMET condition (47) which ensures the consistency between the mean-field
representation and global domain of the system.

$\displaystyle\rho^{mf}_{pq}$ $\displaystyle=$ $\displaystyle
2\sum_{\mu}^{N_{occ}}U_{p\mu}U_{q\mu}^{*}$ (47) $\displaystyle
i\dot{\rho}^{mf}_{pq}$ $\displaystyle=$ $\displaystyle
2\sum_{\mu}^{N_{occ}}\left(i\dot{U}_{p\mu}U_{q\mu}^{*}+U_{p\mu}i\dot{U}_{q\mu}^{*}\right)$
(48)

The time dependence of $U$ is defined in terms of the operator $\hat{G}$, Eq.
(45). Therefore, Eq. (48) can be rewritten as a commutator of $G$ and
$\rho^{mf}$ in the site basis,

$\displaystyle i\dot{\rho}^{mf}_{pq}$ $\displaystyle=$ $\displaystyle
2\sum_{\mu}^{N_{occ}}\sum_{r}\left(G_{pr}{U}_{r\mu}U_{q\mu}^{*}-U_{p\mu}{U}_{r\mu}^{*}G_{rp}\right)$
(49) $\displaystyle=$
$\displaystyle\sum_{r}\left(G_{pr}\rho^{mf}_{rq}-\rho^{mf}_{pr}G_{rp}\right),$

which corresponds to Eq. (22) in the main text.

###  V.3 Derivation of $\dot{\rho}^{emb}$

The time-dependence of the correlated 1-RDM associated with each fragment,
$(\rho_{ab}^{emb})^{A}$, is necessary for the derivation of
$\dot{\rho}^{glob}$ seen below. This 1-RDM is obtained from the fragment
wavefunction as

$(\rho_{ab}^{emb})^{A}=\langle\Psi^{A}|\sum_{\sigma}a^{{\dagger}A}_{b\sigma}a^{A}_{a\sigma}|\Psi^{A}\rangle$
(50)

The time derivative of this expression results in 4 terms, where the explicit
dependence on the fragment $A$ is suppressed for convenience

$\displaystyle
i\dot{\rho}_{ab}=i\langle\dot{\Psi}|\sum_{\sigma}a^{{\dagger}}_{b\sigma}a_{a\sigma}|\Psi\rangle+i\langle\Psi|\sum_{\sigma}a^{{\dagger}}_{b\sigma}a_{a\sigma}|\dot{\Psi}\rangle$
(51)
$\displaystyle+\langle{\Psi}|\sum_{\sigma}i\dot{a}^{{\dagger}}_{b\sigma}a_{a\sigma}|\Psi\rangle+\langle{\Psi}|\sum_{\sigma}ia^{{\dagger}}_{b\sigma}\dot{a}_{a\sigma}|\Psi\rangle$

The first two terms are redefined as a commutator using the time-dependent
Schrödinger equation (TDSE); the second two terms can be written in terms of
the time-dependence of the embedding orbitals by introducing the
transformation of the creation and annihilation operators from the site to the
embedding basis,

$\displaystyle a_{a\sigma}^{{\dagger}}=\sum_{p}R_{pa}a_{p\sigma}^{{\dagger}},$
$\displaystyle a_{a\sigma}=\sum_{p}R^{*}_{pa}a_{p\sigma}$ (52)

$\displaystyle
i\dot{\rho_{ab}}=\langle\Psi|\left[\sum_{\sigma}a^{{\dagger}}_{b\sigma}a_{a\sigma},H\right]|\Psi\rangle+$
(53) $\displaystyle\langle\Psi|\sum_{\sigma
p}\left(i\dot{R}_{pb}a^{{\dagger}}_{p\sigma}a_{a\sigma}+a^{{\dagger}}_{b\sigma}i\dot{R}_{pa}^{*}a_{p\sigma}\right)|\Psi\rangle$

This is further simplified by introducing the matrix $X$ using Eq. (34) and
solving for the commutator using the definitions of the 1- and 2-RDMs, Eqs.
(29) and (30).

$\displaystyle i\dot{\rho_{ab}}$ $\displaystyle=$
$\displaystyle\langle\Psi|\left[\sum_{\sigma}a^{{\dagger}}_{b\sigma}a_{a\sigma},H\right]|\Psi\rangle+\sum_{c}\left(\rho_{ac}X_{cb}-X_{ac}\rho_{cb}\right)$
(54) $\displaystyle=$
$\displaystyle\sum_{c}\left(h_{ac}\rho_{cb}-\rho_{ac}h_{cb}\right)+\sum_{cde}\left(V_{acde}\Gamma_{bcde}-V_{cbde}\Gamma_{cade}\right)+\sum_{c}\left(\rho_{ac}X_{cb}-X_{ac}\rho_{cb}\right)$
(55)

The first line of Eq. (55) can be written in terms of a generalized Fock
matrix in a completely analogous form to what is done both in TD-CASSCFSato
and Ishikawa (2013b) and static CASSCF theoryHelgaker  _et al._ (2012), with
the one caveat that the orbitals are complex-valued.

$\displaystyle i\dot{\rho}_{ab}$ $\displaystyle=$ $\displaystyle
F_{ba}-F_{ab}^{*}+\sum_{c}\left(\rho_{ac}X_{cb}-X_{ac}\rho_{cb}\right)$ (56)
$\displaystyle\equiv$ $\displaystyle
i\dot{\tilde{\rho}}_{ab}+\sum_{c}\left(\rho_{ac}X_{cb}-X_{ac}\rho_{cb}\right)$
(57)

This corresponds to Eq. (25) in the main text.

The generalized Fock matrix is given by

$\displaystyle F_{ba}$ $\displaystyle=$
$\displaystyle\sum_{c}h_{ac}\rho_{cb}+\sum_{cde}V_{acde}\Gamma_{bcde}.$ (58)

The expression for the generalized Fock matrix can be further simplified using
the specific form of the 1- and 2-RDMS for the CAS-like wavefunction for each
fragment. Specifically, $\rho_{cb}=0$ and $\Gamma_{bcde}=0$ if any index
corresponds to a virtual embedding orbital. Furthermore,
$\rho_{cb}=2\delta_{cb}$ if either $c$ or $b$ corresponds to a core embedding
orbital, $\Gamma_{bcde}=4\delta_{bc}\delta_{de}-\delta_{be}\delta_{cd}$ if all
indices correspond to a core orbital, $\Gamma_{bcde}=2\rho_{ed}\delta_{bc}$ if
only $b$ and $c$ correspond to a core orbital, and
$\Gamma_{bcde}=-\rho_{cd}\delta_{be}$ if only $b$ and $e$ correspond to a core
orbital. Using these simplifications and the symmetries of the 1- and 2-RDMs,
the generalized Fock matrix can be written as

$F_{ba}=\begin{cases}2^{I}F_{ab}+2^{A}F_{ab}&\mbox{if $b$ is a core
orbital}\\\ \sum\limits_{c\in
sub/bath}\,^{I}F_{ac}\rho_{cb}+\sum\limits_{cde\in
sub/bath}V_{acde}\Gamma_{bcde}&\mbox{if $b$ is a subsystem or bath orbital}\\\
0&\mbox{if $b$ is a virtual orbital},\\\ \end{cases}$ (59)

where the summations in the second line of Eq. (59) are over the subsystem and
bath orbitals. Finally, the inactive and active Fock matrices in Eq. (59) are
given byHelgaker  _et al._ (2012)

${}^{I}F_{ab}$ $\displaystyle=$ $\displaystyle h_{ab}+\sum_{c\in
core}\left(2V_{abcc}-V_{accb}\right)$ (60) ${}^{A}F_{ab}$ $\displaystyle=$
$\displaystyle\sum_{cd\in
sub/bath}\rho_{cd}\left(V_{abdc}-\frac{1}{2}V_{acdb}\right)$ (61)

###  V.4 Derivation of $\dot{\rho}^{glob}$

The explicit time dependence of $\rho^{glob}$ is necessary to define the
operator $\hat{G}$ that governs the mean-field dynamics, Eq. (46). The global
1-RDM is defined through the democratic partitioning of the correlated 1-RDMs
from the separate fragment calculations, such that

$\displaystyle{\rho}^{glob}_{pq}=\frac{1}{2}\sum_{ab}\left[R_{pa}^{f(p)}{(\rho}_{ab}^{emb})^{f(p)}R^{*f(p)}_{qb}\right.$
(62)
$\displaystyle\left.+R_{pa}^{f(q)}({\rho}^{emb}_{ab})^{f(q)}R^{*f(q)}_{qb}\right]$

To simplify the notation in the following steps, $\rho^{emb}$ will simply be
referred to as $\rho$ and the objects that belong to a fragment $f(q)$ will be
referred with a singular superscript $q$. Using this notation, the time
dependence of $\rho^{glob}$ involves the following six terms,

$\displaystyle i\dot{\rho}^{glob}_{pq}=\frac{1}{2}\sum_{ab}$
$\displaystyle\left[i\dot{R}_{pa}^{p}{\rho}_{ab}^{p}R^{*p}_{qb}+i{R}_{pa}^{p}{\rho}_{ab}^{p}\dot{R}^{*p}_{qb}\right.$
(63)
$\displaystyle+\left.i\dot{R}_{pa}^{q}{\rho}_{ab}^{q}R^{*q}_{qb}+i{R}_{pa}^{q}{\rho}_{ab}^{q}\dot{R}^{*q}_{qb}\right.$
$\displaystyle+\left.R_{pa}^{p}i\dot{\rho}_{ab}^{p}R^{*p}_{qb}+R_{pa}^{q}i\dot{\rho}_{ab}^{q}R^{*q}_{qb}\right].$

Inserting the definition of $\dot{R}$, Eq. (34), and $\dot{\rho}^{emb}$, Eq.
(57), yields

$\displaystyle i\dot{\rho}^{glob}_{pq}=\frac{1}{2}\sum_{abc}$
$\displaystyle\left[{R}_{pc}^{p}X_{ca}^{p}{\rho}_{ab}^{p}R^{*p}_{qb}-{R}_{pa}^{p}{\rho}_{ab}^{p}X_{bc}^{p}{R}^{*p}_{qc}+{R}_{pc}^{q}X_{ca}^{q}{\rho}_{ab}^{q}R^{*q}_{qb}-{R}_{pa}^{q}{\rho}_{ab}^{q}X_{bc}^{q}{R}^{*q}_{qc}\right]$
(64) $\displaystyle+\frac{1}{2}\sum_{ab}$
$\displaystyle\left[R_{pa}^{p}i\dot{\tilde{\rho}}_{ab}^{p}R^{*p}_{qb}+R_{pa}^{q}i\dot{\tilde{\rho}}_{ab}^{q}R^{*q}_{qb}\right]+\frac{1}{2}\sum_{abc}\left[{R}_{pa}^{p}{\rho}_{ac}^{p}X_{cb}^{p}{R}^{*p}_{qb}-{R}_{pa}^{p}X_{ac}^{p}{\rho}_{cb}^{p}R^{*p}_{qb}\right.$
$\displaystyle+\left.{R}_{pa}^{q}{\rho}_{ac}^{q}X_{cb}^{q}{R}^{*q}_{qb}-{R}_{pa}^{q}X_{ac}^{q}{\rho}_{cb}^{q}R^{*q}_{qb}\right]$

All terms that contain a dependence on the matrix $X$ cancel out above, such
that the time-dependence of the global density matrix only depends on the
generalized Fock matrices,

$\displaystyle i\dot{\rho}^{glob}_{pq}=$
$\displaystyle\frac{1}{2}\sum_{ab}\left[R_{pa}^{p}i\dot{\tilde{\rho}}_{ab}^{p}R^{*p}_{qb}+R_{pa}^{q}i\dot{\tilde{\rho}}_{ab}^{q}R^{*q}_{qb}\right]$
(65) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{ab}\left[R_{pa}^{p}i(F_{ba}-F_{ab}^{*})^{p}R^{*p}_{qb}\right.$
$\displaystyle+$
$\displaystyle\left.R_{pa}^{q}i(F_{ba}-F_{ab}^{*})^{q}R^{*q}_{qb}\right].$

This corresponds to Eq. (26) in the main text.

## Computational Details

###  V.5 System parameters

In this work we simulate a 1-dimensional (1D) single-impurity Anderson model
(SIAM) with and without a laser pulse, a set of 1D multi-impurity Anderson
models (MIAM), and a 2-dimensional SIAM. The Hamiltonian that encompasses both
the SIAM and MIAM in 1D and 2D, but without the laser pulse is given by

$H=-t\sum_{<ij>\sigma}(a_{i\sigma}^{\dagger}a_{j\sigma}+h.c.)+U\sum_{i\in
imp}n_{i\uparrow}n_{i\downarrow}$ (66)

which corresponds to Eq. (27) in the main text.

In every calculation throughout the manuscript, the initial state of the
system corresponds to the ground-state of the non-interacting system, where
$U=0$ and $t=1$. As described in the main text, non-equilibrium dynamics are
then initialized by running with a new value of the interaction term
$U^{\prime}\neq U$. The hopping term is fixed at $t=1$ between all sites in
all systems, even during the dynamics.

The 1D SIAM consists of $N=60$ sites where the impurity is located at
$i=\frac{N}{2}$. This system was used in Figure 2 in the main text.

We examine two types of MIAMs, both consisting of $N=120$ sites with
$N_{imp}=6$ impurities. The first type has the impurities evenly distributed
amongst the sites, located at sites $i=\\{41,49,57,65,73,81\\}$. This
corresponds to Figure 3(a) in the main text. The second type has the
impurities unevenly distributed amongst the sites, located at sites
$i=\\{45,51,56,67,74,80\\}$. This corresponds to Figures 3(b) and (c) in the
main text.

We also simulate the same 1D SIAM consisting of $N=60$ sites, but with the
addition of a laser pulse that is restricted to interact with only the central
20 sites of the system. The influence of the external electric field is
included into the Hamiltonian via the Peierls substitution in the hopping
terms between sites 20 to 40.

$a^{\dagger}_{i,\sigma}a_{i+1,\sigma}+h.c.\xrightarrow{}e^{iA(t)}a^{\dagger}_{i,\sigma}a_{i+1,\sigma}+h.c.,$
(67)

where $A(t)$ is the vector potential. In this work we choose $A(t)$ to have
the common form of

$A(t)=A_{0}e^{-(t-t_{0})^{2}/2t^{2}_{d}}cos[\omega_{0}(t-t_{0})],$ (68)

which involves an oscillating field of amplitude $A_{0}$ and central frequency
$\omega_{0}$ with a Gaussian envelope centered in time at $t_{0}$ of width
$t_{d}$. The parameters for the vector potential were chosen in agreement with
previous workShao  _et al._ (2020), such that $A_{o}=0.3$, $\omega=4.0$,
$t_{o}=5$, and $t_{d}=0.5$. In addition to the laser pulse, we still include a
sudden change from $U=0$ to $U^{\prime}\neq U$ at $t=0$, such that the
dynamics are initialized from an exact solution. This corresponds to Figure 5
in the main text.

Lastly, we also simulate a 2D SIAM without a laser pulse consisting of 8 sites
in both the $x$ and $y$ directions (total N=64 sites). The impurity is located
at ($x=3$, $y=4$). This corresponds to Figure 6 in the main text.

###  V.6 Numerical implementation of real-time pDMET

The initial state for every real-time pDMET calculation in the manuscript is
obtained from the ground-state solution of a static pDMET calculation of the
corresponding non-interacting system; the interaction term $U=0$ for all
sites. The static pDMET calculation is initialized using the 1-RDM obtained
from a restricted Hartree-Fock calculation. Exact diagonalization is used as
the solver for each fragment embedding problem. The static pDMET calculation
converges in a single-step as the restricted Hartree-Fock calculation is exact
for the non-interacting system and therefore provides the exact set of bath
states. The static pDMET calculation outputs a mean-field 1-RDM, $\rho^{mf}$,
(Eq. (3) in the main text), a set of embedding orbitals for each fragment,
$R^{emb}$, (Eq. (4) in the main text), a set of CI coefficients for each
fragment, $C$, (Eq. (5) in the main text), a global 1-RDM, $\rho^{global}$,
(Eq. (7) in the main text), and a set of natural orbitals, $U$, (Eq. (9) in
the main text). These are used as the inputs for the subsequent real-time
pDMET calculation.

The real-time pDMET calculation involves the simultaneous propagation of the
CAS-like wavefunction for each fragment and the mean-field 1-RDM. Practically
this implies that we directly propagate the expansion coefficients, the
embedding orbitals, and the mean-field 1-RDM as a function of time; the
equations of motion for these three quantities are given by Eqs. (19), (18),
and (22) in the main text, respectively. In addition, we also directly
propagate the global 1-RDM and natural orbitals as a function of time; the
equations of motion for these two quantities are given by Eqs. (26) and (21)
in the main text, respectively. To solve for these 5 coupled equations, it is
also necessary to solve, at each point in time, for the embedding Hamiltonian
along with the operators $\hat{X}$, $\hat{G}$, and the generalized Fock
matrix, $F$, given by Eqs. (20) and (23) in the main text and Eqs. (32)-(34)
in the SI, respectively. This procedure can be summarized as follows:
All equations referenced in the following algorithm are contained within the
SI unless specified otherwise.

  1. 1.

Given $C(t)$, $R(t)$, $\rho^{mf}(t)$, $\rho^{glob}(t)$, and $U(t)$

  2. 2.

Use $R(t)$ to generate the embedding Hamiltonian.Wouters  _et al._ (2016b)

  3. 3.

Use $C(t)$ and the embedding Hamiltonian to calculate $F(t)$ using Eqs.
(59)-(61)

  4. 4.

Use $F(t)$ and $R(t)$ to calculate $\dot{\rho}^{glob}(t)$ using Eq. (65)

  5. 5.

Use $U(t)$ and $\rho^{glob}(t)$ to obtain the eigenvalues of $\rho^{glob}(t)$
using Eq. (44)

  6. 6.

Use $U(t)$, $\dot{\rho}^{glob}(t)$, and the eigenvalues from step 5 to obtain
$G(t)$ using Eq. (46)

  7. 7.

Use $G(t)$ to obtain $\dot{U}(t)$ using Eq. (45)

  8. 8.

Use $\rho^{mf}(t)$ and $G(t)$ to obtain $\dot{\rho}^{mf}(t)$ using Eq. (49)

  9. 9.

Use $R(t)$ and $\rho^{mf}(t)$ to obtain the eigenvalues of $\rho^{mf}(t)$
using Eq. (33)

  10. 10.

Use $R(t)$, $\dot{\rho}^{mf}(t)$, and the eigenvalues from step 9 to obtain
$X(t)$ using Eq. (43)

  11. 11.

Use $X(t)$ to obtain $\dot{R}$(t) using Eq. (34)

  12. 12.

Use $X(t)$ and the embedding Hamiltonian to obtain $\dot{C}$(t) using Eq. (19)
from the main text

In principle, any possible numerical integrator could be used to solve for the
equations of motion once the desired time-derivatives are calculated following
the procedure above. In the context of this work, the 5 coupled equations of
motions are numerically solved using the fourth order Runge-Kutta method. This
implies that the above procedure is performed at each step within the Runge-
Kutta routine.

While theoretically it is only necessary to propagate $\rho^{mf}$ and $C$,
propagation of the embedding and natural orbitals is implemented to enhance
the numerical stability of the method in a highly degenerate space. All
calculations are initialized from the ground state where Coulomb interactions
are absent (U=0). The solution of the static pDMET is exact for this state
which eliminates any potential discrepancy between methods due to the error in
initial state calculation; this choice allows for a direct comparison of the
accuracy of solely the dynamics. However, this initialization procedure means
that many of the eigenvectors of $\rho^{mf}$ and $\rho^{global}$ are
degenerate. While both embedding and natural orbitals can be generated
directly by diagonalizing previously mentioned matrices, it becomes
numerically challenging to obtain smoothly transitioning eigenvectors at each
time step when eigenvectors are nearly degenerate. It was observed that direct
propagation of $\rho^{global}$ helps to ensure that the condition stated in
Eq. (9) in the main text is maintained throughout the calculation.

The time-dependence of nearly-degenerate embedding and natural orbitals
introduces additional numerically challenges associated with the definition of
corresponding operators $\hat{X}$ and $\hat{G}$ (Eq. (43), (46) in SI ). Both
matrices become undefined when orbitals are strictly degenerate. This gauge
invariance is fixed by setting the corresponding elements of the $X$ or $G$
matrices to 0. However, direct treatment of near-degenerate eigenvectors
results in numerically unstable terms and hence rapid divergence between
previously maintained matching conditions. To address this, regularization
parameters $\varepsilon_{X}$ and $\varepsilon_{G}$ are introduced for each
matrix, such that $X^{A}_{ba}=0$ and $G_{\mu\nu}=0$ if
$|\lambda_{a}-\lambda_{b}|<\varepsilon_{X}$ and
$|\lambda_{\mu}-\lambda_{\nu}|<\varepsilon_{G}$, respectively. Both parameters
define a near-degenerate region around strictly equivalent eigenvalues where
the elements are still set to 0. Similar parameters have been utilized in the
context of TD-CASSCF and MCTDH.Kretchmer and Chan (2018a, b); Manthe  _et al._
(1992); Meyer  _et al._ (1990) It was observed that only a small value of
$\varepsilon_{X}$ was necessary to obtain numerically stable dynamics. By
contrast, the numerical stability was more sensitive to the choice of
$\varepsilon_{G}$ and a larger value was necessary. This, however, may be an
artifact due to the high-degeneracy of the mean-field description of the
system at $t=0$ observed in the systems studied in this work; all calculations
are initialized from non-interacting ground-states, such that the natural
orbitals all start off to be either completely occupied or unoccupied.

###  V.7 Simulation parameters

The static pDMET calculations were initialized with a mean-field 1-RDM
obtained from a Restricted Hartree Fock calculation. Real-time pDMET was
performed using a time step of dt=0.0001 for the majority of calculations.
However, the calculations with strong correlation effects ($U^{\prime}=3,4$)
and a large fragment size ($N_{frag}=5$) required an adaptive time step, where
dt was reduced to dt=0.00001 when any two eigenvalues of the natural orbitals
are within $2*10^{-5}$. The stabilizing parameter for G matrix was set to
$\varepsilon_{G}=10^{-5}$ and the analogous parameter for X matrix was
$\varepsilon_{X}=10^{-9}$.

TD-DMRG calculations were initialized and propagated in time using the ITensor
package Fishman  _et al._ (2020). Initial static DMRG calculations for the
SIAM of N=60 sites and the MIAM of N=120 sites were done using 6 sweeps and a
truncation error of $10^{-10}$. The bond dimension was progressively increased
during the 6 sweeps using the following progression: 10, 20, 100, 100, 200,
500. The TD-DMRG calculation was performed using direct time-evolution of the
matrix product state using an exponentiated matrix product operator.

All calculations were treated with the “Fit” method within ITensor, where the
exponential propagator is approximated to be equal $e^{-itH}$ up to terms on
the order of $t^{2}$. The corresponding parameters are “Cutoff”= $10^{-7}$,
“MaxDim” = 500, and dt=0.001. Convergence with respect to the cutoff, time
step, and bond dimension was observed.

Time-dependent Hartree Fock calculations are initialized from a Restricted
Hartree Fock calculation to form the initial 1-RDM. The 1-RDM is propagated
using a time-dependent Fock operator using the fourth-order Runge-Kutta
method. All calculations were performed with a time step of dt=0.001.

##  VI References

## References

  * Peng  _et al._ (2019) P. Peng, C. Marceau, and D. M. Villeneuve, Nature Reviews Physics 1, 144 (2019).
  * Nisoli  _et al._ (2017) M. Nisoli, P. Decleva, F. Calegari, A. Palacios, and F. Martín, Chemical reviews 117, 10760 (2017).
  * Varillas  _et al._ (2022) R. B. Varillas, M. Lucchini, and M. Nisoli, Reports on Progress in Physics (2022).
  * Gallmann  _et al._ (2012) L. Gallmann, C. Cirelli, and U. Keller, Annual review of physical chemistry 63, 447 (2012).
  * Månsson  _et al._ (2021) E. P. Månsson, S. Latini, F. Covito, V. Wanie, M. Galli, E. Perfetto, G. Stefanucci, H. Hübener, U. De Giovannini, M. C. Castrovilli, _et al._ , Communications Chemistry 4, 1 (2021).
  * de la Cruz Valbuena  _et al._ (2019) G. de la Cruz Valbuena, F. VA Camargo, R. Borrego-Varillas, F. Perozeni, C. D’Andrea, M. Ballottari, and G. Cerullo, The journal of physical chemistry letters 10, 2500 (2019).
  * Krueger  _et al._ (2020) T. D. Krueger, S. A. Boulanger, L. Zhu, L. Tang, and C. Fang, Structural Dynamics 7, 024901 (2020).
  * Fuemmeler  _et al._ (2016) E. G. Fuemmeler, S. N. Sanders, A. B. Pun, E. Kumarasamy, T. Zeng, K. Miyata, M. L. Steigerwald, X.-Y. Zhu, M. Y. Sfeir, L. M. Campos, _et al._ , ACS central science 2, 316 (2016).
  * Ponseca  _et al._ (2017) C. S. Ponseca, P. Chábera, J. Uhlig, P. Persson, and V. Sundström, Chemical Reviews 117, 10940 (2017).
  * Cui  _et al._ (2019) L. Cui, S. Hur, Z. A. Akbar, J. C. Klöckner, W. Jeong, F. Pauly, S.-Y. Jang, P. Reddy, and E. Meyhofer, Nature 572, 628 (2019), number: 7771 Publisher: Nature Publishing Group.
  * Joachim  _et al._ (2000) C. Joachim, J. K. Gimzewski, and A. Aviram, Nature 408, 541 (2000), number: 6812 Publisher: Nature Publishing Group.
  * Nitzan and Ratner (2003) A. Nitzan and M. A. Ratner, Science 300, 1384 (2003).
  * Zhang  _et al._ (2021) L. Zhang, J. Zhou, H. Li, L. Shen, and Y. P. Feng, Applied Physics Reviews 8, 021308 (2021).
  * Dayen  _et al._ (2020) J.-F. Dayen, S. J. Ray, O. Karis, I. J. Vera-Marun, and M. V. Kamalakar, Applied Physics Reviews 7, 011303 (2020).
  * Singh  _et al._ (2019) R. Singh, W. Liu, J. Lim, I. Robel, and V. I. Klimov, Nature Nanotechnology 14, 1035 (2019), number: 11 Publisher: Nature Publishing Group.
  * Bande (2013) A. Bande, The Journal of Chemical Physics 138, 214104 (2013).
  * Jahnke  _et al._ (2020) T. Jahnke, U. Hergenhahn, B. Winter, R. Dörner, U. Frühling, P. V. Demekhin, K. Gokhberg, L. S. Cederbaum, A. Ehresmann, A. Knie, and A. Dreuw, Chemical Reviews 120, 11295 (2020).
  * Demekhin  _et al._ (2009) P. V. Demekhin, Y.-C. Chiang, S. D. Stoychev, P. Kolorenč, S. Scheit, A. I. Kuleff, F. Tarantelli, and L. S. Cederbaum, The Journal of Chemical Physics 131, 104303 (2009).
  * Ouchi  _et al._ (2011) T. Ouchi, K. Sakai, H. Fukuzawa, I. Higuchi, P. V. Demekhin, Y.-C. Chiang, S. D. Stoychev, A. I. Kuleff, T. Mazza, M. Schöffler, K. Nagaya, M. Yao, Y. Tamenori, N. Saito, and K. Ueda, Physical Review A 83, 053415 (2011).
  * Aoki  _et al._ (2014) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Reviews of Modern Physics 86, 779 (2014).
  * Härtle  _et al._ (2008) R. Härtle, C. Benesch, and M. Thoss, Physical Review B 77, 205314 (2008).
  * Rammer (2018) J. Rammer,  _Quantum transport theory_ (CRC Press, 2018).
  * Freericks  _et al._ (2006) J. Freericks, V. Turkowski, and V. Zlatić, Physical review letters 97, 266408 (2006).
  * Li  _et al._ (2013) B. Li, T. J. Levy, D. W. Swenson, E. Rabani, and W. H. Miller, The Journal of Chemical Physics 138, 104110 (2013).
  * Swenson  _et al._ (2011) D. W. Swenson, T. Levy, G. Cohen, E. Rabani, and W. H. Miller, The Journal of chemical physics 134, 164103 (2011).
  * Swenson  _et al._ (2012) D. W. Swenson, G. Cohen, and E. Rabani, Molecular Physics 110, 743 (2012).
  * Cohen  _et al._ (2014a) G. Cohen, D. R. Reichman, A. J. Millis, and E. Gull, Physical Review B 89, 115139 (2014a).
  * Cohen  _et al._ (2014b) G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Physical review letters 112, 146802 (2014b).
  * Gull  _et al._ (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Reviews of Modern Physics 83, 349 (2011).
  * Schiró (2010) M. Schiró, Physical Review B 81, 085126 (2010).
  * Werner  _et al._ (2009) P. Werner, T. Oka, and A. J. Millis, Physical Review B 79, 035320 (2009).
  * Cohen  _et al._ (2015) G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Physical review letters 115, 266802 (2015).
  * Mühlbacher and Rabani (2008) L. Mühlbacher and E. Rabani, Physical review letters 100, 176403 (2008).
  * Weiss  _et al._ (2009) S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Physical Review B 79, 249901 (2009).
  * Segal  _et al._ (2010) D. Segal, A. J. Millis, and D. R. Reichman, Physical Review B 82, 205323 (2010).
  * Li  _et al._ (2020a) X. Li, N. Govind, C. Isborn, A. E. DePrince, and K. Lopata, Chemical Reviews 120, 9951 (2020a).
  * Goings  _et al._ (2018) J. J. Goings, P. J. Lestrange, and X. Li, WIREs Computational Molecular Science 8, e1341 (2018), _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1341.
  * Schlünzen  _et al._ (2019) N. Schlünzen, S. Hermanns, M. Scharnke, and M. Bonitz, Journal of Physics: Condensed Matter 32, 103001 (2019).
  * Church and Rubenstein (2021) M. S. Church and B. M. Rubenstein, The Journal of Chemical Physics 154, 184103 (2021).
  * Thiel (2014) W. Thiel, WIREs Computational Molecular Science 4, 145 (2014), _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1161.
  * Chatterjee and Makri (2019) S. Chatterjee and N. Makri, The Journal of Physical Chemistry B 123, 10470 (2019).
  * Helgaker  _et al._ (2012) T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen, and K. Ruud, Chemical Reviews 112, 543 (2012).
  * Norman and Dreuw (2018) P. Norman and A. Dreuw, Chemical Reviews 118, 7208 (2018).
  * V. Toukach and P. Ananikov (2013) F. V. Toukach and V. P. Ananikov, Chemical Society Reviews 42, 8376 (2013).
  * Baiardi and Reiher (2019) A. Baiardi and M. Reiher, Journal of chemical theory and computation 15, 3481 (2019).
  * Frahm and Pfannkuche (2019) L.-H. Frahm and D. Pfannkuche, Journal of Chemical Theory and Computation 15, 2154 (2019).
  * Baiardi (2021) A. Baiardi, Journal of Chemical Theory and Computation 17, 3320 (2021).
  * Li  _et al._ (2020b) W. Li, J. Ren, and Z. Shuai, The Journal of Physical Chemistry Letters 11, 4930 (2020b).
  * Sonk  _et al._ (2011) J. A. Sonk, M. Caricato, and H. B. Schlegel, The Journal of Physical Chemistry A 115, 4678 (2011).
  * Durden and Levine (2022) A. S. Durden and B. G. Levine, Journal of Chemical Theory and Computation 18, 795 (2022).
  * Peng  _et al._ (2018) W.-T. Peng, B. S. Fales, and B. G. Levine, Journal of chemical theory and computation 14, 4129 (2018).
  * Sato and Ishikawa (2013a) T. Sato and K. L. Ishikawa, Phys. Rev. A 88, 023402 (2013a).
  * Miranda  _et al._ (2011) R. P. Miranda, A. J. Fisher, L. Stella, and A. P. Horsfield, J. Chem. Phys. 134, 244101 (2011).
  * Kochetov and Bokarev (2022) V. Kochetov and S. I. Bokarev, Journal of Chemical Theory and Computation 18, 46 (2022).
  * Isborn  _et al._ (2007) C. M. Isborn, X. Li, and J. C. Tully, The Journal of Chemical Physics 126, 134307 (2007).
  * Knizia and Chan (2012) G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
  * Knizia and Chan (2013) G. Knizia and G. K.-L. Chan, J. Chem. Theory and Comput. 9, 1428 (2013).
  * Wouters  _et al._ (2016a) S. Wouters, C. A. Jiménez-Hoyos, Q. Sun, and G. K.-L. Chan, J. Chem. Theory and Comput. 12, 2706 (2016a).
  * Wu  _et al._ (2019) X. Wu, Z.-H. Cui, Y. Tong, M. Lindsey, G. K.-L. Chan, and L. Lin, The Journal of Chemical Physics 151, 064108 (2019).
  * Zheng and Chan (2016) B.-X. Zheng and G. K.-L. Chan, Physical Review B 93, 035126 (2016).
  * Cui  _et al._ (2020) Z.-H. Cui, T. Zhu, and G. K.-L. Chan, Journal of Chemical Theory and Computation 16, 119 (2020).
  * Faulstich  _et al._ (2022) F. M. Faulstich, R. Kim, Z.-H. Cui, Z. Wen, G. Kin-Lic Chan, and L. Lin, Journal of Chemical Theory and Computation 18, 851 (2022).
  * Reinhard  _et al._ (2019) T. E. Reinhard, U. Mordovina, C. Hubig, J. S. Kretchmer, U. Schollwöck, H. Appel, M. A. Sentef, and A. Rubio, Journal of Chemical Theory and Computation 15, 2221 (2019).
  * Pham  _et al._ (2020) H. Q. Pham, M. R. Hermes, and L. Gagliardi, Journal of Chemical Theory and Computation 16, 130 (2020).
  * Fan and Jie (2015) Z. Fan and Q.-l. Jie, Physical Review B 91, 195118 (2015).
  * Gunst  _et al._ (2017) K. Gunst, S. Wouters, S. De Baerdemacker, and D. Van Neck, Physical Review B 95, 195127 (2017).
  * Tran  _et al._ (2019) H. K. Tran, T. V. Voorhis, and A. J. W. Thom, J. Chem. Phys. 151, 034112 (2019).
  * Mitra  _et al._ (2021) A. Mitra, H. Q. Pham, R. Pandharkar, M. R. Hermes, and L. Gagliardi, The Journal of Physical Chemistry Letters , 11688 (2021).
  * Qiao and Jie (2021) J. Qiao and Q. Jie, Computer Physics Communications 261, 107712 (2021).
  * Pham  _et al._ (2018) H. Q. Pham, V. Bernales, and L. Gagliardi, Journal of Chemical Theory and Computation 14, 1960 (2018).
  * Sandhoefer and Chan (2016) B. Sandhoefer and G. K.-L. Chan, Physical Review B 94, 085115 (2016).
  * Ricke  _et al._ (2017) N. Ricke, M. Welborn, H.-Z. Ye, and T. Van Voorhis, Molecular Physics 115, 2242 (2017).
  * Kretchmer and Chan (2018a) J. S. Kretchmer and G. K.-L. Chan, The Journal of Chemical Physics 148, 054108 (2018a).
  * Wouters  _et al._ (2016b) S. Wouters, C. A. Jiménez-Hoyos, Q. Sun, and G. K.-L. Chan, Journal of Chemical Theory and Computation 12, 2706 (2016b).
  * Sato and Ishikawa (2013b) T. Sato and K. L. Ishikawa, Physical Review A 88, 023402 (2013b).
  * Kretchmer and Chan (2018b) J. S. Kretchmer and G. K.-L. Chan, J. Phys. Chem. Lett. 9, 2863 (2018b).
  * Manthe  _et al._ (1992) U. Manthe, H. Meyer, and L. S. Cederbaum, The Journal of Chemical Physics 97, 3199 (1992).
  * Meyer  _et al._ (1990) H. D. Meyer, U. Manthe, and L. S. Cederbaum, Chemical Physics Letters 165, 73 (1990).
  * Helgaker  _et al._ (2014) T. Helgaker, P. Jorgensen, and J. Olsen,  _Molecular electronic-structure theory_ (John Wiley & Sons, 2014).
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, _et al._ , The Journal of chemical physics 153, 024109 (2020).
  * (81) https://github.com/jskretchmer/real_time_pDMET.
  * Titvinidze  _et al._ (2015) I. Titvinidze, A. Schwabe, and M. Potthoff, The European Physical Journal B 88, 9 (2015).
  * Eickhoff and Anders (2020) F. Eickhoff and F. B. Anders, Physical Review B 102, 205132 (2020), arXiv:2008.09034 [cond-mat].
  * Baiardi and Reiher (2020) A. Baiardi and M. Reiher, The Journal of Chemical Physics 152, 040903 (2020).
  * Shao  _et al._ (2020) C. Shao, T. Tohyama, H.-G. Luo, and H. Lu, Physical Review B 101, 045128 (2020).
  * Fishman  _et al._ (2020) M. Fishman, S. R. White, and E. M. Stoudenmire, “The ITensor software library for tensor network calculations,” (2020), arXiv:2007.14822 .