# The hierarchical and perturbative forms of stochastic Schrödinger equations
and their applications to carrier dynamics in organic materials

Yuchen Wang, Yaling Ke11footnotemark: 1, Yi Zhao11footnotemark: 1
,22footnotemark: 2 State Key Laboratory for Physical Chemistry of Solid
Surfaces, Collaborative Innovation Center of Chemistry for Energy Materials,
Fujian Provincial Key Lab of Theoretical and Computational Chemistry, and
College of Chemistry and Chemical Engineering, Xiamen University, Xiamen
361005, ChinaE-mail<EMAIL_ADDRESS>
#### Article Type:

Advanced Review

#### Abstract

A number of non-Markovian stochastic Schrödinger equations, ranging from the
numerically exact hierarchical form towards a series of perturbative
expressions sequentially presented in an ascending degrees of approximations
are revisited in this short review, aiming at providing a systematic framework
which is capable to connect different kinds of the wavefunction-based
approaches for an open system coupled to the harmonic bath. One can
optimistically expect the extensive future applications of those non-Markovian
stochastic Schrödinger equations in large-scale realistic complex systems,
benefiting from their favorable scaling with respect to the system size, the
stochastic nature which is extremely suitable for parallel computing, and many
other distinctive advantages. In addition, we have presented a few examples
showing the excitation energy transfer in Fenna-Matthews-Olson complex, a
quantitative measure of decoherence timescale of hot exciton, and the study of
quantum interference effects upon the singlet fission processes in organic
materials, since a deep understanding of both mechanisms is very important to
explore the underlying microscopic processes and to provide novel design
principles for highly efficient organic photovoltaics.

## INTRODUCTION

With the burst of newly-fabricated materials and rapid developments of
experimental techniques, it has been demonstrated by numerous studies that the
quantum effects play a fundamental role in ultrafast dynamics, such as the
excitation energy transfer in photosynthetic antenna complex[1, 2, 3, 4, 5, 6,
7, 8] and the carrier dynamics in organic photovoltaics[9, 10, 11, 12, 13, 14,
15].

To elucidate the detailed mechanism behind these processes, it is often
necessary to invoke accurate theoretical descriptions, which is challenging
due to the large (often nearly innumerous) degrees of freedom in realistic
complex systems. Despite the difficulty, remarkable progress has been made in
the past decades. For the purpose of reducing the number of basis vectors with
controllable numerical errors, different brilliant approaches have been
developed, the representative examples are the numerical renormalization
group[16, 17], density matrix renormalization group[18], time evolving density
matrix using orthogonal polynomials algorithms[19] and multi-layer multi-
configuration time-dependent Hartree method (ML-MCTDH)[20]. These approaches
have been proven to be very efficient in the low temperature regime and have
achieved great success in many areas[21, 22, 23, 24, 25, 26].

Another formalism different from the aforementioned full-space approaches is
based on the reduced system language. In organic materials, for instance, the
carriers (usually regarded as the system) can be described by the
corresponding reduced density matrix, whereas the vibrational degrees of
freedom are mapped into a bosonic bath and are traced over. This partial trace
operation can be done in an elegant way using the path integral technique,
leading to the appearance of the famous Feynman-Vernon influence functional in
the path integral expression[27], in which all the influences from the bath
are included. But this operation also introduces new difficulties since the
influence functional is nonlocal both in space and time. The long time results
of the path integral, however, cannot be obtained by directly using the Monte
Carlo technique[28, 29, 30] due to the notorious sign problem. By utilizing
the well-behaved quasiadiabatic propagators, Makri and coworkers[31, 32, 33]
creatively proposed the quasiadiabatic propagator path integral to overcome
this problem. Whereafter, various numerical optimization techniques, like the
discrete value representation[34], tensor multiplication scheme[35, 36, 37,
38, 39, 40] and propagator filtering techniques[41] were proposed successively
to further lower the cost and extend the applicability of the method. Another
efficient and accurate approach is the hierarchical equations of motion,
pioneered by Tanimura and Kubo[42] for the special case of Drude-Lorentzian
spectral densities with the high temperature approximation. Afterwards, the
exact version of hierarchical equations of motion was derived by several
contributors via different starting point[43, 44, 45]. To date, a lot of
efforts have been made to equip hierarchical equations of motion with
different powerful techniques, which can be listed as follows: the Padé
spectrum decomposition[46], better truncation schemes[44, 47, 48], new
decomposition schemes for arbitrary spectral densities[49] and bath
correlation functions[50, 51], highly efficient filtering algorithm
controllable in accuracy[52, 53], exponential integrators for time
evolution[54], high-performance computing platform and efficient parallel
algorithm[55, 56, 57, 58, 59, 60]. Nowadays, hierarchical equations of motion
has become a standard approach in the fields of open quantum systems, and has
made great contributions to many areas of physics and chemistry[61, 62, 63,
47]. A different strategy from quasiadiabatic propagator path integral and
hierarchical equations of motion to treat the influence functional is the
direct stochastic unravelling by introducing two correlated Gaussian complex
stochastic processes, leading to the stochastic Liouville-von Neumann
equation[64, 65, 66, 67, 68], which is a stochastic version of the Liouville-
von Neumann equation without explicit environmental memory effect term in the
equation. This equation can be applied to arbitrary spectral densities, but it
may suffer from severe numerical instability in the long time simulation when
the system-bath interaction is not small. To overcome this problem, the idea
of partial stochastic unraveling combined with partial hierarchical expansion
was proposed, prompting the appearances of mixed deterministic-stochastic
approaches[47, 69, 70, 71, 72], which inherit both merits of the stochastic
Liouville-von Neumann equation and hierarchical equations of motion.

Despite the great success, the aforementioned numerically exact methods are
still limited in small sized systems due to the huge computational cost.
Towards the simulations of large-scale systems, one may turn to suitable
approximate methods. Along with the history of quantum dynamics in open
systems, quantum master equations play a prominent role[73]. Starting from a
formally exact quantum master equation[74, 75, 76], systematic approximations
can be made to obtain computable equations of motion for different parameter
regimes. For example, in the weak exciton-phonon coupling regime, a
perturbative treatment on the exciton-phonon coupling leads to the second-
order time-convoluted or time-convolutionless quantum master equation, the
Markovian limit of which is the famous Bloch-Redfield equation[77, 78]. Beyond
the weak exciton-phonon coupling regime, the polaron transformation
technique[79, 80, 81, 82, 83, 84] and its variational version[85, 86, 87, 88,
89, 90] are usually employed for a more accurate quantum master equation.

As compared with density matrix approaches, the non-Markovian stochastic
Schrödinger equation (NMSSE) has natural superiority in numerical calculations
due to the beneficial linear scaling of Hilbert space. About a few decades
ago, Strunz, Diósi and coworkers[91, 92, 93, 94] pioneered the foundation of
NMSSEs and proposed a formally exact NMSSE, the non-Markovian quantum state
diffusion. However, the explicit expression of the functional derivative term
appearing in the equation is usually unknown, and one often invokes the
version with zeroth order functional expansion approximation[95, 96]. Soon
after the work of Strunz et al., Gaspard and Nagoka[97, 98, 99] proposed a
different NMSSE from the Feshbach projection method, which is only applicable
to the situations when the system-bath interaction is weak and the bath
relaxation is fast. Lately another kind of NMSSE named the time-dependent
wavepacket diffusion method was also proposed by Zhong and Zhao[100], and its
high numerical performance was extensively demonstrated[100, 101, 102, 13,
103, 104, 105, 106, 107]. Very recently, inspired by the spirits of mixed
deterministic-stochastic approaches, several studies showed that there exists
a corresponding numerically exact hierarchical counterpart for a specific
NMSSE. Suess et al.[108, 109, 110] applied the hierarchical expansion
technique to the non-Markovian quantum state diffusion method and developed
the hierarchy of pure states, which is proved to be highly efficient. Song et
al.[111] utilized the same idea and obtained the hierarchical form of time-
dependent wavepacket diffusion method. Ke and Zhao developed a new
hierarchical form of stochastic Schrödinger equation[112], and successively
got its NMSSE version[113]. Various applications have shown the power of these
NMSSEs and their hierarchical form in the simulations of complex systems
dynamics.

In this paper, we are aimed at reviewing recent works about NMSSEs, their
hierarchical forms and the applications. In the next section, we start from a
general partial stochastic unravelling scheme in the path-integral formalism,
disentangling the forward and backward paths in the influence functional,
which paves the way to the wavefunction-based framwork. Then the hierarchical
technique and systematic perturbation expansion are applied to the rest part
of the influence functional, the former leads to the hierarchical form of
stochastic Schrödinger equations, whereas the latter leads to a set of NMSSEs.
In the third section, recent applications are presented to show the power of
these approaches. Concluding remarks are given in the last section.

## Methodology

### Background

We start with a generic model where the total Hamiltonian consists of three
parts:

$\hat{H}_{tot}=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{I}.$ (1)

$\hat{H}_{S}$, $\hat{H}_{B}$, and $\hat{H}_{I}$ denote the system, the
harmonic bath, and their interaction, respectively. Although the boundary
between the system and its environment can be quite flexible, in molecular
aggregates, like organic materials and photosynthetic systems, it is typical
to choose the electronic degrees of freedom as the system part. The explicit
form of $\hat{H}_{B}$ and $\hat{H}_{I}$ can be written as

$\hat{H}_{B}=\sum_{k}\frac{1}{2}(\hat{P}_{k}^{2}+\omega_{k}^{2}\hat{Q}_{k}^{2}),$
(2)

and

$\hat{H}_{I}=\hat{x}\otimes\sum_{k}c_{k}\hat{Q}_{k},$ (3)

where $\hat{P}_{k}$, $\hat{Q}_{k}$ and $\omega_{k}$ are the mass-weighted
momentum, coordinate operator, and the frequency of the $k$-th bath mode,
respectively, $\hat{x}$ is an operator characterizing the system-bath coupling
manner and $c_{k}$ the coupling strength. In many cases, only the system
dynamics is of interest, such that we would forward to obtain the reduced
density operator by tracing over all the bath degrees of freedom. Assuming
that the initial total density operator is factorized between the system and a
thermal-equilibrium bath, i.e.,
$\hat{\rho}_{tot}(0)=\hat{\rho}_{S}(0)\otimes\hat{\rho}_{B}(0)$ with
$\hat{\rho}_{B}(0)=e^{-\beta\hat{H}_{B}}/{\rm Tr}\\{e^{-\beta\hat{H}_{B}}\\}$,
where $\beta$ is the inverse temperature, the time-evolving reduced density
operator of the system is given in the path-integral formalism as ($\hbar$ and
$k_{B}$ are set to be unity for simplicity)

$\begin{split}\rho_{S}(x_{f},x_{i},t)=&\langle x_{f}|{\rm
Tr}_{B}\\{e^{-i\hat{H}_{tot}t}\hat{\rho}_{tot}(0)e^{i\hat{H}_{tot}t}\\}|x_{i}\rangle\\\
=&\int\mathrm{d}x_{0}^{+}\int\mathrm{d}x_{0}^{-}\int_{x_{0}^{+}}^{x_{f}}\mathcal{D}x^{+}(t)\int_{x_{0}^{-}}^{x_{i}}\mathcal{D}x^{-}(t)e^{iS_{0}[x^{+}(t)]-iS_{0}[x^{-}(t)]}\rho_{S}(x_{0}^{+},x_{0}^{-},0)\mathcal{F}[x^{+}(t),x^{-}(t),t],\end{split}$
(4)

where $\rho_{S}(x_{0}^{+},x_{0}^{-},0)=\langle
x_{0}^{+}|\hat{\rho}_{S}(0)|x_{0}^{-}\rangle$, $S_{0}[x(t)]$ is the action
functional corresponding to $\hat{H}_{S}$, and
$\mathcal{F}[x^{+}(t),x^{-}(t),t]$ is the Feynman-Vernon influence
functional[27] incorporating all the dynamical influences from the bath with
the explicit form being

$\mathcal{F}[x^{+}(t),x^{-}(t),t]=\exp\left\\{-\int_{0}^{t}\mathrm{d}\tau\int_{0}^{\tau}\mathrm{d}\tau^{\prime}\Phi(\tau,\tau^{\prime})\right\\},$
(5)

and

$\Phi(\tau,\tau^{\prime})=[x^{+}(\tau)-x^{-}(\tau)][\alpha(\tau-\tau^{\prime})x^{+}(\tau^{\prime})-\alpha^{\ast}(\tau-\tau^{\prime})x^{-}(\tau^{\prime})].$
(6)

Here $\alpha(t)$ is the bath correlation function, and its real and imaginary
parts are denoted as $\alpha_{R}(t)$ and $\alpha_{I}(t)$, respectively.
Defining the spectral density

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

which totally characterizes the system-bath interaction, $\alpha(t)$ can be
expressed as

$\begin{split}\alpha(t)=&\alpha_{R}(t)+i\alpha_{I}(t)\\\
=&\int_{0}^{\infty}\mathrm{d}\omega\frac{J(\omega)}{\pi}[\coth{\frac{\beta\omega}{2}}\cos{\omega
t}-i\sin{\omega t}].\end{split}$ (8)

The complexity in the calculation of $\hat{\rho}_{S}(t)$ originates from the
time-nonlocal nature of the influence functional as well as the entanglement
between forward and backward paths in Eq. (5). Two powerful schemes, the
stochastic unravelling[67, 68, 64, 65, 66] and deterministic hierarchical
technique[42, 43, 44, 45], have been proposed to deal with these problems
without any approximation. The total stochastic unravelling of the influence
functional is straightforward, but its application to the realistic systems is
limited by the convergence performance of stochastic average. Thus a variety
of mixed deterministic-deterministic approaches have been proposed utilizing
the partial stochastic unravelling strategy to improve the convergence
speed[47, 69, 70, 71, 72]. In this review, we will follows this stream but
with an additional requirement of the approaches capable of disentangling the
forward and backward paths, leading to the wavefunction-based framework.

### Partial Stochastic Unravelling

The partial stochastic unravelling requires an artificial partition of the
influence functional, and the key to optimizing the convergence behavior of
the stochastic ensemble average lies in this partition strategy. We first
divide the bath correlation function $\alpha(t)$ into two parts

$\alpha(t)=\alpha^{1}(t)+\alpha_{res}(t),$ (9)

with

$\alpha^{1}(t)=\int_{0}^{\infty}\mathrm{d}\omega\frac{J(\omega)}{\pi}[2\bar{n}(\omega)+g(\omega)]\cos{\omega
t},$ (10a)
$\alpha_{res}(t)=\int_{0}^{\infty}\mathrm{d}\omega\frac{J(\omega)}{\pi}\left\\{\left[1-g(\omega)\right]\cos{\omega
t}-i\sin\omega t\right\\},$ (10b)

where $\bar{n}(\omega)=\frac{1}{e^{\beta\omega}-1}$ is the Bose distribution
function, and $g(\omega)$ is a frequency-dependent function that will be
specified later. Based on Eq. (10), we further rewrite
$\Phi(\tau,\tau^{\prime})$ as the sum of a primary part and a residual part,

$\Phi(\tau,\tau^{\prime})=\Phi_{p}(\tau,\tau^{\prime})+\Phi_{r}(\tau,\tau^{\prime}),$
(11)

with

$\begin{split}\Phi_{p}(\tau,\tau^{\prime})=&x^{+}(\tau)\alpha_{R}^{1}(\tau-\tau^{\prime})x^{+}(\tau^{\prime})+x^{-}(\tau)\alpha_{R}^{1}(\tau-\tau^{\prime})x^{-}(\tau^{\prime})\\\
&-x^{+}(\tau)\alpha^{\ast}(\tau-\tau^{\prime})x^{-}(\tau^{\prime})-x^{-}(\tau)\alpha(\tau-\tau^{\prime})x^{+}(\tau^{\prime})\end{split}$
(12)

and

$\Phi_{r}(\tau,\tau^{\prime})=x^{+}(\tau)\alpha_{res}(\tau-\tau^{\prime})x^{+}(\tau^{\prime})+x^{-}(\tau)\alpha_{res}^{\ast}(\tau-\tau^{\prime})x^{-}(\tau^{\prime}).$
(13)

Note that the cross terms between $x^{+}(t)$ and $x^{-}(t)$ are entirely
included in $\Phi_{p}(\tau,\tau^{\prime})$. Thus if the partial stochastic
unravelling is introduced for the corresponding part of the influence
functional, the explicit entanglement between forward and backward paths can
be resolved. One way towards the realization of stochastic unravelling is to
use the Hubbard-Stratonovich transformation[114, 115]. This results in an
equivalent expression of Eq. (5):

$\mathcal{F}[x^{+}(t),x^{-}(t),t]=\int\mathcal{D}^{2}\chi_{1}\int\mathcal{D}^{2}\chi_{2}^{\ast}P[\chi_{1}(t),\chi_{2}^{\ast}(t),t]\exp^{-i\int_{0}^{t}[\chi_{1}(\tau)x^{+}(\tau)-\chi_{2}^{\ast}(\tau)x^{-}(\tau)]\mathrm{d}\tau-\int_{0}^{t}\mathrm{d}\tau\int_{0}^{\tau}\mathrm{d}\tau^{\prime}\Phi_{r}(\tau,\tau^{\prime})}.$
(14)

Here, two correlated complex stochastic processes $\chi_{1}(t)$ and
$\chi_{2}^{\ast}(t)$ are introduced, and $P[\chi_{1}(t),\chi_{2}^{\ast}(t),t]$
is the corresponding Gaussian probability density functional. By averaging
over $\chi_{1}(t)$ and $\chi_{2}^{\ast}(t)$, Eq. (14) must recover to Eq. (5).
To this end, it is found that $\chi_{1}(t)$ and $\chi_{2}^{\ast}(t)$ should
satisfy the following statistical properties

$\left\\{\begin{array}[]{l}\left\langle\chi_{1}(t)\right\rangle=\left\langle\chi_{2}^{\ast}(t)\right\rangle=0,\\\
\left\langle\chi_{1}(t)\chi_{1}(t^{\prime})\right\rangle=\left\langle\chi_{2}^{\ast}(t)\chi_{2}^{\ast}(t^{\prime})\right\rangle=\alpha_{R}^{1}(t-t^{\prime}),\\\
\left\langle\chi_{1}(t)\chi_{2}^{\ast}(t^{\prime})\right\rangle=\alpha^{\ast}(t-t^{\prime}).\end{array}\right.$
(15)

In the coming subsections, we will provide the explicit generation scheme of
$\chi_{1}(t)$ and $\chi_{2}^{\ast}(t)$ once the concrete form of
$\alpha_{R}^{1}(t)$ is specified.

To go further, we need to specify the initial reduced density operator
$\rho_{S}(0)$. Generally, one can always write it as
$\hat{\rho}_{S}(0)=\sum_{i}|\psi_{1i}\rangle\langle\psi_{2i}|$. For
simplicity, in the following we will assume that
$\hat{\rho}_{S}(0)=|\psi_{1}\rangle\langle\psi_{2}|$, and the extension to the
general cases just requires a simple summation. Substituting the initial
condition and Eq. (14) into Eq. (4), we can recast Eq. (4) into a compact form

$\rho_{S}(x_{f},x_{i},t)=\int\mathcal{D}^{2}\chi_{1}\int\mathcal{D}^{2}\chi_{2}^{\ast}P[\chi_{1}(t),\chi_{2}^{\ast}(t),t]\langle
x_{f}|\psi_{\chi_{1}}(t)\rangle\langle\psi_{\chi_{2}}(t)|x_{i}\rangle.$ (16)

The expressions of the forward stochastic wavefunction
$|\psi_{\chi_{1}}(t)\rangle$ and the backward stochastic wavefunction
$|\psi_{\chi_{2}}(t)\rangle$ are given by

$\begin{split}|\psi_{\chi_{1}}(t)\rangle&=\int\mathcal{D}x(t)e^{iS_{0}[x(t)]-i\int_{0}^{t}x(\tau)\left\\{\chi_{1}(\tau)-i\int_{0}^{\tau}\alpha_{res}(\tau-\tau^{\prime})x(\tau^{\prime})\mathrm{d}\tau^{\prime}\right\\}\mathrm{d}\tau}|\psi_{1}\rangle\\\
&=\hat{\mathcal{U}}_{0}(t)\mathcal{T}_{X}^{\leftarrow}\left[e^{-i\int_{0}^{t}\hat{X}(\tau)\left\\{\chi_{1}(\tau)-i\int_{0}^{\tau}\alpha_{res}(\tau-\tau^{\prime})\hat{X}(\tau^{\prime})\mathrm{d}\tau^{\prime}\right\\}\mathrm{d}\tau}\right]|\psi_{1}\rangle,\end{split}$
(17a)
$\begin{split}|\psi_{\chi_{2}}(t)\rangle&=\int\mathcal{D}x(t)e^{iS_{0}[x(t)]-i\int_{0}^{t}x(\tau)\left\\{\chi_{2}(\tau)-i\int_{0}^{\tau}\alpha_{res}(\tau-\tau^{\prime})x(\tau^{\prime})\mathrm{d}\tau^{\prime}\right\\}\mathrm{d}\tau}|\psi_{2}\rangle\\\
&=\hat{\mathcal{U}}_{0}(t)\mathcal{T}_{X}^{\leftarrow}\left[e^{-i\int_{0}^{t}\hat{X}(\tau)\left\\{\chi_{2}(\tau)-i\int_{0}^{\tau}\alpha_{res}(\tau-\tau^{\prime})\hat{X}(\tau^{\prime})\mathrm{d}\tau^{\prime}\right\\}\mathrm{d}\tau}\right]|\psi_{2}\rangle.\end{split}$
(17b)

The second equality holds by realizing the fact that any product of multiple
time-dependent operators is automatically time-ordered in the path integral
formalism[116, 117, 118]:

$\int_{x_{i}}^{x_{f}}\mathcal{D}x(t)e^{iS_{0}[x(t)]}x(t_{n})\cdots x(t_{1})\\\
=\langle
X_{f}|\mathcal{T}_{X}^{\leftarrow}[\hat{X}(t_{n})\cdots\hat{X}(t_{1})]|X_{i}\rangle,\quad
t_{f}>t_{n},\cdots,t_{1}>t_{i}$ (18)

where
$\hat{X}(t)=\hat{\mathcal{U}}_{0}^{\dagger}(t)\hat{x}\hat{\mathcal{U}}_{0}(t)$
and $|X_{i}\rangle=\hat{\mathcal{U}}_{0}^{\dagger}(t_{i})|x_{i}\rangle$ with
the definition of the pure system propagator
$\hat{\mathcal{U}}_{0}(t)=e^{-i\hat{H}_{S}t}$. $\mathcal{T}_{X}^{\leftarrow}$
represents the chronological time-ordering operation of
$\hat{X}(t_{n})\cdots\hat{X}(t_{1})$.

When an enough number of samplings of $|\psi_{\chi_{1}}(t)\rangle$ and
$|\psi_{\chi_{2}}(t)\rangle$ have be accomplished, we can numerically evaluate
$\hat{\rho}_{S}(t)$ via a statistical ensemble average according to Eq. (16).
But before that, by taking the time derivative of Eq. (17), we obtain the
following equation of motion in the operator form

$i\frac{\partial}{\partial
t}|\psi_{\chi_{j}}(t)\rangle=(\hat{H}_{S}+\chi_{j}(t)\hat{x})|\psi_{\chi_{j}}(t)\rangle+\hat{x}\hat{\mathcal{U}}_{0}(t)\hat{\mathcal{R}}(\chi_{j},t)|\psi_{j}\rangle\quad
j=1,2.$ (19)

with

$\hat{\mathcal{R}}(\chi_{j},t)=\mathcal{T}_{X}^{\leftarrow}\left[-i\int_{0}^{t}\alpha_{res}(t-\tau_{1})\hat{X}(\tau_{1})\mathrm{d}\tau_{1}e^{-i\int_{0}^{t}\hat{X}(\tau_{2})\left\\{\chi_{j}(\tau_{2})-i\int_{0}^{\tau_{2}}\alpha_{res}(\tau_{2}-\tau_{3})\hat{X}(\tau_{3})\mathrm{d}\tau_{3}\right\\}\mathrm{d}\tau_{2}}\right].$
(20)

Obviously, the hardcore issue of solving Eq. (19) resides in the treatment of
the inhomogeneous term
$\hat{x}\hat{\mathcal{U}}_{0}(t)\hat{\mathcal{R}}(\chi_{j},t)|\psi_{j}\rangle$,
which seems cumbersome at the first sight. In the following subsections, we
will provide a few solutions using either the hierarchical technique or the
perturbation expansion. The former leads to the numerically exact hierarchical
form of stochastic Schrödinger equations, while the latter serves as the
starting point to obtain a set of perturbative NMSSEs, which differ in the
choice of $\alpha_{res}(t)$ and the level of approximations.

### Hierarchical Expansion

When considering the cases where

$\alpha_{res}(t)=\sum_{m=1}^{M}d_{m}e^{-\omega_{m}t},$ (21)

one is able to construct a complete group consisting of an infinite number of
auxiliary wavefunctions with the definition as

$|\Psi_{\chi_{j}}^{(\textbf{\emph{l}})}(t)\rangle=\hat{\mathcal{U}}_{0}(t)\mathcal{T}_{X}^{\leftarrow}\left[\prod_{m=1}^{M}\left(\int_{0}^{t}d_{m}e^{-\omega_{m}(t-\tau_{1})}\hat{X}(\tau_{1})\mathrm{d}\tau_{1}\right)^{l_{m}}e^{-i\int_{0}^{t}\hat{X}(\tau_{2})\left\\{\chi_{j}(\tau_{2})-i\int_{0}^{\tau_{2}}\alpha_{res}(\tau_{2}-\tau_{3})\hat{X}(\tau_{3})\mathrm{d}\tau_{3}\right\\}\mathrm{d}\tau_{2}}\right]|\psi_{j}\rangle.$
(22)

which is closed with respect the time-derivative operation within this group.
Besides, every frequency component is associated with an index number, leading
to a vectorial index $\textbf{\emph{l}}=(l_{1},l_{2},\cdots,l_{M})$ with
$l_{m}$ being the non-negative integers. Directly taking the time derivative
of Eq. (22), one obtains the following hierarchical set of differential
equations of motion in a general form[108, 112, 111]:

$i\frac{\partial}{\partial
t}|\Psi_{\chi_{j}}^{(\textbf{\emph{l}})}(t)\rangle=\left(\hat{H}_{S}+\chi_{j}(t)\hat{x}-i\sum_{m=1}^{M}l_{m}\omega_{m}\right)|\Psi_{\chi_{j}}^{(\textbf{\emph{l}})}(t)\rangle+\sum_{m=1}^{M}l_{m}d_{m}\hat{x}|\Psi_{\chi_{j}}^{(\textbf{\emph{l}}_{m}^{-})}(t)\rangle+\hat{x}\sum_{m=1}^{M}|\Psi_{\chi_{j}}^{(\textbf{\emph{l}}_{m}^{+})}(t)\rangle,$
(23)

where $\textbf{\emph{l}}_{m}^{\pm}=(l_{1},\cdots,l_{m}\pm 1,\cdots,l_{M})$. In
the numerical simulations, it is truncated at a certain level
$L=\sum_{m=1}^{M}l_{m}$. But obviously, the larger M is, the more
wavefunctions are needed. In the end, only the zeroth-tier term with
$\textbf{\emph{l}}=(0,\cdots,0)$ are preserved for the stochastic ensemble
average. By inspecting Eq. (22) and setting $t=0$, we obtain the initial
condition for Eq. (23):

$|\Psi_{\chi_{j}}^{(\textbf{\emph{l}})}(0)\rangle=\left\\{\begin{array}[]{l}|\psi_{j}\rangle,\quad\textbf{\emph{l}}=(0,\cdots,0);\\\
0,\quad else.\end{array}\right.$ (24)

In the condensed phase, the bath consisting of nuclear vibrational modes and
solvent degrees of freedom usually exhibits the overdamped features, thereby
the spectral density can be represented by an Ohmic form with the Lorentzian
cutoff[119] (often called the Debye or Drude spectral density in the
literatures)

$J(\omega)=\frac{2\lambda\omega\omega_{c}}{\omega^{2}+\omega_{c}^{2}},$ (25)

where $\lambda$ and $\omega_{c}$ are the reorganization energy and the cutoff
frequency, respectively. In this case, the corresponding bath correlation
function can be analytically obtained as a series of exponential decaying
functions with respect to time[119], satisfying the condition in Eq. (21).

$\alpha(t)=\lambda\omega_{c}\coth{\frac{\beta\omega_{c}}{2}}e^{-\omega_{c}t}-\sum_{v=1}^{\infty}\frac{4\lambda\omega_{c}\omega_{v}}{\beta(\omega_{c}^{2}-\omega_{v}^{2})}e^{-\omega_{v}t}-i\lambda\omega_{c}e^{-\omega_{c}t},$
(26)

where $\omega_{v}=2\pi v/\beta$ are the Matsubara frequencies. In the
numerical implementations, a constraint on the upper bound of $\omega_{v}$ is
indispensable in order to keep a finite index number. Nevertheless, it grows
rapidly with the decreasing temperature. Although a much more efficient
frequency decomposition scheme, called Padé spectrum decomposition has been
proposed to alleviate this problem to some extent[46], the numerical cost is
still expensive, especially when the temperature is pretty low. Upon closer
inspection of Eq. (26), it found that all the temperature effects are included
in the real part. Therefore, we prefer to entirely circumvent this problem by
choosing $g(\omega)=1$ in Eq. (10), as a result, we have

$\alpha_{res}(t)=i\alpha_{I}(t)=-i\lambda\omega_{c}e^{-\omega_{c}t}.$ (27)

In addition, an efficient generation scheme for $\chi_{1}$ and
$\chi_{2}^{\ast}$ are introduced:

$\chi_{1}(t)=\xi_{c}(t)+\xi_{1}(t),$ (28a)
$\chi_{2}^{\ast}(t)=\xi_{c}^{\ast}(t)+\xi_{2}(t)$ (28b)

with the explicit expressions

$\xi_{c}(t)=\sum_{k}\sqrt{\frac{c_{k}^{2}}{2\omega_{k}}}\left[\sqrt{\frac{n_{k}+1}{2}}(\mu_{k}^{1}+i\mu_{k}^{2})e^{i\omega_{k}t}+\sqrt{\frac{n_{k}}{2}}(\mu_{k}^{1}-i\mu_{k}^{2})e^{-i\omega_{k}t}\right],$
(29a)
$\xi_{1}(t)=\sum_{k}\sqrt{\frac{c_{k}^{2}}{2\omega_{k}}}(\sqrt{n_{k}+1}-\sqrt{n_{k}})\times(\mu_{k}^{3}\cos{\omega_{k}t}+\mu_{k}^{4}\sin{\omega_{k}t}),$
(29b)
$\xi_{2}(t)=\sum_{k}\sqrt{\frac{c_{k}^{2}}{2\omega_{k}}}(\sqrt{n_{k}+1}-\sqrt{n_{k}})\times(\mu_{k}^{5}\cos{\omega_{k}t}+\mu_{k}^{6}\sin{\omega_{k}t}),$
(29c)

where $n_{k}=\bar{n}(\omega_{k})$, and $\mu_{k}^{j}$ ($j=1,\cdots,6$) are
independent Gaussian white noises obey the normal distribution $N(0,1)$.
Making use of the fact that $\left\langle\mu_{k}^{i}\right\rangle=0$ and
$\left\langle\mu_{k}^{i}\mu_{k^{\prime}}^{j}\right\rangle=\delta_{ij}\delta_{kk^{\prime}}$,
one can find that the generation scheme Eq. (28) and Eq. (29) indeed satisfies
Eq. (15) with $\alpha_{R}^{1}(t)$ being $\alpha_{R}(t)$. Finally, Eq. (23) is
reduced to the following form:

$i\frac{\partial}{\partial
t}|\Psi_{\chi_{j}}^{(n)}(t)\rangle=\left(\hat{H}_{S}+\chi_{j}(t)\hat{x}-in\omega_{c}\right)|\Psi_{\chi_{j}}^{(n)}(t)\rangle+\hat{x}|\Psi_{\chi_{j}}^{(n+1)}(t)\rangle-
in\lambda\omega_{c}\hat{x}|\Psi_{\chi_{j}}^{(n-1)}(t)\rangle,\quad
n=0,1,2,\cdots$ (30)

Eq. (30) possesses the merits of satisfactory statistical convergence
performance, simplest hierarchical structure, and the favorable scaling
property of Hilbert space, and so on. Thus, it is suitable to numerically
exactly and efficiently explore the excitonic dynamics in Fenna-Matthews-Olson
(FMO) trimer complexes of green sulfur bacteria[112] and peripheral light-
harvesting complex 2 of the purple bacteria[120], and other intermediate-sized
systems. A recent study[121] showed that the method exhibits its strengths in
simulating the system dynamics embedded in an ultraslow bath, where the non-
Markovianity has proven to be extremely strong. In reality, there are a few
numerical calculations[112, 121] demonstrating that the truncation level of
hierarchical stochastic Schrödinger equation is generally less than that of
the hierarchical equations of motion. One reason behind is that a majority of
environmental influence has been taken into account through the introduction
of stochastic fields, such that the higher-order phonon effects contained in
the higher-order auxiliary wavefunctions are comparatively small. This has
inspired us to overcome the limitations inherited in the hierarchy formalism,
for example, the factorial scaling with respect to the system degrees of
freedom and a few available spectral density functions, by resorting to the
perturbative treatment[113].

### Perturbation Expansion

Let’s refocus on Eqs. (19) by reformulating it into a formally exact time-
convolutionless equation of motion:

$i\frac{\partial}{\partial
t}|\psi_{\chi_{j}}(t)\rangle=\left[\hat{H}_{S}+\chi_{j}(t)\hat{x}+\hat{x}\hat{\mathcal{U}}_{0}(t)\hat{\mathcal{K}}(\chi_{j},t)\hat{\mathcal{U}}_{0}^{\dagger}\right]|\psi_{\chi_{j}}(t)\rangle,$
(31)

where the so-called residual dissipative generator
$\hat{\mathcal{K}}(\chi_{j},t)$ reads

$\begin{split}\hat{\mathcal{K}}(\chi_{j},t)=&\mathcal{T}_{X}^{\leftarrow}\left[-i\int_{0}^{t}\alpha_{res}(t-\tau_{1})\hat{X}(\tau_{1})\mathrm{d}\tau_{1}e^{-i\int_{0}^{t}\hat{X}(\tau_{2})\left\\{\chi_{j}(\tau_{2})-i\int_{0}^{\tau_{2}}\alpha_{res}(\tau_{2}-\tau_{3})\hat{X}(\tau_{3})\mathrm{d}\tau_{3}\right\\}\mathrm{d}\tau_{2}}\right]\\\
&\times\mathcal{T}_{X}^{\rightarrow}\left[e^{i\int_{0}^{t}\hat{X}(\tau_{1})\left\\{\chi_{j}(\tau_{1})-i\int_{0}^{\tau_{1}}\alpha_{res}(\tau_{1}-\tau_{2})\hat{X}(\tau_{2})\mathrm{d}\tau_{2}\right\\}\mathrm{d}\tau_{1}}\right],\end{split}$
(32)

where $\mathcal{T}_{X}^{\rightarrow}$ denotes the anti-chronological time-
ordering operation. This expression is the starting point of the systematic
perturbation expansion with respect to the system-bath coupling strength
$c_{k}$. It can be seen from Eq. (7), Eq. (8) and Eq. (15) that $\chi_{j}(t)$
is of the first order of $c_{k}$, while $\alpha(t)$ is of the second order. By
expanding the two exponential functions in Eq. (32) directly and paying
particular attention to the time-ordering operation, we can obtain the
perturbation expansion of Eq. (32) up to arbitrary order as
$\hat{\mathcal{K}}(\chi_{j},t)=\sum_{n=2}^{\infty}\hat{\mathcal{K}}_{n}(\chi_{j},t)$.
The lowest (second) order therein is explicitly expressed as

$\hat{\mathcal{K}}_{2}(t)=-i\int_{0}^{t}\alpha_{res}(t-\tau)\hat{X}(\tau)\mathrm{d}\tau.$
(33)

As a matter of fact, It is advisable to approximate the residual dissipative
generator as $\hat{\mathcal{K}}(\chi_{j},t)\approx\hat{\mathcal{K}}_{2}(t)$,
as we have affirmed in a numerically manner that the contributions of higher
order terms are generally small over a quite broad parameter space[113].
Inserting Eq. (33) into Eq. (31) and then making the substitution
$t-\tau\rightarrow\tau$, we obtain the final result of this subsection:

$i\frac{\partial}{\partial
t}|\psi_{\chi_{j}}(t)\rangle=\left[\hat{H}_{S}+\chi_{j}(t)\hat{x}-i\hat{x}\int_{0}^{t}\alpha_{res}(\tau)e^{-i\hat{H}_{S}\tau}\hat{x}e^{i\hat{H}_{S}\tau}\mathrm{d}\tau\right]|\psi_{\chi_{j}}(t)\rangle.$
(34)

A significant point about Eq. (34) is that it is applicable to arbitrary form
of spectral densities, which is very important since more and more studies
have stressed that the realistic highly structured environment might play a
crucial role in various cases like the enhancement of energy transfer
rate[122, 123, 124, 125] and non-Markovianity[126, 127, 128], the long-
sustained quantum coherence in light-harvesting complexes[129, 130, 131, 132],
and the fine-tuning of the dephasing and relaxation time in quantum
information processors[133, 134]. What’s more, since the concrete form of
$\alpha_{res}(t)$ which exerts a substantial influence on the numerical
performance of Eq. (34), has not yet been specified, this flexibility allows
us to seek for a balance point between high accuracy and affordable
computational cost.

After tons of numerical simulations based on the spin-boson model over a broad
parameter range, it is found that an excellent numerical performance is
guaranteed by the following choice of the pairwise $\chi_{j}(t)$ and
$\alpha_{res}(t)$:

$\chi_{1}(t)=\chi_{2}(t)=\xi_{c}(t),\quad\text{as shown in Eq.\,(\ref{sto-
HSSE-2-a}) };$ (35a)
$\alpha_{res}(t)=\int_{0}^{\infty}\mathrm{d}\omega\frac{J(\omega)}{\pi}[\tanh{\frac{\beta\omega}{4}}\cos{\omega
t}-i\sin{\omega t}].$ (35b)

Specifically, we have denoted Eq. (34) together with the definitions of
$\chi_{1/2}(t)$ and $\alpha_{res}(t)$ shown in Eq. (35) as the perturbative
stochastic Schrödinger equation (PSSE) in the coming numerical section. The
numerical efficiency of the PSSE method has been displayed by comparing with
that of other versions of SSEs, benchmarked against the exact results from Eq.
(23) for a wide range of parameters in the spin-boson model[113], and it was
found to be relatively stable from weak to intermediate system-bath coupling
regimes at different temperatures and bath characteristic frequencies. We will
further affirm its valid regime in a simple dimer model before it is applied
to realistic organic aggregates in this paper. In the past few years, the PSSE
method has proven its strength in simulating charge carrier dynamics in large-
scale realistic organic photovoltaic systems composed of a few hundreds
(thousands) monomers[101, 102, 13, 103, 104, 105, 106, 107]. Some specific
examples will be detailedly presented in the next section.

In some limiting cases, for instance, highly-symmetric systems and extremely
high-temperature conditions, it is reasonable to further invoke a crude
approximation[135], i.e.,

$\alpha(t)\approx\alpha_{R}(t),$ (36)

which actually implies a semiclassical assumption. Under this circumstance,
Eq. (31) is greatly simplified to

$i\frac{\partial}{\partial
t}|\psi_{\chi_{j}}(t)\rangle=\left(\hat{H}_{S}+\widetilde{\xi}_{c}(t)\hat{x}\right)|\psi_{\chi_{j}}(t)\rangle,$
(37)

since $\alpha_{res}(t)=0$ and the stochastic process $\widetilde{\xi}_{c}(t)$
can be reduced to the real one

$\widetilde{\xi}_{c}(t)=\sum_{k}\sqrt{\frac{c_{k}^{2}}{\omega_{k}}\coth{\frac{\beta\omega_{k}}{2}}}\left(\mu_{k}^{1}\cos(\omega_{k}t)+\mu_{k}^{2}\sin(\omega_{k}t)\right),$
(38)

where $\mu_{k}^{1}$ and $\mu_{k}^{2}$ are independent real Gaussian random
variables obeying the standard normal distribution, and it is trivial to
verify that $\left\langle\widetilde{\xi}_{c}(t)\right\rangle=0$ and
$\left\langle\widetilde{\xi}_{c}(t)\widetilde{\xi}_{c}(t^{\prime})\right\rangle=\alpha_{R}(t-t^{\prime})$.
If the bath response time is ultrashort in comparison with the system
dynamics, one can proceed to approximate $\alpha_{R}(t)$ as a delta function
(white noise), such that Eq. (37) recovers the classical Haken-Strobl-Reineker
model (HSR)[136, 137, 138], and out of this reason, Eq. (37) is referred to as
the modified-HSR method in the following numerical simulations. Even though
the neglect of $\alpha_{I}(t)$ leads to the break down of detailed balance
condition, as long as one can take care of its range of application in the
practical simulations, Eq. (37) still constitutes a valuable and reliable
simulation tool widely-used in nanoscale systems. For instance, it is utilized
in evaluating and predicting the charge carrier mobilities in many
representative organic materials[139, 140], and is used to investigate the
effects of spatial as well as temporal correlations between site energies on
the charge transfer (CT) process in flexible molecules[141]. More
interestingly, Eq. (37) toghether with the help of massively parallel
computing platforms[142] and advanced experimental techniques, is capable to
set up a bottom-up strategy[143] promising for the optimal design of advanced
organic functional materials.

## Applications

In this section, we will present a few representative applications of the
aforementioned methods. Various numerical algorithms like the forth-order
Runge-Kutta method and Chebyshev polynomial expansion are applicable to the
propagation of stochastic Schrödinger equations. It should be noted that the
norm of the stochastic wavefunction is actually not conserved during the
propagation, so an artificial renormalization is generally performed in the
practical simulations. However, it is possible to transform the linear
stochastic Schrödinger equation into a norm-conserved nonlinear stochastic
Schrödinger equation by means of Girsanov transformation. Various nonlinear
forms of Eq. (23) and Eq. (34) can be found in the literatures[108, 112, 92,
95, 105]. In the hierarchy of pure state method (HOPS)[108], the non-linear
version is reported to remarkably accelerate the convergence speed, while this
conclusion does not apply to Eq. (30)[112]. Concerning the fact that linear
SSE is computationally beneficial than its non-linear correspondence,
especially for large systems, we are prone to using the linear form in the
following simulations.

Up to now, the hierarchical form of stochastic Schrödinger equations have been
used in simulating the exciation energy transfer in full 24-site FMO
complexes[112], the linear spectra of light-harvesting complex 2[120], as well
as investigating the non-Markovianity of a system embedded in an ultraslow
bath[121]. In the following, we will demonstrate in great details the
superiority of the HSSE method over its deterministic counterpart, the HEOM
method with regards to the simulations of excitation energy transfer in large-
scale photosynthetic systems.

#### Energy Transfer in the FMO Complex

Fenna-Matthews-Olson complex is a light harvesting antennae found in green
sulfur bacteria, bridging a large chromosome supercomplex to the reaction
center. Due to its relatively small size and well-resolved network structure,
it has long served as a prototypical model to explore the highly efficient
energy transfer in photosynthesis and has attracted tremendous attention from
numerous research areas, including various experimental and theoretical
spectroscopic studies [1, 2, 5, 144, 63, 145, 57, 146], atomistic
simulations[147, 148, 149], electronic structure and quantum dynamics
calculations[61, 150, 55, 127, 151, 152, 96], and so on. Two-dimensional
electronic spectroscopy studies of FMO complex during the past decade have
brought an unprecedentedly high time- and frequency-resolution of the
underlying ultrafast physical processes, while the consensus about the
interpretation of observed unexpectedly long-lived peak oscillating behaviors
has not yet been reached[2, 5, 129, 130, 57, 131, 61, 132, 146, 74, 153, 154].
Thus, in the quantum dynamics branch, it would be preferable to use the
methods as accurate as possible. The first non-perturbative quantum dynamics
calculation of the excitation energy transfer in FMO complex was performed by
Ishizaki and Fleming[61] using the HEOM approach, based on a well-established
excitonic Hamiltonian for a seven-bacteriochlorophyll subunit out of three-
fold circularly symmetric assembly. It is further assumed that every pigment
is coupled to its own environment which is characterized by the Debye spectral
density function, as shown in Eq. (25). Based on the same model, we would
exemplify the applicability and efficiency of the HSSE method, in comparison
with the popular HEOM method.

Figure 1: Plots of time-evolving populations at BChl 1, 2, 3, and 4 for a seven-site subunit of FMO complex, under four different parameter sets with the specific values given at the top of every panel. The reorganization energy is fixed at $\lambda=35\,\mathrm{cm}^{-1}$. The initial population is located at BChl 1. Parameter Set | HEOM | HSSE
---|---|---
$(\omega_{c}^{-1},\,\mathrm{T})$ | $(\mathrm{L_{HEOM}},\,\mathrm{K})$ | (#AM)$\times\mathrm{N}^{2}$ | $(\mathrm{L_{HSSE}},\,\mathrm{N_{ran}})$ | (#AW)$\times$ 2N
$(50\mathrm{fs},\,77\mathrm{K})$ | $(4,\,3)$ | 619850 | $(3,\,20000)$ | 1680
$(50\mathrm{fs},\,300\mathrm{K})$ | $(5,\,1)$ | 38808 | $(3,\,10000)$ | 1680
$(166\mathrm{fs},\,77\mathrm{K})$ | $(6,\,2)$ | 1899240 | $(3,\,10000)$ | 1680
$(166\mathrm{fs},\,300\mathrm{K})$ | $(10,\,1)$ | 952950 | $(3,\,5000)$ | 1680
Table 1: The computational requirements for the HEOM and HSSE methods to obtain the converged results as shown Fig. 1. $\mathrm{L_{HEOM}}$ and $\mathrm{L_{HSSE}}$ are the truncation levels of hierarchy for HEOM and HSSE method, respectively. $\mathrm{K}$ is the number of Matsubara frequencies needed and #AM denotes the number of auxiliary matrices in HEOM method. $\mathrm{N_{ran}}$ is the random trajectory number to obtain smooth curves and #AW denotes the number of auxiliary wavefunctions used in HSSE method. Besides, $\mathrm{N}=7$ in this case. Parameter Set | HEOM | HSSE
---|---|---
$(\omega_{c}^{-1},\,\mathrm{T})$ | $(\mathrm{L_{HEOM}},\,\mathrm{K})$ | (#AM)$\times\mathrm{N}^{2}$ | $(\mathrm{L_{HSSE}},\,\mathrm{N_{ran}})$ | (#AW)$\times$ 2N
$(50\mathrm{fs},\,77\mathrm{K})$ | $(4,\,3)$ | 738993600 | $(3,\,20000)$ | 140400
$(50\mathrm{fs},\,300\mathrm{K})$ | $(5,\,1)$ | 68402880 | $(3,\,10000)$ | 140400
$(166\mathrm{fs},\,77\mathrm{K})$ | $(6,\,2)$ | 14876447040 | $(3,\,10000)$ | 140400
$(166\mathrm{fs},\,300\mathrm{K})$ | $(10,\,1)$ | 75529808640 | $(3,\,5000)$ | 140400
Table 2: The computational requirements for the HEOM and HSSE methods to
obtain the converged results for full 24-site FMO trimer. The meaning of the
characters are the same as those in Table. 1 except that $\mathrm{N}=24$ in
this case.
Figure 2: Plots of time-evolving populations for a full 24-site FMO complex,
under four different parameter sets with the specific values given at the top
of every panel. The reorganization energy is fixed at
$\lambda=35\,\mathrm{cm}^{-1}$. The initial population is located at BChl 1
within one subunit out of a homotrimer. For the purpose of a better
visualization, only those corresponding to the most populated eight BChls are
displayed.

Shown in Fig. 1 are the population evolutions of the first four
bacteriochlorophylls (BChl) obtained through both the HEOM and HSSE methods,
under four different parameter sets of the temperature $\beta$ and bath cutoff
frequency $\omega_{c}$, while the bath reorganization energy is fixed at
$\lambda=35\,\mathrm{cm}^{-1}$. The initial excitation is prepared at BChl 1.
The numerical conditions for convergence are listed in Table. 1. Note that the
higher the temperature or the smaller $\omega_{c}$ is, the deeper the
truncation level $\mathrm{L_{HEOM}}$ is needed in HEOM method. $\mathrm{K}$
will be larger as well at the lower temperature. In contrast,
$\mathrm{L_{HSSE}}$ seems less susceptible to the temperature and bath cutoff
frequency[121]. Besides, as mentioned in previous section about HSSE (Eq.
(30)), all the temperature effects are included in the stochastic noises and
consequently neither Matsubara nor Padé frequency decomposition schemes are
needed. In the end, the number of elements needed to be propagated in HSSE is
far smaller than that in HEOM method, which indicates a substantial saving of
memory requirements. Typically in photosynthetic systems or organic
aggregates, the electron-phonon interaction is comparatively small or in close
proximity to the excitonic couplings, such that the random trajectory number
$\mathrm{N_{ran}}$ in most cases remains manageable. In passing, the
stochastic nature renders HSSE method well-suited for parallel computing and
allows for the inclusion of static disorders in a trivial fashion[155]. The
advantages will be more evident in larger complexes.

Very recently, it is found that there exits in vivo an additional BChl, termed
as BChl 8, residing in a cleft at the surface of the protein to connect two
neighboring subunits[156, 157, 147]. As such, one has to consider the intact
trimer structure composed of 24 BChls. The computational requirements grow a
lot for both HEOM and HSSE methods, but obviously more demanding for the
former, as illustrated in Table. 2. The results of population evolutions of
full trimer at the bath condition $\omega_{c}^{-1}=50\,\mathrm{fs}$ and
$\mathrm{T}=300\,\mathrm{K}$ are calculated by Wilkins and Dattani[54] through
the HEOM method with the usage of a more efficient numerical propagation
algorithm. Nevertheless, applications to the rest three parameter sets are far
more challenging. For the most straightforward example, at the case
$\omega_{c}^{-1}=166\,\mathrm{fs}$ and $\mathrm{T}=300\,\mathrm{K}$, up to 75
billions elements are needed to be propagated in a single time step. While,
complementarily, the numerical expenses for the HSSE method in this parameter
regime appear to be the cheapest. The converged results acquired by the HSSE
method under four different sets of bath parameters are shown in Fig. 2, with
the excitonic Hamiltonian give by

$\hat{H}_{ex}=\left[\begin{array}[]{ccc}\hat{H}^{1}_{A}&\hat{H}_{B}&\hat{H}_{B}^{T}\\\
\hat{H}_{B}^{T}&\hat{H}^{1}_{A}&\hat{H}_{B}\\\
\hat{H}_{B}&\hat{H}_{B}^{T}&\hat{H}^{1}_{A}\\\ \end{array}\right],$ (39)

$H_{A}^{1}$ covers the intra-subunit excitonic elements[158] and $H_{B}$
represents all the inter-subunit couplings[147], and they are specifically
tabulated in Ref. [112]. In addition, we have also plotted the results
obtained through the PSSE method, i.e., Eq. (34), in Fig. 2. Benchmarked by
the exact results, the PSSE method has proven to be remarkably accurate for
delineating the population relaxation in FMO complexes. In fact, the
simulations using modified-HSR method (Eq. (37)) are also performed, but the
results are qualitatively wrong and then not shown here. Since Eq. (34) and
Eq. (37) are applicable to nanoscale realistic systems consisting of several
thousands and tens of thousands of system degrees of freedom, respectively. It
is quite essential to systematically investigate their ranges of validity
beforehand.

#### Validity Regime of Two Approximate SSEs

Figure 3: Plots of the time-evolving populations at site 1 in a dimer model
system, subjected to various reorganization energies $\lambda$
$\left[(1)-(4)\right]$ and two different energy biases $\big{[}\Delta
E=E_{1}-E_{2}=0$ and $\Delta E=100\,\mathrm{cm}^{-1}\big{]}$ at both high and
low temperatures $\big{[}\mathrm{T}=300\,\mathrm{K}$ and
$\mathrm{T}=50\,\mathrm{K}\big{]}$. The initial population is put at site 1.
Other parameters are set as: the excitonic coupling $V=100\,\mathrm{cm}^{-1}$,
$\omega_{c}^{-1}=100\,\mathrm{fs}$ in the Debye spectral density function Eq.
(25).

The perturbative forms of stochastic Schrödinger equations, including Eq. (34)
and Eq. (37), due to their favorable linear scaling to the system size, well-
behaved stochastic convergence property, and the time-convolutionless
structure, have proven to be extremely suitable for describing the quantum
dynamics of various large-scale organic systems and have achieved remarkable
success, for example, successfully illuminating the temperature dependence of
carrier diffusion coefficients in one-dimensional chain models including both
static and dynamic disorders for a wide range of disorder strength, and
bridging the band-like and hopping-type transport[101], demonstrating the
effects of the laser pulse, the Coulomb interaction and the temperature in the
energy relaxation as well as the separation of a photon-induced hot electron-
hole pair in organic aggregates[102]. Besides, it can be straightforwardly
extended to investigate the mixed Frenkel exciton and CT exciton states, and
then demonstrated the important role of hole delocalization on the Frenkel-CT
decoherence time[103].

It is recommendable to testify the validity of approaches on the basis of a
simple model subjected to various parameter conditions. The simplest molecular
aggregate, namely, the dimer, is adopted here, while some key features should
be shared by larger analogous systems. As was found previously in Ref. [113]
and also demonstrated in the FMO complex shown above, the accuracy of the PSSE
and modified-HSR methods seems barely influenced by the cutoff frequency. As
such, we have fixed it at the value $\omega_{c}^{-1}=100\,\mathrm{fs}$ in the
calculations. The excitonic coupling is chosen as a reference and set as
$V=100\,\mathrm{cm}^{-1}$. In a homogeneous sample with the energy biases
being zeros, as shown in Fig. 3 (a) and (b), the results of both PSSE and
modified-HSR methods are in good agreement with the exact ones, even when the
reorganization energy $\lambda$ is comparable to or slightly larger than the
excitonic couplings. In fact, the results remain qualitatively solid till
$\lambda=3\,V$. Nevertheless, this range is greatly reduced when the finite
energy bias $(\Delta E=E_{1}-E_{2}\neq 0)$ is introduced. It could be
concluded that the modified-HSR method is completely not suited for an
inhomogeneous system. As can clearly seen in Fig. 3 (c) and (d), it fails even
under very weak system-bath coupling condition $(\lambda/V=0.1)$, and the
reason behind is that the imaginary part of bath correlation function is
completely neglected, that is, all the feedbacks from the environment to the
system are missed and thus the detailed balance is broken down. The PSSE
method has partially repaired this disability, especially at the low
temperatures. In short, the PSSE method might overestimate the relaxation
rates when applied to an inhomogenous system, but generally speaking, the
results should be trustworthy when $\lambda/V<1$.

The rest two subsections mainly introduce the applications of the PSSE method
(Eq. (34)) to unveil the multiple time scales hot exciton relaxation process
in a neat organic aggregate[13], and study the CT states, exciton migrations,
and the quantum interference effects in singlet fission (SF) dynamics[106].

#### Hot Exciton Relaxation

The energy relaxation of a hot exciton (an excitonic state with energy
exceeding $k_{B}T$ considerably) is believed to be critical to many dynamic
processes in organic photovoltaics like charge separation and multiple exciton
generation. Recent experiments found that the excess energy of the hot exciton
goes through an initial fast followed by a slow relaxation process[159, 2, 9,
160, 161, 162, 163, 164]. Further studies indicate that the fast one, which
may be relevant to the coherent motion of the exciton, can largely enhance the
generation of separated charges at the donor-acceptor interface[12]. However,
the detailed mechanism of the hot exciton energy relaxation is still not
clear. A recent work[13] have given an unambiguous picture of this process and
clarified the ultimate origin of the multiple time scales processes.

Here we consider the Frenkel exciton model for an one-dimensional aggregate
chain. The Hamiltonian of the exciton part is

$\hat{H}_{ex}=\sum_{n}\epsilon_{n}\hat{B}_{n}^{\dagger}\hat{B}_{n}+\sum_{n>m}V_{nm}(\hat{B}_{n}^{\dagger}\hat{B}_{m}+\hat{B}_{m}^{\dagger}\hat{B}_{n}),$
(40)

where $\hat{B}_{n}^{\dagger}$ and $\hat{B}_{n}$ are the creation and
annihilation operator of a local exciton at the nth site, and $\epsilon_{n}$
is the corresponding energy. $V_{nm}$ is the excitonic coupling between mth
and nth sites, which can be calculated through the Coulomb force of transition
electron density[165]. Here we characterize the influence of the bath on the
excitonic system by an Ohmic spectral density function
$J(\omega)=\frac{1}{2}\pi\xi\omega e^{-\omega/\omega_{c}}$, where $\xi$ and
$\omega_{c}$ is the Kondo parameter and the cutoff frequency, respectively,
and their product $\xi\omega_{c}$ gives the reorganization energy $\lambda$.
For the initial state of the exciton, we assume a normalized Gaussian
wavepacket to mimic the laser induced hot exciton:

$f(E_{n})=\frac{1}{Z}e^{-\frac{(E_{n}-\overline{E}_{0})^{2}}{2\Delta_{E}}},$
(41)

where $E_{n}$ is the energy of the nth eigenstate of $H_{ex}$,
$\overline{E}_{0}$ and $\Delta_{E}$ represent the average energy and the
magnitude of the energy delocalization of the initial exciton, corresponding
to the pulse energy and pulse duration, respectively,
$Z^{2}=\sum_{n}e^{-\frac{(E_{n}-\overline{E}_{0})^{2}}{\Delta_{E}}}$ is the
normalization constant. From Eq. (41), we can obtain the initial wavefunction
in the site representation by a simple unitary transformation, as is shown in
Fig. 4. Several experiments indicate that the coherence length of the initial
exciton may vary from $10^{2}$ to $10^{4}$[166, 167, 168, 169], depending on
the system of interest. Thus it is reasonable to assume an initial wavepacket
with a coherence length of 100 sites in the simulation. Specifically speaking,
we consider an organic aggregate chain with 200 sites, and choose the 100
sites in the center of the chain to form initial excitonic wavepacket. We set
$\epsilon_{n}=0$ and only consider the nearest-neighbouring coupling,
$V_{i,i\pm 1}=1$ for simplicity.

Figure 4: Transform of initial distribution from the energy representation to
the site representation. The energy band ranges from $-2V_{i,i\pm 1}$ to
$2V_{i,i\pm 1}$ as the excitonic energies $\\{\epsilon_{n}\\}$ are all set as
zero. Reprinted with permission from Ref. [13]. Copyright 2015 American
Chemical Society.

When there is no coupling between the exciton and the bath, the dynamics is
purely coherent, exhibiting an interference pattern of the exciton population,
and the average energy of the exciton, defined as $\overline{E}(t)={\rm
Tr}{\\{\hat{\rho}_{ex}(t)\hat{H}_{ex}\\}}$ ($\hat{\rho}_{ex}(t)$ is the
reduced density operator of the exciton, and the energy relaxation curves of
the hot exciton under different exciton-phonon interaction are shown in Fig. 5
(a)), will keep constant in this situation. However, the interference pattern
quickly disappears when the exciton-phonon interaction is turned on, and a
larger $\lambda$ leads to a faster localization of the exciton. An interesting
feature is that now the energy relaxation shows a decay trend with multiple
time scales, which has been observed by many experiments[9, 170, 171, 172,
173, 174]. A possible two-step mechanism is proposed to understand this
phenomenon, an initial excitonic state with energy $E_{h}$ first relaxes to an
intermediate state with energy $E_{d}$ through the fast depasing dynamics, and
then relaxes to the quasi-thermal equilibrium state with energy $E_{e}$ by
slow hopping. This mechanism can be expressed as[13]

$\overline{E}(t)=E_{e}+Ae^{-k_{1}t}+Be^{-k_{2}t},$ (42)

where $A=\frac{((E_{d}-E_{h})k_{1}+(E_{h}-E_{e})k_{2})C}{k_{2}-k_{1}}$,
$B=\frac{(E_{e}-E_{d})k_{1}C}{k_{2}-k_{1}}$, $C$ is the initial population of
the exciton, $k_{1}$ and $k_{2}$ are the rates of fast dephasing and slow
hopping respectively. One thing should be noted is that when the ratio
$k_{1}/k_{2}$ is larger than $\frac{E_{e}-E_{h}}{E_{d}-E_{h}}$, $A$ becomes
negative, which has also been observed in experiments[9, 170].

To validate the above model, we further investigate the decoherence dynamics
of the exciton by calculating the coherence-length sequence$L_{k}(t)$, defined
as

$L_{k}(t)=\sum_{i=1}^{N+1-k}|\langle i|\hat{\rho}_{ex}|i+k-1\rangle|,$ (43)

where $|i\rangle$ is a brief notation representing the localized exciton state
at the ith site, and $k$ varies from $1$ to $N$, where $N$ is the number of
sites. $L_{k}(t)$ actually reflects the coherence between two sites separated
from $k-1$ sites. We calculate $L_{k}(t)$ for two different $\lambda$ and the
results are shown in Fig. 5(b). As can be seen, the dephasing is faster with a
larger $\lambda$, consistent with experiments[175]. The quantitative dephasing
time $t_{d}$ can be obtained from the averaged coherent length
$\overline{L}(t)=\sum_{k=1}^{N}L_{k}(t)/(N-k+1)$ by fitting $\overline{L}(t)$
exponentially. Fig. 5(c) shows $\overline{L}(t)$ for the two different
$\lambda$ cases, and one can find that $t_{d}$ is very close to the fast
relaxation time $t_{1}=1/k_{1}$. To further identify their relationship, we
calculate $t_{d}$ and $t_{1}$ for a wide range of $\lambda$, and the results
are shown in Fig. 5 (c). It can be seen that these two time constants are
indeed consistent with each other, validating the assumption that the fast
energy relaxation is corresponding to the dephasing of the exciton.

Figure 5: (a) The energy relaxation dynamics corresponding to two different
exciton-phonon interaction (weak: $\lambda=0.4$ and strong: $\lambda=4.0$,
respectively). (b) The dephasing process of the averaged coherent length
$\overline{L(t)}$. (c) The dephasing time ($t_{d}$) and the fast energy
relaxation time ($t_{1}$) with respect to the varying reorganization energies.
The initial wavepacket is generated by a laser pulse with the parameters
$\overline{E}_{0}=0$ ($2V_{i,i\pm 1}$ above the bottom of the energy band),
$\Delta_{E}=0.008$, and $\omega_{c}=1$, $\beta=10$. Adapted with permission
from Ref. [13]. Copyright 2015 American Chemical Society.

As an illustrative example of realistic materials, we calculate the ultrafast
exciton dynamics in a one dimensional chain of PBDTTPD copolymer, which is a
promising electron donor in organic solar cells[176, 177, 178, 179, 180, 181,
182, 183]. The structure of the PBDTTPD copolymer and one of its unit are
shown in Fig. 6, the latter is regarded as a site in the Frenkel exciton
model. All the parameters needed for the quantum dynamics calculations are
obtained from ab initio calculations, based on the Gaussian09 and Q-chem
program packages. Quantum chemistry calculations are performed on a four-unit
oligomer, and alkyl sides are replaced by methyl groups to lower the
computational cost. The semiempirical quantum chemical PM6 method is used to
optimize the structure of the oligomer. We consider the excitation energy
transfer as a reaction of $DD^{\ast}\rightarrow D^{\ast}D$. The total
reorganization energy is calculated via the vibrational mode method[184] as
the sum of the reorganization energy from the ground state $D$ to the excited
state $D^{\ast}$ and that from $D^{\ast}$ to $D$, and can be expressed as
$\lambda=\sum_{k}\lambda_{k}=\sum_{k}\frac{1}{2}\omega_{k}^{2}\Delta
Q_{k}^{2}$, where $\omega_{k}$ is the the frequency corresponding to the kth
normal mode, $\Delta Q_{k}$ represents the coordinate shift between the
optimized geometries of $D$ and $D^{\ast}$ along the kth mode. Based on the
first-order approximation in the single excitation theory, the excitonic
coupling can be calculated as[185]

$\begin{split}V_{nm}=&\langle
M_{n}^{\ast}M_{m}|\hat{H}|M_{n}M_{m}^{\ast}\rangle\\\
\approx&\int\mathrm{d}\vec{r}\int\mathrm{d}\vec{r^{\prime}}\rho_{n}^{i}(\vec{r})\left[\frac{1}{|\vec{r}-\vec{r^{\prime}}|}+g_{xc}(\vec{r},\vec{r^{\prime}})\right]\rho_{m}^{j}(\vec{r^{\prime}})-\omega_{0}\int\mathrm{d}\vec{r}\rho_{n}^{j}(\vec{r})\rho_{m}^{j}(\vec{r}),\end{split}$
(44)

where $\rho_{n}^{j}(\vec{r})$ is the transition density of the jth excited
state of the molecule m, $\omega_{0}$ is the transition frequency, and
$g_{xc}(\vec{r},\vec{r^{\prime}})$ is the exchange-correlation potential. All
other calculations, like the first excited states, are based on the long-range
corrected CAM-B3LYP functional and $6-31G^{\ast}$ basis set from Gaussian09
package.

Figure 6: The frequency dependence of mode-specific reorganization energy for
PBDTTPD (frequencies between 3000$cm^{-1}$ and 3300$cm^{-1}$ are not shown
because the peaks are too weak to display). The molecular structure of PBDTTPD
and it’s component unit is embedded. Reprinted with permission from Ref. [13].
Copyright 2015 American Chemical Society.

The mode-specific reorganization energies of PBDTTPD is shown in Fig. 6. In
the quantum dynamics simulation, an ultrafast laser pulse with a duration of
100 fs induces an initial exciton wavepacket with excess energy about 1.056 eV
higher than the bottom of the exciton band. At room temperature (298 K), after
1ps evolution we obtain the two time scales of the energy relaxation, 36 fs
and 491 fs, respectively. The contours of the population and coherence length
sequence evolution are shown in Fig. 7, from where one can find that the fast
time (36 fs) is consistent with the dephasing process, while the slow one (491
fs) corresponds to the hopping motion, although partial coherence is still
maintained. Using the femtosecond-resolved fluorescence up-conversion
technique coupled with global analysis, Banerji and coworkers[170, 176, 177]
have shown that more than 90% of the Stoke shift takes place within the
instrumental time resolution (200 fs), and a characteristic time of about 500
fs was attributed to a single adjacent hopping (0.5-1 ps). Our simulation
results are consistent with the experiments, thus we suggest that the fast
time scale observed experimentally is essentially the dephasing time of
exciton, whereas the slow one is corresponding to hopping motion.

Figure 7: (a) The time-evolving exciton population of PBDTTPD at room
temperature. (b) The corresponding contour map of the coherence length
sequence $\\{L_{k}(t)\\}$. Reprinted with permission from Ref. [13]. Copyright
2015 American Chemical Society.

#### Quantum Interferences in SF

Much attention has been given to the SF process recently[186, 187, 175, 188,
189, 190, 191, 192, 193, 194, 195, 196, 197] for the potential of overcoming
the Shockley-Queisser limit in single-junction photovoltaic devices[198, 199,
200]. Several works have found that the quantum interference between different
pathways of the SF process, such as the direct and indirect CT-mediated
pathways, can have strong impact on SF rates[201, 202, 203, 204].
Interestingly, there is an analogous process in aggregate spectra[205, 206,
207, 208], where a simple dimer model are denoted as H (J)-aggregates when the
direction of transition dipole moments of monomer are perpendicular (parallel)
to the line connecting the the centers of monomer, with the higher (lower)
excitonic states optically bright as revealed in the molecular spectra. Note
that H-aggregates usually have a positive exciton-exciton coupling strength,
while the J-aggregates have a negative one. Very recently, Zang and
coworkers[106] have investigated the quantum interferences among multiple
pathways in SF dynamics and found that they are closely related to the
properties of J- and H-aggregates.

We begin with a heterodimer model AB consisting of five electronic states
denoted as $A^{\ast}B$ ($AB^{\ast}$), $A^{+}B^{-}$ ($A^{-}B^{+}$) and
$A^{T}B^{T}$, corresponding to two localized singlet excited (SE) states, two
CT states and a spin-correlated TT state respectively, where $A^{*}(B^{*})$,
$A^{+}(B^{+})$, and $A^{-}(B^{-})$ represent the excited, positively charged,
and negatively charged states of chromophore $A$ ($B$), respectively. For
convenience, in the following we will refer to them as $S_{1}S_{0}$
($S_{0}S_{1}$), CA (AC) and TT states, respectively. The electronic
Hamiltonian can be expressed as

$\hat{H}_{e}=\left(\begin{array}[]{ccccc}E_{S_{1}S_{0}}&V_{ex}&V_{LL}&V_{HH}&0\\\
V_{ex}&E_{S_{0}S_{1}}&V_{HH}&V_{LL}&0\\\ V_{LL}&V_{HH}&E_{CA}&0&V_{LH}\\\
V_{HH}&V_{LL}&0&E_{AC}&V_{HL}\\\ 0&0&V_{LH}&V_{HL}&E_{TT}\\\
\end{array}\right).$ (45)

Fig. 8 presents the relationship and the couplings among those electronic
states.

Figure 8: A schematic diagram for the five electronic states and the couplings
among them. The signs of couplings marked with red are alterable, whereas the
couplings marked with black are fixed to be positive. Reprinted with
permission from Ref. [106]. Copyright 2017 American Chemical Society.

By changing the signs of the five couplings $V_{ex}$, $V_{LL}$, $V_{HH}$,
$V_{LH}$ and $V_{HL}$, various different quantum interference effects may
appear. However, it has been shown[106] that the results can be sorted into
several different groups. Under the assumption that the energies of the two CT
states are much higher than other states, the Hamiltonian Eq. (45) can be
simplified as an effective three-levels Hamiltonian via a partitioning
technique[209, 210]:

$\left(\begin{array}[]{ccc}\widetilde{E}_{S_{1}S_{0}}&\widetilde{J}&\widetilde{V}_{1}\\\
\widetilde{J}&\widetilde{E}_{S_{0}S_{1}}&\widetilde{V}_{2}\\\
\widetilde{V}_{1}&\widetilde{V}_{2}&\widetilde{E}_{TT}\\\ \end{array}\right),$
(46)

where the effective energies are

$\displaystyle\widetilde{E}_{S_{1}S_{0}}=E_{S_{1}S_{0}}+\frac{V^{2}_{LL}}{E_{S_{1}S_{0}}-E_{CA}}+\frac{V^{2}_{HH}}{E_{S_{1}S_{0}}-E_{AC}},$
$\displaystyle\widetilde{E}_{S_{0}S_{1}}=E_{S_{0}S_{1}}+\frac{V^{2}_{HH}}{E_{S_{0}S_{1}}-E_{CA}}+\frac{V^{2}_{LL}}{E_{S_{0}S_{1}}-E_{AC}},$
$\displaystyle\widetilde{E}_{TT}=E_{TT}+\frac{V^{2}_{LH}}{E_{TT}-E_{CA}}+\frac{V^{2}_{HL}}{E_{TT}-E_{AC}},$
(47)

and the effective couplings are

$\displaystyle\widetilde{J}=V_{ex}+\frac{2V_{LL}V_{HH}}{E_{S_{1}S_{0}}+E_{S_{0}S_{1}}-2E_{CA}}+\frac{2V_{HH}V_{LL}}{E_{S_{1}S_{0}}+E_{S_{0}S_{1}}-2E_{AC}},$
$\displaystyle\widetilde{V}_{1}=\frac{2V_{LL}V_{LH}}{E_{S_{1}S_{0}}+E_{TT}-2E_{CA}}+\frac{2V_{HH}V_{HL}}{E_{S_{1}S_{0}}+E_{TT}-2E_{AC}},$
$\displaystyle\widetilde{V}_{2}=\frac{2V_{HH}V_{LH}}{E_{S_{0}S_{1}}+E_{TT}-2E_{CA}}+\frac{2V_{LL}V_{HL}}{E_{S_{0}S_{1}}+E_{TT}-2E_{AC}}.$
(48)

$\widetilde{V}_{1}$ ($\widetilde{V}_{2}$) represents the effective coupling
between $S_{1}S_{0}$ ($S_{0}S_{1}$) state and TT state. The expression of
$\widetilde{V}_{1}$ includes two terms corresponding to two pathways
$S_{1}S_{0}$$\to$CA$\to$TT and $S_{1}S_{0}$$\to$AC$\to$TT respectively, and
when terms $V_{LL}V_{LH}$ and $V_{HH}V_{HL}$ have the same (opposite) signs,
the interference between the two pathways is constructive (destructive).

In the following, we will investigate what effects the effective exciton-
exciton coupling $\widetilde{J}$ will bring in. In this situation, an
interesting analogy between the SF dynamics and the aggregate emission
process[205, 206, 207] can be drawn if we regard the TT state as the ground
state and regard $\widetilde{V}_{1}$ and $\widetilde{V}_{1}$ as the transition
dipoles of the two monomers. Therefore, to explain the interference effects in
SF dynamics, we can directly borrow the concepts of J- and H-type aggregates
by partially diagonalize the effective three-levels Hamiltonian Eq. (45) (see
the upper panel of Fig. 9 for the schematic process):

$\left(\begin{array}[]{ccc}\widetilde{E}_{S_{1}S_{0}}&\widetilde{J}&\widetilde{V}_{1}\\\
\widetilde{J}&\widetilde{E}_{S_{0}S_{1}}&\widetilde{V}_{2}\\\
\widetilde{V}_{1}&\widetilde{V}_{2}&\widetilde{E}_{TT}\\\
\end{array}\right)\to\left(\begin{array}[]{ccc}E_{S_{+}}&0&V_{+}\\\
0&E_{S_{-}}&V_{-}\\\ V_{+}&V_{-}&\widetilde{E}_{TT}\\\ \end{array}\right),$
(49)

where
$E_{S_{+}}=\widetilde{E}_{S_{1}S_{0}}\mathrm{sin}^{2}\theta+\widetilde{E}_{S_{0}S_{1}}\mathrm{cos}^{2}\theta+\widetilde{J}\mathrm{sin}2\theta$,
$E_{S_{-}}=\widetilde{E}_{S_{1}S_{0}}\mathrm{cos}^{2}\theta+\widetilde{E}_{S_{0}S_{1}}\mathrm{sin}^{2}\theta-\widetilde{J}\mathrm{sin}2\theta$,
and the two new introduced couplings are
$V_{+}=\widetilde{V}_{1}\mathrm{sin}\theta+\widetilde{V}_{2}\mathrm{cos}\theta$
and
$V_{-}=\widetilde{V}_{1}\mathrm{cos}\theta-\widetilde{V}_{2}\mathrm{sin}\theta$,
with
$\mathrm{tan}2\theta=\frac{2\widetilde{J}}{\widetilde{E}_{S_{0}S_{1}}-\widetilde{E}_{S_{1}S_{0}}}$.

Figure 9: Schematic SF processes in J- and H-aggregates in terms of Eq. (49).
Reprinted with permission from Ref. [106]. Copyright 2017 American Chemical
Society.

It is noted that through artificial design strategies such as changing the
symmetry of covalent dimer[203], breaking local crystal structures[211],
controlling side chains or elements[212], one can effectively tune the
couplings $\widetilde{V}_{1}$ and $\widetilde{V}_{2}$. The consequence of J-
and H-aggregates in aggregate spectra can be further borrowed to guide the
tuning of $\widetilde{V}_{1}$ and $\widetilde{V}_{2}$. Under special
assumptions that $\widetilde{E}_{S_{1}S_{0}}=\widetilde{E}_{S_{0}S_{1}}=E$,
one obtains $E_{S_{\pm}}=E\pm\widetilde{J}$ and
$|V_{\pm}|=\frac{1}{\sqrt{2}}|\widetilde{V}_{1}\pm\widetilde{V}_{2}|$. If the
absolute value of $\widetilde{V}_{1}$ and $\widetilde{V}_{2}$ are further
assumed to be identical, one can conclude that no matter the type of the
aggregates, in the cases that
$\widetilde{J}\widetilde{V}_{1}\widetilde{V}_{2}<0$, the low-energy pathway is
favored and the SF rate is enhanced, while in the cases that
$\widetilde{J}\widetilde{V}_{1}\widetilde{V}_{2}>0$, only the high-energy
pathway is opened and the SF rate is suppressed. A schematic diagram for the
above analysis is shown in the lower panel of Fig. 9.

Figure 10: Population evolution of SF with (a) $V_{ex}=20$ meV, (b)
$V_{ex}=-40$ meV and (c) $V_{ex}=90$ meV, as well as TT state population from
only one SE state pathway (dashed line). (d) The SF rates in terms of $V_{ex}$
in dimer (solid line) and ten monomer (dashed line) models. Reprinted with
permission from Ref. [106]. Copyright 2017 American Chemical Society.

This effective three-levels picture is only valid when the energies of CT
states are much higher than the energies of other electronic states. Beyond
this regime, it is not easy to find the interference effects analytically.
Besides, The complex electron-phonon interaction should also be incorporated
in the analysis. In the following the perturbative non-Markovian stochastic
Schrödinger equation is adopted to investigate the interference effects behind
J- and H-aggregates, where the Debye spectral density Eq. (25) is adopted for
all the five electronic states with an identical set of parameters,
$\lambda=50$ meV and $\omega_{c}=1450$ cm-1[213], and the temperature being
300K. The electronic Hamiltonian elements are set as as
$E_{S_{1}S_{0}}-E_{TT}=200$ meV, $V_{LL}=V_{LH}=V_{HL}=V_{HH}=50$ meV, and
$E_{CA}(E_{AC})-E_{TT}=400$ meV[214, 201, 215, 216], whereas $V_{ex}$ is a
variable that controlling the type of the aggregates. The initial wavefunction
of the exciton is set as $(|S_{1}S_{0}\rangle+|S_{0}S_{1}\rangle)/\sqrt{2}$.
Fig. 10(a)-(c) show the population evolution with $V_{ex}=20,-40,90$ meV,
corresponding to the cases when $\widetilde{J}=0$, $<0$, and $>0$
respectively. It is clear that when $\widetilde{J}=0$, the interference
between two SE states disappears[208, 217], therefore the SF process should be
similar to that of the null-aggregate case, as can be seen in Fig. 10(a). It
is also clear in Fig. 10(b)-(c) that the negative (positive) $\widetilde{J}$
speeds up (suppresses) the SF processes as expected. Via an exponential
fitting procedure $P_{SE}(t)=e^{-kt}$, we can obtain the SF rate $k$, and the
results in terms of $V_{ex}$ are shown in Fig. 10(d). Focusing on the case of
$\widetilde{V}_{1}=\widetilde{V}_{2}$, it is clear that the SF rates when
$V_{ex}<20$ meV are larger than those when $V_{ex}>20$ meV, totally consistent
with the explanation borrowed from J- and H-aggregates. However, in the region
of $V_{ex}>20$ meV, an interesting behavior is that the SF rates first
increase and then decrease with the increase of $V_{ex}$. In this region, the
pathway from the lower-energy state $S_{-}$ (high-efficient one) to the TT
state is prohibited, while that from the higher-energy state $S_{+}$ (low-
efficient one) is allowed. When $\widetilde{J}$ is small, the energy splitting
between $S_{-}$ and $S_{+}$ states is not so high that the low-efficient
pathway may still enhance the SF rates as compared to the null aggregate case,
since $\widetilde{V}_{+}$ is larger than $\widetilde{V}_{1}$ and
$\widetilde{V}_{2}$. When $\widetilde{J}$ is high enough, the transition from
$S_{+}$ state to TT state becomes inefficient due to the large energy gap,
therefore the SF rates decrease.

In the above analysis, the J-aggregates with negative $\widetilde{J}$ allow
both the non-radiative SF and radiative emission pathways, which is not
favored from a practical view since the two processes may have similar time
scales[218]. One can suppress one of the pathways by controlling the values of
$\widetilde{V}_{1}$ and $\widetilde{V}_{2}$. For example, let $V_{LL}\to-
V_{LL}$ and $V_{LH}\to-V_{LH}$, we can get
$\widetilde{V}_{1}=-\widetilde{V}_{2}$ (see Eq. (Quantum Interferences in
SF)), and the roles of J- and H-aggregates are totally inversed, as is seen
from the purple line in Fig. 10(d). To show that these results are general,
the same calculations are also done for a one-dimensional chain model with 10
monomers[104], and the same relationship between SF rates and $V_{ex}$ is
obtained as shown in Fig. 10(d). The only discrepancy is that the SF rates are
larger than those of the dimer model, which is due to the exciton
delocalization effect[219, 104].

## CONCLUSIONS

Based on the path integral formalism of an open system coupled to harmonic
baths, we have outlined a group of stochastic wavefunction-based approaches,
i.e., non-Markovian stochastic Schrödinger equations, which are presented in a
progressive manner in terms of the level of approximations so as to set up an
integrated platform for its versatile applications in various physical and
chemical scenarios. First of all, we have introduced a mixed strategy
combining the stochastic processes and the hierarchical techniques to obtain a
hierarchical set of stochastic differential equations, which are numerically
exact and suitable for the cases where at least part of the bath correlation
functions exhibit exponential decaying features. We have further demonstrated
that within the Debye bath, by choosing the residual bath correlation function
to be a pure imaginary function, one can obtain a specific set of equations
which require the smallest number of the auxiliary wavefunctions, effectively
avoiding the perplexity along with the low temperature conditions in other
hierarchical quantum master equations for density matrix. Then, considering
the cases where the system-bath coupling is not quite strong, the bath is
highly structured, or the system size is very large, we further summarized a
systematic perturbation expansion scheme with respect to the system-bath
coupling. Under the lowest order approximation, the perturbative stochastic
Schrödinger equations are found to be reliable in a very broad parameter
range, which is highly commendable since it arrives at an optimal balance
between accuracy and computational cost. Of course, the equations can be
further simplified if the systems are highly-symmetric or the bath memory is
extremely short-lived. In the last part, a few applications of the above-
mentioned approaches are presented, including the calculations of excitation
energy transfer in the FMO complex, quantum transport properties of hot
exciton in a copolymer chain and the studies about the quantum interference
effects on the the promising SF processes in organic materials. All in all,
the non-Markovian stochastic Schrödinger equations due to its distinctive
advantages, are expected to play a crucial role in the nanoscale realistic
systems in the near future.

## ACKNOWLEDGEMENTS

This work is supported by the National Science Foundation of China (Grant No.
21573175 and 21773191).

## References

  * [1] Tobias Brixner, Jens Stenger, Harsha M Vaswani, Minhaeng Cho, Robert E Blankenship, and Graham R Fleming.  Two-dimensional spectroscopy of electronic couplings in photosynthesis.  Nature, 434(7033):625–628, 2005.
  * [2] Gregory S. Engel, Tessa R. Calhoun, Elizabeth L. Read, Tae-Kyu Ahn, Tomáš Mančal, Yuan-Chung Cheng, Robert E. Blankenship, and Graham R. Fleming.  Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems.  Nature, 446(7137):782–786, 2007.
  * [3] Hohjai Lee, Yuan-Chung Cheng, and Graham R Fleming.  Coherence dynamics in photosynthesis: protein protection of excitonic coherence.  Science, 316(5830):1462–1465, 2007.
  * [4] Elisabetta Collini, Cathy Y Wong, Krystyna E Wilk, Paul MG Curmi, Paul Brumer, and Gregory D Scholes.  Coherently wired light-harvesting in photosynthetic marine algae at ambient temperature.  Nature, 463(7281):644–647, 2010.
  * [5] Gitt Panitchayangkoon, Dugan Hayes, Kelly A Fransted, Justin R Caram, Elad Harel, Jianzhong Wen, Robert E Blankenship, and Gregory S Engel.  Long-lived quantum coherence in photosynthetic complexes at physiological temperature.  Proc. Natl. Acad. Sci. USA, 107(29):12766–12770, 2010.
  * [6] Gregory D Scholes, Graham R Fleming, Lin X Chen, Alán Aspuru-Guzik, Andreas Buchleitner, David F Coker, Gregory S Engel, Rienk van Grondelle, Akihito Ishizaki, David M Jonas, et al.  Using coherence to enhance function in chemical and biophysical systems.  Nature, 543(7647):647–656, 2017.
  * [7] Carles Curutchet and Benedetta Mennucci.  Quantum chemical studies of light harvesting.  Chem. Rev., 2016.
  * [8] Aurélia Chenu and Gregory D Scholes.  Coherence in energy transfer and photosynthesis.  Annu. Rev. Phys. Chem., 66:69–96, 2015.
  * [9] Natalie Banerji, Sarah Cowan, Eric Vauthey, and Alan J Heeger.  Ultrafast relaxation of the poly (3-hexylthiophene) emission spectrum.  J. Phys. Chem. C, 115(19):9726–9739, 2011.
  * [10] Nicholas R Monahan, Kristopher W Williams, Bharat Kumar, Colin Nuckolls, and X-Y Zhu.  Direct observation of entropy-driven electron-hole pair separation at an organic semiconductor interface.  Phys. Rev. Lett., 114(24):247003, 2015.
  * [11] Simon Gélinas, Akshay Rao, Abhishek Kumar, Samuel L Smith, Alex W Chin, Jenny Clark, Tom S van der Poll, Guillermo C Bazan, and Richard H Friend.  Ultrafast long-range charge separation in organic semiconductor photovoltaic diodes.  Science, 343(6170):512–516, 2014.
  * [12] Loren G Kaake, Daniel Moses, and Alan J Heeger.  Coherence and uncertainty in nanostructured organic photovoltaics.  J. Phys. Chem. Lett., 4(14):2264–2268, 2013.
  * [13] Yaling Ke, Yuxiu Liu, and Yi Zhao.  Visualization of hot exciton energy relaxation from coherent to diffusive regimes in conjugated polymers: A theoretical analysis.  J. Phys. Chem. Lett., 6(9):1741–1747, 2015.
  * [14] Jean-Luc Brédas, Edward H Sargent, and Gregory D Scholes.  Photovoltaic concepts inspired by coherence effects in photosynthetic systems.  Nat. Mater., 16(1):35–44, 2017.
  * [15] Gordon J Hedley, Arvydas Ruseckas, and Ifor DW Samuel.  Light harvesting for organic photovoltaics.  Chem. rev., 117(2):796–837, 2016.
  * [16] Ralf Bulla, Ning-Hua Tong, and Matthias Vojta.  Numerical renormalization group for bosonic systems and application to the sub-ohmic spin-boson model.  Phys. Rev. Lett., 91(17):170601, 2003.
  * [17] Ralf Bulla, Hyun-Jung Lee, Ning-Hua Tong, and Matthias Vojta.  Numerical renormalization group for quantum impurities in a bosonic bath.  Phys. Rev. B, 71(4):045122, 2005.
  * [18] Steven R White and Adrian E Feiguin.  Real-time evolution using the density matrix renormalization group.  Phys. Rev. Lett., 93(7):076401, 2004.
  * [19] Alex W Chin, Ángel Rivas, Susana F Huelga, and Martin B Plenio.  Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials.  J. Math. Phys., 51(9):092109, 2010.
  * [20] Haobin Wang and Michael Thoss.  Multilayer formulation of the multiconfiguration time-dependent hartree theory.  J. Chem. Phys., 119(3):1289–1299, 2003.
  * [21] Mathias Duckheim and Joachim Ankerhold.  Influence of low-frequency noise on macroscopic quantum tunneling in superconducting circuits.  Phys. Rev. B, 71(13):134501, 2005.
  * [22] Frithjof B. Anders, Ralf Bulla, and Matthias Vojta.  Equilibrium and nonequilibrium dynamics of the sub-ohmic spin-boson model.  Phys. Rev. Lett., 98(21), 2007.
  * [23] André Winter, Heiko Rieger, Matthias Vojta, and Ralf Bulla.  Quantum phase transition in the sub-ohmic spin-boson model: Quantum monte carlo study with a continuous imaginary time cluster algorithm.  Phys. Rev. Lett., 102(3):030601, 2009.
  * [24] Haobin Wang and Michael Thoss.  From coherent motion to localization: Ii. dynamics of the spin-boson model with sub-ohmic spectral density at zero temperature.  Chem. Phys., 370(1):78–86, 2010.
  * [25] Denis Kast and Joachim Ankerhold.  Bipartite entanglement dynamics of two-level systems in sub-ohmic reservoirs.  Phys. Rev. B, 90(10):100301, 2014.
  * [26] Denis Kast and Joachim Ankerhold.  Persistence of coherent quantum dynamics at strong dissipation.  Phys. Rev. Lett., 110(1):010402, 2013.
  * [27] Richard Phillips Feynman and Frank Lee Vernon.  The theory of a general quantum system interacting with a linear dissipative system.  Ann. Phys., 24:118–173, 1963.
  * [28] MF Herman, EJ Bruskin, and BJ Berne.  On path integral monte carlo simulations.  J. Chem. Phys., 76(10):5150–5155, 1982.
  * [29] R Egger and U Weiss.  Quantum monte carlo simulation of the dynamics of the spin-boson model.  Zeitschrift für Physik B Condensed Matter, 89(1):97–107, 1992\.
  * [30] Reinhold Egger and C. H. Mak.  Low-temperature dynamical simulation of spin-boson systems.  Phys. Rev. B, 50(20):15210–15220, 1994.
  * [31] Nancy Makri.  Improved feynman propagators on a grid and non-adiabatic corrections within the path integral framework.  Chem. Phys. Lett., 193(5):435–445, 1992.
  * [32] Maria Topaler and Nancy Makri.  Multidimensional path integral calculations with quasiadiabatic propagators: Quantum dynamics of vibrational relaxation in linear hydrocarbon chains.  J. Chem. Phys, 97(12):9001–9015, 1992.
  * [33] Maria Topaler and Nancy Makri.  Quasi-adiabatic propagator path integral methods. exact quantum rate constants for condensed phase reactions.  Chem. Phys. Lett., 210(1):285–293, 1993.
  * [34] Maria Topaler and Nancy Makri.  System-specific discrete variable representations for path integral calculations with quasi-adiabatic propagators.  Chem. Phys. Lett., 210(4-6):448–457, 1993.
  * [35] Dmitrii E Makarov and Nancy Makri.  Path integrals for dissipative systems by tensor multiplication. condensed phase quantum dynamics for arbitrarily long time.  Chem. Phys. Lett., 221(5):482–491, 1994.
  * [36] Nancy Makri and Dmitrii E Makarov.  Tensor propagator for iterative quantum time evolution of reduced density matrices. i. theory.  J. Chem. Phys., 102(11):4600–4610, 1995.
  * [37] Nancy Makri and Dmitrii E Makarov.  Tensor propagator for iterative quantum time evolution of reduced density matrices. ii. numerical methodology.  J. Chem. Phys., 102(11):4611–4618, 1995.
  * [38] Jiushu Shao and Nancy Makri.  Iterative path integral calculation of quantum correlation functions for dissipative systems.  Chem. Phys., 268(1-3):1–10, 2001.
  * [39] Jiushu Shao and Nancy Makri.  Iterative path integral formulation of equilibrium correlation functions for quantum dissipative systems.  J. Chem. Phys., 116(2):507–514, 2002.
  * [40] Nancy Makri.  Numerical path integral techniques for long time dynamics of quantum dissipative systems.  J. Math. Phys., 36(5):2430–2457, 1995.
  * [41] Eunji Sim and Nancy Makri.  Filtered propagator functional for iterative dynamics of quantum dissipative systems.  Comput. Phys. Commun., 99(2):335–354, 1997.
  * [42] Yoshitaka Tanimura and Ryogo 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.
  * [43] Yun-An Yan, Fan Yang, Yu Liu, and Jiushu Shao.  Hierarchical approach based on stochastic decoupling to dissipative systems.  Chem. Phys. Lett., 395(4):216–221, 2004.
  * [44] Rui-Xue Xu, Ping Cui, Xin-Qi Li, Yan Mo, and YiJing Yan.  Exact quantum master equation via the calculus on path integrals.  J. Chem. Phys., 122(4):041103, 2005.
  * [45] Akihito Ishizaki and Yoshitaka 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.
  * [46] Jie Hu, Rui-Xue Xu, and Yijing Yan.  Communication: Padé spectrum decomposition of fermi function and bose function.  J. Chem. Phys., 133(10):101106, 2010.
  * [47] Yoshitaka Tanimura.  Stochastic liouville, langevin, fokker–planck, and master equation approaches to quantum dissipative systems.  J. Phys. Soc. Jpn., 75(8):082001, 2006.
  * [48] Markus Schröder, Michael Schreiber, and Ulrich Kleinekathöfer.  Reduced dynamics of coupled harmonic and anharmonic oscillators using higher-order perturbation theory.  J. Chem. Phys., 126(11):114102, 2007.
  * [49] Christoph Meier and David J Tannor.  Non-markovian evolution of the density operator in the presence of strong laser fields.  J. Chem. Phys., 111(8):3365–3376, 1999.
  * [50] Zhoufei Tang, Xiaolong Ouyang, Zhihao Gong, Haobin Wang, and Jianlan Wu.  Extended hierarchy equation of motion for the spin-boson model.  J. Chem. Phys., 143(22):224112, 2015.
  * [51] Chenru Duan, Zhoufei Tang, Jianshu Cao, and Jianlan Wu.  Zero-temperature localization in a sub-ohmic spin-boson model investigated by an extended hierarchy equation of motion.  Phys. Rev. B, 95(21):214308, 2017.
  * [52] Qiang Shi, Liping Chen, Guangjun Nan, Rui-Xue Xu, and YiJing Yan.  Efficient hierarchical liouville space propagator to quantum dissipative dynamics.  J. Chem. Phys., 130(8):084105, 2009.
  * [53] Hao Liu, Lili Zhu, Shuming Bai, and Qiang 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):134106, 2014.
  * [54] David M Wilkins and Nikesh S Dattani.  Why quantum coherence is not important in the fenna–matthews–olsen complex.  J. Chem. Theory Comput., 11(7):3411–3419, 2015.
  * [55] Christoph Kreisbeck, Tobias Kramer, Mirta Rodriguez, and Birgit Hein.  High-performance solution of hierarchical equations of motion for studying energy transfer in light-harvesting complexes.  J. Chem. Theory Comput., 7(7):2166–2174, 2011.
  * [56] Christoph Kreisbeck, Tobias Kramer, and Alan Aspuru-Guzik.  Disentangling electronic and vibronic coherences in two-dimensional echo spectra.  J. Phys. Chem. B, 117(32):9380–9385, 2013.
  * [57] Christoph Kreisbeck and Tobias Kramer.  Long-lived electronic coherence in dissipative exciton dynamics of light-harvesting complexes.  J. Phys. Chem. Lett., 3(19):2828–2833, 2012.
  * [58] Masashi Tsuchimoto and Yoshitaka Tanimura.  Spins dynamics in a dissipative environment: Hierarchal equations of motion approach using a graphics processing unit (gpu).  J. Chem. Theory Comput., 11(8):3859–3865, 2015.
  * [59] Johan Strümpfer and Klaus Schulten.  Open quantum dynamics calculations with the hierarchy equations of motion on parallel computers.  J. Chem. Theory Comput., 8(8):2808–2816, 2012.
  * [60] Christoph Kreisbeck, Tobias Kramer, and Alán 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.
  * [61] Akihito Ishizaki and Graham R Fleming.  Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature.  Proc. Natl. Acad. Sci. U.S.A., 106(41):17255–17260, 2009.
  * [62] Liping Chen, Renhui Zheng, Qiang Shi, and YiJing Yan.  Optical line shapes of molecular aggregates: Hierarchical equations of motion method.  J. Chem. Phys., 131(9):094502, 2009.
  * [63] Liping Chen, Renhui Zheng, Yuanyuan Jing, and Qiang Shi.  Simulation of the two-dimensional electronic spectra of the fenna-matthews-olson complex using the hierarchical equations of motion method.  J. Chem. Phys., 134(19):194508, 2011.
  * [64] Jürgen T Stockburger and CH Mak.  Dynamical simulation of current fluctuations in a dissipative two-state system.  Phys. Rev. Lett., 80(12):2657, 1998.
  * [65] Jürgen T Stockburger and Chi H Mak.  Stochastic liouvillian algorithm to simulate dissipative quantum dynamics with arbitrary precision.  J. Chem. Phys., 110(11):4983–4985, 1999.
  * [66] Jürgen T Stockburger and Hermann Grabert.  Non-markovian quantum state diffusion.  Chem. Phys., 268(1):249–256, 2001.
  * [67] Jürgen T Stockburger and Hermann Grabert.  Exact c-number representation of non-markovian quantum dissipation.  Phys. Rev. Lett., 88(17):170407, 2002.
  * [68] Jiushu Shao.  Decoupling quantum dissipation interaction via stochastic fields.  J. Chem. Phys., 120(11):5053–5056, 2004.
  * [69] Y Zhou, Y Yan, and J Shao.  Stochastic simulation of quantum dissipative dynamics.  Europhys. Lett., 72(3):334, 2005.
  * [70] Yun Zhou and Jiushu Shao.  Solving the spin-boson model of strong dissipation with flexible random-deterministic scheme.  J. Chem. Phys., 128(3):034106, 2008.
  * [71] Jeremy M Moix and Jianshu Cao.  A hybrid stochastic hierarchy equations of motion approach to treat the low temperature dynamics of non-markovian open quantum systems.  J. Chem. Phys., 139(13):134106, 2013.
  * [72] Lili Zhu, Hao Liu, and Qiang Shi.  A new method to account for the difference between classical and quantum baths in quantum dissipative dynamics.  New J. Phys., 15(9):095020, 2013.
  * [73] Inés de Vega and Daniel Alonso.  Dynamics of non-markovian open quantum systems.  Rev. Mod. Phys., 89(1), 2017.
  * [74] Masoud Mohseni, Yasser Omar, Gregory S Engel, and Martin B Plenio.  Quantum effects in biology.  Cambridge University Press, 2014.
  * [75] Volkhard May and Oliver Kühn.  Charge and energy transfer dynamics in molecular systems.  John Wiley & Sons, 2008.
  * [76] Heinz-Peter Breuer and Francesco Petruccione.  The theory of open quantum systems.  Oxford University Press on Demand, 2002.
  * [77] Alfred G Redfield.  On the theory of relaxation processes.  IBM J. Res. Dev., 1(1):19–31, 1957.
  * [78] AG Redfield.  The theory of relaxation processes.  In Advances in Magnetic and Optical Resonance, volume 1, pages 1–32. Elsevier, 1965.
  * [79] Seogjoo Jang, Yuan-Chung Cheng, David R Reichman, and Joel D Eaves.  Theory of coherent resonance energy transfer, 2008.
  * [80] Seogjoo Jang.  Theory of coherent resonance energy transfer for coherent initial condition.  J. Chem. Phys., 131(16):164101, 2009.
  * [81] Ahsan Nazir.  Correlation-dependent coherent to incoherent transitions in resonant energy transfer dynamics.  Phys. Rev. Lett., 103(14):146404, 2009.
  * [82] Seogjoo Jang.  Theory of multichromophoric coherent resonance energy transfer: A polaronic quantum master equation approach.  J. Chem. Phys., 135(3):034105, 2011.
  * [83] Avinash Kolli, Ahsan Nazir, and Alexandra Olaya-Castro.  Electronic excitation dynamics in multichromophoric systems described via a polaron-representation master equation.  J. Chem. Phys., 135(15):154112, 2011.
  * [84] Dara PS McCutcheon and Ahsan Nazir.  Coherent and incoherent dynamics in excitonic energy transfer: Correlated fluctuations and off-resonance effects.  Phys. Rev. B, 83(16):165101, 2011.
  * [85] Robert Silbey and Robert A Harris.  Variational calculation of the dynamics of a two level system interacting with a bath.  J. Chem. Phys., 80(6):2615–2617, 1984.
  * [86] Robert A Harris and Robert Silbey.  Variational calculation of the tunneling system interacting with a heat bath. ii. dynamics of an asymmetric tunneling system.  J. Chem. Phys., 83(3):1069–1074, 1985.
  * [87] Dara PS McCutcheon and Ahsan Nazir.  Consistent treatment of coherent and incoherent energy transfer dynamics using a variational master equation.  J. Chem. Phys., 135(11):114501, 2011.
  * [88] Dara PS McCutcheon, Nikesh S Dattani, Erik M Gauger, Brendon W Lovett, and Ahsan Nazir.  A general approach to quantum dynamics using a variational master equation: Application to phonon-damped rabi rotations in quantum dots.  Phys. Rev. B, 84(8):081305, 2011.
  * [89] Chee Kong Lee, Jeremy Moix, and Jianshu Cao.  Accuracy of second order perturbation theory in the polaron and variational polaron frames.  J. Chem. Phys., 136(20):204120, 2012.
  * [90] Felix A Pollock, Dara PS McCutcheon, Brendon W Lovett, Erik M Gauger, and Ahsan Nazir.  A multi-site variational master equation approach to dissipative energy transfer.  New J. Phys., 15(7):075018, 2013.
  * [91] Lajos Diósi and Walter T Strunz.  The non-markovian stochastic schrödinger equation for open systems.  Phys. Lett. A, 235(6):569–573, 1997.
  * [92] Lajos Diósi, Nicolas Gisin, and Walter T Strunz.  Non-markovian quantum state diffusion.  Phys. Rev. A, 58(3):1699, 1998.
  * [93] Walter T Strunz, Lajos Diósi, and Nicolas Gisin.  Open system dynamics with non-markovian quantum trajectories.  Phys. Rev. Lett., 82(9):1801, 1999.
  * [94] Ting Yu.  Non-markovian quantum trajectories versus master equations: Finite-temperature heat bath.  Phys. Rev. A, 69(6), 2004.
  * [95] Ting Yu, Lajos Diósi, Nicolas Gisin, and Walter T Strunz.  Non-markovian quantum-state diffusion: Perturbation approach.  Phys. Rev. A, 60(1):91, 1999.
  * [96] Gerhard Ritschel, Jan Roden, Walter T Strunz, and Alexander Eisfeld.  An efficient method to calculate excitation energy transfer in light-harvesting systems: application to the fenna–matthews–olson complex.  New J. Phys., 13(11):113034, 2011.
  * [97] Pierre Gaspard and Masataka Nagaoka.  Non-markovian stochastic schrödinger equation.  J. Chem. Phys., 111(13):5676–5690, 1999.
  * [98] Inés de Vega, Daniel Alonso, Pierre Gaspard, and Walter T Strunz.  Non-markovian stochastic schrödinger equations in different temperature regimes: A study of the spin-boson model.  J. Chem. Phys., 122(12):124106, 2005.
  * [99] Robert Biele and Roberto D’Agosta.  A stochastic approach to open quantum systems.  J. Phys.: Condens. Matter, 24(27):273201, 2012.
  * [100] Xinxin Zhong and Yi Zhao.  Non-markovian stochastic schrödinger equation at finite temperatures for charge carrier dynamics in organic crystals.  J. Chem. Phys., 138(1):014111, 2013.
  * [101] Xin Xin Zhong, Yi Zhao, and Jian Shu Cao.  Coherent quantum transport in disordered systems: Ii. temperature dependence of carrier diffusion coefficients from the time-dependent wavepacket diffusion method.  New J. Phys., 16(4):045009, 2014.
  * [102] Lu Han, Xin Xin Zhong, Wan Zhen Liang, and Yi Zhao.  Energy relaxation and separation of a hot electron-hole pair in organic aggregates from a time-dependent wavepacket diffusion method.  J. Chem. Phys., 140(21):214107, 2014.
  * [103] Takatoshi Fujita, Sule Atahan-Evrenk, Nicolas P. D. Sawaya, and Alán Aspuru-Guzik.  Coherent dynamics of mixed frenkel and charge-transfer excitons in dinaphtho[2,3-b:2$\prime$3$\prime$-f]thieno[3,2-b]-thiophene thin films: The importance of hole delocalization.  J. Phys. Chem. Lett., 7(7):1374–1380, apr 2016.
  * [104] Hang Zang, Yaling Ke, Yi Zhao, and WanZhen Liang.  Effects of charge transfer state and exciton migration on singlet fission dynamics in organic aggregates.  J. Phys. Chem. C, 120(25):13351–13359, jun 2016.
  * [105] Thomas Plehn and Volkhard May.  Charge and energy migration in molecular clusters: A stochastic schrödinger equation approach.  J. Chem. Phys., 146(3):034107, 2017.
  * [106] Hang Zang, Yi Zhao, and WanZhen Liang.  Quantum interference in singlet fission: J-and h-aggregate behavior.  J. Phys. Chem. Lett., 8(20):5105–5112, 2017.
  * [107] ZhiYe Zhu, Hang Zang, Yi Zhao, and WanZhen Liang.  Charge carrier mobilities and singlet fission dynamics in thienoquinoidal compounds.  J. Phys. Chem. C, 121(41):22587–22596, 2017.
  * [108] D Suess, A Eisfeld, and WT Strunz.  Hierarchy of stochastic pure states for open quantum system dynamics.  Phys. Rev. Lett., 113(15):150403, 2014.
  * [109] Daniel Suess, Walter T Strunz, and Alexander Eisfeld.  Hierarchical equations for open system dynamics in fermionic and bosonic environments.  J. Stat. Phys., 159(6):1408–1423, 2015.
  * [110] Gerhard Ritschel, Daniel Suess, Sebastian Möbius, Walter T Strunz, and Alexander Eisfeld.  Non-markovian quantum state diffusion for temperature-dependent linear spectra of light harvesting aggregates.  J. Chem. Phys., 142(3):034115, 2015.
  * [111] Kai Song, Linze Song, and Qiang Shi.  An alternative realization of the exact non-markovian stochastic schrödinger equation.  J. Chem. Phys., 144(22):224105, 2016.
  * [112] Yaling Ke and Yi Zhao.  Hierarchy of forward-backward stochastic schrödinger equation.  J. Chem. Phys., 145(2):024101, 2016.
  * [113] Yaling Ke and Yi Zhao.  Perturbation expansions of stochastic wavefunctions for open quantum systems.  J. Chem. Phys., 147(18):184103, 2017.
  * [114] J Hubbard.  Calculation of partition functions.  Phys. Rev. Lett., 3(2):77, 1959.
  * [115] RL Stratonovich.  On a method of calculating quantum distribution functions.  In Soviet Physics Doklady, volume 2, page 416, 1957.
  * [116] Heinz-Peter Breuer and Francesco Petruccione.  Destruction of quantum coherence through emission of bremsstrahlung.  Phys. Rev. A, 63(3), 2001.
  * [117] Akihito Ishizaki and Yoshitaka Tanimura.  Nonperturbative non-markovian quantum master equation: Validity and limitation to calculate nonlinear response functions.  Chem. Phys., 347(1):185–193, 2008.
  * [118] Masashi Ban, Sachiko Kitajima, and Fumiaki Shibata.  Reduced dynamics and the master equation of open quantum systems.  Phys. Lett. A, 374(23):2324–2330, 2010.
  * [119] Heinz-Peter Breuer and Francesco Petruccione.  The Theory of Open Quantum Systems.  Oxford University Press, 2002.
  * [120] Yaling Ke and Yi Zhao.  Hierarchy of stochastic schrödinger equation towards the calculation of absorption and circular dichroism spectra.  J. Chem. Phys., 146(17):174105, 2017.
  * [121] Yaling Ke and Yi Zhao.  Quantum dynamics simulations in an ultraslow bath using hierarchy of stochastic schrödinger equations.  Mol. Phys., pages 1–10, 2018.
  * [122] Marco del Rey, Alex W Chin, Susana F Huelga, and Martin B Plenio.  Exploiting structured environments for efficient energy transfer: the phonon antenna mechanism.  J. Phys. Chem. Lett., 4(6):903–907, 2013.
  * [123] Edward J O’Reilly and Alexandra Olaya-Castro.  Non-classicality of the molecular vibrations assisting exciton energy transfer at room temperature.  Nat. Commun., 5, 2014.
  * [124] Arend G Dijkstra, Chen Wang, Jianshu Cao, and Graham R Fleming.  Coherent exciton dynamics in the presence of underdamped vibrations.  J. Phys. Chem. Lett., 6(4):627–632, 2015.
  * [125] Jake Iles-Smith, Arend G Dijkstra, Neill Lambert, and Ahsan Nazir.  Energy transfer in structured and unstructured environments: Master equations beyond the born-markov approximations.  J. Chem. Phys, 144(4):044110, 2016.
  * [126] Patrick Rebentrost and Alán Aspuru-Guzik.  Communication: Exciton–phonon information flow in the energy transfer process of photosynthetic complexes.  J. Chem. Phys., 134(10):101103, 2011.
  * [127] C. A. Mujica-Martinez, P. Nalbach, and M. Thorwart.  Quantification of non-markovian effects in the fenna-matthews-olson complex.  Phys. Rev. E, 88(6), 2013.
  * [128] Felix Vaughan, Noah Linden, and Frederick R Manby.  How markovian is exciton dynamics in purple bacteria?  J. Chem. Phys., 146(12):124113, 2017.
  * [129] Niklas Christensson, Harald F Kauffmann, Tonu Pullerits, and Tomáš Mančal.  Origin of long-lived coherences in light-harvesting complexes.  J. Phys. Chem. B, 116(25):7449–7454, 2012.
  * [130] Vivek Tiwari, William K Peters, and David M Jonas.  Electronic resonance with anticorrelated pigment vibrations drives photosynthetic energy transfer outside the adiabatic framework.  Proc. Natl. Acad. Sci. U.S.A., 110(4):1203–1208, 2013.
  * [131] Aurélia Chenu, Niklas Christensson, Harald F Kauffmann, and Tomáš Mančal.  Enhancement of vibronic and ground-state vibrational coherences in 2d spectra of photosynthetic complexes.  Sci. Rep., 3:2029, 2013.
  * [132] AW Chin, J Prior, R Rosenbach, F Caycedo-Soler, SF Huelga, and MB Plenio.  The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment-protein complexes.  Nat. Phys., 9(2):113, 2013.
  * [133] F.K. Wilhelm, S. Kleff, and J. von Delft.  The spin-boson model with a structured environment: a comparison of approaches.  Chem. Phys., 296(2-3):345–353, 2004.
  * [134] M Thorwart, E Paladino, and M Grifoni.  Dynamics of the spin-boson model with a structured environment.  Chem. Phys., 296(2-3):333–344, 2004.
  * [135] Xin Xin Zhong and Yi Zhao.  Charge carrier dynamics in phonon-induced fluctuation systems from time-dependent wavepacket diffusion approach.  J. Chem. Phys., 135(13):134110, 2011.
  * [136] H Haken and P Reineker.  The coupled coherent and incoherent motion of excitons and its influence on the line shape of optical absorption.  Z. Phys., 249(3):253–268, 1972.
  * [137] H Haken and G Strobl.  An exactly solvable model for coherent and incoherent exciton motion.  Z. Phys., 262(2):135–148, 1973.
  * [138] P Reineker, VM Kenkre, and R Kühne.  Drift mobility of photo-electrons in organic molecular crystals: Quantitative comparison between theory and experiment.  Phys. Lett. A, 84(5):294–296, 1981.
  * [139] WeiWei Zhang, XinXin Zhong, and Yi Zhao.  Electron mobilities of n-type organic semiconductors from time-dependent wavepacket diffusion method: pentacenequinone derivatives.  J. Phys. Chem. A, 116(46):11075–11082, 2012.
  * [140] Yuqian Jiang, Xinxin Zhong, Wen Shi, Qian Peng, Hua Geng, Yi Zhao, and Zhigang Shuai.  Nuclear quantum tunnelling and carrier delocalization effects to bridge the gap between hopping and bandlike behaviors in organic semiconductors.  Nanoscale Horizons, 1(1):53–59, 2016.
  * [141] Chaoren Liu, David N Beratan, and Peng Zhang.  Coarse-grained theory of biological charge transfer with spatially and temporally correlated noise.  J. Phys. Chem. B, 120(15):3624–3633, 2016.
  * [142] Nicolas P. D. Sawaya, Joonsuk Huh, Takatoshi Fujita, Semion K. Saikin, and Alán Aspuru-Guzik.  Fast delocalization leads to robust long-range excitonic transfer in a large quantum chlorosome model.  Nano Lett., 15(3):1722–1729, mar 2015.
  * [143] Étienne Boulais, Nicolas PD Sawaya, Rémi Veneziano, Alessio Andreoni, James L Banal, Toru Kondo, Sarthak Mandal, Su Lin, Gabriela S Schlau-Cohen, Neal W Woodbury, et al.  Programmed coherent coupling in a synthetic dna-based excitonic circuit.  Nat. Mater., 17(2):159, 2018.
  * [144] Minhaeng Cho, Harsha M Vaswani, Tobias Brixner, Jens Stenger, and Graham R Fleming.  Exciton analysis in 2d electronic spectroscopy.  J. Phys. Chem. B, 109:10542–10556, 2005.
  * [145] Birgit Hein, Christoph Kreisbeck, Tobias Kramer, and Mirta Rodríguez.  Modelling of oscillations in two-dimensional echo-spectra of the fenna–matthews–olson complex.  New J. Phys., 14(2):023018, 2012.
  * [146] Hong-Guang Duan, Valentyn I Prokhorenko, Richard J Cogdell, Khuram Ashraf, Amy L Stevens, Michael Thorwart, and RJ Dwayne Miller.  Nature does not rely on long-lived electronic quantum coherence for photosynthetic energy transfer.  Proc. Natl. Acad. Sci. U.S.A., 114(32):8493–8498, 2017.
  * [147] Carsten Olbrich, Thomas LC Jansen, Jörg Liebers, Mortaza Aghtar, Johan Strümpfer, Klaus Schulten, Jasper Knoester, and Ulrich Kleinekathöfer.  From atomistic modeling to excitation transfer and two-dimensional spectra of the fmo light-harvesting complex.  J. Phys. Chem. B, 115(26):8609–8621, 2011.
  * [148] Julia Adolphs and Thomas Renger.  How proteins trigger excitation energy transfer in the fmo complex of green sulfur bacteria.  Biophys. J., 91(8):2778–2797, 2006.
  * [149] Junkuo Gao, Wu-Jun Shi, Jun Ye, Xiaoqing Wang, Hajime Hirao, and Yang Zhao.  Qm/mm modeling of environmental effects on electronic transitions of the fmo complex.  J. Phys. Chem. B, 117(13):3488–3495, 2013.
  * [150] Jeremy Moix, Jianlan Wu, Pengfei Huo, David Coker, and Jianshu Cao.  Efficient energy transfer in light-harvesting systems, iii: The influence of the eighth bacteriochlorophyll on the dynamics and efficiency in fmo.  J. Phys. Chem. Lett., 2(24):3045–3052, 2011.
  * [151] P Nalbach, CA Mujica-Martinez, and M Thorwart.  Vibronically coherent speed-up of the excitation energy transfer in the fenna-matthews-olson complex.  Phys. Rev. E, 91(2):022706, 2015.
  * [152] Jan Schulze, Mohamed F Shibl, Mohammed J Al-Marri, and Oliver Kühn.  Multi-layer multi-configuration time-dependent hartree (ml-mctdh) approach to the correlated exciton-vibrational dynamics in the fmo complex.  J. Chem. Phys., 144(18):185101, 2016.
  * [153] Neill Lambert, Yueh-Nan Chen, Yuan-Chung Cheng, Che-Ming Li, Guang-Yin Chen, and Franco Nori.  Quantum biology.  Nat. Phys., 9(1):10–18, 2013.
  * [154] Gregory D Scholes, Graham R Fleming, Alexandra Olaya-Castro, and Rienk van Grondelle.  Lessons from nature about solar light harvesting.  Nat. Chem., 3(10):763–774, 2011.
  * [155] Alexander S Fokas, Daniel J Cole, Nicholas DM Hine, Stephen A Wells, Michael C Payne, and Alex W Chin.  Evidence of correlated static disorder in the fenna–matthews–olson complex.  J. Phys. Chem. Lett., 8:2350–2356, 2017.
  * [156] Dale E Tronrud, Jianzhong Wen, Leslie Gay, and Robert E Blankenship.  The structural basis for the difference in absorbance spectra for the fmo antenna protein from various green sulfur bacteria.  Photosynth. Res., 100(2):79–87, 2009.
  * [157] Marcel Schmidt am Busch, Frank Müh, Mohamed El-Amine Madjet, and Thomas Renger.  The eighth bacteriochlorophyll completes the excitation energy funnel in the fmo protein.  J. Chem. Theory Comput., 2(2):93–98, 2010.
  * [158] X Jia, Y Mei, JZ Zhang, and Y Mo.  Hybrid qm/mm study of fmo complex with polarized protein-specific charge.  Sci. Rep., 5:17096, 2015.
  * [159] Xiu Juan Yang, Tieneke E Dykstra, and Gregory D Scholes.  Photon-echo studies of collective absorption and dynamic localization of excitation in conjugated polymers and oligomers.  Phys. Rev. B, 71(4):045203, 2005.
  * [160] Ivan G Scheblykin, Arkady Yartsev, Tonu Pullerits, Vidmantas Gulbinas, and Villy Sundström.  Excited state and charge photogeneration dynamics in conjugated polymers.  J. Phys. Chem. B, 111(23):6303–6321, 2007.
  * [161] Tak W Kee.  Femtosecond pump–push–probe and pump–dump–probe spectroscopy of conjugated polymers: New insight and opportunities.  J. Phys. Chem. Lett., 5(18):3231–3240, 2014.
  * [162] D C Dai and A P Monkman.  Femtosecond hot-exciton emission in a ladder-type $\pi$-conjugated rigid-polymer nanowire.  Phys. Rev. B, 87(4):045308, 2013.
  * [163] Nathan P Wells and David A Blank.  Correlated exciton relaxation in poly (3-hexylthiophene).  Phys. Rev. Lett., 100(8):086403, 2008.
  * [164] Julia Kempe.  Quantum random walks: an introductory overview.  Contemp. Phys., 44(4):307–327, 2003.
  * [165] Wichard JD Beenken and Tonu Pullerits.  Excitonic coupling in polythiophenes: comparison of different calculation methods.  J. Chem. Phys., 120(5):2490–2495, 2004.
  * [166] François Dubin, Romain Melet, Thierry Barisien, Roger Grousson, Laurent Legrand, Michel Schott, and Valia Voliotis.  Macroscopic coherence of a single exciton state in an organic quantum wire.  Nat. Phys., 2(1):32–35, 2006.
  * [167] Wai Lun Chan, Manuel Ligges, and X Y Zhu.  The energy barrier in singlet fission can be overcome through coherent coupling and entropic gain.  Nat. Chem., 4(10):840–845, 2012.
  * [168] Frank Würthner, Theo E Kaiser, and Chantu R Saha-Möller.  J-aggregates: From serendipitous discovery to supramolecular engineering of functional dye materials.  Angew. Chem. Int. Ed., 50(15):3376–3410, 2011.
  * [169] K Birgitta Whaley, Aleksey A Kocherzhenko, and Abraham Nitzan.  Coherent and diffusive time scales for exciton dissociation in bulk heterojunction photovoltaic cells.  J. Phys. Chem. C, 118(47):27235–27244, 2014.
  * [170] Natalie Banerji, Sarah Cowan, Mario Leclerc, Eric Vauthey, and Alan J Heeger.  Exciton formation, relaxation, and decay in pcdtbt.  J. Am. Chem. Soc., 132(49):17459–17470, 2010.
  * [171] Tieneke E Dykstra, Emmanuelle Hennebicq, David Beljonne, Johannes Gierschner, Gil Claudio, Eric R Bittner, Jasper Knoester, and Gregory D Scholes.  Conformational disorder and ultrafast exciton relaxation in ppv-family conjugated polymers.  J. Phys. Chem. B, 113(3):656–667, 2009.
  * [172] Arvi Freiberg, Su Lin, Kou Timpmann, and Robert E Blankenship.  Exciton dynamics in fmo bacteriochlorophyll protein at low temperatures.  J. Phys. Chem. B, 101(37):7211–7220, 1997.
  * [173] Kai Chen, Alex J Barker, Matthew E Reish, Keith C Gordon, and Justin M Hodgkiss.  Broadband ultrafast photoluminescence spectroscopy resolves charge photogeneration via delocalized hot excitons in polymer: Fullerene photovoltaic blends.  J. Am. Chem. Soc., 135(49):18502–18512, 2013.
  * [174] Nathan P Wells, Bryan W Boudouris, Marc A Hillmyer, and David A Blank.  Intramolecular exciton relaxation and migration dynamics in poly (3-hexylthiophene).  J. Phys. Chem. C, 111(42):15404–15414, 2007.
  * [175] Ti Wang and Wai-Lun Chan.  Dynamical localization limiting the coherent transport range of excitons in organic crystals.  J. Phys. Chem. Lett., 5(11):1812–1818, 2014.
  * [176] Arun Aby Paraecattil, Serge Beaupré, Mario Leclerc, Jacques-E Moser, and Natalie Banerji.  Intensity dependent femtosecond dynamics in a pbdttpd-based solar cell material.  J. Phys. Chem. Lett., 3(20):2952–2958, 2012.
  * [177] Arun Aby Paraecattil and Natalie Banerji.  Charge separation pathways in a highly efficient polymer: Fullerene solar cell material.  J. Am. Chem. Soc., 136(4):1472–1482, 2014.
  * [178] Zhi Guo, Doyun Lee, Richard D Schaller, Xiaobing Zuo, Byeongdu Lee, TengFei Luo, Haifeng Gao, and Libai Huang.  Relationship between interchain interaction, exciton delocalization, and charge separation in low-bandgap copolymer blends.  J. Am. Chem. Soc., 136(28):10024–10032, 2014.
  * [179] Inchan Hwang, Serge Beaupré, Mario Leclerc, and Gregory D Scholes.  Ultrafast relaxation of charge-transfer excitons in low-bandgap conjugated copolymers.  Chem. Sci., 3(7):2270–2277, 2012.
  * [180] Chad Risko, Michael D McGehee, and Jean-Luc Brédas.  A quantum-chemical perspective into low optical-gap polymers for highly-efficient organic solar cells.  Chem. Sci., 2(7):1200–1218, 2011.
  * [181] Yong Zhang, Steven K Hau, Hin-Lap Yip, Ying Sun, Orb Acton, and Alex K-Y Jen.  Efficient polymer solar cells based on the copolymers of benzodithiophene and thienopyrroledione.  Chem. Mater., 22(9):2696–2698, 2010.
  * [182] Yingping Zou, Ahmed Najari, Philippe Berrouard, Serge Beaupré, Badrou Réda Aïch, Ye Tao, and Mario Leclerc.  A thieno [3, 4-c] pyrrole-4, 6-dione-based copolymer for efficient solar cells.  J. Am. Chem. Soc., 132(15):5330–5331, 2010.
  * [183] N Van den Brande, G Van Lier, F Da Pieve, G Van Assche, B Van Mele, F De Proft, and P Geerlings.  A time dependent dft study of the efficiency of polymers for organic photovoltaics at the interface with pcbm.  RSC Adv., 4(95):52658–52667, 2014.
  * [184] WeiWei Zhang, WanZhen Liang, and Yi Zhao.  Non-condon effect on charge transport in dithiophene-tetrathiafulvalene crystal.  J. Chem. Phys., 133(2):024501, 2010.
  * [185] Fei Pan, Fang Gao, WanZhen Liang, and Yi Zhao.  Nature of low-lying excited states in h-aggregated perylene bisimide dyes: results of td-lrc-dft and the mixed exciton model.  J. Phys. Chem. C, 113(44):14581–14587, 2009.
  * [186] Verner K Thorsmølle, Richard D Averitt, Jure Demsar, DL Smith, S Tretiak, RL Martin, X Chi, BK Crone, AP Ramirez, and AJ Taylor.  Morphology effectively controls singlet-triplet exciton relaxation and charge transport in organic semiconductors.  Phys. Rev. Lett., 102(1):017401, 2009.
  * [187] Justin C Johnson, Arthur J Nozik, and Josef Michl.  High triplet yield from singlet fission in a thin film of 1, 3-diphenylisobenzofuran.  J. Am. Chem. Soc., 132(46):16302–16303, 2010.
  * [188] Mark WB Wilson, Akshay Rao, Jenny Clark, R Sai Santosh Kumar, Daniele Brida, Giulio Cerullo, and Richard H Friend.  Ultrafast dynamics of exciton fission in polycrystalline pentacene.  J. Am. Chem. Soc., 133(31):11830–11833, 2011.
  * [189] Brian J Walker, Andrew J Musser, David Beljonne, and Richard H Friend.  Singlet exciton fission in solution.  Nat. Chem., 5(12):1019, 2013.
  * [190] Mark WB Wilson, Akshay Rao, Kerr Johnson, Simon Gélinas, Riccardo Di Pietro, Jenny Clark, and Richard H Friend.  Temperature-independent singlet exciton fission in tetracene.  J. Am. Chem. Soc., 135(44):16680–16688, 2013.
  * [191] Samuel W Eaton, Leah E Shoer, Steven D Karlen, Scott M Dyar, Eric A Margulies, Brad S Veldkamp, Charusheela Ramanan, Daniel A Hartzler, Sergei Savikhin, Tobin J Marks, et al.  Singlet exciton fission in polycrystalline thin films of a slip-stacked perylenediimide.  J. Am. Chem. Soc., 135(39):14701–14712, 2013.
  * [192] Kolja Kolata, Tobias Breuer, Gregor Witte, and Sangam Chatterjee.  Molecular packing determines singlet exciton fission in organic semiconductors.  ACS nano, 8(7):7377–7383, 2014.
  * [193] Eric A Margulies, Claire E Miller, Yilei Wu, Lin Ma, George C Schatz, Ryan M Young, and Michael R Wasielewski.  Enabling singlet fission by controlling intramolecular charge transfer in $\pi$-stacked covalent terrylenediimide dimers.  Nat. Chem., 8(12):1120, 2016.
  * [194] Christopher J Bardeen.  The structure and dynamics of molecular excitons.  Annu. Rev. Phys. Chem., 65:127–148, 2014.
  * [195] Andrew J Musser, Matz Liebel, Christoph Schnedermann, Torsten Wende, Tom B Kehoe, Akshay Rao, and Philipp Kukura.  Evidence for conical intersection dynamics mediating ultrafast singlet exciton fission.  Nat. Phys., 11(4):352, 2015.
  * [196] Artem A Bakulin, Sarah E Morgan, Tom B Kehoe, Mark WB Wilson, Alex W Chin, Donatas Zigmantas, Dassia Egorova, and Akshay Rao.  Real-time observation of multiexcitonic states in ultrafast singlet fission using coherent 2d electronic spectroscopy.  Nat. Chem., 8(1):16, 2016.
  * [197] You-Dan Zhang, Yishi Wu, Yanqing Xu, Qiang Wang, Ke Liu, Jian-Wei Chen, Jing-Jing Cao, Chunfeng Zhang, Hongbing Fu, and Hao-Li Zhang.  Excessive exoergicity reduces singlet exciton fission efficiency of heteroacenes in solutions.  J. Am. Chem. Soc., 138(21):6739–6745, 2016.
  * [198] William Shockley and Hans J Queisser.  Detailed balance limit of efficiency of p-n junction solar cells.  J. Appl. Phys., 32(3):510–519, 1961.
  * [199] MC Hanna and AJ Nozik.  Solar conversion efficiency of photovoltaic and photoelectrolysis cells with carrier multiplication absorbers.  J. Appl. Phys., 100(7):074510, 2006.
  * [200] Murad JY Tayebjee, Angus A Gray-Weale, and Timothy W Schmidt.  Thermodynamic limit of exciton fission solar cell efficiency.  J. Phys. Chem. Lett., 3(19):2749–2754, 2012.
  * [201] Guohua Tao.  Electronically nonadiabatic dynamics in singlet fission: A quasi-classical trajectory simulation.  J. Phys. Chem. C, 118(31):17299–17305, 2014.
  * [202] Fatemeh Mirjani, Nicolas Renaud, Natalie Gorczak, and Ferdinand C Grozema.  Theoretical investigation of singlet fission in molecular dimers: The role of charge transfer states and quantum interference.  J. Phys. Chem. C, 118(26):14192–14199, 2014.
  * [203] Niels H Damrauer and Jamie L Snyder.  Symmetry-directed control of electronic coupling for singlet fission in covalent bis–acene dimers.  J. Phys. Chem. Lett., 6(22):4456–4462, 2015.
  * [204] Maria A Castellanos and Pengfei Huo.  Enhancing singlet fission dynamics by suppressing destructive interference between charge-transfer pathways.  J. Phys. Chem. Lett., 8(11):2480–2488, 2017.
  * [205] Michael Kasha.  Energy transfer mechanisms and the molecular exciton model for molecular aggregates.  Radiat. Res., 20(1):55–70, 1963.
  * [206] M Kasha, HR Rawls, and M Ashraf El-Bayoumi.  The exciton model in molecular spectroscopy.  Pure Appl. Chem., 11(3-4):371–392, 1965.
  * [207] Frank C Spano and Carlos Silva.  H-and j-aggregate behavior in polymeric semiconductors.  Annu. Rev. Phys. Chem., 65:477–500, 2014.
  * [208] Nicholas J Hestand and Frank C Spano.  Molecular aggregate photophysics beyond the kasha model: Novel design principles for organic materials.  Acc. Chem. Res., 50(2):341–350, 2017.
  * [209] Sven Larsson.  Electron transfer in chemical and biological systems. orbital rules for nonadiabatic transfer.  J. Am. Chem. Soc., 103(14):4034–4040, 1981.
  * [210] Abraham Nitzan.  Chemical dynamics in condensed phases: relaxation, transfer and reactions in condensed molecular systems.  Oxford university press, 2006.
  * [211] Piotr Petelenz and Mateusz Snamina.  Locally broken crystal symmetry facilitates singlet exciton fission.  J. Phys. Chem. Lett., 7(10):1913–1916, 2016.
  * [212] Catherine M Mauck, Patrick E Hartnett, Eric A Margulies, Lin Ma, Claire E Miller, George C Schatz, Tobin J Marks, and Michael R Wasielewski.  Singlet fission via an excimer-like intermediate in 3, 6-bis (thiophen-2-yl) diketopyrrolopyrrole derivatives.  J. Am. Chem. Soc., 138(36):11749–11761, 2016.
  * [213] Timothy C. Berkelbach, Mark S. Hybertsen, and David R. Reichman.  Microscopic theory of singlet exciton fission. II. application to pentacene dimers and the role of superexchange.  J. Chem. Phys., 138(11):114103, mar 2013.
  * [214] Eric C Greyson, Josh Vura-Weis, Josef Michl, and Mark A Ratner.  Maximizing singlet fission in organic dimers: theoretical investigation of triplet yield in the regime of localized excitation and fast coherent electron transfer.  J. Phys. Chem. B, 114(45):14168–14177, 2010.
  * [215] Wai-Lun Chan, Timothy C Berkelbach, Makenzie R Provorse, Nicholas R Monahan, John R Tritsch, Mark S Hybertsen, David R Reichman, Jiali Gao, and X-Y Zhu.  The quantum coherent mechanism for singlet fission: experiment and theory.  Acc. Chem. Res., 46(6):1321–1329, 2013.
  * [216] Linjun Wang, Yoann Olivier, Oleg V Prezhdo, and David Beljonne.  Maximizing singlet fission by intermolecular packing.  J. Phys. Chem. Lett., 5(19):3345–3353, 2014.
  * [217] Nicholas J Hestand and Frank C Spano.  Interference between coulombic and ct-mediated couplings in molecular aggregates: H-to j-aggregate transformation in perylene-based $\pi$-stacks.  J. Chem. Phys, 143(24):244707, 2015.
  * [218] Shane R Yost, Jiye Lee, Mark WB Wilson, Tony Wu, David P McMahon, Rebecca R Parkhurst, Nicholas J Thompson, Daniel N Congreve, Akshay Rao, Kerr Johnson, et al.  A transferable model for singlet-fission kinetics.  Nat. Chem., 6(6):492, 2014.
  * [219] Ryan D Pensack, Andrew J Tilley, Sean R Parkin, Tia S Lee, Marcia M Payne, Dong Gao, Ashlee A Jahnke, Daniel G Oblinsky, Peng-Fei Li, John E Anthony, et al.  Exciton delocalization drives rapid singlet fission in nanoparticles of acene derivatives.  J. Am. Chem. Soc., 137(21):6790–6803, 2015.