# Exciton dynamics from the mapping approach to surface hopping: Comparison
with Förster and Redfield theories

Johan E. Runeson<EMAIL_ADDRESS>Department of Chemistry,
University of Oxford, Physical and Theoretical Chemistry Laboratory, South
Parks Road, Oxford, OX1 3QZ, UK  Thomas P. Fay  Department of Chemistry,
University of California, Berkeley, California 94720, USA  David E.
Manolopoulos  Department of Chemistry, University of Oxford, Physical and
Theoretical Chemistry Laboratory, South Parks Road, Oxford, OX1 3QZ, UK

###### Abstract

We compare the recently introduced multi-state mapping approach to surface
hopping (MASH) with the Förster and Redfield theories of excitation energy
transfer. Whereas Förster theory relies on weak coupling between chromophores,
and Redfield theory assumes the electronic excitations to be weakly coupled to
fast chromophore vibrations, MASH is free from any perturbative or Markovian
approximations. We illustrate this with an example application to the rate of
energy transfer in a Frenkel-exciton dimer, showing that MASH interpolates
correctly between the opposing regimes in which the Förster and Redfield
results are reliable. We then compare the three methods for a realistic model
of the Fenna–Matthews–Olson complex with a structured vibrational spectral
density and static disorder in the excitation energies. In this case there are
no exact results for comparison so we use MASH to assess the validity of
Förster and Redfield theories. We find that Förster theory is the more
accurate of the two on the picosecond timescale, as has been shown previously
for a simpler model of this particular light-harvesting complex. We also
explore various ways to sample the initial electronic state in MASH and find
that they all give very similar results for exciton dynamics.

##  I Introduction

Revealing the mechanisms of excitation energy transfer is fundamental to the
study of biological photosynthetic systems and has motivated methodological
development for several decades.Mirkovic  _et al._ (2017); Segatta  _et al._
(2019) Excitonic systems typically consist of an assembly of chromophores
(sites) embedded in a protein or solvent environment. Traditionally,
excitation energy transfer has been divided into two regimes, depending on the
relative magnitudes of the vibrational relaxation time and the inverse of the
inter-site coupling strength. May and Kühn (2011); Ishizaki  _et al._ (2010)
For chromophores that are weakly coupled to each other (but have rapid
vibrational relaxation, which usually means strong coupling to the
environment) the transfer is referred to as ‘incoherent’ and proceeds via
hopping between sites. This regime can be well described by Förster theory. In
the opposite limit of slow vibrational relaxation (or strong coupling between
the chromophores), the transfer is referred to as ‘coherent’ and proceeds
between delocalized exciton states. This regime can be described by a
perturbative master equation such as Redfield theory. We emphasize that in
either regime, the transfer of population is often well described by a kinetic
rate equation, and that the electronic coherences are both small and short-
lived in an appropriate electronic basis.Alvertis  _et al._ (2020); Cao  _et
al._ (2020)

The main challenge is that many systems of interest tend to appear in the
intermediate regime where the two timescales are comparable.Scholes (2018);
Sneyd, Beljonne, and Rao (2022) For simple system–bath models one can solve
the dynamics non-perturbatively using fully quantum methods such as the
hierarchical equations of motion (HEOM).Tanimura and Kubo (1989) Such a study
of a two-site model has shown that exciton transfer is fastest in the
intermediate regime, where the rate is not adequately captured by either of
the perturbative theories mentioned above.Ishizaki and Fleming (2009a)
However, fully quantum methods like HEOM and quasi-adiabatic path
integralsMakri and Makarov (1995); Makri (2020) are only practical for
harmonic models and they can be hard to converge for strong system–environment
couplings. It is therefore important to develop accurate methods that can
describe nonadiabatic transitions more generally, including in systems with
anharmonic atomistic potentials.

So far, the most successful strategy to model nonadiabatic transitions has
been Tully’s fewest switches surface hopping (FSSH).Tully (1990) Since it was
first proposed in 1990, this stochastic algorithm has become immensely popular
in photochemistry and is widely implemented in open software. However, it
continues to suffer from long-standing issues often referred to as
‘overcoherence’, despite much effort on the development of more or less ad hoc
‘decoherence corrections’.Wang  _et al._ (2020)

Recently, Mannouch and Richardson have proposed a different strategy based on
a phase-space mapping of two-state systems onto a spin degree of
freedom.Mannouch and Richardson (2023) This so-called ‘mapping approach to
surface hopping’ (MASH) uses a deterministic algorithm that, in contrast to
FSSH, hops to the adiabatic surface with the largest instantaneous population.
This strategy has many appealing features. Firstly, it removes all ambiguity
about the need for velocity rescaling/reversal for successful/unsuccessful
hops. Secondly, it replaces ad hoc decoherence corrections with a rigorous
‘quantum jump’ procedure, even without which it has been found to be more
accurate than FSSH in a range of benchmark applications.Mannouch and
Richardson (2023) Thirdly, unlike FSSH, it can correctly describe the
transition between adiabatic and non-adiabatic rates in the spin-boson model
and it recovers Marcus theory in the limit of a perturbative inter-state
coupling.Lawrence, Mannouch, and Richardson (2023) Finally, MASH has been
proven to relax to the correct quantum–classical equilibrium distribution for
ergodic systems, a feature that is not shared by any other nonadiabatic
trajectory method.Amati, Mannouch, and Richardson (2023)

In its original formulation, MASH was limited to systems with two electronic
states, but an adaptation to multiple states has recently been
proposed.Runeson and Manolopoulos (2023) Compared to the original approach,
the multi-state formulation differs in the way it calculates electronic
observables. The multi-state estimators are constructed to be equivariant
under unitary basis transformations, meaning that populations and coherences
are treated on the same footing and observables can be directly evaluated in
any basis. With these estimators, multi-state MASH provably relaxes to the
correct quantum–classical equilibrium for a general $N$-state system, in any
basis.Runeson and Manolopoulos (2023)

Unfortunately, the equivariant multi-state formulation of MASHRuneson and
Manolopoulos (2023) does not reduce to the original formulationMannouch and
Richardson (2023) in the two-state case. The two methods have nevertheless
been shown to be of comparable accuracy for a wide variety of two-state
systems.Runeson and Manolopoulos (2023) A noteworthy exception was found for a
spin-boson model in the Marcus inverted regime, for which the original
formulation was more accurate than the multi-state formulation.Runeson and
Manolopoulos (2023) However, we shall show below that this conclusion does not
hold generally, and that for a two-site exciton transfer model the two
formulations lead to essentially the same rates, even in the inverted regime.

In this article, we apply multi-state MASH to a more challenging set of
exciton systems than previously considered. In particular, we investigate the
transition between the regimes of ‘coherent’ and ‘incoherent’ rate theories
and compare with fully quantum benchmark results. We then go on to investigate
exciton transfer in the Fenna–Matthews–Olson (FMO) complex, using an
experimentally measured spectral density to account for coupling to
vibrations, and including static disorder in the site energies. Comparison to
simple rate theories indicates that FMO is better described by Förster theory
than by Redfield theory, despite the popularity of the latter in the exciton
literature.

We also analyse in detail the choice of initial conditions in multi-state
MASH. The choice used in Ref. Runeson and Manolopoulos, 2023 is not unique
and, more importantly, not basis-equivariant, in contrast to the treatment of
electronic observables. We exemplify why this could become a problem for a
system with three states, and consider alternative sets of initial conditions.
In particular, we present an approach that overcomes the objection for three
states and restores the equivariance of the initial distribution. Upon
comparison for exciton transfer in the dimer and FMO models, however, we find
that the initial distribution has little practical influence on the dynamics.
Based on these results, we conclude that using the simplest set of initial
conditions is well-justified for the kind of system we consider here.

##  II Perturbative rate theories

As a prototypical model for exciton energy transfer, we consider the Frenkel-
exciton Hamiltonian,

$H=H_{\rm S}+H_{\rm B}+H_{\rm SB}.$ (1)

The first term is the system Hamiltonian in the basis of localized pigment
(‘site’) excitations

$H_{\rm S}=\sum_{n=1}^{N}\epsilon_{n}|n\rangle\langle
n|+\sum_{n>m}J_{nm}(|n\rangle\langle m|+|m\rangle\langle n|),$ (2)

where $\\{\epsilon_{n}\\}$ and $\\{J_{nm}\\}$ are the site energies and the
inter-site couplings, respectively. To model interaction with vibrational and
solvent degrees of freedom, the on-diagonal site energies are linearly coupled
to a harmonic bath,

$H_{\rm
B}=\sum_{n=1}^{N}\sum_{j=1}^{f}\left(\frac{p_{j,n}^{2}}{2}+\frac{1}{2}\omega_{j}^{2}q_{j,n}^{2}\right),$
(3) $H_{\rm SB}=\sum_{j=1}^{f}\sum_{n=1}^{N}\kappa_{j}q_{j,n}|n\rangle\langle
n|.$ (4)

Here, $p_{j,n}$ and $q_{j,n}$ are mass-scaled ($m_{j,n}=1$) momentum and
coordinate variables for the vibrational modes. The bath frequencies and
vibrational couplings are specified through the spectral density

$J(\omega)=\frac{\pi}{2}\sum_{j}\frac{\kappa_{j}^{2}}{\omega_{j}}\delta(\omega-\omega_{j}),$
(5)

and all sites are assumed to be coupled to identical and independent baths. As
a measure of the overall system–bath coupling strength, we define the bath
reorganization energy

$\lambda=\frac{1}{\pi}\int_{0}^{\infty}\frac{J(\omega)}{\omega}\mathrm{d}\omega=\sum_{j}\frac{\kappa_{j}^{2}}{2\omega_{j}^{2}}.$
(6)

The system–bath interaction involves the bath operator
$B_{n}=\sum_{j}\kappa_{j}q_{j,n}$, which is the energy gap between a localized
excitation on site $n$ and the ground state. In the following, we drop the
index $n$ since the baths are identical. To characterize the dynamics of the
bath, it is useful to define the autocorrelation function of the bath
operator,

$C(t)=\mathrm{Tr}_{\rm B}[\mathrm{e}^{-\beta{H}_{\rm B}}B(0)B(t)],$ (7)

where $\beta=1/(k_{\rm B}T)$ and the time-dependence refers to dynamics under
${H}_{\rm B}$. Using known expressions for the thermal correlation functions
of a harmonic oscillator, one can show that

$\displaystyle C(t)$
$\displaystyle=\frac{1}{\pi}\int_{0}^{\infty}\mathrm{d}\omega\,J(\omega)\left[\coth\left(\frac{\beta\omega}{2}\right)\cos\omega
t+\mathrm{i}\sin\omega t\right]$ (8)
$\displaystyle=\frac{1}{\pi}\int_{0}^{\infty}\mathrm{d}\omega\,J(\omega)\left[(1+n(\omega))\mathrm{e}^{\mathrm{i}\omega
t}+n(\omega)\mathrm{e}^{-\mathrm{i}\omega t}\right]$ (9)

where $n(\omega)=1/(\mathrm{e}^{\beta\omega}-1)$ is the Bose–Einstein
distribution and throughout we use units where $\hbar=1$.

###  II.1 Redfield theory

In the limit of weak system–bath coupling, the effect of the bath on the
excitonic system can be described by a second-order perturbative master
equation. The perturbation expansion is usually truncated in the eigenbasis of
$H_{\rm S}$, which is called the _exciton_ basis. Diagonalizing the system
Hamiltonian gives $H_{\rm S}=\sum_{\mu}\omega_{\mu}|\mu\rangle\langle\mu|$,
where $\omega_{\mu}$ is an eigenenergy and $|\mu\rangle=\sum_{n}U^{-1}_{\mu
n}|n\rangle=\sum_{n}U_{n\mu}|n\rangle$ is an exciton state. We also assume
that $C(t)$ decays faster than the timescale of the system dynamics. Under
these two assumptions (a weak and fast bath), a standard textbook
derivationMay and Kühn (2011) leads to a Markovian master equation for the
reduced density matrix of the system,

$\frac{\partial}{\partial
t}\rho_{\mu\nu}(t)=-\mathrm{i}\omega_{\mu\nu}\rho_{\mu\nu}(t)+\sum_{\mu^{\prime}\nu^{\prime}}R_{\mu\nu\mu^{\prime}\nu^{\prime}}\rho_{\mu^{\prime}\nu^{\prime}}(t),$
(10)

where $\omega_{\mu\nu}=\omega_{\mu}-\omega_{\nu}$ is the energy gap between
the corresponding eigenstates of $H_{\rm S}$. The second term involves the
Redfield tensor

$\displaystyle R_{\mu\nu\mu^{\prime}\nu^{\prime}}$
$\displaystyle=\Gamma_{\nu^{\prime}\nu\mu\mu^{\prime}}+\Gamma^{*}_{\mu^{\prime}\mu\nu\nu^{\prime}}$
$\displaystyle-\delta_{\nu\nu^{\prime}}\sum_{\kappa}\Gamma_{\mu\kappa\kappa\mu^{\prime}}-\delta_{\mu\mu^{\prime}}\sum_{\kappa}\Gamma^{*}_{\nu\kappa\kappa\nu^{\prime}},$
(11)

which is expressed in terms of the damping tensor

$\Gamma_{\mu\nu\mu^{\prime}\nu^{\prime}}=\sum_{n}\langle\mu|n\rangle\langle
n|\nu\rangle\langle\nu^{\prime}|n\rangle\langle
n|\mu^{\prime}\rangle\tilde{C}(\omega_{\mu\nu}).$ (12)

Here, $\tilde{C}(\omega)$ is the Fourier–Laplace transform of the bath
correlation function,

$\tilde{C}(\omega)=\int_{0}^{\infty}\mathrm{d}t\;\mathrm{e}^{-\mathrm{i}\omega
t}C(t).$ (13)

Inserting the expression for $C(t)$ from Eq. (8) gives

$\displaystyle\operatorname{Re}\tilde{C}(\omega)$
$\displaystyle=J(\omega)(1+n(\omega))+J(-\omega)n(-\omega)$ (14)
$\displaystyle\operatorname{Im}\tilde{C}(\omega)$
$\displaystyle=\frac{1}{\pi}\mathrm{P}\int_{-\infty}^{\infty}\mathrm{d}\omega^{\prime}\;\frac{\operatorname{Re}\tilde{C}(\omega)}{\omega^{\prime}-\omega},$
(15)

where $\mathrm{P}$ denotes principal value and, with the notation used in this
paper, $J(\omega<0)=0$, so only one term on the right-hand side of Eq. (14) is
non-zero. Finally, if the diagonal elements of ${\rho}$ (the exciton
populations) are only weakly influenced by the off-diagonal elements (the
exciton coherences), one can replace Eq. (10) by a kinetic rate equation for
the populations (the secular approximation),

$\frac{\partial}{\partial
t}\rho_{\mu\mu}(t)=-\sum_{\mu\neq\nu}(k_{\nu\to\mu}\rho_{\nu\nu}-k_{\mu\to\nu}\rho_{\mu\mu})$
(16)

with the Redfield rate constants

$k_{\mu\to\nu}=2\sum_{n}\langle\mu|n\rangle\langle
n|\nu\rangle\langle\nu|n\rangle\langle
n|\mu\rangle\operatorname{Re}\tilde{C}(\omega_{\mu\nu}).$ (17)

###  II.2 Förster theory

In the opposite limit of strong system–bath coupling, one may instead pick the
perturbation parameter to be the intersite coupling $J_{nm}$. This is the
Förster-type incoherent hopping limit, in which the subsystem follows a
kinetic rate equation in the _site_ basis. The Förster rate constants areYang
and Fleming (2002)

$k_{n\to
m}=2|J_{nm}|^{2}\operatorname{Re}\int_{0}^{\infty}\mathrm{d}t\;F_{n}^{*}(t)A_{m}(t)$
(18)

where $F_{n}(t)$ and $A_{m}(t)$ are the flourescence and absorption lineshape
functions, which based on the cumulant expansion technique can be written as

$\displaystyle F_{n}^{*}(t)$
$\displaystyle=\mathrm{e}^{+\mathrm{i}(\epsilon_{n}-\lambda)t-g(t)}$ (19)
$\displaystyle A_{m}(t)$
$\displaystyle=\mathrm{e}^{-\mathrm{i}(\epsilon_{m}+\lambda)t-g(t)}.$ (20)

The function $g(t)$ is

$g(t)=\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\,C(t_{2}),$
(21)

where $C(t)$ is the bath correlation function in Eq. (8). For the Frenkel-
exciton model with identical baths, the rate expression in Eq. (18) reduces to

$k_{n\to
m}=2|J_{nm}|^{2}\operatorname{Re}\int_{0}^{\infty}\mathrm{d}t\;\mathrm{e}^{\mathrm{i}(\epsilon_{n}-\epsilon_{m})t-2g(t)},$
(22)

where

$g(t)=\sum_{j}\frac{\kappa_{j}^{2}}{2\omega_{j}^{3}}\left(\coth\left(\frac{\beta\omega_{j}}{2}\right)(1-\cos\omega_{j}t)+\mathrm{i}\sin\omega_{j}t\right).$
(23)

##  III Multi-state MASH

An alternative strategy is to simulate the bath dynamics explicitly with
surface hopping. To this end, we rewrite the Frenkel-exciton Hamiltonian in
Eq. (1) as

$\hat{H}(p,q)=\sum_{j,n}\frac{p_{j,n}^{2}}{2}+\hat{V}(q),$ (24)

where

$\hat{V}(q)=\sum_{nm}V_{nm}(q)|n\rangle\langle m|$ (25)

is a (diabatic) potential operator with matrix elements

$V_{nm}(q)=\bigg{(}\epsilon_{n}+\sum_{j}\kappa_{j}q_{j,n}+\sum_{j,l}\frac{1}{2}\omega_{j}^{2}q_{j,l}^{2}\bigg{)}\delta_{nm}+J_{nm}.$
(26)

Surface hopping is almost always run on adiabatic surfaces, i.e., in the
eigenbasis of $\hat{V}(q)$,

$\hat{V}(q)=\sum_{a}V_{a}(q)|a(q)\rangle\langle(q)|.$ (27)

We will refer to $p,q$ as _nuclear_ variables (physically, they represent
intramolecular vibrations as well as collective modes of the solvent). Their
dynamics is coupled to the electronic wavefunction $|\psi\rangle$, which can
be expanded in the diabatic or the adiabatic basis as

$|\psi\rangle=\sum_{n}c_{n}|n\rangle=\sum_{a}c_{a}|a(q)\rangle.$ (28)

The coefficients obey the Schrödinger equation

$\dot{c}_{n}=-\mathrm{i}\sum_{m}V_{nm}(q)c_{m}$ (29)

or (equivalently)

$\dot{c}_{a}=-\mathrm{i}V_{a}(q)c_{a}-\sum_{j}\dot{q}_{j}\sum_{b}d_{ab}^{j}(q)c_{b},$
(30)

where $d_{ab}^{j}(q)=\langle a(q)|\nabla_{j}|b(q)\rangle$ is a nonadiabatic
coupling matrix element. (Here and in the rest of this section, the index $j$
runs over the nuclear degrees of freedom of all sites, not just one.)

In MASH, the nuclear trajectories evolve according to the classical equations
of motion

$\displaystyle\dot{q}_{j}$ $\displaystyle=p_{j}/m_{j}$ (31a)
$\displaystyle\dot{p}_{j}$ $\displaystyle=-\langle
a(q)|\nabla_{j}\hat{V}|a(q)\rangle$ (31b)

where $a$ is the adiabatic state with the largest instantaneous population
$|c_{a}|^{2}$. Effectively, this means that the nuclei evolve on a potential
with abrupt steps. By introducing the classical state projectors,

$\Theta_{n}=\begin{cases}1&\text{if }|c_{n}|^{2}>|c_{k}|^{2}\;\forall k\neq
n\\\ 0&\text{otherwise},\end{cases}$ (32)

one can express the effective potential perceived by the nuclei as

$V_{\rm eff}=\sum_{a}V_{a}(q)\Theta_{a}(c).$ (33)

At each instant, precisely one of the $\Theta_{a}(c)$ factors in the sum of
Eq. (33) is non-zero. Whenever a new state $b$ reaches a higher population
than the current state $a$, the nuclei meet a potential step $V_{b}-V_{a}$.
When crossing such a potential step, the momentum is rescaled so as to
conserve energy, and if there is insufficient kinetic energy to overcome the
step, the momentum is instead reversed. In multi-state MASH, the momentum
component subject to rescaling/reversal is the projection onto the direction
of a vector $v$ with elementsRuneson and Manolopoulos (2023)

$v_{j}=\sum_{a^{\prime}}\operatorname{Re}\left[c^{*}_{a^{\prime}}(d_{a^{\prime}a}^{j}c_{a}-d_{a^{\prime}b}^{j}c_{b})\right].$
(34)

Because the hops occur deterministically, it is straightforward to adjust the
timestep if necessary to better resolve a given hopping event (whereas in a
stochastic algorithm, changing the timestep would change the locations of the
hops).

###  III.1 Estimators

Having defined the dynamics of each trajectory, we next address how to measure
observables. Consider a process starting in a pure electronic state
$|i\rangle\langle i|$ with a (normalized) nuclear density $\rho_{\rm B}(p,q)$.
The time-dependent expectation value of an observable $O$ is then

$\langle O(t)\rangle=\mathrm{Tr}[\rho_{\rm B}(p,q)|i\rangle\langle
i|(0)\,\hat{O}(t)],$ (35)

where the trace runs over nuclear as well as electronic degrees of freedom.
The corresponding expression in multi-state MASH is the phase-space integral

$\langle
O(t)\rangle\approx\frac{\int\mathrm{d}p\,\mathrm{d}q\int_{|c|=1}\mathrm{d}c\,\rho_{\rm
B}(p,q)\rho_{i}(c)\,O(p_{t},q_{t},c_{t})}{\int\mathrm{d}p\,\mathrm{d}q\int_{|c|=1}\mathrm{d}c\,\rho_{\rm
B}(p,q)\rho_{i}(c)}$ (36)

where $\int_{|c|=1}\mathrm{d}c$ is an integral over all normalized electronic
wavefunctions. A simple way to sample this integral is generate $2N$ normal
deviates $\\{x_{n},y_{n}\\}_{n=1}^{N}$ and set
$c_{n}=(x_{n}+\mathrm{i}y_{n})/\sqrt{\sum_{k}(x_{k}^{2}+y_{k}^{2})}$.

In the following, we consider the case where $\hat{O}_{n}=|n\rangle\langle n|$
is an electronic population, for which the time-dependent estimator is simply
$O_{n}(c_{t})$. There are multiple ways to construct this estimator.Runeson
and Manolopoulos (2023) If we are interested in an adiabatic population, then
the state projector $\Theta_{a}(c_{t})$ is a natural choice that is consistent
with the adiabatic surface the nuclei are evolving on. However, the state
projector is not a good estimator for diabatic populations because it does not
transform correctly under unitary basis transformations.Runeson and
Manolopoulos (2023)

Another estimator that would transform correctly between bases is the
Ehrenfest population $|c_{n}|^{2}$, but inserting this choice into Eq. (36)
leads to the wrong long-time equilibrium populations. In Ref. Runeson and
Manolopoulos, 2023 it was shown that a simple estimator that fulfils both
criteria (equivariance under unitary basis transformations and consistency
with the quantum-classical equilibrium populations) is

$O_{n}(c)=\alpha_{N}|c_{n}|^{2}+\beta_{N}$ (37)

where

$\alpha_{N}=\frac{N-1}{\sum_{k=1}^{N}(1/k)-1}\qquad\beta_{N}=\frac{1-\alpha_{N}}{N}$
(38)

are two scalars that require no more information than the number of states.
This is the population estimator that is used in multi-state MASH
calculations.Runeson and Manolopoulos (2023)

###  III.2 Initial conditions

What remains to be defined is the choice of initial electronic distribution
$\rho_{i}(c)$ in Eq. (36). This distribution is not unique, in the sense that
many (quasi)probability distributions will fulfil the initial condition

$\frac{\int_{|c|=1}\mathrm{d}c\,\rho_{i}(c)O_{j}(c)}{\int_{|c|=1}\mathrm{d}c\,\rho_{i}(c)}=\delta_{ij}.$
(39)

In the following we consider a few options.

Figure 1: Schematical overview of the MASH initial conditions considered in
this article for the case of two (top row) and three (bottom row) states. In
the two-state case, $\sigma_{x}=2\operatorname{Re}c_{1}^{*}c_{2}$,
$\sigma_{y}=2\operatorname{Im}c_{1}^{*}c_{2}$, and the vertical axis
$\sigma_{z}=|c_{1}|^{2}-|c_{2}|^{2}$ corresponds to polarization in a given
diabatic basis, while the tilted axis in the right column corresponds to
polarization in the adiabatic basis. Dark blue shading indicates a higher
weight and red shading a negative weight. The initial conditions in the first
three columns are used together with the equivariant time-dependent observable
in Eq. (37). The original MASH method uses an alternative prescription that is
(so far) limited to two states.

####  III.2.1 Cap initial condition

A simple choice is

$\rho_{i}(c)=\Theta_{i}(c),$ (40)

which was used in Ref. Runeson and Manolopoulos, 2023. What this equation
implies is that the initial state is chosen randomly from the region where
$|c_{i}|^{2}$ is the largest population. The left column in Fig. 1 visualizes
this region for two and three-level systems. We will refer to these regions as
‘caps’ on the sphere/simplex and consequently to Eq. (40) as the ‘cap’ initial
distribution.

Although the cap distribution has been found to be accurate in a variety of
benchmark calculations,Runeson and Manolopoulos (2023) it is not basis-
equivariant, unlike the population estimator in Eq. (37). A related issue is
that the diabatic basis is not unique, which makes the diabatic state
projector ambigously defined. As a simple example, consider a three-state
system in which we want to start from state 1. Suppose we sample the vector
$c=\frac{1}{\sqrt{3}}(\sqrt{1+2\epsilon},\sqrt{1-\epsilon},\sqrt{1-\epsilon})$,
where $\epsilon$ is some number in the range $0<\epsilon<1$. This vector has
$\Theta_{1}(c)=1$ and would therefore contribute to the dynamics. But if we
define a new set of diabatic basis vectors $|\tilde{1}\rangle=|1\rangle$,
$|\tilde{2}\rangle=(|2\rangle+|3\rangle)/\sqrt{2}$,
$|\tilde{3}\rangle=(|2\rangle-|3\rangle)/\sqrt{2}$, then for
$\epsilon<\frac{1}{4}$ the maximally populated state is $|\tilde{2}\rangle$.
Thus, even though the ket of the initial state is unchanged, the same vector
$c$ would no longer contribute to the dynamics.

####  III.2.2 Focused initial condition

Perhaps the most intuitive choice for the initial wavefunction corresponding
to $|i\rangle\langle i|$ would be to set $c_{i}=1$ and $c_{j\neq i}=0$,
corresponding to a pole on the sphere or a corner of the simplex. This
approach is standard in Ehrenfest dynamics and (provided $i$ is an adiabat) in
FSSH. However, for MASH it would violate the constraint in Eq. (39). The
reason is that in the initial corner of the simplex, $O_{i}(c)>1$ and
$O_{j\neq i}(c)<0$.

The analogous initial condition for MASH with the correct initial value would
be to start from wavefunctions $c$ for which $O_{i}(c)=1$ and $O_{j\neq
i}(c)=0$. Such wavefunctions are confined to circles on the sphere and
isolated points on the simplex, as shown in the second column of Fig. 1. The
circle (point) is defined by $|c_{i}|^{2}=(1-\beta_{N})/\alpha_{N}$ and
$|c_{j\neq i}|^{2}=-\beta_{N}/\alpha_{N}$. These conditions fix the magnitudes
of all components of $c$, leaving the phases to be sampled uniformly from
$[0,2\pi)$. The resulting ‘focused’ initial distribution can be written as

$\rho_{i}(c)=\delta\left(|c_{i}|^{2}-\frac{1-\beta_{N}}{\alpha_{N}}\right)\prod_{j\neq
i}\delta\left(|c_{j}|^{2}+\frac{\beta_{N}}{\alpha_{N}}\right).$ (41)

An advantage of this choice is that each trajectory is initialized with
physical population observables, so it may be possible to use fewer
trajectories than with the cap initial condition to reach statistical
convergence. Nevertheless, the focused distribution is also not basis-
equivariant.

####  III.2.3 Equivariant initial condition

In this section, we derive a quasiprobability distribution $\rho_{i}$ that
transforms correctly under unitary basis transformations. To satisfy the
condition in Eq. (39), the simplest approach is to try the same functional
form as the time-dependent observable in Eq. (37), i.e.

$\rho_{i}(c)=a_{N}|c_{i}|^{2}+b_{N}$ (42)

with some constants $a_{N},b_{N}$ that need not be the same as
$\alpha_{N},\beta_{N}$. The resulting $\rho_{i}(c)$ need not be positive
definite since we can multiply each $c$ sampled from the $|c|=1$ sphere by a
weight that is positive or negative.

To evaluate the integral in Eq. (39), we make use of the following moments:

$\langle|c_{n}|^{2}\rangle=\frac{1}{N},\qquad\langle|c_{n}|^{2}|c_{m}|^{2}\rangle=\frac{\delta_{nm}+1}{N(N+1)}$
(43)

where

$\langle
f\rangle\equiv\frac{\int_{|c|=1}\mathrm{d}c\,f}{\int_{|c|=1}\mathrm{d}c}.$
(44)

These expectation values can be derived using standard formulas for integrals
over a sphere.Folland (2001) With the help of the moments in Eq. (43), Eq.
(39) reduces to two equations in two unknowns, with the solution

$a_{N}=\frac{N+1}{\alpha_{N}},\qquad b_{N}=-1-a_{N}\beta_{N}.$ (45)

Note that there is some similarity between the Roman constants and the Greek
ones in Eq. (38). For example, $\alpha_{N}$ and $\beta_{N}$ are related
through $N\beta_{N}=1-\alpha_{N}$, and likewise $a_{N}$ and $b_{N}$ are
related through $Nb_{N}=1-a_{N}$. This means that not only $O_{n}(c)$ but also
$\rho_{i}(c)$ involves a scaling relative to the centre of the simplex: after
inserting $\beta_{N}$ and $b_{N}$ we get

$O_{n}(c)=\frac{1}{N}+\alpha_{N}\left(|c_{n}|^{2}-\frac{1}{N}\right)$ (46)

and

$\rho_{i}(c)=\frac{1}{N}+a_{N}\left(|c_{i}|^{2}-\frac{1}{N}\right).$ (47)

The special status of the centre of the simplex is analogous to the special
role of the identity operator in phase-space mapping methods.Saller, Kelly,
and Richardson (2019); Gao  _et al._ (2020)

####  III.2.4 Original MASH

For two-level systems, the original MASH method of Mannouch and
RichardsonMannouch and Richardson (2023) uses an alternative prescription in
which diabatic observables are first converted to the adiabatic basis, where
populations and coherences are then measured with different estimators.
Explicitly, for $a\neq a^{\prime}$ and $b\neq b^{\prime}$, their prescription
is (in our notation)

$\displaystyle\mathrm{Tr}[|a\rangle\langle a|(0)|b\rangle\langle b|(t)]$
$\displaystyle\mapsto\frac{2\int\mathrm{d}c\,W_{a}(c)\Theta_{a}(c)\Theta_{b}(c_{t})}{\int\mathrm{d}c}$
(48) $\displaystyle\mathrm{Tr}[|a\rangle\langle a|(0)|b\rangle\langle
b^{\prime}|(t)]$
$\displaystyle\mapsto\frac{4\int\mathrm{d}c\,\Theta_{a}(c)c_{b,t}^{*}c_{b^{\prime},t}}{\int\mathrm{d}c}$
(49) $\displaystyle\mathrm{Tr}[|a\rangle\langle a^{\prime}|(0)|b\rangle\langle
b^{\prime}|(t)]$
$\displaystyle\mapsto\frac{6\int\mathrm{d}c\,c_{a}^{*}c_{a^{\prime}}\,c_{b,t}^{*}c_{b^{\prime},t}}{\int\mathrm{d}c}$
(50)

where $W_{a}(c)=4|c_{a}|^{2}-2$ is a weight that goes to zero when the two
adiabats have equal populations. Note that this approach differs from the
others above not only in the initial distribution but also in the construction
of the time-dependent observable.

##  IV Results

Figure 2: Rate of intersite population transfer in a Frenkel-exciton dimer as
a function of the bath reorganization energy. MASH agrees closely with the
HEOM benchmark across the entire parameter range (left panel), regardless of
the particular choice of initial condition (right panel).

###  IV.1 Exciton dimer

To investigate the transition from Redfield to Förster-like transfer, we
consider a two-site exciton model with
$\epsilon_{1}-\epsilon_{2}=$100\text{\,}\mathrm{c}\mathrm{m}^{-1}$$ and
$J_{12}=$20\text{\,}\mathrm{c}\mathrm{m}^{-1}$$. Each site is coupled to a
bath at $T=$300\text{\,}\mathrm{K}$$ with the Debye spectral density

$J(\omega)=2\lambda\frac{\omega\omega_{\rm c}}{\omega^{2}+\omega^{2}_{\rm c}}$
(51)

where $\omega_{\rm c}=$53\text{\,}\mathrm{c}\mathrm{m}^{-1}$$. The quantity of
interest is the forward intersite rate $k_{1\to 2}$ as a function of
$\lambda$. This model is well-studied in the literatureIshizaki and Fleming
(2009a, b); Huo and Miller III (2015); Runeson (2022) and therefore allows
comparison with a wide range of methods. Quantum mechanical (HEOM) benchmark
results have been computed by Ishizaki and Fleming,Ishizaki and Fleming
(2009a) who observed that Redfield theory is accurate for small $\lambda$ but
qualitatively wrong for large $\lambda$, whereas Förster theory is only valid
for large $\lambda$. In their calculations, the forward and backward rate
constants were obtained by fitting the population dynamics to the kinetic
model

$\displaystyle\frac{d}{dt}\langle P_{1}(t)\rangle$ $\displaystyle=-k_{1\to
2}\langle P_{1}(t)\rangle+k_{2\to 1}\langle P_{2}(t)\rangle$
$\displaystyle\frac{d}{dt}\langle P_{2}(t)\rangle$
$\displaystyle=\phantom{+}k_{1\to 2}\langle P_{1}(t)\rangle-k_{2\to 1}\langle
P_{2}(t)\rangle$ (52)

where the site density matrix was initialized as $|1\rangle\langle 1|$. No
secular approximation was applied in the Redfield theory.

To assess the performance of MASH, we have calculated the population dynamics
using all three of the initial conditions considered in Section III, and
performed additional calculations with the original version of MASH for
comparison. The bath was discretized into 100 modes per site using a standard
discretization scheme Hele (2011) and the nuclei were initialized from the
classical Boltzmann distribution of an uncoupled bath. The dynamics was
averaged over $10^{5}$ trajectories.

When extracting the rate, we observed that fitting the population difference
$\langle\sigma_{z}(t)\rangle=\langle P_{1}(t)-P_{2}(t)\rangle$ to a single
exponential,

$\langle\sigma_{z}(t)\rangle=(\langle\sigma_{z}(0)\rangle-\langle\sigma_{z}\rangle_{\rm
eq})\mathrm{e}^{-k_{\rm tot}t}+\langle\sigma_{z}\rangle_{\rm eq},$ (53)

with $k_{\rm tot}=k_{1\to 2}+k_{2\to 1}$ as a fitting parameter, was more
stable than using the two-parameter fit in Eqs. (IV.1). Since the nuclear
statistics is essentially classical ($k_{\rm B}T<\omega_{\rm c}$), MASH is
guaranteed to recover the correct equilibrium value
$\langle\sigma_{z}\rangle_{\rm eq}=\langle P_{1}-P_{2}\rangle_{\rm eq}$, as do
HEOM, Redfield and Förster theory. So there is no need for an additional free
parameter. Once $k_{\rm tot}$ has been extracted from Eq. (53), the forward
rate constant can be calculated as

$k_{1\to 2}=k_{\rm tot}\langle P_{2}\rangle_{\rm eq}=k_{\rm
tot}\frac{1}{2}\left(1-\langle\sigma_{z}\rangle_{\rm eq}\right),$ (54)

which follows from the detailed balance relation $k_{1\to 2}/k_{2\to
1}=\langle P_{2}\rangle_{\rm eq}/\langle P_{1}\rangle_{\rm eq}$. To ensure a
fair comparison, we have also recalculated the Redfield and HEOM population
dynamics (using the Pyrho open source software packageBerkelbach (2020)), and
applied the same fitting procedure to those. This was found to lead to slighly
($<10$%) different rates compared to Ref. Ishizaki and Fleming, 2009a.

Figure 2 shows our results together with the Förster theory rates from Ref.
Ishizaki and Fleming, 2009a. We find that MASH agrees closely with HEOM across
the entire parameter range (left panel), including the Redfield and Förster-
type regimes. Moreover, all four versions of MASH lead to essentially the same
rates (right panel). This is interesting because the cap initial condition has
previously been found to be less accurate than the original version of MASH
for a spin-boson model in the Marcus inverted regime.Runeson and Manolopoulos
(2023) In the present calculations, the region
$2\lambda<\epsilon_{1}-\epsilon_{2}=$100\text{\,}\mathrm{c}\mathrm{m}^{-1}$$
is formally in the inverted regime, but the different initial conditions
nevertheless lead to similar behaviour. We reach the same conclusion when we
convert the Frenkel-exciton model into a spin-boson model with matching bias
and total reorganization energy, which we find does not noticeably change
either the HEOM or the MASH rates. Hence, all four versions of MASH can be
regarded as reliable in the present inverted regime. The inverted regime
considered in Ref. Runeson and Manolopoulos, 2023 was more challenging owing
to its larger bias (20 times the diabatic coupling matrix element rather than
the 5 times considered here), and in that regime the original version of MASH
is to be preferred.

Figure 3: Population dynamics in FMO at $300\text{\,}\mathrm{K}$ comparing
MASH with different initial conditions to well-established rate theories.
Left: dynamics in the site basis after an initial excitation of site 1. The
three MASH initial conditions lead to identical results and agree
qualitatively with Förster theory at long times. The inset shows the site
labels using the same colouring as for the data curves. Right: dynamics in the
exciton basis after an initial excitation of exciton 8. The three MASH initial
conditions lead to similar results and predict notably slower transfer than
(secular) Redfield theory. The inset depicts qualitatively the spatial extent
of the exciton states and their labels in order of increasing energy.

###  IV.2 Eight-site Fenna–Matthews–Olson complex

Another well-known benchmark system for exciton energy transfer is the
Fenna–Matthews–Olson complex found in green sulfur bacteria. We have
previously demonstrated that MASH (with cap initial conditions) agrees closely
with HEOM for a standard seven-site FMO model with a Debye spectral
density.Runeson and Manolopoulos (2023) Here, we consider a more challenging
(and realistic) eight-site model with a structured spectral density extracted
from fluorescence line narrowing experiments.Wendling  _et al._ (2000) The
resulting bath has a reorganization energy of
$\lambda=$45\text{\,}\mathrm{c}\mathrm{m}^{-1}$$. The intersite couplings and
average site energies (shown in Table 1) were obtained from electrostatic
calculations for the FMO complex of _Prosthecochloris aestuarii_.Schmidt am
Busch  _et al._ (2011) Static disorder was included by sampling the site
energies with the Gaussian widths shown in Table 2, which were calculated by
Müh _et al_.Müh  _et al._ (2007)

We are not aware of any fully quantum benchmarks for this model, so instead we
compare MASH to Förster and Redfield theory. These methods have a long history
in modelling the dynamics of FMO.Milder  _et al._ (2010); Cao  _et al._
(2020); Renger (2021) Here, we calculate the dynamics in the site basis using
Förster theory and the dynamics in the exciton basis using Redfield theory
within the secular approximation. In each basis, the dynamics is therefore
simply a propagation of the populations with a constant rate matrix. For
simplicity, we start from an excitation localized on a single site (for the
site basis calculation) or on a single exciton (for the exciton basis
calculation). The Förster and Redfield dynamics were averaged over 1000
samples of the site energies to account for static disorder. In the MASH
calculations, the bath was discretized into 100 modes per site using an
equally spaced grid up to $\omega_{\rm
max}=$500\text{\,}\mathrm{c}\mathrm{m}^{-1}$$, and the modes were initialized
from the classical Boltzmann distribution of an uncoupled bath at
$300\text{\,}\mathrm{K}$. The dynamics were averaged over $10^{6}$
trajectories for the cap and equivariant initial conditions and $10^{5}$ for
the focused initial condition to ensure tight convergence.

The left panel of Fig. 3 shows the site populations after an initial
excitation of site 1. All three MASH initial conditions give indistinguishable
results. Apart from a transient ($<$0.5\text{\,}\mathrm{p}\mathrm{s}$$)
coherence between sites 1 and 2, the dynamics is essentially rate-like.
Although Förster theory does not capture the coherence and differs from MASH
at short times, in particular for site 8, it agrees qualitatively with MASH at
longer times. This observation is consistent with previous studies for simpler
FMO models,Wu  _et al._ (2012); Wilkins (2015) where Förster theory was found
to be qualitatively reliable in comparison with exact benchmark results,
despite several site couplings being as strong as
$90\text{\,}\mathrm{c}\mathrm{m}^{-1}$. The reason is likely that the strong
couplings only matter for the first $\sim$100\text{\,}\mathrm{f}\mathrm{s}$$,
whereas on the $\sim$1\text{\,}\mathrm{p}\mathrm{s}$$ timescale the population
transfer is controlled by the weaker couplings for which Förster theory is
accurate.

The right panel of Fig. 3 shows the exciton populations after an initial
excitation of exciton 8. This exciton state is spatially located on sites 8
and 1, and has been identified as one of the dominant pathways when captured
photon energy enters the FMO complex from the baseplate of the chlorosomes.Cao
_et al._ (2020) Again, all MASH initial conditions lead to similar dynamics up
to a slight difference that washes out within
$1\text{\,}\mathrm{p}\mathrm{s}$. Notably, the overall transfer is
significantly slower than in Redfield theory, by roughly a factor of 2. This
observation is consistent with a previous study using a phase-space mapping of
the electronic states,Runeson  _et al._ (2022a) where it was shown that even
though the Markovian approximation is valid for the present bath, the
system–bath coupling is too large for second-order perturbative approaches
like Redfield theory to be reliable (see Figs. S2 and S5 of Ref. Runeson  _et
al._ , 2022a). Since MASH has the additional advantage of relaxing to the
correct long-time limit, we expect it to be more accurate than those previous
mapping calculations. Note, however, that MASH can experience negative
populations for intermediate times. In the present calculations, exciton state
7 becomes slightly negative between 0.1 and 0.3 ps with the ‘cap’ and
‘equivariant’ initial conditions. This could be a real effect or due to
insufficient sampling. Currently, neither of the versions of MASH guarantees
complete positivity of the system density matrix except in the long-time
limit. (For two states, the original MASH gives strictly non-negative
populations only in the adiabatic basis.) The quantum-jump correctionMannouch
and Richardson (2023); Lawrence, Mannouch, and Richardson (2023) may help to
alleviate this deficiency in future work.

Table 1: Average site energies and couplings for FMOSchmidt am Busch  _et al._
(2011) in units of ${\rm cm}^{-1}$.
Site | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
1 | 310 | $-94.8$ | $5.5$ | $-5.9$ | $7.1$ | $-15.1$ | $-12.2$ | $39.5$
2 |  | 230 | $29.8$ | $7.6$ | $1.6$ | $13.1$ | $5.7$ | $7.9$
3 |  |  | 0 | $-58.9$ | $-1.2$ | $-9.3$ | $3.4$ | $1.4$
4 |  |  |  | 180 | $-64.1$ | $-17.4$ | $-62.3$ | $-1.6$
5 |  |  |  |  | 405 | $89.5$ | $-4.6$ | $4.4$
6 |  |  |  |  |  | 320 | $35.1$ | $-9.1$
7 |  |  |  |  |  |  | 270 | $-11.1$
8 |  |  |  |  |  |  |  | 505
Table 2: Gaussian widths (full width at half maximum) of the site energiesMüh
_et al._ (2007) in ${\rm cm}^{-1}$.
Site | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
FWHM | 60 | 100 | 60 | 60 | 120 | 120 | 120 | 100

##  V Conclusions

In this article, we have shown by comparison with exact results that MASH
correctly captures the transition from the Redfield to Förster regimes for an
exciton dimer. This is the case no matter if one uses the original two-state
observables or the equivariant estimators in multi-state MASH. In conjuction
with the recent finding that MASH recovers Marcus theory in the diabatic
limit,Lawrence, Mannouch, and Richardson (2023) our results further establish
MASH as a generally reliable rate theory across several relevant parameter
regimes. Since it additionally relaxes to the correct equilibrium populations
for excitonic systems in classical environments, in contrast to any other
nonadiabatic dynamics method we are aware of, and since it is applicable to
systems described by general anharmonic interaction potentials, we would argue
that MASH is a practical tool that is capable of capturing almost all of the
relevant ingredients of exciton transfer. (It has yet to be generalized to
include quantum mechanical effects in the nuclear motion, which is a work in
progress.)

For a challenging model of FMO including static disorder and an experimental
spectral density, we find that MASH agrees qualitatively with Förster theory
(apart from a short transient coherence in the site basis), even though
several inter-site couplings are expected to be beyond the range of
applicability of the golden rule. In the exciton basis, MASH differs from
Redfield theory in its slower energy transfer timescale, confirming findings
from spin-mapping methodsRuneson and Richardson (2020); Runeson  _et al._
(2022b) that the system–bath coupling is too large to treat as a
perturbation.Runeson  _et al._ (2022a)

We have also described and resolved an important issue regarding the initial
conditions in multi-state MASH. For the present systems, we find that the
results are virtually identical for various different choices of the initial
conditions. Although the situation would likely be different for applications
in excited-state photochemistry, we conclude that for the condensed-phase
environments considered here one may use whichever initial condition is more
practical. A previous calculation for a spin-boson model in the Marcus
inverted regime has found that the ‘cap’ initial condition in multi-state MASH
and the original MASH method give different relaxation timescales,Runeson and
Manolopoulos (2023) but for the present dimer model there is no noticeable
difference between the two methods even for model parameters that correspond
to the inverted regime.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

J.E.R. was funded by a mobility fellowship from the Swiss National Science
Foundation. T.P.F. was supported by the U.S. Department of Energy, Office of
Science, Basic Energy Sciences, CPIMS Program Early Career Research Program
under Award DE-FOA0002019.

## References

  * Mirkovic  _et al._ (2017) T. Mirkovic, E. E. Ostroumov, J. M. Anna, R. van Grondelle, Govindjee, and G. D. Scholes, “Light absorption and energy transfer in the antenna complexes of photosynthetic organisms,” Chem. Rev. 117, 249–293 (2017).
  * Segatta  _et al._ (2019) F. Segatta, L. Cupellini, M. Garavelli, and B. Mennucci, “Quantum chemical modeling of the photoinduced activity of multichromophoric biosystems,” Chem. Rev. 119, 9361–9380 (2019).
  * May and Kühn (2011) V. May and O. Kühn,  _Charge and Energy Transfer Dynamics in Molecular Systems_ , 3rd ed. (Wiley-VCH, Weinheim, 2011).
  * Ishizaki  _et al._ (2010) A. Ishizaki, T. R. Calhoun, G. S. Schlau-Cohen, and G. R. Fleming, “Quantum coherence and its interplay with protein environments in photosynthetic electronic energy transfer,” Phys. Chem. Chem. Phys. 12, 7319–7337 (2010).
  * Alvertis  _et al._ (2020) A. M. Alvertis, W. Barford, S. Bourne Worster, I. Burghardt, A. Datta, A. Dijkstra, T. Fay, S. Ghosh, T. Grünbaum, S. Habershon, P. J. Hore, D. Hutchinson, S. Iyengar, A. R. Jones, G. Jones, K. Komarova, J. Lawrence, J. Léonard, Y. Litman, J. Mannouch, D. Manolopoulos, C. Martens, M. Mondelo-Martell, D. Picconi, D. Plant, K. Sakaushi, M. A. C. Saller, A. Schile, G. D. Scholes, J. Segarra-Martí, F. Segatta, A. Troisi, and G. Worth, “Quantum coherence in complex environments: general discussion,” Faraday Discuss. 221, 168–201 (2020).
  * Cao  _et al._ (2020) J. Cao, R. J. Cogdell, D. F. Coker, H.-G. Duan, J. Hauer, U. Kleinekathöfer, T. L. C. Jansen, T. Mančal, R. J. D. Miller, J. P. Ogilvie, V. I. Prokhorenko, T. Renger, H.-S. Tan, R. Tempelaar, M. Thorwart, E. Thyrhaug, S. Westenhoff, and D. Zigmantas, “Quantum biology revisited,” Sci. Adv. 6, eaaz4888 (2020).
  * Scholes (2018) G. D. Scholes, “Coherence from light harvesting to chemistry,” J. Phys. Chem. Lett. 9, 1568–1572 (2018).
  * Sneyd, Beljonne, and Rao (2022) A. J. Sneyd, D. Beljonne, and A. Rao, “A new frontier in exciton transport: Transient delocalization,” J. Phys. Chem. Lett. 13, 6820–6830 (2022).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath,” J. Phys. Soc. Jpn. 58, 101–114 (1989).
  * Ishizaki and Fleming (2009a) A. Ishizaki and G. R. Fleming, “Unified treatment of quantum coherent and incoherent hopping dynamics in electronic energy transfer: Reduced hierarchy equation approach,” J. Chem. Phys. 130, 234111 (2009a).
  * Makri and Makarov (1995) N. Makri and D. E. Makarov, “Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory,” J. Chem. Phys. 102, 4600–4610 (1995).
  * Makri (2020) N. Makri, “Small matrix disentanglement of the path integral: Overcoming the exponential tensor scaling with memory length,” J. Chem. Phys. 152, 041104 (2020).
  * Tully (1990) J. C. Tully, “Molecular dynamics with electronic transitions,” J. Chem. Phys. 93, 1061–1071 (1990).
  * Wang  _et al._ (2020) L. Wang, J. Qiu, X. Bai, and J. Xu, “Surface hopping methods for nonadiabatic dynamics in extended systems,” WIREs Comput. Mol. Sci. 10, e1435 (2020).
  * Mannouch and Richardson (2023) J. R. Mannouch and J. O. Richardson, “A mapping approach to surface hopping,” J. Chem. Phys. 158, 104111 (2023).
  * Lawrence, Mannouch, and Richardson (2023) J. E. Lawrence, J. R. Mannouch, and J. O. Richardson, “Recovering Marcus theory rates and beyond without the need for decoherence corrections: The mapping approach to surface hopping,” (2023), arXiv:2311.08802 [physics.chem-ph].
  * Amati, Mannouch, and Richardson (2023) G. Amati, J. R. Mannouch, and J. O. Richardson, “Detailed balance in mixed quantum-classical mapping approaches,” (2023), arXiv:2309.04686 [quant-ph].
  * Runeson and Manolopoulos (2023) J. E. Runeson and D. E. Manolopoulos, “A multi-state mapping approach to surface hopping,” J. Chem. Phys. 159, 094115 (2023).
  * Yang and Fleming (2002) M. Yang and G. R. Fleming, “Influence of phonons on exciton transfer dynamics: comparison of the Redfield, Förster, and modified Redfield equations,” Chem. Phys. 275, 355–372 (2002).
  * Folland (2001) G. B. Folland, “How to integrate a polynomial over a sphere,” Amer. Math. Monthly 108, 446–448 (2001).
  * Saller, Kelly, and Richardson (2019) M. A. C. Saller, A. Kelly, and J. O. Richardson, “On the identity of the identity operator in nonadiabatic linearized semiclassical dynamics,” J. Chem. Phys. 150, 071101 (2019).
  * Gao  _et al._ (2020) X. Gao, M. A. C. Saller, Y. Liu, A. Kelly, J. O. Richardson, and E. Geva, “Benchmarking quasiclassical mapping Hamiltonian methods for simulating electronically nonadiabatic molecular dynamics,” J. Chem. Theory Comput. 16, 2883–2895 (2020).
  * Ishizaki and Fleming (2009b) A. Ishizaki and G. R. Fleming, “On the adequacy of the redfield equation and related approaches to the study of quantum dynamics in electronic energy transfer,” J. Chem. Phys. 130, 234110 (2009b).
  * Huo and Miller III (2015) P. Huo and T. F. Miller III, “Electronic coherence and the kinetics of inter-complex energy transfer in light-harvesting systems,” Phys. Chem. Chem. Phys. 17, 30914–30924 (2015).
  * Runeson (2022) J. E. Runeson,  _Spin-mapping approaches for mixed quantum-classical dynamics_ , Ph.D. thesis, ETH Zurich (2022).
  * Hele (2011) T. J. H. Hele,  _An Electronically Non-Adiabatic Generalization of Ring Polymer Molecular Dynamics_ , Master’s thesis, University of Oxford (2011).
  * Berkelbach (2020) T. C. Berkelbach, “Pyrho: A Python package for reduced density matrix techniques,”  (2020).
  * Wendling  _et al._ (2000) M. Wendling, T. Pullerits, M. A. Przyjalgowski, S. I. Vulto, T. J. Aartsma, R. van Grondelle, and H. van Amerongen, “Electron–vibrational coupling in the Fenna–Matthews–Olson complex of _Prosthecochloris aestuarii_ determined by temperature-dependent absorption and fluorescence line-narrowing measurements,” J. Phys. Chem. B 104, 5825–5831 (2000).
  * Schmidt am Busch  _et al._ (2011) M. Schmidt am Busch, F. Müh, M. El-Amine Madjet, and T. Renger, “The eighth bacteriochlorophyll completes the excitation energy funnel in the FMO protein,” J. Phys. Chem. Lett. 2, 93–98 (2011).
  * Müh  _et al._ (2007) F. Müh, M. E.-A. Madjet, J. Adolphs, A. Abdurahman, B. Rabenstein, H. Ishikita, E.-W. Knapp, and T. Renger, “$\alpha$-Helices direct excitation energy flow in the Fenna–Matthews–Olson protein,” P. Natl. Acad. Sci. USA 104, 16862–16867 (2007).
  * Milder  _et al._ (2010) M. T. W. Milder, B. Brüggemann, R. van Grondelle, and J. L. Herek, “Revisiting the optical properties of the FMO protein,” Photosynth. Res. 104, 257–274 (2010).
  * Renger (2021) T. Renger, “Semiclassical modified Redfield and generalized Förster theories of exciton relaxation/transfer in light-harvesting complexes: The quest for the principle of detailed balance,” J. Phys. Chem. B 125, 6406–6416 (2021).
  * Wu  _et al._ (2012) J. Wu, F. Liu, J. Ma, R. J. Silbey, and J. Cao, “Efficient energy transfer in light-harvesting systems: Quantum-classical comparison, flux network, and robustness analysis,” J. Chem. Phys. 137, 174111 (2012).
  * Wilkins (2015) D. M. Wilkins,  _A Theoretical Investigation Into Energy Transfer In Photosynthetic Open Quantum Systems_ , Master’s thesis, University of Oxford (2015).
  * Runeson  _et al._ (2022a) J. E. Runeson, J. E. Lawrence, J. R. Mannouch, and J. O. Richardson, “Explaining the efficiency of photosynthesis: Quantum uncertainty or classical vibrations?” J. Phys. Chem. Lett. 13, 3392–3399 (2022a).
  * Runeson and Richardson (2020) J. E. Runeson and J. O. Richardson, “Generalized spin mapping for quantum-classical dynamics,” J. Chem. Phys. 152, 084110 (2020).
  * Runeson  _et al._ (2022b) J. E. Runeson, J. R. Mannouch, G. Amati, M. R. Fiechter, and J. O. Richardson, “Spin-mapping methods for simulating ultrafast nonadiabatic dynamics,” Chimia 76, 582 (2022b).