# Constructing Tensor Network Influence Functionals for General Quantum
Dynamics

Erika Ye<EMAIL_ADDRESS>Division of Engineering and Applied Sciences,
California Institute of Technology, Pasadena, USA 91125  Garnet Kin-Lic Chan
<EMAIL_ADDRESS>Division of Chemistry and Chemical Engineering, California
Institute of Technology, Pasadena, USA 91125

###### Abstract

We describe an iterative formalism to compute influence functionals that
describe the general quantum dynamics of a subsystem beyond the assumption of
linear coupling to a quadratic bath. We use a space-time tensor network
representation of the influence functional and investigate its approximability
in terms of its bond dimension and time-like entanglement in the tensor
network description. We study two numerical models, the spin-boson model and a
model of interacting hard-core bosons in a 1D harmonic trap. We find that the
influence functional and the intermediates involved in its construction can be
efficiently approximated by low bond dimension tensor networks in certain
dynamical regimes, which allows the quantum dynamics to be accurately computed
for longer times than with direct time evolution methods. However, as one
iteratively integrates out the bath, the correlations in the influence
functional can first increase before decreasing, indicating that the final
compressibility of the influence functional is achieved via non-trivial
cancellation.

##  I Introduction

Obtaining the long-time dynamics of a large quantum system is in general
intractable due to the exponential scaling of Hilbert space with respect to
system size and the associated exponential growth of spatial entanglement with
time. Fortunately, in many cases, one is most interested in the dynamics of
observables defined on a small subset of the full system, and can thus reframe
the dynamics of the observable from the viewpoint of a subsystem coupled to a
bath Breuer and Petruccione (2002). If one is able to determine the dynamics
of the bath and its influence on the subsystem, then the original dynamics
problem is reduced to obtaining the dynamics of the subsystem Feynman (1948);
Nakajima (1958); Zwanzig (1960); Breuer  _et al._ (2001).

The influence functional (IF) Feynman (1948) method provides an exact
framework for computing the dynamics of an arbitrary bath and its interactions
with the subsystem. However, the cost of computing the IF without
approximation is comparable to determining the dynamics of the full system and
the size of the IF scales exponentially with the number of time steps. Thus it
is not usually possible to use the IF method without additional
approximations.

The IF can be viewed as reweighting the path integral of the subsystem. In the
case of linear coupling to a harmonic bath, Feynman and Vernon derived an
analytical form of the IF Feynman and Vernon (1963); Vernon (1959), which
takes the form of the Boltzmann weight of a complex valued Hamiltonian defined
in the time direction with pairwise interactions between time points.

For many physical bath spectral densities, it is natural to assume that the
pairwise time interaction is short-ranged in time, corresponding to a finite
”memory” in the influence of the bath, and many numerical approximations have
successfully taken advantage of this short-range temporal nature Makri and
Makarov (1995); Makri (1999); Strathearn  _et al._ (2017); Weiss  _et al._
(2008); Segal  _et al._ (2010); Cohen  _et al._ (2015); Chen  _et al._ (2017);
Cerillo and Cao (2014); Rosenbach  _et al._ (2016); Shi and Geva (2003); Kelly
_et al._ (2016); Cohen and Rabani (2011); Kidon  _et al._ (2018); Chatterjee
and Makri (2019). For IF methods, the assumption of limited memory allows one
to remove the exponential growth of cost of the quantum dynamics with
simulation time, thus making long time-scale quantum dynamical simulations
possible. For example, in QuAPI Makri and Makarov (1995), one approximates the
analytical harmonic bath IF by only including terms acting on time steps
within a finite time window. Alternatively one can construct an ansatz for the
IF; a natural choice is a matrix product state (MPS) in the time direction,
which compactly encodes short-range time correlations. This approach has been
used in a number of recent works, which, although they do not necessarily use
the language of IFs, all proceed by constructing a compressed version of the
IF Lerose  _et al._ (2021) or a closely related object such as the auxiliary
density operator as defined in QuAPI Strathearn  _et al._ (2018); Jorgensen
and Pollock (2019), the process tensor Pollock  _et al._ (2018); Cygorek  _et
al._ (2021), or other variants Bañuls  _et al._ (2009); Müller-Hermes  _et
al._ (2012); Luchnikov  _et al._ (2019, 2020).

In this work, we are interested in using the IF method for computing the
dynamics of a subsystem within a general quantum system. Such a subsystem may
arise as part of a larger interacting problem (in which case, the subsystem
might not be different from other parts of the system) or it may arise from a
system-bath model. In either case, the couplings and bath cannot be assumed to
be linear and quadratic respectively and thus the analytical form of the IF is
not known. Instead, the IF is simply a particular integral of the space-time
dynamics that must be obtained numerically. To do this concretely, we can use
a tensor network description of the space-time dynamics. For a 1-dimensional
representation of the system and bath, the tensor network is thus defined in
1+1 dimensions. Traditional time evolution corresponds to contracting this
network first along the time axis, and if done exactly leads to exponential
cost with time. Computing the IF corresponds to the contraction of the network
in the spatial direction to yield a final matrix product state IF defined
along the time direction, and if done exactly has an exponential cost with the
spatial length. A similar framework appears in the recently proposed modular
path integral (MPI) procedure Makri (2018a, 2021); Kundu and Makri (2019,
2021); Makri (2018b); Kundu and Makri (2020). However, to avoid exponential
computational cost, the contractions must be performed approximately. There
has been some discussion regarding using filtering methods to reduce the cost
of the MPI algorithm Makri (2018a). In contrast, we utilize tensor network
compression and propose a tensor network contraction procedure similar to that
presented by Banuls et al. Bañuls  _et al._ (2009); Müller-Hermes  _et al._
(2012) and Lerose et al. Lerose  _et al._ (2021), assuming a translationally
invariant infinite problem. However, we will describe a general tensor network
procedure to construct the IF without such assumptions and explore the
numerical feasibility of doing so to compute quantum dynamics in different
regimes beyond the commonly considered linear coupling to quadratic baths.

The paper is organized as follows. We first translate the IF into a space-time
tensor network language and describe an iterative algorithm to compute it. We
then investigate the compressibility of the IF and its ability to produce
long-time dynamics, first for the canonical spin-boson model where the
analytical IF is known, and then an interacting hard-core boson model where
there is no analytical expression, which corresponds to the case of general
quantum dynamics. We analyze the time-like entanglement both in the IF itself
as well as the intermediates that arise as the bath is numerically integrated
out. We end with a brief discussion of the implications of this work for
future studies.

##  II Theory

###  II.1 Definition of the Influence Functional

To introduce notation, we first recall the definition of the influence
functional (IF) Feynman (1948). The influence functional describes how the
path integral of a subsystem is reweighted, under the influence of dynamical
coupling to a bath. To obtain an explicit form, we define a full system
composed of the subsystem of interest and the coupled bath. At time $T$, we
denote the subsystem density matrix by $\rho_{s}(s_{T})$, where $s_{T}$ is a
basis for the density matrix, and the bath density matrix is analogously
written as $\rho_{b}(b_{T})$. The basis of the full system is spanned by the
product space $\\{s\\}\otimes\\{b\\}$. The evolution of the density matrix is
given by a linear operator, the Liouville operator $L$, which we partition as
$L=L_{s}+L_{bs}$ where $L_{s}$ contains the component operating only on the
subsystem and $L_{bs}$ contains the component on the bath and interactions
between the subsystem and bath. If we further assume the system dynamics obeys
Hamiltonian evolution, then the Liouville action can be written as
$L\,\cdot=[H,\cdot]$.

Formally, $\rho_{s}(s_{T})$ is obtained by time evolving the entire system and
tracing out the bath degrees of freedom. The path integral expression,
assuming a second-order Trotter decomposition of the time evolution operator
into $N$ timesteps of length $\epsilon$, is

$\displaystyle\rho_{s}(s_{T})$ $\displaystyle=$
$\displaystyle\operatorname{Tr}_{b_{T}}\biggl{[}$
$\displaystyle\sum_{s_{t_{N-1}}}\sum_{b_{t_{N-1}}}\dots\sum_{s_{t_{0}}}\sum_{b_{t_{0}}}$
$\displaystyle\mathinner{\langle\langle{s_{T}}|}e^{-\frac{i}{2}L_{s}\epsilon}e^{-iL_{bs}\epsilon}e^{-iL_{s}\epsilon}\mathinner{|{s_{t_{N-1}},b_{t_{N-1}}}\rangle\rangle}\times$
$\displaystyle\mathinner{\langle\langle{s_{t_{N-1}},b_{t_{N-1}}}|}e^{-iL_{bs}\epsilon}e^{-iL_{s}\epsilon}\mathinner{|{s_{t_{N-2}},b_{t_{N-2}}}\rangle\rangle}\times\dots$
$\displaystyle\times$
$\displaystyle\mathinner{\langle\langle{s_{t_{1}},b_{t_{1}}}|}e^{-iL_{bs}\epsilon}e^{-\frac{i}{2}L_{s}\epsilon}\mathinner{|{s_{t_{0}},b_{t_{0}}}\rangle\rangle}\times$
$\displaystyle\langle\langle
s_{t_{0}},b_{t_{0}}|\rho(s_{t_{0}},b_{t_{0}})\rangle\rangle\,\biggr{]}$ (1)

where $\rho(s_{t_{0}},b_{t_{0}})$ is the initial state of the system, and the
double bra/ket notation indicates we are working in Liouville space, with the
density matrix being a vector in this space. For simplicity, we assume there
are no correlations between the subsystem and bath initially such that
$\mathinner{|{\rho(s_{t_{0}},b_{t_{0}})}\rangle\rangle}=\mathinner{|{\rho_{s}(s_{t_{0}})}\rangle\rangle}\mathinner{|{\rho_{b}(b_{t_{0}})}\rangle\rangle}$.
Furthermore, $L_{bs}$ is typically assumed to be diagonal in the basis
$\\{s\\}$ (we lift both these restrictions below). Then Eq. (1) becomes

$\displaystyle\rho_{s}(s_{T})$
$\displaystyle=\sum_{s_{t_{N-1}}}\,\dots\sum_{s_{t_{0}}}\,\,\mathinner{\langle\langle{s_{T}}|}e^{-\frac{i}{2}L_{s}\epsilon}\mathinner{|{s_{t_{N}}}\rangle\rangle}\times$
$\displaystyle\mathinner{\langle\langle{s_{t_{N}}}|}e^{-iL_{s}\epsilon}\mathinner{|{s_{t_{N-1}}}\rangle\rangle}\times\dots$
$\displaystyle\times\mathinner{\langle\langle{s_{t_{1}}}|}e^{-\frac{i}{2}L_{s}\epsilon}\mathinner{|{s_{t_{0}}}\rangle\rangle}\langle\langle
s_{t_{0}}|\rho_{s}(s_{t_{0}})\rangle\rangle$ $\displaystyle\times
I(s_{t_{1}},s_{t_{2}},...,s_{t_{N}})$ (2)

where $I(s_{t_{1}},s_{t_{2}},...,s_{t_{N}})$ is the influence functional

$I(\ldots)=\operatorname{Tr}_{b_{T}}\left[e^{-iL_{bs}(s_{t_{N}})\epsilon}\dots
e^{-iL_{bs}(s_{t_{1}})\epsilon}\mathinner{|{\rho_{b}(b_{t_{0}})}\rangle\rangle}\right]$
(3)

with $L_{bs}(s)=\langle\langle s|L_{bs}|s\rangle\rangle$. The IF assigns a
complex weight to each configuration of the system path integral.
Consequently, the storage of the IF grows exponentially with number of time
steps $N$.

###  II.2 Generalized Influence Functional in Tensor Network Language

####  II.2.1 Influence Functional Structure

Translating Eq. (1) and Eq. (2) into the tensor network language is
straightforward, and is shown in diagrams Fig. 1(a) and (b), respectively. The
time evolution operators of the system and system-bath dynamics respectively
appear as boxes with two and four legs in Fig. 1(a), with the legs labelled by
the bra and ket basis states. The elements in the tensors are
$\langle\langle{s_{t}}|e^{-iL_{s}\epsilon}|s_{t^{\prime}}\rangle\rangle$ and
$\langle\langle{s_{t},b_{t}}|e^{-iL_{bs}\epsilon}|s_{t^{\prime}},b_{t^{\prime}}\rangle\rangle$,
respectively. In the case that $L_{bs}$ is diagonal with respect to the
subsystem basis states $\\{s\\}$, then the elements of the system-bath time
evolution operator are
$\langle\langle{s_{t},b_{t}}|e^{-iL_{bs}\epsilon}|s_{t},b_{t^{\prime}}\rangle\rangle$
and can be depicted as a box with three legs labelled by $s_{t},b_{t}$, and
$b_{t^{\prime}}$, as in Fig. 1(b). The final trace operation over the bath
degrees of freedom at the last time step can be written as a vector that is
contracted with the corresponding leg of the time evolution operator tensor.

The influence functional element $I(s_{t_{1}},s_{t_{2}},\ldots,s_{t_{N}})$ is
the object within the blue rectangle of Fig. 1(b), obtained after performing
tensor contractions over all bath tensors as denoted by the connected lines.
Within this diagrammatic picture it is easy to depict the generalization of
the influence functional to a correlated initial state. In this case, the
dotted line indicates a correlated initial state with entanglement between the
subsystem and bath,

$\displaystyle\rho(s_{t_{0}},b_{t_{0}})=\sum_{\alpha}\rho_{s\alpha}(s_{t_{0}})\rho_{b\alpha}(b_{t_{0}})$
(4)

and the influence functional is defined with an additional index,
$I(s_{t_{1}},s_{t_{2}},\ldots,s_{t_{N}};\alpha)$. Similarly, if $L_{bs}$
cannot be diagonalized in the subsystem basis $\\{s\\}$ then we can generalize
the influence functional to contain two subsystem indices at each intermediate
time,
$I(s_{t_{1}}s_{t^{\prime}_{1}},s_{t_{2}}s_{t^{\prime}_{2}},\ldots,s_{t_{N}}s_{t^{\prime}_{N}})$,
and it is the object within the blue rectangle of Fig. 1(a). Given the
influence functional, arbitrary time-correlation functions can be computed as
shown in Fig. 1(c).

Because the influence functional has a one-dimensional structure along the
time axis, it is natural to rewrite it as a matrix product state of $N$
tensors (see Fig. 1(d)), i.e.

$\displaystyle I(s_{t_{1}},s_{t_{2}},\ldots,s_{t_{N}})=$ (5)
$\displaystyle\sum_{\\{i\\}}A^{(1)}_{i_{1}}(s_{t_{1}})$ $\displaystyle
A^{(2)}_{i_{1},i_{2}}(s_{t_{2}})\ldots A^{(N)}_{i_{N-1}}(s_{t_{N}})$

where $A(s_{t})$ denotes a matrix of dimension $D\times D$ for each element of
the basis $s_{t}$, except for $A^{(1)}(s_{t_{1}})$ and $A^{(N)}(s_{t_{N}})$
which are $D$ dimensional row and column vectors respectively. In MPS
language, $D$ is referred to as the virtual bond dimension. Although any IF
can be represented as a MPS for sufficiently large $D$, the MPS of small bond
dimension naturally capture sums of exponentially decaying time-correlations
along the time axis. The key system-specific questions to understand are thus
(i) is the IF itself representable by an MPS of low bond dimension, in
physically relevant dynamical and interaction regimes, and (ii) can the IF be
constructed with manageable cost in those regimes. It is important to note
that an affirmative answer to (i) does not imply an affirmative answer to
(ii).

Figure 1: (a) Time evolution of
$\rho(s_{t_{0}})=\sum_{\alpha}\rho_{s,\alpha}(s_{t_{0}})\otimes\rho_{b,\alpha}(s_{t_{0}})$
in Liouville space with second order Trotter decomposition between system and
interaction dynamics. (b) The same as (a) but assuming $L_{bs}$ is diagonal
with respect to subsystem basis. The boxed regions are the generalized and
traditional definitions of the influence functional, respectively. (c)
Measurement of the time-correlation
$\langle\hat{O}_{2}(t_{1})\,\hat{O}_{1}(t_{0})\rangle_{\rho}$. (d) Matrix
product state representation of influence functional. The labels
$\\{s_{t_{m}}\\}$ and $\\{b_{t_{m}}\\}$ index the system and bath states at
time step $m$, respectively. The labels $\\{i_{m}\\}$ index the virtual bonds.
Lines that connect two tensors (blocks) represent tensor contraction over the
labeled indices.

####  II.2.2 Space-Time Tensor Network Representation

To define an approximate procedure to construct the influence functional for
complex bath dynamics, we first write down a space-time representation of the
full system dynamics. We first assume that the bath Hilbert space is a product
space over $K$ modes,

$\displaystyle\\{b\\}=\\{b^{1}\\}\otimes\\{b^{2}\\}\otimes\ldots\otimes\\{b^{K}\\}$
(6)

We can then formally express the system density matrix at any time as a matrix
product state (MPS)

$\displaystyle\rho(s_{t},b_{t})=\sum_{\\{i\\}}C^{(0)}_{i_{0}}({s_{t}})\,C^{(1)}_{i_{0},i_{1}}(b^{1})\ldots
C^{(K)}_{i_{K-1}}(b^{K})$ (7)

with a bond dimension denoted $D_{\rho}$.

Figure 2: (a) Space-time tensor network representing time evolution of a
system represented as a 1D MPS (row of circles) with second order Trotter
decomposition between system (single square on left most column) and
interaction dynamics (row of squares). Pictured are two time steps applied to
a system coupled to five bath sites. (b) Tensor network representation of
influence functional containing $N$ time steps. The top-most row represents
the time evolution operator with the trace operation applied.

Similarly, the Liouville evolution operator can be written as a matrix product
operator (MPO) with the elements

$\displaystyle\langle\langle$ $\displaystyle
s_{t},b_{t}|e^{-iL_{bs}\epsilon}|s_{t^{\prime}},b_{t^{\prime}}\rangle\rangle=$
(8)
$\displaystyle\sum_{\\{i\\}}M^{(0)}_{i_{0}}({s_{t},s_{t^{\prime}}})\,M^{(1)}_{i_{0},i_{1}}(b^{1}_{t},b^{1}_{t^{\prime}})\ldots
M^{(K)}_{i_{K-1}}(b^{K}_{t},b^{K}_{t^{\prime}})$

and has a bond dimension $D_{L}$. Note that since the Liouville operator is
assumed time-independent, $D_{L}$ is fixed.  At last the time step, the bath
degrees of freedom are traced out so the bath sites are now in MPS form,

$\displaystyle\langle\langle s_{t_{N}}$
$\displaystyle|\,\text{Tr}_{b_{T}}\left[e^{-iL_{bs}\epsilon}\right]\,|\,s_{t_{N}^{\prime}},b_{t_{N-1}}\rangle\rangle=$
(9)
$\displaystyle\sum_{\\{i\\}}M^{(0)}_{i_{0}}(s_{t_{N}},s^{\prime}_{t_{N}})\tilde{M}^{(1)}_{i_{0},i_{1}}(b^{1}_{t_{N-1}})\ldots\tilde{M}^{(K)}_{i_{K-1}}(b^{K}_{t_{N-1}})$

where $\tilde{M}=\text{Tr}_{b_{T}}(M)$ with $M$ as defined in Eq. (8).

If the Hamiltonian consists of only nearest neighbor interactions, one can
obtain the operator using a Trotter-Suzuki decomposition of nearest neighbour
gates and then directly map them onto a matrix product operator. Otherwise,
for more general interactions one can use a 4th-order Runge-Kutta expansion
García-Ripoll (2006). In this case, the matrix product operator can have large
bond dimension, but can be compressed by allowing for truncation errors of
$~{}\mathcal{O}(\epsilon^{5})$. In the cases studied here, the subsystem is
small enough such that $e^{-iL_{s}\epsilon}$ can be obtained exactly.

The full time evolution of the system with $K$ bath modes and $N$ time steps
thus corresponds to the two-dimensional tensor network diagram shown in Fig.
2(a). Correspondingly, the space-time representation of the influence
functional is shown in Fig. 2(b).

####  II.2.3 Transverse Contraction Scheme

The most common way to contract a 2D space-time tensor network is from bottom
to top, i.e. in the direction of increasing time Vidal (2004); Daley  _et al._
(2004); Schollwöck (2011); White and Feiguin (2004); White  _et al._ (2018);
Schollwöck and White (2006); Ronca  _et al._ (2017); Haegeman  _et al._
(2016). We refer to this as direct time evolution. For example, contracting
the network in Fig. 2(a) row by row yields the system density matrix at each
time step as an MPS. The cost of contracting two rows is
$\mathcal{O}(KD_{L}^{2}D_{\rho}^{2}d_{\rho}^{2})$ where $d_{\rho}$ is the
dimension of the density matrices represented at each MPS site, and $D_{\rho}$
is the bond dimension of the density matrix MPS. For exact time evolution
$D_{\rho}$ grows by a factor of $D_{L}$ at each time step. This means that in
the worst case, $D_{\rho}$ grows exponentially with time.

However, if the time-correlations in the influence functional decay with long
time, then this implies that the influence functional ultimately can be
represented by a matrix product with low bond dimension along the time axis.
This suggests that a more efficient contraction strategy is to contract column
wise (in the transverse direction to time). The cost of contracting columns
together in the process of construct the final influence functional is
$\mathcal{O}(ND_{L}^{2}D_{I}^{2}d_{\rho}^{2})$ where $D_{I}$ is the bond
dimension of the column bath MPS (defined along the time direction). For exact
contraction, $D_{I}$ will grow by $d_{\rho}$ at each contraction step and thus
in the worst case scales exponentially with $K$. But, it need not have a
dependence on the total simulation time.

In practice, due to the exponential growth of the bond dimension, exact
contraction of the 2D tensor network (in either direction) is often too
expensive. In Banuls et al Bañuls  _et al._ (2009); Müller-Hermes  _et al._
(2012), the explicit transverse contraction of the 2D tensor network was
avoided by assuming that the system is infinite and translationally invariant,
in which case the result of the infinite contraction of columns is
proportional to the maximal eigenvector of the column transfer operator.

Figure 3: Iteration of transverse contraction scheme to compute IF tensor
network. The time evolution operators $e^{-iL_{bs}\epsilon}$ at each time step
are the rows of the grid and are each represented as an MPO, and the subsystem
of interest is at the left-most site. Before contraction, we first
canonicalize each row into left canonical form as indicated by the right
pointing triangles along the rows. The rightmost two columns are then
contracted and compressed to fixed bond dimension $D_{I}$ using the standard
MPS compression algorithm, where the column is first converted into a
canonical form (here, top canonical form) and then compressed by singular
value decomposition in the reverse direction (leaving it in bottom canonical
form). The canonical form implies that the tensors satisfy an isometric
condition (see diagram on the right); e.g. the right pointing arrow implies
contraction of a tensor with its complex conjugate over the left, up, and down
indices yields the identity matrix. The procedure is repeated until all
columns have been contracted.

Alternatively, one may use standard matrix product state techniques to
compress the intermediates that arise during the contraction to restrict bond
dimensions $D_{\rho}$ or $D_{I}$ to some constant value Vidal (2004); Daley
_et al._ (2004); Schollwöck (2011). In this case, the cost of the algorithm is
dominated by the cost of the MPS compression (which requires performing a
series of singular value decompositions), which scales like
$\mathcal{O}(KD_{\rho}^{3}D_{L}^{3}d_{\rho})$ for the direct time evolution
case (upward contraction by row), and
$\mathcal{O}(ND_{I}^{3}d_{\rho}^{3}D_{L})$ for the IF time evolution procedure
(sideways contraction by column).  In this paper, we use such an approximate
transverse contraction scheme (compressing to a fixed bond dimension of the
contraction intermediates) to compute the IF (Fig. 2(b)) for systems with
arbitrary baths. The algorithm is to iteratively contract and compress the
columns from the edges of the bath inwards to the sites connected to the
subsystem (Fig. 3). Assuming the subsystem is defined as the leftmost site, we
start from the rightmost boundary column, and then the column is absorbed
leftward to make a new boundary column, which is compressed using standard MPS
compression to a pre-specified maximum bond dimension $D_{I}$. This is the
standard ”boundary contraction” algorithm of 2D tensor networks Murg  _et al._
(2007). Key to the success of the algorithm and the quality of the compression
is the choice in gauge (a redundant degree of freedom in all tensor networks).
We choose the gauge as shown in Fig. 3. In the case where the bath extends
both to the left and right of the site of interest, we compute the IFs
corresponding to the left and right bath sites separately (see Appendix for
details). This contraction scheme was implemented within the Quimb tensor
network library Gray (2018).

For purposes of comparison, we will also present reference dynamics generated
by standard MPS time evolution (ie. contracting the space-time tensor network
in the usual time direction) Vidal (2004); Zwolak and Vidal (2004); García-
Ripoll (2006); Schollwöck (2011). Because the underlying full system dynamics
is governed by Hamiltonian evolution in the problems that we study, we have
the option to apply $e^{-iH\epsilon}$ as a commutator to the square root of
the density matrix
($[e^{-iH\epsilon}\rho^{1/2}][\rho^{1/2{\dagger}}e^{iH\epsilon}]$) or via
$e^{-iL\epsilon}$ directly. We refer to the former as Hilbert time evolution
(HTE) and the latter as Liouville time evolution (LTE). In HTE, the compressed
tensor network dynamics is carried out for the pseudowavefunction
$\psi=e^{-iH\epsilon}\rho^{1/2}$ Verstraete  _et al._ (2004); Barthel  _et
al._ (2009); Karrasch  _et al._ (2012). HTE has the advantage that the
compressed density matrix is always positive definite, although correlations
between the bra and ket sides of the density matrix are less compressible. In
existing literature, this is sometimes referred to as purification-based time
evolution. In the case that $\rho$ is a pure state, this method is equivalent
to traditional MPS Hilbert space time evolution.

##  III Results

###  III.1 Spin Boson Model

First, we consider the well-studied spin-boson model, in which a single spin
is linearly coupled to a bath of non-interacting harmonic oscillators,

$\displaystyle H_{SB}=\Delta
S_{X}+\int\\!d\omega\left[S_{Z}(g(\omega)a_{\omega}+g^{*}(\omega)a^{\dagger}_{\omega})+\omega
a^{\dagger}_{\omega}a_{\omega}\right]$ (10)

where $\Delta$ is the tunneling strength between the two subsystem states, and
the system-bath coupling strength $g(\omega)$ is determined from the bath
spectral density function $J(\omega)$ by

$\displaystyle|g(\omega)|^{2}=\frac{1}{\pi}J(\omega)$ (11)

In the case of an Ohmic bath with exponential cut-off,

$J(\omega)=\frac{\pi}{2}\alpha\omega e^{-\omega/\omega_{c}}$ (12)

where $\alpha$ is the Kondo parameter and $\omega_{c}$ is the cut-off
frequency. Typically one computes the dynamics from a factorized initial state
$\mathinner{|{\rho_{s}(s_{t_{0}})}\rangle\rangle}\mathinner{|{\rho_{b}(\beta)}\rangle\rangle}$,
where $\mathinner{|{\rho_{b}(\beta)}\rangle\rangle}$ is the Gibbs thermal
state of the isolated bath at finite temperature $\beta$. In this paper, we
set $\Delta=1.0$, $\omega_{c}=7.5$, and $\beta=5.0$.

The spin-boson model exhibits a dynamical phase transition from thermalizing
to localizing behavior at $\alpha=1.0+\mathcal{O}(\Delta/\omega_{c})$
Strathearn  _et al._ (2018), and is often cited as an example of physically
relevant non-Markovian dynamics Breuer and Petruccione (2002); Vega and Alonso
(2017). Because of the linear coupling and harmonic bath, the IF may be
computed via an analytical expression. There already exist several methods of
obtaining accurate dynamics for various bath coupling strengths and spectral
densities  Makri and Makarov (1995); Makri (2014, 2020a); Wang and Thoss
(2008); Strathearn  _et al._ (2018) . We thus use this model as a benchmark to
understand the properties of the influence functional, its compressibility,
and the accuracy of the tensor network contraction approximation.

Figure 4:  (a, d) Dynamics obtained using the analytical IF with
$m_{\max}=\infty$ and compressed to bond dimension $D_{I}=256$ for various
coupling strengths $\alpha$. (b, e) R.m.s. error for the analytical IF with
respect to $D_{I}$. (c, f) R.m.s error for capped IFs with respect to
$m_{\max}$ computed using $D_{I}=256$. The errors are obtained using dynamics
from the IF with no cap ($m_{\max}=\infty$) and $D_{I}=256$ as reference.
These calculations are performed using $N=100$ time steps. The system
parameters are $\Delta=1.0$, bath inverse temperature $\beta=5.0$, and an
Ohmic bath spectral density with exponential cut-off $\omega_{c}=7.5$ for
(a-c) and $\omega_{c}=2.0$ for (d-f).

####  III.1.1 Compressibility of Analytical IF

We first investigate the compressibility of the analytical expression for the
influence functional for the spin-boson model. Denoting the system basis
$\mathinner{|{s}\rangle\rangle}\equiv\mathinner{|{s^{+}}\rangle}\mathinner{\langle{s^{-}}|}$
where
$\mathinner{|{s^{\pm}}\rangle}\in\\{\mathinner{|{1}\rangle},\mathinner{|{-1}\rangle}\\}$
and $\mathinner{|{1}\rangle},\mathinner{|{-1}\rangle}$ are the eigenstates of
the $S_{Z}$ operator, then the influence functional can be written as

$\displaystyle I_{SB}$
$\displaystyle=\exp\left\\{-\sum_{k=1}^{N}\sum_{k^{\prime}=1}^{k}(s_{t_{k}}^{+}-s_{t_{k}}^{-})(\eta_{kk^{\prime}}s_{t_{k^{\prime}}}^{+}-\eta_{kk^{\prime}}^{*}s_{t_{k^{\prime}}}^{-})\right\\}$

This explicitly shows the form of the influence functional as the Boltzmann
weight of a complex spin Hamiltonian with the spins interacting along the time
axis via the long-range pairwise ”interaction” $\eta_{kk^{\prime}}$. We can
further factorize the weights into contributions for times $(t_{1})$,
$(t_{1},t_{2}),\ldots,(t_{1},t_{N})$, giving

$\displaystyle I_{SB}=$
$\displaystyle\prod_{k=1}^{N}I_{0}(s_{t_{k}}^{\pm})\prod_{k=1}^{N-1}I_{1}(s_{t_{k}}^{\pm},s_{t_{k+1}}^{\pm})\,\dots$
$\displaystyle\prod_{k=1}^{N-m}I_{m}(s_{t_{k}}^{\pm},s_{t_{k+m}}^{\pm})\>\dots\prod_{k=1}^{1}I_{N}(s_{t_{k}}^{\pm},s_{t_{k+N-1}}^{\pm})$
(14)

where

$I_{m}=\exp\left\\{-(s_{t_{k+m}}^{+}-s_{t_{k+m}}^{-})(\eta_{k+m,k}s_{t_{k}}^{+}-\eta_{k+m,k}s_{t_{k}}^{-})\right\\}$
(15)

The $\eta_{kk^{\prime}}$ interaction terms can be derived from the spectral
density of the bosonic bath. Expressions for $\eta_{kk^{\prime}}$ are given in
Eq. (12) of Ref. Makri and Makarov, 1995, where they use $\Delta t$ instead of
$\epsilon$ to denote the timestep.

In QuAPI, one considers interactions $\eta_{kk^{\prime}}$ for times $k$ and
$k^{\prime}$ within $m_{\text{max}}$ of each other, where $m_{\text{max}}$ is
treated as a convergence parameter Makri and Makarov (1995). Then, one can
evaluate the influence functional (or its effect on the dynamics Makri (2020a,
b)) with a computational cost exponential in $m_{\text{max}}$. Since QuAPI
often converges rapidly with $m_{\text{max}}$, one also expects the matrix
product state representation of the IF (Eq. (5)) to be compressible to small
bond dimension. One way to verify this would be to construct the large
influence functional object as an exact tensor, and then compress it into a
matrix product state. Because of the exponential storage of the tensor with
time, this is possible only for a small number of time points $N$.
Alternatively, one could build the influence functional iteratively (i.e.
piece by piece in Eq. (14)) and compress at each step. This is the idea behind
TEMPO and related methods Strathearn  _et al._ (2018); Jorgensen and Pollock
(2019) which exploit the compressibility of the augmented density matrix, the
influence functional applied to the subsystem density matrix, i.e.

$\displaystyle
A(s_{t_{1}},\ldots,s_{t_{m_{\max}}})=\sum_{s_{t_{0}}}I(s_{t_{0}},\ldots,s_{t_{m_{\max}}})\rho_{s}(s_{t_{0}})\,\,.$
(16)

Note that because $I$ is composed of commuting pieces there are many possible
decompositions and thus sequences of iterative constructions.

To verify the compressibility of the IF itself, we approximate $I_{SB}$ with
no memory cutoff ($m_{\text{max}}=\infty)$ as an MPS of bond dimension $D_{I}$
using an iterative scheme (see Appendix) and determine the error in the
resulting on-site dynamics, using the $D_{I}=256$ result as reference. Here
and throughout the paper, the error is computed as the r.m.s. deviation of the
dynamics of an observable with respect to some reference over the time
interval of the plot. The $D_{I}=256$ results are used as the reference
because most of the $D_{I}=128$ dynamics calculations are already converged to
within an r.m.s. error of 0.001 with respect to it.

Fig. 4 compares convergence with respect to $D_{I}$ for infinite
$m_{\text{max}}$, as well as convergence with respect to $m_{\text{max}}$ for
fixed $D_{I}=256$. As expected, the influence functional generally becomes
less compressible with increased Kondo parameter $\alpha$. However, while
convergence of the IF with respect to $m_{\text{max}}$ and $D_{I}$ seem
comparable for baths with larger cut-off frequencies ($\omega_{c}=7.5$), the
IF converges more quickly with respect to $D_{I}$ for smaller $\omega_{c}$.
The key advantage of using a compressed matrix representation, as opposed to
truncating $\eta_{k,k+m}$ at some $m_{\max}$ as in QUAPI, is that this does
not eliminate the effects of long-range memory Strathearn  _et al._ (2018).
Overall, this result confirms that for certain spectral densities, the
influence functional can be efficiently written as a low-rank MPS, consistent
with the findings from TEMPOStrathearn  _et al._ (2018); Jorgensen and Pollock
(2019).

####  III.1.2 Finite Size Harmonic Bath

We now use the spin-boson model to examine if the IF can be constructed
efficiently from our tensor network contraction scheme. To do so, we consider
a finite size harmonic bath with $K$ sites, where the bosons are capped to
some finite number of states. For small baths and small boson cap, the
transverse contraction can be performed without compression, allowing a
numerical test of the compression procedure. The analytical IF, assuming
bosons with infinite boson cap, can also be computed. The bath is
characterized by the discretized spectral density,

$J_{D}(\omega)=\sum_{j}\frac{J(\omega)}{\rho(\omega)}\delta(\omega-\omega_{j})$

We use a linear discretization of the bath sites, such that the bath density
is $\rho(\omega)=K/\omega_{m}$, where $\omega_{m}$ is the maximum boson
frequency used. Here, we set $\omega_{m}=10$.

Figure 5: Comparisons of average error in
$\mathinner{\langle{S_{Z}(t)}\rangle}$ with respect to reference direct time
evolution results for the IF where the IF is (left) constructed using exact
time evolution, contracted exactly, and finally compressed to bond dimension
$D_{I}$ at the end, and where the IF is constructed with (middle) exact time
evolution and (right) RK4 time evolution but iteratively contracted and
compressed using the transverse compression scheme. In these calculations, we
use $\Delta=1.0$, bath inverse temperature $\beta=5.0$, cut-off frequency
$\omega_{c}=7.5$, and assume a discrete Ohmic bath with 2 modes at
$\omega=[5.0,10.]$. The bosonic bath sites are approximated to have only a
maximum boson number of 2. Time evolution is performed using $N=100$ time
steps with a time step of $\epsilon=0.05$. The plots show that for small bath
sizes, the error in the iterative compression scheme is dominated by the lack
of compressibility of the final IF.

First, we consider a system with only 2 bath modes (in the number basis) each
with a maximum boson number of 2, for $N=100$ time steps of size
$\epsilon=0.05$. For this small system, we use the exact time evolution
operator of $L_{bs}$ and compute the IF by exact transverse contraction,
applying compression only to the final IF object. Fig. 5(a) shows the error of
the exact IF compressed to bond dimension $D_{I}$. As expected, the IF is much
less compressible than with the continuous bath density in the last section,
due to the small bath size. The error decreases only slightly until it drops
suddenly once the bond dimension is large enough to capture the IF exactly,
and further, the compression error increases with $\alpha$. We then perform
the same analysis but for IFs computed by transverse contraction with
compression, for both exact and RK4 time evolution. As seen in Fig. 5, the
errors are comparable to those obtained when compressing the final exact IF.
This indicates that for this small problem, there is little additional error
added by the iterative contraction, and that the time-step error is
negligible: the error is dominated by the compressibility of the final IF
itself (which is low when the bath size is small).

Figure 6: (top) Expectation values $\mathinner{\langle{S_{Z}(t)}\rangle}$
obtained from the analytical IF (thinner colored lines) and IF from transverse
contraction (thicker colored lines) for bath size $K=11$ with $\alpha=$ 0.5,
1.0, and 1.5. The thick dashed line corresponds to reference dynamics from
direct time evolution of the density matrix. The transverse contraction scheme
introduces additional error with respect to the analytical result, which
increases with $\alpha$. (bottom) Time-averaged error in IF dynamics with
respect to $D_{I}=128$ results, obtained using (left) the transverse
contraction scheme and (right) the analytical IF for the finite bath of size
$K$, respectively. The $K=\infty$ bath size corresponds to a continuous bath.
Figure 7: Expectation values $\mathinner{\langle{N_{i}(t)}\rangle}$ for the
hard-core boson model with $U=0$ for system sizes (top) $L=9$ and (bottom)
$L=43$ for sites $i=\\{0,L/2,L/4\\}$ obtained using the iterative contraction
scheme with bond dimension $D_{I}$ as labeled. The thick dashed line is the
dynamics obtained by direct Hilbert-space time evolution, compressed to bond
dimension $D_{\rho}=256$, shown only up to times for which the results are
converged. The thin dashed line, shown only in the left-most plots but with
similar behavior for all, is the dynamics obtained by direct time evolution of
the full density matrix. Unphysical behavior suggests loss of positivity of
the density matrix. For these calculations, the initial state is a pure
product state with alternating spins, $\mathinner{|{0,1,0,1,...}\rangle}$.
Time evolution is performed using Trotter steps with a $N=100$ time steps of
$\Delta t=0.1$. Consistent with earlier observations, larger $D_{I}$ is needed
to accurately capture the IF for smaller bath sizes. However, as shown in the
$L=43,i=21$ simulation, the IF of comparable bond dimension can simulate
dynamics for longer times than direct HTE.

Next we investigate systems with larger bath sizes in Fig. 6(a). We first
examine the time-dynamics of the analytical IF (i.e. without any boson cap,
and without transverse contraction) for a discretized spectral density with 11
bath sites, as well as the IF computed by transverse contraction, using a
boson cap of 2. Because the system size is so small, the exact reference
dynamics for a boson cap of 2 can be generated by direct MPS time evolution
(here we use $D_{\rho}=64$ and Hamiltonian time evolution). From this
comparison we observe two things. First, compared to using the continuous bath
density, the error of the analytical IF dynamics is increased, although it is
still somewhat compressible. For the same $D_{I}$, the errors using transverse
contraction are larger, suggesting that at intermediate points in the
transverse contraction, there is more time-like entanglement than in the final
IF itself. In Fig. 6(b) we show the time-averaged error of the IF dynamics as
a function of the number of bath sites. We see that this error decreases as
the number of bath sites increases, both for the analytical IF and the
transverse contraction. This is consistent with the idea that smoother bath
densities are more ”compressible”.

###  III.2 1D Hard-core Boson Model

We next study dynamics of a 1D hard-core boson (HCB) lattice model, in which
each lattice site is either unoccupied ($\mathinner{|{0}\rangle}$) or occupied
($\mathinner{|{1}\rangle}$) by a single bosonic particle. We consider the
Hamiltonian

$H_{HCB}=\sum_{j}\left[-J(a_{j}^{\dagger}a_{j+1}+h.c.)+Un_{j}n_{j+1}+\frac{K}{2}n_{j}j^{2}\right]$
(17)

where $a_{j}^{\dagger}=\mathinner{|{1}\rangle}\mathinner{\langle{0}|}$,
$a_{j}=\mathinner{|{0}\rangle}\mathinner{\langle{1}|}$ are hard-core boson
creation and annihilation operators at the $j^{\text{th}}$ lattice site, and
$n_{j}=a^{\dagger}_{j}a_{j}$ is the number operator. This Hamiltonian is
intended to mimic the Bose-Hubbard model dynamics often simulated by cold-atom
experiments, which was shown to be difficult to compute using direct time
evolution methods Trotzky  _et al._ (2012), where the on-site interaction term
is replaced by a nearest-neighbor interaction term. We also include a harmonic
potential term to emulate a cold atom trap. We assume a pure initial state
$\mathinner{|{0,1,0,...}\rangle}$ such that there is one particle at every
other lattice site, and set the parameters $J=1$ and $K=10^{-2}$ while varying
$U$. For non-zero interaction term $U$, there is no analytical form for the
IF.

Figure 8: (a,b,c) Expectation values of the site-occupancy
$\mathinner{\langle{N_{i}(t)}\rangle}$ for the hard-core boson model of length
$L=43$ at lattice site $i=21$ for different coupling strengths $U=$ 0.5, 1.0,
and 2.0. Lines become dotted after divergence of ¿0.03 with respect to the
$D_{I}=256$ IF results. (d) The r.m.s. errors with respect to the $D_{I}=256$
IF results as a function of $U$. Simulation parameters are otherwise the same
as in Fig. 7. Compared to the $U=0$ case, a larger bond dimension is needed,
particularly at around $U=3.0$ where the r.m.s. error peaks. In contrast, the
HTE dynamics converge more quickly with increasing $U$.

We compute the dynamics of $\mathinner{\langle{N_{i}(t)}\rangle}$ at lattice
sites $i=\\{0,L/4,L/2\\}$, where $L$ is the length of the 1-D chain. For
longer chains, the rapid growth of entanglement means that direct MPS time
evolution (either using Hamiltonian evolution, denoted HTE or Liouvillian
evolution, denoted LTE) with a finite $D_{\rho}$ can only obtain converged
dynamics up to a finite time. We consider two chain lengths: $L=9$ where
converged HTE MPS dynamics can be used as a reference, and $L=43$, where the
HTE MPS dynamics appear to be not fully converged (no longer within 0.03 of
$D_{\rho}=512$ results) for the full simulated time. To obtain the dynamics
using the influence functional method, we partition the lattice such that site
$i$ is the subsystem of interest and the remaining sites are the bath.

The $U=0$ dynamics for $L=9$ and $L=43$ is shown in Fig. 7. Direct MPS LTE
shows unphysical behaviour for large system sizes, presumably because of the
loss of positivity at some point in the dynamics. In contrast, the dynamics
obtained using the iteratively contracted IF are more stable, highlighting the
innate compressibilty of time evolution tensor network along the time axis as
opposed to the spatial axis. For $L=9$, the IF dynamics only appears to begin
to converge by $D_{I}=128$ with respect to the (exact) HTE dynamics, having
less than 0.03 r.m.s error in $\mathinner{\langle{N_{i}(t)}\rangle}$ over the
simulated time interval and deviations within 0.08. Thus, there appears to be
no significant advantage to using the the IF method over direct HTE for small
system sizes. Conversely, for the $L=43$ system, the IF dynamics are converged
with respect to $D_{I}=256$ results by $D_{I}=64$, with less than 0.02 r.m.s
error and a maximum deviation of 0.04 for $i=10$ and less than 0.001 r.m.s
error and a maximum deviation of 0.004 for $i=21$. Note that the $i=10$
dynamics converge more slowly because effectively the site is coupled to two
separate baths, one of which is small. However, the IF method still
outperforms direct HTE.

For $U>0$, as shown in Fig. 8, the $D_{I}=64$ results appear less converged
than the non-interacting case, but the $D_{I}=128$ results are converged
(r.m.s. errors of the $D_{I}=128$ observable dynamics with respect to the
$D_{I}=256$ results are less than 0.03), for times longer than that accessible
by direct time evolution. This shows that the IF-based dynamics can produce
the correct oscillatory behaviour of the density as a function of time, which
is not captured by the direct MPS time evolution despite using a larger bond
dimension (this difficulty with the long time oscillatory tail has previously
been noted in other cold atom simulations Trotzky  _et al._ (2012)). However,
while the IF method notably outperforms direct HTE at small $U$, the two
methods become comparable at larger $U\approx 1.5$ once the oscillatory tail
is sufficiently dampened.

###  III.3 Entanglement Spectrum

Figure 9: (Top left) Entanglement entropy and (top right) spectrum of
normalized singular values at the middle of the ”bath” boundary column MPS
after each contraction and compression iteration for the SB model
($\alpha=1.0$). As bath size increases, the EE of the IF decreases, converging
to some finite value. (Bottom left) Entanglement entropy of the hard-core
boson model with $L=43$ for the ”bath” boundary MPS with contractions starting
from the right edge, plotted for different values of $U$. The decrease in EE
with respect to iteration shows that EE decreases with system size. (Bottom
right) Normalized singular values for the hardcore-boson model with $L=43$ and
$U=0.5$. Surprisingly, for insufficient $D_{I}$, the singular values take on
large and incorrect values, yielding an artificially large EE.

The accuracy of the transverse contraction scheme depends on the entanglement
in the time-like direction. Recall that our contraction algorithm starts with
the farthest column (an MPS), and at each iteration another column is
contracted into this boundary. Thus, as the iteration number increases, the
boundary column represents more of the bath. For both the spin-boson model and
HCB model, we measure the singular values at the middle of the boundary ”bath”
MPS during each step of the iterative contraction scheme. The entanglement
entropy (EE) and spectrum of the singular values (normalized so that
$\sum_{i}|s_{i}|^{2}=1$) are plotted in Fig. 9. For the spin-boson model, only
results for $\alpha=1.0$ are shown; the only notable difference for other
$\alpha$ is that the EE increases with $\alpha$. Consistent with the
observations in our simulations above, the EE decreases as one increases bath
size. For the SB model, the EE decreases with increased number of bath sites
in the discretization until convergence. For the HCB model, only if
sufficiently large enough bond dimension is used does the EE decrease with
increasing iteration number. Otherwise, the EE stays at a large value
throughout the contraction scheme and the gap between the dominant and non-
dominant singular values decreases; this makes the EE of the smaller $D_{I}$
approximation larger than that of the larger $D_{I}$ approximation. Overall,
this suggests that the final compressibility of the IF emerges from the
cancellation of many different correlations as one iteratively contracts out
the bath.

##  IV Conclusions

In this work, similar to some other recent contributions Lerose  _et al._
(2021); Cygorek  _et al._ (2021), we have used the representation of the
influence functional within the tensor network language, motivated by the
limitations of modeling spatial entanglement growth in quantum dynamics. We
have discussed a transverse tensor network contraction algorithm that allows
us to compute the influence functional in cases where the analytical form is
not known. We have applied this algorithm to study both the canonical spin-
boson model as well as an interacting hard-core boson chain where the bath is
not quadratic (i.e. interacting). We find that the compressibility of the
influence functional is controlled by several factors, principally the size of
the bath, as well as the nature of the interactions. In addition, although the
time-like correlations may ultimately be short-ranged in the final influence
functional, during the transverse tensor network contraction to construct it,
it is possible to proceed through intermediate quantities with larger time-
like entanglement. This suggests a complicated picture where time-like
correlations first accumulate as the bath is integrated out before finally
cancelling in the influence functional itself. In the regimes where the
influence functional and all intermediate quantities are compressible, as in
some interaction regimes in the interacting hardcore boson model we have
studied, it is possible to outperform conventional tensor network time
evolution methods at longer times.

There are many possible directions for further investigation. For example,
there are natural extensions to higher-dimensional interacting problems and
fermionic systems, as well as more complicated correlation functions. Also, a
better theoretical understanding of how correlations grow and cancel out in
the transverse contraction scheme may lead to a deeper understanding of the
generation of memory in quantum dynamics, the master equation formalism Shi
and Geva (2003), improved contraction schemes, and ultimately new algorithms
to carry out longer time dynamical simulations.

##  V Acknowledgments

GKC was supported by the Center for Molecular Magnetic Quantum Materials, an
Energy Frontier Research Center funded by the U.S. Department of Energy,
Office of Science, Basic Energy Sciences under Award No. DE-SC0019330. EY was
primarily supported by the Google PhD fellowship program with supplemental
salary support from the Center for Molecular Magnetic Quantum Materials.
Support for industrial mentorship for EY was provided via the QISE-NET
program, funded via NSF award DMR-1747426.

##  VI Data and Code Availability

Data and code is available upon request.

##  Appendix A Transverse Contraction Procedure

Figure 10: Tensor network diagrams depicting the contraction and compression
of an MPO-MPS product. (a) The MPO and MPS columns to be contracted together.
(b) The resulting MPS after exactly contracting the network along the
horizontal bonds indexed by $\\{i\\}$. We use a double line to represent the
merging of the two original vertical lines, as in Eq. (23). (c) The new MPS
put into right canonical form. (d) The MPS is compressed to a smaller bond
dimension starting from the left-most tensor and moving to the right, and the
final MPS is now in left canonical form. The double line is reduced to a
single line to signify the size reduction of the large vertical bonds down to
the desired bond dimension.

We describe a single step of the transverse contraction procedure, which in
essence is performing an MPO-MPS contraction and then compressing the
resulting MPS to the desired bond dimension $D$, described in detail in Ref.
Schollwöck, 2011. We assume that all tensors have been prepared in the
appropriate gauge (discussed later) as shown in Fig. 3. This algorithm is also
the standard boundary contraction algorithm for 2D tensor networks, as
described in Ref. Murg  _et al._ , 2007.

For simplicity, consider the influence functional tensor network for $K$ bath
sites, as discussed in the main text. We start with the farthest boundary
column MPS (the column for the $K^{\text{th}}$ bath site) which is given by

$\displaystyle P_{K}$
$\displaystyle(i_{0},...,i_{N})=\sum_{\\{a\\}}A^{(0)}_{a_{0}}({i_{0}})A^{(1)}_{a_{0},a_{1}}(i_{1})\ldots
A^{(N)}_{a_{N-1}}(i_{N})$ (18)

with

$\displaystyle A^{(0)}_{a}({i})=C^{(K)}_{i}(a)$ (19a) $\displaystyle
A^{(n)}_{a,a^{\prime}}(i)=M^{(K)}_{i}(a^{\prime},a)\qquad\text{for}\,\,n=1,...,N-1$
(19b) $\displaystyle A^{(N)}_{a}(i)=\tilde{M}^{(K)}_{i}(a)$ (19c)

where $C^{(K)}$ and $M^{(K)}$ are the right-most tensors of the density matrix
MPS (Eq. (7)) and Liouville evolution operator (Eq. (8)), respectively, and
$\tilde{M}^{(K)}$ denotes the Liouville evolution operator with traced out
bath degrees of freedom.

The next column can be interpreted as an MPO, given by

$\displaystyle P_{K-1}(j_{0},...,j_{N};i_{0},...,i_{N})$ $\displaystyle\hskip
1.42271pt=\sum_{\\{b\\}}B^{(0)}_{b_{0}}({j_{0},i_{0}})\,B^{(1)}_{b_{0},b_{1}}(j_{1},i_{1})\ldots
B^{(N)}_{b_{N-1}}(j_{N},i_{N})\,\,.$ (20)

with

$\displaystyle B^{(0)}_{b}(j,i)=C^{(K-1)}_{j,i}(b)$ (21a) $\displaystyle
B^{(n)}_{b,b^{\prime}}(j,i)=M^{(K-1)}_{j,i}(b^{\prime},b)\qquad\text{for}\,\,n=1,...,N-1$
(21b) $\displaystyle B^{(N)}_{b}(j,i)=\tilde{M}^{(K-1)}_{j,i}(b)$ (21c)

Contracting the two columns yields a new MPS,

$\displaystyle\sum_{\\{i_{K-1}\\}}P_{K-1}(j_{0},...,j_{N};i_{0},...,i_{N})\,\,\times\,\,P_{K}(i_{0},...,i_{N})$
$\displaystyle\hskip
28.45274pt=\sum_{\\{c\\}}F^{(0)}_{c_{0}}({i_{0}})F^{(1)}_{c_{0},c_{1}}(i_{1})\ldots
F^{(N)}_{c_{N-1}}(i_{N})$ (22)

where the sum over $c$ indicates sums over all $(a,b)$, and

$\displaystyle
F^{(n)}_{c,c^{\prime}}(i)=\sum_{i}B^{(n)}_{b,b^{\prime}}(j,i)\,A^{(n)}_{a,a^{\prime}}(i)$
(23)

for all $n$ indexing the tensors in the MPO and MPS. Thus, if the MPS was of
bond dimension $D_{K}$ and the MPO was of bond dimension $D_{K-1}$, the new
MPS has a bond dimension of $D_{K}D_{K-1}$. We then must compress this new MPS
back to one of a lower bond dimension $D$ (again, see Ref. Schollwöck, 2011
for a more in depth discussion). To do so, we first put the MPS in right or
left canonical form.

In left canonical form, all tensors in the MPS (except for the rightmost one)
satisfies

$\displaystyle\sum_{l,u}L_{l,r}(u)\,\,\,$ $\displaystyle
L^{*}_{l,r^{\prime}}(u)=\delta_{r,r^{\prime}}$ (24)

where $\delta$ is the Kronecker delta. The tensor network diagram depicting
Eq. (24) is

In our tensor network diagrams, we denote tensors in canonical form using
triangles, pointing in the direction of the uncontracted leg.

Similarly, in right canonical form, all tensors (except for the leftmost one)
satisfies

$\displaystyle\sum_{r,u}R_{l,r}(u)\,\,\,$ $\displaystyle
R^{*}_{l^{\prime},r}(u)=\delta_{l,l^{\prime}}$ (25)

We are able to define canonical forms because tensor networks have a gauge
degree of freedom. This means that the choice of tensors in the network is not
unique. One can see this by introducing a set of matrices $X$, $X^{-1}$, which
clearly satisfy $XX^{-1}=\mathbb{I}$, along any line connecting two tensors.

Canonicalization can be performed using singular value decomposition (SVD).
Suppose that we are interested in written the MPS in left canonical form, and
that all tensors left of the $n^{\text{th}}$ tensor are already in left
canonical form. We then take the SVD of the $n^{\text{th}}$ tensor,

$\displaystyle
F^{(n)}_{c,c^{\prime}}(i)=\sum_{\sigma}U_{c,\sigma}(i)\,\Sigma_{\sigma}V^{\dagger}_{\sigma,c^{\prime}}$
(26)

where $\Sigma$ are the singular values from the decomposition of the tensor.
Note that by definition, $U$ is left canonical, as desired, and thus will be
used as the new $n^{\text{th}}$ tensor. The remaining matrices are then pushed
into the $(n+1)^{\text{th}}$ tensor,

$\displaystyle F^{(n)}_{c,\sigma}(i)\leftarrow U_{c,\sigma}(i)$ (27)
$\displaystyle
F^{(n+1)}_{\sigma,c^{\prime\prime}}(i)\leftarrow\sum_{c^{\prime}}\Sigma_{\sigma}V^{\dagger}_{\sigma,c^{\prime}}F^{(n+1)}_{c^{\prime},c^{\prime\prime}}(i)\,\,.$
(28)

By iteratively performing this operation starting from the left-most tensor
all the way to the right end of the MPS, the MPS is put into left canonical
form. The procedure for expressing the MPS in right canonical form is
analogous.

MPS compression is performed in the same way, except only the largest $D$
singular values are retained, generating some error. For minimal compression
errors, the MPS must be in left (right) canonical form prior to performing the
iterative compression procedure starting from the right (left) end.

In Fig. 3 we depict the contraction of the columns of the (1+1)D influence
functional tensor network. The rows are initialized in left canonical form.
After the two right-most columns are contracted, the product is canonicalized
and the compressed using the procedure discussed above. Because of the
vertical orientation, the left and right canonical forms are depicted by
triangles pointing upwards and downwards along the column.

##  Appendix B Matrix Product Form of Analytical IF

As discussed in the text, the expression for the analytical IF in discretized
time steps is given by Eq. (LABEL:eq:exactSB) and we wish to write it in the
MPS form with physical bonds that index the states of the density matrix at
each timestep,

$\displaystyle I(s_{t_{1}},$ $\displaystyle s_{t_{2}},\ldots,s_{t_{N}})=$
$\displaystyle\sum_{i_{1},...,i_{N-1}}A(s_{t_{1}})_{i_{1}}A(s_{t_{2}})_{i_{1},i_{2}}\ldots
A(s_{t_{N}})_{i_{N-1}}$ (29)

One possible way to construct the IF is to take the product of factors in Eq.
(LABEL:eq:exactSB) in the order $I=I_{0}I_{1}\ldots I_{N-1}$. We start by
using the $I_{0}$ terms which are in the form of a product state (MPS with
bond dimension $1$). We then multiply by each of the subsequent $I_{m}$ and
compress into an MPS after each $I_{m}$ is applied. Multiplying by $I_{m}$ can
be viewed as multiplication by an MPO where the tensors are very sparse. The
two-body terms $I_{m}$ for $m>1$ are long-range operators, and must be padded
with identities to skip over the times in the middle. More explicitly, the
MPOs are

$\displaystyle
I_{m}(s_{t_{k}},s_{t_{k+m}})\,=\sum_{v}\left[I_{m}^{(0)}(s_{t_{k}})\right]_{v}\left[I_{m}^{(1)}(s_{t_{k+m}})\right]_{v}$
(30) $\displaystyle\hskip 14.22636pt\rightarrow\hskip
12.80365pt\leavevmode\hbox to157.46pt{\vbox
to64.18pt{\pgfpicture\makeatletter\hbox{\hskip 78.7305pt\lower-32.09074pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}} {}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{-68.28731pt}{0.0pt}\pgfsys@lineto{-12.80386pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{8.5359pt}{0.0pt}\pgfsys@lineto{68.28731pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-6.31065pt}{-2.5pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\ldots$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{-68.28731pt}{-17.07182pt}\pgfsys@lineto{-68.28731pt}{17.07182pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{-34.14365pt}{-17.07182pt}\pgfsys@lineto{-34.14365pt}{17.07182pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{34.14365pt}{-17.07182pt}\pgfsys@lineto{34.14365pt}{17.07182pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{68.28731pt}{-17.07182pt}\pgfsys@lineto{68.28731pt}{17.07182pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{-78.5305pt}{-8.5359pt}\pgfsys@moveto{-78.5305pt}{-8.5359pt}\pgfsys@lineto{-78.5305pt}{8.5359pt}\pgfsys@lineto{-58.0441pt}{8.5359pt}\pgfsys@lineto{-58.0441pt}{-8.5359pt}\pgfsys@closepath\pgfsys@moveto{-58.0441pt}{8.5359pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{-44.38684pt}{-8.5359pt}\pgfsys@moveto{-44.38684pt}{-8.5359pt}\pgfsys@lineto{-44.38684pt}{8.5359pt}\pgfsys@lineto{-23.90044pt}{8.5359pt}\pgfsys@lineto{-23.90044pt}{-8.5359pt}\pgfsys@closepath\pgfsys@moveto{-23.90044pt}{8.5359pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{17.07182pt}{-8.5359pt}\pgfsys@moveto{17.07182pt}{-8.5359pt}\pgfsys@lineto{17.07182pt}{8.5359pt}\pgfsys@lineto{51.21548pt}{8.5359pt}\pgfsys@lineto{51.21548pt}{-8.5359pt}\pgfsys@closepath\pgfsys@moveto{51.21548pt}{8.5359pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{58.0441pt}{-8.5359pt}\pgfsys@moveto{58.0441pt}{-8.5359pt}\pgfsys@lineto{58.0441pt}{8.5359pt}\pgfsys@lineto{78.5305pt}{8.5359pt}\pgfsys@lineto{78.5305pt}{-8.5359pt}\pgfsys@closepath\pgfsys@moveto{78.5305pt}{8.5359pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-71.25403pt}{-26.12303pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$s_{t_{k}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-38.54149pt}{-26.0297pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$s_{t_{k+1}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.18565pt}{-27.35774pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{\scriptsize$\ldots$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{28.3891pt}{-26.0297pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$s_{t_{k+m-1}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{63.4661pt}{-26.0297pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$s_{t_{k+m}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-70.46692pt}{24.548pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$s^{\prime}_{t_{k}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-36.32326pt}{24.64134pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$s^{\prime}_{t_{k+1}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.18565pt}{23.85774pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{\scriptsize$\ldots$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{31.96404pt}{24.64134pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$s^{\prime}_{t_{k+m-1}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{66.1077pt}{24.64134pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$s^{\prime}_{t_{k+m}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-73.62212pt}{-3.12666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{\scriptsize$B^{(0)}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-39.47845pt}{-3.12666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{\scriptsize$B^{(1)}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{26.4346pt}{-3.12666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\scriptsize$B^{(m-1)}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{62.2116pt}{-3.12666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{\scriptsize$B^{(m)}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$

where we first decompose $I_{m}$ into two tensors (eg. via SVD or QR
decomposition) and then define the tensors in the MPO as

$\displaystyle B^{(0)}_{v}(s^{\prime},s)$ $\displaystyle=\hskip
4.26773pt\sum_{j}\left(\delta_{s,s^{\prime}}\delta_{s,v}\right)\left[I^{(0)}_{m}(s)\right]_{v}$
(31) $\displaystyle B^{(m)}_{v}(s^{\prime},s)$ $\displaystyle=\hskip
4.26773pt\sum_{j}\left(\delta_{s,s^{\prime}}\delta_{s,v}\right)\left[I^{(1)}_{m}(s^{\prime})\right]_{v}$
(32) $\displaystyle B^{(i)}_{v,v^{\prime}}(s^{\prime},s)$
$\displaystyle=\hskip
4.26773pt\delta_{v,v^{\prime}}\delta_{s,s^{\prime}}\hskip
14.22636pt\forall\,i\in[1,...,m-1]\,\,.$ (33)

In diagrammatic form, the MPS for an IF with 6 time steps is

In our contraction scheme, we start from the bottom row and contract upwards.
However, because each $I_{m}$ factor commutes with the rest, other choices of
ordering are possible and are the basis of algorithms such as TEMPO Strathearn
_et al._ (2018); Jorgensen and Pollock (2019).

##  Appendix C Influence functional transverse contraction around an arbitrary
site

Sometimes the site whose dynamics we are interested in may be at the middle of
the MPS representation of the system (e.g. in the hard core boson model).
Thus, we need to generalize the tensor network diagrams presented in the main
text to consider IFs for subsystems at arbitrary lattice site $i$.

Figure 11: (left) Tensor network showing generalized IF isolating dynamics at
the $i=3$ lattice site. Triangles denote gauging of tensors along the row, as
defined in Fig. 3 in the main text. (right) Tensor network computing
expectation value using [blue] left and right environment columns obtained
separately using the iterative contraction scheme described above, [white]
original tensors at site $i$ dictating interactions of both environment
columns with site itself and the environments with each other, and [gray] on-
site terms including [circle] the initial state and [square] time evolution
operators and the observable of interest (see Fig. 1). It is cheapest to
contract this network vertically from the row at one end and continuing to the
other end. Note that in using this method one does not explicitly compute the
IF itself.

If the Hamiltonian only consists of nearest-neighbor interactions, the IFs
from the two sides of site $i$ are separable and can be computed
independently. Otherwise, the tensor network can be initialized as shown in
Fig. 11, and one contracts inwards from the outer columns separately. Once
only the column corresponding to the isolated site is left, one can now
include on-site terms (initial state, on-site time evolution operators,
observable) such that the network now corresponds to the expectation value of
the observable at the desired time step (a scalar). The cost of contracting
this network scales like $\mathcal{O}(D_{I}^{3})$, which is much cheaper than
explicitly computing the full IF first and then computing the observable
expectation values.

## References

  * Breuer and Petruccione (2002) H. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, 2002).
  * Feynman (1948) R. P. Feynman, Rev. Mod. Phys. 20, 367 (1948).
  * Nakajima (1958) S. Nakajima, Progress Theor. Phys. 20, 948 (1958).
  * Zwanzig (1960) R. Zwanzig, J. Chem. Phys. 33, 1338 (1960).
  * Breuer  _et al._ (2001) H. Breuer, B. Kappler, and F. Petruccione, Ann. Phys. 291, 36 (2001).
  * Feynman and Vernon (1963) R. P. Feynman and F. L. Vernon, Annals Phys. 24, 118 (1963).
  * Vernon (1959) F. L. Vernon,  _The Theory of a General Quantum System Interacting with a Linear Dissipative System_ , Ph.D. thesis, California Institute of Technology (1959).
  * Makri and Makarov (1995) N. Makri and D. E. Makarov, J. Chem. Phys 102, 4600 (1995).
  * Makri (1999) N. Makri, J. Chem. Phys. 111, 6164 (1999).
  * Strathearn  _et al._ (2017) A. Strathearn, B. W. Lovett, and P. Kirton, New J. Phys. 19, 093009 (2017).
  * Weiss  _et al._ (2008) S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Phys. Rev. B 77, 195316 (2008).
  * Segal  _et al._ (2010) D. Segal, A. J. MIllis, and D. R. Reichman, Phys. Rev. B 82, 205323 (2010).
  * Cohen  _et al._ (2015) G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Phys. Rev. Lett. 115, 266802 (2015).
  * Chen  _et al._ (2017) H.-T. Chen, G. Cohen, and D. R. Reichman, J. Chem. Phys. 146, 054105 (2017).
  * Cerillo and Cao (2014) J. Cerillo and J. Cao, Phys. Rev. Lett. 112, 110401 (2014).
  * Rosenbach  _et al._ (2016) R. Rosenbach, J. Cerrillo, S. F. Huelga, J. Cao, and M. B. Plenio, New J. Phys. 18, 023035 (2016).
  * Shi and Geva (2003) Q. Shi and E. Geva, J. Chem. Phys. 119, 12063 (2003).
  * Kelly  _et al._ (2016) A. Kelly, A. Montoya-Castillo, L. Wang, and T. E. Markland, J. Chem. Phys. 144, 184105 (2016).
  * Cohen and Rabani (2011) G. Cohen and E. Rabani, Phys. Rev. B 84, 075150 (2011).
  * Kidon  _et al._ (2018) L. Kidon, H. Wang, M. Thoss, and E. Rabani, J. Chem. Phys. 149, 104105 (2018).
  * Chatterjee and Makri (2019) S. Chatterjee and N. Makri, J. Phys. Chem. B 123, 10470 (2019).
  * Lerose  _et al._ (2021) A. Lerose, M. Sonner, and D. A. Abanin, Phys. Rev. X , 021040 (2021).
  * Strathearn  _et al._ (2018) A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, Nature Communications 9, 3322 (2018).
  * Jorgensen and Pollock (2019) M. R. Jorgensen and F. A. Pollock, Phys. Rev. Lett. 123, 240602 (2019).
  * Pollock  _et al._ (2018) F. A. Pollock, C. Rodriguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Phys. Rev. A 97, 012127 (2018).
  * Cygorek  _et al._ (2021) M. Cygorek, M. Cosacchi, A. Vagov, V. M. Axt, B. W. Lovett, J. Keeling, and E. M. Gauger, arXiv quant-ph , 2101.01653 (2021).
  * Bañuls  _et al._ (2009) M. C. Bañuls, M. B. Hastings, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 102, 240603 (2009).
  * Müller-Hermes  _et al._ (2012) A. Müller-Hermes, J. I. Cirac, and M. C. Bañuls, New J. Phys. 14, 075003 (2012).
  * Luchnikov  _et al._ (2019) I. A. Luchnikov, S. V. Vintskevich, H. Ouerdane, and S. N. Filippov, Phys. Rev. Lett. 122, 160401 (2019).
  * Luchnikov  _et al._ (2020) I. A. Luchnikov, S. V. Vintskevich, D. A. Grigoriev, and S. N. Filippov, Phys. Rev. Lett. 124, 140502 (2020).
  * Makri (2018a) N. Makri, J. Chem. Phys. 149, 214108 (2018a).
  * Makri (2021) N. Makri, Phys. Chem. Chem. Phys. 23, 12537 (2021).
  * Kundu and Makri (2019) S. Kundu and N. Makri, J. Chem. Phys. 151, 074110 (2019).
  * Kundu and Makri (2021) S. Kundu and N. Makri, J. Phys. Chem. C 125, 201 (2021).
  * Makri (2018b) N. Makri, J. Chem. Phys. 148, 101101 (2018b).
  * Kundu and Makri (2020) S. Kundu and N. Makri, Molecular Physics (2020), https://doi.org/10.1080/00268976.2020.1797200.
  * García-Ripoll (2006) J. J. García-Ripoll, New J. Phys. 8, 305 (2006).
  * Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
  * Daley  _et al._ (2004) A. J. Daley, C. Kollath, U. Schollwock, and G. Vidal, J. Stat. Mech.: Theor. Exp. 2004, P04005 (2004).
  * Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
  * White and Feiguin (2004) S. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
  * White  _et al._ (2018) C. D. White, M. Zalatel, R. S. K. Mong, and G. Refael, Phys. Rev. B 97, 035127 (2018).
  * Schollwöck and White (2006) U. Schollwöck and S. R. White, AIP Conference Proc. 816, 155 (2006).
  * Ronca  _et al._ (2017) E. Ronca, Z. Li, C. Jimenez-Hoyos, and G. K.-L. Chan, J. Chem. Theory Comput. 13, 5560 (2017).
  * Haegeman  _et al._ (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
  * Murg  _et al._ (2007) V. Murg, F. Verstraete, and J. I. Cirac, Phys. Rev. A 75, 033605 (2007).
  * Gray (2018) J. Gray, Journal of Open Source Software 3, 819 (2018).
  * Zwolak and Vidal (2004) M. Zwolak and G. Vidal, Phys. Rev. Lett. 93, 207205 (2004).
  * Verstraete  _et al._ (2004) F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
  * Barthel  _et al._ (2009) T. Barthel, U. Schollwöck, and S. R. WHite, Phys. Rev. B 79, 245101 (2009).
  * Karrasch  _et al._ (2012) C. Karrasch, J. H. Bardarson, and J. E. Moore, Phys. Rev. Lett. 108, 227206 (2012).
  * Vega and Alonso (2017) I. d. Vega and D. Alonso, Rev. Mod. Phys. 89, 015001 (2017).
  * Makri (2014) N. Makri, J. Chem. Phys. 141, 134117 (2014).
  * Makri (2020a) N. Makri, J. Chem. Theory Comput. 16, 4038 (2020a).
  * Wang and Thoss (2008) H. Wang and M. Thoss, New J. Phys. 10, 115005 (2008).
  * Makri (2020b) N. Makri, J. Chem. Phys. 152, 041104 (2020b).
  * Trotzky  _et al._ (2012) S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwock, J. Eisert, and I. Bloch, Nature Phys. 8, 325 (2012).