A stochastic reorganizational bath model for electronic energy transfer

The fluctuations of optical gap induced by the environment play crucial roles in electronic energy transfer dynamics. One of the simplest approaches to incorporate such fluctuations in energy transfer dynamics is the well known Haken-Strobl-Reineker model, in which the energy-gap fluctuation is approximated as a white noise. Recently, several groups have employed molecular dynamics simulations and excited-state calculations in conjunction to take the thermal fluctuation of excitation energies into account. Here, we discuss a rigorous connection between the stochastic and the atomistic bath models. If the phonon bath is treated classically, time evolution of the exciton-phonon system can be described by Ehrenfest dynamics. To establish the relationship between the stochastic and atomistic bath models, we employ a projection operator technique to derive the generalized Langevin equations for the energy-gap fluctuations. The stochastic bath model can be obtained as an approximation of the atomistic Ehrenfest equations via the generalized Langevin approach. Based on the connection, we propose a novel scheme to correct reorganization effects within the framework of stochastic models. The proposed scheme provides a better description of the population dynamics especially in the regime of strong exciton-phonon coupling. Finally, we discuss the effect of the bath reorganization in the absorption and fluorescence spectra of ideal J-aggregates in terms of the Stokes shifts. For this purpose, we introduce a simple relationship that relates the reorganization contribution to the Stokes shifts - the reorganization shift - to three parameters: the monomer reorganization energy, the relaxation time of the optical gap, and the exciton delocalization length. This simple relationship allows one to classify the origin of the Stokes shifts in molecular aggregates.


I. INTRODUCTION
Natural photosynthesis starts with light energy absorption by an assembly of photosynthetic pigments, after which this excitation energy is transferred to a reaction center, [1][2][3][4][5][6][7][8][9] where charge transfer is carried out. In most organisms, various pigment-protein complexes are responsible for this light harvesting. The highly ordered structures of light-harvesting complexes have motivated researchers to design artificial light-harvesting antenna systems such as self-assembled supramolecular systems, [10][11][12] quantum dots, 13,14 and metalorganic frameworks. 15,16 Experimental studies of natural and artificial lightharvesting systems have spurred theoretical descriptions of electronic energy transfer (EET) in complex systems. EET has been modeled using Förster theory [17][18][19] which describes exciton transport as incoherent hopping between chromophores. However, this approach is applicable only in the strong exciton-phonon coupling regime, where the excitonic couplings between chromophores are small relative to the exciton-phonon couplings. In the opposite limit, a quantum master equation can be derived by treating the excitonphonon coupling perturbatively. 20 The most commonly a) Electronic addresses: tfujita@fas.harvard.edu and aspuru@chemistry.harvard.edu used theory from this limit is the Redfield approach. 21 In order to more accurately model EET in the intermediate regime between these two limiting cases, one could consider using the hierarchy equations of motions (HEOM), [22][23][24][25][26][27][28][29][30] path integral Monte Carlo, 31, 32 a polaron-transformed master equation, [33][34][35] non-Markovian quantum state diffusion, [36][37][38] or density matrix renormalization group methods. 39 Different methodologies for simulating energy transfer have been widely reviewed. [40][41][42][43][44][45][46][47] The most expensive and accurate methods are unified models that combine electronic structure directly with bath models. 48 The details of exciton-phonon coupling are crucial in light-harvesting systems because they strongly modify the excited-state dynamics. The environmental effect is usually characterized as a two-time bath correlation function or a bath spectral density. The Haken-Strobl-Reinker (HSR) 49, 50 model approximates bath fluctuation as white noise. The model can be extended to treat colored noise. 51,52 However, HSR-like approaches are incapable of describing bath reorganization and finite temperature effects. A lack of a bath reorganization process causes the dephasing rate to be underestimated 20 and leads to the same peak positions for absorption and fluorescence spectra without the Stokes shift. As mentioned previously, in recent years, exciton-phonon interactions have been treated at the atomistic level by combining molecular dynamics simulations and excited-state calculations. [53][54][55][56][57][58] From these simulations, one can perform ensemble-averaged wavepacket dynamics, 29, 59-61 also called Ehrenfest dynamics, or extract a spectral density 62 for its use in a quantum master equation.
In earlier work on the Fenna-Matthes-Olson complex, 57 we compared the atomistic Ehrenfest dynamics -the atomistic bath model -with the stochastic HSR model. Surprisingly, these two methods gave similar exciton population dynamics regardless of the different descriptions of the phonon baths. Motivated by these results, we discuss a rigorous connection between the atomistic and stochastic bath models and develop a novel correction scheme for the stochastic approaches. Atomistic or parameterized Ehrenfest dynamics require the time evolution of many bath degrees of freedom. The advantage of the stochastic approaches lies in the use of collective bath variables. This computational simplicity of the stochastic bath model will be useful to simulate EET in large systems, for example, biological scales. To relate the atomistic and stochastic bath models, we exploit the projection operator technique 63 to derive the generalized Langevin equations for energy-gap fluctuations. The stochastic model can be derived directly from the atomistic Ehrenfest equations via the generalized Langevin approach. Based on the connection between the stochastic and the atomic approaches, we propose a novel scheme to correct a reorganization effect in the stochastic model. Finally, our scheme will be compared with the HEOM through the population dynamics of a model system.
The remainder of this paper is organized as follows: In Sec. II, we introduce the exciton-phonon Hamiltonian and the Ehrenfest equations. Next, we exploit the projection operator technique in the atomistic Ehrenfest equations to derive stochastic bath models. Furthermore, we propose a scheme to correct reorganization effects within the stochastic bath models. In Sec. III, we present numerical results for the population dynamics and Stokes shifts in idealized J-aggregates. In Sec. IV, we present the paper with concluding remarks.

A. Excitonic Hamiltonian and Ehrenfest equations
The excitonic Hamiltonian in the single exciton manifold of a molecular aggregate can be written as follows: 20 where |m denotes the state where an electronic excitation is localized at mth molecule (site) and all other molecules are in the ground states. R and P refer to the nuclear coordinates and momentum, respectively. m (R) represents the excitation energy of mth site in the nuclear configuration R, and V mn (R) is an excitonic coupling constant between mth and nth molecules. Here, T (P) and V G (R) are the kinetic energy and the ground-state potential energy for the nuclear coordinates. Hereafter, the dependence of V mn on R is neglected. Then, we decompose the total Hamiltonian into the form of the system-bath Hamiltonian as follows: Here, m (R) = m (R) − m , and a bracket denotes the ensemble average over H B .
If we treat the nuclear degree of freedoms classically, the wave function of the exciton system is described with the time-dependent Schrödinger equation, The nuclear degrees of freedom follow from Hamilton's equation in the mean-field interaction, with the Hamiltonian for the nuclei on the excited-state potential energy surface (PES) In Eq. (8), H (t) gives the Ehrenfest mean potential; that is, the average of site energies weighted by exciton populations, | m|ψ(t) | 2 . This term shifts the excited-state PES with respect to the ground-state one, as illustrated in Fig. 1. The PES shift is essential to describe reorganization processes after the photoexcitation. This effect is also referred to as a back reaction, 29 because through this term the dynamics of the exciton system can affect the dynamics of the phonon bath.
The density matrix of the excitonic system ρ is obtained as the average of an ensemble of unitary evolutions FIG. 1. Schematic of H (t), which shifts the excited-state PES compared with the ground-state one. V E refers to the excited-state PES: Equations (5)-(10) form the basis for mixed quantumclassical approaches for exciton dynamics. [53][54][55][56][57] In these applications, the Ehrenfest mean force (i.e., the second term in Eq. (8)) is neglected, because one needs to run MD simulations in excited-state PES, which is computationally more expensive to determine. This term may be ignored when the reorganization energy is small relative to the excitonic couplings. Later, we discuss the condition under which this term can be ignored.

B. Generalized Langevin equations for energy-gap fluctuations
In this section, we derive stochastic bath models from the atomistic Ehrenfest equations. To relate the atomistic and stochastic bath models, we apply the projection operator technique [63][64][65] first developed in classical statistical mechanics. The projection operator techniques can provide a microscopic derivation of a phenomenological equation such as the Langevin equation. We briefly review Mori's projection operator formalism 63 and the derivation of the generalized Langevin equations. The standard technique is applied to the exciton-phonon system to derive stochastic equations for the energy-gap fluctuations.
Here, we consider the time evolution of the site energy fluctuations, In classical mechanics, any physical variable A is a function of R and P, and the time evolution of A is given by a Poisson bracket with the Hamiltonian. Accordingly, the equations of motion of (t) are given by Here, we introduce the Liouville operators of H B and H (t) as follows: First, we neglect the iL term in such a way that the time evolution of can be determined solely from iL 0 . This approximation leads to the exciton dynamics ignoring the reorganization effects. We then apply the projection operator technique 63 to derive the generalized Langevin equations for . We define the projection operator P from the site-energy fluctuations In this paper, we will use a simplified notation (A, B) ≡ BA* . By applying the projection operator to Eq. (12), the time evolution for (t) is obtained as the generalized Langevin Here, the memory matrix M(t) and the random force vector F(t) are defined as where Q = 1 − P.
Given the generalized Langevin equations for the siteenergy fluctuations, deriving stochastic models for them is straightforward. 49,50,66,67 We introduce the Markov approximation to the memory functions and neglect cross correlations between different sites as follows: By using this approximation, Eq. (16) becomes a set of Langevin equations with Gaussian fluctuations F m (t), where The inverse of m is related to the decay time 1/ m = τ m , which characterizes the relaxation time of the bath fluctuations. As the correlation function induced by Eq. (20) is an exponentially decaying function, i.e., this method is identical to the Kubo-Anderson (KA) stochastic model. 66,67 Furthermore, if we consider the limit of τ →0, the correlation function becomes a delta function So far, we have discussed that the stochastic bath models, such as the KA or HSR models, are the limiting case of the atomistic Ehrenfest equations. The stochastic bath models are known to result in high-temperature dynamics and are incapable of describing reorganization processes. The hightemperature dynamics results from the classical treatment of the phonon bath. Classical correlation functions do not satisfy the detailed balance condition. Within the stochastic Schrödinger equations, the non-Markovian terms and the detailed balance condition are important for describing the thermal relaxation process. 44 The lack of reorganization is due to the neglect of the H (t) term, which plays an essential role in the reorganization process that follows photoexcitation. The parameters m and 2 m for the dynamical equations can be obtained from MD simulations and excited-state calculations. The proposed simplified descriptions of the bath fluctuations allow us to simulate EET in large systems. For example, in earlier works, we used the KA and HSR approaches to simulate the light-harvesting apparatus of green-sulfur bacteria, which consists of thousands of chromophores. 58,68,69

C. Reorganization correction to stochastic bath model
In this section, we modify the KA model such that the reorganization effect can be included via the Ehrenfest mean potential. We begin by the exact time evolution of the site-energy fluctuations in the excited-state manifold. The dynamics of in the excited-state can be given by the time-evolution op- The definition of m has not been changed, i.e., m = − g . Here, g and e refer to the ensemble average over ground-and excited-state PESs, respectively. The first-order expansion of the time-evolution operator with respect to iL (t) gives The first term in the right-hand side is identical to Eq. (16), i.e., the time evolution in the ground-state. The second term in the right-hand side is the first-order correction of the reorganization effects. We redefine the time-local second term by K (t) that is a characteristic force induced by H (t). Finally, we can obtain the following equation for : Within the first-order perturbation with respect to iL , we can correct the reorganization effects by adding the time-local term to the Langevin equations. After applying the same approximation as Eq. (19), the time evolution of the in the excited-state follows the Langevin equation with the white noise, The formal solution of Eq. (28) is (28) After taking the ensemble average of Eq. (28) and using the definition of and F(t), we have Here, we consider a situation where the excitonic system is excited at t = 0, resulting in the same nuclear configuration which implies that m (0) e = m g . Intuitively, K m (t) describes the Stokes shifts and thus the reorganization process.
In the following, we use the linear-response theory to obtain K m (t). The first-order perturbation of site energies due to H (t) 70 is where P n (t) = n|ρ|n , and the inverse temperature is denoted by β = 1/(k B T). We approximate the correlation function as an exponential (Eq. (22)) in such a way that it coincides with the previous derivations. Thus, the difference in average site energies is given by This is the main result leading to a new term in the stochastic equations and the description of the reorganization effects. The steady-state solution of Eq. (27) is m e − m g = −β 2 m P m . Note that, in the high-temperature limit, the variance and reorganization energy λ are related by 2 g = 2k B Tλ. If a system is composed of a single pigment (P=1), then it is easy to see that one gets e − g = −2λ. This is a well-known result that the Stokes shift is twice the reorganization energy. As expected, K m (t) shifts the excited-state PES and induces the Stokes shifts. In Fig. 2, we use Eqs. (20) and (27) to show the distribution of site energies for the groundand excited-state of a two-level molecule. As expected, the distributions are shifted by 2λ.
By solving Eq. (5) with Eq. (27), we have shown that the reorganization effects can be incorporated within the framework of a stochastic bath model. For the remainder of this paper, we refer to our scheme as the reorganization-corrected Kubo-Anderson (RECKA) model. The exciton populations P m are necessary for input into RECKA, and they can be obtained after taking an ensemble average of unitary evolutions. Therefore, Eqs. (5), (10) and (27) must be solved selfconsistently until P m converge.
We note that the newly proposed term, K m (t), can be determined by the variance and the inverse of relaxation times, so additional computational time is not required. The lack of additional computational expense is a huge advantage of the RECKA over others methods that evaluate Ehrenfest potentials accurately; these methods demand the derivative calculation of the excited-state PES, which is a computationally demanding task that limits such approaches to small systems. 71 Noticeably, the low cost of employing linear-response theory makes the applicability of the approximation reasonable for a certain range of parameters. We obtain the reorganization correction within the linear response theory. Firstorder perturbation provides a reasonable approximation when the Ehrenfest mean potential is small. This is true for small reorganization energy or for excitons delocalized over a large number of chromophores, N. The first condition is readily observed; if the reorganization energy is small, the response of the site energy is small enough to be treated by the first-order perturbation. The second condition suggests that the response becomes small for large systems irrespective of the magnitude of the reorganization energy. Because the site population is of the order 1/N owing to the normalization condition of the delocalized exciton populations, the response from secondorder perturbation goes as 1/N 2 . This analysis indicates that the higher-order perturbations diminish with increasing N.
Here, we mention possible generalizations of this approach. The present method can be extended to treat correlated bath fluctuations by keeping the off-diagonal terms in the memory matrix. The correlated site-energy fluctuations can be obtained by solving coupled Langevin equations. Another possible extension to include non-Markovian effects of the memory function is continued fractal expansion. [72][73][74] The exponential correlation functions given by the KA or RECKA models correspond to the Lorentz-Drude spectral densities. By introducing non-delta-function memory kernels, 75 we can treat a structured spectral density and vibronic coupling. 76,77 The effect of the structured spectral density in optical and transport properties of a biological-scale light-harvesting system is being studied in our group. 78

A. Population dynamics
To demonstrate the reorganization scheme, we begin by considering an homogeneous dimer, i.e., 1 = 2 . We fix the excitonic couplings V 12 = V 21 ≡ V to 100 cm −1 and the temperature to 300 K. The initial conditions of site-energy fluctuations were sampled from a Gaussian distribution of variance 2 ≡ 2 = 2k B Tλ. Numerical results were obtained by averaging 50 000 trajectories with a timestep of 1.0 fs. To examine the accuracy of our approach, we essentially compare the KA and RECKA methods with the numerically exact HEOM approach pioneered by Tanimura and coworkers. [22][23][24][25][26][27][28][29][30] with the Lorentz-Drude spectral density of the same parameters. We chose the inverse relaxation times /V to be 0.1, 1, and 10, which approximately correspond to 531, 53.1, and 5.3 fs, respectively. The range of suggested time depends on the system studied and the assumptions made by the authors. For example, Ishizaki and Fleming 24 proposed a slow bath correlation time greater than 100 fs for bacteriochlorophyll a in the photosynthetic Fenna-Matthews-Olson complex. However, atomistic simulations have predicted much faster relaxation times 55, 57 of 5-10 fs. In Fig. 3, we compare the Ehrenfest-based methods with the HEOM for λ = 100 cm −1 and different relaxation times. The population dynamics from the KA model agrees well with HEOM, which is consistent with earlier works. 29,59,61 They show a slight deviation for the fast relaxation time of 5.3 fs. Comparing the KA and RECKA models, we have found that the effects of reorganization corrections are negligible for slow and intermediate bath relaxation times. However, the relaxation included in the RECKA approach for fast relaxation times yields results that are nearly identical to those of the benchmark HEOM approach.
Next, we turn to the case of a stronger exciton-phonon coupling of λ = 500 cm −1 . The population dynamics are presented in Fig. 4 for the same series of relaxation times as those in Fig. 3. KA and RECKA still agree well with HEOM for slow and intermediate bath relaxation times. However, they deteriorate for the faster relaxation time: KA predicts a much faster time scale for inter-site energy transfer. Owing to the lack of reorganization effects, the KA approach provides the same energy-level distributions for absorption and fluorescence, which results in larger estimated hopping rates. This deficiency can be improved by the RECKA approach that exhibits a population decay similar to that of HEOM. In terms of population dynamics, we have obtained encouraging results for the set of parameters considered here. Previous studies using the atomistic Ehrenfest-based methods neglect the Ehrenfest mean-field potential in their description of the exciton dynamics. [53][54][55][56][57] This approximation can be justified when the reorganization energies are small compared to the excitonic couplings. We have found that, in addition to the small reorganization energies, slow bath fluctuation is an important condition for justifying the neglect of the reorganization effect.
Similar to the case for the HSR, KA, and Ehrenfest approaches, our correction scheme does not provide a Boltzmann distribution in the long-time. As discussed earlier, this deficiency originates from the classical treatment of the phonon bath and cannot be improved by reorganization correction. Based on previous work by Bastida et al., 79 Aghtar et al. 29 have introduced a correction factor that leads to the fulfillment of detailed balance. Although their scheme gives the correct thermal equilibrium distribution, it overestimates dephasing rates. An explicit quantum correction for the thermal limit may be derived by considering multiconfigurational Ehrenfest dynamics, 80,81 which lies outside the scope of this paper.

B. Stokes shifts in ideal J-aggregates
J-aggregates, 82-84 named after Jelley, are relatively ordered molecular arrays whose transition dipoles are arranged in such a way that large exciton delocalization is present in the low-energy excitation bands. They have technologically relevant optical properties such as a narrow absorption band that exhibits a large absorption intensity at the expense of a large number of quasi-dark excited states. An additional optical property of J-aggregates is a small Stokes shift. The Stokes shifts of aggregates include internal energy relaxation among exciton states and the reorganization shifts of the exciton level. Because these two contributions cannot be separated, information on exciton-phonon coupling in aggregates cannot be obtained directly from experimental Stokes shifts.
Here, we discuss exciton-phonon coupling in ideal Jaggregates in terms of the Stokes shift magnitude. The Hamiltonian of ideal J-aggregates is H S = N m { |m m| + V (|m m + 1| + |m + 1 m|)}, where the cyclic boundary condition is applied, and the nearest-neighbor coupling, V , is negative. In this ideal J-aggregate, only the optical transition to the lowest exciton state is allowed, and hence the Stokes shift is identical to the reorganization shift.
We calculate the absorption and emission spectra based on the RECKA models. Essentially, the absorption spectra can be calculated from the Fourier transform of the dipole autocorrelation function, 20 In the Ehrenfest-based method, the dipole autocorrelation function can be obtained from the ensemble average of the time-evolution operator, 53 where U mn is the (m,n)th element of the time-evolution operator for H S + H SB (t), and μ m denotes a transition dipole moment of mth site. Relaxed fluorescence spectra can be obtained from the dipole autocorrelation function averaged with initial conditions in the electronic excited-state. 85 When we propagated wavefunctions and calculated Eq. (34) for the fluorescence spectra, we sampled (t = 0) from a Gaussian distribution of same variance and mean value of −2λ/N. The lowest exciton state was selected as the initial condition for the exciton propagation.
To compare our results with those previously published, we focus on a monomer spectra where analytical solutions are available for the Lorentz-Drude spectral density. 86,87 For slow relaxation time, the Stokes shift is known to become 2λ, and the lineshapes become Gaussian. In the opposite limit (i.e., fast relaxation), the absorption and fluorescence spectra become the same without the Stokes shift. In Fig. 5, we present absorption and fluorescence spectra obtained from the RECKA approach. The lineshape is Gaussian for the fast relaxation time of = 10V , whereas it becomes Lorentzian for slow relaxation times. To describe the Stokes shifts as a simple function, we introduce the dimensionless parameter κ = / . 87 Fig. 6(a) illustrates the Stokes shifts of the monomer as a function of κ. For the KA model where the twotime correlation is given by an exponentially decaying function, we have found that the Stokes shifts are exponentially Figs. 6(b)-6(d) illustrates the fits for α for the dimer with V = −10 , −100, and −1000 cm −1 , respectively. We have ob-tained the same exponential dependence of κ, and the Stokes shifts of the dimer become half of those of the monomer. This result suggests that the delocalization of the exciton over the aggregate can lead to smaller reorganization shifts. Here, α differs slightly for the dimer with different excitonic couplings.
Finally, in Fig. 7, we have obtained the Stokes shifts as a function of the number of sites. As expected from the correction term, K m ∝1/N; Stokes shifts are inversely proportional to N. From the above results for the Stokes shifts, we have obtained a simple expression for the Stokes shifts for the ideal J-aggregates Here, α = 0.7-0.8. Stokes shifts attributed to the reorganization vanish when the bath correlation time is fast or when N is large. In macroscopic systems where N is sufficiently large, the reorganization shifts essentially becomes zero. This result strongly suggests that Stokes shifts observed in J-aggregates result from the internal energy relaxation within the exciton manifold. Three-pulse photon echo peak shift experiments 88 have indicated that in cyclic dye aggregates the exciton-phonon coupling decreases as the ring size increases, which is consistent with the present result. We stress that Eq. (35) has been obtained by using the ideal Jaggregates, where the exciton can delocalize over the entire aggregate. In a inhomogeneous system, exciton delocalization length would be limited by static disorder. In this case, N in Eq. (35) would be replaced by the inverse participation ratio. 89

IV. CONCLUSION
In this paper, we have discussed the connection between the atomistic Ehrenfest equation and stochastic bath models. By exploiting Mori's projection operator technique, we have discussed that stochastic models such as HSR can be derived as approximations of the atomistic Ehrenfest equations for exciton-phonon systems. Next, we have proposed a reorganization correction scheme to the KA model. The proposed RECKA scheme produces more accurate population dynamics than the KA model, especially in the strong excitonphonon coupling regime, and provides similar time scales for inter-site energy transfer as compared to the HEOM. The short-time dynamics in the intermediate regime are also more accurate when compared to the KA model. Finally, we have obtained a simple expression for predicting the reorganization shift in ideal J-aggregates. For the stochastic model with exponential correlation functions, the reorganization shifts depend exponentially on the inverse relaxation times. From our findings, the reorganization shifts can be described by three parameters: the monomer reorganization energy, bath relaxation time, and the exciton delocalization length. This simple relationship allows the physical origin of the Stokes shift to be understood and will be useful to separate the reorganization shift from internal energy relaxation.