#  Tree tensor network state approach for solving hierarchical equations of
motions

Yaling Ke  Institute of Physics, University of Freiburg, Hermann-Herder-
Strasse 3, 79104 Freiburg, Germany

###### Abstract

The hierarchical equations of motion (HEOM) method is a numerically exact open
quantum system dynamics approach. The method is rooted in an exponential
expansion of the bath correlation function, which in essence strategically
reshapes a continuous environment into a set of effective bath modes that
allow for more efficient cutoff at finite temperatures. Based on this
understanding, one can map the HEOM method into a Schrödinger-like equation
with a non-Hermitian super Hamiltonian for an extended wavefunction being the
tensor product of the central system wave function and the Fock state of these
effective bath modes. Recognizing that the system and these effective bath
modes form a star-shaped entanglement structure, in this work, we explore the
possibility of representing the extended wave function as an efficient tree
tensor network state (TTNS), the super Hamiltonian as a tree tensor network
operator of the same structure, as well as the application of a time
propagation algorithm using the time-dependent variational principle. Our
benchmark calculations based on the spin-boson model with a slow-relaxing bath
show that, the proposed HEOM+TTNS approach yields consistent results with that
of the conventional HEOM method, while the computation is considerably sped up
by a factor of a few orders of magnitude. Besides, the simulation with a
genuine TTNS is four times faster than a one-dimensional matrix product state
decomposition scheme.

The efficient simulation of open quantum system dynamics lies at the heart of
a great variety of chemical and physical applications,Breuer, Petruccione  _et
al._ (2002); May and Kühn (2008); Weimer, Kshetrimayum, and Orús (2021)
ranging from efficient energy and charge transfer in natural and artificial
photosynthetic systems to quantum transport in nanoscale electronic
devices.Ishizaki and Fleming (2009); Brédas, Sargent, and Scholes (2017);
Thoss and Evers (2018); Cuevas and Scheer (2010) Over the past few decades,
the hierarchical equations of motion (HEOM) method has grown into a mainstream
non-perturbative and non-Markovian open quantum system dynamics
approach.Tanimura and Kubo (1989); Yan  _et al._ (2004); Xu and Yan (2007);
Jin, Zheng, and Yan (2008); Yan (2014); Schinabeck, Härtle, and Thoss (2018);
Hsieh and Cao (2018); Ye  _et al._ (2016); Tanimura (2020) The core idea of
the method is based on an exponential series expansion of the bath correlation
function as well as successively taking time-derivative of the reduced system
density operator in the influence functional formalism, so as to construct a
group of auxiliary density operators (ADOs) that obey a hierarchical set of
differential equations.

While the HEOM method has gained tremendous success aided by abundant
optimization schemes and advances in computational architecture,Ishizaki and
Tanimura (2005); Shi  _et al._ (2009); Hu, Xu, and Yan (2010); Strümpfer and
Schulten (2012); Tsuchimoto and Tanimura (2015); Kramer  _et al._ (2018); Chen
_et al._ (2022); Xu  _et al._ (2022a); Ikeda and Nakayama (2022) there remains
a computational bottleneck when the central system size is large and many
exponential terms are required to converge the original bath correlation
function. Recently, these limitations have been effectively broached by
employing the matrix product states (MPS),Schollwöck (2011) also known as
tensor trains (TT) decompositionOseledets (2011) of these ADOs.Shi  _et al._
(2018); Borrelli (2019); Yan, Xing, and Shi (2020); Ke, Borrelli, and Thoss
(2022) MPS/TT provides a compact and ideal way to encode one-dimensional (1D)
short-range quantum many-body system entanglement,Vidal (2004) by representing
the multi-dimensional wave function as the product of a chain of low-rank
tensors. However, it may not be the optimal choice for the HEOM method, given
that the system degree of freedom (DoF) and effective bath modes form a star-
shaped entanglement structure, i.e., the system DoF is coupled to every bath
mode and these bath modes are not directly coupled to each other. Therefore,
it is important to go beyond the MPS/TT ansatz and probe other higher-
dimensional tensor network state structures that resemble the inherent
entanglement topology between the system and effective bath modes.

In this work, we will demonstrate the applicability of tree tensor network
state (TTNS) which is a generalization of the MPS/TT into a tree-shaped
network of tensors,Shi, Duan, and Vidal (2006); Tagliacozzo, Evenbly, and
Vidal (2009); Murg  _et al._ (2010); Li, von Delft, and Xiang (2012);
Changlani  _et al._ (2013); Nakatani and Chan (2013); Murg  _et al._ (2015);
Gunst  _et al._ (2018); Schröder  _et al._ (2019); Larsson (2019); Ferrari,
Magnifico, and Montangero (2022); Seitz  _et al._ (2022); Milsted  _et al._
(2019) and a time propagation algorithm based on the time-dependent
variational principle (TDVP)Bauernfeind and Aichhorn (2020); Kloss, Reichman,
and Bar Lev (2020); Ceruti, Lubich, and Walach (2021) for solving the HEOM
method.

In general, an open quantum system describes a central system of interest
coupled to a macroscopic environment, and the Hamiltonian reads

$H=H_{S}+H_{B}+H_{SB},$ (1)

where $H_{S}$, $H_{B}$, and $H_{SB}$ correspond to the system, bath
Hamiltonian, and their interaction, respectively. As a demonstration of the
concept, we start by considering a simple but paradigmatic open quantum
system, the spin-boson model, where the system Hamiltonian is given by

$H_{S}=\epsilon\sigma_{z}+\Delta\sigma_{x}.$ (2)

Here, the energy bias between two spin states is given by $2\epsilon$ and
$\Delta$ denotes the coupling between two states. The spin is coupled to a
dissipative bosonic environment, modeled as a phonon bath comprised of an
infinite number of harmonic oscillators,

$H_{B}=\sum_{j}\left(\frac{p_{j}^{2}}{2}+\frac{1}{2}\omega_{j}^{2}q_{j}^{2}\right),$
(3)

where $q_{j}$ and $p_{j}$ are the mass-weighted position and conjugated
momentum operator of the $j$th harmonic oscillator, and $\omega_{j}$ is the
corresponding frequency. The system-bath coupling Hamiltonian $H_{SB}$ takes a
linear form with respect to the bath coordinate $q_{j}$,

$H_{SB}=\sigma_{z}\sum_{j}c_{j}q_{j},$ (4)

and the system-bath coupling operator in the system subspace is given by
$\sigma_{z}$, the coupling strength specified by $c_{j}$.

The composite system-bath dynamics are described by the density operator
$\rho(t)$. Note that, we assume the system and bath are factorized at the
initial moment and the bath is in its own thermal equilibrium state at inverse
temperature $\beta$. Then, the initial density operator $\rho(0)$ is given by

$\rho(0)=\rho_{S}(0)\otimes\frac{e^{-\beta H_{B}}}{\mathrm{Tr}_{B}\\{e^{-\beta
H_{B}}\\}}.$ (5)

In this case, in the reduced description of the system dynamics
$\rho_{S}(t)=\mathrm{tr}_{B}\\{\rho(t)\\}$ where all environmental DoFs are
traced out, the influence of the environment on the system is characterized
statistically in a time correlation function

$C(t)=\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{e^{-i\omega
t}}{1-e^{-\beta\omega}}J(\omega)\mathrm{d}\omega.$ (6)

The spectral density function $J(\omega)$ encodes the coupling-weighted
density of states of the bath in frequency space, defined by

$J(\omega)=\frac{\pi}{2}\sum_{j}\frac{c_{j}^{2}}{\omega_{j}}\delta(\omega-\omega_{j}).$
(7)

For some specific forms of $J(\omega)$, such as Drude-Lorentz spectral density
function,

$J(\omega)=2\lambda\frac{\omega\Omega}{\omega^{2}+\Omega^{2}},$ (8)

where $\Omega$ is the bath characteristic frequency and $\lambda$ quantifies
the system-bath coupling strength, the bath correlation function in Eq. (6)
can be expanded analytically into an exponential summation,

$C(t)=\sum_{p=0}^{\infty}\lambda\eta_{p}e^{-i\gamma_{p}t}.$ (9)

The explicit expressions of $\eta_{p}$ and $\gamma_{p}$ depend on the sum-
over-pole decomposition scheme of the Bose distribution function
$f(\omega)=1/(1-e^{-\beta\omega})$. Throughout this work, we adopt the Padé
pole decomposition scheme Hu, Xu, and Yan (2010); Hu  _et al._ (2011) and the
explicit expressions for $\eta_{p}$ and $\gamma_{p}$ can be found in the
supporting information (SI). At finite temperatures, the infinite summation in
Eq. (9) can be cut off at a finite $P$, which is chosen large enough to well
reproduce the original correlation function. For an arbitrary form of the
spectral density function, the exponential expansion analogous to Eq. (9) can
be implemented numerically.

In fact, the exponential expansion in Eq. (9) lays the groundwork for the
derivations of the HEOM method. For a detailed derivation of the method, we
refer readers to a review Ref. Tanimura, 2020 and the references therein.
Here, we only briefly introduce the main concept of the method. The key is to
introduce a group of ADOs, $\rho^{\bm{n}}(t)$, where the superscript
$\bm{n}=(n_{0},\cdots,n_{p},\cdots,n_{P})$ is an ordered set of $P+1$ non-
negative integers and each $n_{p}$ is associated with a frequency component
$\gamma_{p}$ in Eq. (9). These ADOs are closed with respect to the time-
derivative operation, which yields

a)

b)

c)

Figure 1: Examples of the MPS/TT (a) and two TTNS decompositions (b and c) in
the graphical notation of the extended wave function $|\Psi(t)\rangle$ (see
Eq. (14)) with $P=3$ for the spin-boson model. The circle nodes represent the
tensors, black dangling lines denote the physical indices, and shared links
correspond to virtual indices. The orange dotted lines with arrows specify the
ordered sequence sweeping through all nodes in the entire TTNS, which start
and end at the root node.

$\begin{split}i\frac{d\rho^{\bm{n}}(t)}{dt}=&H_{S}\rho^{\bm{n}}(t)-\rho^{\bm{n}}(t)H_{S}-i\sum_{p=0}^{P}n_{p}\gamma_{p}\rho^{\bm{n}}(t)\\\
&+\sum_{p=0}^{P}\sqrt{\lambda(n_{p}+1)}\left(\sigma_{z}\rho^{\bm{n}^{+}_{p}}(t)-\rho^{\bm{n}_{p}^{+}}(t)\sigma_{z}\right)\\\
&+\sum_{p=0}^{P}\sqrt{\lambda
n_{p}}\left(\eta_{p}\sigma_{z}\rho^{\bm{n}_{p}^{-}}(t)-\eta_{p}^{*}\rho^{\bm{n}^{-}_{p}}(t)\sigma_{z}\right),\end{split}$
(10)

where $\bm{n}_{p}^{\pm}=(n_{0},\cdots,n_{p}\pm 1,\cdots,n_{P})$. Directly
propagating Eq. (10) with a certain hierarchy truncation scheme is termed as
the conventional HEOM method. The simplest truncation scheme is to set
$\rho^{\bm{n}}(t)=0$ when $\sum_{p}n_{p}>L$ and $L$ is the hierarchical
truncation tier. Using this truncation scheme, the number of elements to be
stored and propagated is $\frac{(P+1+L)!}{(P+1)!L!}d^{2N_{s}}$, where $N_{s}$
is the number of system DoFs and $d$ is the size of basis set per system DoF
(for the spin-boson model, we have $N_{s}=1$ and $d=2$). As such, the
conventional HEOM method is formidably expensive in the cases of many bath
poles, a deep truncation tier, and a large central system. An alternative way
to reduce the computational cost and memory requirement is to employ tensor
network states. To this end, it is opportune to reformulate the above
hierarchical set of differential equations into a Schrödinger-like
equation,Borrelli (2019); Ke, Borrelli, and Thoss (2022)

$i\frac{d|\Psi(t)\rangle}{dt}=\mathcal{H}|\Psi(t)\rangle,$ (11)

via the introduction of an effective phonon bath and the density matrix
purification scheme.Schmutz (1978); Suzuki (1991); Arimitsu and Umezawa
(1987); Feiguin and White (2005); Verstraete, Garcia-Ripoll, and Cirac (2004);
Borrelli and Gelin (2021)

First, for each $n_{p}$, a virtual dissipative harmonic oscillator is
introduced with the basis $|n_{p}\rangle$ in the occupation number
representation as well as a pair of creation and annihilation operators,
$b^{\dagger}_{p}$ and $b_{p}$, acting on $|n_{p}\rangle$ to yield

$\displaystyle b^{\dagger}_{p}$
$\displaystyle|n_{p}\rangle=\sqrt{n_{p}+1}|n_{p}+1\rangle,$ (12a)
$\displaystyle b_{p}$
$\displaystyle|n_{p}\rangle=\sqrt{n_{p}}|n_{p}-1\rangle.$ (12b)

All these virtual modes constitute an effective phonon bath, whose influence
on the system dynamics is exactly the same as that of the original bath (see
Eqs. (3) and (4)). When all these effective modes are placed in their ground
states, i.e., for the particular configuration of the Fock state
$|\bm{n}\rangle=|0,\cdots,0\rangle$, the reduced density operator is
reproduced, $\rho_{S}(t)=\rho^{(0,\cdots,0)}(t)$.

Second, each ADO is usually expressed as a density matrix in the Hilbert
space, but it can also be recast into a rank-$N_{s}$ tensor ,

$\rho^{\bm{n}}(t)\equiv\sum_{s_{1}\cdots s_{N_{s}}}C^{\bm{n}}_{s_{1}\cdots
s_{N_{s}}}(t)|s_{1}\cdots s_{N_{s}}\rangle.$ (13)

where $s_{i}$ ranges from $1$ to $d^{2}$. For more theoretical and technical
details with regard to this space transformation, we refer readers to Refs.
Schmutz, 1978; Suzuki, 1991; Arimitsu and Umezawa, 1987; Feiguin and White,
2005; Verstraete, Garcia-Ripoll, and Cirac, 2004; Borrelli and Gelin, 2021.

All the ADOs combined comprise an extended wave function in the enlarged space

$|\Psi(t)\rangle=\sum_{\begin{subarray}{c}s_{1}\cdots s_{N_{s}}\\\ n_{0}\cdots
n_{P}\end{subarray}}C^{\bm{n}}_{\bm{s}}(t)|s_{1}\cdots
s_{N_{s}}\rangle\otimes|n_{0}\cdots n_{P}\rangle.$ (14)

The reduced system operator is obtained by the partial product
$\rho_{S}(t)\equiv|\Psi_{S}(t)\rangle=\langle\bm{n}=\bm{0}|\Psi(t)\rangle$.
The super Hamiltonian $\mathcal{H}$ in Eq. (11) can be written out explicitly
with the operators introduced in Eq. (12) as

$\begin{split}\mathcal{H}=&\mathcal{H}_{S}-i\sum_{p=0}^{P}\gamma_{p}b^{\dagger}_{p}b_{p}+\sqrt{\lambda}\hat{\sigma}_{z}\sum_{p=0}^{P}\left(b_{p}+\eta_{p}b^{\dagger}_{p}\right)\\\
&-\sqrt{\lambda}\tilde{\sigma}_{z}\sum_{p=0}^{P}\left(b_{p}+\eta_{p}^{*}b^{\dagger}_{p}\right).\end{split}$
(15)

where
$\mathcal{H}_{S}=\epsilon\left(\hat{\sigma}_{z}-\tilde{\sigma}_{z}\right)+\Delta\left(\hat{\sigma}_{x}-\tilde{\sigma}_{x}\right)$
with $\hat{\sigma}_{z/x}=\sigma_{z/x}\otimes\mathbbm{1}_{2}$,
$\tilde{\sigma}_{z/x}=\mathbbm{1}_{2}\otimes\sigma_{z/x}$, and
$\mathbbm{1}_{2}$ being an $2\times 2$ unit matrix. It is worth noting that
the super Hamiltonian $\mathcal{H}$ is non-Hermitian. As such, the norm of the
extended wave function $|\Psi(t)\rangle$ is not conserved, but the norm
conservation of reduced system wave function $|\Psi_{S}(t)\rangle$ is in
principle guaranteed.

The direct time integration of Eq. (11) along with the explicit definition of
the extended wave function $|\Psi(t)\rangle$ and super Hamiltonian
$\mathcal{H}$ given in Eq. (14) and Eq. (15), respectively, is usually an
intractable task as the size of the coefficient tensor
$C^{\bm{n}}_{\bm{s}}(t)$ scales exponentially with the number of extended
system DoFs, $N=N_{s}+P+1$. This huge amount of data can be compressed in the
format of a tensor network state, which is a structured product of low-rank
tensors.

The well-established MPS/TT format is a special and widely-used instance of
tensor network states, where the tensors are arranged into a one-dimensional
chain.Schollwöck (2011); Oseledets (2011); Cirac and Verstraete (2009); Chan
_et al._ (2016); Jaschke, Wall, and Carr (2018) The high-rank coefficient
tensor $C^{\bm{n}}_{\bm{s}}(t)$ in the extended wave function
$|\Psi(t)\rangle$ can be approximated in the MPS/TT formalism as

$\begin{split}C_{\bm{s}}^{\bm{n}}(t)&\approx\sum_{r_{1}\cdots
r_{N-1}}A^{[1]}_{r_{1},s_{1}}\cdots
A^{[N_{s}+1]}_{r_{N_{s}},r_{N_{s}+1},n_{0}}\cdots A^{[N]}_{r_{N-1},n_{P}}\\\
&\approx\mathrm{Contr}\left\\{A^{[1]}_{s_{1}}\cdots
A^{[N_{s}+1]}_{n_{0}}\cdots A^{[N]}_{n_{P}}\right\\}.\end{split}$ (16)

Here, $A^{[i]}$ are rank-3 tensors with one physical index ($s_{i}$ or
$n_{p}$) and two virtual indices $r_{i-1}$ and $r_{i}$, except for $A^{[1]}$
and $A^{[N]}$ which are rank-2 tensors and thus have only one virtual index.
The virtual index $r_{i}$ runs from 1 to $D_{i}$, where the bond dimension
$D_{i}$ is a controllable parameter that can be systematically increased to
reduce the degree of approximation in Eq. (16). The maximum value among
$\\{D_{i}\\}$ is designated as the maximal bond dimension $D_{\mathrm{max}}$.
An example of the MPS/TT representation of $|\Psi(t)\rangle$ with $P=3$ is
illustrated graphically in Fig. 1 (a). The circle nodes correspond to the
tensors. The dangling legs denote the physical indices, and the connected link
between two neighboring nodes is assigned with a virtual index $r_{i}$. The
summation over a virtual index shared between two tensors is called the
contraction. Contracting all virtual indices (denoted in Eq. (16) as
$\mathrm{Contr}$) reproduces approximately $C_{\bm{s}}^{\bm{n}}(t)$.

The TTNS is an extension of the MPS/TT beyond 1D geometry and thus can
accommodate a more complicated entanglement network.Shi, Duan, and Vidal
(2006); Tagliacozzo, Evenbly, and Vidal (2009); Murg  _et al._ (2010); Li, von
Delft, and Xiang (2012); Changlani  _et al._ (2013); Nakatani and Chan (2013);
Murg  _et al._ (2015); Gunst  _et al._ (2018); Schröder  _et al._ (2019);
Larsson (2019); Ferrari, Magnifico, and Montangero (2022); Seitz  _et al._
(2022); Milsted  _et al._ (2019) It is also called hierarchical Tucker
format,Lubich  _et al._ (2013) which is the ansatz behind an important quantum
dynamics approach, multi configuration time-dependent Hartree and its multi-
layer variant.Beck and Meyer (1997); Wang and Thoss (2003); Meyer, Gatti, and
Worth (2009); Manthe (2015); Wang (2015); Lindoy, Kloss, and Reichman (2021) A
general TTNS connects all tensors into a tree structure and the network is
loop-free,Montangero, Rico, and Silvi (2022) which means that there exists a
unique path between any two tensors. There are many different ways to connect
the tensors into a tree-shaped network. Optimally, the tree shape should be in
alignment with the natural entanglement topology of the whole system, which is
determined by the super Hamiltonian $\mathcal{H}$ in Eq. (15). On the one
hand, all effective bath modes are not directly coupled to each other but to
the system DoF. On the other, the influence of a specific effective bath mode
on the system dynamics decreases with an increasing pole number $p$. Based on
these considerations, we construct a TTNS representation for
$C_{\bm{s}}^{\bm{n}}(t)$ in the extended wave function $|\Psi(t)\rangle$,
which is directly analogous to Eq. (16),

$\begin{split}C_{\bm{s}}^{\bm{n}}(t)\approx\mathrm{Contr}\left\\{T^{[1]}_{s_{1}}T^{[2]}_{n_{0}}T^{[3]}_{X}\cdots
T^{[k]}\cdots T^{[K]}_{n_{P}}\right\\},\end{split}$ (17)

except that $T^{[k]}$ is now a rank-$z$ tensor where $z$ can be higher than 3.
The useful graphical representation for an example with $P=3$ is shown in Fig.
1 (b). Note in passing that, in the cases where the environment consists of
discrete vibrational modes and each mode contributes equally to the impact on
the system dynamics, a more balanced structure as shown in Fig. 1 (c) can be
taken.Yan  _et al._ (2021)

a)

b)

c)

Figure 2: (a) An example of the TTNO decomposition of the super Hamiltonian
$\mathcal{H}$ (see Eq. (15)), which has the same network structure as the TTNS
in Fig. 1 (b). To distinguish the TTNO from the TTNS in the graphical
notation, each tensor in the TTNO is represented as a rectangle node. (b) and
(c) display an example of the effective Hamiltonian $H_{[k]}^{\mathrm{eff}}$
(see Eq. (25)) for node $T^{[7]}$ and $G_{[km]}^{\mathrm{eff}}$ (see Eq. (26))
for the link between nodes $T^{[5]}$ and $T^{[7]}$, respectively. The nodes
with an asterisk denote taking their complex conjugation.

For convenience, the circle nodes in Fig. 1 (b) are enumerated and some
notations that will be used in what follows are introduced here. A node with
$z$ legs corresponds to a rank-$z$ tensor. Similar to an MPS/TT, there are
nodes assigned with one dangling leg denoting the physical index
$s_{1}/n_{p}$, such as nodes $T^{[1]}_{s_{1}}$, $T^{[2]}_{n_{0}}$,
$T^{[4]}_{n_{1}}$, $T^{[6]}_{n_{2}}$, and $T^{[7]}_{n_{3}}$. However, a TTNS
can also have connecting nodes, which have exclusively connected legs, i.e.,
virtual indices, like nodes $T^{[3]}_{X}$ and $T^{[5]}_{X}$. Therefore, the
number of nodes $K$ can be larger than $N$. We assign the red-shaded node with
one leg corresponding to the first system DoF as the root node, $T^{[1]}$, as
shown in Fig. 1 (b). Every node can have none or multiple children nodes. For
example, the root node has two children nodes, $T^{[2]}$ and $T^{[3]}$.
$T^{[2]}$ has no children node. Leaves are the nodes (green-shaded) without
children nodes. Other than the root node, every node has a parent node, which
is the adjacent node on the path pointing toward the root node. For each node,
all the legs are enumerated in a counter-clockwise manner with the bond
connected to the parent node designated as the first leg, and the dangling leg
(if it has) as the last leg, as shown exemplarily for node $T^{[5]}$ at the
top-right corner of Fig. 1 (b).

The super Hamiltonian $\mathcal{H}$ in Eq. (15) can also be decomposed into
the product of tree tensor network operators (TTNO), which have the same
network structure as the TTNS,

$\mathcal{H}=\mathrm{Contr}\left\\{W^{[1]}_{s_{1},s^{\prime}_{1}}W^{[2]}_{n_{0},n^{\prime}_{0}}\left(\prod_{p=1}^{P-1}W_{X}^{[2p+1]}W^{[2p+2]}_{n_{p},n_{p}^{\prime}}\right)W^{[2P+1]}_{n_{P},n^{\prime}_{P}}\right\\},$
(18)

as exhibited graphically in Fig. 2 (a) for $P=3$. The TTNO tensors
$W^{[1]}_{s_{1},s^{\prime}_{1}}$ and $W^{[2p+2]}_{n_{p},n^{\prime}_{p}}$ in
Eq. (18) read

$W^{[1]}_{s_{1},s^{\prime}_{1}}=\left(\begin{array}[]{cccc}\mathcal{H}_{S}&\hat{\sigma}_{z}&-\tilde{\sigma}_{z}&\mathbbm{1}_{4}\\\
\hat{\sigma}_{z}&&&\\\ -\tilde{\sigma}_{z}&&&\\\ \mathbbm{1}_{4}&&&\\\
\end{array}\right),$ (19)
$W^{[2p+2]}_{n_{p},n^{\prime}_{p}}=\left(\begin{array}[]{c}\mathbbm{1}_{L}\\\
\sqrt{\lambda}\left(b_{p}+\eta_{p}b^{\dagger}_{p}\right)\\\
\sqrt{\lambda}\left(b_{p}+\eta^{*}_{p}b^{\dagger}_{p}\right)\\\
-i\gamma_{p}b^{\dagger}_{p}b_{p}\\\ \end{array}\right),$ (20)

where the terms equal to zero have been left blank in Eq. (19).
$\mathbbm{1}_{L}$ is an $L\times L$ unit matrix. $W_{X}$ is a rank-3
connecting tensor of size $4\times 4\times 4$ with the elements

$W_{X}[i,j,k]=\frac{\delta_{ij}\delta_{k=1}+\delta_{ik}\delta_{j=1}}{2\delta_{i=1}+\delta_{i\neq
1}}.$ (21)

With the TTNS and TTNO expression for the extended wave function
$|\Psi(t)\rangle$ and the super Hamiltonian $\mathcal{H}$ explicitly given in
Eqs. 17 and 18 respectively, we can use the TDVP-based time propagation
algorithm to solve Eq. (11). TDVP algorithm has been proven to be a stable and
efficient time propagation scheme for the MPS/TT.Haegeman, Osborne, and
Verstraete (2013); Lubich, Oseledets, and Vandereycken (2015); Haegeman  _et
al._ (2016); Schröder and Chin (2016); Baiardi and Reiher (2019); Mendl
(2018); Paeckel  _et al._ (2019) The feasibility of the TDVP algorithm on the
TTNS, the detailed derivation and related numerical analyses have been
reported in Refs. Bauernfeind and Aichhorn, 2020; Kloss, Reichman, and Bar
Lev, 2020; Ceruti, Lubich, and Walach, 2021. Here, we only concisely discuss
the implementation of the algorithm for our model and provide a pseudocode.

In the one-site version of the TDVP algorithm, the Schrödinger equation in Eq.
(11) is solved by projecting the wave function into the so-called tangent
space, i.e., a manifold of all TTNSs with fixed bond dimensions, Bauernfeind
and Aichhorn (2020); Kloss, Reichman, and Bar Lev (2020); Ceruti, Lubich, and
Walach (2021)

$i\frac{|\Psi[T]\rangle}{\partial
t}=\mathcal{P}_{\mathcal{M}\\{\Psi[T]\\}}\mathcal{H}|\Psi[T]\rangle.$ (22)

The explicit expression of the tangent space projection operator
$\mathcal{P}_{\mathcal{M}\\{\Psi[T]\\}}$ can be found in Refs. Bauernfeind and
Aichhorn, 2020; Kloss, Reichman, and Bar Lev, 2020; Ceruti, Lubich, and
Walach, 2021. By using Trotter breakups, the projected time-dependent
Schrödinger equation in Eq. (22) can be integrated locally, where each tensor
$T^{[k]}$ and the link $V^{[km]}$ between node $k$ and one of its neighboring
node $m$ in a specified sequence are updated according to the following
equations

$\frac{\partial T^{[k]}(t)}{\partial t}=-iH_{[k]}^{\mathrm{eff}}T^{[k]}(t),$
(23) $\frac{\partial V^{[km]}(t)}{\partial
t}=iG_{[km]}^{\mathrm{eff}}V^{[km]}(t),$ (24)

while all other tensors are fixed. The effective local Hamiltonian
$H_{[k]}^{\mathrm{eff}}$ in Eq. (23) and $G_{[km]}^{\mathrm{eff}}$ in Eq. (24)
are given by

$\begin{split}H_{[k]}^{\mathrm{eff}}=\mathrm{Contr}^{\prime}_{k}&\left\\{T^{[1]}_{s_{1}}W^{[1]}_{s_{1},s^{\prime}_{1}}T^{[1]*}_{s^{\prime}_{1}}\cdots
T^{[k-1]}W^{[k-1]}T^{[k-1]*}W^{[k]}\right.\\\
&\left.T^{[k+1]}W^{[k+1]}T^{[k+1]*}\cdots
T^{[K]}_{n_{P}}W^{[K]}_{n_{P},n^{\prime}_{P}}T^{[K]*}_{n^{\prime}_{P}}\right\\},\end{split}$
(25)
$\begin{split}G_{[km]}^{\mathrm{eff}}=\mathrm{Contr}{{}^{\prime\prime}}_{km}&\left\\{T^{[1]}_{s_{1}}W^{[1]}_{s_{1},s^{\prime}_{1}}T^{[1]*}_{s^{\prime}_{1}}\cdots
T^{[k]}W^{[k]}T^{[k]*}\right.\\\ &\left.\cdots
T^{[K]}_{n_{P}}W^{[K]}_{n_{P},n^{\prime}_{P}}T^{[K]*}_{n^{\prime}_{P}}\right\\}.\end{split}$
(26)

Here, we use the Einstein convention that the duplicate indices are summed
over. The asterisk denotes taking the complex conjugation of the tensor.
$\mathrm{Contr}^{\prime}_{k}$ denotes contracting over all physical and
virtual indices in the TTNS decomposition of $C_{\bm{s}}^{\bm{n}}(t)$ and
$C_{\bm{s}}^{\bm{n}*}(t)$ as well as the TTNO expression for $\mathcal{H}$,
except for leaving the indices connected to $T^{[k]}$ and $T^{[k]*}$ and
physical indices to $W^{[k]}$ open. For
$\mathrm{Contr}{{}^{\prime\prime}}_{[km]}$, the bond connecting $T^{[k]}$ and
$T^{[m]}$ as well as the bond connecting $T^{[k]*}$ and $T^{[m]*}$ are left
open. As an example, we show the graphical representation of
$H_{[7]}^{\mathrm{eff}}$ and $G_{[57]}^{\mathrm{eff}}$ in Fig. 2 (b) and (c),
respectively.

A full TDVP step starts from a canonicalized TTNS where all tensors except for
the root node are orthonormalized,
$\bm{T}^{[1]}(t),\mathcal{T}^{[2]}(t),\cdots,\mathcal{T}^{[K]}(t)$, and is
accomplished when all tensors are propagated to $\bm{T}^{[1]}(t+\Delta
t),\mathcal{T}^{[2]}(t+\Delta t),\cdots,\mathcal{T}^{[K]}(t+\Delta t)$ after a
sweeping along the entire network, as illustrated by the orange dotted lines
in Fig. 1. A full sweeping procedure starts and ends at the root node, and
every node is entered and exited $n_{c}^{k}+1$ times where $n_{c}^{k}$ is the
number of the children nodes to node $k$. For example, a sweep over the TTNS
in Fig. 1 (b) goes through the nodes sequentially, and the path is
$1\rightarrow 2\rightarrow 1\rightarrow 3\rightarrow 4\rightarrow 3\rightarrow
5\rightarrow 6\rightarrow 5\rightarrow 7\rightarrow 5\rightarrow 3\rightarrow
1$. To keep track of the time evolution, we introduce a variable $h_{k}$ and
it is initialized as $h_{k}=0$ before the sweeping procedure starts. Moving
from node $k$ to a neighboring node $m$ along the direction of the sweeping
path, the following steps are implemented:

  1. 1.

Evolve $\bm{T}^{[k]}(t+h_{k})$ forward in time according to

$\bm{T}^{[k]}\left(t+h_{k}+\frac{\Delta
t}{n_{c}^{k}+1}\right)=e^{-iH^{\mathrm{eff}}_{[k]}\frac{\Delta
t}{n_{c}^{k}+1}}\bm{T}^{[k]}(t+h_{k}),$ (27)

and update $h_{k}=h_{k}+\frac{\Delta t}{n_{c}^{k}+1}$.
$H^{\mathrm{eff}}_{[k]}$ is computed with
$\mathcal{T}^{[1]}(t+h_{1}),\cdots,\mathcal{T}^{[k-1]}(t+h_{k-1}),\mathcal{T}^{[k+1]}(t+h_{k+1}),\cdots,\mathcal{T}^{[K]}(t+h_{K})$.

  2. 2.

Permute the indices of $\bm{T}^{[k]}(t+h_{k})$ and then reshape the tensor
into a matrix such that the index $r_{k,m}$ corresponding to the bond
connecting nodes $k$ and $m$ now serves as the column of the matrix, while all
other indices are grouped into one index as the row of the matrix; Conduct the
QR decomposition of the matricized $\bm{T}^{[k]}(t+h_{k})$ into
$\mathcal{T}^{[k]}(t+h_{k})$ and $V^{[km]}(t+h_{k})$; Retensorize
$\mathcal{T}^{[k]}(t+h_{k})$ as a new local tensor on node $k$.

  3. 3.

Evolve $V^{[km]}(t+h_{k})$ in time according to

$V^{[km]}\left(t+h_{m}\right)=e^{iG^{\mathrm{eff}}_{[km]}(h_{k}-h_{m})}\,V^{[km]}(t+h_{k}),$
(28)

where $G^{\mathrm{eff}}_{[km]}$ is computed with
$\mathcal{T}^{[1]}(t+h_{1}),\cdots,\mathcal{T}^{[k]}(t+h_{k}),\cdots,\mathcal{T}^{[K]}(t+h_{K})$.

  4. 4.

Absorb $V^{[km]}\left(t+h_{m}\right)$ into $\mathcal{T}^{[m]}(t+h_{m})$ to
obtain $\bm{T}^{[m]}(t+h_{m})$.

After implementing the above steps, the orthogonality center in the TTNS is
also moved from node $k$ to the adjacent node $m$.

To assess the accuracy and performance of the proposed HEOM+TTNS method
against the conventional HEOM method, we first conduct a benchmark calculation
on the spin-boson model. The spin is initially placed in the spin-up state.
The parameters can be found in the caption of Fig. 3. Note that, we adopt a
much smaller value for $\Omega$ than other parameters, which indicates that
the bath relaxation is slow as compared to other system dynamical processes,
and the non-Markovian feature in this case can be strong. We inspected the
spin dynamics $\langle\sigma_{z}(t)\rangle=\mathrm{tr}\\{\sigma_{z}\rho(t)\\}$
and found that the converged result using the conventional HEOM method is
obtained with a very deep truncation tier $L=\sum_{p=0}^{P}n_{p}=50$ and
$P=3$. Using the HEOM+TTNS method with the TTNS displayed in Fig. 1 (b),
$\langle\sigma_{z}(t)\rangle$ is calculated by

$\begin{split}\langle\sigma_{z}(t)\rangle=&\mathrm{Contr}\left\\{T^{[1]}_{s_{1}}\sigma_{z}[s_{1},s^{\prime}_{1}]T^{[1]*}_{s^{\prime}_{1}}T^{[2]}_{0}T^{[2]*}_{0}\right.\\\
&\left.\left(\prod_{p=1}^{P-1}T^{[2p+1]}_{X}T^{[2p+1]*}_{X}T^{[2p+2]}_{0}T^{[2p+2]*}_{0}\right)T^{[K]}_{0}T^{[K]*}_{0}\right\\},\end{split}$
(29)

where $T^{[2p+2]}_{0}$ denotes $T^{[2p+2]}_{n_{p}=0}$. For every effective
bath mode, we take $0\leq n_{p}\leq L$, which means that actually more ADOs
are taken into account than in the conventional HEOM method. We found that
with a very low maximal bond dimension $D_{max}=5$, the result obtained
through the HEOM+TTNS method is already in excellent agreement with that
obtained through the conventional HEOM method, as shown in Fig. 3.

For a fair comparison, all the codes are written using Julia programming
language and our simulations are run on a single CPU core of an Apple MacBook
pro laptop with M1 chip. For the conventional HEOM method, we use an
exponential integrator as proposed in Ref. Wilkins and Dattani, 2015, which
requires less memory usage as compared to Runge-Kutta methods. Every
simulation step takes on an average 6.7 seconds with a time step $\Delta
t=0.025\text{ fs}$ and requires a memory of 7.6 GB. Using the HEOM+TTNS method
with the previously stated TDVP integration scheme, we gain a significant
speed-up in computational time and save in memory usage. Every simulation step
takes only 0.028 seconds in CPU time and memory usage of 10 MB. Besides, it is
also found that the maximally allowed time step in the HEOM+TTNS method is
$\Delta t=0.5\text{ fs}$, much larger than that used in the conventional HEOM
method. In a word, even for the relatively simple spin-boson model, the
simulation using the HEOM+TTNS method can be a few orders of magnitude faster
than the conventional HEOM method.

We also tested the MPS/TT decomposition of the extended wave function
$|\Psi(t)\rangle$ for the above model (see the tensor network structure in
Fig. 1 (a)) and the matrix product operator format of the super Hamiltonian
$\mathcal{H}$ is obtained using ITensor Package.Fishman, White, and
Stoudenmire (2022) The same time propagation algorithm and code as in the
HEOM+TTNS method are used. For this model, although we found that the result
of the HEOM+MPS/TT method (data not shown) also converges with a maximal bond
dimension of $D_{max}=5$, every simulation step is four times slower than that
of the HEOM+TTNS method. This implies the pure influence of a tree shape in
determing the practical computational cost.

Figure 3: Dynamics of the spin-boson model calculated using the conventional
HEOM with a hierarchical truncation tier $L=50$ and the HEOM+TTNS method with
a maximal bond dimension $D_{max}=5$. The following parameters are used:
$\epsilon=50\text{ cm}^{-1}$, $\Delta=200\text{ cm}^{-1}$, $\beta^{-1}=300K$,
$\lambda=100\text{ cm}^{-1}$, and $\Omega=10\text{ cm}^{-1}$.

As a proof-of-concept, we demonstrated in the main text the algorithm and
results specific to the spin-boson model. The method can be generalized to
much more complicated systems. In the SI, we show the applicability of the
HEOM+TTNS method for the Frenkel exciton model which is widely used in
studying the excitation energy transfer in photosynthetic systems,Ishizaki and
Fleming (2009) as well as a spin-boson-like chain model where each spin is
coupled to its own bosonic environment and the nearest-neighbor spins.

In summary, we presented the feasibility of adopting a genuine TTNS and TDVP-
based time propagation algorithm for solving the HEOM method, and compared its
performance to the conventional HEOM method as well as the HEOM-MPS/TT scheme.
Because the TTNS inherits many properties from the MPS/TT formalism, we expect
that the proposed HEOM+TTNS method can also benefit from algorithmic
optimizations and advances that have been developed in the context of the
MPS/TT representation.Yang and White (2020); Dunnett and Chin (2021); Borrelli
and Dolgov (2021); Li, Gleis, and Von Delft (2022); Xu  _et al._ (2022b);
Ceruti, Lubich, and Sulz (2023) Besides, we also hope that the idea presented
in this work may inspire more applications of the tensor network states in
combination with the HEOM method to enable a more efficient simulation of a
much broader class of chemical and physical systems. For example, we will
extend the current approach to study open fermionic system problems in our
future work.Evers  _et al._ (2022)

## Acknowledgements

The author thanks M. Thoss for helpful discussions. This work was supported by
the German Research Foundation (DFG).

## Supplementary Material

See the supplementary material for the implementation details of the HEOM+TTNS
method for the Frenkel exciton model and a spin-boson-like chain model.

## Data Availability

The data and code that support the findings of this study are available from
the corresponding author upon reasonable request.

## References

  * Breuer, Petruccione  _et al._ (2002) H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * May and Kühn (2008) V. May and O. Kühn,  _Charge and energy transfer dynamics in molecular systems_ (John Wiley & Sons, 2008).
  * Weimer, Kshetrimayum, and Orús (2021) H. Weimer, A. Kshetrimayum, and R. Orús, “Simulation methods for open quantum many-body systems,” Rev. Mod. Phys. 93, 015008 (2021).
  * Ishizaki and Fleming (2009) A. Ishizaki and G. R. Fleming, “Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature,” Proceedings of the National Academy of Sciences 106, 17255–17260 (2009).
  * Brédas, Sargent, and Scholes (2017) J.-L. Brédas, E. H. Sargent, and G. D. Scholes, “Photovoltaic concepts inspired by coherence effects in photosynthetic systems,” Nat. Mater. 16, 35–44 (2017).
  * Thoss and Evers (2018) M. Thoss and F. Evers, “Perspective: Theory of quantum transport in molecular junctions,” J. Chem. Phys. 148, 030901 (2018).
  * Cuevas and Scheer (2010) J. C. Cuevas and E. Scheer,  _Molecular electronics: an introduction to theory and experiment_ (World Scientific, Singapore, 2010).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly gaussian-markoffian noise bath,” J. Phys. Soc. Jpn. 58, 101–114 (1989).
  * Yan  _et al._ (2004) Y.-a. Yan, F. Yang, Y. Liu, and J. Shao, “Hierarchical approach based on stochastic decoupling to dissipative systems,” Chem. Phys. Lett. 395, 216–221 (2004).
  * Xu and Yan (2007) R.-X. Xu and Y. Yan, “Dynamics of quantum dissipation systems interacting with bosonic canonical bath: Hierarchical equations of motion approach,” Phys. Rev. E 75, 031107 (2007).
  * Jin, Zheng, and Yan (2008) J. Jin, X. Zheng, and Y. Yan, “Exact dynamics of dissipative electronic systems and quantum transport: Hierarchical equations of motion approach,” J. Chem. Phys. 128, 234703 (2008).
  * Yan (2014) Y. Yan, “Theory of open quantum systems with bath of electrons and phonons and spins: Many-dissipaton density matrixes approach,” J. Chem. Phys. 140, 054105 (2014).
  * Schinabeck, Härtle, and Thoss (2018) C. Schinabeck, R. Härtle, and M. Thoss, “Hierarchical quantum master equation approach to electronic-vibrational coupling in nonequilibrium transport through nanosystems: Reservoir formulation and application to vibrational instabilities,” Phys. Rev. B 97, 235429 (2018).
  * Hsieh and Cao (2018) C.-Y. Hsieh and J. Cao, “A unified stochastic formulation of dissipative quantum dynamics. i. generalized hierarchical equations,” J. Chem. Phys. 148, 014103 (2018).
  * Ye  _et al._ (2016) L. Ye, X. Wang, D. Hou, R.-X. Xu, X. Zheng, and Y. Yan, “Heom-quick: a program for accurate, efficient, and universal characterization of strongly correlated quantum impurity systems,” WIREs Comput Mol Sci 6, 608–638 (2016).
  * Tanimura (2020) Y. Tanimura, “Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (heom),” J. Chem. Phys. 153, 020901 (2020).
  * Ishizaki and Tanimura (2005) A. Ishizaki and Y. Tanimura, “Quantum dynamics of system strongly coupled to low-temperature colored noise bath: Reduced hierarchy equations approach,” J. Phys. Soc. Jpn. 74, 3131–3134 (2005).
  * Shi  _et al._ (2009) Q. Shi, L. Chen, G. Nan, R.-X. Xu, and Y. Yan, “Efficient hierarchical liouville space propagator to quantum dissipative dynamics,” J. Chem. Phys. 130, 084105 (2009).
  * Hu, Xu, and Yan (2010) J. Hu, R.-X. Xu, and Y. Yan, “Communication: Padé spectrum decomposition of fermi function and bose function,” J. Chem. Phys. 133, 101106 (2010).
  * Strümpfer and Schulten (2012) J. Strümpfer and K. Schulten, “Open quantum dynamics calculations with the hierarchy equations of motion on parallel computers,” J. Chem. Theory Comput. 8, 2808–2816 (2012).
  * Tsuchimoto and Tanimura (2015) M. Tsuchimoto and Y. Tanimura, “Spins dynamics in a dissipative environment: Hierarchal equations of motion approach using a graphics processing unit (gpu),” J. Chem. Theory Comput. 11, 3859–3865 (2015).
  * Kramer  _et al._ (2018) T. Kramer, M. Noack, A. Reinefeld, M. Rodríguez, and Y. Zelinskyy, “Efficient calculation of open quantum system dynamics and time-resolved spectroscopy with distributed memory heom (dm-heom),” J. Comput. Chem. 39, 1779–1794 (2018).
  * Chen  _et al._ (2022) Z.-H. Chen, Y. Wang, X. Zheng, R.-X. Xu, and Y. Yan, “Universal time-domain prony fitting decomposition for optimized hierarchical quantum master equations,” J. Chem. Phys. 156, 221102 (2022).
  * Xu  _et al._ (2022a) M. Xu, Y. Yan, Q. Shi, J. Ankerhold, and J. Stockburger, “Taming quantum noise for efficient low temperature simulations of open quantum systems,” Phys. Rev. Lett. 129, 230601 (2022a).
  * Ikeda and Nakayama (2022) T. Ikeda and A. Nakayama, “Collective bath coordinate mapping of “hierarchy” in hierarchical equations of motion,” J. Chem. Phys. 156, 104104 (2022).
  * Schollwöck (2011) U. Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Ann. Phys. (NY) 326, 96–192 (2011).
  * Oseledets (2011) I. V. Oseledets, “Tensor-train decomposition,” SIAM J. Sci. Comput. 33, 2295–2317 (2011).
  * Shi  _et al._ (2018) Q. Shi, Y. Xu, Y. Yan, and M. Xu, “Efficient propagation of the hierarchical equations of motion using the matrix product state method,” J. Chem. Phys. 148, 174102 (2018).
  * Borrelli (2019) R. Borrelli, “Density matrix dynamics in twin-formulation: An efficient methodology based on tensor-train representation of reduced equations of motion,” J. Chem. Phys. 150, 234102 (2019).
  * Yan, Xing, and Shi (2020) Y. Yan, T. Xing, and Q. Shi, “A new method to improve the numerical stability of the hierarchical equations of motion for discrete harmonic oscillator modes,” J. Chem. Phys. 153, 204109 (2020).
  * Ke, Borrelli, and Thoss (2022) Y. Ke, R. Borrelli, and M. Thoss, “Hierarchical equations of motion approach to hybrid fermionic and bosonic environments: Matrix product state formulation in twin space,” J. Chem. Phys. 156, 194102 (2022).
  * Vidal (2004) G. Vidal, “Efficient simulation of one-dimensional quantum many-body systems,” Phys. Rev. Lett. 93, 040502 (2004).
  * Shi, Duan, and Vidal (2006) Y.-Y. Shi, L.-M. Duan, and G. Vidal, “Classical simulation of quantum many-body systems with a tree tensor network,” Phys. Rev. A 74, 022320 (2006).
  * Tagliacozzo, Evenbly, and Vidal (2009) L. Tagliacozzo, G. Evenbly, and G. Vidal, “Simulation of two-dimensional quantum systems using a tree tensor network that exploits the entropic area law,” Phys. Rev. B 80, 235127 (2009).
  * Murg  _et al._ (2010) V. Murg, F. Verstraete, Ö. Legeza, and R. M. Noack, “Simulating strongly correlated quantum systems with tree tensor networks,” Phys. Rev. B 82, 205105 (2010).
  * Li, von Delft, and Xiang (2012) W. Li, J. von Delft, and T. Xiang, “Efficient simulation of infinite tree tensor network states on the bethe lattice,” Phys. Rev. B 86, 195137 (2012).
  * Changlani  _et al._ (2013) H. J. Changlani, S. Ghosh, C. L. Henley, and A. M. Läuchli, “Heisenberg antiferromagnet on cayley trees: Low-energy spectrum and even/odd site imbalance,” Phys. Rev. B 87, 085107 (2013).
  * Nakatani and Chan (2013) N. Nakatani and G. K.-L. Chan, “Efficient tree tensor network states (ttns) for quantum chemistry: Generalizations of the density matrix renormalization group algorithm,” J. Chem. Phys. 138, 134113 (2013).
  * Murg  _et al._ (2015) V. Murg, F. Verstraete, R. Schneider, P. R. Nagy, and O. Legeza, “Tree tensor network state with variable tensor order: An efficient multireference method for strongly correlated systems,” J. Chem. Theory Comput. 11, 1027–1036 (2015).
  * Gunst  _et al._ (2018) K. Gunst, F. Verstraete, S. Wouters, O. Legeza, and D. Van Neck, “T3ns: Three-legged tree tensor network states,” J. Chem. Theory Comput. 14, 2026–2033 (2018).
  * Schröder  _et al._ (2019) F. A. Schröder, D. H. Turban, A. J. Musser, N. D. Hine, and A. W. Chin, “Tensor network simulation of multi-environmental open quantum dynamics via machine learning and entanglement renormalisation,” Nat. Commun. 10, 1–10 (2019).
  * Larsson (2019) H. R. Larsson, “Computing vibrational eigenstates with tree tensor network states (ttns),” J. Chem. Phys. 151, 204102 (2019).
  * Ferrari, Magnifico, and Montangero (2022) G. Ferrari, G. Magnifico, and S. Montangero, “Adaptive-weighted tree tensor networks for disordered quantum many-body systems,” Phys. Rev. B 105, 214201 (2022).
  * Seitz  _et al._ (2022) P. Seitz, I. Medina, E. Cruz, Q. Huang, and C. B. Mendl, “Simulating quantum circuits using tree tensor networks,” arXiv preprint arXiv:2206.01000 (2022).
  * Milsted  _et al._ (2019) A. Milsted, M. Ganahl, S. Leichenauer, J. Hidary, and G. Vidal, “Tensornetwork on tensorflow: A spin chain application using tree tensor networks,” arXiv preprint arXiv:1905.01331 (2019).
  * Bauernfeind and Aichhorn (2020) D. Bauernfeind and M. Aichhorn, “Time dependent variational principle for tree tensor networks,” SciPost Physics 8, 024 (2020).
  * Kloss, Reichman, and Bar Lev (2020) B. Kloss, D. Reichman, and Y. Bar Lev, “Studying dynamics in two-dimensional quantum lattices using tree tensor network states,” SciPost Physics 9, 070 (2020).
  * Ceruti, Lubich, and Walach (2021) G. Ceruti, C. Lubich, and H. Walach, “Time integration of tree tensor networks,” SIAM J. Numer. Anal. 59, 289–313 (2021).
  * Hu  _et al._ (2011) J. Hu, M. Luo, F. Jiang, R.-X. Xu, and Y. Yan, “Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems,” J. Chem. Phys. 134, 244106 (2011).
  * Schmutz (1978) M. Schmutz, “Real-time green’s functions in many body problems,” Zeitschrift für Physik B Condensed Matter 30, 97–106 (1978).
  * Suzuki (1991) M. Suzuki, “Density matrix formalism, double-space and thermo field dynamics in non-equilibrium dissipative systems,” Internat. J. Modern Phys. B 5, 1821–1842 (1991).
  * Arimitsu and Umezawa (1987) T. Arimitsu and H. Umezawa, “Non-equilibrium thermo field dynamics,” Prog. Theor. Phys. 77, 32–52 (1987).
  * Feiguin and White (2005) A. E. Feiguin and S. R. White, “Finite-temperature density matrix renormalization using an enlarged hilbert space,” Phys. Rev. B 72, 220401 (2005).
  * Verstraete, Garcia-Ripoll, and Cirac (2004) F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, “Matrix product density operators: Simulation of finite-temperature and dissipative systems,” Phys. Rev. Lett. 93, 207204 (2004).
  * Borrelli and Gelin (2021) R. Borrelli and M. F. Gelin, “Finite temperature quantum dynamics of complex systems: Integrating thermo-field theories and tensor-train methods,” WIREs Comput Mol Sci , e1539 (2021).
  * Cirac and Verstraete (2009) J. I. Cirac and F. Verstraete, “Renormalization and tensor product states in spin chains and lattices,” J. Phys. A Math. Theor. 42, 504004 (2009).
  * Chan  _et al._ (2016) G. K.-L. Chan, A. Keselman, N. Nakatani, Z. Li, and S. R. White, “Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms,” J. Chem. Phys. 145, 014102 (2016).
  * Jaschke, Wall, and Carr (2018) D. Jaschke, M. L. Wall, and L. D. Carr, “Open source matrix product states: Opening ways to simulate entangled many-body quantum systems in one dimension,” Comput. Phys. Commun. 225, 59–91 (2018).
  * Lubich  _et al._ (2013) C. Lubich, T. Rohwedder, R. Schneider, and B. Vandereycken, “Dynamical approximation by hierarchical tucker and tensor-train tensors,” SIAM J. Matrix Anal. Appl. 34, 470–494 (2013).
  * Beck and Meyer (1997) M. H. Beck and H.-D. Meyer, “An efficient and robust integration scheme for the equations of motion of the multiconfiguration time-dependent hartree (mctdh) method,” Zeitschrift für Physik D Atoms, Molecules and Clusters 42, 113–129 (1997).
  * Wang and Thoss (2003) H. Wang and M. Thoss, “Multilayer formulation of the multiconfiguration time-dependent hartree theory,” J. Chem. Phys. 119, 1289–1299 (2003).
  * Meyer, Gatti, and Worth (2009) H.-D. Meyer, F. Gatti, and G. A. Worth,  _Multidimensional quantum dynamics: MCTDH theory and applications_ (John Wiley & Sons, 2009).
  * Manthe (2015) U. Manthe, “The multi-configurational time-dependent hartree approach revisited,” J. Chem. Phys. 142, 244109 (2015).
  * Wang (2015) H. Wang, “Multilayer multiconfiguration time-dependent hartree theory,” J. Phys. Chem. A 119, 7951–7965 (2015).
  * Lindoy, Kloss, and Reichman (2021) L. P. Lindoy, B. Kloss, and D. R. Reichman, “Time evolution of ml-mctdh wavefunctions. i. gauge conditions, basis functions, and singularities,” J. Chem. Phys. 155, 174108 (2021).
  * Montangero, Rico, and Silvi (2022) S. Montangero, E. Rico, and P. Silvi, “Loop-free tensor networks for high-energy physics,” Philos. Trans. R. Soc. A 380, 20210065 (2022).
  * Yan  _et al._ (2021) Y. Yan, M. Xu, T. Li, and Q. Shi, “Efficient propagation of the hierarchical equations of motion using the tucker and hierarchical tucker tensors,” J. Chem. Phys. 154, 194104 (2021).
  * Haegeman, Osborne, and Verstraete (2013) J. Haegeman, T. J. Osborne, and F. Verstraete, “Post-matrix product state methods: To tangent space and beyond,” Phys. Rev. B 88, 075133 (2013).
  * Lubich, Oseledets, and Vandereycken (2015) C. Lubich, I. V. Oseledets, and B. Vandereycken, “Time integration of tensor trains,” SIAM J. Numer. Anal. 53, 917–941 (2015).
  * Haegeman  _et al._ (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, “Unifying time evolution and optimization with matrix product states,” Phys. Rev. B 94, 165116 (2016).
  * Schröder and Chin (2016) F. A. Schröder and A. W. Chin, “Simulating open quantum dynamics with time-dependent variational matrix product states: Towards microscopic correlation of environment dynamics and reduced system evolution,” Phys. Rev. B 93, 075105 (2016).
  * Baiardi and Reiher (2019) A. Baiardi and M. Reiher, “Large-scale quantum dynamics with matrix product states,” J. Chem. Theory Comput. 15, 3481–3498 (2019).
  * Mendl (2018) C. B. Mendl, “Time evolution of matrix product operators with energy conservation,” arXiv preprint arXiv:1812.11876 (2018).
  * Paeckel  _et al._ (2019) S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, “Time-evolution methods for matrix-product states,” Ann. Phys. (NY) 411, 167998 (2019).
  * Wilkins and Dattani (2015) D. M. Wilkins and N. S. Dattani, “Why quantum coherence is not important in the fenna–matthews–olsen complex,” J. Chem. Theory Comput. 11, 3411–3419 (2015).
  * Fishman, White, and Stoudenmire (2022) M. Fishman, S. R. White, and E. M. Stoudenmire, “The ITensor Software Library for Tensor Network Calculations,” SciPost Phys. Codebases , 4 (2022).
  * Yang and White (2020) M. Yang and S. R. White, “Time-dependent variational principle with ancillary krylov subspace,” Phys. Rev. B 102, 094315 (2020).
  * Dunnett and Chin (2021) A. J. Dunnett and A. W. Chin, “Efficient bond-adaptive approach for finite-temperature open quantum dynamics using the one-site time-dependent variational principle for matrix product states,” Phys. Rev. B 104, 214302 (2021).
  * Borrelli and Dolgov (2021) R. Borrelli and S. Dolgov, “Expanding the range of hierarchical equations of motion by tensor-train implementation,” J. Phys. Chem. B 125, 5397–5407 (2021).
  * Li, Gleis, and Von Delft (2022) J.-W. Li, A. Gleis, and J. Von Delft, “Time-dependent variational principle with controlled bond expansion for matrix product states,” arXiv preprint arXiv:2208.10972 (2022).
  * Xu  _et al._ (2022b) Y. Xu, Z. Xie, X. Xie, U. Schollwöck, and H. Ma, “Stochastic adaptive single-site time-dependent variational principle,” JACS Au 2, 335–340 (2022b).
  * Ceruti, Lubich, and Sulz (2023) G. Ceruti, C. Lubich, and D. Sulz, “Rank-adaptive time integration of tree tensor networks,” SIAM J. Numer. Anal. 61, 194–222 (2023).
  * Evers  _et al._ (2022) F. Evers, A. Aharony, N. Bar-Gill, O. Entin-Wohlman, P. Hedegård, O. Hod, P. Jelinek, G. Kamieniarz, M. Lemeshko, K. Michaeli, _et al._ , “Theory of chirality induced spin selectivity: Progress and challenges,” Adv. Mater. 34, 2106629 (2022).