# Quantum dynamics of coupled excitons and phonons in chain-like systems:
tensor train approaches and higher-order propagators

Patrick Gelß  Institut für Mathematik, Freie Universität Berlin
Arnimallee 3–9, D-14195 Berlin, Germany  Zuse-Institut Berlin, Takustraße 7,
D-14195 Berlin, Germany  Sebastian Matera  Institut für Mathematik, Freie
Universität Berlin
Arnimallee 3–9, D-14195 Berlin, Germany  Fritz-Haber-Institut der Max-Planck-
Gesellschaft, Faradayweg 4–6, D-14195 Berlin, Germany  Rupert Klein  Institut
für Mathematik, Freie Universität Berlin
Arnimallee 3–9, D-14195 Berlin, Germany  Burkhard Schmidt
<EMAIL_ADDRESS>Institut für Mathematik, Freie Universität
Berlin
Arnimallee 3–9, D-14195 Berlin, Germany  Weierstraß-Institut für Angewandte
Analysis und Stochastik, Mohrenstraße 39, D-10117 Berlin, Germany

###### Abstract

We investigate the use of tensor-train approaches to the solution of the time-
dependent Schrödinger equation for chain-like quantum systems with on-site and
nearest-neighbor interactions only. Using the efficient SLIM representation
for low-rank tensor train representations of quantum-mechanical Hamiltonians,
we aim at reducing the memory consumption as well as the computation costs, in
order to mitigate the curse of dimensionality as much as possible. As an
example, coupled excitons and phonons modeled in terms of Fröhlich-Holstein
type Hamiltonians are studied here. By comparing our tensor-train based
results with semi-analytical results, we demonstrate the key role of the ranks
of tensor-train representations for quantum state vectors. Both the
computational effort of the propagations and the accuracy that can be reached
crucially depend on the maximum number of ranks chosen. Typically, an
excellent quality of the solutions is found only when the ranks exceeds a
certain value. That threshold, however, is very different for excitons,
phonons, and coupled systems. One class of propagation schemes used in the
present work builds on splitting the Hamiltonian into two groups of
interleaved nearest-neighbor interactions which commutate within each of the
groups. In addition to the first order Lie-Trotter and the second order
Strang-Marchuk splitting schemes, we have also implemented the 4-th order
Yoshida-Neri and the 8-th order Kahan-Li symplectic compositions. Especially
the latter two are demonstrated to yield very accurate results, close to
machine precision. However, due to the computational costs, currently their
use is restricted to rather short chains. Another class of propagators
involves explicit, time-symmetrized Euler integrators. Building on the
traditional second order differencing method, we have also implemented higher
order methods. Especially the 4-th order variant is recommended for quantum
simulations of longer chains, even though the high precision of the splitting
schemes cannot be reached. Moreover, the scaling of the computational effort
with the dimensions of the local Hilbert spaces is much more favorable for the
differencing than for the splitting schemes.

##  I Introduction

Excitonic energy transport (EET), i.e., the transport of electronic energy
without net charge is of paramount importance in many biological and
technological systems. In particular, the efficiency of the conversion of
light into chemical or electric energy crucially depends on it [1, 2]. For
example, in photovoltaic devices the transport of energy from the region of
primary excitation to the charge separating interfaces has to be fast enough
to compete with radiative and non-radiative decay channels, most prominently
the recombination of electrons and holes. The efficiency of EET is strongly
affected by the coupling of electronic excitations to (intramolecular and/or
intermolecular) vibrational degrees of freedom [3, 4]. The simplest
theoretical approaches to this exciton-phonon coupling (EPC) are based on the
assumption of chain-like, quasi-1D arrangements. Examples relevant for
semiconducting materials are lattice vibrations in 1D molecular crystal models
or various stretching, bending, and/or torsional motions in conjugated polymer
chains. The ”dressing” of excitons by vibrations in these classes of systems
is often described in terms of Holstein or Fröhlich type Hamiltonians.
Originally developed for polarons [5], these model Hamiltonians, or extensions
thereof, are routinely used in a variety of applications, ranging from
vibrational energy transport in helical proteins [6, 7, 8] up to electronic
energy transport in conjugated polymer chains [9, 10]. Despite of the apparent
simplicity of the Fröhlich-Holstein type systems, their quantum dynamics still
poses many open questions, as noted in the abstract of Ref. [5]: ”It is
remarkable how the Fröhlich polaron … has resisted full analytical or
numerical solution … since $\sim$1950, when its Hamiltonian was first written.
The field has been a testing ground for analytical, semi-analytical and
numerical techniques, such as path integrals, strong-coupling perturbation
expansion, advanced variational, exact diagonalization (ED) and quantum Monte
Carlo (QMC) techniques.” Hence, the study of Fröhlich-Holstein type systems as
models for exciton dynamics in the presence of EPC remains an intriguing and
timely research field.

Essentially, the problems of theoretical studies of EPC are due to the quantum
dynamical treatment required for the above-mentioned models, i.e., solutions
of time-dependent Schrödinger equations (TDSE) in many dimensions. Quantum
dynamical simulations to solve the TDSE with conventional methods eventually
become unfeasible when increasing the number of dimensions [11]. One way to
circumvent this so-called curse of dimensionality is to resort to mixed (or
hybrid) quantum-classical methods. There, only the electronic degrees of
freedom are treated quantum-mechanically, while the nuclear (ionic) degrees of
freedom are treated classically, typically justified by the small electron to
ion mass ratio [12]. The most basic of these hybrid methods is the quantum-
classical Ehrenfest approach which plays an important role, e.g., in the
Davydov soliton theory for EET in helical proteins [6, 7, 8]. Apart from the
classical approximation, the main limitation of Ehrenfest dynamics is that it
builds on a mean-field coupling of the two sub-systems [13, 14, 15]. This
bottleneck is overcome in surface hopping trajectory approaches, where non-
adiabatic transitions are approximately treated by stochastic hopping of
trajectories between different (diabatic or adiabatic) energy states [16].
Even though surface hopping trajectory techniques have lately been used a lot
in studying the photo-induced physics and chemistry of EET [17, 18, 19], the
validity of the underlying approximations of quantum-classical dynamics is not
always clear. In particular, for EET in organic semiconductors, there is
sometimes no clear distinction between slow and fast degrees of freedom [4].

The focus of the present work is on methods for fully quantum-mechanical
approaches to treat coupled excitons and phonons. In the literature, reports
of such simulations are scarce. One example are quantum dynamical
investigations of the related problem of formation of polarons in crystals
where variational techniques for representing quantum state vectors were used
[20, 21]. In more recent work, structure and dynamics of excitons in
$\pi$-conjugated polymers were investigated, with emphasis on coupling to the
torsional dynamics of the chain [9, 10]. In those studies, the multi-
configurational time-dependent Hartree (MCTDH) method and its multi-layer (ML)
extension were used, which are known as highly efficient algorithms for
quantum dynamics in high dimensions [22, 23].

In recent years, the interest in low-rank tensor decompositions has been
growing steadily. The potential of different tensor formats like the canonical
format [24], the Tucker format [25], and particularly the tensor-train (TT,
aka matrix product states, MPS) format [26] for mitigating the curse of
dimensionality has been shown in various application areas. So far, tensor
decomposition techniques have been used in the fields of, e.g., quantum
mechanics [27, 28, 29, 30], neuro-imaging [31, 32, 33], system identification
[34, 35, 36], and supervised learning [37, 38, 39], just to name a few. In
previous work, we proposed TT-based techniques for numerically solving the
time-independent Schrödinger equations (TISE) for coupled excitons and phonons
[40]. There, it has been demonstrated that high-dimensional quantum systems
comprising only NN interactions may be simulated and analyzed without an
exponential scaling of the memory consumption and/or computational effort by
using so-called SLIM decompositions [41]. Note that the underlying TT
representations are a special case of the hierarchical Tucker format [42],
which is equivalent to the wave function representation in the multi-layer
(ML) variant of MCTDH. In recent comparisons of the efficiency of both methods
for calculating ground states of chains of molecular rotors, it has been
found, that TT/MPS techniques are less time- and memory-consuming than the
state-of-the-art implementation of ML-MCTDH [43, 44].

In the present work, we want to exploit low-rank TT decompositions of quantum
state vectors and Hamiltonian matrices for the case of the time-dependent
Schrödinger equation (TDSE). In addition to previously published time
integrators for TTs using projectors onto the tangent space of the tensor
manifold [45, 46, 47], also symmetric Euler differencing schemes [48] and
split operator techniques [49, 50] have been suggested recently for use with
TT representations [27, 51, 52]. Being explicit and easy to implement, these
two integrators have been in use since the early days of numerical quantum
dynamics [53]. Our goal is to investigate them along with their higher-order
variants, see Sec. III. There, the various time integrators are subjected to
rigorous tests for realistic models of excitons, phonons, and exciton-phonon-
coupling. To that end, we are not only considering their ability to conserve
norm and energy but also how accurately they can reproduce semi-analytic
propagations which are available for special cases, see Sec. IV. On the one
hand, our focus will be on the scaling of the numerical effort with the chain
length and the Hilbert space dimensions, thus assessing the propagators for
use in simulations of large systems. On the other hand, we will elucidate on
the use of higher-order methods in comparison to lower-order ones in the
context of low-rank TT representations of high-dimensional wave functions. In
particular, we will address the question whether the increased costs for each
propagation step will be compensated by the possibility to use larger step
sizes, and, hence, a smaller number of steps, while obtaining solutions with
the same accuracy. Note that implicit schemes such as the trapezoidal rule or
the midpoint rule are not considered here. Even though they offer the
advantage of conserving energy, their computational effort is unfavorable
because they involve solutions of large-scale linear systems of equations.
While this is in principle possible for TTs by virtue of the alternating
linear scheme [54], which is also used in our previous TISE study [40], it is
computationally too expensive to do this in every time step.

##  II Model Hamiltonians

In the present work, we use Fröhlich-Holstein type Hamiltonians to describe
the dynamics of excitons and phonons. With the coupling of the two subsystems
assumed to be linear, these can be considered as the simplest models of EPC in
organic semiconductors [3, 4]. Originally developed for electrons in a
polarizable lattice [5], these models, partly with minor modifications, are
also known as Huang-Rhys, or Peierls models. There is also a close analogy to
Davydov models describing the interaction of the amide I vibrations and the
hydrogen bond stretching in $\alpha$-helices of proteins [7, 8].

Here, we restrict ourselves to the consideration of effectively one-
dimensional, linear chains comprising $N$ exciton-supporting sites assuming NN
interactions only, see also our previous work [40]. Such models are suitable
for the description of two classes of systems [1, 55] in the context of
photovoltaics applications. First, in conjugated polymer chains without major
kinks or turns, the excitonic sites interact with each other mainly via a
Dexter mechanism (”through bond”) [56]. Second, in crystals of polycyclic
aromatic molecules, the coupling of excitons occurs mainly in the direction
perpendicular to the planes of the stacked molecules, governed by the Förster
mechanism (”through space”) [57].

The total Hamiltonians for the coupled excitons and phonons can be written as

$H=H^{\mathrm{(ex)}}\otimes\mathbb{I}^{\mathrm{(ph)}}+\mathbb{I}^{\mathrm{(ex)}}\otimes
H^{\mathrm{(ph)}}+H^{\mathrm{(epc)}}$ (1)

with the superscripts (ex), (ph), and (epc) standing for excitons, phonons and
their coupling, respectively, while the $\mathbb{I}$ are the respective
identity operators. Note that atomic units with $\hbar=1$ are used throughout
this work.

Restricting ourselves to a chain of two-state-systems, the excitonic
Hamiltonian for a linear or cyclic system of $N$ (identical or not) sites is
given in terms of (bosonic) exciton raising, $b_{i}^{\dagger}$, and lowering,
$b_{i}$, operators

$H^{\mathrm{(ex)}}=\sum_{i=1}^{N}\alpha_{i}b_{i}^{\dagger}b_{i}+\sum_{i=1}^{N}\beta_{i}\left(b_{i}^{\dagger}b_{i+1}+b_{i}b_{i+1}^{\dagger}\right)$
(2)

with local (”on site”) excitation energies $\alpha_{i}$ for site $i$. The
nearest-neighbor (NN) coupling energies $\beta_{i}$ between site $i$ and
$i+1$, also known as transfer or hopping integrals, determine the exciton
delocalization and mobility. Here and throughout the following, the last
summand ($i=N$) of the NN coupling term (with indices $i+1$ replaced by 1) is
intended for systems with periodic boundary conditions only and is deleted
otherwise.

The vibrational (phononic) Hamiltonian for a one-dimensional lattice can be
written in terms of masses $m_{i}$, displacement coordinates $R_{i}$, and
conjugated momenta $P_{i}$ in harmonic approximation

$H^{\mathrm{(ph)}}=\frac{1}{2}\sum_{i=1}^{N}\frac{P_{i}^{2}}{m_{i}}+\frac{1}{2}\sum_{i=1}^{N}m_{i}\nu_{i}^{2}R_{i}^{2}+\frac{1}{2}\sum_{i=1}^{N}\mu_{i}\omega_{i}^{2}\left(R_{i}-R_{i+1}\right)^{2}$
(3)

where each site $i$ is restrained to oscillate around its equilibrium position
by a harmonic potential with frequency $\nu_{i}$ and where the NN interactions
between sites $i$ and $i+1$ are modeled by harmonic oscillators with frequency
$\omega_{i}$ and corresponding reduced masses
$\mu_{i}=m_{i}m_{i+1}/(m_{i}+m_{i+1})$. Note that the frequencies $\nu_{i}$
and $\omega_{i}$ are assumed to be independent of the excitonic state of the
system.

In parallel with the treatment of the excitons, we introduce second
quantization also for the phononic Hamiltonian of Eq. (3) resulting in

$H^{\mathrm{(ph)}}=\sum_{i=1}^{N}\tilde{\nu}_{i}\left(c_{i}^{\dagger}c_{i}+\frac{1}{2}\right)-\sum_{i=1}^{N}\tilde{\omega}_{i}\left(c_{i}^{\dagger}+c_{i}\right)\left(c_{i+1}^{\dagger}+c_{i+1}\right)$
(4)

where raising and lowering operators of (local) vibrations of site $i$ are
indicated by $c_{i}^{\dagger}$ and $c_{i}$, respectively, and where for the
effective frequencies of single site and NN pair vibrations the following
expressions

$\displaystyle\tilde{\nu}_{i}$ $\displaystyle=$
$\displaystyle\sqrt{\nu_{i}^{2}+\frac{m_{i-1}}{m_{i}+m_{i-1}}\omega_{i-1}^{2}+\frac{m_{i+1}}{m_{i}+m_{i+1}}\omega_{i}^{2}}$
(5) $\displaystyle\tilde{\omega}_{i}$ $\displaystyle=$
$\displaystyle\frac{\mu_{i}\omega_{i}^{2}}{2\sqrt{m_{i}\tilde{\nu}_{i}m_{i+1}\tilde{\nu}_{i+1}}}$
(6)

have been found where for linear systems without periodic boundary conditions
the second or third term under the square root of Eq. (5) is omitted for the
first ($i=1$) or last ($i=N$) site, respectively.

Linear models for the coupling of excitons and phonons can be formulated in
terms of a Fröhlich-Holstein type Hamiltonian [5]

$H^{\mathrm{(epc)}}=\sum_{i=1}^{N}\sigma_{i}b_{i}^{\dagger}b_{i}\otimes\left[R_{i+1}-R_{i-1}\right]=\sum_{i=1}^{N}b_{i}^{\dagger}b_{i}\otimes\left[\bar{\sigma}_{i}\left(c_{i+1}^{\dagger}+c_{i+1}\right)-\bar{\bar{\sigma}}_{i}\left(c_{i-1}^{\dagger}+c_{i-1}\right)\right]$
(7)

where the coupling constants
$\sigma_{i}=\mathrm{d}\alpha_{i}/\mathrm{d}(R_{i+1}-R_{i-1})$ give the
linearized dependence of the exciton energies $\alpha_{i}$ on the sum of the
NN distances to the left and to the right of each site. The bar notation is
used here to convert the EPC constants to second quantization,
$\bar{\sigma}_{i}=\sigma_{i}/\sqrt{2m_{i+1}\tilde{\nu}_{i+1}}$ and
$\bar{\bar{\sigma}}_{i}=\sigma_{i}/\sqrt{2m_{i-1}\tilde{\nu}_{i-1}}$. Note
that in our previous work [40], the distinction between EPC constants with
bars and double bars was missing, which, however, is not required for the
homogeneous systems mainly investigated there.

This symmetric formulation of the EPC is suitable in case of a mirror symmetry
of the individual sites frequently occurring in semiconductor materials
whereas one-sided expressions are typically used for excitons in protein
helices where an amide I vibration couples much more strongly to the directly
adjacent H-bond [58, 7, 8]. For other model assumptions, including Peierls
type models based on an $R$–dependence of the NN coupling parameters
$\beta_{i}$, see our previous work [40].

##  III Tensor train decompositions

###  III.1 Statement of the problem

Quantum dynamics of non-relativistic, closed systems (i.e., without
dissipation or dephasing) is given in terms of the time-dependent Schrödinger
equation (TDSE)

$\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}t}\Psi(t)=H\Psi(t),$ (8)

where $t$ is the time and where we restrict ourselves to time-independent
Hamiltonians $H$, in our case the Hamiltonian given in Eq. (1). Note that
atomic units with $\hbar=1$ are used throughout this work. The initial quantum
state is given by $\Psi(t=0)=\Psi_{0}$, which in typical photophysical or
photochemical simulation scenarios describes a state prepared by the
interaction with light. Let us assume that a coupled excitonic-phononic state
vector for each of the constituent sites will be given by a vector in the
product Hilbert space

$\psi_{i}\in\mathcal{H}_{i}:=\mathcal{H}_{i}^{\mathrm{(ex)}}\otimes\mathcal{H}_{i}^{\mathrm{(ph)}}\cong\mathbb{C}^{d_{i}}$
(9)

with dimension $d_{i}:=d_{i}^{\mathrm{(ex)}}d_{i}^{\mathrm{(ph)}}$, where
$d_{i}^{\mathrm{(ex)}}$ and $d_{i}^{\mathrm{(ph)}}$ are the dimensions of the
Hilbert spaces of excitonic and vibrational state vectors, respectively. Then,
time-dependent quantum state vectors $\Psi(t)$ of a complete chain (or ring)
comprising $N$ sites can be understood as tensors of order $N$,

$\Psi\in\mathcal{H}:=\mathcal{H}_{1}\otimes\mathcal{H}_{2}\otimes\cdots\otimes\mathcal{H}_{N}\cong\mathbb{C}^{D}$
(10)

with a total dimension of $D:=\prod_{i=1}^{N}d_{i}$, i.e., the numerical
discretization of $\Psi$ is a tensor in $\mathbb{C}^{d_{1}\times\dots\times
d_{N}}$. In analogy, the tensorized representation $H$ of the (real-valued,
symmetric) Hamiltonian given in Eq. (1) is a tensor in
$\mathbb{R}^{(d_{1}\times d_{1})\times\dots\times(d_{N}\times d_{N})}$.

Because the formal time evolution of quantum states

$\Psi(t)=\exp(-\mathrm{i}tH)\Psi_{0}$ (11)

is not computationally available for large $N$, our aim is to numerically
integrate the TDSE (8) in terms of tensor representations by employing time
integrators. In general, however, the storage consumption of (non-sparse)
tensors of the form (1) or (10) grows exponentially with the order $N$. Due to
this so-called curse of dimensionality, storing the considered tensors in full
format may be not feasible for large $N$. Therefore, special tensor
representations are required. As in [40], we will focus on the tensor train
format (TT format) [59, 60]. This type of tensor decomposition is not only
tailor-made for chain-like topologies such as our Hamiltonians with NN
interactions, but it is also a preferred choice for representing high-
dimensional tensors in terms of storage consumption and computational
robustness. We will also introduce properly designed numerical methods based
on tensor trains which we will use to solve evolution equations such as the
TDSE (8) in high-dimensional tensor spaces which otherwise would be extremely
expensive or even impossible with today’s computational resources.

###  III.2 Tensor trains and SLIM decomposition

In this work, the TT format is used to represent a quantum state $\Psi$. That
is, we decompose a tensor in $\mathbb{C}^{d_{1}\times\dots\times d_{N}}$ into
$N$ component tensors which are coupled in a chain, see Fig. 1 and Def. 1.

$r_{1}$$r_{2}$$r_{3}$$r_{4}$$d_{1}$$d_{2}$$d_{3}$$d_{4}$$d_{5}$ (a)

$r_{1}$$r_{2}$$r_{3}$$r_{4}$$d_{1}$$d_{2}$$d_{3}$$d_{4}$$d_{5}$$d_{1}$$d_{2}$$d_{3}$$d_{4}$$d_{5}$
(b)

Figure 1: Graphical representation of tensor trains: Tensors are depicted as
circles with different arms indicating the set of free indices. (a) Tensor of
order 5 in TT format, the first and the last core are matrices, the other
cores are tensors of order 3. (b) Linear operator in TT format, the first and
the last core are tensors of order 3, the other cores are tensors of order 4.
In both cases, arms corresponding to first and last TT ranks are omitted since
$r_{0}=r_{5}=1$.

###### Definition 1

A tensor $\Psi\in\mathbb{C}^{d_{1}\times d_{2}\times\dots\times d_{N}}$ is
said to be in the TT format if

$\Psi=\sum_{k_{0}=1}^{r_{0}}\cdots\sum_{k_{N}=1}^{r_{N}}\Psi^{(1)}_{k_{0},:,k_{1}}\otimes\cdots\otimes\Psi^{(N)}_{k_{N-1},:,k_{N}}.$
(12)

The tensors $\Psi^{(i)}\in\mathbb{C}^{r_{i-1}\times d_{i}\times r_{i}}$ of
order $3$ are called TT cores and the numbers $r_{i}$ are called TT ranks. It
holds that $r_{0}=r_{N}=1$ and $r_{i}\geq 1$ for $i=1,\dots,N-1$.

To reduce the memory consumption as well as the computational costs, we are
especially interested in low-rank TT decompositions. In Ref. [40], we used the
TT format to represent stationary quantum states, i.e., solutions of the TISE.
Here, we are interested in the numerical solution of Eq. (8) in form of a time
series $(\Psi(t_{0}),\Psi(t-1),\dots,\Psi(T))$ of TT representations. That is,
if the quantum states at the considered time points can be approximated
accurately enough by tensor trains with ranks of manageable size, we are able
to mitigate the curse of dimensionality. For more details about the TT format,
we refer to Refs. [59, 60, 26].

Since the efficiency of the tensor-based integrators described in the
following sections also depends on the ranks of the operator
$H\in\mathbb{R}^{(d_{1}\times d_{1})\times\dots\times(d_{N}\times d_{N})}$, we
want to construct low-rank TT representations of the Hamiltonian tensor as
well. All terms of the Hamiltonians introduced in Sec. II either act locally
on a single site or on two neighboring sites. Thus, as we explained in Ref.
[40], $H$ can be expressed in the TT format by using a so-called SLIM
decomposition [41]:

$\begin{split}H=\left\llbracket\begin{matrix}S_{1}&L_{1}&I_{1}&M_{1}\end{matrix}\right\rrbracket\otimes\left\llbracket\begin{matrix}I_{2}&0&0&0\\\
M_{2}&0&0&0\\\ S_{2}&L_{2}&I_{2}&0\\\
0&0&0&J_{2}\end{matrix}\right\rrbracket\otimes\dots\otimes\left\llbracket\begin{matrix}I_{N-1}&0&0&0\\\
M_{N-1}&0&0&0\\\ S_{N-1}&L_{N-1}&I_{N-1}&0\\\
0&0&0&J_{N-1}\end{matrix}\right\rrbracket\otimes\left\llbracket\begin{matrix}I_{N}\\\
M_{N}\\\ S_{N}\\\ L_{N}\end{matrix}\right\rrbracket.\end{split}$ (13)

Above, we use the core notation of the TT format, cf. Refs. [61, 41, 40]. The
single-site components $S_{i}$ for $i=1,\dots,N$, see Eqs. (2,4), are given by

$S_{i}=\alpha_{i}b_{i}^{\dagger}b_{i}+\tilde{\nu}_{i}\left(c_{i}^{\dagger}c_{i}+\frac{1}{2}\right)$
(14)

The NN interactions are represented by the components $L_{i}$ and $M_{i+1}$
defined as

$\displaystyle\llbracket L_{i}\rrbracket$
$\displaystyle=\left\llbracket\begin{matrix}L_{i,1}&\dots&L_{i,\xi_{i}}\end{matrix}\right\rrbracket$
$\displaystyle\in\mathbb{R}^{1\times d_{i}\times d_{i}\times\xi_{i}},$
$\displaystyle\llbracket M_{i+1}\rrbracket$
$\displaystyle=\left\llbracket\begin{matrix}M_{i+1,1}&\dots&M_{i+1,\xi_{i}}\end{matrix}\right\rrbracket^{\mathbb{T}}$
$\displaystyle\in\mathbb{R}^{\xi_{i}\times d_{i+1}\times d_{i+1}\times 1},$

for $i=1,\dots,N-1$, and

$\displaystyle\llbracket L_{N}\rrbracket$
$\displaystyle=\left\llbracket\begin{matrix}L_{N,1}&\dots&L_{N,\xi_{N}}\end{matrix}\right\rrbracket^{\mathbb{T}}$
$\displaystyle\in\mathbb{R}^{\xi_{N}\times d_{N}\times d_{N}\times 1},$ (15)
$\displaystyle\llbracket M_{1}\rrbracket$
$\displaystyle=\left\llbracket\begin{matrix}M_{1,1}&\dots&M_{1,\xi_{d}}\end{matrix}\right\rrbracket$
$\displaystyle\in\mathbb{R}^{1\times d_{1}\times d_{1}\times\xi_{N}},$

for the case of cyclic systems only.

The two-site contributions to the total Hamiltonian, see Eqs. (2,4,7), are
expressible in terms of five operators $L_{i,\lambda}$ acting on site $i$ and
five corresponding operators $M_{i+1,\lambda}$ acting on the neighboring site
$i+1$:

$L_{i,1}=\beta_{i}b_{i}^{\dagger}$, |  |  $M_{i+1,1}=b_{i+1}$,
---|---|---
$L_{i,2}=\beta_{i}b_{i}$, |  |  $M_{i+1,2}=b_{i+1}^{\dagger}$,
$L_{i,3}=-\tilde{\omega}_{i}\left(c_{i}^{\dagger}+c_{i}\right)$, |  |  $M_{i+1,3}=c_{i+1}^{\dagger}+c_{i+1}$,
$L_{i,4}=\bar{\sigma}_{i}b_{i}^{\dagger}b_{i}$, |  |  $M_{i+1,4}=c_{i+1}^{\dagger}+c_{i+1}$,
$L_{i,5}=-\left(c_{i}^{\dagger}+c_{i}\right)$, |  |  $M_{i+1,5}=\bar{\bar{\sigma}}_{i+1}b_{i+1}^{\dagger}b_{i+1}$.
(16)

If certain coefficients vanish, e.g., if we consider purely excitonic or
phononic Hamiltonians, the numbers $\xi_{i}$ of two-site interactions may even
be smaller, leading to reduced TT ranks of the Hamiltonian. A detailed
derivation of the above components was given in Ref. [40]. Besides the fact
that the storage consumption scales only linearly with the number of sites if
we consider systems with a fixed/bounded number of NN contributions, the
repeating pattern of the SLIM decomposition allows us to easily increase or
decrease the number of sites. We only have to insert or remove certain TT
cores, respectively.

###  III.3 Higher-order differencing schemes

The simplest way to integrate the TDSE (8) forward by one time step $\Delta t$
is to employ an explicit Euler method

$\Psi(t+\Delta t)=\Psi(t)-\mathrm{i}\Delta
t\mathbf{H}\Psi(t)+\mathcal{O}(\Delta t^{2})$ (17)

which is known to be unstable because neither the time reversal symmetry nor
the symplecticity of the TDSE is conserved. However, when combining Euler
steps forward and backward in time, one arrives at the second order symmetric
Euler (S2) method

$\Psi(t+\Delta t)=\Psi(t-\Delta t)-2\mathrm{i}\Delta
t\mathbf{H}\Psi(t)+\mathcal{O}(\Delta t^{3})$ (18)

which has been in use for more than 40 years in numerical quantum dynamics
where it is known as second order differencing scheme [62, 53]. Not requiring
a special structure of the Hamiltonian (such as the separability, see Sec.
III.4), the S2 scheme is universal and very simple to implement. In passing,
we note the analogy with the Verlet integrator often used for numerical
integration of Newton’s classical equation of motion [63]. Even though this
explicit propagation scheme has the advantages of being time-reversible and
conditionally stable, it is neither symplectic nor strictly unitary nor energy
conserving [48]. However, in practice norm and energy conservation are
approximately realized for small enough time steps $\Delta t$. Since the above
method is a two step method, it requires two initial conditions: (i)
$\Psi(t)$, which is given, and (ii) $\Psi(t-\Delta t)$, which we obtain by
propagating backwards in time. A common way to achieve this is to start by a
first-order scheme (17) for half a time step, and then use the S2 scheme (18)
to propagate backward another half step [53].

Higher order schemes are based on Euler methods of order $\ell_{\mathrm{max}}$
forward ($t+\Delta t$) and backward ($t-\Delta t$) in time

$\Psi(t\pm\Delta
t)=\sum_{\ell=0}^{\ell_{\mathrm{max}}}\frac{1}{\ell!}\left(\mp\mathrm{i}\Delta
t\mathbf{H}\right)^{\ell}\Psi(t)+\mathcal{O}(\Delta
t)^{\ell_{\mathrm{max}}+1}$ (19)

Combining the forward and backward equation then yields time reversible
schemes

$\Psi(t+\Delta t)=\Psi(t-\Delta
t)+\underbrace{\sum_{k=1}^{k_{\mathrm{max}}}\frac{2}{(2k-1)!}\left(-\mathrm{i}\Delta
t\mathbf{H}\right)^{2k-1}}_{=:\mathbf{A}}\Psi(t)+\mathcal{O}(\Delta
t)^{2k_{\mathrm{max}}+1}$ (20)

which offer the advantage of eliminating all even powers of $\Delta t$, in
particular the term with $\mathcal{O}(\Delta t)^{2k_{\mathrm{max}}}$. Hence,
the resulting schemes are of order $2k_{\mathrm{max}}$, and they will be
referred to as S2, S4, S6 for $k_{\mathrm{max}}=1,2,3$ throughout Sec. IV.

In our practical work using TT representations of state vectors and
(Hamiltonian) operators, solving Eq. (18) or Eq. (20) is straightforward:
given $\Psi(t)$, $\Psi(t-\Delta t)$ and $\mathbf{H}$ in the TT format, we can
obtain $\Psi(t+\Delta t)$ by simple tensor additions and multiplications.
However, since TT ranks add under addition and multiply under multiplication,
we have

$\text{rank}(\Psi(t+\Delta t))=\text{rank}(\Psi(t-\Delta
t))+\text{rank}(\mathbf{H})\cdot\text{rank}(\Psi(t))$ (21)

already for the second order scheme of Eq. (18). For the higher-order schemes
of (20), the rank of $\mathbf{H}$ in Eq. (21) has to be replaced by
$\text{rank}(\mathbf{A})\leq\sum_{k=1}^{k_{\mathrm{max}}}\text{rank}(\mathbf{H})^{2k-1}$
where $\mathbf{A}$ is defined in the underbrace of Eq. (20). Hence, the rank
of our approximate $\Psi(t)$ would grow exponentially with the number of time
steps and would become exceedingly large already after very few time steps.
Therefore, one needs to perform a rank reduction after every time step, i.e.,
we approximate the high-rank $\Psi(t+\Delta t)$ resulting from the propagation
step with a low-rank TT. That is, given a core of $\Psi(t+\Delta t)$, we
reshape it into a matrix and then apply a truncated singular value
decomposition (SVD). Depending on the direction of the orthonormalization, the
left- or right-singular vectors then define the updated cores and the non-
orthonormal part is contracted with the next core, see [26, 40] for details.
Applying this procedure successively to each core from one end to the other,
we obtain a low-rank approximation which is then used for the next time step.
The rank reduction requires $N-1$ singular value decompositions on matrices of
size $s\times d\cdot s^{\prime}$ where $d$ is local Hilbert space dimension
and $s$ and $s^{\prime}\leq s$ are the maximum TT ranks of $\Psi(t+\Delta t)$
before and after the orthonormalization. The total cost in each time step can
be estimated by $\mathcal{O}(NR^{2}r^{2}d^{2}+NR^{3}r^{3}d)$, where $R$ and
$r$ are the maximum TT ranks of $\mathbf{A}$ and $\Psi(t)$, respectively. The
first term results from estimating the contraction of $\mathbf{A}$, see Eq.
(20), with $\Psi(t)$, whereas the second term corresponds to the
truncation/orthonormalization of the tensor cores of $\Psi(t-\Delta
t)+\mathbf{A}\Psi(t)$.

###  III.4 Higher order splitting schemes

Since the early days of numerical quantum dynamics, splitting approaches have
been widely used for the matrix exponential (11) governing the formal solution
of the TDSE [64, 49]. These methods are suitable wherever the quantum-
mechanical Hamiltonian can be decomposed into (two or more) exactly (or, at
least, easily) solvable parts, e.g., the kinetic and potential energy
operator. Because these parts of the Hamiltonian typically don’t commute,
deriving approximations to the Baker–Campbell–Hausdorff formula has lead to a
hierarchy of exponential splitting methods for use as integrators in solving
the time discretized TDSE which are explicit and easy to implement. Based on
the classical first order Lie-Trotter (LT) and second order Strang-Marchuk
(SM) schemes, higher-order methods can be obtained as compositions of the
basic methods [65, 66, 67, 68]. For a systematic overview of such composition
methods, the reader is referred to section V.3 of Ref. [69]. Common advantages
of the exponential splitting methods are that they typically are unitary,
symplectic and time-reversible [70, 50].

In what follows, we construct various splitting schemes for Hamiltonians
restricted to NN interactions in quasi one-dimensional chains. Moreover, we
assume that all state vectors and Hamiltonian operators are given in TT format
as SLIM decompositions. Splitting methods have been considered for operators
in TT format, see, e.g., [27], where the case of Lie-Trotter splitting is
explained in detail. In the closely related context of density matrix
renormalization group (DMRG), TT/MPS decompositions are used to construct
time-evolving block decimation (TEBD) schemes for propagations of quantum
state vectors [51, 52], which resemble the techniques described below. For
other closely related splitting methods on TTs, we refer to [45, 46].

The key to splitting the exponential in the evolution (11) is to separate the
operators constituting the Hamiltonian (1) into groups as large as possible.
At the same time, all the operators belonging to one group should commute with
each other such that corresponding computations can be carried out in parallel
[51]. Thus, for an implementation of splitting methods for TTs using SLIM
decomposition (13), all single site operators $S$ commute by construction, and
operators $L$ and $M$ acting on odd/even NN pairs commute within their oddity
groups. More specificially, we can split our Hamiltonian tensor operator
$H\in\mathbb{R}^{(d_{1}\times d_{1})\times\dots\times(d_{N}\times d_{N})}$
into two parts $H_{\text{odd}}$ and $H_{\text{even}}$ which are defined by

$\begin{split}H_{\text{odd}}&=\sum_{\begin{subarray}{c}i=1,\\\
i\,\text{odd}\end{subarray}}^{N}I^{1:i-1}\otimes S^{(i)}\otimes
I^{i+1:N}+\sum_{\begin{subarray}{c}i=1,\\\
i\,\text{odd}\end{subarray}}^{N-1}I^{1:i-1}\otimes\left\llbracket
L^{(i)}\right\rrbracket\otimes\left\llbracket M^{(i+1)}\right\rrbracket\otimes
I^{i+2:N}\\\ &=\sum_{i=1}^{\lceil\frac{N}{2}\rceil}I^{1:2i-2}\otimes
S^{(2i-1)}\otimes
I^{2i:N}+\sum_{i=1}^{\lceil\frac{N-1}{2}\rceil}I^{1:2i-2}\otimes\left\llbracket
L^{(2i-1)}\right\rrbracket\otimes\left\llbracket
M^{(2i)}\right\rrbracket\otimes I^{2i+1:N}\end{split}$

and

$\begin{split}H_{\text{even}}&=\sum_{\begin{subarray}{c}i=1,\\\
i\,\text{even}\end{subarray}}^{N}I^{1:i-1}\otimes S^{(i)}\otimes
I^{i+1:N}+\sum_{\begin{subarray}{c}i=1,\\\
i\,\text{even}\end{subarray}}^{N-1}I^{1:i-1}\otimes\left\llbracket
L^{(i)}\right\rrbracket\otimes\left\llbracket M^{(i+1)}\right\rrbracket\otimes
I^{i+2:N}\\\ &=\sum_{i=1}^{\lfloor\frac{N}{2}\rfloor}I^{1:2i-1}\otimes
S^{(2i)}\otimes
I^{2i+1:N}+\sum_{i=1}^{\lfloor\frac{N-1}{2}\rfloor}I^{1:2i-1}\otimes\left\llbracket
L^{(2i)}\right\rrbracket\otimes\left\llbracket
M^{(2i+1)}\right\rrbracket\otimes I^{2i+2:N},\end{split}$

respectively. Next, we define $H_{i}$, $i=1,\dots,N-1$, to denote the
components

$\begin{split}H_{i}&=I^{1:i-1}\otimes S^{(i)}\otimes
I^{i+1:N}+I^{1:i-1}\otimes L^{(i)}\otimes M^{(i+1)}\otimes I^{i+2:N}\\\
&=I^{1:i-1}\otimes\left(S^{(i)}\otimes I^{(i+1)}+L^{(i)}\otimes
M^{(i+1)}\right)\otimes I^{i+2:N},\end{split}$

whereas the last component with $i=N$ is simply given by

$H_{N}=I^{1:N-1}\otimes S^{(N)}.$

Note that we here only consider non-cyclic systems because defining
$H_{N}=I^{1:N-1}\otimes S^{(N)}+M_{1}\otimes J^{2:N_{1}}\otimes L_{N}$, cf.
(13), would either lead to a significant increase of the computational
complexity of the following schemes or would involve cyclic TT representations
[42] of quantum states, both of which are beyond the scope of the present
work. For the use of SLIM decompositions of cyclic systems, we refer to [41,
40]. With this definition of the $H_{i}$, the tensors $H_{\text{odd}}$ and
$H_{\text{even}}$ can then be written as

$H_{\text{odd}}=\sum_{i\in\mathcal{I}_{\text{odd}}}H_{i},\quad\mathcal{I}_{\text{odd}}=\\{1,3,\dots,2\lceil
N/2\rceil-1\\}$ (22)

and

$H_{\text{even}}=\sum_{i\in\mathcal{I}_{\text{even}}}H_{i},\quad\mathcal{I}_{\text{even}}=\\{2,4,\dots,2\lfloor
N/2\rfloor\\}.$ (23)

The components of $H_{\text{odd}}$ as given in (22) commute pairwise. The same
holds for the components of $H_{\text{even}}$ as given in (23). This implies
that the components of the propagator in Eq. (11) can be expressed as

$\begin{split}\exp(-\mathrm{i}tH_{\text{odd}})&=\prod_{i\in\mathcal{I}_{\text{odd}}}\exp(-\mathrm{i}tH_{i})\\\
&=\prod_{i\in\mathcal{I}_{\text{odd}}}I^{1:i-1}\otimes\exp\left(-\mathrm{i}t\left(S^{(i)}\otimes
I^{(i+1)}+L^{(i)}\otimes M^{(i+1)}\right)\right)\otimes I^{i+2:N}\end{split}$

with $I_{N+1}=1$ and $M_{N+1}=0$, see above. An analogous statement holds for
$\exp(-\mathrm{i}tH_{\text{even}})$. That is, the tensor exponentials of
$-\mathrm{i}tH_{\text{odd}}$ and $-\mathrm{i}tH_{\text{even}}$ can be written
as the product of $\lceil N/2\rceil$ and $\lfloor N/2\rfloor$, respectively,
commuting tensor operators acting only on two sites at most. Thus, we can
directly implement splitting schemes such as first-order Lie-Trotter (LT) or
second-order Strang-Marchuk (SM), see Fig. 2, given by

$\exp(-\mathrm{i}tH)\overset{\text{LT}}{\approx}\exp(-\mathrm{i}tH_{\text{odd}})\exp(-\mathrm{i}tH_{\text{even}})$

and

$\exp(-\mathrm{i}tH)\overset{\text{SM}}{\approx}\exp\left(-\mathrm{i}\frac{t}{2}H_{\text{odd}}\right)\exp(-\mathrm{i}tH_{\text{even}})\exp\left(-\mathrm{i}\frac{t}{2}H_{\text{odd}}\right),$

respectively.

The SM split-operator approach offers the advantage of being unitary,
symplectic, stable, symmetric, and time-reversible [49], However, to obtain
highly accurate results, this algorithm requires using small time steps
because it has only second-order accuracy. Hence, we consider TT-based
splitting schemes with higher order by using the fact that these methods can
be obtained as suitable compositions of SM steps [50], i.e.,

$\exp(-\mathrm{i}tH)=\text{SM}(-\mathrm{i}\gamma_{s}tH)\cdots\text{SM}(-\mathrm{i}\gamma_{1}tH)$

with

$\text{SM}(-\mathrm{i}\gamma tH):=\exp\left(-\mathrm{i}\frac{\gamma
t}{2}H_{\text{odd}}\right)\exp(-\mathrm{i}\gamma
tH_{\text{even}})\exp\left(-\mathrm{i}\frac{\gamma
t}{2}H_{\text{odd}}\right),$

where $s$ is the number of stages and where consistency implies that
$\sum_{j=1}^{s}\gamma_{j}=1$. For more information see Ref. [70], where an
overview of splitting methods with different order is given. In the present
work, we consider the palindromic Yoshida-Neri (YN, order $4$ with 3 stages)
[65] and Kahan-Li (KL, order $8$ with 17 stages) [71] compositions with
coefficients

$\gamma_{1}=\gamma_{3}=1.3512,\quad\gamma_{2}=-1.7024,$

and

$\displaystyle\gamma_{1}=\gamma_{17}=\phantom{-}0.1302,\quad\gamma_{2}=\gamma_{16}=0.5612,\quad\gamma_{3}=\gamma_{15}=-0.3895,\quad\gamma_{4}=\gamma_{14}=0.1588,$
$\displaystyle\gamma_{5}=\gamma_{13}=-0.3959,\quad\gamma_{6}=\gamma_{12}=0.1845,\quad\gamma_{7}=\gamma_{11}=\phantom{-}0.2584,\quad\gamma_{8}=\gamma_{10}=0.2950,$
$\displaystyle\gamma_{9}=-0.6055,$

respectively. Because some of the coefficients are negative, both the YN and
the KL scheme can be considered as combinations of backward and forward
propagations in time.

$\exp\left(H_{\text{odd}}\frac{\Delta
t}{2}\right)$$\exp\left(H_{\text{even}}\Delta
t\right)$$\exp\left(H_{\text{odd}}\frac{\Delta t}{2}\right)$$\Psi(t)$ Figure
2: Strang splitting for SLIM operators for even number $N$ of sites: In each
stage of the splitting scheme, only pairs of cores or single cores are
altered. If no rank reduction is performed between these stages, the ranks of
the resulting TT representation of $\Psi(t+\Delta t)$ are bounded by $rd^{3}$,
where $r$ is the maximum TT rank of $\Psi(t)$.

Note that the composition schemes described above are time-reversible,
symplectic and unitary, i.e., norm-preserving. In analogy to the higher-order
differencing schemes, see Sec. III.3, the advantages of TT-based splitting
methods are that, in contrast to implicit time-integration methods, we do not
have to rely on iterative solvers for linear systems like the _alternating
linear scheme_ (ALS) [54] and we can, in principle, parallelize the tensor
contractions and rank reductions in each stage. In particular, homogeneous
systems can be efficiently simulated by splitting methods since the two-site
tensor exponentials of $H_{\text{odd}}$ and $H_{\text{even}}$ are identical
for all sites. In general, each two-site operator $\exp(\mathrm{i}tH_{i})$ has
a TT representation with rank bounded by $d^{2}$. Thus, by orthonormalizing
(without truncation, see Sec. III.3) between the stages, the ranks of the
resulting tensor after one SM step are alternately bounded by $rd^{2}$ and
$rd^{3}$. In this case, the computational complexity for sufficiently large
$N$ can be estimated by $\mathcal{O}(Nr^{3}d^{8})$, where we assume that in
each stage the matricizations of the two-site operators are multiplied by the
contractions of the corresponding cores and the resulting matrix is reshaped
and decomposed by using an SVD. Note that a subsequent truncation for reducing
the ranks of $\Psi(t+\Delta t)$ to $r$ would take $O(Nr^{3}d^{5})$ floating-
point operations. However, if we only allow TT ranks bounded by some number
$s<rd$ between the SM stages, we can reduce the computational costs to
$O(Ns^{3}d^{4})$. Truncating the ranks of the resulting tensor train to $r$
then only takes $O(Nrs^{2}d)$ operations. In this way, we can significantly
reduce the computational complexity but we have to accept losses in accuracy
since our experiments showed that higher TT ranks are needed particularly in
between the SM stages. Choosing $s=2r$, however, turned out to yield nearly
the same numerical results in our test cases as when allowing arbitrarily
large TT ranks in between the SM stages.

###  III.5 Software

All simulations presented in this work are carried out using our recently
developed Python package WaveTrain which is publicly available via the GitHub
platform [72]. It encompasses various numerical solvers for TISE and TDSE for
Hamiltonians in TT form using a SLIM representation, among them all the
integrators mentioned in the present work. WaveTrain builds on Scikit-TT [73],
an open-source tensor train toolbox for Python based on NumPy and SciPy. For
short chains, $N\leq 5$, our results obtained with WaveTrain are in excellent
agreement with values obtained from the conventional, grid-based solvers for
general Hamiltonians available within the WavePacket software package [74, 75,
76].

##  IV Results and Discussion

In this Section, we aim at assessing the performance and accuracy of the TT
techniques, in particular the various propagators to numerically solve the
TDSE for excitons, for phonons, as well as for the coupled systems. First of
all, our comparisons are based on the conservation of norm and energy over
time. Where possible, we also compare our numerical results with semi-
analytical solutions which are available for short chains, see Sec. IV.3,
and/or for harmonic phonons, see Sec. IV.4.

###  IV.1 Model parameters

Unless stated otherwise, in our computer experiments we restrict ourselves to
the case of non-periodic, homogeneous chains, with the model parameters
$\alpha,\beta,m,\nu,\omega,\sigma$ being equal for all sites, with their
numerical values adapted from our previous work [40]. For the excitonic
Hamiltonian of (2), we choose a local excitation energy $\alpha=0.1E_{h}$
($\approx 2.7$ eV) which is of the order of typical band gaps in organic
semiconductors, whereas the NN coupling energy is chosen one order of
magnitude smaller, $\beta=-0.01E_{h}$. The phononic Hamiltonian of (3) is
parametrized in terms of mass-weighted displacement coordinates,
$\tilde{R}=\sqrt{m}R$, here with unit masses, $m=1$, and harmonic frequencies
$\nu=10^{-3}E_{h}/\hbar$ ($\approx 220$ cm-1) and $\omega=\sqrt{2}\times
10^{-3}E_{h}/\hbar$ for the restraining and for the NN oscillators,
respectively. The coupling constant occurring in the (symmetric) EPC mechanism
of Eq. (7) is chosen even smaller than the vibrational frequencies,
$\sigma=2\times 10^{-4}E_{h}/a_{0}$. Note that for this value of the EPC
constant, the stationary states obtained as solutions of the TISE display
mutual trapping of excitons and phonons, the former ones localized with a full
width at half maximum of about 7 sites, see our previous work [40]. Note that
the use of mass-weighted coordinates implies a scaling of the EPC constant
with $\tilde{\sigma}=\sigma/\sqrt{m}$. Thus, it is straight-forward to
transfer all our results described below to arbitrary values of the particle
masses involved. Except where stated otherwise, we employ a basis set
consisting of $d^{\mathrm{(ex)}}=2$ (electronic two-state model) and
$d^{\mathrm{(ph)}}=8$ vibrational basis functions (from the second
quantization introduced in Sec. II).

###  IV.2 Time scales

For a homogeneous, non-periodic chain of infinite length, the TDSE for the
purely excitonic two-state Hamiltonian of Eq. (2) can be solved analytically:
For the case of an initial excitation localized at site $i_{0}$, the
probability of finding the exciton at time $t$ on site $i$ is given by

$\mathcal{P}_{i}(t)=J_{i-i_{0}}^{2}(2|\beta|t)$ (24)

where the $J_{n}$ are $n$-th order Bessel functions of the first kind [77, 7].
In practice, these solutions are of limited use as a benchmark for numerical
quantum dynamics because typically the excitations quickly reach the end of a
finite length chain. However, the above solution is instrumental in defining a
time scale for the excitonic energy transport as the time between the first
maximum of $J_{0}$ (at $t=0$) and the first maximum of $J_{1}$ (at
$2|\beta|t/\approx 1.8412$). Hence, for our example we suggest a rounded value
for the time scale, $\tau^{\mathrm{(ex)}}\equiv 1/|\beta|=100$ in atomic
units, or approximately 2.4 fs. In the absence of analytical results for
phonon dynamics given by the Hamiltonian of Eq. (4), we simply define an
analogous time scale for phononic energy transport,
$\tau^{\mathrm{(ph)}}\equiv 1/\nu=1000$, or approximately 24 fs. This value
approximately coincides with our empirical observations of a vibrational
excitation moving from one site to its nearest neighbors. Similarly, the time
scale of the exciton-phonon coupling of (7) is defined as
$\tau^{\mathrm{(epc)}}\equiv 1/|\sigma|=5000$, or approximately 120 fs,
characterizing the redistribution of energy between excitonic and phononic
subsystems.

###  IV.3 Comparison with quantum dynamics

For short chains, it is straight-forward to obtain semi-analytical solutions
to the TDSE (8) by calculating the matrix exponential in the formal solution
(11) directly. Assuming a maximum matrix size of $4096\times 4096$, the semi-
analytical approach is limited to $N=12$, $N=4$, or $N=3$ for excitons with
$d^{\mathrm{(}ex)}=2$, phonons with $d^{\mathrm{(}ph)}=8$, or coupled systems
with $d=d^{\mathrm{(}ex)}d^{\mathrm{(}ph)}=16$, respectively.

First, we run semi-analytical propagations with 100 main time steps of length
$\tau^{\mathrm{(ex)}}/2=50$ for excitons with and without EPC and
$\tau^{\mathrm{(ph)}}/2=500$ for phonons only. As expected, these simulations
display preservation of norm and energy to machine precision. Subsequently, we
use those results as a reference in assessing the quality of TT-based
propagations of the same duration where we split the main time steps into 1,
2, 5, 10, 20, 50, 100, 200, 500, or 1000 sub-steps each. Propagations are
carried out using the second order symmetric Euler (S2) method, see Sec.
III.3, as well as four different splitting schemes, i.e., Lie-Trotter (LT),
Strang-Marchuk (SM), Yoshida-Neri (YN), and Kahan-Li (KL), see Sec. III.4. In
all simulations of exciton dynamics with or without EPC, the initial state
corresponds to a single site near the center of the chain being excited. Our
simulations of phonon dynamics were initialized with the vibration of that
site being in a coherent state where a mean displacement $\langle\Delta
R\rangle=1$ has been chosen [78]. Note that this value is small enough such
that the truncation of the vibrational basis (here $d^{\mathrm{(ph)}}=8$) does
not have a notable influence on the precision of the results. All other sites
were prepared in their excitonic and/or vibrational ground states.

Our results for the three classes of systems investigated are displayed in
Figs. 3, 4, and 5. The columns of these figures correspond to different values
of $r$, the maximum ranks of the solutions, characterizing the flexibility of
the solutions. Note that for given dimension $d$ of the Hilbert spaces for
each of the $N$ sites, the maximum rank possible, $r_{\mathrm{max}}$, of a
tensor train (after left- and right-orthonormalization) is given by
$d^{\left\lfloor N/2\right\rfloor}$ where the values of $d$ are chosen as 2,
8, or 16 for excitons, phonons, or coupled systems, respectively. Hence, we
have $r_{\mathrm{max}}=64$ both for pure excitons ($N=12$) and pure phonons
($N=4$) whereas we have $r_{\mathrm{max}}=16$ for coupled systems ($N=3$). The
CPU times for a single thread of a Xeon Skylake 6130 processor, shown in the
upper rows of the figures, increase linearly with the number of time sub-
steps. Owing to the simplicity of the integrator (18), the S2 is the fastest
of the 5 integrators compared here. It is followed by the first order (LT) and
second order (SM) splitting schemes with slightly higher computational effort.
Note that results for S2, LT, SM propagations are partly not available for the
larger sub-step sizes where they become unstable. The elapsed CPU times are
considerable higher for the 4-th order (YN) and 8-th order (KL) schemes.

The panels in the second rows of Figs. 3–5 show the root mean squared
deviation (RMSD, averaged over the 100 time steps considered) of the norm of
the quantum-mechanical states from unity. The S2 method reproduces the
preservation of norm moderately well with RMSD values typically between
$10^{-2}$ and $10^{-11}$, steeply descending with decreasing step size. RMSD
values between $10^{-9}$ and $10^{-13}$ are often found for the splitting
methods, which is expected because of the unitarity of those schemes. However,
this favorable behavior is only found where the ranks of the solutions, $r$,
are high enough. While this is already the case for $r=2$ for the exciton
results shown in Fig. 3, one has to go close to the maximum ranks for the
other two systems, i.e., $r=r_{\mathrm{max}}=64$ for the phonons shown in Fig.
4 or $r=r_{\mathrm{max}}=16$ for the coupled systems, see Fig. 5. A few
outliers are found for the YN scheme for excitons, see Fig. 3, which are
apparently caused by the rank reduction in between the SM stages, see Section
III.4.

A similar behavior is found in our (relative) RMSDs of energies w.r.t. their
initial values and for the (absolute) RMSDs of the TT-based solutions from the
semi-analytical solutions, as can be seen in the third and fourth rows of
Figs. 3–5. Once the ranks of the solutions, $r$, exceed the above-mentioned
thresholds, these two quantities can be brought down to $10^{-11}$ or even
lower for small enough time steps. While such a high quality of the numerical
solutions requires rather small time steps for the lower order propagators, it
is reached already for much larger time steps when using the higher order
schemes. As expected, our data in the last rows (RMSDs of the state vectors)
of the right columns of Figs. 3–5 are well approximated by power laws with
exponents of approximately 1 (LT), 2 (S2, SM), 4 (YN), and 8 (KL). The
deviating behavior of the RMSDs when applying YN and KL for small time steps
can be explained by the high number of SVDs used for
truncating/orthonormalizing the intermediate solutions, see Sec. III.4. In
this case, the (numerical) errors induced by the core decompositions appear to
outweigh the time-discretization error.

To summarize, Fig. 6 shows a synopsis of our results for excitons, phonons,
and coupled systems. Again, the respective chain lengths of $N=12$, $N=4$, and
$N=3$ are equivalent in the sense that the dimensions of the underlying
Hilbert spaces,
$(d^{\mathrm{(}ex)})^{12}=(d^{\mathrm{(}ph)})^{4}=(d^{\mathrm{(}ex)}d^{\mathrm{(}ph)})^{3}=4096$,
are equal. In this figure, the accuracy, calculated as RMSD between semi-
analytical and TT-based numerical results, is shown as a function of the total
CPU time. Essentially, the comparison shows the superiority of the higher-
order splitting integrators when high precision is required. Only the 4-th
order (YN) and 8-th order schemes (KL) are able to reach a precision of
$10^{-8}$ and below. The most precise results are obtained for the KL
integrator, with RMSDs below $10^{-11}$ for suitable choice of the time step
size. In particular, the CPU times required to reach the maximal precision are
about $10^{2}$ s for excitons and phonons, whereas about $10^{3}$ s are found
for the coupled systems. Since we here consider rather short chains,
especially for phonons and coupled systems, we observe only a weak influence
of the maximum TT rank on the computational costs. But the complexity of the
simulations employing the splitting schemes still depends on $d^{4}$, see Sec.
III.4, which we identify as the reason for the higher CPU times required for
computing accurate solutions for the coupled exciton-phonon systems.

Figure 3: Quantum dynamics of purely excitonic chains with $d=2$ and $N=12$
sites. From left to right: Maximum number of ranks of state vectors
increasing. From top to bottom: CPU-time versus size of temporal sub-steps,
deviation of norm from unity, relative deviation of energy from initial value,
deviation of state vectors from semi-analytical reference data. For second
order symmetric Euler (S2), Lie-Trotter (LT), Strang-Marchuk (SM), Yoshida-
Neri (YN), and Kahan-Li (KL) integration methods Figure 4: Quantum dynamics of
purely phononic chains with $d=8$ and $N=4$ sites. For details, see caption of
Fig. 3 Figure 5: Quantum dynamics of coupled excitons and phonons for $d=16$
and $N=3$ sites. For details, see caption of Fig. 3

Figure 6: Accuracy of solutions versus computational effort. From left to
right: excitons ($N=12$, $r=8$, $d=2$), phonons ($N=4$, $r=64$, $d=8$), and
coupled systems ($N=3$, $r=16$, $d=16$).

###  IV.4 Comparison with classical dynamics

For the dynamics of phonons, reference data for assessing our TT-based results
can also be generated on the basis of the quantum-classical correspondence.
According to the Ehrenfest theorem, the quantum-mechanical expectation values
of positions and momenta coincide with results from classical trajectories, as
long as the vibrational Hamiltonian is a polynomial of order not higher than
two [79]. This is indeed the case for our vibrational model Hamiltonian (3)
which can be re-formulated as
$H^{\mathrm{(ph)}}=R^{T}\mathbf{F}R+P^{T}\mathbf{G}P$ where $\mathbf{F}$ and
$\mathbf{G}$ stand for the Hessians of the Hamiltonian function with respect
to positions, $R$, and momenta, $P$, respectively. Under these assumptions,
Hamilton’s classical equations of motion can be solved in a semi-analytical
manner using a matrix exponential

$\left(\begin{array}[]{c}R(t)\\\
P(t)\end{array}\right)=\exp\left[\left(\begin{array}[]{cc}\mathbf{0}&\mathbf{G}\\\
-\mathbf{F}&\mathbf{0}\end{array}\right)t\right]\left(\begin{array}[]{c}R(0)\\\
P(0)\end{array}\right)$ (25)

where a symplectic phase space formulation [69] is used and where $\mathbf{0}$
stands for a zero matrix of size $N\times N$. Our classical simulations are
initialized with a single site (at the center of the chain) displaced by
$\Delta R=1$ while all other sites are at their equilibrium positions. We
generate a set of classical simulations with a duration of 100 time steps of
size $\tau^{\mathrm{(ph)}}/2=500$ for $N\in\\{4,8,16\\}$. As expected for
these system sizes, the semi-analytical solutions display excellent
preservation of energy within machine accuracy.

Next, we conduct TT-based quantum propagations of the same duration
($50\tau^{\mathrm{(ph)}}=50000$) using both splitting (LT, SM) and
differencing (S2, S4) propagators, again sub-dividing the main time steps into
1, 2, 5, 10, 20, 50, 100, 200, 500, 1000 sub-steps each. Moreover, we restrict
our analysis to the results for $r=32$. For shorter chains with $N=4$, this
value has already been shown to yield reasonable results for the RMSDs of the
state vectors, see Fig. 4. Note that we here exclude the higher-order
splitting schemes (YN, KL) due to the exceedingly high computational effort.
Furthermore, we find that S6 as well as S8 splitting methods yield results
that are almost identical to those obtained for the S4 method which is the
reason why we do not show them anywhere in Fig. 7. We explain this effect as a
consequence of the used orthonormalizations and rank reductions which
eliminate the information stored in the terms corresponding to small singular
values and expect S6 to show better results than S4 when allowing for higher
TT ranks.

In analogy to our classical simulations of phonon dynamics, we initialize our
quantum simulations by using a coherent state for a single site with the mean
position matching the displacement given above. Our results are shown in Fig.
7 where RMSD comparisons of quantum versus classical position values of the
vibrational coordinates are shown in the last row. In contrast to our
simulations of phonon dynamics presented in Sec. IV.3, we use a larger basis
set with $d^{\mathrm{(ph)}}=16$ here which has turned out to be mandatory when
aiming at high accuracy for the comparison of vibrational coordinates.

The CPU times for the four integrators used here are shown in the top panels
of Fig. 7 for chains of 4, 8, and 16 repeat units. We see an increase in the
CPU times of roughly one order of magnitude upon doubling the chain length
both for the differencing schemes (S2, S4) and for the splitting (LT, SM)
schemes. As stated in Sec. III.3 and Sec. III.4, the computational
complexities of the differencing and splitting schemes, respectively, scale
linearly with the chain length. The reason why this behavior cannot be
observed in Fig. 7 is the comparatively small number of sites and the fact
that the ranks $r_{1}$ and $r_{N-1}$ of a TT representation of a quantum state
vector with dimensions $d$ are always bounded by $d$. Thus, in this
experiment, every tensor train with $N=4$ has only one rank bounded by $32$,
whereas it has five when $N=8$ and $13$ when $N=16$. Hence, due to the
dominant computational costs of contractions and orthonormalizations of TT
cores with maximum ranks, the CPU times in Fig. 7 are not expected to show a
linear scaling when doubling the chain length. However, by bounding the TT
ranks of the intermediate solutions in between the SM stages to $2r$, cf.
Section III.4, we are able to reduce the increase of the CPU times for
splitting methods significantly while observing almost no loss in accuracy.

The behavior of the norm of the state vectors is shown in the second row of
Fig. 7. The norm is preserved up to machine precision for the unitary LT and
SM splitting schemes only for $N=4$ whereas notable deviations from unity are
found for $N=8$ and $N=16$. This behavior is attributed to the finite ranks,
$r=32$, used here, see also Fig. 4 and our discussion of rank effects in Sec.
IV.3. In contrast, for the S2 and S4 differencing schemes, the norm is only
approximately conserved for short time steps. For $N=4$ and $N=8$, the
deviation from unity scales is approximately proportional to $\tau^{-4}$ for
the S2 scheme as expected [62], whereas it scales as $\tau^{-3}$ for the S4
scheme. Note that the RMSDs for longer chains ($N=16$) are much larger due to
finite-rank effects.

The third row of Fig. 7 shows that the energies are not preserved exactly, as
expected for our choice of integrators which are known not to conserve
quadratic invariants rigorously. The comparison of our results does not show
any characteristic differences between the different integrators for longer
chains ($N=8$ and $N=16$). Only for short chains ($N=4$), the double-
logarithmic representations show an approximately linear behavior with slopes
of 1.0 (LT), 2.0 (SM), 2.5 (S2) and 3.0 (S4).

Finally, we consider the deviations (RMSD) of our numerical TT quantum results
for the calculated positions from our (semi-analytical) classical results
described at the beginning of this section. Again, only for short chains we
find straight lines in the bottom panels of Fig. 7, this time for $N=4$ with
slopes of 1.0 (LT), 2.0 (SM), 2.0 (S2) and 4.0 (S4) as expected. Also for the
longer chains with $N=8$ and $N=16$ we find that the higher order differencing
schemes are superior to the lower order ones. Except for the largest time
steps, the SM scheme clearly beats the LT scheme, and the S4 scheme is often
several orders of magnitude more precise than the S2 scheme which, however, is
the fastest.

Figure 7: Quantum dynamics of phonons for $N=4$, $N=8$, and $N=16$ sites (from
left to right) for $d=16$ and $r=32$. Note that the last row shows the RMSDs
of the displacements (position coordinates) of the oscillators

##  V Conclusions and Prospects

In the present work, we have shown how to apply tensor-trains (TT) to solve
the time-dependent Schrödinger equation for chain-like systems with nearest
neighbor (NN) interactions only. The aim of this approach is to mitigate the
curse of dimensionality, i. e., to reduce the memory consumption as well as
the computational costs, as much as possible. Throughout this study, we have
employed memory-efficient SLIM representations to construct TT representations
of the quantum-mechanical Hamiltonians and state vectors.

Two classes of explicit integrators to solve the Schrödinger equation are
investigated here. First, we consider time-symmetrized Euler integrators, see
Sec. III.3, based on the second order differencing method regularly used since
the early days of numerical quantum dynamics. Second, we study splitting
propagators such as the well-known first order Lie-Trotter and the second-
order Strang-Marchuk schemes, see Sec. III.4, originally based on a splitting
of the Hamiltonian into kinetic and potential energy terms. However, in the
present context of chain-like systems we build on a splitting of the
Hamiltonian into groups of interleaved NN interactions commutating within each
group. Both for the differencing and splitting methods, we have extended the
traditional schemes to higher order ones.

This work presents rigorous tests of these approaches for three different
classes of Hamiltonians, i.e., for two-state excitons, for harmonic acoustic
phonons, and for coupled exciton-phonon systems modeled in terms of Frenkel-
Fröhlich-Holstein type Hamiltonians. In particular, we have conducted
extensive investigations of the influence of the TT ranks, the time step size,
and the orders of the individual propagators for homogeneous and non-periodic
chains of various lengths.

The focus of our studies in Sec. IV.3 is on the dependence of the conservation
properties as well as the approximation quality on the maximum number of
ranks, $r$, of the solutions. Typically we find a threshold behavior, i.e.,
once the number $r$ exceeds a certain threshold, the quality of the solutions
increases drastically. For the three classes of physical systems investigated
here, this threshold is found at very different values of $r$. For example,
for an excitonic chain of length $N=12$, we find excellent agreement with
semi-analytic results already for $r=2$. This may be due to the fact that the
exciton dynamics is confined to the Fock space of states with one excitation
which also contains our initial state. (Note that the number of excitons is a
conserved quantity because it commutes with the Hamiltonian (2)). In marked
contrast, for a phononic chain of length $N=4$, this threshold is already
beyond $r=50$. Nevertheless, we summarize that the TT-based propagators
provide excellent solutions once the ranks are chosen large enough. In
particular the higher order splitting methods such as the YN (4-th order) and
KL (8-th order) schemes are recommended where very high precision (close to
machine precision) is required.

The main objective of our studies presented in Sec. IV.4 is the dependence of
the approximation quality on the number of sites, $N$. Here, splitting methods
are not our preferred choice. Instead, differencing methods are more efficient
for longer chains, i.e. the S2 scheme is the faster whereas the S4 scheme is
more accurate than LT or SM splitting. However, one has to be aware that even
S4 achieves only moderate precision.

As an outlook toward future work, we expect our approach to be extended to
more complex scenarios than the ones investigated in the present work. For
example, considering more than one vibrational degree of freedom per site into
account is necessary, e.g., for realistically modeling the conformational
dynamics of $\pi$-conjugated polymer chains [9, 10]. Obviously, this requires
larger dimensions, $d$, of the local Hilbert spaces. In the practice of our
approaches using TT based propagators, in that case differencing methods are
expected to be superior to splitting methods because of the more favorable
scaling of the computational costs with $d$ for the former ones, see our
estimates of the complexities at the end of Secs. III.3 and III.4.

Finally we remark that we have not discussed implicit integrators. For
example, the implicit trapezoidal scheme offers the intriguing advantage of
conserving quadratic invariants such as the energy associated with the
Hamiltonians proposed in Sec. II. Unfortunately, in our test calculations this
integrator was found to be considerably slower than any of the explicit
schemes considered. This is because it involves the solution of a large set of
equations for every time step. In the context of TT techniques, this is
typically done by means of the alternating linear scheme [54] which is rather
time-consuming, see also our previous TISE work [40]. However, implicit
schemes may become more attractive when combined with spatially and/or
temporally adaptive schemes, similar to the trapezoidal rule for adaptive
integration of Liouville equations (TRAIL) [80, 81].

###### Acknowledgements.

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy – The Berlin Mathematics
Research Center MATH+ (EXC-2046/1, project ID: 390685689) and by the CRC 1114
“Scaling Cascades in Complex Systems” funded by the Deutsche
Forschungsgemeinschaft (project ID: 235221301, project B06). Felix Henneke (FU
Berlin) is acknowledged for insightful discussions and Jerome Riedel for
valuable help with the implementation of the WaveTrain Python codes. Moreover,
the authors would like to thank the HPC Service of ZEDAT, FU Berlin, for
generous allocation of computing resources.

## Author declarations

### Conflict of Interest

No potential conflict of interest was reported by the authors.

## Data Availability

The Python and Bash scripts used to generate the simulation results shown in
Figs. 3–7 are openly available in the Zenodo repository at
https://doi.org/10.5281/zenodo.7351812. Moreover, we used version 1.0 of our
open-source WaveTrain software [72], based on version 1.2 of the scikit_tt
tensor train package, both of which are freely available from the GitHub
platform at https://github.com/PGelss/.

## References

  * Mikhnenko  _et al._ [2015] O. V. Mikhnenko, P. W. M. Blom, and T.-Q. Nguyen, Exciton diffusion in organic semiconductors, Energy Environ. Sci. 8, 1867 (2015).
  * Kranz and Elstner [2016] J. J. Kranz and M. Elstner, Simulation of Singlet Exciton Diffusion in Bulk Organic Materials, Journal of Chemical Theory and Computation 12, 4209 (2016).
  * Schröter  _et al._ [2015] M. Schröter, S. D. Ivanov, J. Schulze, S. P. Polyutov, Y. Yan, T. Pullerits, and O. Kühn, Exciton-vibrational coupling in the dynamics and spectroscopy of Frenkel excitons in molecular aggregates, Physics Reports 567, 1 (2015).
  * Zhugayevych and Tretiak [2015] A. Zhugayevych and S. Tretiak, Theoretical Description of Structural and Electronic Properties of Organic Photovoltaic Materials, Annual Review of Physical Chemistry 66, 305 (2015).
  * Devreese and Alexandrov [2009] J. T. Devreese and A. S. Alexandrov, Fröhlich polaron and bipolaron: Recent developments, Reports on Progress in Physics 72, 066501 (2009).
  * Davydov [1985] A. S. Davydov,  _Solitons in Molecular Systems_ (Reidel, 1985).
  * Scott [1992] A. Scott, Davydov’s soliton, Physics Reports 217, 1 (1992).
  * Georgiev and Glazebrook [2019] D. D. Georgiev and J. F. Glazebrook, On the quantum dynamics of Davydov solitons in protein $\alpha$-helices, Physica A: Statistical Mechanics and its Applications 517, 257 (2019).
  * Binder  _et al._ [2018] R. Binder, D. Lauvergnat, and I. Burghardt, Conformational Dynamics Guides Coherent Exciton Migration in Conjugated Polymer Materials: First-Principles Quantum Dynamical Study, Physical Review Letters 120, 227401 (2018).
  * Di Maiolo  _et al._ [2020] F. Di Maiolo, D. Brey, R. Binder, and I. Burghardt, Quantum dynamical simulations of intra-chain exciton diffusion in an oligo ( para -phenylene vinylene) chain at finite temperature, The Journal of Chemical Physics 153, 184107 (2020).
  * Kosloff [1988] R. Kosloff, Time-dependent quantum-mechanical methods for molecular dynamics, The Journal of Physical Chemistry 92, 2087 (1988).
  * Lenz [1951] F. Lenz, The Ratio of Proton and Electron Masses, Physical Review 82, 554 (1951).
  * Bornemann  _et al._ [1996] F. A. Bornemann, P. Nettesheim, and C. Schütte, Quantum-classical molecular dynamics as an approximation to full quantum dynamics, The Journal of Chemical Physics 105, 1074 (1996).
  * Nettesheim  _et al._ [1996] P. Nettesheim, F. A. Bornemann, B. Schmidt, and C. Schütte, An explicit and symplectic integrator for quantum-classical molecular dynamics, Chemical Physics Letters 256, 581 (1996).
  * Choi and Vaníček [2021] S. Choi and J. Vaníček, High-order geometric integrators for representation-free Ehrenfest dynamics, The Journal of Chemical Physics 155, 124104 (2021).
  * Tully [1990] J. C. Tully, Molecular dynamics with electronic transitions, The Journal of Chemical Physics 93, 1061 (1990).
  * Xie  _et al._ [2020] W. Xie, D. Holub, T. Kubaš, and M. Elstner, Performance of Mixed Quantum-Classical Approaches on Modeling the Crossover from Hopping to Bandlike Charge Transport in Organic Semiconductors, Journal of Chemical Theory and Computation 16, 2071 (2020).
  * Nelson  _et al._ [2020] T. R. Nelson, A. J. White, J. A. Bjorgaard, A. E. Sifain, Y. Zhang, B. Nebgen, S. Fernandez-Alberti, D. Mozyrsky, A. E. Roitberg, and S. Tretiak, Non-adiabatic Excited-State Molecular Dynamics: Theory and Applications for Modeling Photophysics in Extended Molecular Materials, Chemical Reviews 120, 2215 (2020).
  * Freixas  _et al._ [2021] V. M. Freixas, A. J. White, T. Nelson, H. Song, D. V. Makhov, D. Shalashilin, S. Fernandez-Alberti, and S. Tretiak, Nonadiabatic Excited-State Molecular Dynamics Methodologies: Comparison and Convergence, Journal of Physical Chemistry Letters 12, 2970 (2021).
  * Trugman  _et al._ [2004] S. A. Trugman, L.-C. Ku, and J. Bonča, Jahn–Teller and the Dynamics of Polaron Formation, Journal of Superconductivity 17, 193 (2004).
  * Ku and Trugman [2007] L.-C. Ku and S. A. Trugman, Quantum dynamics of polaron formation, Physical Review B 75, 014307 (2007).
  * Beck  _et al._ [2000] M. H. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer, The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets, Physics Reports 324, 1 (2000).
  * Meyer  _et al._ [2009] H. D. Meyer, F. Gatti, and G. A. Worth, eds.,  _Multidimensional quantum dynamics: MCTDH theory and applications_ (Wiley-VCH, 2009).
  * Hitchcock [1927] F. L. Hitchcock, The expression of a tensor or a polyadic as a sum of products, Journal of Mathematics and Physics 6, 164 (1927).
  * Tucker [1964] L. R. Tucker, The extension of factor analysis to three-dimensional matrices, in  _Contributions to mathematical psychology_ , edited by H. Gulliksen and N. Frederiksen (Holt, Rinehart and Winston, 1964) pp. 110–127.
  * Oseledets [2011] I. V. Oseledets, Tensor-train decomposition, SIAM Journal on Scientific Computing 33, 2295 (2011).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117 (2014).
  * Borrelli and Gelin [2016] R. Borrelli and M. F. Gelin, Quantum electron-vibrational dynamics at finite temperature: Thermo field dynamics approach, Journal of Chemical Physics 145, 10.1063/1.4971211 (2016).
  * Bose and Walters [2021] A. Bose and P. L. Walters, A Multisite Decomposition of the Tensor Network Path Integrals, The Journal of Chemical Physics 156, 024101 (2021).
  * Gelß  _et al._ [2022] P. Gelß, S. Klus, Z. Shakibaei, and S. Pokutta, Low-rank tensor decompositions of quantum circuits (2022), arXiv:arXiv:2205.09882 .
  * Karahan  _et al._ [2015] E. Karahan, P. A. Rojas-López, M. L. Bringas-Vega, P. A. Valdés-Hernández, and P. A. Valdes-Sosa, Tensor analysis and fusion of multimodal brain images, Proceedings of the IEEE 103, 1531 (2015).
  * Chatzichristos  _et al._ [2019] C. Chatzichristos, E. Kofidis, M. Morante, and S. Theodoridis, Blind fMRI source unmixing via higher-order tensor decompositions, Journal of Neuroscience Methods 315, 17 (2019).
  * Erol and Hunyadi [2022] A. Erol and B. Hunyadi, Chapter 12 - tensors for neuroimaging: A review on applications of tensors to unravel the mysteries of the brain, in  _Tensors for Data Processing_, edited by Y. Liu (Academic Press, 2022) pp. 427–482.
  * Klus  _et al._ [2018] S. Klus, P. Gelß, S. Peitz, and C. Schütte, Tensor-based dynamic mode decomposition, Nonlinearity 31, 3359 (2018).
  * Gelß  _et al._ [2019] P. Gelß, S. Klus, J. Eisert, and C. Schütte, Multidimensional approximation of nonlinear dynamical systems, Journal of Computational and Nonlinear Dynamics 14, 061006 (2019).
  * Kargas and Sidiropoulos [2020a] N. Kargas and N. D. Sidiropoulos, Nonlinear system identification via tensor completion, in  _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 34 (2020) pp. 4420–4427.
  * Stoudenmire and Schwab [2016] E. Stoudenmire and D. J. Schwab, Supervised learning with tensor networks, in  _Advances in Neural Information Processing Systems_ , Vol. 29, edited by D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (Curran Associates, Inc., 2016).
  * Klus and Gelß [2019] S. Klus and P. Gelß, Tensor-based algorithms for image classification, Algorithms 12, 240 (2019).
  * Kargas and Sidiropoulos [2020b] N. Kargas and N. D. Sidiropoulos, Supervised learning via ensemble tensor completion, in  _2020 54th Asilomar Conference on Signals, Systems, and Computers_ (2020) pp. 196–199.
  * Gelß  _et al._ [2022] P. Gelß, R. Klein, S. Matera, and B. Schmidt, Solving the time-independent Schrödinger equation for chains of coupled excitons and phonons using tensor trains, The Journal of Chemical Physics 156, 024109 (2022).
  * Gelß  _et al._ [2017] P. Gelß, S. Klus, S. Matera, and C. Schütte, Nearest-neighbor interaction systems in the tensor-train format, Journal of Computational Physics 341, 140 (2017).
  * Hackbusch [2012] W. Hackbusch,  _Tensor Spaces and Numerical Tensor Calculus_, Vol. 42 (2012).
  * Mainali  _et al._ [2021] S. Mainali, F. Gatti, D. Iouchtchenko, P.-N. Roy, and H.-D. Meyer, Comparison of the multi-layer multi-configuration time-dependent Hartree (ML-MCTDH) method and the density matrix renormalization group (DMRG) for ground state properties of linear rotor chains, The Journal of Chemical Physics 154, 174106 (2021).
  * Serwatka and Roy [2022] T. Serwatka and P.-N. Roy, Ground state of asymmetric tops with DMRG: water in one dimension, The Journal of Chemical Physics 156, 044116 (2022).
  * Lubich and Oseledets [2014] C. Lubich and I. V. Oseledets, A projector-splitting integrator for dynamical low-rank approximation, BIT Numerical Mathematics 54, 171 (2014).
  * Lubich  _et al._ [2015] C. Lubich, I. V. Oseledets, and B. Vandereycken, Time Integration of Tensor Trains, SIAM Journal on Numerical Analysis 53, 917 (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, Physical Review B 94, 165116 (2016).
  * Choi and Vaníček [2019] S. Choi and J. Vaníček, Efficient geometric integrators for nonadiabatic quantum dynamics. I. The adiabatic representation, The Journal of Chemical Physics 150, 204112 (2019).
  * Feit  _et al._ [1982] M. D. Feit, J. A. Fleck, and A. Steiger, Solution of the Schrödinger equation by a spectral method, Journal of Computational Physics 47, 412 (1982).
  * Roulet  _et al._ [2019] J. Roulet, S. Choi, and J. Vaníček, Efficient geometric integrators for nonadiabatic quantum dynamics. II. the diabatic representation, Journal of Chemical Physics 150, 204113 (2019).
  * Volokitin  _et al._ [2019] V. Volokitin, I. Vakulchyk, E. Kozinov, A. Liniov, I. Meyerov, M. Ivanchenko, T. Laptyeva, and S. Denisov, Propagating large open quantum systems towards their asymptotic states: cluster implementation of the time-evolving block decimation scheme, Journal of Physics: Conference Series 1392, 012061 (2019).
  * 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, Annals of Physics 411, 167998 (2019).
  * Leforestier  _et al._ [1991] C. Leforestier, R. Bisseling, C. Cerjan, M. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, O. Roncero, and R. Kosloff, A comparison of different propagation schemes for the time dependent Schrödinger equation, Journal of Computational Physics 94, 59 (1991).
  * Holtz  _et al._ [2012] S. Holtz, T. Rohwedder, and R. Schneider, The Alternating Linear Scheme for Tensor Optimization in the Tensor Train Format, SIAM Journal on Scientific Computing 34, A683 (2012).
  * Hedley  _et al._ [2017] G. J. Hedley, A. Ruseckas, and I. D. W. Samuel, Light Harvesting for Organic Photovoltaics, Chemical Reviews 117, 796 (2017).
  * Dexter [1953] D. L. Dexter, A theory of sensitized luminescence in solids, The Journal of Chemical Physics 21, 836 (1953).
  * Dexter  _et al._ [1969] D. L. Dexter, R. S. Knox, and T. Förster, The radiationless transfer of energy of electronic excitation between impurity molecules in crystals, physica status solidi (b) 34, K159 (1969).
  * Scott [1991] A. C. Scott, Davydov’s soliton revisited, Physica D: Nonlinear Phenomena 51, 333 (1991).
  * Oseledets [2009] I. V. Oseledets, A new tensor decomposition, Doklady Mathematics 80, 495 (2009).
  * Oseledets and Tyrtyshnikov [2009] I. V. Oseledets and E. E. Tyrtyshnikov, Breaking the curse of dimensionality, or how to use SVD in many dimensions, SIAM Journal on Scientific Computing 31, 3744 (2009).
  * Kazeev  _et al._ [2013] V. Kazeev, O. Reichmann, and C. Schwab, Low-rank tensor structure of linear diffusion operators in the TT and QTT formats, Linear Algebra and its Applications 438, 4204 (2013).
  * Askar and Cakmak [1978] A. Askar and A. S. Cakmak, Explicit integration method for the time-dependent Schrödinger equation for collision problems, The Journal of Chemical Physics 68, 2794 (1978).
  * Verlet [1967] L. Verlet, Computer ”Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules, Physical Review 159, 98 (1967).
  * Fleck  _et al._ [1976] J. A. Fleck, J. R. Morris, and M. D. Feit, Time-dependent propagation of high energy laser beams through the atmosphere, Applied Physics 10, 129 (1976).
  * Yoshida [1990] H. Yoshida, Construction of higher order symplectic integrators, Physics Letters A 150, 262 (1990).
  * Suzuki [1990] M. Suzuki, Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations, Physics Letters A 146, 319 (1990).
  * Bandrauk and Shen [1992] A. D. Bandrauk and H. Shen, Higher order exponential split operator method for solving time-dependent Schrödinger equations, Canadian Journal of Chemistry 70, 555 (1992).
  * Blanes  _et al._ [2006] S. Blanes, F. Casas, and A. Murua, Symplectic splitting operator methods for the time-dependent Schrödinger equation., The Journal of Chemical Physics 124, 234105 (2006).
  * Hairer  _et al._ [2006] E. Hairer, C. Lubich, and G. Wanner,  _Geometric Numerical Integration_ (Springer, New York, 2006).
  * Lubich [2008] C. Lubich,  _From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis_ (European Mathematical Society, Zürich, 2008).
  * Kahan and Li [1997] W. Kahan and R.-C. Li, Composition constants for raising the orders of unconventional schemes for ordinary differential equations, Math. Comput. 66, 1089 (1997).
  * Riedel  _et al._ [2023] J. Riedel, P. Gelß, R. Klein, and B. Schmidt, WaveTrain: A Python package for numerical quantum mechanics of chain-like systems based on tensor trains (2023), arXiv:2302.03725 .
  * Gelß  _et al._ [2021] P. Gelß, S. Klus, M. Scherer, F. Nüske, and M. Lücke, Scikit-TT, https://github.com/PGelss/scikit_tt (2021).
  * Schmidt and Lorenz [2017] B. Schmidt and U. Lorenz, WavePacket: A Matlab package for numerical quantum dynamics. I: Closed quantum systems and discrete variable representations, Computer Physics Communications 213, 223 (2017).
  * Schmidt and Hartmann [2018] B. Schmidt and C. Hartmann, WavePacket: A Matlab package for numerical quantum dynamics.II: Open quantum systems, optimal control, and model reduction, Computer Physics Communications 228, 229 (2018).
  * Schmidt  _et al._ [2019] B. Schmidt, R. Klein, and L. Cancissu Araujo, WavePacket: A Matlab package for numerical quantum dynamics. III. Quantum‐classical simulations and surface hopping trajectories, Journal of Computational Chemistry 40, 2677 (2019).
  * Kenkre and Phatak [1984] V. Kenkre and S. Phatak, Exact probability propagators for motion with arbitrary degree of transport coherence, Physics Letters A 100, 101 (1984).
  * Schleich [2001] W. P. Schleich,  _Quantum Optics in Phase Space_ (Wiley–VCH, Berlin, 2001).
  * Cohen-Tannoudji  _et al._ [1977] C. Cohen-Tannoudji, B. Liu, and F. Laloe,  _Quantum Mechanics_ (Wiley, New York, 1977).
  * Horenko and Weiser [2003] I. Horenko and M. Weiser, Adaptive integration of molecular dynamics., Journal of Computational Chemistry 24, 1921 (2003).
  * Horenko  _et al._ [2004] I. Horenko, M. Weiser, B. Schmidt, and C. Schütte, Fully adaptive propagation of the quantum-classical Liouville equation., The Journal of Chemical Physics 120, 8913 (2004).