2021

[4]Michèle Desouter-Lecomte

1]MSME, Université Gustave Eiffel, UPEC, CNRS, Marne-La-Vallée, F-77454,
France

2]ECE - Paris, Graduate School of Engineering, Paris, 75015, France

3]Institut des Nanosciences de Paris, Sorbonne Université, CNRS, Paris, 75015,
France

[4]Institut de Chimie Physique UMR8000, Université Paris-Saclay, CNRS, Orsay,
91405, France

# Survey of the Hierarchical Equations of Motion in Tensor-Train format for
non-Markovian quantum dynamics

Etienne Mangaud<EMAIL_ADDRESS>Amine Jaouadi<EMAIL_ADDRESS>Alex Chin<EMAIL_ADDRESS>michele.desouter-lecomte@universite-paris-
saclay.fr [  [  [  *

###### Abstract

This work is a pedagogical survey about the hierarchical equations of motion
and their implementation with the tensor-train format. These equations are a
great standard in non-perturbative non-Markovian open quantum systems. They
are exact for harmonic baths in the limit of relevant truncation of the
hierarchy. We recall the link with the perturbative second order time
convolution equations also known as the Bloch-Redfield equations.  Some
theoretical tools characterizing non-Markovian dynamics such as the non-
Markovianity measures or the dynamical map are also briefly discussed in the
context of HEOM simulations. The main points of the tensor-train expansion are
illustrated in an example with a qubit interacting with a bath described by a
Lorentzian spectral density. Finally, we give three illustrative applications
in which the system-bath coupling operator is similar to that of the
analytical treatment. The first example revisits a model in which population-
to-coherence transfer via the bath creates a long-lasting coherence between
two states. The second one is devoted to the computation of stationary
absorption and emission spectra. We illustrate the link between the spectral
density and the Stokes shift in situations with and without nonadiabatic
interaction. Finally, we simulate an excitation transfer when the spectral
density is discretized by undamped modes to illustrate a situation in which
the TT formulation is more efficient than the standard one.

###### keywords:

Open quantum sytems, Non-Markovian Quantum Dynamics, HEOM, Tensor-Train,
Coherence, Linear spectroscopy

##  1 Introduction

Simulating quantum dynamics of complex systems with a large number of degrees
of freedom (DoF) remains a computational challenge. However, measured
observables often depend on a limited number of DoFs. Thus, the full system
can be described as an active subsystem embedded in an environment, which
makes fluctuate the energy levels of the subsystem. The latter is often
described quantum mechanically while the extended environment is treated by a
wide range of possibilities based on semi-classical or quantum or statistical
approaches adequately chosen with respect to the choice of system-environment
partitioning. When the number of DoFs increases, standard methods become
computationnally untractable, which is known as the "curse of dimensionality".
In this context, the low-rank tensor decomposition has aroused a constantly
growing interest. When the surrounding is modelled by an ensemble of discrete
modes, Multi Configuration Time-Dependent Hartree (MCTDH) and the multi-layer
version (ML-MCTDH) Manthe2008 ; Thoss2009 ; Meyer2011 ; Reichman2021  are
based on tensor network algorithms, mainly the Tucker and hierarchical Tucker
tensors Bader2009 ; Tobler2013 . Similarly, the expression of the time-
dependent multi-mode wave functions in Matrix Product State (MPS) also called
Tensor Train (TT) Oseledets2011 ; Oseledets2015 ; Oseledets2016 ; Orus2014
expansion has revealed its efficiency in many applications Chin2016 ;
SchroderChin2019 ; Chin2019 ; Dunnett2021 ; Dunnett2021bis ; Baiardi2019 ;
Ma2019 ; Plenio2019 ; Batista2022 ; Schmidt2022 .

In the usual approach of open quantum systems based on statistical mechanics
with a surrounding at thermal equilibrium, the environment is described less
explicitly and the active system is treated by a reduced density matrix by
tracing over the bath degrees of freedom. Path integrals derived from the
Feynman–Vernon influence functional Feynman1963  and the hierarchical
equations of motion (HEOM) Kubo1989 ; Tanimura2006 ; Tanimura_rev_20  are a
priori exact methods for harmonic baths and are closely related Yan2007 ;
Yan2009 . HEOM may also be derived from the Nakajima–Zwanzig Nakajima1958 ;
Zwanzig1960  partition of the Liouville equation by using the cumulant
expansion of the reduced propagator and properties of the Gaussian
distribution of the bath linked to the harmonic approximation
Ishizaki_Fleming2009 . Both Path Integral and HEOM formalisms have also been
recently treated by the tensor formalism in the Tensor Network Path Integral
Walters2022  and in the MPS Shi2018 ; Shi2020 ; Borrelli2021 ;
Borrelli_Gelin2021  or Tucker and hierarchical Tucker tensor format for HEOM
Shi2021 . HEOM has been applied to describe many physico-chemical processes
(see the recent review by Y. Tanimura Tanimura_rev_20 ), for instance,
excitation transfer in photosynthetical complexes Ishizaki_Fleming2009 ;
Fleming2010 ; Kreisbeck2012 ; Fleming2015 ; Zhao2015 ; Limmer2022  or in other
devices Tanimura2021 , non-adiabatic interactions and electron transfer
Tanaka_Tanimura2010 ; MangaudDMP2015 ; Thiago2016 , dynamics via conical
intersections Gelin2016 ; Duan2017 ; Dijkstra2017 ; MDL_19 ; Breuil2021 ;
Jaouadi2022 , proton tranfer BorelliTanimura2020 , laser optimal control
MDL2018 ; Jaouadi2022 , non-equilibrium fluxes Kato2016 ; Shi2017 ; Thoss2021
, and non-linear spectroscopies Tanimura_rev_20 ; Fleming2019 ; Zhu2011 ;
Tanimura2012 .

It is worth noting that accounting for temperature is handled in a similar way
in the HEOM system of equations with discrete undamped modes Shi2014  and in a
recent MPS implementation with wave functions fchem2021  in the context of
T-TEDOPA formalism Tamascelli2019  . In the latter, the wave function approach
uses discrete vibrational bath transformed into a chain and introduces finite
temperature by sampling two baths representing absorption and emission
respectively, with the latter being described by a bath of oscillators with
negative frequencies. Emission into the vacuum of these negative modes mimics
the absorption of environmental quanta that would be present in a physical
(mixed-state) thermal bath, allowing a pure wavefunction description to
capture the physics of a mixed-state initial condition without the need for
thermal sampling Tamascelli2019 ; fchem2021 . On the other hand, the
particular implementation of HEOM with undamped discrete modes also samples
baths with positive and negative frequencies Shi2014 . By comparing these
methods, we also emphasize that even if the reduced density matrix in open
quantum systems is obtained by tracing over the bath modes, this does not mean
that the information about the environment disappears. Relevant information
about the baths may be extracted, for instance, the time-dependent
distribution of the collective bath modes in each electronic state, which may
be seen as the square modulus of a dissipative wave packet Shi2012 ; Shi2014 ;
ChinChevet2019 ; MDL_19  or projection of the coherence among electronic
states along the collective modes MDL_19  and fluxes Shi2017 .

In this work, we present a pedagogical survey about HEOM and their
implementation with the TT format. We recall the main lines and in particular
the link with the perturbative second order time convolution equations also
known as the Bloch-Redfield equations, which are an illuminating step to
understand the complicated structure of HEOM. We briefly discuss how HEOM may
be used to compute some theoretical tools (non-Markovianity measure or
dynamical map) related to non-Markovian dynamics.  We then present the
principal points of the TT expansion by giving the expressions related to a
simple example where a qubit interacts with a bath. Finally, we give three
applications based on models in which the system-bath coupling operator is
similar to the one used in the survey. Finally, the appendix explains how to
encode the main expressions with a Python package.

##  2 Open quantum systems and HEOM

The standard starting point in open quantum system Breuer2002 ; Weiss2012 ;
Kuhn2011  is the partition of the DoFs of the full complex system. The active
subsystem may be only electronic DoFs like in the usual spin-boson model or
include some Brownian coordinates coupled to residual baths Chenel2014 ;
MangaudDMP2015 ; Nazir2016 ; Gelin2016 . The generic partitioning of the full
Hamiltonian in three parts is written as:

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

where ${H}_{S}$ and ${H}_{B}$ are the Hamiltonians of the active system and of
the vibrational or phonon bath(s) respectively. The system Hamiltonian may be
time dependent if it contains interaction with external fields. When the
system interacts with $N_{bath}$, the system-bath coupling is
${{H}_{SB}}=\sum\nolimits_{\alpha=1}^{{{N}_{bath}}}{{{S}_{\alpha}}}{{B}_{\alpha}}$.
${S}_{\alpha}$ and ${{B}_{\alpha}}$ are operators in the space of the system
and in the complementary space respectively.

In the electronic-nuclear partition, the system operators $S_{\alpha}$ are
$n\times n$ matrices with $n$ the number of electronic states. They act as
projectors on some states when the baths tune the electronic energies (i.e.
are diagonally coupled) or transition matrices between some of them when the
baths make fluctuate the off-diagonal electronic coupling such as in the case
of conical intersections Duan2017 ; Dijkstra2017 ; MDL_19 ; Breuil2021 ;
Jaouadi2022 . In other system-bath partition cases, the system operators are
the Brownian coordinates included in the active subspace Chenel2014 ;
MangaudDMP2015 ; Nazir2016 ; Gelin2016 . Furthermore, each bath operator
$B_{\alpha}$ is a linear combination of the position operators $q_{j}$ of the
oscillators
${{B}_{\alpha}}=\sum\nolimits_{j}^{{{N}_{\alpha}}}{c_{j}^{(\alpha)}}{{q}_{j}}$
in the discrete bath representation with $N$ oscillators. The expression of
the coupling coefficients $c_{j}^{(\alpha)}$ depends on the partition and on
the choice of the coordinates. In the following, we adopt mass-weighted
coordinates and the electronic-nuclear partition.

The initial total density operator
$\rho_{tot}(0)=\rho_{S}(0)\otimes\rho_{B,eq}$ is assumed to be factorized and
is the product of the system density operator $\rho_{S}(0)$ and a thermally
equilibrated bath density operator $\rho_{B,eq}=e^{-\beta
H_{B}}/Tr_{B}\left[e^{-\beta H_{B}}\right]$. Extension to correlated initial
conditions have been proposed in Refs. Pomyalov2010 ; Dijkstra2010 ;
Tanimura2014 ; Shi2015 .

###  2.1 Second order auxiliary operators

It is very instructive to first examine the non-Markovian perturbative
equations, which contain all the crucial tools occurring in HEOM. To do so, we
consider the simplest case with only one bath, i.e. ${H}_{SB}=SB$. The exact
formal Nakajima-Zwanzig equation given the evolution of the reduced density
matrix
$\dot{\rho_{S}}(t)=-\frac{i}{\hbar}T{{r}_{B}}\left[H,{{\rho}_{tot}}(t)\right]$
may be written as:

$\dot{\rho_{S}}(t)={{L}_{S}}\rho_{S}(t)+\int_{0}^{t}{K(t,\tau)\rho_{S}(\tau)d\tau}+I(t)$
(2)

with ${{L}_{S}}\centerdot=-i\left[{{H}_{S}},\centerdot\right]$ the system
Liouvillian and $\hbar=1$. $K(t,\tau)$ is the memory kernel which embeds the
bath influence on the system and $I(t)$ is an initial correlation term which
cancels when the system and bath can be initially factorized (i.e. $I(t)=0$
when $\rho_{tot}(0)=\rho_{S}\otimes\rho_{B,eq}$).

At the second order in the ${{H}_{SB}}$ coupling, the memory kernel becomes
(Meier1999, ):

$\displaystyle\int_{0}^{t}{d\tau}K(t,\tau){{\rho}_{S}}(\tau)=\int_{0}^{t}{d\tau}K^{(2)}(t,\tau){{\rho}_{S}}(\tau)$
$\displaystyle=i\int_{0}^{t}{d\tau}\left[S,\left\\{iC(t-\tau){{U}_{S}}(t-\tau)S{{\rho}_{S}}(\tau)\right.\right.$
$\displaystyle\left.\left.\times
U_{S}^{\dagger}(t-\tau)\right\\}+\left\\{hc\right\\}\right]$ (3)

where $C(t)$ is the correlation function of the collective bath mode
$C(t)={{\left\langle B(t)B(0)\right\rangle}_{eq}}$, $B(t)$ is the Heisenberg
representation of the operator with Hamiltonian $H_{B}$ and
${{\left\langle\centerdot\right\rangle}_{eq}}$ denotes the average over a
Boltzmann distribution at a given temperature $T$. $U_{S}(t)=e^{-iH_{S}t}$ is
the propagator of the system. The treatment up to the fourth order may be
found in Ref. Cao2002  and the extension by the generalized master equation
method in Ref. Geva2019 .

The quantum bath correlation function is the main descriptor of the bath
Mukamel1995 . The usual step to go towards the second order auxiliary density
operator (ADO) or HEOM is the representation of the correlation function as a
sum of damped decaying functions:

$C\left(t-\tau\right)={{\left\langle
B(t)B(\tau)\right\rangle}_{eq}}=\underset{k=1}{\overset{K}{\mathop{\sum}}}\,{{\alpha}_{k}}{{e}^{i{{\gamma}_{k}}\left(t-\tau\right)}}$
(4)

where ${{\alpha}_{k}}$ and ${{\gamma}_{k}}$ are complex parameters. The sum is
a priori infinite but truncated to $K$ modes, which have been called the
"artificial bath modes" in Ref. Meier1999 . Some extensions to different
analytical forms or arbitrary correlation functions have been proposed
recently in Refs. Kleinekat2019 ; Yan2019 ; Ikeda2020 ; Yan2022 ; Lambert2020
.

By inserting Eq.(4) in Eqs.(2) and (3), each artificial mode corresponds to a
particular memory integral of the time dependent integro-differential
equation. Each integral is set equal to an ADO :

$\displaystyle{\rho}_{k}(t)/(i{\alpha}_{k})$
$\displaystyle=\int_{0}^{t}{d\tau}\left[S,\left\\{i{{e}^{-{{\gamma}_{k}}(t-\tau)}}{{U}_{S}}(t-\tau)S{{\rho}_{S}}(\tau)\right.\right.$
$\displaystyle\left.\left.\times
U_{S}^{\dagger}(t-\tau)\right\\}+\left\\{hc\right\\}\right].$ (5)

The dimension of ${{\rho}_{k}}(t)$ is that of the $S$ operator and thus of
${{H}_{S}}$. A time local system of coupled equations may be obtained by
taking the first derivative Meier1999 ; Kleinekat2004 ; Pomyalov2010 .
Different choices are possible to define the ADOs. From Eq.(5), we obtain the
operational equations Pomyalov2010 :

$\displaystyle{{\dot{\rho}}_{S}}(t)={{L}_{S}}{{\rho}_{S}}(t)+i\underset{k=1}{\overset{K}{\mathop{\sum}}}\,{\left[S,{{\rho}_{k}}(t)\right]}$
$\displaystyle{{\dot{\rho}}_{k}}(t)=\left(i{{\gamma}_{k}}+{{L}_{S}}\right){{\rho}_{k}}(t)+i\left[{{\alpha}_{k}}S{{\rho}_{S}}(t)\right.$
$\displaystyle\left.-{{{\tilde{\alpha}}}_{k}}{{\rho}_{S}}(t)S\right]$ (6)

where ${{L}_{S}}$ is the system Liouvillian. The ${{\tilde{\alpha}}_{k}}$
parameters will be discussed below.

Each ADO being associated only to one decay mode, it may be considered as
resulting from a single excitation in this mode and denoted by an array of $K$
indexes with one in the $k^{th}$ position and zero everywhere else as shown in
figure 1. These ADOs will constitute the first level of the HEOM hierarchy.

Figure 1: Schematic representation of the auxiliary operators associated to
each artificial bath mode $k=1,K$ and involved in the second order master
equation (Eq.( 6). Each ADO corresponds to a single excitation in the mode.

A real classical correlation function may be obtained by molecular dynamics
Valleau2012 ; Fleming2015 ; Thiago2016 ; Dunnett2021bis  or directly from
experimental results (Marcus2002, ) and corrected to get the complex quantum
correlation function satisfying the fluctuation-dissipation theorem
(Mukamel1995, ). Direct parametrization of $C(t)$ fitted by different
approaches Lambert2020 , Prony method Yan2022  or expansion on Chebyshev or
Bessel functions Kleinekat2019  have also been proposed. However, the relation
with the bath spectral density, defined as follows:

$C(t)=\int_{-\infty}^{\infty}{d\omega{{J}}\left(\omega\right){{\left({{e}^{\beta\omega}}-1\right)}^{-1}}}{{e}^{-i\omega
t}}$ (7)

where $\beta=1/{{k}_{B}}T$ is the Boltzmann factor, is currently used and the
$\left\\{{{\alpha}_{k}},{{\gamma}_{k}}\right\\}$ parameters are obtained from
the parametrization of
$C(\omega)={{J}}\left(\omega\right){{\left({{e}^{\beta\omega}}-1\right)}^{-1}}$
where ${{J}}\left(\omega\right)$ is independent of the temperature and
${{\left({{e}^{\beta\omega}}-1\right)}^{-1}}$ is the Bose function related to
the quantum fluctuation-dissipation theorem. In the discrete case, the
spectral density is defined from the system-bath coupling coefficients
${{c}_{j}}$ related to $N_{disc}$ discrete modes

$J(\omega)=\frac{\pi}{2}\underset{j=1}{\overset{N_{disc}}{\mathop{\sum}}}\,{\frac{c_{j}^{2}}{{{\omega}_{j}}}}\delta(\omega-{{\omega}_{j}}).$
(8)

In the continuous representation, the spectral density is often approximated
by an Ohmic function ${{J}_{Ohm}}\approx c\omega{{f}_{cutoff}}(\omega)$ or a
super Ohmic function ${{J}_{SuperOhm}}\approx
c{{\omega}^{3}}{{f}_{cutoff}}(\omega)$, with $c$ a constant and $f_{cutoff}$
an exponential or Lorentzian cutoff. The Ohmic function is used for solvents
while the super Ohmic one is relevant for the solid phase and phonon baths. In
order to get a parametrization of the bath correlation function (Eq.(4))
through Eq.(7), one has to perform a fitting procedure of the spectral
functions with relevant analytical functions. A special care in fitting the
spectral density low frequency behavior is often necessary to get accurate
results (Eisfeld2014, ). Here, we will discuss only the cases of Ohmic and
super-Ohmic Lorentzian functions. We adopt the Tannor-Meier parametrization
Meier1999  in which the spectral density is fitted by  ${n}_{lor}$ two-pole
Lorentzian functions depending on three parameters
$\left({{p}_{l}},{{\Omega}_{l}},{{\Gamma}_{l}}\right)$ where ${{p}_{l}}$ is
the associated weight, ${\Omega}_{l}$ the central frequency, and
${\Gamma}_{l}$ the bandwidth.

${{J}_{Ohm}}\left(\omega\right)=\sum\limits_{l=1}^{{{n}_{lor}}}{\frac{{{p}_{l}}\omega}{\Upsilon\left({{\Omega}_{l}},{{\Gamma}_{l}}\right)}}$
(9)

with

${\Upsilon\left({{\Omega}_{l}},{{\Gamma}_{l}}\right)=\left[{{\left(\omega+{{\Omega}_{l}}\right)}^{2}}+\Gamma_{l}^{2}\right]\left[{{\left(\omega-{{\Omega}_{l}}\right)}^{2}}+\Gamma_{l}^{2}\right]}$

All the $\left\\{{{\alpha}_{k}},{{\gamma}_{k}}\right\\}$ parameters may be
obtained from the integral (7) after substitution of Eq.(9). Each Lorentzian
corresponds to two artificial decay modes. The Bose function generates an
infinite series of terms coming from its poles. They are known as the
Matsubara terms. In practice, the number of Matsubara terms is very small at
high temperature but may become numerous at low temperature rendering the ADO
method computationally more demanding. In the low temperature regime, Padé
approximants of the Bose function (Yan2011, ),  fitting procedure of $C(t)$
used to capture the Matsubara terms Lambert2019  or a recent correction scheme
Fay2022  can be used. The analytical expressions of
$\left\\{{{\alpha}_{k}},{{\gamma}_{k}}\right\\}$ (Eq.(4)) as functions of the
$\left({{p}_{l}},{{\Omega}_{l}},{{\Gamma}_{l}}\right)$ parameters of
$J_{Ohm}(\omega)$ and those related to the Matsubara terms are given in Refs.
(Pomyalov2010, ) or (Chenel2014, ). One also could find in these references
the expression of the ${{\tilde{\alpha}}_{k}}$ parameters (see Eq.(6)). They
come from a particular expression of the complex conjugate of $C(t)$ that may
be written as:
${{C}^{*}}\left(t\right)=\underset{k=1}{\overset{K}{\mathop{\sum}}}\,{{\tilde{\alpha}}_{k}}{{e}^{i{{\gamma}_{k}}t}}$
with the same ${{\gamma}_{k}}$ as in Eq.(4). In the super Ohmic case, the
fitting functions have four poles leading to four artificial decay channels

${{J}_{SuperOhm}}\left(\omega\right)=\sum\limits_{l=1}^{{{n}_{lor}}}{\frac{{{p}_{l}}{{\omega}^{3}}}{\Upsilon\left({{\Omega}_{l,1}},{{\Gamma}_{l,1}}\right)\Upsilon\left({{\Omega}_{l,2}},{{\Gamma}_{l,2}}\right)}}.$
(10)

The analytical expressions to get the
$\left\\{{{\alpha}_{k}},{{\gamma}_{k}}\right\\}$ parameters of $C(t)$ (Eq.(4))
from
$\left({{p}_{l}},{{\Omega}_{l,1}},{{\Gamma}_{l,1}},{{\Omega}_{l,2}},{{\Gamma}_{l,2}}\right)$
parameters of $J_{SuperOhm}(\omega)$ are gathered in Ref. MangaudMeier2017 .

###  2.2 HEOM

First, we discuss the case with a single bath, i.e. we consider only one $S$
operator. The generalization will be discussed at the end of the section. When
the bath correlation time becomes long with respect to the system
characteristic timescale due to a strongly peaked spectral density, the high
non-Markovianity generally involves a non-perturbative regime. HEOM are one of
the reference dynamical methods for open quantum systems modelled with a
harmonic bath. The equations giving the evolution of the reduced density
matrix were originally derived for a Drude-Lorentz spectral density in the
high-temperature limit from the Kubo stochastic Liouville equation
(kubo_stochastic_1963, ; Kubo1989, ; Tanimura2006, ; Tanimura_rev_20, ) and
the Feynman-Vernon influence functional formalism Yan2007 ; Yan2009 . Another
approach to derive these equations is by referring to the remarkable property
that the cumulant expansion limited to second order is exact when the bath
statistics are Gaussian (Ishizaki_2005, ; Ishizaki_Fleming2009, ). This is
based on the Wick theorem Wick1950 . The reduced density matrix
${{\rho}_{S,I}}(t)$ for the system in interaction representation
$\rho_{S,I}(t)=e^{iH_{0}t}\rho_{S}(t)e^{-iH_{0}t}$ with $H_{0}=H_{S}+H_{B}$,
is given by the partial trace over the bath of the time evolution of the total
density matrix :

${{\rho}_{S,I}}(t)=T{{r}_{B}}\left[{\mathcal{T}^{(+)}}{{e}^{\int\limits_{0}^{t}{d\tau{{L}_{SB,I}}(\tau)}}}\rho_{B}^{eq}\right]{{\rho}_{S,I}}(0)$
(11)

where $\mathcal{T}^{(+)}$ is a time-ordering operator and a factorization is
assumed at the initial time $t=0$. ${{L}_{SB,I}}(t)$ is the Liouvillian in
interaction representation,
${{L}_{SB,I}}(t)\centerdot=-i[{{S}}(t){{B}}(t),\centerdot]$, where
$S(t)={{e}^{i{{H}_{S}}t}}S{{e}^{-i{{H}_{S}}t}}$ and
$B(t)={{e}^{i{{H}_{B}}t}}S{{e}^{-i{{H}_{B}}t}}$.

At the second order in the cumulant expansion, Eq.(11) becomes

${{\rho}_{S,I}}(t)={\mathcal{T}^{(+)}}{{e}^{\int\limits_{0}^{t}{d\tau\mathcal{K}_{I}(\tau)}}}{{\rho}_{S,I}}(0).$
(12)

$\mathcal{K}_{I}(\tau)=\int\limits_{0}^{\tau}dt^{\prime}K_{I}^{(2)}(\tau,t^{\prime})$
corresponds to the second-order memory term occurring in the second-order
perturbation theory (Eq.(3)) but here, the second order is an exact expression
for the cumulant expansion :

$\displaystyle\int_{0}^{\tau}dt^{\prime}K_{I}^{(2)}(\tau,t^{\prime})\centerdot$
$\displaystyle=-[S(\tau),\int\limits_{0}^{\tau}dt^{\prime}C(\tau-t^{\prime})S(t^{\prime})\centerdot-\\{h.c.\\}].$
(13)

By inserting the correlation function parametrization (Eq.(4)), equation (13)
is separable into a sum of $K$ operators $\mathcal{K}_{Ik}$ where $K$ is the
number of artificial decay modes. In addition, each $\mathcal{K}_{Ik}$ is
decomposed as :

$\mathcal{K}_{Ik}(\tau)=\Phi_{I}(\tau)W_{Ik}(\tau)$ (14)

where

$\Phi_{I}(\tau)\centerdot=-[S(\tau),\centerdot]$ (15)

and

$W_{Ik}(\tau)=\int\limits_{0}^{\tau}dt^{\prime}e^{i{{\gamma}_{k}}\left(\tau-t^{\prime}\right)}{{\Theta}_{Ik}}\left(t^{\prime}\right)$
(16)

with

$\Theta_{Ik}(t^{\prime})\centerdot={{\alpha}_{k}}S\left(t^{\prime}\right)\centerdot-{{{\tilde{\alpha}}_{k}}}\centerdot
S\left(t^{\prime}\right).$ (17)

The exact solution takes the following form:

$\displaystyle{{\rho}_{S,I}}(t)={\mathcal{T}^{(+)}}\prod\nolimits_{k}^{K}{{{e}^{\int_{0}^{t}{d\tau}{\mathcal{K}_{Ik}}(\tau)}}}{{\rho}_{S,I}}(0)$
$\displaystyle={\mathcal{T}^{(+)}}\prod\nolimits_{k}^{K}{{{e}^{\int_{0}^{t}{d\tau}\int_{0}^{\tau}{dt^{\prime}\Phi_{I}(t^{\prime}){{e}^{i{{\gamma}_{k}}(\tau-t^{\prime})}}{{\Theta}_{Ik}}(t^{\prime})}}}}$
$\displaystyle\times{{\rho}_{S,I}}(0).$ (18)

The master equation can then be derived as

${{\dot{\rho}}_{S,I}}(t)={{\mathcal{T}}^{(+)}}\sum\nolimits_{k}^{K}{{{\mathcal{K}}_{Ik}}(t){{\rho}_{S,I}}(t)}.$
(19)

This equation is time non-local since $\mathcal{K}_{Ik}(t)$ contains an
integral. A time local system of coupled equations is obtained by defining the
ADOs by the same strategy as in the second order approach :

$\displaystyle{{\rho}_{\mathbf{n},I}}(t)$
$\displaystyle={\mathcal{T}^{(+)}}\prod\nolimits_{k}^{K}{{{W}_{Ik}}}{{(t)}^{{{n}_{k}}}}{{e}^{\int_{0}^{t}{d\tau{\mathcal{K}_{Ik}}(\tau)}}}{{\rho}_{S,I}}(0)$
(20)

where

$\mathbf{n}=\left\\{{{n}_{1}},...,{{n}_{k}},...{{n}_{K}}\right\\}$ (21)

in a vector of nonnegative integers giving the occupation number in each
artificial decay mode. The case $\mathbf{n}=\left\\{0,...,0,...0\right\\}$
corresponds to $\rho_{S,I}$. Time local coupled equations among the ADOs may
be derived by working with the Fourier-Laplace transforms of Eqs. (18) and
(20) and using integration by parts Fleming2010 ; MDL_19 . One recovers the
relations established in Ref.Kubo1989 . After the inverse Laplace-Fourier
transform and the return in the Schrödinger representation, the HEOM read :

$\displaystyle{{{\dot{\rho}}}_{\mathbf{n}}}(t)=$
$\displaystyle{{L}_{S}}{{\rho}_{\mathbf{n}}}(t)+i\sum\limits_{k=1}^{K}{{{n}_{k}}{{\gamma}_{k}}{{\rho}_{\mathbf{n}}}}(t)$
$\displaystyle-i\left[S,\sum\limits_{k=1}^{K}{{{\rho}_{\mathbf{n}_{k}^{+}}}(t)}\right]$
$\displaystyle-i\sum\limits_{k=1}^{K}{{{n}_{k}}\left({{\alpha}_{k}}{{S}}{{\rho}_{\mathbf{n}_{k}^{-}}}(t)-{{{\tilde{\alpha}}}_{k}}{{\rho}_{\mathbf{n}_{k}^{-}}}(t)S\right)}$
(22)

or more shortly by using definitions (15) and (17) in Schrödinger
representation (i.e., $\Phi(\tau)\centerdot=-[S,\centerdot]$ and
$\Theta_{k}\centerdot={{\alpha}_{k}}S\centerdot-{{{\tilde{\alpha}}_{k}}}\centerdot
S$):

$\displaystyle{{\dot{\rho}}_{\mathbf{n}}(t)}={{L}_{S}}{{\rho}_{\mathbf{n}}}$
$\displaystyle+i\sum\limits_{k=1}^{K}{\left({{n}_{k}}{{\gamma}_{k}}{{\rho}_{\mathbf{n}}}-\Phi{{\rho}_{\mathbf{n}_{k}^{+}}}-{{n}_{k}}{{\Theta}_{k}}{{\rho}_{\mathbf{n}_{k}^{-}}}\right)}.$
(23)

The subscripts

$\mathbf{n}_{k}^{+}=\left\\{{{n}_{1}},\cdots,{{n}_{k}}+1,\ldots,{{n}_{{{n}_{K}}}}\right\\}$

and

$\mathbf{n}_{k}^{-}=\left\\{{{n}_{1}},\cdots,{{n}_{k}}-1,\ldots,{{n}_{{{n}_{K}}}}\right\\}$

denote the matrices for which one occupation number differs by one unit in the
hierarchy ${{n}_{k}}\to{{n}_{k}}\pm 1$. The sum of the occupation numbers
defines the level of the hierarchy $L=\sum n_{k}$. The total number of
matrices when the hierarchy is limited at level $L$ is $(L+K)!/L!K!$. Each
matrix is connected to the matrices of the lower or the upper levels. Figure 2
illustrates the case with two artificial decay channels (i.e. one Ohmic
Lorentzian in the spectral density and no Matsubara term). Initially,
$\rho_{00...0}=\rho_{S}$ and all the ADOs are zero. They represent the initial
bath at equilibrium. When the final relaxed state is reached, the converged
ADOs represent the new equilibrium and may be taken as the initial condition
to describe another process affecting the equilibrated system. For instance to
describe the stationary fluorescence. All the terms occurring in the
derivative of each ADO are schematized in figure 2 for a simple example with
two bath modes.

Figure 2: Schematic representation of the auxiliary operators of the HEOM
hierarchy up to level $L=2$ for the case with two artificial bath modes $K=2$.
$\Phi$ acts on matrices for which one occupation number increases by one unit.
$\Theta_{k}$ is applied on matrices for which one occupation number decreases
by one unit. The sum of the decay rates $n_{k}\gamma_{k}$ multiplying each
matrix is indicated. The schematic representation of the terms involved in the
derivative of one matrix is complete for $\rho_{00}$ (orange); $\rho_{10}$
(blue) and $\rho_{01}$ (green).

The generalization is straightforward when there are $N_{bath}$ uncorrelated
baths, each associated with $K_{b}$ artificial bath modes. Then,
$K=\sum\nolimits_{b=1}^{{{N}_{bath}}}{{{K}_{b}}}$. Each bath is linked to the
system by an operator $S_{b}$ and the equations become :

$\displaystyle{{\dot{\rho}}_{\mathbf{n}}(t)}={{L}_{S}}{{\rho}_{\mathbf{n}}}+i\sum\limits_{b}^{{{N}_{bath}}}{\sum\limits_{k=1}^{{{K}_{b}}}{\left({{n}_{b,k}}{{\gamma}_{b,k}}{{\rho}_{\mathbf{n}}}\right.}}$
$\displaystyle{{\left.-{{\Phi}_{b}}{{\rho}_{\mathbf{n}_{k}^{+}}}-{{n}_{b,k}}{{\Theta}_{b,k}}{{\rho}_{\mathbf{n}_{k}^{-}}}\right)}}.$
(24)

An alternative to the continuous representation of the spectral density with
$K$ artificial modes is the discretization with $N_{disc}$ undamped modes in
the spirit of MCTDH or ML-MCTDH computations Burghardt2012 ; Burghardt2020 ;
Burghardt2021 ; kuhn2016  or MPS simulations in the Hilbert space
SchroderChin2019 ; Chin2019 ; Dunnett2021 ; Dunnett2021bis . Discretization
has also been illustrated with HEOM Shi2014  and checked in TT format Shi2018
; Shi2021 . The expression of the corresponding HEOM equations are given in
Refs. Shi2014 ; Shi2018 . The active system is then coupled to two identical
baths with positive or negative frequencies describing emission and absorption
of energy. Since there are two baths, the number of terms in the correlation
function is large $K=2N_{disc}$ but there are no Matsubara terms associated
with the poles of the Bose function in the continuous case. The discrete
couplings are smaller than those of the artificial modes. We will illustrate
in one application that the TT implementation may be more efficient than the
standard one in that case.

The system of the coupled differential equations (24) or those related to
discrete modes given in refs. Shi2014 ; Shi2018  can then be solved using
standard numerical methods such as Runge-Kutta 4 (RK4), Cash-Karp (RK4-5)
adaptative step-size or Arnoldi algorithms. Although HEOM are an infinite
system of differential equations, they are truncated at a maximum value of
$n_{k}=n_{max}$. Thus, several simulations with increasing values of $n_{max}$
should be carried out until the results (density matrix or observables) are
properly converged. For comparison between the standard and TT formulations,
we have used a home-made fortran code parallelized with OpenMP. Other
implementations in python or on new architectures such as GPU Kreisbeck2014
could be useful. See a review of different softwares in Ref. Lambert2020 .

###  2.3 HEOM and non-Markovianity

As mentioned above, HEOM are a standard method to tackle non-perturbative and
non-Markovian regimes. Non-Markovianity obviously depends on the partition
since it is roughly predicted by the characteristic timescales of the active
system and the corresponding bath. These typical dynamical times may be very
different with respect to the system definition. Inserting some effective bath
coordinates into the system part (to reduce the coupling towards the residual
bath and modify the corresponding timescales) may be an efficient strategy to
change the Markovianity regime Burghardt2010 ; Chenel2014 ; Nazir2014 ;
MangaudDMP2015 ; Nazir2016 ; MDL_19 . When the partition leads to a non-
Markovian regime, an abundant literature has been devoted to the
characterization of the non-Markovianity by different measures and from a
fundamental perspective, to a mathematical definition of non-Markovian quantum
dynamical maps, which is still an open problem. These fundamental questions
are reviewed for instance in Refs. RevBreuer2016 ; deVega2017 . A detailed
presentation of these items is beyond the scope of this paper and we summarize
here only the main points. A Markovian behavior is linked to a continuous loss
of information from the open system to the surrounding while a flow from the
environment back to the open system is the signature of a non-Markovian
effect. The return of information from the bath could modify the system
dynamics. The measures aim at quantifying this backward flow. One may cite,
among others, the trace distance measure Laine2009 ; RevBreuer2016 , the
entanglement measure Rivas2010 , negativity of time-dependent canonical rates
Anderson2014 , and the Bloch volume measure Paternostro2013 . They are
compared for instance in Refs. Anderson2014 ; Yi2015 . In the last case, the
measure is based on an estimation of the volume of the accessible states in
the generalized Bloch sphere for a $n-$state system. A non-monotonic decrease
of this volume is a signature of non-Markovianity. Any density matrix may be
expanded in the basis set of ${{n}^{2}}$ operators : the normalized identity
$G_{0}=I/\sqrt{n}$ and $G_{i}(i=1,...,n^{2}-1)$: the $n^{2}-1$ generators of
$SU(n)$ Kimura2003 ; Aerts2014 . They are the Pauli matrices for $n=2$ and the
Gell-Mann matrices for $n=3$. The volume of accessible states
$V(t)=\det(\mathbf{F}(t))$ is obtained from the determinant of the matrix
${{F}_{m,n}}(t)=Tr\left({{G}_{m}(0)}{{G}_{n}}(t)\right)$. This requires
$n^{2}$ propagations of the basis operators and this is easily obtained with
HEOM Atabek2017 ; Chevet2018 ; MDL2018 ; MangaudMeier2017 ; ChinChevet2019 ;
MDL_19 . Similarly, the $n^{2}$ basis operators may be used to build the
$n^{2}\times n^{2}$ decoherence matrix Anderson2014 ; MangaudMeier2017
${{D}_{ij}}(t)=\sum\nolimits_{m=0}^{{{n}^{2}}-1}{Tr\left[{{G}_{m}}{{G}_{i}}{{\Lambda}_{t}}\left[{{G}_{m}(t)}\right]{{G}_{j}}\right]}$
where ${{\Lambda}_{t}}\left[{{G}_{m}(t)}\right]$ is the right member of the
master equation ${{\dot{G}}_{m}}(t)={{\Lambda}_{t}}\left[{{G}_{m}(t)}\right]$,
in other words, the expression of ${{\dot{\rho}}_{00..0}}(t)$ when the initial
state is $G_{m}$ (see Eq.(22)). The eigenvalues of this decoherence matrix are
called the canonical decay rates ${{\gamma}_{k}}(t)$ Anderson2014 ;
MangaudMeier2017 . They witness non-Markovianity when some of them become
negative signalizing an information return towards the system.

The concept of dynamical map in the theory of open quantum is relevant when
the initial state of the total system is a factorized product state but it
remains a debated point for initial system-bath entangled state RevBreuer2016
; Aspuru2010 ; Rosario2013 . While the master equation is related to the time
derivative ${{\dot{\rho}}}(t)={{\Lambda}_{t}}\left[{{\rho}(t)}\right]$, the
corresponding dynamical map $\phi_{t}$ transformed any initial state $\rho(0)$
to $\rho(t)$, i.e., $\rho(t)={{\phi}_{t}}\left[\rho(0)\right]$. The map is
expected to be completely positive or at least positive to ensure that
$\phi_{t}$ maps physical states to physical states. This preserves the
Hermiticity and the trace of operators. The mathematical properties of the map
are discussed mainly in the quantum information community. By expanding both
$\phi_{t}\left[\rho(0)\right]$ and $\rho(0)$ in the basis of the $n^{2}$
operators $G_{k}$, the map may be expressed in matrix notation
$\rho(t)=\sum\nolimits_{k,l=0}^{{{n}^{2}}-1}{{{F}_{lk}}}(t)Tr\left[{{G}_{l}}\rho(0)\right]{{G}_{k}}$,
i.e., as a function of the $\mathbf{F}(t)$ matrix discussed above and easily
computed by HEOM. In order to analyze the positivity, it is rather the
$n^{2}\times n^{2}$ Choi matrix choi1975 ; Kye2022 ; Anderson2014  that is
used. It corresponds to the expansion of the map in the basis set of the
$n^{2}$ projector matrices $\lvert j\rangle\langle k\rvert$ related to a basis
of the system Hilbert space in place of the generators $G_{n}$. The map is
completely positive if and only if the eigenvalues of the Choi matrix are
positive. This Choi matrix would be easily obtained by HEOM by a procedure
similar to that providing the $\mathbf{F}(t)$ matrix via the propagation of
the $n^{2}$ projector matrices $\lvert j\rangle\langle k\rvert$. To our
knowledge, this systematic analysis has not been carried out with HEOM and
could be interesting. In our previous works, we have mainly computed with HEOM
$\mathbf{F}(t)$ Chevet2018 ; Atabek2017  and the canonical rates
$\gamma_{k}(t)$ MangaudMeier2017 ; MDL2018 . Anyhow, we did not detect a loss
of conservation of the trace of $\rho_{S}(t)$ even when dynamics is non-
Markovian as shown by the volume of accessible states or by the canonical
rates. It is well known that numerical instabilities with non-conservation of
the norm may occur mainly for long time in the standard formulation of HEOM
Reichman2019  or in the TT approach due to the variational approach Shi2020
and limitation of the tensor ranks Borrelli2021 .

The control of open quantum systems Koch2016  has aroused a renewed interest
mainly in the context of quantum technologies that rely on the coherent
manipulation and transfer of information, encoded in quantum states
Maniscalco2013 ; Mirkin2019 ; Sugny2022qt . Non-Markovianity is expected to be
a resource to improve the control by exploiting the transitory flow back. The
HEOM formalism has proven its efficiency when it is coupled with different
control strategies, in particular with optimal control protocols MDL2018 ;
Jaouadi2022 .

##  3 HEOM in Tensor-Train format

In this section, we summarize the main relations of the tensor train formalism
(also named matrix product state (MPS)) that is an interesting way for
representing a high-dimensional tensor as the one we are using for HEOM. The
idea of this TT format is to decompose the tensor into a network of low-
dimensional tensors called TT cores coupled in a chain. MPS has received a
growing interest in the quantum physics community Chin2016 ; SchroderChin2019
; Chin2019 ; Dunnett2021 ; Dunnett2021bis ; Baiardi2019 ; Ma2019 ; Plenio2019
; Batista2022 ; Schmidt2022 . The application in HEOM was already suggested by
Shi Shi2018 ; Shi2020  who also uses the Tucker representation Shi2021  and
later by Borelli Borrelli2021 ; Borrelli_Gelin2021 . We present a pedagogical
survey showing the way the super-operators involved in the TT-HEOM formalism
are written and by expliciting them in a simple case of a two-level system
coupled to a bath with a spectral density fitted by an Ohmic Lorentzian
(Eq.(9)) leading to two artificial decay modes. The Appendix gathers the main
steps for the encoding in python.

###  3.1 Representation of the ADOs

When the system is a $n$-state case, each ADO is a $n\times n$ matrix that may
be reshaped in a super-vector with $n^{2}$ elements
$\bar{\rho}^{\alpha}_{\textbf{n}}$ where $\alpha\in[1,n^{2}]$ stands for a
($a,b$) element of the system density matrix ($a,b\in[1,n]$). The global index
$\mathbf{n}$ corresponds as in Eq.(21) to the occupation number in each decay
mode. When the TT format is adopted, each occupation number $n_{k}$ runs from
0 to $n_{max}$. The total number of matrices is then larger than when the
hierarchy is truncated at a given level $L$ in the standard formulation. Each
element of the high-dimensional array $\mathbf{\bar{\rho}}$ is written in TT-
format as :

$\displaystyle\bar{\rho}^{\alpha}_{\textbf{n}}\approx\sum_{j_{0}}\sum_{j_{1}}\cdots\sum_{j_{k}}\cdots\sum_{j_{K+1}}A_{0}(j_{0},\alpha,j_{1})$
$\displaystyle\times A_{1}(j_{1},n_{1},j_{2})\cdots
A_{k}(j_{k},n_{k},j_{k+1})$ $\displaystyle\times\cdots
A_{K}(j_{K},n_{K},j_{K+1}).$ (25)

The summation index $j_{k}$ goes from 1 to $r_{k}$ where $r_{k}$ is the rank
also called the bond ($r_{0}=r_{K+1}=1$ for dimensionality consistency). $K$
is the number of decay modes. $A_{k}$ are the cores, i.e. arrays of dimension
$r_{k}\times n_{dim}\times r_{k+1}$ where $n_{dim}=n^{2}$ for $k=0$ and
$n_{dim}=n_{HEOM}=n_{max}+1$ for $k\neq 0$ ($n_{HEOM}=L+1$ where $L$ is the
hierarchy level). The choice of the rank $r_{k}$ is crucial and convergence
must be carefully checked. Complex baths (and realistic ones) often need many
decay modes and high hierarchy level: $K$ and $L$ increase dramatically which
leads to heavy simulations.

TT decomposition allows a priori to deal with this dimensionality curse.
Indeed, instead of using an array with $n^{2}\times{n_{HEOM}}^{K}$ in a
standard approach, the tensor train decomposition stores $n^{2}r_{1}\times
n_{HEOM}\prod_{k=1}^{K}r_{k}r_{k+1}$ tensor elements. For instance, if the
ranks for every core are all the same ($\forall k\in[1,K],r_{k}=r$), the
number of elements is $r(n^{2}+n_{HEOM})+Kr^{2}n_{HEOM}$ in the TT approach.
With the TT formalism, the increase with the number of modes is linear with
$K$. This might allow to deal with large values of $K$ and thus extend the
capabilities of the HEOM method to describe more complex environments.
However, the TT decomposition approximates the real tensor and is exact only
if the ranks $r_{k}$ grow up to infinity. In practise, they are truncated to a
sufficiently high value to ensure the computation convergence. Choosing the
best ranks is not a trivial task and goes beyond the scope of this article. We
discuss this point in section 3.3.

Tensor representations are schematized in figure 3.

Figure 3: Upper panel : Schematic representation of a tensor train with the
HEOM formalism for a two-level system ($n=2$) with a bath modelled by two
single Tannor-Meier Lorentzian functions with no Matsubara frequencies ($K=4$)
at a hierarchy level of $n_{HEOM}=3$. Individual matrices or vectors are
called cores of the tensor (and often denoted as $A_{k}$). Arbitrary values
have been chosen here for the tensor ranks: $r_{1}=4$, $r_{2}=2$, $r_{3}=3$
and $r_{4}=4$. For instance, the tensor element $\rho_{2,3,2,2,1}$, which
corresponds to $\rho^{12}_{3,2,2,1}$, i.e. the element 1,2 of the ADO for
occupation number $\mathbf{n}=3,2,2,1$, is computed by performing vector-
matrix, matrix-matrix or matrix-vector products with the blue cores on the
figure. As the initial and final tensor ranks ($r_{0}=r_{K+1}=1$) are equal to
one, the final result is a scalar number. Lower panel : Tensor train
representation of the ADOs. $A_{k}$ are the cores of the tensor. The circles
represent the physical legs. $\alpha$ runs from 1 to $n^{2}$ where $n$ is the
number of states in the system and $n_{k}$ is the index (occupation number)
for each decay mode that runs from 0 to $n_{max}$. The rectangles are matrices
$r_{k-1}\times r_{k}$ where the rank $r_{k}$ is also called the bond.

###  3.2 HEOM super-Liouvillian

The HEOM equations (23) and (24) contain a Liouvillian operator related to the
Hamitonian of the system, a damping term, a term coupling to matrices of a
higher level in the hierarchy and a term coupling to a lower level. In super-
operator notation, one has:

$\mathcal{L}=\mathcal{L}_{S}+\sum_{k^{\prime}=1}^{K}\left(\mathcal{L}_{k^{\prime}}+\mathcal{L}_{k^{\prime}+}+\mathcal{L}_{k^{\prime}-}\right)$
(26)

where, $k^{\prime}=(k,b)$ is a collective index which adresses both the index
of the correlation function terms ($k\in\left[1,n_{cor,b}\right]$) and the
bath $b$.

For pedagogical purpose, in this section, we give the expression of the
different contributions to the super-operator of the forward propagation with
TT for a simple example. We consider a two-level system ($n=2$). The excited
state is coupled to a single bath that tunes the energies and is coupled
diagonally to the system. The spectral density is a two-pole Lorentzian
leading to two decay modes ($K=2$). We assume a high temperature so there is
no Matsubara term. The corresponding coupling operator is

${{S}_{k}}=\left(\begin{matrix}0&0\\\ 0&1\\\ \end{matrix}\right)$ (27)

for decay mode $k=1,2$. HEOM is treated at order two $n_{HEOM}=2$, i.e., for
each mode, the occupation number $n_{k}$ may take the value 0 or 1
($n_{max}=1$). The hierarchy contains 4 $(2\times 2)$ matrices related to
occupation numbers 00, 01, 10 and 11. The super-operator is then a $(16\times
16)$ matrix in this example.

All the elements of all the ADOs ($n\times n$ matrices) are reshaped in a
supervector that contains $n^{2}$ groups corresponding to a given element of
the system density matrix $\alpha=(a,b)$ when the ADO indexes of
$\left\\{{{n}_{1}},..,{{n}_{k}},.,{{n}_{K}}\right\\}$ run from $0$ to
$n_{max}$ beginning by the last one. In the example with a two-level system,
two decay modes ($K=2$) and $n_{max}=1$, there are 4 groups labelled 00, 01,
10 and 11 (see figure 4.

####  3.2.1 Matrix $\mathcal{L}_{S}$

The expression of the $n^{2}\times n^{2}$ super-operator for the system
Liouvillian, which must act on the system density matrix and on each auxiliary
density operator (ADO) is straightforward ($\otimes$ denotes here the
Kronecker product):

$\displaystyle{\mathcal{L}_{S(ADO)}}=-i(H\otimes{{I}_{n}}-{{I}_{n}}\otimes H)$
$\displaystyle=-i\left(\begin{matrix}0&-{{H}_{12}}&{{H}_{12}}&0\\\
-{{H}_{21}}&{{H}_{11}}-{{H}_{22}}&0&{{H}_{12}}\\\
{{H}_{21}}&0&{{H}_{22}}-{{H}_{11}}&-{{H}_{12}}\\\
0&{{H}_{21}}&-{{H}_{21}}&0\\\ \end{matrix}\right).$ (28)

The expression of the superoperator $\mathcal{L}_{S}$ is then

$\mathcal{L}_{S}=\mathcal{L}_{S(ADO)}\otimes\prod_{k^{\prime\prime}=1}^{K}I_{n_{HEOM}}$
(29)

becoming $\mathcal{L}_{S}=\mathcal{L}_{S(ADO)}\otimes I_{2}\otimes I_{2}$ in
the example with $K=2$ and $n_{HEOM}=2$. The corresponding matrix is given in
figure 4 where we have adopted the concise notation
$L_{ij}=(\mathcal{L}_{S(ADO)})_{ij}$.

Figure 4: Matrix $\mathcal{L}_{S}$ (Eq.(29)). The $L_{ij}$ elements are the
concise notation for $(\mathcal{L}_{S(ADO)})_{ij}$ defined in Eq.(28). The
empty regions correspond to zero.

####  3.2.2 Damping Liouvillian $\mathcal{L}_{damp}$

The damping term is
$\mathcal{L}_{damp}=\sum_{k^{\prime}=1}^{K}\left(\mathcal{L}_{k^{\prime}}\right)$
where

$\mathcal{L}_{k^{\prime}}=i\gamma_{k^{\prime}}I_{n^{2}}\otimes\prod_{k^{\prime\prime}=1}^{K}M_{k^{\prime\prime}}$
(30)

with $M_{k^{\prime\prime}}=I_{n_{HEOM}}$ if $k^{\prime\prime}\neq k^{\prime}$
and $M_{k^{\prime\prime},lm}=(l-1)$ $\delta_{l,m}$ if
$k^{\prime\prime}=k^{\prime}$ ($l,m\in\left[1,n_{HEOM}\right]$). It is
diagonal in the superoperator representation. For each decay mode
$k^{\prime}$, all the elements of the matrices ${{\rho}_{\mathbf{n}}}$ are
multiplied by the decay rate ${{\gamma}_{k^{\prime}}}$ times the occupation
number ${{n}_{k^{\prime}}}$. In the two-level and two-decay case with
$n_{HEOM}=n_{max}$, one has:

$\displaystyle{\mathcal{L}_{1}}=i{{\gamma}_{1}}{{I}_{{{n}^{2}}}}\otimes
M\otimes{{I}_{{{n}_{NHEOM}}}}$
$\displaystyle=i{{\gamma}_{1}}\left(\begin{matrix}1&{}&{}&{}\\\ {}&1&{}&{}\\\
{}&{}&1&{}\\\ {}&{}&{}&1\\\
\end{matrix}\right)\otimes\left(\begin{matrix}0&{}&{}&{}\\\ {}&1&{}&{}\\\
{}&{}&\ddots&{}\\\ {}&{}&{}&{{n}_{max}}\\\ \end{matrix}\right)$
$\displaystyle\otimes\left(\begin{matrix}1&{}&{}&{}\\\ {}&1&{}&{}\\\
{}&{}&\ddots&{}\\\ {}&{}&{}&1\\\ \end{matrix}\right)$ (31)

and

$\displaystyle{\mathcal{L}_{2}}=i{{\gamma}_{2}}{{I}_{{{n}^{2}}}}\otimes{{I}_{{{n}_{NHEOM}}}}\otimes
M$ $\displaystyle=i{{\gamma}_{2}}\left(\begin{matrix}1&{}&{}&{}\\\
{}&1&{}&{}\\\ {}&{}&1&{}\\\ {}&{}&{}&1\\\
\end{matrix}\right)\otimes\left(\begin{matrix}1&{}&{}&{}\\\ {}&1&{}&{}\\\
{}&{}&\ddots&{}\\\ {}&{}&{}&1\\\ \end{matrix}\right)$
$\displaystyle\otimes\left(\begin{matrix}0&{}&{}&{}\\\ {}&1&{}&{}\\\
{}&{}&\ddots&{}\\\ {}&{}&{}&{{n}_{max}}\\\ \end{matrix}\right)$ (32)

When $n_{max}=1$ and thus $n_{NHEOM}=2$, ${{M}}=\left(\begin{matrix}0&{}\\\
{}&1\\\ \end{matrix}\right)$ and ${{I}_{{{n}_{NHEOM}}}}=I_{2}$. The
corresponding $\mathcal{L}_{damp}$ matrix divided by $i$ is given in figure 5.

|
---|---
Figure 5: Color on line. Damping superoperator
$(\mathcal{L}_{1}+\mathcal{L}_{2})/i$ (Eqs.(31) and (32)). (a) Example of a
two-state system with two artificial decay modes and $n_{max}=1$.
$\mathcal{L}_{1}$ in blue and $\mathcal{L}_{2}$ in red. (b) Block related to a
particular element $\rho^{\alpha}_{\mathbf{n}}$ where $\alpha=ij$, when
$n_{max}>1$. All the blocks $(ij)$ have the same structure

####  3.2.3 Matrix $\mathcal{L}_{+}$

Each term $\mathcal{L}_{k^{\prime}+}$ addressing the upper level of the
hierarchy in Eqs.(23) or (24) corresponds to the superoperator

$\mathcal{L}_{k^{\prime}+}=-i(S_{k^{\prime}}\otimes I_{n}-I_{n}\otimes
S_{k^{\prime}})\otimes\prod_{k^{\prime\prime}=1}^{K}M^{\prime}_{k^{\prime\prime}}.$
(33)

where $M^{\prime}_{k^{\prime\prime}}=I_{n_{HEOM}}$ if $k^{\prime\prime}\neq
k^{\prime}$ and $M^{\prime}_{k^{\prime\prime},lm}=\delta_{l+1,m}$ if
$k^{\prime\prime}=k^{\prime}$ ($l,m\in\left[1,n_{HEOM}\right]$). We consider
the case with a single bath with a system-bath coupling operator $S$
(Eq.(27)). Then $S_{k^{\prime}}$ is independant of $k^{\prime}$. The factor
$Q_{+}=(S\otimes I_{n}-I_{n}\otimes S)$ in the two-state case with $n=2$ is:

$\displaystyle Q_{+}=\left(\begin{matrix}0&0&0&0\\\ 0&-1&0&0\\\ 0&0&1&0\\\
0&0&0&0\\\ \end{matrix}\right).$ (34)

The two contributions $\mathcal{L}_{k^{\prime}+}$ for two decay modes and
$n_{max}>1$ are:

$\displaystyle{\mathcal{L}_{1+}}=-iQ_{+}\otimes
M^{\prime}\otimes{{I}_{{{n}_{NHEOM}}}}$
$\displaystyle=-i\left(\begin{matrix}0&0&0&0\\\ 0&-1&0&0\\\ 0&0&1&0\\\
0&0&0&0\\\ \end{matrix}\right)\otimes\left(\begin{matrix}0&1&{}&{}\\\
{}&0&\ddots&{}\\\ {}&{}&\ddots&1\\\ {}&{}&{}&0\\\ \end{matrix}\right)$
$\displaystyle\otimes\left(\begin{matrix}1&{}&{}&{}\\\ {}&1&{}&{}\\\
{}&{}&\ddots&{}\\\ {}&{}&{}&1\\\ \end{matrix}\right)$ (35)

and

$\displaystyle{\mathcal{L}_{2+}}=-iQ_{+}\otimes{{I}_{{{n}_{NHEOM}}}}\otimes
M^{\prime}$ $\displaystyle=-i\left(\begin{matrix}0&0&0&0\\\ 0&-1&0&0\\\
0&0&1&0\\\ 0&0&0&0\\\
\end{matrix}\right)\otimes\left(\begin{matrix}1&{}&{}&{}\\\ {}&1&{}&{}\\\
{}&{}&\ddots&{}\\\ {}&{}&{}&1\\\ \end{matrix}\right)$
$\displaystyle\otimes\left(\begin{matrix}0&1&{}&{}\\\ {}&0&\ddots&{}\\\
{}&{}&\ddots&1\\\ {}&{}&{}&0\\\ \end{matrix}\right).$ (36)

When $n_{max}=1$, ${{I}_{{{n}_{NHEOM}}}}=I_{2}$ and
$M^{\prime}=\left(\begin{matrix}0&1\\\ 0&0\\\ \end{matrix}\right)$ The
corresponding matrix of the sum $(\mathcal{L}_{1+}+\mathcal{L}_{2+})/(-i)$ is
displayed in figure 6. By comparing with Eq.(22), the contribution to
${{\dot{\rho}}_{00}}$ of the term related to the upper ADOs with only single
excitation is $-i\left[S,{{\rho}_{10}}+{{\rho}_{01}}\right]$, i.e.,
$\dot{\rho}_{00}^{12}\to-\rho_{01}^{12}-\rho_{10}^{12}$,
$\dot{\rho}_{00}^{21}\to\rho_{01}^{21}+\rho_{10}^{21}$ and
$\dot{\rho}_{00}^{11}=\dot{\rho}_{00}^{22}\to 0$. Matrix-vector products of
lines 5 and 9 with the column vector $\rho^{\alpha}_{\textbf{n}}$ in figure 6
provide these expressions for the contribution to $\dot{\rho}_{00}^{12}$ and
$\dot{\rho}_{00}^{21}$ respectively. The results for $\dot{\rho}_{00}^{11}$
and $\dot{\rho}_{00}^{22}$ can be obtained in the same way at lines 1 and 13.

Figure 6: Color online. Matrix $(\mathcal{L}_{1+}+\mathcal{L}_{2+})/(-i)$ in
the case of a two-state system with two artificial decay modes and
$n_{max}=1$. $\mathcal{L}_{1+}$ in blue and $\mathcal{L}_{2+}$ in red.

####  3.2.4 Matrix $\mathcal{L}_{-}$

The superoperator connecting the ADOs with a lower layer in the hierarchy
involves the operator:

$\displaystyle\mathcal{L}_{k^{\prime}-}=$
$\displaystyle-i(\alpha_{k^{\prime},t/c}S_{k^{\prime}}\otimes
I_{n}-\tilde{\alpha}_{k^{\prime},t/c}I_{n}\otimes S_{k^{\prime}})$
$\displaystyle\otimes\prod_{k^{\prime\prime}=1}^{K}M^{\prime\prime}_{k^{\prime\prime}}$
(37)

where $M^{\prime\prime}_{k^{\prime\prime}}=I_{n_{HEOM}}$ if
$k^{\prime\prime}\neq k^{\prime}$ and
$M^{\prime\prime}_{k^{\prime\prime},lm}=(l-1)$ $\delta_{l-1,m}$ if
$k^{\prime\prime}=k^{\prime}$ ($l,m\in\left[1,n_{HEOM}\right]$). In the two-
level case with $S_{k^{\prime}}$ independent of $k^{\prime}$, for each mode,
one has the factor:

$\displaystyle{{Q}_{-k^{\prime}}}={{\alpha}_{k^{\prime}}}{{S}}\otimes{{I}_{n}}-{{{\tilde{\alpha}}}_{k^{\prime}}}{{I}_{n}}\otimes{{S}}$
$\displaystyle={{\alpha}_{k^{\prime}}}\left(\begin{matrix}0&0\\\ 0&1\\\
\end{matrix}\right)\otimes\left(\begin{matrix}1&0\\\ 0&1\\\
\end{matrix}\right)-{{{\tilde{\alpha}}}_{k^{\prime}}}\left(\begin{matrix}1&0\\\
0&1\\\ \end{matrix}\right)\otimes\left(\begin{matrix}0&0\\\ 0&1\\\
\end{matrix}\right)$ $\displaystyle=\left(\begin{matrix}0&0&0&0\\\
0&{{{-\tilde{\alpha}}}_{k^{\prime}}}&0&0\\\ 0&0&{{\alpha}_{k^{\prime}}}&0\\\
0&0&0&{{\alpha}_{k^{\prime}}}-{{{\tilde{\alpha}}}_{k^{\prime}}}\\\
\end{matrix}\right).$ (38)

The two contributions to the super-Liouvillian when $n_{max}>1$ are:

$\displaystyle{\mathcal{L}_{1-}}={{Q}_{-1}}\otimes
M^{\prime\prime}\otimes{{I}_{{{n}_{HEOM}}}}$
$\displaystyle=\left(\begin{matrix}0&0&0&0\\\
0&{{{-\tilde{\alpha}}}_{1}}&0&0\\\ 0&0&{{\alpha}_{1}}&0\\\
0&0&0&{{\alpha}_{1}}-{{{\tilde{\alpha}}}_{1}}\\\
\end{matrix}\right)\otimes\left(\begin{matrix}0&{}&{}&{}&{}\\\ 1&0&{}&{}&{}\\\
{}&{2}&\ddots&{}&{}\\\ {}&{}&\ddots&\ddots&{}\\\ {}&{}&{}&{{n}_{\max}}&0\\\
\end{matrix}\right)$ $\displaystyle\otimes\left(\begin{matrix}1&{}&{}&{}\\\
{}&1&{}&{}\\\ {}&{}&\ddots&{}\\\ {}&{}&{}&1\\\ \end{matrix}\right)$ (39)

and

$\displaystyle{\mathcal{L}_{2-}}={{Q}_{-2}}\otimes{{I}_{{{n}_{HEOM}}}}\otimes
M^{\prime\prime}$ $\displaystyle=\left(\begin{matrix}0&0&0&0\\\
0&{{{-\tilde{\alpha}}}_{2}}&0&0\\\ 0&0&{{\alpha}_{2}}&0\\\
0&0&0&{{\alpha}_{2}}-{{{\tilde{\alpha}}}_{2}}\\\
\end{matrix}\right)\otimes\left(\begin{matrix}1&{}&{}&{}\\\ {}&1&{}&{}\\\
{}&{}&\ddots&{}\\\ {}&{}&{}&1\\\ \end{matrix}\right)$
$\displaystyle\otimes\left(\begin{matrix}0&{}&{}&{}&{}\\\ 1&0&{}&{}&{}\\\
{}&{2}&\ddots&{}&{}\\\ {}&{}&\ddots&\ddots&{}\\\ {}&{}&{}&{{n}_{\max}}&0\\\
\end{matrix}\right).$ (40)

In the example of the two decay modes with $n_{max}=1$, one has
$M^{\prime\prime}=\left(\begin{matrix}0&0\\\ 1&0\\\ \end{matrix}\right)$. The
matrices for $\mathcal{L}_{1-}$ and $\mathcal{L}_{2-}$ are given in figure 7.

Figure 7: Color on line. Matrix $(\mathcal{L}_{1-}+\mathcal{L}_{2-})/(-i)$ in
the case of a two-state system with two artificial decay modes and
$n_{max}=1$. $\mathcal{L}_{1+}$ in blue and $\mathcal{L}_{2+}$ in red.

###  3.3 Dynamics with Tensor-Train format

Dynamics is driven by solving

$\mathbf{\dot{\bar{\rho}}}(t)=\mathcal{L}\mathbf{\bar{\rho}}(t)$ (41)

where $\bar{\rho}$ is the full TT-converted vector of the elements of all the
ADOs and $\mathcal{L}$ is the super-operator described in Sec.3.2. We use the
projector-splitting KSL scheme (Lubich2014, ; Lubich2015, ; Oseledets2016, ;
Batista2022, ) implemented in the ttpy package (`tt.ksl.ksl`) ttpy . The
method is based on the dynamical low-rank approximation which is equivalent to
the Dirac-Frenkel time-dependent variational principle used in MCTDH. It
consists in using an approximate low-rank tensor with fixed ranks instead of
getting a solution with a high rank tensor and then truncate it with singular
value decomposition (SVD). To comply with this goal, the derivative of the
approximate low-rank tensor is obtained by orthogonally projecting the
derivative of the tensor on the tangent space of the approximate low-rank
tensor at its current position. Time-integration is then obtained by a
splitting scheme (second order in this work) of the projector (see Refs.
(Lubich2014, ; Lubich2015, ; Oseledets2016, ; Batista2022, ) for more
details). An adaptative rank may be necessary during the propagation as
proposed in Refs.DunnettChin2021 ; Dolgov19 . We have adopted a mixed
strategy. The standard Runge-Kutta integrator (written with TT algebra
available with the ttpy package) is run after some time steps to allow the
increase of the ranks during the propagation. In the application, we use a
Runge-Kutta run after 10 timesteps.

##  4 Illustrative Applications

We give three examples for which dynamics is driven by the TT method. In the
first two cases, the spectral density is continuous and is an Ohmic two-pole
Lorentzian function (Eq.(9)) leading to two artificial bath modes. The third
application uses undamped discrete modes.

###  4.1 Population to coherence transfer via a bath

The zero-order model and the eigenstates are schematized in figure 8. This
model was introduced to analyze one of the first applications of an
experimental ‘quantum simulator’ for molecular quantum dynamics Potocnik2018 .
As we show, using a circuit of qubits to represent a network of chromophores
gives access to strongly non-Markovian regimes of open dynamics, including
situations in which strong dissipation actually induces coherent dynamics.
Indeed, one of the original motivations for the experiment and analysis of
Ref.(ChinPRA2018, ) was to explore the existence, robustness and uses of
coherent transport in photosynthetic light-harvesting proteins, which are
described by analogous models. However, they are much harder to control
compared to superconducting circuits. The Hamiltonian of the active subsystem
is:

$\displaystyle H_{S}(t)=H_{s}+H_{ren}$
$\displaystyle-\sum_{n=0}^{3}\sum_{m\neq{n}=0}^{3}\mu_{nm}|n\rangle\langle
m|\mathcal{E}(t)$ (42)

with

$H_{s}=\sum_{n=0}^{3}\varepsilon_{n}|n\rangle\langle
n|+\sum_{n=0}^{3}\sum_{m\neq{n}=0}^{3}H_{nm}|n\rangle\langle m|$ (43)

where $H_{ren}$ is the renormalization term given below and $\mu_{ij}$ is the
dipolar coupling. The ground state is coupled to the excited states only
radiatively, i.e. $H_{0j}=0$ for $j\neq 0$ and only $\mu_{02}$ and $\mu_{03}$
induce the radiative coupling. Two degenerate excited bright states
($\left|2\right\rangle$ and $\left|3\right\rangle$) strongly interact by
interstate coupling $H_{23}$ and state $\left|2\right\rangle$ is weakly
coupled to a dark lower state $\left|1\right\rangle$ and to a tuning bath that
makes fluctuate the energy. The strong $H_{23}$ coupling leads to eigenstates
that are mainly the bright in phase and the dark out of phase superpositions
denoted $\left|B\right\rangle$ and $\left|D_{+}\right\rangle$ respectively.
Both eigenstates $\left|D_{+}\right\rangle$ and $\left|D_{-}\right\rangle$
form a dark doublet. Their dipole transition moments $\mu_{0D_{\pm}}$ are very
weak, being two orders of magnitude smaller than the transition moment
$\mu_{0B}$ to the $\left|B\right\rangle$ state. Spontaneous radiative decay is
not taken into account in the simulation.

Figure 8: Schematic representation of the energy levels, the interstate
couplings $H_{ij}$, the dipolar couplings $\mu_{ij}$ and the system-bath
couplings of a device in which a long-lasting coherence in a dark doublet may
be created by interaction with a bath after excitation of a bright state. (a)
zero-order basis, (b) eigenstates.

The bath is coupled to state $\left|2\right\rangle$ only. The corresponding
system-bath operator takes the form:

$\displaystyle S=\left(\begin{matrix}0&0&0&0\\\ 0&0&0&0\\\ 0&0&1&0\\\
0&0&0&0\\\ \end{matrix}\right).$ (44)

The spectral density is highly structured and peaks nearly at the mean
($\Delta_{B-D}=\langle B|H_{s}|B\rangle-\langle D|H_{s}\left|D\right\rangle$)
transition energy, $\langle D|H_{s}\left|D\right\rangle$ being the average
energy of states $\langle D_{+}|H_{s}\left|D_{+}\right\rangle$ and $\langle
D_{-}|H_{s}\left|D_{-}\right\rangle$. It is a single Ohmic Lorentzian (Eq.(9))
with parameters $p=2.0\times 10^{-12}a.u.$, $\Omega=4.5\times 10^{-3}a.u.$ and
$\Gamma=4.0\times 10^{-4}a.u.$. It is represented with the correlation
function at $T$ = 298K in figure 9. The renormalisation energy is
$\lambda=(1/\pi)\int_{0}^{\infty}{d\omega J(\omega)/\omega}$. The
corresponding renormalization term is
$H_{ren}=\lambda\left|2\right\rangle\left\langle 2\right|$. The system-bath
coupling is weak. It is estimated by the ratio
$\eta=\lambda/{\Delta_{{(\left|B\right\rangle-\left|D\right\rangle)}}}$ where
$\Delta_{{(\left|B\right\rangle-\left|D\right\rangle)}}$ is the energy gap
corresponding here to the cutoff of the spectral density. It is a perturbative
regime, however, the correlation time is long (about 250 fs) due to the peaked
shape of the spectral density. It is longer than the Rabi period of the
($\left|B\right\rangle-\left|D\right\rangle$) transition that is 34 fs.
Dynamics is then non-Markovian with $n_{HEOM}=5$.

Figure 9: (a) Spectral density of the model presented in figure 8 to
illustrate the population-to-coherence transfer. The spectral density maximum
is close to the mean $BD$ energy gap; (b) corresponding normalized correlation
function $C(t)/C(0)$ (Eq.(4)) at $T$ = 298K with real part (dashed line),
imaginary part (dots) and modulus (solid line).

In the eigenstate representation (Fig.8(b)), the $S$ operator becomes:

$\displaystyle{{S}_{eigen}}=\left(\begin{matrix}0&0&0&0\\\ 0&0.26&0.25&0.36\\\
0&0.25&0.24&0.34\\\ 0&0.36&0.34&0.50\\\ \end{matrix}\right).$ (45)

In this basis set, the bath is coupled diagonally and off-diagonally to the
system. This is an interesting device leading to a population-to-coherence
transfer from the bright $\left|B\right\rangle$ state to the dark doublet.
This process requires exact non-Markovian dynamics taking terms that are
neglected in the secular approximation of a Redfield treatment (Lovett2016, ;
ChinPRA2018, ). From the analysis in Ref.(ChinPRA2018, ), we choose a coupling
strength $\eta=0.013$ relevant to illustrate an efficient population-to-
coherence process. The populations in the eigenstates and the modulus of the
coherence $(\rho_{S})_{D_{+}D_{-}}$ between the two dark states are displayed
in figures 10(a) and 10(c) when the initial state is the bright state
$\left|B\right\rangle$. In this ideal case, the population-to-coherence is
complete in 250 fs. The coherence is long-lasting and slowly decays in about
25 ps. Figures 10(b) and 10(d) compare this ideal preparation with the results
obtained by an excitation from the ground state by a laser field of 24 fs. To
respect the condition that the area of the oscillating field must be equal to
zero Brabec2000 ; Milosevic2006 , the field is then given by
$\mathcal{E}(t)=-\frac{\partial\mathcal{A}(t)}{\partial t}$ with the vector
potential
$\mathcal{A}(t)=\left(\frac{{{\mathcal{E}}_{0}}}{\omega}\right){{\sin}^{2}}\left(\frac{\pi(t-{{t}_{i}})}{\tau}\right)\sin\left(\omega(t-{{t}_{i}})\right)$,
where $t_{i}$ is the initial time of the pulse, $t_{i}$ = 0 here, $\tau$ is
the pulse duration and $\mathcal{E}_{0}$ is the field maximum amplitude. The
carrier frequency is in resonance with the $|0\rangle\rightarrow|B\rangle$
transition. When the number of cycles is large for a pulse duration longer
than about 20 fs, the expression becomes similar to a pulse with a sine square
envelope ${\mathcal{E}}(t)={{\mathcal{E}}_{0}}{{\sin}^{2}}\left(\frac{\pi
t}{\tau}\right)\cos\left({{\omega}_{0B}}t\right)$. We use this expression to
estimate $\mathcal{E}_{0}$ providing a $\pi$ pulse for which the integral of
the Rabi frequency $\Omega(t)=\mathcal{E}_{0}{{\sin}^{2}}\left(\frac{\pi
t}{\tau}\right)\mu_{0B}/\hbar$ is equal to $\pi$ (Thomas1983, ; Holthaus1994,
). $\mu_{0B}$ is the transition dipole. Then
${\mathcal{E}_{0}}=2\pi/({{\mu}_{0B}}{{\tau}})$. This field would induce a
complete transfer towards the bright state in the absence of coupling to the
bath. A slight modification of the amplitude should be necessary for short
pulses with few cycles for which the envelope slighly differs from a sine
square.

Figure 10: Upper panels: populations in the eigenstates $B$, $D_{-}$ and
$D_{+}$, (a) ideal preparation in the bright state $B$, (b) excitation of the
$B$ state from the ground state by a $\pi$ pulse of 24 fs. lower panels:
modulus of the coherence $(\rho_{S})_{D_{+}D_{-}}$, (c) ideal preparation in
the bright state, (d) excitation by a $\pi$ pulse.

Figure 11 shows the influence of the pulse duration $\tau$ on the coherence
generation in the dark doublet. The result is close to the ideal case for
pulses with $\tau$ in the range 20-50 fs. The decrease of the yield comes from
the environment that makes fluctuate the energies.

Simulations are made with a timestep of 0.24 fs and the maximum rank is 10
with $\epsilon=10^{-15}$.  The temperature is 298 K. No Matsubara term is
needed due to the high temperature and the narrow spectral density. No
significant change has been observed when carrying out this simulation with 3
additional Matsubara terms.

Figure 11: Excitation of the bright state $B$ by a $\pi$ pulse of different
durations. (a) Population in the ground state, (b) Modulus of the coherence
$(\rho_{S})_{D_{+}D_{-}}$ generated in the dark doublet .

###  4.2 Simulation of absorption and emission spectra

HEOM have already been used to simulate stationary absorption $\sigma_{abs}$
and emission $\sigma_{em}$ spectra Tanimura2006 ; Tanimura_rev_20 ; Shi2013 ;
Cao2014 ; Yan2017 ; Kramer2018 . They are computed by the linear response
theory as:

${{\sigma}_{abs}}(\omega)={\textrm{Re}\int_{0}^{\infty}{dt{{e}^{i\omega
t}}T{{r}_{S}}}\left[\rho_{{{\mu}^{-}}}^{\dagger}(t){{\rho}_{{{\mu}^{-}}}}(0)\right]}$
(46) ${{\sigma}_{em}}(\omega)={\textrm{Re}}\int_{0}^{\infty}{dt{{e}^{i\omega
t}}T{{r}_{S}}}\left[\rho_{{{\mu}^{+}}}^{\dagger}(t){{\rho}_{{{\mu}^{+}}}}({{t}_{eq}})\right]$
(47)

where ${{\mu}^{-}}=\sum\nolimits_{k\neq
0}{{{\mu}_{k}}}\left|0\right\rangle\left\langle k\right|$ and
${{\mu}^{+}}=\sum\nolimits_{k\neq
0}{{{\mu}_{k}}}\left|k\right\rangle\left\langle 0\right|$. The initial
conditions for the absorption is
${{\rho}_{{{\mu}^{-}}}}(0)={{\mu}^{-}}{{\rho}_{S}}(0)$ where ${{\rho}_{S}}(0)$
is the reduced density matrix of the system in its ground state and all the
ADOs are zero. The stationary emission is computed at the thermal equilibrium
of the emitting state. The first strategy is to propagate an arbitrary system
density matrix with population in the excited states towards the thermally
equilibrated state, which is independent of the chosen initial state
(Dijkstra2010, ; Shi2013, ; Kramer2018, ). Then, the initial conditions for
the emission is
${{\rho}_{{{\mu}^{+}}}}(t_{eq})={{\mu}^{+}}{{\rho}_{S}}(t_{eq})$ where
${{\rho}_{S}}(t_{eq})$ is the reduced system density matrix at equilibrium and
the ADOs take their asymptotic values at time $t_{eq}$. $t_{eq}$ is estimated
by verify the population and the average value of the collective bath mode.
Another possibility to reach the equilibrium involves propagation with an
imaginary time Pomyalov2010 ; Tanimura2014 ; Shi2015 .

Figure 12:  Schematic representation of the diabatic potential energy surface
of the excited states used in the simulation of absorption and emission
spectra. $q$ is a bath nuclear coordinate and $\Delta_{q_{e}}$ is the
difference of equilibrium positions. The bath is at equilibrium in the ground
state and in the dark state $\left|D\right\rangle$ since the equilibrium
positions are assumed to be the same. The electronic energies of the system
correspond to Franck-Condon vertical energies represented by full circles.
$\lambda$ is the renormalization energy. Two cases are simulated, the first
one with $H_{12}=0$ (the $\left|D\right\rangle$ state is spectator) and the
second one with ${{H}_{12}}\neq 0$

.

Figure 12 schematizes the model. According to the value of the interstate
coupling $H_{12}$, we consider a single bright state when ${{H}_{12}}=0$ or a
two-excited-state case with a nonadiabatic coupling with the dark state when
${{H}_{12}}\neq 0$. The spectra are computed by making the electronic-nuclear
partition. The system is composed of two or three electronic states
corresponding to vertical energies at the ground state equilibrium geometry.
The bath consists in all the nuclear intermolecular and solvent vibrators.
Within this partition, the system-bath couplings
$c_{j}=\omega^{2}_{j}\Delta_{eq,j}$ depend on the difference
${{\Delta}_{{{q}_{e}}}}$ between the equilibrium position of all the
vibrational modes in the ground and excited state. In our example, only state
$\left|B\right\rangle$ is displaced and therefore coupled to the bath. In the
first example ($H_{12}=0$), the system-bath coupling operator is that given in
Eq.(27). In the second case, in the diabatic representation, it becomes

$\displaystyle S=\left(\begin{matrix}0&0&0\\\ 0&1&0\\\ 0&0&0\\\
\end{matrix}\right).$ (48)

We intend to illustrate the strong influence of the spectral density on the
Stokes shift by referring to Mukamel’s $\kappa$ parameter Mukamel1995

$\kappa=\frac{\Lambda}{\Delta}.$ (49)

$\Lambda$ corresponds to the cutoff of the spectral density.
${{\Lambda}^{-1}}$ is an estimation of the bath fluctuation timescale.
${{\Delta}^{2}}=C(t=0)$ is related to the initial value of the bath
correlation function, which is a real value obtained from Eq.(7). $\Delta$ is
an estimation of the amplitude of the fluctuations.

Figure 13 displays the four spectral densities adopted in the simulation and
the corresponding $\kappa$.

(i) The case with a single excited bright state is the basic example for which
the behavior is well predicted by the $\kappa$ value.

Figure 13: Spectral densities used in the simulation of the absorption and
emission spectra. $\kappa$ is defined in Eq.(49). The $\Delta$ parameter is
given by the initial value of the correlation function (Eq.(4)). The cutoff
$\Lambda$ is estimated by the energy at the maximum. The renormalization
energy $\lambda$ for decreasing values of $\kappa$ are 0.034 eV, 0.17 eV, 0.34
eV and 0.52 eV.

When $\kappa>1$, the bath dynamics become fast and this situation is known as
a homogeneous dephasing case leading to Lorentzian profiles for linear
absorption or relaxed emission spectra and no Stokes shift. On the contrary,
when $\kappa<1$, bath dynamics become slow, this is the static limit or
inhomogeneous case for which the profiles become Gaussian and the maximum
Stokes shift is given by $2\lambda$. The four simulations are presented in
figure 14. One could observe the expected behavior passing from no Stokes
shift when $\kappa$ = 1.45 to a large shift when $\kappa$ = 0.21 for which
$\lambda$ = 0.52 eV. In this example, the relaxation is simply the evolution
towards the thermal mixture at the new equilibrium geometry of the excited
bright state. This stationary state is obtained in about 150 fs. The ADOs are
taken after a propagation of 180 fs. Generally, the computation of the
emission spectrum is more demanding than for the absorption. Convergence of
the correlation function requires higher HEOM level (up to level 55 for
$\kappa$ = 0.21). Propagation with the TT method has required a small time
step of 0.012 fs. The maximum tensor rank remains below 10 with
$\epsilon=10^{-20}$. The spectra are computed at 298 K. We have verified for
the case with $\kappa=0.21$ that no Matsubara term is necessary.

Figure 14: Normalized absorption (solid lines) and emission (dashed lines)
spectra for the case with a single bright state ($H_{12}=0$) represented in
figure 12. The panels correspond to the different spectral densities given in
figure 13. $\kappa$ is defined in Eq.(49).

(ii) In the second example presented in figure 15, dynamics is more
complicated since the emission occurs from the eigen vibronic states having a
component on the bright state. The propagation duration to reach the
asymptotic populations and prepare the ADOs depends on $\lambda$. It is 500 or
600 fs for the different examples. The absorption spectrum is more affected by
the nonadiabatic interaction. For the emission spectrum, one recovers the
shape corresponding to the relaxed bright state in this example since the
second excited state is assumed to be dark. The evolution of the Stokes shift
with respect to the $\kappa$ parameter is the same and corresponds to the
expected behavior.

These results highlight that tools like HEOM can now predict optical spectra
with great precision, given the spectral density. Therefore, attention must
now be given to extracting high quality spectral densities to match
experimental results, particulary in condensed phases where there can be very
distinct timescales in the environment, i.e. fast intramolecular motions and
slow solvent/lattice/protein reorganisation. The importance of including the
latter and a way of obtaining them from first principles was given in
Ref.Dunnett2021bis .

Figure 15: Normalized absorption (solid lines) and emission (dashed lines)
spectra for the case with a bright state coupled to a dark state
(${{H}_{12}}\neq 0$) represented in figure 12. The panels correspond to the
different spectral densities given in figure 13. $\kappa$ is defined in
Eq.(49).

###  4.3 Excitation transfer with discrete bath modes

In this example, the bath is formed by discrete undamped vibrational modes. We
revisit the excitation transfer in a dimer oligothiophene (OT4)-fullerene(C60)
investigated by I. Burghardt and coworkers Tamura2011 ; Burghardt2012 ;
Chenel2014 ; Burghardt2021 . The excited states may be schematized as in
Fig.12 with $D$ being the excitonic state XT and $B$ becoming the charge
transfer (CT) state. We consider the parameters for an inter-fragment distance
of 0.35 nm Chenel2014 . The energy gap is 0.07 eV and the electronic coupling
is 0.007 eV. The continuous spectral density has been obtained from the
difference of the equilibrium position $\Delta q_{eq}$ of the fragment normal
modes Tamura2011 . The reference position of the equilibrated bath is taken at
the middle point $\Delta q_{eq}/2$ so that the coupling operator is

$\displaystyle S=\left(\begin{matrix}-0.5&0\\\ 0&0.5\\\ \end{matrix}\right).$
(50)

The smooth spectral density Burghardt2012  is fitted by five two-pole
Lorentzians (Eq.(9)) and displayed in Fig.16(a). The parameters are given in
the supplementary material of Ref.Chenel2014 . This leads to $K=10$ artificial
decay modes without Matsubara terms at 298 K. We have used the discretization
procedure of Ref.Makri2017  providing unequally spaced modes that correspond
to equal fractions of the reorganization energy. The discrete spectral density
at the $N_{disc}=40$ selected points $k$ ($\pi
c_{k}^{2}/(2\omega_{k}\Delta(\omega))$ where $\Delta(\omega)$ is the local
state density) is shown in Fig.16(b). Since two baths are involved in the
discrete procedure, this corresponds to $K=80$ modes. Figure 16(c) gives the
occupation probability of the XT and CT states when the initial state is XT.
The dashed lines are obtained with the continuous spectral density and the
five artificial modes, at level 5 of the hierarchy, with the tolerance
parameter $\epsilon=10^{-10}$ and the maximum rank (rmax) is 20 (see
Appendix). This example requires only 3003 matrices and the TT formulation is
not more advantageous than the standard one. However, the discrete case may
emphasize the utility of the TT propagation when compared with a standard
treatment. The corresponding population are drawn in solid lines in Fig.16(c).
We may estimate the total number of elements in the full matrix array. One has
$n=2$, $K=80$, $L=5$ leading to
${{N}_{standard}}={{n}^{2}}(K+L)!/\left(K!L!\right)=$ 131,206,068 elements.
This is a computationally heavy simulation ($\approx$ 21 Go for the density
matrix only) and intractable with our current Fortran implementation. In the
TT implementation, one has $n=2$, $K=80$, $L=5$ and $r_{max}=80$. The number
of stored elements might reach a maximum value of
${{N}_{TT}}=({{n}^{2}}+L)*{{r}_{max}}+{{r}_{max}}^{2}*(K-1)*L$ = 2,528,720.
The storage is obviously more attractive ($\approx$ 400 Mo). The tolerance is
$\epsilon=10^{-8}$. The norm is well conserved up to about 250 fs. The result
is promising but it will be necessary to improve the adaptation of the maximum
rank to further extend the performance Borrelli2021 ; DunnettChin2021 .

Figure 16: Excitation energy transfer dynamics with a continuous or
discretized spectral density by the TT implemetation. (a) : continuous
spectral density fitted by five two-pole Lorentzian functions (Eq.(9)), see
Refs.Burghardt2012 ; Chenel2014 , (b) : discretization in $N_{disc}=40$ modes
by the Makri procedure (see Ref.Makri2017 ) leading to $K=2N_{disc}$ decay
modes, (c) : population in the XT (black) and CT (red) states when the initial
state is XT. Dashed lines: TT implementation with the continuous spectral
density leading to $K=10$ decay modes without Matsubara terms at 298 K; solid
lines : TT implementation with $K=80$ discrete modes .

##  5 Conclusion

We have given a detailed description of the TT formulation of HEOM and how it
is implemented which we hope will be a useful guide to newcomers to the field.
This is an exciting development, linking the widely used HEOM method with
rapidly developing advances in tensor networks across physics and theoretical
chemistry. One could connect to them the recent explosion of physical science
applications of machine learning, for example: tensor network simulation of
multi-environmental open quantum dynamics via machine learning
AspuruGuzik_2017 ; Zhao2018  and entanglement renormalisation Schroder2019 .

We have shown three examples where the TT method allows us to obtain highly
accurate results for complex phenomena, such as ‘noisy’ generation of
coherences (population to coherence transfer) and the dramatic impact of bath
relaxation times on optical spectra. The latter is important, as Stokes shifts
are often used as a direct measure of system-bath coupling. This is clearly
only one part of the story, as the relaxation time is also important. Hence
suitable techniques that can handle multiple timescales (long and short-lived
ADOs) are essential. These effects described above will be even more important
in larger multistate systems. The third example simulates an ultrafast
excitation transfer by using only undamped decay modes. This example
illustrates the efficiency of the TT formulation when compared to the standard
method, which would involve a huge number of matrices. All the applications
involve non-Markovian dynamics, either due to long correlation times or high
level of hierachy.

It would be interesting to consider extensions to time-dependent
spectroscopies, 2D-spectroscopies Tanimura_rev_20 ; Fleming2019 ; Zhu2011 ;
Tanimura2012  and optical control. Adaptation of the TT formulation for
optimal control has already been given in Ref.Jaouadi2022 . Another important
prospect is the consideration of non-harmonic environments, for instance by
including some modes in the active system Chenel2014 ; MangaudDMP2015 ;
Nazir2016 ; Gelin2016  and the treatment of fermionic baths Thoss2021 ;
Lambert2020 .

## Appendix: Numerical implementation of HEOM with ttpy

Most of the TT algebra is carried out by the ttpy package developed by
Oseledets and coworkers ttpy . In this appendix, we show a minimal code to
build the system Liouvillian and time-integrate a given system density matrix
with initial system-bath factorization in TT format with Python3, Numpy and
ttpy packages.

The system Liouvillian is defined as
${\mathcal{L}_{S(ADO)}}=-i(H\otimes{{I}_{n}}-{{I}_{n}}\otimes H)$ where $H$ is
the system Hamiltonian, $I_{n}$ the identity matrix with $n$ the number of
system states. Thus, ${\mathcal{L}_{S(ADO)}}$ is a $n^{2}\times n^{2}$ matrix
which can be built from standard numpy functions (`np.eye` returns the
identity matrix and `np.kron` the Kronecker product of both matrices) :

    ids = np.eye(n)
    Ls = -j * (np.kron(hamil,ids)  \
    - np.kron(ids,hamil))

where `ids` is the identity matrix of size $n\times n$ and ` j` the imaginary
unit. To convert this numpy array to a TT format, `tt.matrix` routine performs
an approximation of ${\mathcal{L}_{S(ADO)}}$ for a maximal rank (`rmax`) and
an accuracy `eps` :

    Ls = tt.matrix(Ls, eps=eps, rmax=rmax)

where `Ls` is the system Liouvillian superoperator. At this point, we are
still only dealing with a representation of an array of $n^{2}\times n^{2}$
dimensions. The HEOM superoperator $\mathcal{L}_{S}$ which spans over the
whole Liouville space is defined as
$\mathcal{L}_{S}=\mathcal{L}_{S(ADO)}\otimes\prod_{k^{\prime\prime}=1}^{K}I_{n_{HEOM}}$.
To avoid memory issues due to the high dimensionality of the array, one must
work with the Kronecker products ( `tt.kron`) of ttpy packages instead of the
one of numpy (`np.kron`). Indeed, numpy will build the full tensor which can
be very large and thus might suffer from the dimensionality curse. To carry
out this task, a single loop iterates over the number of artificial decay
modes $K$ with successive Kronecker products:

    idr = np.eye(nheom)
    idr = tt.matrix(idr)
    for i in range(K):
        Ls = tt.kron(Ls,idr)

The total Liouvillian super-operator $\mathcal{L}$ is a sum of several other
super-operators, i.e. $\mathcal{L}_{k^{\prime}}$, $\mathcal{L}_{k^{\prime}+}$,
$\mathcal{L}_{k^{\prime}-}$ (see Eq.(26)). Addition can be performed directly
with the usual algebraic symbol (+) on TT objects. The ttpy library carries
out automatically the correct tensor operations. However, tensor ranks
increase at each iteration. Thus, rounding operations are regularly performed
to reduce the rank for a given accuracy and maximum rank with the following
command:

    L = L.round(eps,rmax)

where `L` is the total HEOM Liouvillian super-operator.

The initial vectorized density matrix is defined directly from its cores. The
first core is filled with the initial system density matrix. As all auxiliary
density matrices are vectors of zeros when assuming system-bath factorization
all the other cores are vectors of dimensions $n_{HEOM}$ with only the first
index equal to 1.

For a given system density matrix `rho` defined as a numpy array and r the
initial core rank (equal to 1) in this example, we compute `rho` the initial
density matrix in TT format with the following algorithm:

    rho = np.zeros((n**2), \
          dtype=np.complex128)
    rho = np.reshape(rhos,n**2)
    cores =[]
    a = np.zeros((1,n**2,r))
    a[0,:,0] = rho.copy()
    cores.append(a)
    vec = np.zeros((nheom), \
          dtype=np.complex128)
    vec[0] = 1.
    for i in range(K-1):
        a = np.zeros((r,nheom,r))
        a[0,:,0] = vec[:]
        cores.append(a)
    a = np.zeros((r,nheom,1))
    a[0,:,0] = vec[:]
    cores.append(a)
    rho = tt.vector().from_list(cores)}

Time-integration is performed with the KSL second-order splitting algorithm
Lubich2014 . For each time step `dt`, the density matrix in TT format is
updated by using the function ksl implemented in tt.ksl routines :

    rho = tt.ksl.ksl(L,rho,dt)

By iterating over the desired number of timesteps, we compute the full density
matrix at time $t$. In order to extract the system density matrix, projection
techniques (by expressing projectors on the full Liouville space in the TT
format) or core manipulations can be used.

Acknowledgments

This work was performed within the French GDR 686 3575 THEMS. We want to
dedicate this work to the memory of Christoph Meier and Osman Atabek. Both
authors made significant contributions to the research areas described above,
such as the Meier-Tannor spectral density and the introduction of the
auxiliary matrices. Both have also greatly contributed to many strategies of
quantum control in a wide range of processes such as molecular orientation,
dynamics of excited states, isomerisation, molecular cooling to cite only few.
Dominik Domin is warmly acknowledged for his efficient technical support.

## Declarations

  * •

Funding: No funding was received for conducting this study.

  * •

Conflict of interest/Competing interests: The authors have no relevant
financial or non-financial interests to disclose.

  * •

Ethics approval: Not applicable.

  * •

Availability of data and materials: This manuscript has no associated data or
the data will not be deposited. [Authors’ comment: The datasets generated
during and/or analysed during the current study are available from the
corresponding author on reasonable request.].

  * •

Authors’ contributions: All the authors were involved in the preparation of
the manuscript. All the authors have read and approved the final manuscript.

## References

  *   * (1) U. Manthe, A multilayer multiconfigurational time-dependent hartree approach for quantum dynamics on general potential energy surfaces.  J. Chem. Phys. 128(16), 164,116 (2008).  10.1063/1.2902982
  * (2) H. Wang, M. Thoss, Numerically exact quantum dynamics for indistinguishable particles: The multilayer multiconfiguration time-dependent hartree theory in second quantization representation.  J. Chem. Phys. 131(2), 024,114 (2009).  10.1063/1.3173823
  * (3) O. Vendrell, H.D. Meyer, Multilayer multiconfiguration time-dependent hartree method: Implementation and applications to a henon–heiles hamiltonian and to pyrazine.  J. Chem. Phys. 134(4), 044,135 (2011).  10.1063/1.3535541
  * (4) L.P. Lindoy, B. Kloss, D.R. Reichman, Time evolution of ml-mctdh wavefunctions. ii. application of the projector splitting integrator.  The Journal of Chemical Physics 155(17), 174,109 (2021).  10.1063/5.0070043.  URL https://doi.org/10.1063/5.0070043
  * (5) T.G. Kolda, B.W. Bader, Tensor decompositions and applications.  SIAM Review 51(3), 455–500 (2009).  10.1137/07070111X
  * (6) L. Grasedyck, D. Kressner, C. Tobler, A literature survey of low-rank tensor approximation techniques.  GAMM-Mitteilungen 36(1), 53–78 (2013).  https://doi.org/10.1002/gamm.201310004
  * (7) I.V. Oseledets, Tensor-train decomposition.  SIAM Journal on Scientific Computing 33(5), 2295–2317 (2011).  10.1137/090752286
  * (8) C. Lubich, I.V. Oseledets, B. Vandereycken, Time integration of tensor trains.  SIAM Journal on Numerical Analysis 53(2), 917–941 (2015).  10.1137/140976546
  * (9) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, F. Verstraete, Unifying time evolution and optimization with matrix product states.  Phys. Rev. B 94, 165,116 (2016).  10.1103/PhysRevB.94.165116
  * (10) R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  Annals of Physics 349, 117–158 (2014).  https://doi.org/10.1016/j.aop.2014.06.013.  URL https://www.sciencedirect.com/science/article/pii/S0003491614001596
  * (11) F.A.Y.N. Schröder, 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, 075,105 (2016).  10.1103/PhysRevB.93.075105.  URL https://link.aps.org/doi/10.1103/PhysRevB.93.075105
  * (12) F.A.Y.N. Schröder, D.H.P. Turban, A.J. Musser, N.D.M. Hine, A.W. Chin, Tensor network simulation of multi-environmental open quantum dynamics via machine learning and entanglement renormalisation.  Nature Communications 10, 1062 (2019)
  * (13) A.M. Alvertis, F.A.Y.N. Schröder, A.W. Chin, Non-equilibrium relaxation of hot states in organic semiconductors: Impact of mode-selective excitation on charge transfer.  J. Chem. Phys. 151(8), 084,104 (2019)
  * (14) T. Lacroix, A. Dunnett, D. Gribben, B.W. Lovett, A. Chin, Unveiling non-markovian spacetime signaling in open quantum systems with long-range tensor network dynamics.  Phys. Rev. A 104, 052,204 (2021).  10.1103/PhysRevA.104.052204
  * (15) A.J. Dunnett, D. Gowland, C.M. Isborn, A.W. Chin, T.J. Zuehlsdorff, Influence of non-adiabatic effects on linear absorption spectra in the condensed phase: Methylene blue.  J. Chem. Phys. 155(14), 144,112 (2021).  10.1063/5.0062950
  * (16) A. Baiardi, M. Reiher, Large-scale quantum dynamics with matrix product states.  Journal of Chemical Theory and Computation 15(6), 3481–3498 (2019).  10.1021/acs.jctc.9b00301.  PMID: 31067052
  * (17) X. Xie, Y. Liu, Y. Yao, U. Schollwöck, C. Liu, H. Ma, Time-dependent density matrix renormalization group quantum dynamics for realistic chemical systems.  J. Chem. Phys. 151(22), 224,101 (2019).  10.1063/1.5125945
  * (18) A.D. Somoza, O. Marty, J. Lim, S.F. Huelga, M.B. Plenio, Dissipation-assisted matrix product factorization.  Phys. Rev. Lett. 123, 100,502 (2019).  10.1103/PhysRevLett.123.100502.  URL https://link.aps.org/doi/10.1103/PhysRevLett.123.100502
  * (19) N. Lyu, M.B. Soley, V.S. Batista, Tensor-train split-operator ksl (tt-soksl) method for quantum dynamics simulations.  Journal of Chemical Theory and Computation 18(6), 3327–3346 (2022).  10.1021/acs.jctc.2c00209.  PMID: 35649210
  * (20) P. Gelß, R. Klein, S. Matera, B. Schmidt, Solving the time-independent schrödinger equation for chains of coupled excitons and phonons using tensor trains.  J. Chem. Phys. 156(2), 024,109 (2022).  10.1063/5.0074948.  URL https://doi.org/10.1063/5.0074948.  https://doi.org/10.1063/5.0074948
  * (21) R. Feynman, F. Vernon, The theory of a general quantum system interacting with a linear dissipative system.  Annals of Physics 24, 118–173 (1963).  https://doi.org/10.1016/0003-4916(63)90068-X.  URL https://www.sciencedirect.com/science/article/pii/000349166390068X
  * (22) Y. Tanimura, R. Kubo, Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath.  J. Phys. Soc. Jpn. 58(1), 101–114 (1989).  10.1143/JPSJ.58.101
  * (23) Y. Tanimura, Stochastic liouville, langevin, fokker–planck, and master equation approaches to quantum dissipative systems.  Journal of the Physical Society of Japan 75(8), 082,001 (2006).  10.1143/JPSJ.75.082001
  * (24) Y. Tanimura, Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (heom).  J. Chem. Phys. 153(2), 020,901 (2020).  10.1063/5.0011599
  * (25) R.X. Xu, Y. Yan, Dynamics of quantum dissipation systems interacting with bosonic canonical bath: Hierarchical equations of motion approach.  Phys. Rev. E 75, 031,107 (2007)
  * (26) Q. Shi, L. Chen, G. Nan, R.X. Xu, Y. Yan, Efficient hierarchical Liouville space propagator to quantum dissipative dynamics.  J. Chem. Phys. 130(8), 084,105 (2009).  10.1063/1.3077918
  * (27) S. Nakajima, On Quantum Theory of Transport Phenomena: Steady Diffusion.  Progress of Theoretical Physics 20(6), 948–959 (1958).  10.1143/PTP.20.948
  * (28) R. Zwanzig, Ensemble method in the theory of irreversibility.  The Journal of Chemical Physics 33(5), 1338–1341 (1960).  10.1063/1.1731409
  * (29) A. Ishizaki, G.R. Fleming, Unified treatment of quantum coherent and incoherent hopping dynamics in electronic energy transfer: Reduced hierarchy equation approach.  J. Chem. Phys. 130(23), 234,111 (2009).  10.1063/1.3155372
  * (30) A. Bose, P.L. Walters, A multisite decomposition of the tensor network path integrals.  J.Chem. Phys. 156(2), 024,101 (2022).  10.1063/5.0073234
  * (31) Q. Shi, Y. Xu, Y. Yan, M. Xu, Efficient propagation of the hierarchical equations of motion using the matrix product state method.  J. Chem. Phys. 148(17), 174,102 (2018).  10.1063/1.5026753
  * (32) Y. Yan, T. Xing, 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(20), 204,109 (2020).  10.1063/5.0027962
  * (33) R. Borrelli, S. Dolgov, Expanding the range of hierarchical equations of motion by tensor-train implementation.  J. Phys. Chem. B 125(20), 5397–5407 (2021).  10.1021/acs.jpcb.1c02724
  * (34) R. Borrelli, M.F. Gelin, Finite temperature quantum dynamics of complex systems: Integrating thermo-field theories and tensor-train methods.  WIREs Computational Molecular Science 11(6), e1539 (2021).  https://doi.org/10.1002/wcms.1539
  * (35) Y. Yan, M. Xu, T. Li, Q. Shi, Efficient propagation of the hierarchical equations of motion using the tucker and hierarchical tucker tensors.  J. Chem. Phys. 154(19), 194,104 (2021).  10.1063/5.0050720
  * (36) A. Ishizaki, T.R. Calhoun, G.S. Schlau-Cohen, G.R. Fleming, Quantum coherence and its interplay with protein environments in photosynthetic electronic energy transfer.  Phys. Chem. Chem. Phys. 12, 7319–7337 (2010).  10.1039/C003389H.  URL http://dx.doi.org/10.1039/C003389H
  * (37) C. Kreisbeck, T. Kramer, Long-lived electronic coherence in dissipative exciton dynamics of light-harvesting complexes.  J. Phys. Chem. Lett. 3(19), 2828–2833 (2012).  10.1021/jz3012029.  URL http://pubs.acs.org/doi/abs/10.1021/jz3012029.  arXiv:1203.1485
  * (38) S. Saito, M. Higashi, G.R. Fleming, Site-dependent fluctuations optimize electronic energy transfer in the fenna–matthews–olson protein.  The Journal of Physical Chemistry B 123(46), 9762–9772 (2019).  10.1021/acs.jpcb.9b07456.  URL https://doi.org/10.1021/acs.jpcb.9b07456.  PMID: 31657928
  * (39) L. Chen, P. Shenai, F. Zheng, A. Somoza, Y. Zhao, Optimal energy transfer in light-harvesting systems.  Molecules 20(8), 15,224–15,272 (2015).  10.3390/molecules200815224.  URL https://www.mdpi.com/1420-3049/20/8/15224
  * (40) T.P. Fay, D.T. Limmer, Coupled charge and energy transfer dynamics in light harvesting complexes from a hybrid hierarchical equations of motion approach.  The Journal of Chemical Physics 157(17), 174,104 (2022).  10.1063/5.0117659.  URL https://doi.org/10.1063/5.0117659.  https://doi.org/10.1063/5.0117659
  * (41) M. Cainelli, Y. Tanimura, Exciton transfer in organic photovoltaic cells: A role of local and nonlocal electron–phonon interactions in a donor domain.  J. Chem. Phys. 154(3), 034,107 (2021).  10.1063/5.0036590
  * (42) M. Tanaka, Y. Tanimura, Multistate electron transfer dynamics in the condensed phase: Exact calculations from the reduced hierarchy equations of motion approach.  J. Chem. Phys. 132(21), 214,502 (2010).  10.1063/1.3428674
  * (43) E. Mangaud, A. de la Lande, C. Meier, M. Desouter-Lecomte, Electron transfer within a reaction path model calibrated by constrained dft calculations: application to mixed-valence organic compounds.  Phys. Chem. Chem. Phys. 17, 30,889–30,903 (2015).  10.1039/C5CP01194A.  URL http://dx.doi.org/10.1039/C5CP01194A
  * (44) T. Firmino, E. Mangaud, F. Cailliez, A. Devolder, D. Mendive-Tapia, F. Gatti, C. Meier, M. Desouter-Lecomte, A. de la Lande, Quantum effects in ultrafast electron transfers within cryptochromes.  Phys. Chem. Chem. Phys. 18, 21,442–21,457 (2016).  10.1039/C6CP02809H.  URL http://dx.doi.org/10.1039/C6CP02809H
  * (45) L. Chen, M.F. Gelin, V.Y. Chernyak, W. Domcke, Y. Zhao, Dissipative dynamics at conical intersections: Simulations with the hierarchy equations of motion method.  Faraday Discuss. 194(0), 61–80 (2016).  10.1039/C6FD00088F
  * (46) H.G. Duan, V.I. Prokhorenko, R.J. Cogdell, K. Ashraf, A.L. Stevens, M. Thorwart, R.D. Miller, Nature does not rely on long-lived electronic quantum coherence for photosynthetic energy transfer.  Proceedings of the National Academy of Sciences 114(32), 8493–8498 (2017)
  * (47) A.G. Dijkstra, V.I. Prokhorenko, Simulation of photo-excited adenine in water with a hierarchy of equations of motion approach.  J. Chem. Phys. 147(6), 064,102 (2017).  10.1063/1.4997433
  * (48) E. Mangaud, B. Lasorne, O. Atabek, M. Desouter-Lecomte, Statistical distributions of the tuning and coupling collective modes at a conical intersection using the hierarchical equations of motion.  J. Chem. Phys. 151(24), 244,102 (2019).  10.1063/1.5128852
  * (49) G. Breuil, E. Mangaud, B. Lasorne, O. Atabek, M. Desouter-Lecomte, Funneling dynamics in a phenylacetylene trimer: Coherent excitation of donor excitonic states and their superposition.  J. Chem. Phys. 155(3), 034,303 (2021).  10.1063/5.0056351
  * (50) A. Jaouadi, J. Galiana, E. Mangaud, B. Lasorne, O. Atabek, M. Desouter-Lecomte, Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control.  Phys. Rev. A 106, 043,121 (2022).  10.1103/PhysRevA.106.043121.  URL https://link.aps.org/doi/10.1103/PhysRevA.106.043121
  * (51) J. Zhang, R. Borrelli, Y. Tanimura, Proton tunneling in a two-dimensional potential energy surface with a non-linear system–bath interaction: Thermal suppression of reaction rate.  The Journal of Chemical Physics 152(21), 214,114 (2020).  10.1063/5.0010580.  URL https://doi.org/10.1063/5.0010580.  https://doi.org/10.1063/5.0010580
  * (52) E. Mangaud, R. Puthumpally-Joseph, D. Sugny, C. Meier, O. Atabek, M. Desouter-Lecomte, Non-markovianity in the optimal control of an open quantum system described by hierarchical equations of motion.  New J. Phys. 20, 043,050 (2018)
  * (53) A. Kato, Y. Tanimura, Quantum heat current under non-perturbative and non-markovian conditions: Applications to heat machines.  The Journal of Chemical Physics 145(22), 224,105 (2016).  10.1063/1.4971370.  URL https://doi.org/10.1063/1.4971370.  https://doi.org/10.1063/1.4971370
  * (54) L. Song, Q. Shi, Hierarchical equations of motion method applied to nonequilibrium heat transport in model molecular junctions: Transient heat current and high-order moments of the current operator.  Phys. Rev. B 95(6), 064,308 (2017).  10.1103/PhysRevB.95.064308
  * (55) J. Bätge, Y. Ke, C. Kaspar, M. Thoss, Nonequilibrium open quantum systems with multiple bosonic and fermionic environments: A hierarchical equations of motion approach.  Phys. Rev. B 103, 235,413 (2021).  10.1103/PhysRevB.103.235413.  URL https://link.aps.org/doi/10.1103/PhysRevB.103.235413
  * (56) E.C. Wu, Q. Ge, E.A. Arsenault, N.H.C. Lewis, N.L. Gruenke, M.J. Head-Gordon, G.R. Fleming, Two-dimensional electronic-vibrational spectroscopic study of conical intersection dynamics: An experimental and electronic structure study.  Phys. Chem. Chem. Phys. 21(26), 14,153–14,163 (2019).  10.1039/C8CP05264F
  * (57) K.B. Zhu, R.X. Xu, H.Y. Zhang, J. Hu, Y.J. Yan, Hierarchical dynamics of correlated system-environment coherence and optical spectroscopy.  The Journal of Physical Chemistry B 115(18), 5678–5684 (2011).  10.1021/jp2002244.  URL https://doi.org/10.1021/jp2002244.  PMID: 21452824.  https://doi.org/10.1021/jp2002244
  * (58) Y. Tanimura, Reduced hierarchy equations of motion approach with drude plus brownian spectral distribution: Probing electron transfer processes by means of two-dimensional correlation spectroscopy.  The Journal of Chemical Physics 137(22), 22A550 (2012).  10.1063/1.4766931.  URL https://doi.org/10.1063/1.4766931.  https://doi.org/10.1063/1.4766931
  * (59) H. Liu, L. Zhu, S. Bai, Q. Shi, Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes.  J. Chem. Phys. 140(13), 134,106 (2014).  10.1063/1.4870035
  * (60) A.J. Dunnett, A.W. Chin, Simulating quantum vibronic dynamics at finite temperatures with many body wave functions at 0 k.  Frontiers in Chemistry 8 (2021).  10.3389/fchem.2020.600731.  URL https://www.frontiersin.org/articles/10.3389/fchem.2020.600731
  * (61) D. Tamascelli, A. Smirne, J. Lim, S.F. Huelga, M.B. Plenio, Efficient simulation of finite-temperature open quantum systems.  Phys. Rev. Lett. 123, 090,402 (2019).  10.1103/PhysRevLett.123.090402.  URL https://link.aps.org/doi/10.1103/PhysRevLett.123.090402
  * (62) L. Zhu, H. Liu, W. Xie, Q. Shi, Explicit system-bath correlation calculated using the hierarchical equations of motion method.  J. Chem. Phys. 137(19), 194,106 (2012).  10.1063/1.4766358
  * (63) A. Chin, E. Mangaud, V. Chevet, O. Atabek, M. Desouter-Lecomte, Visualising the role of non-perturbative environment dynamics in the dissipative generation of coherent electronic motion.  Chem. Phys. 525, 110,392 (2019).  https://doi.org/10.1016/j.chemphys.2019.110392
  * (64) H.P. Breuer, F. Petruccione, _The Theory of Open Quantum Systems_ (Oxford University Press, 2002)
  * (65) U. Weiss, _Quantum Dissipative Systems_ (Singapore, World Scientific, 2012)
  * (66) V. May, O. Kühn, _Charge and Energy Transfer Dynamics in Molecular Systems_ (John Wiley and Sons, Ltd, 2011)
  * (67) A. Chenel, E. Mangaud, I. Burghardt, C. Meier, M. Desouter-Lecomte, Exciton dissociation at donor-acceptor heterojunctions: Dynamics using the collective effective mode representation of the spin-boson model.  J. Chem. Phys. 140(4), 044,104 (2014).  10.1063/1.4861853
  * (68) J. Iles-Smith, A.G. Dijkstra, N. Lambert, A. Nazir, Energy transfer in structured and unstructured environments: Master equations beyond the Born-Markov approximations.  J. Chem. Phys. 144(4), 044,110 (2016).  10.1063/1.4940218
  * (69) A. Pomyalov, C. Meier, D.J. Tannor, The importance of initial correlations in rate dynamics: A consistent non-Markovian master equation approach.  Chem. Phys. 370(1–3), 98–108 (2010).  10.1016/j.chemphys.2010.02.017
  * (70) A.G. Dijkstra, Y. Tanimura, Non-markovian entanglement dynamics in the presence of system-bath coherence.  Phys. Rev. Lett. 104, 250,401 (2010).  10.1103/PhysRevLett.104.250401
  * (71) Y. Tanimura, Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities.  J. Chem. Phys. 141(4), 044,114 (2014).  10.1063/1.4890441
  * (72) L. Song, Q. Shi, Calculation of correlated initial state in the hierarchical equations of motion method using an imaginary time path integral approach.  J. Chem. Phys. 143(19), 194,106 (2015).  10.1063/1.4935799
  * (73) C. Meier, D.J. Tannor, Non-Markovian evolution of the density operator in the presence of strong laser fields.  J. Chem. Phys. 111(8), 3365–3376 (1999).  10.1063/1.479669
  * (74) S. Jang, J. Cao, R.J. Silbey, Fourth-order quantum master equation and its markovian bath limit.  J. Chem. Phys. 116(7), 2705–2717 (2002).  10.1063/1.1445105
  * (75) E. Mulvihill, A. Schubert, X. Sun, B.D. Dunietz, E. Geva, A modified approach for simulating electronically nonadiabatic dynamics via the generalized quantum master equation.  The Journal of Chemical Physics 150(3), 034,101 (2019).  10.1063/1.5055756.  URL https://doi.org/10.1063/1.5055756.  https://doi.org/10.1063/1.5055756
  * (76) S. Mukamel, _Principles of nonlinear optical spectroscopy_ (Oxford University Press, 1995)
  * (77) H. Rahman, U. Kleinekathöfer, Chebyshev hierarchical equations of motion for systems with arbitrary spectral densities and temperatures.  J. Chem. Phys. 150(24), 244,104 (2019)
  * (78) L. Cui, H.D. Zhang, X. Zheng, R.X. Xu, Y. Yan, Highly efficient and accurate sum-over-poles expansion of Fermi and Bose functions at near zero temperatures: Fano spectrum decomposition scheme.  J. Chem. Phys. 151(2), 024,110 (2019)
  * (79) T. Ikeda, G.D. Scholes, Generalization of the hierarchical equations of motion theory for efficient calculations with arbitrary correlation functions.  J. Chem. Phys. 152(20), 204,101 (2020).  10.1063/5.0007327
  * (80) Z.H. Chen, Y. Wang, X. Zheng, R.X. Xu, Y. Yan, Universal time-domain prony fitting decomposition for optimized hierarchical quantum master equations.  J. Chem. Phys. 156(22), 221,102 (2022).  10.1063/5.0095961
  * (81) N. Lambert, T. Raheja, S. Cross, P. Menczel, S. Ahmed, A. Pitchford, D. Burgarth, F. Nori.  Qutip-bofin: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics (2020).  10.48550/ARXIV.2010.10806.  URL https://arxiv.org/abs/2010.10806
  * (82) U. Kleinekathöfer, Non-markovian theories based on a decomposition of the spectral density.  J. Chem. Phys. 121(6), 2505–2514 (2004).  10.1063/1.1770619
  * (83) S. Valleau, A. Eisfeld, A. Aspuru-Guzik, On the alternatives for bath correlators and spectral densities from mixed quantum-classical simulations.  The Journal of Chemical Physics 137(22), 224,103 (2012).  10.1063/1.4769079.  URL https://doi.org/10.1063/1.4769079.  https://doi.org/10.1063/1.4769079
  * (84) T. Renger, R.A. Marcus, On the relation of protein dynamics and exciton relaxation in pigment–protein complexes: An estimation of the spectral density and a theory for the calculation of optical spectra.  The Journal of Chemical Physics 116(22), 9997–10,019 (2002).  10.1063/1.1470200.  URL https://doi.org/10.1063/1.1470200.  https://doi.org/10.1063/1.1470200
  * (85) G. Ritschel, A. Eisfeld, Analytic representations of bath correlation functions for ohmic and superohmic spectral densities using simple poles.  The Journal of Chemical Physics 141(9), 094,101 (2014).  10.1063/1.4893931.  URL https://doi.org/10.1063/1.4893931.  https://doi.org/10.1063/1.4893931
  * (86) J. Hu, M. Luo, F. Jiang, R.X. Xu, Y. Yan, Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems.  J. Chem. Phys. 134(24), 244,106 (2011).  10.1063/1.3602466
  * (87) N. Lambert, S. Ahmed, M. Cirio, F. Nori, Modelling the ultra-strongly coupled spin-boson model with unphysical modes.  Nature Communications 10(1), 3721 (2019).  10.1038/s41467-019-11656-1.  https://doi.org/10.1038/s41467-019-11656-1
  * (88) T.P. Fay, A simple improved low temperature correction for the hierarchical equations of motion.  The Journal of Chemical Physics 157(5), 054,108 (2022).  10.1063/5.0100365.  URL https://doi.org/10.1063/5.0100365.  https://doi.org/10.1063/5.0100365
  * (89) E. Mangaud, C. Meier, M. Desouter-Lecomte, Analysis of the non-markovianity for electron transfer reactions in an oligothiophene-fullerene heterojunction.  Chem. Phys. 494, 90–102 (2017).  https://doi.org/10.1016/j.chemphys.2017.07.011.  URL https://www.sciencedirect.com/science/article/pii/S0301010417304548
  * (90) R. Kubo, Stochastic Liouville Equations.  Journal of Mathematical Physics 4(2), 174 (1963).  10.1063/1.1703941
  * (91) A. Ishizaki, Y. Tanimura, Quantum Dynamics of System Strongly Coupled to Low-Temperature Colored Noise Bath: Reduced Hierarchy Equations Approach.  J. Phys. Soc. Jpn. 74(12), 3131–3134 (2005).  10.1143/JPSJ.74.3131
  * (92) G.C. Wick, The evaluation of the collision matrix.  Phys. Rev. 80, 268–272 (1950).  10.1103/PhysRev.80.268.  URL https://link.aps.org/doi/10.1103/PhysRev.80.268
  * (93) H. Tamura, R. Martinazzo, M. Ruckenbauer, I. Burghardt, Quantum dynamics of ultrafast charge transfer at an oligothiophene-fullerene heterojunction.  J. Chem. Phys. 137(22), 22A540 (2012).  10.1063/1.4751486
  * (94) F. Di Maiolo, D. Brey, R. Binder, I. Burghardt, Quantum dynamical simulations of intra-chain exciton diffusion in an oligo (para-phenylene vinylene) chain at finite temperature.  J. Chem. Phys. 153(18), 184,107 (2020).  10.1063/5.0027588
  * (95) W. Popp, D. Brey, R. Binder, I. Burghardt, Quantum dynamics of exciton transport and dissociation in multichromophoric systems.  Annual Review of Physical Chemistry 72(1), 591–616 (2021).  10.1146/annurev-physchem-090419-040306.  URL https://doi.org/10.1146/annurev-physchem-090419-040306
  * (96) J. Schulze, M.F. Shibl, M.J. Al-Marri, O. Kühn, Multi-layer multi-configuration time-dependent hartree (ml-mctdh) approach to the correlated exciton-vibrational dynamics in the fmo complex.  The Journal of Chemical Physics 144(18), 185,101 (2016).  10.1063/1.4948563.  URL https://doi.org/10.1063/1.4948563
  * (97) C. Kreisbeck, T. Kramer, A. Aspuru-Guzik, Scalable high-performance algorithm for the simulation of exciton dynamics. application to the light-harvesting complex II in the presence of resonant vibrational modes.  J. Chem. Theory Comput. 10(9), 4045–4054 (2014).  10.1021/ct500629s.  https://doi.org/10.1021/ct500629s
  * (98) R. Martinazzo, K.H. Hughes, F. Martelli, I. Burghardt, Effective spectral densities for system-environment dynamics at conical intersections: S2–S1 conical intersection in pyrazine.  Chem. Phys. 377(1), 21 – 29 (2010)
  * (99) J. Iles-Smith, N. Lambert, A. Nazir, Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems.  Phys. Rev. A 90(3), 032,114 (2014).  10.1103/PhysRevA.90.032114
  * (100) H.P. Breuer, E.M. Laine, J. Piilo, B. Vacchini, Colloquium: Non-markovian dynamics in open quantum systems.  Rev. Mod. Phys. 88, 021,002 (2016).  10.1103/RevModPhys.88.021002.  URL https://link.aps.org/doi/10.1103/RevModPhys.88.021002
  * (101) I. de Vega, D. Alonso, Dynamics of non-Markovian open quantum systems.  Rev. Mod. Phys. 89(1), 015,001 (2017).  10.1103/RevModPhys.89.015001
  * (102) H.P. Breuer, E.M. Laine, J. Piilo, Measure for the degree of non-markovian behavior of quantum processes in open systems.  Phys. Rev. Lett. 103, 210,401 (2009).  10.1103/PhysRevLett.103.210401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.103.210401
  * (103) A. Rivas, S.F. Huelga, M.B. Plenio, Entanglement and non-markovianity of quantum evolutions.  Phys. Rev. Lett. 105, 050,403 (2010).  10.1103/PhysRevLett.105.050403.  URL https://link.aps.org/doi/10.1103/PhysRevLett.105.050403
  * (104) M.J.W. Hall, J.D. Cresser, L. Li, E. Andersson, Canonical form of master equations and characterization of non-markovianity.  Phys. Rev. A 89, 042,120 (2014).  10.1103/PhysRevA.89.042120.  URL https://link.aps.org/doi/10.1103/PhysRevA.89.042120
  * (105) S. Lorenzo, F. Plastina, M. Paternostro, Geometrical characterization of non-markovianity.  Phys. Rev. A 88, 020,102 (2013).  10.1103/PhysRevA.88.020102.  URL https://link.aps.org/doi/10.1103/PhysRevA.88.020102
  * (106) S.C. Hou, S.L. Liang, X.X. Yi, Non-markovianity and memory effects in quantum open systems.  Phys. Rev. A 91, 012,109 (2015).  10.1103/PhysRevA.91.012109.  URL https://link.aps.org/doi/10.1103/PhysRevA.91.012109
  * (107) G. Kimura, The Bloch vector for N-level systems.  Phys. Lett. A 314(5), 339 – 349 (2003).  https://doi.org/10.1016/S0375-9601(03)00941-1
  * (108) D. Aerts, M.S. de Bianchi, The extended Bloch representation of quantum mechanics and the hidden-measurement solution to the measurement problem.  Ann. Phys. 351, 975 – 1025 (2014).  https://doi.org/10.1016/j.aop.2014.09.020
  * (109) R. Puthumpally-Joseph, O. Atabek, E. Mangaud, M. Desouter-Lecomte, D. Sugny, Towards laser control of open quantum systems: memory effects.  Molecular Physics 115(15-16), 1944–1954 (2017).  10.1080/00268976.2017.1319085.  URL https://doi.org/10.1080/00268976.2017.1319085.  https://doi.org/10.1080/00268976.2017.1319085
  * (110) R. Puthumpally-Joseph, E. Mangaud, V. Chevet, M. Desouter-Lecomte, D. Sugny, O. Atabek, Basic mechanisms in the laser control of non-markovian dynamics.  Phys. Rev. A 97, 033,411 (2018).  10.1103/PhysRevA.97.033411.  URL https://link.aps.org/doi/10.1103/PhysRevA.97.033411
  * (111) C.A. Rodríguez-Rosario, K. Modi, A. Aspuru-Guzik, Linear assignment maps for correlated system-environment states.  Phys. Rev. A 81, 012,313 (2010).  10.1103/PhysRevA.81.012313.  URL https://link.aps.org/doi/10.1103/PhysRevA.81.012313
  * (112) A. Brodutch, A. Datta, K. Modi, A. Rivas, C.A. Rodríguez-Rosario, Vanishing quantum discord is not necessary for completely positive maps.  Phys. Rev. A 87, 042,301 (2013).  10.1103/PhysRevA.87.042301.  URL https://link.aps.org/doi/10.1103/PhysRevA.87.042301
  * (113) M.D. Choi, Completely positive linear maps on complex matrices.  Linear Algebra and its Applications 10(3), 285–290 (1975).  https://doi.org/10.1016/0024-3795(75)90075-0.  URL https://www.sciencedirect.com/science/article/pii/0024379575900750
  * (114) S.H. Kye, Choi matrices revisited.  Journal of Mathematical Physics 63(9), 092,202 (2022).  10.1063/5.0107646.  URL https://doi.org/10.1063/5.0107646.  https://doi.org/10.1063/5.0107646
  * (115) I.S. Dunn, R. Tempelaar, D.R. Reichman, Removing instabilities in the hierarchical equations of motion: Exact and approximate projection approaches.  The Journal of Chemical Physics 150(18), 184,109 (2019).  10.1063/1.5092616.  URL https://doi.org/10.1063/1.5092616.  https://doi.org/10.1063/1.5092616
  * (116) C.P. Koch, Controlling open quantum systems: tools, achievements, and limitations.  Journal of Physics: Condensed Matter 28(21), 213,001 (2016).  10.1088/0953-8984/28/21/213001.  URL https://doi.org/10.1088/0953-8984/28/21/213001
  * (117) B. Bylicka, D. Chruściński, S. Maniscalco.  Non-markovianity as a resource for quantum technologies (2013).  10.48550/ARXIV.1301.2585.  URL https://arxiv.org/abs/1301.2585
  * (118) N. Mirkin, P. Poggi, D. Wisniacki, Information backflow as a resource for entanglement.  Physical Review A 99(6) (2019).  10.1103/physreva.99.062327.  URL https://doi.org/10.1103%2Fphysreva.99.062327
  * (119) C.P. Koch, U. Boscain, T. Calarco, G. Dirr, S. Filipp, S.J. Glaser, R. Kosloff, S. Montangero, T. Schulte-Herbrüggen, D. Sugny, F.K. Wilhelm, Quantum optimal control in quantum technologies. strategic report on current status, visions and goals for research in europe.  EPJ Quantum Technology 9(1), 19 (2022).  10.1140/epjqt/s40507-022-00138-x.  URL https://doi.org/10.1140/epjqt/s40507-022-00138-x
  * (120) C. Lubich, I.V. Oseledets, A projector-splitting integrator for dynamical low-rank approximation.  Bit Numer Math 54, 171–188 (2014).  https://doi.org/10.1007/s10543-013-0454-0
  * (121) C. Lubich, I.V. Oseledets, B. Vandereycken, Time integration of tensor trains.  SIAM Journal on Numerical Analysis 53(2), 917–941 (2015).  10.1137/140976546.  URL https://doi.org/10.1137%2F140976546
  * (122) https://github.com/oseledets/ttpy
  * (123) A.J. Dunnett, 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, 214,302 (2021).  10.1103/PhysRevB.104.214302.  URL https://link.aps.org/doi/10.1103/PhysRevB.104.214302
  * (124) S.V. Dolgov, A tensor decomposition algorithm for large odes with conservation laws.  Comput. Methods Appl. Math. 19(1), 23–38 (2019).  10.1515/cmam-2018-0023.  URL https://doi.org/10.1515/cmam-2018-0023
  * (125) A. Potočnik, A. Bargerbos, F.A.Y.N. Schröder, S.A. Khan, M.C. Collodo, S. Gasparinetti, Y. Salathé, C. Creatore, C. Eichler, H.E. Türeci, A.W. Chin, A. Wallraff, Studying light-harvesting models with superconducting circuits.  Nature Communications 9(1), 1–7 (2018).  10.1038/s41467-018-03312-x.  URL https://doi.org/10.1038/s41467-018-03312-x.  https://doi.org/10.1038/s41467-018-03312-x
  * (126) A.W. Chin, E. Mangaud, O. Atabek, M. Desouter-Lecomte, Coherent quantum dynamics launched by incoherent relaxation in a quantum circuit simulator of a light-harvesting complex.  Phys. Rev. A 97, 063,823 (2018).  10.1103/PhysRevA.97.063823.  URL https://link.aps.org/doi/10.1103/PhysRevA.97.063823
  * (127) P.R. Eastham, P. Kirton, H.M. Cammack, B.W. Lovett, J. Keeling, Bath-induced coherence and the secular approximation.  Phys. Rev. A 94, 012,110 (2016).  10.1103/PhysRevA.94.012110.  URL https://link.aps.org/doi/10.1103/PhysRevA.94.012110
  * (128) T. Brabec, F. Krausz, Intense few-cycle laser fields: Frontiers of nonlinear optics.  Rev. Mod. Phys. 72, 545–591 (2000).  10.1103/RevModPhys.72.545.  URL https://link.aps.org/doi/10.1103/RevModPhys.72.545
  * (129) D.B. Milošević, G.G. Paulus, D. Bauer, W. Becker, Above-threshold ionization by few-cycle pulses.  J. Phys. B: At. Mol. Opt. Phys. 39(14), R203–R262 (2006).  10.1088/0953-4075/39/14/r01.  URL https://doi.org/10.1088/0953-4075/39/14/r01
  * (130) G.F. Thomas, Validity of the rosen-zener conjecture for gaussian-modulated pulses.  Phys. Rev. A 27, 2744 (1983)
  * (131) M. Holthaus, B. Just, Generalized $\pi$ pulse.  Phys. Rev. A 49, 1950 (1994)
  * (132) Y. Jing, L. Chen, S. Bai, Q. Shi, Equilibrium excited state and emission spectra of molecular aggregates from the hierarchical equations of motion approach.  The Journal of Chemical Physics 138(4), 045,101 (2013).  10.1063/1.4775843.  URL https://doi.org/10.1063/1.4775843.  https://doi.org/10.1063/1.4775843
  * (133) J. Ma, J. Cao, Förster resonance energy transfer, absorption and emission spectra in multichromophoric systems. i. full cumulant expansions and system-bath entanglement.  The Journal of Chemical Physics 142(9), 094,106 (2015).  10.1063/1.4908599.  URL https://doi.org/10.1063/1.4908599.  https://doi.org/10.1063/1.4908599
  * (134) H.D. Zhang, Q. Qiao, R.X. Xu, X. Zheng, Y. Yan, Efficient steady-state solver for hierarchical quantum master equations.  J. Chem. Phys. 147(4), 044,105 (2017).  10.1063/1.4995424.  URL https://doi.org/10.1063/1.4995424.  https://doi.org/10.1063/1.4995424
  * (135) T. Kramer, M. Noack, A. Reinefeld, M. Rodríguez, Y. Zelinskyy, Efficient calculation of open quantum system dynamics and time-resolved spectroscopy with distributed memory heom (dm-heom).  J. Comput. Chem. 39(22), 1779–1794 (2018).  https://doi.org/10.1002/jcc.25354.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.25354.  https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.25354
  * (136) H. Tamura, I. Burghardt, M. Tsukada, Exciton dissociation at thiophene/fullerene interfaces: The electronic structures and quantum dynamics.  The Journal of Physical Chemistry C 115(20), 10,205–10,210 (2011).  10.1021/jp203174e.  URL https://doi.org/10.1021/jp203174e
  * (137) P.L. Walters, T.C. Allen, N. Makri, Direct determination of discrete harmonic bath parameters from molecular dynamics simulations.  Journal of Computational Chemistry 38(2), 110–115 (2017).  https://doi.org/10.1002/jcc.24527.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.24527.  https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.24527
  * (138) F. Häse, C. Kreisbeck, A. Aspuru-Guzik, Machine learning for quantum dynamics: deep learning of excitation energy transfer properties.  Chem. Sci. 8, 8419–8426 (2017).  10.1039/C7SC03542J.  URL http://dx.doi.org/10.1039/C7SC03542J
  * (139) S. Bandyopadhyay, Z. Huang, K. Sun, Y. Zhao, Applications of neural networks to the simulation of dynamics of open quantum systems.  Chemical Physics 515, 272–278 (2018).  https://doi.org/10.1016/j.chemphys.2018.05.019.  URL https://www.sciencedirect.com/science/article/pii/S0301010418304336.  Ultrafast Photoinduced Processes in Polyatomic Molecules:Electronic Structure, Dynamics and Spectroscopy (Dedicated to Wolfgang Domcke on the occasion of his 70th birthday)
  * (140) F.A.Y.N. Schröder, D.H.P. Turban, A.J. Musser, N.D.M. Hine, A.W. Chin, Tensor network simulation of multi-environmental open quantum dynamics via machine learning and entanglement renormalisation.  Nature Communications 10(1), 1–10 (2019).  10.1038/s41467-019-09039-7.  URL https://doi.org/10.1038/s41467-019-09039-7.  https://doi.org/10.1038/s41467-019-09039-7