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

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.^1,2 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.^3,4 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.^5,6

The main challenge is that many systems of interest tend to appear in the intermediate regime where the two timescales are comparable.^7,8 For simple system-bath models one can solve the dynamics non-perturbatively using fully quantum methods such as the hierarchical equations of motion (HEOM).⁹ 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.¹⁰ However, fully quantum methods like HEOM and quasi-adiabatic path integrals^11,12 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).¹³ 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’.¹⁴

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.¹⁵ 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.¹⁵ 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.¹⁶ 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.¹⁷

In its original formulation, MASH was limited to systems with two electronic states, but an adaptation to multiple states has recently been proposed.¹⁸ 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.¹⁸

Unfortunately, the equivariant multi-state formulation of MASH¹⁸ does not reduce to the original formulation¹⁵ 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.¹⁸ 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.¹⁸ 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. 18 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_S + H_B + H_SB. (1)

The first term is the system Hamiltonian in the basis of localized pigment (‘site’) excitationswhere 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,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 densityand 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 energyThe system-bath interaction involves the bath operator , 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) = Tr_B[e^−βH_BB(0)B(t)], (7)

where β = 1/(k_BT) and the time-dependence refers to dynamics under H_B. Using known expressions for the thermal correlation functions of a harmonic oscillator, one can show thatwhere n(ω) = 1/(e^βω − 1) is the Bose–Einstein distribution and throughout we use units where ħ = 1.

A. 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_S, which is called the exciton basis. Diagonalizing the system Hamiltonian gives where ω_μ is an eigenenergy and 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 derivation³ leads to a Markovian master equation for the reduced density matrix of the system,where ω_μν = ω_μ − ω_ν is the energy gap between the corresponding eigenstates of H_S. The second term involves the Redfield tensorwhich is expressed in terms of the damping tensorHere, C̃(ω) is the Fourier-Laplace transform of the bath correlation function,Inserting the expression for C(t) from eqn (8) gives

Re C̃(ω) = J(ω)(1 + n(ω)) + J(−ω)n(−ω) (14)

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

B. 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 system follows a kinetic rate equation in the site basis. The Förster rate constants are¹⁹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 asThe function g(t) iswhere C(t) is the bath correlation function in eqn (8). For the Frenkel-exciton model with identical baths, the rate expression in eqn (18) reduces towhere

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 eqn (1) aswhereis a (diabatic) potential operator with matrix elementsSurface hopping is almost always run on adiabatic surfaces, i.e., in the eigenbasis of V̂(q),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 |ψ〉, which can be expanded in the diabatic or the adiabatic basis asThe coefficients obey the Schrödinger equationor (equivalently)where d^j_ab(q) = 〈a(q)|∇_j|b(q)〉 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

q̇_j = p_j/m_j (31a)

ṗ_j = −〈a(q)|∇_jV̂|a(q)〉 (31b)

where a is the adiabatic state with the largest instantaneous population |c_a|². Effectively, this means that the nuclei evolve on a potential with abrupt steps. By introducing the classical state projectors,one can express the effective potential perceived by the nuclei asAt each instant, precisely one of the Θ_a(c) factors in the sum of eqn (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 elements¹⁸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).

A. 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〉〈i| with a (normalized) nuclear density ρ_B(p,q). The time-dependent expectation value of an observable O is then

〈O(t)〉 = Tr[ρ_B(p,q)|i〉〈i|(0)Ô(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 integralwhere 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_n=1 and set .

In the following, we consider the case where Ô_n = |n〉〈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.¹⁸ If we are interested in an adiabatic population, then the state projector Θ_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.¹⁸

Another estimator that would transform correctly between bases is the Ehrenfest population |c_n|², but inserting this choice into eqn (36) leads to the wrong long-time equilibrium populations. In ref. 18 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) = α_N|c_n|² + β_N (37)

whereare 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.¹⁸

B. Initial conditions

What remains to be defined is the choice of initial electronic distribution ρ_i(c) in eqn (36). This distribution is not unique, in the sense that many (quasi)probability distributions will fulfil the initial conditionIn the following we consider a few options.

1. Cap initial condition. A simple choice is

ρ_i(c) = Θ_i(c), (40)

which was used in ref. 18. What this equation implies is that the initial state is chosen randomly from the region where |c_i|² 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 eqn (40) as the ‘cap’ initial distribution.

Fig. 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, , , and the vertical axis σ_z = |c₁|² − |c₂|² 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 eqn (37). The original MASH method uses an alternative prescription that is (so far) limited to two states.

Although the cap distribution has been found to be accurate in a variety of benchmark calculations,¹⁸ it is not basis-equivariant, unlike the population estimator in eqn (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 where is some number in the range . This vector has Θ₁(c) = 1 and would therefore contribute to the dynamics. But if we define a new set of diabatic basis vectors |1̃〉 = |1〉, , then for the maximally populated state is |2̃〉. Thus, even though the ket of the initial state is unchanged, the same vector c would no longer contribute to the dynamics.

2. Focused initial condition. Perhaps the most intuitive choice for the initial wavefunction corresponding to |i〉〈i| would be to set c_i = 1 and c_j≠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 eqn (39). The reason is that in the initial corner of the simplex, O_i(c) > 1 and O_j≠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≠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|² = (1 − β_N)/α_N and |c_j≠i|² = −β_N/α_N. These conditions fix the magnitudes of all components of c, leaving the phases to be sampled uniformly from [0,2π). The resulting ‘focused’ initial distribution can be written as

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.

3. Equivariant initial condition. In this section, we derive a quasiprobability distribution ρ_i that transforms correctly under unitary basis transformations. To satisfy the condition in eqn (39), the simplest approach is to try the same functional form as the time-dependent observable in eqn (37), i.e.

ρ_i(c) = a_N|c_i|² + b_N (42)

with some constants a_N, b_N that need not be the same as α_N, β_N. The resulting ρ_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 eqn (39), we make use of the following moments:

whereThese expectation values can be derived using standard formulas for integrals over a sphere.²⁰ With the help of the moments in eqn (43), eqn (39) reduces to two equations in two unknowns, with the solutionNote that there is some similarity between the Roman constants and the Greek ones in eqn (38). For example, α_N and β_N are related through Nβ_N = 1 − α_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 ρ_i(c) involves a scaling relative to the centre of the simplex: after inserting β_N and b_N we getand

The special status of the centre of the simplex is analogous to the special role of the identity operator in phase-space mapping methods.^21,22

4. Original MASH. For two-level systems, the original MASH method of Mannouch and Richardson¹⁵ 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 ≠ a′ and b ≠ b′, their prescription is (in our notation)where W_a(c) = 4|c_a|² − 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

A. Exciton dimer

To investigate the transition from Redfield to Förster-like transfer, we consider a two-site exciton model with and J₁₂ = 20 cm⁻¹. Each site is coupled to a bath at T = 300 K with the Debye spectral densitywhere ω_c = 53 cm⁻¹. The quantity of interest is the forward intersite rate k_1→2 as a function of λ. This model is well-studied in the literature^10,23–25 and therefore allows comparison with a wide range of methods. Quantum mechanical (HEOM) benchmark results have been computed by Ishizaki and Fleming,¹⁰ who observed that Redfield theory is accurate for small λ but qualitatively wrong for large λ, whereas Förster theory is only valid for large λ. In their calculations, the forward and backward rate constants were obtained by fitting the population dynamics to the kinetic modelwhere the site density matrix was initialized as |1〉〈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 3, 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²⁶ and the nuclei were initialized from the classical Boltzmann distribution of an uncoupled bath. The dynamics was averaged over 10⁵ trajectories.

When extracting the rate, we observed that fitting the population difference 〈σ_z(t)〉 = 〈P₁(t) − P₂(t)〉 to a single exponential,

〈σ_z(t)〉 = (〈σ_z(0)〉 − 〈σ_z〉_eq)e^−k_tott + 〈σ_z〉_eq, (53)

with k_tot = k_1→2 + k_2→1 as a fitting parameter, was more stable than using the two-parameter fit in eqn (52). Since the nuclear statistics is essentially classical (k_BT < ω_c), MASH is guaranteed to recover the correct equilibrium value 〈σ_z〉_eq = 〈P₁ − P₂〉_eq, as do HEOM, Redfield and Förster theory. So there is no need for an additional free parameter. Once k_tot has been extracted from eqn (53), the forward rate constant can be calculated aswhich follows from the detailed balance relation k_1→2/k_2→1 = 〈P₂〉_eq/〈P₁〉_eq. To ensure a fair comparison, we have also recalculated the Redfield and HEOM population dynamics (using the Pyrho open source software package²⁷), and applied the same fitting procedure to those. This was found to lead to slighly (<10%) different rates compared to ref. 10.

Fig. 2 shows our results together with the Förster theory rates from ref. 10. 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.¹⁸ In the present calculations, the region 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. 18 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.

Fig. 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).

B. 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.¹⁸ Here, we consider a more challenging (and realistic) eight-site model with a structured spectral density extracted from fluorescence line narrowing experiments.²⁸ The resulting bath has a reorganization energy of λ = 45 cm⁻¹. The intersite couplings and average site energies (shown in Table 1) were obtained from electrostatic calculations for the FMO complex of Prosthecochloris aestuarii.²⁹ 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.³⁰

Table 1 Average site energies and couplings for FMO²⁹ in units of cm⁻¹

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 energies³⁰ in cm⁻¹

Site	1	2	3	4	5	6	7	8
FWHM	60	100	60	60	120	120	120	100

---

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.^6,31,32 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 ω_max = 500 cm⁻¹, and the modes were initialized from the classical Boltzmann distribution of an uncoupled bath at 300 K. The dynamics were averaged over 10⁶ trajectories for the cap and equivariant initial conditions and 10⁵ 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 ps) 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,^33,34 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 cm⁻¹. The reason is likely that the strong couplings only matter for the first ∼100 fs, whereas on the ∼1 ps timescale the population transfer is controlled by the weaker couplings for which Förster theory is accurate.

Fig. 3 Population dynamics in FMO at 300 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.

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.⁶ Again, all MASH initial conditions lead to similar dynamics up to a slight difference that washes out within 1 ps. 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,³⁵ 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 Fig. S2 and S5 of ref. 35). 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 correction^15,16 may help to alleviate this deficiency in future work.

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,¹⁶ 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 methods^36,37 that the system-bath coupling is too large to treat as a perturbation.³⁵

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,¹⁸ 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 under Award P500PN_206641/1. 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

1. T. Mirkovic, E. E. Ostroumov, J. M. Anna, R. van Grondelle, Govindjee and G. D. Scholes, Chem. Rev., 2017, 117, 249–293.
2. F. Segatta, L. Cupellini, M. Garavelli and B. Mennucci, Chem. Rev., 2019, 119, 9361–9380.
3. V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, Wiley-VCH, Weinheim, 3rd edn, 2011.
4. A. Ishizaki, T. R. Calhoun, G. S. Schlau-Cohen and G. R. Fleming, Phys. Chem. Chem. Phys., 2010, 12, 7319–7337.
5. 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, Faraday Discuss., 2020, 221, 168–201.
6. J. Cao, R. J. Cogdell, D. F. Coker, H.-G. Duan, J. Hauer, U. Kleinekathöfer, T. L. C. Jansen, T. Mancal, 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, Sci. Adv., 2020, 6, eaaz4888.
7. G. D. Scholes, J. Phys. Chem. Lett., 2018, 9, 1568–1572.
8. A. J. Sneyd, D. Beljonne and A. Rao, J. Phys. Chem. Lett., 2022, 13, 6820–6830.
9. Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn., 1989, 58, 101–114.
10. A. Ishizaki and G. R. Fleming, J. Chem. Phys., 2009, 130, 234111.
11. N. Makri and D. E. Makarov, J. Chem. Phys., 1995, 102, 4600–4610.
12. N. Makri, J. Chem. Phys., 2020, 152, 041104.
13. J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071.
14. L. Wang, J. Qiu, X. Bai and J. Xu, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2020, 10, e1435.
15. J. R. Mannouch and J. O. Richardson, J. Chem. Phys., 2023, 158, 104111.
16. 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].
17. G. Amati, J. R. Mannouch and J. O. Richardson, Detailed balance in mixed quantum-classical mapping approaches, 2023, arXiv:2309.04686 [quant-ph].
18. J. E. Runeson and D. E. Manolopoulos, J. Chem. Phys., 2023, 159, 094115.
19. M. Yang and G. R. Fleming, Chem. Phys., 2002, 275, 355–372.
20. G. B. Folland, Amer. Math. Monthly, 2001, 108, 446–448.
21. M. A. C. Saller, A. Kelly and J. O. Richardson, J. Chem. Phys., 2019, 150, 071101.
22. X. Gao, M. A. C. Saller, Y. Liu, A. Kelly, J. O. Richardson and E. Geva, J. Chem. Theory Comput., 2020, 16, 2883–2895.
23. A. Ishizaki and G. R. Fleming, J. Chem. Phys., 2009, 130, 234110.
24. P. Huo and T. F. Miller III, Phys. Chem. Chem. Phys., 2015, 17, 30914–30924.
25. J. E. Runeson, PhD thesis, ETH Zurich, 2022.
26. T. J. H. Hele, MSc thesis, University of Oxford, 2011, arXiv:1308.3950 [physics.chem-ph].
27. T. C. Berkelbach, Pyrho: A Python package for reduced density matrix techniques, 2020, https://github.com/berkelbach-group/pyrho.
28. M. Wendling, T. Pullerits, M. A. Przyjalgowski, S. I. Vulto, T. J. Aartsma, R. van Grondelle and H. van Amerongen, J. Phys. Chem. B, 2000, 104, 5825–5831.
29. M. Schmidt am Busch, F. Müh, M. El-Amine Madjet and T. Renger, J. Phys. Chem. Lett., 2011, 2, 93–98.
30. F. Müh, M. E.-A. Madjet, J. Adolphs, A. Abdurahman, B. Rabenstein, H. Ishikita, E.-W. Knapp and T. Renger, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 16862–16867.
31. M. T. W. Milder, B. Brüggemann, R. van Grondelle and J. L. Herek, Photosynth. Res., 2010, 104, 257–274.
32. T. Renger, J. Phys. Chem. B, 2021, 125, 6406–6416.
33. J. Wu, F. Liu, J. Ma, R. J. Silbey and J. Cao, J. Chem. Phys., 2012, 137, 174111.
34. D. M. Wilkins, MSc thesis, University of Oxford, 2015, arXiv:1503.03277 [physics.chem-ph].
35. J. E. Runeson, J. E. Lawrence, J. R. Mannouch and J. O. Richardson, J. Phys. Chem. Lett., 2022, 13, 3392–3399.
36. J. E. Runeson and J. O. Richardson, J. Chem. Phys., 2020, 152, 084110.
37. J. E. Runeson, J. R. Mannouch, G. Amati, M. R. Fiechter and J. O. Richardson, Chimia, 2022, 76, 582.

---