# Simulation of absorption spectra of molecular aggregates: a Hierarchy of
Stochastic Pure States approach

Lipeng Chen  Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Str 38, Dresden, Germany<EMAIL_ADDRESS>Doran I. G. Bennett
Department of Chemistry, Southern Methodist University, PO Box 750314, Dallas,
TX, USA<EMAIL_ADDRESS>Alexander Eisfeld  Max Planck Institute for the
Physics of Complex Systems, Nöthnitzer Str 38, Dresden, Germany
<EMAIL_ADDRESS>
###### Abstract

The simulation of spectroscopic observables for molecular aggregates with
strong and structured coupling of electronic excitation to vibrational degrees
of freedom is an important but challenging task. The hierarchy of pure states
(HOPS) provides a formally exact solution based on local, stochastic
trajectories. Exploiting the localization of HOPS for the simulation of
absorption spectra in large aggregares requires a formulation in terms of
normalized trajectories. Here we provide a normalized dyadic equation where
the ket\- and bra-states are propagated in different electronic Hilbert
spaces. This work opens the door to apply adaptive HOPS methods for the
simulation of absorption spectra and also to a formulation for non-linear
spectroscopy that is perturbative with respect to interactions with the
electric field.

##  I Introduction

The photophysics of molecular aggregates plays a central role in light
harvesting by both artificial materials and photosynthetic organisms.
ReviewPS1 ; ReviewPS2 ; ReviewPS3 ; ReviewPS4  Optical absorption is a simple,
yet powerful, probe of the delocalized excited states (excitons) formed in
molecular aggregates. The basic molecular exciton theory treats each
chromophore as an electronic two-level (or few-level) system coupled to
molecular vibrations and a continuum of vibrational degrees of freedom which
are responsible for energy dissipation and electronic dephasing (described by
the spectral density).ReviewPS2 ; SD1 ; SD2 ; SD3 ; SD4  When the coupling
between the electronic states of the chromophore and the system-environment
coupling are of similar magnitudes, then perturbative treatments are
inappropriate; when the vibrational modes are highly structured it is
essential to account for non-Markovian effects. From a theoretical point of
view, it is thus of great importance to develop approaches which are capable
of treating exciton dynamics at finite temperature in a non-perturbative and
non-Markovian manner.

While open quantum system methods appropriate to molecular excitons have
typically used density matrices,NPNMRMP ; OQS1 ; OQS2 ; Tanimura2006JPSJ ;
Tanimura2020JCP ; Makri1995JCP1 ; Makri1995JCP2 ; VarMaster1 ; VarMaster2
wave function methods have been gaining popularity. Beside methods that
simultaneously propagate electronic and vibrational degrees of freedom (see
e.g Refs. MCTDH1 ; MCTDH2 ; MCTDH3 ; GMCTDH ; MCG ; DA1 ; DA2 ; MCE1 ; MCE2
), it is also possible to propagate only the relevant electronic degrees of
freedom (or occasionally vibronic) subject to a stochastic force term that is
constructed to ensure the ensemble reproduces the correct density matrix. One
of these stochasic approaches is the non-Markovian quantum state diffusion
(NMQSD).NMQSD1 ; NMQSD2 ; NMQSD3  The evolution equations in NMQSD come in two
flavours: a linear equation, where individual wave functions are not
normalized, and a non-linear one, where each wave function is properly
normalized ensuring more efficient convergence of the ensemble.

The hierarchy of pure states (HOPS) HOPS1 ; HOPS2  is a formally exact way to
handle the general NMQSD equation. In HOPS the non-Markovian influence of the
bath degrees of freedom is captured by a hierarchy of vectors in the system
Hilbert space that evolve subject to a Gaussian stochastic noise that shares a
correlation function with the environment. Averaging over the stochastic noise
trajectories reproduces the complete system dynamics. One appealing feature of
NMQSD and HOPS is that dynamics are localized by interaction with the thermal
environment which allows adaptive algorithms to calculate dynamics in
molecules aggregates composed of thousands of particles.AdpGao ; AdpDoran

The ability to connect exciton dynamics as obtained by HOPS with spectroscopic
observables is essential. For HOPS calculations, the linear absorption
spectrum can be calculated exactly using a single trajectory by setting the
noise trajectory equal to zero. The full cancellation of noise terms, however,
does not directly extend to non-linear spectroscopic calculations.
Furthermore, the zero-noise trajectory does not show dynamic localization and
therefore cannot be solved using efficient adaptive algorithms which exploit
the localization of single trajectories. These constraints greatly limits the
applicability of HOPS for simulating non-linear spectroscopy or even linear
spectra of large molecular aggregates.

Figure 1: Sketch of the dyadic scheme to calculate linear absorption.
Initially the aggregate is in the electronic ground state g, represented by
the density matrix
$\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}$.
Interactions with the electromagnetic field are indicated by incoming and
outgoing arrows. The dipole operator then brings the ket into the excited
state $\bigl{|}\psi_{\mathrm{ex}}\bigr{\rangle}$; the density matrix is in the
state
$\bigl{|}\psi_{\mathrm{ex}}\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}$.
Finally, the dipole operator brings the excited state back to the electronic
ground state, and the trace is taken. For further details about the formalism
see Ref. MukamelTextBook,

Here, we provide a framework for simulating linear absorption using stochastic
HOPS equations. One particular goal is to connect the HOPS approach, where
pure states are propagated to the common formulation of spectroscopy in terms
of dyadic states, where the bra and ket can be in different electronic states
and are propagated separately. Such a dyadic propagation scheme is
symbolically displayed in the form of a so-called double-sided Feynman
diagram, as displayed in Fig. 1. We demonstrate that by expressing the density
matrix after the first interaction with an electric field as a sum over pure
states, we can write a simple equation that is equivalent to propagating the
bra and ket side of the initial density matrix separately according to HOPS.
For a dyadic propagation scheme it is not obvious how to construct the
normalization of the non-linear NMQSD equation. We find an correct
normalization contains contributions from both the bra and ket, coupling the
two time-evolution equations. We study in detail the convergence of the
dipole-dipole correlation function with the number of stochastic trajectories
using both the linear and non-linear HOPS equations. We find that both, the
linear and the non-linear equation are well behaved at the single trajectory
level and show good convergence properties. Remarkably, this hold even in
regimes where the linear HOPS fails catastrophically for the calculation of
site populations.

The paper is organized as follows: In section II we introduce the open quantum
system formulation of the molecular aggregate and provide the NMQSD framework
and the linear and the non-linear HOPS equations. In section III, which
contains our central results, we formulate absorption in terms of stochastic
NMQSD trajectories and we establish a dyadic HOPS formulation. In section IV
we demonstrate the applicability of our formulas and investigate the
convergence with respect to the number of trajectories. Finally, in section V
we summarize our findings and we conclude with an outlook.

Throughout the work we use $\hbar=1$.

##  II Background

###  II.1 Open quantum system description of the aggregate

We consider a molecular aggregate consisting of $N$ interacting chromophores
embedded in a condensed phase environment such as a solvent, protein, or solid
matrix. In the language of open quantum system, the Hamiltonian for the
aggregate can be written as ReviewPS2

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

with $\hat{H}_{\mathrm{S}}$, $\hat{H}_{\mathrm{B}}$, and
$\hat{H}_{\mathrm{SB}}$ denoting the system Hamiltonian, the environment
(bath) Hamiltonian, and the system-bath interaction Hamiltonian, respectively.
The system Hamiltonian $\hat{H}_{\mathrm{S}}$ describing the Frenkel-exciton
reads

$\hat{H}_{\mathrm{S}}=\hat{H}_{\mathrm{g}}+\hat{H}_{\mathrm{ex}}$ (2)

with the ground state and excited state Hamiltonians

$\displaystyle\hat{H}_{\mathrm{g}}$ $\displaystyle=$
$\displaystyle\epsilon_{\mathrm{g}}|\mathrm{g}\rangle\langle{\mathrm{g}}|$ (3)
$\displaystyle\hat{H}_{\mathrm{ex}}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\mathrm{N}}\epsilon_{n}|n\rangle\langle{n}|+\sum_{n\neq{m}}^{\mathrm{N}}V_{nm}|n\rangle\langle{m}|,$
(4)

where $\epsilon_{\mathrm{g}}$ is the energy in the ground state,
$\epsilon_{n}$ is the site energy of the $n$th molecule, $|n\rangle$ is the
state describing the $n$th excited molecule, and $V_{nm}$ is the electronic
coupling between excited states of molecules $n$ and $m$. Each chromophore is
coupled to inter- and intra-molecular vibrations which can be modeled as a
collection of harmonic oscillators. In what follows, we assume each molecule
is coupled to an independent collection of vibrations. We thus have the bath
Hamiltonian

$\hat{H}_{\mathrm{B}}=\sum_{n=1}^{N}\sum_{q}\omega_{nq}\hat{b}_{nq}^{\dagger}\hat{b}_{nq}$
(5)

where $\hat{b}_{nq}^{\dagger}\textrm{ }(\hat{b}_{nq})$ are the creation
(annihilation) operator of $q$th mode of chromophore $n$ with frequency
$\omega_{nq}$. The system-bath coupling Hamiltonian $\hat{H}_{\mathrm{SB}}$ is
expressed as

$\hat{H}_{\mathrm{SB}}=-\sum_{n=1}^{N}\hat{L}_{n}\sum_{q}g_{nq}\left(\hat{b}_{nq}^{\dagger}+\hat{b}_{nq}\right)$
(6)

with system coupling operators

$\hat{L}_{n}=\bigl{|}n\bigr{\rangle}\bigl{\langle}n\bigr{|}$ (7)

and $g_{nq}$ is the exciton-bath coupling strength of the $q$th bath mode for
chromophore $n$. The influence of the vibrational modes on the dynamics of the
electronic system is described by the bath-correlation function

$\alpha_{n}(\tau)=\int_{0}^{\infty}\\!\\!\\!\mathrm{d}\omega\,J_{n}(\omega)\Big{(}\coth\big{(}\frac{\beta\omega}{2}\big{)}\cos(\omega\tau)-i\sin(\omega\tau)\Big{)}$
(8)

that contains the spectral density
$J_{n}(\omega)=\sum_{q}|g_{nq}|^{2}\delta(\omega-\omega_{nq})$ and the inverse
temperature $\beta=1/T$.

###  II.2 Non-Markovian quantum state diffusion (NMQSD) and the Hierarchy of
Pure States (HOPS)

For a separable initial state

$\hat{\rho}_{\mathrm{ini}}=\bigl{|}\psi_{\mathrm{ini}}\bigr{\rangle}\bigl{\langle}\psi_{\mathrm{ini}}\bigr{|}\otimes\hat{\rho}_{\mathrm{B}}$
(9)

the expectation value of a system operator $\hat{F}$ can be expressed as

$\langle\hat{F}(t)\rangle={\mathrm{Tr}}[\hat{F}\hat{\rho}(t)]=\mathcal{M}[\bigl{\langle}\psi(t,\mathbf{z}^{*})\bigr{|}\hat{F}\bigl{|}\psi(t,\mathbf{z}^{*})\bigr{\rangle}]$
(10)

where $\mathcal{M}[\cdots]$ denotes ensemble average over realizations of
stochastic trajectories which obey the (linear) non-Markovian quantum state
diffusion (NMQSD) equation NMQSD1 ; NMQSD2

$\begin{split}\partial_{t}|\psi(t,\mathbf{z}^{*})\rangle=&-i\hat{H}_{\mathrm{S}}|\psi(t,\mathbf{z}^{*})\rangle\\\
&+\sum_{n}\hat{L}_{n}z_{t,n}^{*}|\psi(t,\mathbf{z}^{*})\rangle\\\
&-\sum_{n}\hat{L}_{n}^{\dagger}\int_{0}^{t}\mathrm{d}s\,\alpha_{n}(t-s)\frac{\delta{|}\psi(t,\mathbf{z}^{*})\rangle}{\delta{z}_{s,n}^{*}}\end{split}$
(11)

where $|\psi(t,\mathbf{z}^{*})\rangle$ is a vector in the system Hilbert space
and $\mathbf{z}$ comprises complex Gaussian stochastic processes $z_{t,n}^{*}$
with mean $\mathcal{M}[z_{t,n}]=0$, and correlations
$\mathcal{M}[z_{t,n}z_{s,m}]=0$, and
$\mathcal{M}[z_{t,n}z_{s,m}^{*}]=\alpha_{n}(t-s)\delta_{nm}$. The important
non-linear NMQSD equation is obtained by making in the above equation the
following replacements:NMQSD2
$\hat{L}_{n}^{\dagger}\rightarrow\hat{L}_{n}^{\dagger}-\langle\hat{L}_{n}^{\dagger}\rangle_{t}$
and
$\tilde{z}_{t,n}=z_{t,n}^{*}+\int_{0}^{t}\mathrm{d}s\,\alpha_{n}^{*}(t-s)\langle{\hat{L}}_{n}^{\dagger}\rangle_{s}$.
Expectation values $\langle\cdot\rangle_{t}$ are calculated using the
normalized state.

The hierarchy of pure states (HOPS) is a formally exact solution to the open
quantum system model discussed in the previous sections within the NMQSD
framework. HOPS requires the bath correlation function $\alpha_{n}(\tau)$ to
be expanded as a sum of exponentials

$\alpha_{n}(\tau)\approx\sum_{j=1}^{J}p_{nj}e^{-w_{nj}\tau}$ (12)

where $w_{nj}=\gamma_{nj}+i\Omega_{nj}$. Assuming, each molecule has the same
bath correlation function, then the linear HOPS equation is

$\displaystyle\begin{aligned}
\partial_{t}|\psi^{(\mathbf{k})}&(t,\mathbf{z}^{*})\rangle\\\
=&\Big{(}-i\hat{H}_{\mathrm{S}}-\mathbf{k}\cdot\mathbf{w}+\sum_{n}\hat{L}_{n}z_{t,n}^{*}\Big{)}|\psi^{(\mathbf{k})}(t,\mathbf{z}^{*})\rangle\\\
&+\sum_{n}\hat{L}_{n}\sum_{j}k_{nj}p_{nj}|\psi^{(\mathbf{k}-\mathbf{e}_{nj})}(t,\mathbf{z}^{*})\rangle\\\
&-\sum_{n}\hat{L}_{n}^{\dagger}\sum_{j}|\psi^{(\mathbf{k}+\mathbf{e}_{nj})}(t,\mathbf{z}^{*})\rangle\end{aligned}$
(13)

where $\mathbf{w}=\left\\{{w}_{1,1},\cdots,w_{N,J}\right\\}$, and
$\mathbf{k}=\left\\{{k}_{1,1},\cdots,k_{N,J}\right\\}$ with non-negative
integers $k_{nj}$. The appearance of $\mathbf{z}$ in
$|\psi(t,\mathbf{z}^{*})\rangle$ indicates that the state depends on the
complete set of stochastic processes up to time $t$. With ∗ we denote complex
conjugate. The physical wave function is given by
$|\psi^{(\mathbf{0})}(t,\mathbf{z}^{*})\rangle$ while the remaining terms are
‘auxiliary wave functions’ and represent the influence of the finite memory
time of the bath on the time evolution of the electronic state. The time
evolution of the reduced density matrix associated with the system degrees of
freedom is given by
$\hat{\rho}(t)=\mathcal{M}\big{[}|\psi^{(\mathbf{0})}(t,\mathbf{z}^{*})\rangle\langle\psi^{(\mathbf{0})}(t,\mathbf{z}^{*})|\big{]}$.
HOPS consists of infinite set of coupled equations, which can be truncated at
a finite number of hierarchy elements, for example using the triangular
truncation condition for the hierarchy where only auxiliary states with
$|\mathbf{k}|\leq\mathcal{K}$ are taken into account, more flexible scheme to
truncate the HOPS can be found in Ref. TrunHOPS, .

The non-linear version of the HOPS equation reads

$\displaystyle\begin{aligned}
\partial_{t}|\tilde{\psi}^{(\mathbf{k})}(t,\mathbf{z}^{*})\rangle=&\Big{(}-i\hat{H}_{\mathrm{S}}-\mathbf{k}\cdot\mathbf{w}+\sum_{n}\hat{L}_{n}\tilde{z}_{t,n}\Big{)}|\tilde{\psi}^{(\mathbf{k})}(t,\mathbf{z}^{*})\rangle\\\
&+\sum_{n}\hat{L}_{n}\sum_{j}k_{nj}p_{nj}|\tilde{\psi}^{(\mathbf{k}-\mathbf{e}_{nj})}(t,\mathbf{z}^{*})\rangle\\\
&-\sum_{n}\left(\hat{L}_{n}^{\dagger}-\langle{\hat{L}}_{n}^{\dagger}\rangle_{t}\right)\sum_{j}|\tilde{\psi}^{(\mathbf{k}+\mathbf{e}_{nj})}(t,\mathbf{z}^{*})\rangle.\end{aligned}$
(14)

Here,
$\tilde{z}_{t,n}=z_{t,n}^{*}+\int_{0}^{t}\mathrm{d}s\,\alpha_{n}^{*}(t-s)\langle{\hat{L}}_{n}^{\dagger}\rangle_{s}$.
Expectation values $\langle\cdot\rangle_{t}$ are calculated using the
normalized physical state

$\bigl{|}\Psi^{(\mathbf{0})}(t,\mathbf{z}^{*})\bigr{\rangle}=\frac{\bigl{|}\tilde{\psi}^{(\mathbf{0})}(t,\mathbf{z}^{*})\bigr{\rangle}}{\sqrt{\langle\tilde{\psi}^{(\mathbf{0})}(t,\mathbf{z}^{*})|\tilde{\psi}^{(\mathbf{0})}(t,\mathbf{z}^{*})\rangle}}.$
(15)

We emphasize that in particular for strong coupling the non-linear equation is
essential to ensure convergence. Also for adaptivity the non-linear normalized
equation is very beneficial.

##  III Linear spectra from stochastic trajectories

To calculate absorption spectra of the aggregate, we need to extend the
Hamiltonian to include the interaction of the electronic system with the
external electric field MukamelTextBook

$\displaystyle\hat{H}_{\mathrm{L}}$ $\displaystyle=$
$\displaystyle-\sum_{n}^{N}\big{(}(\bm{\mu}_{n}\cdot\bm{\varepsilon})E(t)\big{)}|n\rangle\langle{\mathrm{g}}|+h.c.$
(16)

where $\bm{\mu}_{n}$ is the transition dipole moment of chromophore $n$, and
$\bm{\varepsilon}$ and $E(t)$ are the polarization and the envelope of the
laser pulse, which we have taken to be constant over the size of the
aggregate. Our results can be easily extended to electric field that vary from
molecule to molecule, as needed for example for circular dichroism FTODTFD  or
for near-field spectroscopy.NFS1 ; NFS2

The linear absorption spectra ($\mathcal{A}(\omega)$) is given by the Fourier
transformation

$\mathcal{A}(\omega)=\mathrm{Re}\int_{0}^{\infty}\mathrm{d}t\,e^{i\omega{t}}C(t)$
(17)

of the dipole-dipole auto-correlation function MukamelTextBook

$\displaystyle
C(t)=\mathrm{Tr}\left\\{\hat{\mu}_{\mathrm{eff}}\,{e}^{-i\hat{H}t}\big{(}\hat{\mu}_{\mathrm{eff}}{|}\mathrm{g}\rangle\langle{\mathrm{g}}|\otimes\hat{\rho}_{\mathrm{B}}\big{)}e^{i\hat{H}t}\right\\}$
(18)

where $\hat{\mu}_{\mathrm{eff}}$ is the scalar, collective dipole moment
operator

$\hat{\mu}_{\mathrm{eff}}=\sum_{n}^{N}(\bm{\mu}_{n}\cdot\bm{\varepsilon})\,|n\rangle\langle{\mathrm{g}}|+h.c.$
(19)

and
$\hat{\rho}_{0}=|\mathrm{g}\rangle\langle{\mathrm{g}}|\otimes\hat{\rho}_{\mathrm{B}}$
is a factorized initial total density matrix with
$\hat{\rho}_{\mathrm{B}}=e^{-\beta{\hat{H}}_{\mathrm{B}}}/\mathrm{Tr}_{\mathrm{B}}\left\\{{e}^{-\beta{\hat{H}}_{\mathrm{B}}}\right\\}$
being the density matrix of the thermal bath.

###  III.1 Pure state decomposition of initial density matrix

The NMQSD framework cannot directly applied to the expression (Eq. 18), since
the initial state
$\hat{\mu}_{\mathrm{eff}}\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}$
is not pure. To construct the correlation function in terms of pure states, we
follow the approach of Hartmann and Strunz.HOPSEDRichard  We find the
decomposition

$\hat{\mu}_{\mathrm{eff}}\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}=\frac{1}{2}\,\mu_{\mathrm{tot}}\sum_{\eta\in\\{\pm
1,\pm i\\}}\eta\ \bigl{|}v_{\eta}\bigr{\rangle}\bigl{\langle}v_{\eta}\bigr{|}$
(20)

into pure states given by

$\bigl{|}v_{\eta}\bigr{\rangle}=\frac{1}{\sqrt{2}}(\bigl{|}\psi_{\mathrm{ex}}\bigr{\rangle}+\eta\bigl{|}\mathrm{g}\bigr{\rangle}),\quad\quad\eta\in\\{\pm
1,\pm i\\}$ (21)

where

$\bigl{|}\psi_{\mathrm{ex}}\bigr{\rangle}=\frac{1}{{\mu}_{\mathrm{tot}}}\hat{\mu}_{\mathrm{eff}}\bigl{|}\mathrm{g}\bigr{\rangle}=\frac{1}{{\mu}_{\mathrm{tot}}}\sum_{n=1}^{N}(\bm{\mu}_{n}\cdot\bm{\varepsilon})\bigl{|}n\bigr{\rangle}$
(22)
$\mu_{\mathrm{tot}}=\sqrt{\sum_{n=1}^{N}(\bm{\mu}_{n}\cdot\bm{\varepsilon})^{2}}$
(23)

and $\mu_{\mathrm{tot}}$ is the total transition strength.

With this pure state decomposition, the dipole-dipole auto-correlation
function can be written as a sum over contributions from each pure state and
accordingly each pure state contribution can be treated within the NMQSD
framework. Using Eq. (10) one finds

$C(t)=\frac{1}{2}\mu_{\mathrm{tot}}\sum_{\eta\in\\{\pm 1,\pm
i\\}}\eta\,\mathcal{M}[\bigl{\langle}v_{\eta}(t,z^{*})\bigr{|}\hat{\mu}_{\mathrm{eff}}\bigl{|}v_{\eta}(t,z^{*})\bigr{\rangle}]$
(24)

where $\bigl{|}v_{\eta}(t,z^{*})\bigr{\rangle}$ denotes the time evolved
states $\bigl{|}v_{\eta}\bigr{\rangle}$ according to the NMQSD equation in the
full Hilbert space.

###  III.2 Direct propagation of the dyadic equation

The correlation function can be calculated by directly propagating a single
initial state defined in the one exciton manifold. We will demonstrate that
this method is equivalent to propagating a dyadic equation where the bra and
ket states are independently propagated to construct the correlation function.

First, we will simplify the time evolution of the pure state components of the
correlation function to separate the excited state and ground state dynamics.
We note that the Hamiltonian (Eq. 1) does not couple the electronic ground
state $\bigl{|}\mathrm{g}\bigr{\rangle}$ and the single exciton states
$\bigl{|}n\bigr{\rangle}$: the system Hamiltonian has the form
$\hat{H}_{\mathrm{S}}=\hat{H}_{\mathrm{g}}+\hat{H}_{\mathrm{ex}}$ and the
system-bath coupling Hamiltonian acts only in the one-exciton space.
Therefore, initial states of the form of Eq. 21 will evolve under the NMQSD
equation to give

$\bigl{|}v_{\eta}(t,z^{*})\bigr{\rangle}=\frac{1}{\sqrt{2}}\big{(}\bigl{|}\psi_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}+\eta
e^{-i\epsilon_{\mathrm{g}}t}\bigl{|}\mathrm{g}\bigr{\rangle})$ (25)

where $\bigl{|}\psi_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}$ denotes the state
that is obtained upon evolving the initial state
$\bigl{|}\psi_{\mathrm{ex}}\bigr{\rangle}$ with the NMQSD equation restricted
to the one-exciton space (i.e. replacing $\hat{H}_{\mathrm{S}}$ by
$\hat{H}_{\mathrm{ex}}$ in Eq. (11), and correspondingly in the HOPS
equations).

Using this result together with the explicit form of the collective dipole
operator (Eq. (19)) we find that the terms appearing in Eq. (24) have the form

$\bigl{\langle}v_{\eta}(t,z^{*})\bigr{|}\hat{\mu}_{\mathrm{eff}}\bigl{|}v_{\eta}(t,z^{*})\bigr{\rangle}=\frac{1}{2}\eta^{*}\mu_{\mathrm{tot}}\,\langle\psi_{\mathrm{ex}}|\psi_{\mathrm{ex}}(t,z^{*})\rangle
e^{i\epsilon_{\mathrm{g}}t}+h.c.$ (26)

Recognizing that

$\begin{split}\eta\,\bigl{\langle}v_{\eta}(t,z^{*})\bigr{|}&\hat{\mu}_{\mathrm{eff}}\bigl{|}v_{\eta}(t,z^{*})\bigr{\rangle}\\\
=&\phantom{+}\frac{1}{2}\mu_{\mathrm{tot}}\,|\eta|^{2}\langle\psi_{\mathrm{ex}}|\psi_{\mathrm{ex}}(t,z^{*})\rangle
e^{i\epsilon_{\mathrm{g}}t}\\\
&+\frac{1}{2}\mu_{\mathrm{tot}}\,\eta^{2}\langle\psi_{\mathrm{ex}}(t,z^{*})|\psi_{\mathrm{ex}}\rangle
e^{-i\epsilon_{\mathrm{g}}t}\end{split}$ (27)

where $|\eta|^{2}=1$ for all $\eta=\pm 1,\pm i$ then the four contributions to
the correlation function simplify (due to cancellation) to give

$C(t)=\mu_{\mathrm{tot}}^{2}\mathcal{M}[\langle\psi_{\mathrm{ex}}|\psi_{\mathrm{ex}}(t,z^{*})\rangle]e^{i\epsilon_{\mathrm{g}}t}.$
(28)

Eq. 28 is our central result and shows that the dipole correlation function
can be calculated by propagating the initial state
$\bigl{|}\psi_{\mathrm{ex}}\bigr{\rangle}=\sum_{n=1}^{N}\frac{(\bm{\mu}_{n}\cdot\bm{\varepsilon})}{\mu_{\mathrm{tot}}}\bigl{|}n\bigr{\rangle}$
with the NMQSD equation up to time $t$ and then project on the same initial
state again. The average over many such trajectories will reproduce the
ensemble correlation function. Note that the form of Eq. (28) is very similar
to that of the noiseless equation.FTODTFD

It is easy to see (details can be found in appendix B.2) that Eq. 28 can be
written as

$C(t)={\mu}_{\mathrm{tot}}\,\mathcal{M}\Big{[}\mathrm{Tr}_{\mathrm{S}}\Big{\\{}\hat{\mu}_{\mathrm{eff}}\bigl{|}\psi_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}\bigl{\langle}\mathrm{g}(t)\bigr{|}\Big{\\}}\Big{]}$
(29)

with
$\bigl{|}\mathrm{g}(t)\bigr{\rangle}=e^{-i\hat{H}_{\mathrm{g}}t}\bigl{|}\mathrm{g}\bigr{\rangle}$,
which can also be written as

$C(t)=\mathcal{M}\Big{[}\mathrm{Tr}_{\mathrm{S}}\\{\hat{\mu}_{\mathrm{eff}}\,(\hat{U}_{\mathrm{ex}}(t,z^{*})\hat{\mu}_{\mathrm{eff}}\bigl{|}\mathrm{g}\bigr{\rangle})\bigl{\langle}\mathrm{g}\bigr{|}\hat{U}^{\dagger}_{\mathrm{g}}(t)\\}\Big{]}$
(30)

where $\hat{U}_{\mathrm{ex}}(t,z^{*})$ denotes the stochastic propagator in
the excited state manifold, i.e. evolution with the respective NMQSD equation
NMQSDfootnote  We see that this equation has now exactly the form of a dyadic
scheme as sketched in Fig. 1: The bra is during the whole time evolution
propagated in the electronic ground state. The ket is first lifted via
$\hat{\mu}_{\mathrm{eff}}$ to the excited state manifold where it is then
propagated during the time period $t$. Finally, $\hat{\mu}_{\mathrm{eff}}$ is
applied again and the trace is taken. One has to be careful about the formal
meaning of the stochastic state
$\bigl{|}\psi_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}$ in the case of the non-
linear NMQSD equation. We will discuss this in detail in the next section.

Figure 2:  Example calculations for the case of a bath-correlation function
described by one exponential with $\Omega=1.0$ and $\gamma=0.25\Omega$. Two
cases of the coupling strength are shown: intermediate coupling $p=0.5$ (left
box) and strong coupling $=2$ (right box). The dipole-dipole coupling between
the monomer is $V=\Omega$. In each box we show the following: (a) Monomer
spectrum. (b) real and imaginary parts of $C(t)$, (c) absorption spectrum
calculated by the linear HOPS equation. (d) real and imaginary parts of
$C(t)$, (e) absorption spectrum calculated by the non-linear HOPS equation.
The truncation depth of the hierarchy for HOPS are $\mathcal{K}=6$ ($p=0.5$)
and $\mathcal{K}=12$ ($p=2$). The noise-free HOPS calculations are performed
using the formulas provided in appendix A. For the stochastic calculations Eq.
(28) and Eq. (33) are used for the linear and non-linear HOPS, respectively.

###  III.3 Dyadic non-linear HOPS equation

While we have demonstrated that the NMQSD equations can be used to directly
propagate a dyadic equation for the dipole-dipole auto-correlation function,
but care is required to establish the corresponding non-linear HOPS equation.
In the non-linear form of the starting equation Eq. (28) one propagates the
pure state vectors $\bigl{|}\tilde{v}_{\eta}(t,z^{*})\bigr{\rangle}$ according
to Eq. (14) in the total Hilbert space spanned by the singly excited states
$\bigl{|}n\bigr{\rangle}$ and the ground state
$\bigl{|}\mathrm{g}\bigr{\rangle}$. It has the same form as the Eq. (25)
$\bigl{|}\tilde{v}_{\eta}(t,z^{*})\bigr{\rangle}=\frac{1}{\sqrt{2}}\big{(}\bigl{|}\tilde{\psi}_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}+\eta
e^{-i\epsilon_{\mathrm{g}}t}\bigl{|}\mathrm{g}\bigr{\rangle})$, where
$\bigl{|}\tilde{\psi}_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}$ is evolved
according to the non-linear equation Eq. (14). To perform expectation values
one uses the normalized states

$\bigl{|}\tilde{V}_{\eta}(t,z^{*})\bigr{\rangle}=\frac{\bigl{|}\tilde{v}_{\eta}(t,z^{*})\bigr{\rangle}}{||\tilde{v}_{\eta}(t,z^{*})||}$
(31)

with

$\begin{split}||\tilde{v}_{\eta}(t,z^{*})||^{2}&\equiv\langle\tilde{v}_{\eta}(t,z^{*})|\tilde{v}_{\eta}(t,z^{*})\rangle\\\
&=\frac{1}{2}\left[\langle\tilde{\psi}_{\mathrm{ex}}(t,z^{*})|\tilde{\psi}_{\mathrm{ex}}(t,z^{*})\rangle+1\right].\end{split}.$
(32)

Note, that the last expression is independent of $\eta$, which means that for
all four initial state one has the same normalization factor. Therefore,
repeating the steps that leads form Eq. (24) to Eq. (28) we now arrive at

$C(t)=\mu_{\mathrm{tot}}^{2}\mathcal{M}\Big{[}\frac{\langle\psi_{\mathrm{ex}}|\tilde{\psi}_{\mathrm{ex}}(t,z^{*})\rangle}{\frac{1}{2}\left(||\tilde{\psi}_{\mathrm{ex}}(t,z^{*})||^{2}+1\right)}\Big{]}e^{i\epsilon_{\mathrm{g}}t}.$
(33)

We emphasize, that here $\tilde{\psi}_{\mathrm{ex}}(t,z^{*})$ is propagated
with the nonlinear Eq. (14) in the excited state manifold only, but using
expectation values $\langle L_{n}^{\dagger}\rangle_{t}$ calculated with
respect to the normalized state
$\bigl{|}\tilde{\psi}_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}/\sqrt{({||\tilde{\psi}_{\mathrm{ex}}(t,z^{*})||^{2}+1})}$,
which contains a ground state contribution in the normalization.

##  IV Numerical calculations

In this section, we investigate the numerical performance of our dyadic
equation (Eq. (28)) using both the linear (Eq. (13)) and non-linear (Eq. (14))
HOPS equations and compare their convergence with respect to the number of
trajectories. In the following, we will calculate the linear absorption
spectra for a homodimer ($\epsilon_{n}=\epsilon$) with electronic coupling
$V=1$, where the transition dipoles of both chromophores are parallel. We use
a simple bath correlation function
$\alpha_{n}(\tau)=\alpha(\tau)=p\,e^{-i\Omega\tau-\gamma|\tau|}$ where
$\gamma=0.25$. Here, and in the following, we choose $\Omega$ as the unit of
energy and express $p$ in units of $(\Omega)^{2}$. The truncation constant for
the HOPS, $\mathcal{K}$, defined after Eq. (13), is provided in the respective
plots and is always chosen large enough to be well converged in the hierarchy.

Figure 3: Analysis of the error of the stochastic equations compared to the
noiseless reference. We consider the same parameters as in Fig. 2. Top: The
distribution of errors for four different numbers of trajectories, for both
the linear and the non-linear HOPS. Bottom: The mean error of
$\mathcal{A}(\omega)$ as a function of $\mathrm{N_{traj}}$. The straight lines
are guides for the eye and indicate a scaling of the mean error as
$r_{0}\sqrt{N_{\mathrm{traj}}}$

In Fig. 2 we demonstrate that both linear and non-linear HOPS are reasonably
converged using a rather small number of trajectories for both $p=0.5$ (left
box) and $p=2$ (right box). For both cases the spectrum of the monomer is
shown in the top row. In the second row the numerically calculated correlation
function of the dimer are shown, for both the linear and non-linear HOPS along
with the reference spectrum calculated using the noise-free HOPS algorithm.
The bottom row shows the corresponding absorption spectra.

We quantify the convergence with respect to the number of trajectories using
non-parametric error estimation and find the stochastic HOPS calculations
converge roughly with $1/\sqrt{N_{\mathrm{traj}}}$, with the non-linear HOPS
showing overall faster convergence. We quantify the difference between the
stochastic ($\mathcal{A}(\omega)$) and reference
($\mathcal{A}_{\mathrm{ref}}(\omega)$) spectrum using

$\mathrm{error}=\frac{1}{\omega_{\mathrm{max}}-\omega_{\mathrm{min}}}\int_{\omega_{\mathrm{min}}}^{\omega_{\mathrm{max}}}\big{|}\mathcal{A}(\omega)-\mathcal{A}_{\mathrm{ref}}(\omega)\big{|}\mathrm{d}\omega$
(34)

where the integration extend over the region of $\omega$ in which one is
interested in (we take the $\omega$ range shown in Fig. 2). Note that
according to this definition the error is proportional to
$\mu_{\mathrm{tot}}^{2}$, which in the case shown is equal to 2. In Fig. 3
(upper panels) we compare the distribution errors calculated from average
spectra for different ensembles with fixed $N_{\mathrm{traj}}=100$, $1000$,
$5000$, and $10000$ (see appendix C for details). In all cases the non-linear
HOPS calculation has error that is approximately half that of the linear HOPS
and both have a slightly asymmetric distribution with a full width at half
maximum roughly half of their mean. In the bottom row of Fig. 3 we show the
mean values as a function of $N_{\mathrm{traj}}$ on a log-log scale, where the
solid line provides a guide for the expected scaling of a mean error
proportional to $1/\sqrt{N_{\mathrm{traj}}}$.

Figure 4: Population dynamics of site 1 for $p=0.5$ (a) and $p=2$ (b). The
other parameters of the model are the same as in Fig. 2 and Fig. 3. The number
of trajectories used in both linear and non-linear HOPS is 10000. One sees
clearly that for population dynamics the linear HOPS does not converge for the
parameters considered. The truncation depth of the hierarchy for HOPS and HEOM
are $\mathcal{K}=6$ and $7$ for (a), $\mathcal{K}=16$ and $18$ for (b),
respectively.

These results demonstrate that both the linear and non-linear HOPS can produce
accurate correlation function $C(t)$ and absorption spectra
$\mathcal{A}(\omega)$. We note that for the strong coupling regime presented
here, the linear HOPS converges very slowly for population transfer. Fig. 4
compares the population dynamics of site 1 calculated by linear (red line) and
non-linear (blue line) HOPS for the two values $p=0.5$ and $p=2$ with the
exact HEOM calculations (black line). While non-linear HOPS with 10000
trajectories produces very accurate population dynamics, the linear HOPS
starts to deviate from the exact results at short time and completely fails to
describe the population dynamics at longer time. This inability of linear HOPS
to describe the strong coupling regime is well known for population dyanmics.
That linear HOPS can reproduce absorption spectra is a consequence of the
cancellation of noise that can be clearly understood in the context of the
pure state decomposition.

##  V Conclusions

We derived stochastic propagation schemes for the calculation of absorption
spectra within the NMQSD approach, using both the linear and the non-linear
NMQSD equation. Beside a scheme that relies on formally propagating pure
states in the full Hilbert-space, we derived a scheme, where individual
propagations in the electroinc ground and electronic excited state are used to
obtain the desired spectrum. This scheme directly resembles the diagrammatic
perturbation theory where dyadic matrices are propagated. While for the linear
NMQSD equation the two propagations are completely independent, for the
nonlinear NMQSD they become coupled via a common normalization factor that
contains the norm of the two wavefunctions in the different Hilbert-spaces.
Besides its favorable convergence properties with respect to the number of
trajectories, this non-linear normalized version will allow a direct
implementation of adaptive algorithms. We investigated the convergence with
respect to trajectories $N_{\mathrm{traj}}$in detail and found that the error
decreases as $1/\sqrt{N_{\mathrm{traj}}}$ with the linear version needing
roughly twice the number of trajectories as the non-linear one to achieve the
same accuracy. This behaviour of linear equation is remarkable, since we
worked in a parameter regime where population dynamics does not converge for
the linear NMQSD and our decomposition into pure states does not lead to a
cancellation as was found in Ref. HOPSEDRichard, .

We have formulated the absorption spectrum in terms of propagating a single
initial state in the electronic excited state. This state contains all the
information of the dipoles (magnitude and orientation) of the individual
molecules as well as the local fields at these molecules (for ease of notation
we used the same electric field for all molecules). Since arbitrary
distributions of the electromagnetic field can be treated, the formalism can
be directly applied to the calculation of circular dichroism or near-field
spectra. We note that similar to the treatment of Ref. FTODTFD, , one can also
construct the spectrum from an summation of initially localized states, with
the correct wheighting factors. Such a scheme will be beneficial for an
adaptive treatment.

We would like to remark, that our derivations relies on the assumption that
decay processes from excited to ground state are negligible (as we have
already used in the starting Hamiltonian). This is the case in many relevant
situations, where these processes are on the nanosecond timescale, compared to
a few femtoseconds that are needed to find well resolved spectra. While such
coupling between ground and excited electronic states can be directly treated
using the decomposition into pure states and propagating via NMQSD in the full
electronic Hilbert space, presumably the dyadic equation will need major
modifications.

We demonstrated the applicability of our schemes explicitly using the HOPS
formulation of NMQSD. It should be noted that our equation are formulated for
arbitrary temperatures. Within HOPS there are different ways how to
incooperate temperature either in the hierarchy FTODTFD  or in additional
stochastic processes.HOPS2  Because our schemes propagate single trajectories
subject to noise, it is straightforward to account for the effect of the
static disorder induced by the inhomogeneity of the environment without much
additional computational cost. Since the HOPS provides efficient treatment of
the environmental degrees of freedom and furthermore one propagates vectors
instead of matrices the propagation schemes offer promising techniques to
simulate absorption spectra of large molecular aggregates with complicated
structured environments. Our dyadic HOPS equation is also an important step
towards application of HOPS for non-linear spectroscopy.

###### Acknowledgements.

LPC acknowledges support from the Max-Planck Gesellschaft via the MPI-PKS
visitors program. AE acknowledges support from the DFG via a Heisenberg
fellowship (Grant No EI 872/5-1). DIGB acknowledges support from Robert A.
Welch Foundation (Grant N-2026-20200401).

## Data availability

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

##  Appendix A Linear absortion using HOPS without noise

It is demonstrated in Ref. FTODTFD,  that for the calculation of absorption
spectra, one can propagate the linear HOPS equation (Eq. 13) without noise,
i.e. with all $z_{t,n}=0$. The dipole-correlation function is calculated
accoring to

$C(t)=\mu_{\mathrm{tot}}^{2}\langle\psi_{\mathrm{ex}}|\psi^{(\mathbf{0})}_{\mathrm{ex}}(t)\rangle$
(35)

with

$\begin{split}\partial_{t}|\psi_{\mathrm{ex}}^{(\mathbf{k})}(t)\rangle=\left(-i\hat{H}_{\mathrm{S}}-\mathbf{k}\cdot\mathbf{w}\right)|\psi_{\mathrm{ex}}^{(\mathbf{k})}(t)\rangle\\\
+\sum_{n}\hat{L}_{n}\sum_{j}k_{nj}p_{nj}|\psi_{\mathrm{ex}}^{(\mathbf{k}-\mathbf{e}_{nj})}(t)\rangle\\\
-\sum_{n}\hat{L}_{n}^{\dagger}\sum_{j}|\psi_{\mathrm{ex}}^{(\mathbf{k}+\mathbf{e}_{nj})}(t)\rangle\end{split}$
(36)

and the initial condition
$|\psi_{\mathrm{ex}}^{(\mathbf{0})}(t\\!=\\!0)\rangle=\bigl{|}\psi_{\mathrm{ex}}\bigr{\rangle}$.

##  Appendix B Calculations for the Hartmann-Strunz approach

###  B.1 Writing the initial operator as a sum of pure states

To obtain equation (20) we first rewrite
$\hat{\mu}_{\mathrm{eff}}\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}$
as the sum of two hermitian matrices,

$\hat{\mu}_{\mathrm{eff}}\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}=\frac{1}{2}\big{(}\hat{A}+\hat{B}\big{)}$
(37)

with

$\displaystyle\hat{A}=$
$\displaystyle\big{\\{}\hat{\mu}_{\mathrm{eff}}\,,\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}\big{\\}}=$
$\displaystyle\sum_{n=1}^{N}(\bm{\mu}_{n}\cdot\bm{\varepsilon})\Big{(}\bigl{|}n\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}+\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}n\bigr{|}\Big{)}$
(38) $\displaystyle\hat{B}=$
$\displaystyle\big{[}\hat{\mu}_{\mathrm{eff}}\,,\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}\big{]}=$
$\displaystyle\sum_{n=1}^{N}(\bm{\mu}_{n}\cdot\bm{\varepsilon})\Big{(}\bigl{|}n\bigr{\rangle}\bigl{\langle}\mathrm{g}\bigr{|}-\bigl{|}\mathrm{g}\bigr{\rangle}\bigl{\langle}n\bigr{|}\Big{)}$
(39)

The eigenvalues and eigenvectors of these two operators can be calculated
analytically

$\displaystyle\hat{A}\bigl{|}v_{\pm 1}\bigr{\rangle}$ $\displaystyle=$
$\displaystyle\pm\mu_{\mathrm{tot}}\bigl{|}v_{\pm 1}\bigr{\rangle}$ (40)
$\displaystyle\hat{B}\bigl{|}v_{\pm i}\bigr{\rangle}$ $\displaystyle=$
$\displaystyle\pm{i}\mu_{\mathrm{tot}}\bigl{|}v_{\pm i}\bigr{\rangle}$ (41)

where the eigenvectores are that defined in Eq. (21) of the main text.

###  B.2 Derivation of the dyadic equation

Here we show, that the dyadic equation (29) is identical to Eq. (28). We start
with Eq. (29)

$\begin{split}{\mu}_{\mathrm{tot}}&\,\mathcal{M}\Big{[}\mathrm{Tr}_{\mathrm{S}}\Big{\\{}\hat{\mu}_{\mathrm{eff}}\bigl{|}\psi_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}\bigl{\langle}\mathrm{g}(t)\bigr{|}\Big{\\}}\Big{]}\\\
&={\mu}_{\mathrm{tot}}\,\mathcal{M}\Big{[}\mathrm{Tr}_{\mathrm{S}}\Big{\\{}\hat{\mu}_{\mathrm{eff}}\bigl{|}\psi_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}e^{i\epsilon_{\mathrm{g}}t}\bigl{\langle}\mathrm{g}\bigr{|}\Big{\\}}\Big{]}\\\
&={\mu}_{\mathrm{tot}}\,\mathcal{M}\Big{[}\bigl{\langle}\mathrm{g}\bigr{|}\hat{\mu}_{\mathrm{eff}}\bigl{|}\psi_{\mathrm{ex}}(t,z^{*})\bigr{\rangle}e^{i\epsilon_{\mathrm{g}}t}\Big{]}\\\
&={\mu}_{\mathrm{tot}}^{2}\,\mathcal{M}\Big{[}\langle\psi_{\mathrm{ex}}|\psi_{\mathrm{ex}}(t,z^{*})\rangle
e^{i\epsilon_{\mathrm{g}}t}\Big{]}\end{split}$ (42)

where in the first step we used
$\bigl{\langle}\mathrm{g}(t)\bigr{|}=e^{i\epsilon_{\mathrm{g}}t}\bigl{\langle}\mathrm{g}\bigr{|}$,
in the second step we evaluated the trace and in the last step we have used
Eq. (22). We see that the last line is identical to Eq. (28).

##  Appendix C Analysis of the statistic error by the bootstrapping method

To get a quantitative description of the convergence properties of linear and
non-linear HOPS with respect to the number of trajectories, we conduct a
detailed analysis of the statistic error due to a finite number of
trajectories by using the bootstrapping technique BootStrap . To this end, we
first calculate $3\times{10}^{4}$ trajectories. For each value of
$\mathrm{N_{traj}}$, we construct $10^{4}$ ensembles by selecting
$\mathrm{N_{traj}}$ trajectories randomly from the original $3\times{10}^{4}$
trajectories (the same trajectory can appear multiple times within each
ensemble). For each ensemble, we quantify the error by evaluating the average
absolute difference of $\mathcal{A}(\omega)$.

## References

  * (1) H. van Amerongen, R. van Grondelle, and L. Valkunas, Photosynthetic Excitons. (World Scientific, 2000).
  * (2) V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 3rd ed. (Wiley-VCH, 2011).
  * (3) D. Abramavicius, B. Palmieri, D. V. Voronine, F. Šanda, and S. Mukamel, Chem. Rev. 109, 2350-2408 (2009).
  * (4) L. Chen, P. Shenai, F. Zheng, A. Somoza, and Y. Zhao, Molecules. 20, 15224-15272 (2015).
  * (5) C. Olbrich, and U. Kleinekathöfer, J. Phys. Chem. B. 114, 12427-12437 (2010).
  * (6) C. Olbrich, J. Strümpfer, K. Schulten, and U. Kleinekathöfer, J. Phys. Chem. Lett. 2, 1771-1776 (2011).
  * (7) S. Valleau, A. Eisfeld, and A. Aspuru-Guzik, J. Chem. Phys. 137, 224103 (2012).
  * (8) S. Maity, V. Daskalakis, M. Elstner, and U. Kleinekathöfer, Phys. Chem. Chem. Phys. 23, 7407 (2021).
  * (9) S. Jang, and B. Mennucci, Rev. Mod. Phys. 90, 035003 (2018).
  * (10) U. Weiss, Quantum Dissipative Systems, Series in Modern Condensed Matter Systems. (World Scientific, 2008).
  * (11) H. P. Breuer, and F. Petruccione, The Theory of Open Quantum Systems. (Oxford University Press, 2007).
  * (12) Y. Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006).
  * (13) Y. Tanimura, J. Chem. Phys. 153, 020901 (2020).
  * (14) N. Makri, and D. E. Makarov, J. Chem. Phys. 102, 4600-4610 (1995).
  * (15) N. Makri, and D. E. Makarov, J. Chem. Phys. 102, 4611-4618 (1995).
  * (16) D. P. S. McCutcheon, and A. Nazir, J. Chem. Phys. 135, 114501 (2011).
  * (17) A. Kimura, and Y. Fujihashi, J. Chem. Phys. 141, 194110 (2014).
  * (18) M. H. Beck, A. Jäckle, G. A. Worth, and H. D. Meyer, Phys. Rep. 324, 1-105 (2000).
  * (19) H. Wang, and M. Thoss, J. Chem. Phys. 119, 1289-1299 (2003).
  * (20) U. Manthe, J. Chem. Phys. 128, 164116 (2008).
  * (21) I. Burghardt, M. Nest, and G. A. Worth, J. Chem. Phys. 119, 5364 (2003).
  * (22) G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt, and B. Lasorne, Int. Rev. Phys. Chem. 34, 269-308 (2015).
  * (23) L. Chen, R. Borrelli, and Y. Zhao, J. Phys. Chem. A. 121, 8757-8770 (2017).
  * (24) L. Chen, M. F. Gelin, and Y. Zhao, Chem. Phys. 515, 108-118 (2018).
  * (25) D. V. Shalashilin, J. Chem. Phys. 130, 244101 (2009).
  * (26) L. Chen, M. F. Gelin, and D. V. Shalashilin, J. Chem. Phys. 151, 244116 (2019).
  * (27) L. Diósi, and W. T. Strunz, Phys. Lett. A. 235, 569-573 (1997).
  * (28) L. Diósi, N. Gisin, and W. T. Strunz, Phys. Rev. A. 58, 1699 (1998).
  * (29) T. Yu, L. Diósi, N. Gisin, and W. T. Strunz, Phys. Rev. A. 60, 91 (1999).
  * (30) D. Suess, A. Eisfeld, and W. T. Strunz, Phys. Rev. Lett. 113, 150403 (2014).
  * (31) R. Hartmann, and W. T. Strunz, J. Chem. Theory. Comput. 13, 5834-5845 (2017).
  * (32) X. Gao, and A. Eisfeld, J. Chem. Phys. 150, 234115 (2019).
  * (33) L. Varvelo, J. K. Lynd, and D. I. G. Bennett, Chem. Sci. 12, 9704 (2021).
  * (34) S. Mukamel, Principles of Nonlinear Optical Spectroscopy. (Oxford University Press, 1995).
  * (35) P. Zhang, C. D. B. Bentley, and A. Eisfeld, J. Chem. Phys. 148, 134103 (2018).
  * (36) G. Ritschel, D. Suess, S. Möbius, W. T. Strunz, and A. Eisfeld, J. Chem. Phys. 142, 034115 (2015).
  * (37) X. Gao, and A. Eisfeld, J. Phys. Chem. Lett. 9, 6003-6010 (2018).
  * (38) F. Zheng, X. Gao, and A. Eisfeld, Phys. Rev. Lett. 123, 163202 (2019).
  * (39) R. Hartmann, and W. T. Strunz, J. Phys. Chem. A. 125, 7066-7079 (2021).
  * (40) This equation is equivalent to $C(t)=\mathcal{M}\Big{[}\mathrm{Tr}_{\mathrm{S}}\\{\hat{\mu}_{\mathrm{eff}}\,(\hat{U}(t,\mathbf{z}^{*})\hat{\mu}_{\mathrm{eff}}\bigl{|}\mathrm{g}\bigr{\rangle})\bigl{\langle}\mathrm{g}\bigr{|}\hat{U}^{\dagger}(t,\mathbf{z}^{*})\\}\Big{]}$ where $\hat{U}(t,\mathbf{z}^{*})$ denotes the stochastic propagator in the complete Hilbert space.
  * (41) A. C. Davison and D. V. Hinkley, Bootstrap Methods and Their Application, Cambridge Series in Statistical and Probabilistic Mathematics. (Cambridge University Press, 1997).