## Abstract

Recent experiments on light-harvesting complexes have shown clear indication of coherent transport of excitations in these aggregates. We discuss the theoretical models that have been used to study energy transfer in molecular aggregates, beginning with the early models of Förster and Davydov and ending with the theoretical models of the present day.

## 1. Introduction

Just after the second world war, Förster [1,2] introduced a model for excitation transfer between donor and acceptor molecules and showed that the rate of energy transfer, *k*_{ET}, (which we would now classify as incoherent energy transfer) depends on the distance between donor and acceptor, *R*, as *k*_{ET}∼(*R*_{0}/*R*)^{6}, where *R*_{0} can be calculated from the overlap of the donor emission spectrum and the acceptor absorption spectrum. This theory has been immensely useful in the intervening years and is now used routinely in single-molecule biophysical spectroscopy experiments. Recently, a number of theoretical generalizations of the theory make it applicable in situations that are more complex than the original one that Förster discussed [3].

At the same time, Davydov, in the USSR, following earlier work by Frenkel [4,5], introduced the idea of molecular excitons, in order to discuss the low-temperature spectroscopy of molecular crystals and aggregates. His work was translated and popularized by Kasha, appearing in 1962 [6,7]. Davydov's description of the electronic (and vibronic) states of molecular crystals emphasized the effects of coherence: the states were linear combinations of excitations localized on different sites. The theory predicted that in molecular crystals in which there were two translationally inequivalent molecular sites (but connected by a symmetry element), each spectral line would be split into two, giving rise to Davydov splitting. If the electronic states were not coherent superpositions, but were localized on a site, there would be no splitting. Experiments found the Davydov splittings in a number of crystals at low temperature [8]. The theoretical model of Davydov has been generalized (e.g. to include charge transfer and other two site states) [9], but the basic ideas have continued to be fertile for experimenters and theorists alike.

Förster's model only considered incoherent energy transfer between localized excitations; Davydov's model only considered excitations that were coherently coupled among many sites, although he discussed the effects of phonon–exciton interaction on line shapes. A more general model allows for both coherence and interactions with phonons and other excitations in the environment of the excitations that lead to decoherence, so that the Davydov and Förster pictures emerge as limiting cases.

As experimental methods improved, it was possible to experimentally determine the degree of coherence (and decoherence) in certain systems. For example, Schwoerer & Wolf [10,11] examined isotopically mixed naphthalene crystals (small amount of naphthalene-h8 in a crystal of naphthalene-d8). The naphthalene-h8 excited triplet state is approximately 100 cm^{−1} below the same state of the deuterated molecule. At sufficiently high concentrations of the h8, both monomers and dimers of h8 (made up of the nearest-neighbour translationally inequivalent molecules) can be seen in optical spectroscopy, electron paramagnetic resonance and electron nuclear double resonance experiments at very low temperature. High-resolution optical spectra showed two lines, one 1.25 cm^{−1} above and one 1.25 cm^{−1} below the monomer transition. Thus, the eigenstates of the dimer at *T*=0 *K* were taken to be
1.1
The matrix element of the Hamiltonian between the two terms is *J*=〈*A*_{ex}*B*_{g}|*H*|*A*_{g}*B*_{ex}〉, and the splitting between the two eigenvalues is 2*J*, so that *J*=1.25 cm^{−1} [*Φ*_{−} was assumed to be the lower state (*J*>0)]. This is the first example of the ‘mini-exciton.’

The question arose as to how to describe the experimental observations, in particular, whether the energy transfer between the two localized states should be considered to be incoherent hopping (the Förster picture) or coherent (the Davydov picture). The interaction with phonons in the crystal lattice will cause scattering from *Φ*_{+} to *Φ*_{−} and vice versa. If the scattering rate from the upper state to the lower is *τ*^{−1}, then by detailed balance, the scattering from the lower to the upper state is given by . If *W*>2*J*, then the better description is hopping, and if *W*<2*J*, the better description is coherent transfer. Botter *et al.* [12] measured the phase memory time of the triplet spins and the zero magnetic field resonance frequencies (using optical detection) between 1.2 *K* and 4.2 *K*, from which they were able to conclude that *W*/2*J*=10^{−3} at 1.2 *K* and 10^{−2} at 4.2 *K*, and thus the coherent description of the states and the energy transfer is preferred.

As optical experiments improved, it was possible to perform photon echo experiments on isolated molecules and isolated pairs of molecules in a host crystal [13] at liquid He temperatures. This showed that the coherent description of the pair states was correct, at least at liquid He temperatures. There were a number of other groups that contributed to the study of dimers (and one-dimensional molecular excitons [14–16]) in this period.

At very low temperatures, the coherent description of the electronic states of aggregates and crystals was basically correct, so energy transfer was dominated by coherent interactions; that is, the electronic (or vibronic) states of the system were wavelike and were scattered weakly by phonons (and defects). At higher temperatures, most experiments continued to be described using the Förster mechanism, so the energy transfer was considered to be incoherent; that is, the scattering of the electronic states was large enough to localize the excitation, which then hopped from site to site owing to the fluctuation of the interactions. A number of theories were proposed in order to understand how the description changes as one goes from low temperature to high temperature.

Experimental measurements continue to improve. Recently, the Fleming [17–19], Engel [20] and Scholes [21] groups have been able to see coherent oscillations over a few hundred femtoseconds in multi-chromophoric aggregates, after which time, decoherence sets in. These light-harvesting aggregates are part of biological systems that use the energy of the light to drive chemical reactions. Some of the questions that many groups are trying to answer include the following. Does this long-lived coherence play a role in the efficiency of this process? How does the coherence and decoherence depend on the structure of the system? Does the protein in which these chromophoric molecules are arrayed help or hinder the coherence? Does decoherence play a role in optimizing the efficiency of the transport of the excitation to the reaction centre?

Theory also continues to improve, although the basic ideas that are used were set down by Förster and Davydov. In the rest of the paper, we will discuss the theoretical descriptions briefly.

## 2. Theoretical descriptions

In order to theoretically describe the dynamics of excitations in aggregates and crystals, one must use a density matrix formalism because it is impossible to control the other degrees of freedom (i.e phonons, librations, two-level defects, etc.) in the crystal or protein environment of the chromophores; all one can do is set the temperature. The equation of motion for the density matrix for the entire system, *ρ*, is
2.1
where the bracket represents a commutator. The Hamiltonian is given in general as
2.2
where *H*_{s} is the Hamiltonian of the system under study (e.g the dimer of h8 molecules mentioned above), *H*_{b} is the Hamiltonian of the bath, i.e the other degrees of freedom (phonons, etc.) and *V* is the interaction between the system and the bath. Normally the interaction term *V* is considered to be a sum of system operators multiplied by phonon (or vibrational) operators. In the site basis for the system (i.e. *not* the eigenbasis), one can characterize the diagonal terms as site-energy fluctuations owing to phonons and the off-diagonal terms as transfer matrix element (i.e the *J* above) fluctuations owing to phonons. In the eigenstate basis for the system, the fluctuations in eigenenergies are due to both site energy fluctuations and transfer matrix element fluctuations.

Since we only measure properties of the system, we take the trace of *ρ* over the bath degrees of freedom to find the *reduced* density matrix of the system, *σ*(*t*)=Tr_{b}[*ρ*(*t*)]. Using a variety of methods, one can derive an exact equation for the dynamics of *σ*(*t*), called the stochastic Liouville equation (SLE) [22–24],
2.3

There are very few (almost none!) models for which the exact form of *K*(*t*) can be found; therefore, one must use sensible approximations in order to proceed. The most common (and perhaps the most important) of these is due to Redfield [25], who studied the relaxation of nuclear spin states in condensed phases.

Redfield made the approximations (quite appropriate for nuclear spins) that (i) the coupling of the system to the bath was weak and (ii) the dynamics of *σ*(*t*) was slow compared with the time scale of relaxation of the bath, which governs the time scale of relaxation of *K*(*t*). Therefore, one can replace *σ*(*t*−*τ*) by *σ*(*t*) in equation (2.3). In addition, if one restricts the time of examination of *σ*(*t*) to times long compared with the relaxation of the bath, one can take the upper limit of the integral to infinity. The Redfield equation for the reduced density matrix becomes
2.4

Here, *R* is the relaxation matrix (it is complex in general). The matrix elements of *R* in the eigenstate basis are Fourier transforms of phonon operator time correlation functions. In the eigenbasis of *H*_{s}, {|*n*〉}, we can write this as
2.5

The diagonal elements, *σ*_{nn}(*t*), are the populations of the eigenstates, and the off-diagonal elements, *σ*_{nm}(*t*), are called the coherences. Note that the coherences oscillate with frequencies given by the difference between the eigenenergies. The populations do not oscillate unless they are coupled to the coherences and the coherences are non-zero. Of course, the *site* populations can oscillate, even if the eigenstate populations do not.

Owing to the assumption of weak coupling, the form of *R*_{mn,pq} is equivalent to second-order perturbation theory. The equation for the populations exhibits this clearly,
2.6
If one can neglect the coupling of the populations to the coherences, the equation for the populations becomes a simple kinetic (i.e master) equation.

There are some immediate consequences of the Redfield model. (i) The scattering rates between two states satisfy detailed balance. (ii) The decay of the populations is governed by the scattering among the states, but the decay of the coherences has an added term, the pure dephasing rate, which is related to the time-dependent correlation function of the fluctuations of the system eigenenergy levels owing to the interaction with the bath. (iii) Since the pure dephasing rate is a zero-frequency Fourier transform, it normally goes to zero as the temperature goes to zero because the spectral density of phonon (and other) modes goes to zero at zero frequency. (iv) In general, there will be coupling between populations and coherences [26–28].

The Redfield model has been applied to many spin resonance and optical experiments with great success. One must be careful, however, because it is inherently a weak coupling model. One consequence of this is that for some values of the parameters, the predicted populations can become negative or larger than 1. This indicates a breakdown of the model. If this happens, we have to go beyond second order or define a different system Hamiltonian, by, for example, making a unitary transformation that puts some of the coupling to bath modes into the system Hamiltonian [27,28]. Another approach is to use the Lindblad equation, which guarantees positivity of the reduced density matrix [29].

Haken and co-workers [30–33] introduced a model, following earlier work by Sewell [34], that has many attractive features. The model replaces the time dependence of the phonon operators in the Redfield model by classical stochastic Gaussian white noise, i.e. it assumes that all time correlation functions have a *δ*(*t*) time dependence and all fluctuations are Gaussian. The ensuing equations are the same as the Redfield equations given above, but with the matrix elements of *R* calculated with the Haken–Strobl–Reinecker (HSR) assumptions. And these equations are *exact* within these assumptions because these assumptions cause the higher than second-order cumulants to be zero, and so the weak coupling equations (i.e. Redfield) that only consider the second cumulant become exact.

There are a number of concerns about the HSR model. (i) Assuming a *δ*(*t*) time dependence means that the relaxation of the bath modes is instantaneous and that the bath has energies on all scales, or equivalently, the temperature is very large compared with any energy differences in the system. Thus, detailed balance is lost and all equilibrium populations are equal. (ii) the model predicts that the pure dephasing rate is independent of temperature (although one can always put in a *T* dependence by hand). It is known that the pure dephasing rate should vanish at low temperatures. (iii) The model predicts that all the Fourier transforms are independent of frequency, which neglects the differences in the phonon spectral density of real systems. (iv) the original HSR model neglected the correlations between fluctuations on different sites and on the fluctuations between *J* and site energies. Neighbouring molecules in a crystal or protein environment should exhibit correlations among these quantities and putting them into the model can qualitatively change the results [35].

In spite of these concerns, the HSR model, including all correlations, is useful and important, especially at higher temperatures, where (i)–(iii) are less troubling.

Let us now turn to explicit phonon models in which the phonons are assumed to be harmonic, , where *b*_{k} and are the usual Boson operators. The most common form of interaction is a site diagonal excitation–phonon interaction, i.e. the energy at site *n* in the system varies with the *k* phonon coordinate *Q*_{k} as ) [where ]. This form is usually interpreted as the change in site energy owing to a quasi-localized oscillator made up by a linear combination of the normal harmonic modes. For the asymmetric dimer, the Hamiltonian then becomes
2.7
Note, that if *E*_{1}=*E*_{2} and *g*_{1k}=−*g*_{2k} (the symmetric dimer), this is the famous spin–boson problem [36,37].

In the limit that *J*=0, the eigenstates are easily found by displacing the coordinates an amount depending on which molecules are excited (polaron transformation). The eigenenergies become .

The second term is the reorganization energy of the lattice when the excitation is on site *n*, and is related to the Huang–Rhys factor used in spectroscopy [38]. Also note the resemblance to the theory of electron transfer.

With this form for the exciton–phonon interaction, we can calculate the population relaxation rate, *Γ*, in the Redfield model to lowest order (i.e. second-order perturbation theory),
2.8
where is the Boson occupation number and *β*=1/(*kT*). To this order in perturbation theory, only one-phonon absorption or emission processes are allowed. In order to go further, we need explicit forms for the frequencies and the coupling constants. Or we can make the following *crude* approximations and assume that (i) the sum over the harmonic modes yields a decaying exponential function, , where *ν* is proportional to the phonon bandwidth, (ii) , and (iii) (*g*_{1k}−*g*_{2k})^{2}≈*g*^{2}(1−*c*), where *c* is the correlation coefficient between the local modes on sites 1 and 2. Then,
2.9
where 〈*ωg*^{2}〉(1−*c*) is the average of *ω*_{k}(*g*_{1k}−*g*_{2k})^{2} over the phonon band(s). This is crude, but gives insight into the process of decoherence. If , then the decoherence is weak (this is the motional narrowing limit, in which the excitation transfer between sites is so rapid that the phonons cannot keep up). In the opposite limit (), the phonon bandwidth is so large that there are so few decoherence processes at 2*J* and again the decoherence is weak. Calculations going to higher order confirm this qualitative picture. We note the similarity to the results of the high-temperature Brownian oscillator model introduced by Mukamel for nonlinear spectroscopy [38].

Almost all theoretical calculations have been made with the assumption of site-energy fluctuations only. However, the Haken–Strobl model shows that the fluctuations of the transfer matrix element, *J*, are important. In particular, the long time diffusion constant at high temperatures is dominated by these fluctuations. It is important for a full understanding of the problem of energy transfer that we take these fluctuations (and the correlation of fluctuations) into account in theoretical models.

We can deal with the Hamiltonian of equation (2.7) in another way, by making the polaron transformation to remove the linear displacements of the harmonic phonons for the states |1〉 and |2〉. The transformed Hamiltonian is *H*′=*UHU*^{†} where , so
2.10
where . Next, we average the term in *J* over the thermal phonon bath and add and subtract that average from *H*′ [27,28,39–41] to find
2.11

Note that the perturbation term can be interpreted as the fluctuation of the transformed transfer matrix element from its average value. This perturbation is small in both the weak coupling limit () and the strong coupling limit (), suggesting that this is a particularly good starting point for studying excitation transfer. In the strong coupling limit, second-order perturbation theory reproduces the Förster picture exactly. In the weak coupling limit, it reproduces the Davydov picture with weak scattering of the exciton states.

Another approach is to keep the upper limit of the integral in equation (2.3) as *t*, and keep the second-order form of *K*(*t*), to get (after some mathematical manipulations),
2.12

This is called the time-convolutionless second-order equation. Jang *et al.* [42] and McCutcheon & Nazir [43] used this equation (with careful consideration of the initial conditions) to discuss excitation transfer in a two-molecule system. They first applied a polaron transformation to the Hamiltonian [27,28,39–41] that, as we have seen, leads to a small perturbation in both the weak coupling and strong coupling regimes of exciton–phonon coupling. The results show qualitative differences from those of the Redfield approach in the moderate to strong coupling regime.

Cao [44] introduced a phase-space method for going beyond the Redfield picture for the SLE, and Tanimura & Kubo [45,46] and Ishizaki & Fleming [47,48] have discussed hierarchical approaches to improving the approximations to the kernel, *K*(*t*). All of these approaches have assumed that the fluctuations are of the Gaussian–Markov type with exponential time dependence in the correlation function. Since the harmonic phonon bath is Gaussian, this is a good assumption. Ishizaki and Fleming show that for the dimer system, their approach yields results that are valid in both the weak and strong coupling regime.

Finally, the Reineker group [49,50] has treated the dimer with dichotomic noise, i.e fluctuations that take only two values with a time correlation function that is exponential. With this assumption, the hierarchical equations close [51,52]; thus, this procedure is exact in principle, but requires much computation. This is an extension of the Haken–Strobl ideas for bath correlation times that are non-zero.

In the last few years, a number of theoretical groups have examined the effect of decoherence (owing to noise or fluctuations) on the efficiency of energy transfer, particularly in the Fenna–Matthews–Olson protein [53–64] because of the exciting recent experimental results [17–21].

## 3. Variational polaron transform

### (a) Introduction

It is particularly desirable to develop a theory for energy transfer that can correctly describe the dynamics over a wide range of parameter values. Förster theory and related theories discussed earlier are perturbative in the electronic coupling and cannot correctly describe the coherent dynamics that result from strong coupling. On the other hand, many of the master equation approaches developed recently instead use the system–linear phonon coupling as a perturbation; these theories are of course unsuitable when the system–bath coupling becomes large, unless one can go to high order in the perturbation [44–48].

In order to provide a second-order perturbation method that will give good results in both the weak coupling and strong coupling limits, we will use the variational polaron method for such problems that was introduced by Yarkony & Silbey [65–67] and used by Silbey & Harris [41]. Recently, we learned that McCutcheon and Nazir [68] have also done a calculation with this method. Their results and ours coincide for the same parameter set.

### (b) Model

We use the two-state Hamiltonian introduced earlier in equation (2.7) as our starting point. This Hamiltonian assumes linear coupling between the electronic states and the bath, but is otherwise quite general. In particular, we have allowed both electronic states to interact with the same phonon modes; we can also treat the case of independent baths as a special case by requiring either *g*_{1k} or *g*_{2k} to be zero for each *k*, such that each mode only couples to one of the electronic states, or have correlation between the modes on different sites.

Before proceeding to look in more detail at this system, we re-write the Hamiltonian in the following form:
3.1
where
3.2
and
3.3
We have simply written the Hamiltonian in a form where the coupling to the bath is treated in a symmetric manner. It is always possible to write a Hamiltonian of the form of equation (2.7) in this manner by performing a polaron transform,
3.4
and by redefining the zero of energy if necessary. A nice feature of this Hamiltonian is that it does not depend on *g*_{1k} and *g*_{2k} individually, but only on their difference.

At this point, one could solve for the dynamics of this Hamiltonian by perturbing in the system–bath interaction, as was done by Redfield [25]. Alternatively, if the system–bath coupling is large, one could perform the polaron transform introduced in equation (2.10) and instead perturb in *Jθ*−*J*〈*θ*〉, as in equation (2.11). For ohmic spectral densities, *J*〈*θ*〉 is always zero, and there is no coherence in the zeroth-order Hamiltonian; as a result, the polaron-transform technique always yields completely incoherent dynamics and is qualitatively incorrect at small system–bath coupling. At large system–bath coupling, a polaron-transform is required since directly perturbing in the large system–bath coupling will not yield good results.

In the case of intermediate system–bath coupling, it is not always clear which of these two approaches will yield better results. A method that can interpolate between the two limiting cases would thus be very useful, particularly at intermediate system–bath coupling. It is such a method that is presented here.

### (c) Partial polaron transform

Instead of choosing between a complete polaron transform and no polaron transform at all, we introduce the partial polaron transform [65,66],
3.5
for some set of parameters *f*_{k}. We note that choosing *f*_{k}=*g*_{k} corresponds to performing a complete polaron transform as introduced earlier, while choosing *f*_{k}=0 corresponds to performing no polaron transform at all. In general, the *f*_{k} will be between these two limits and we will only partially transform the phonon modes.

We then obtain 3.6 where 3.7 We note that there will be a polaron shift added to both states, but it will be of the same magnitude and sign for both states so we have not included it.

Our Hamiltonian is now in a suitable form for a perturbative expansion, *H*=*H*_{0}+*V* . We take as our zeroth-order Hamiltonian
3.8
where *J*_{p}=*J*〈*e*^{2G}〉 and the angle brackets indicate an average over the bath in thermal equilibrium.

Our perturbation is then 3.9 We note that the thermal average of our perturbation is zero by explicit construction.

Our choice of *H*_{0} and *V* depends on our choice of the parameters *f*_{k}; it will of course be desirable to choose the *f*_{k} in a manner that will minimize the perturbation *V* . We will follow the suggestion of Silbey & Harris [41,67] and choose the *f*_{k} such that the Gibbs–Bogoliubov free energy of *H*_{0} is a minimum,
3.10
We note that *J*_{p} itself depends on our choice of *f*_{k}, and so this equation is not an explicit formula for *f*_{k}, but rather a self-consistent equation. It may be interesting in the future to look at other criteria for minimizing the perturbation.

For our further treatment, we will express the perturbation as
3.11
where the *u*_{i} act only on the bath subspace, and the *V* _{i}=*σ*_{i} are Pauli matrices that act only in the system subspace. We will refer to the system operators as *V* _{i} rather than *σ*_{i} to avoid confusion with the reduced density matrix *σ*. Specifically, we have
3.12
3.13

and 3.14

### (d) Perturbative solution of dynamics

Now that we have a zeroth-order Hamiltonian and a perturbation, we can solve perturbatively for the dynamics. We choose to apply the second-order time-convolutionless master equation introduced in equation (2.12). This equation is well known and has been used previously to study excitonic energy transfer [42,43,69]. We thus obtain an equation of motion for the reduced density matrix of the system in the interaction representation of *H*_{0}, as [69]
3.15
We have introduced the notation *V* ^{×}≡[*V*,⋅] and *V* ^{°}≡{*V*,⋅}. The correlation functions *χ*(*t*) and *S*(*t*) are defined as
3.16

In order to derive this form of the time-convolutionless master equation, one must assume that the bath is initially in equilibrium. Jang *et al.* have looked in detail at the effect of including the initial condition in the time-convolutionless master equation for a similar system; including the initial condition effects made some small changes to their results, but did not qualitatively change the energy transfer dynamics [42,43].

Many of the correlation functions needed to solve the time-convolutionless master equation are already known. In particular, the correlation functions of *u*_{z} are known from treating the problem without a polaron transform [69], while the correlation functions involving *u*_{x} and *u*_{y} are known from the theories that perform a full polaron transform [42,43]. In this case, we also need the correlation functions, which are a cross between the polaron and non-polaron terms. They are straightforwardly calculated as
3.17
3.18
and *S*_{xz}=*S*_{zx}=*χ*_{xz}=*χ*_{zx}=0, where
3.19

The operators, *V* _{m}(*t*)^{×} and *V* _{n}(*t*)^{°} are just commutators and anticommutators of Pauli matrices transformed to the interaction picture and are likewise easily computed as needed. We thus have all of the required terms to compute the time dependence of the density matrix via equation (3.15).

### (e) Results

We will now apply the theory outlined here to a few cases of interest. Recently, Ishizaki *et al*. [69] have looked at the rate of energy transfer as a function of system–bath coupling using several of the existing theories. We will use the same parameters as they did so that our results can be compared with theirs. In particular, we will choose *Δ*=100 cm^{−1}, *J*=20 cm^{−1} and *T*=300 *K*. The two sites are coupled with independent bath modes with the same spectral density,
3.20
for *γ*=53 cm^{−1}. We will examine the dynamics for various values of *λ*, which is the reorganization energy of the bath.

For these parameters, the dynamics are mainly incoherent and thus it makes sense to define a rate of population transfer. In figure 1, we plot the rate of population transfer from site 1 to site 2 as a function of reorganization energy, calculated using the variational polaron-transform method outlined here. For comparison, we also plot the results obtained by performing a full polaron transform and by performing no polaron transform. (The case of no polaron transform was considered by Ishizaki *et al*. [69] and our results agree with theirs.) As expected, the results obtained without performing a polaron transform are qualitatively incorrect at large reorganization energies, but the polaron-transform result correctly shows a decreasing rate as the reorganization energy becomes very large.

More interestingly, the variational polaron transform interpolates between these results and obtains reasonable results for all values of the reorganization energy. Our variational polaron transfer rates are very similar to those obtained using the hierarchical equation approach of Ishizaki *et al*. [69]. We also point out that for these parameters, the full polaron transform actually predicts the transfer rates quite well over the whole range of reorganization energies considered, but the variational technique does introduce some minor corrections at small *λ*.

Finally, we note that for , the variational polaron results are identical to the full polaron results. This is because for large enough reorganization energy, the solution to the variational condition, equation (3.10), is *J*_{p}=0 and *f*_{k}=*g*_{k}; thus, a full polaron transform is performed, even in the variational calculation. As was shown analytically by Silbey and Harris, *J*_{p} does not go continuously to zero, but does so abruptly when a critical value of *λ* is reached. It is for this reason that we see a discontinuous jump around *λ*=5 cm^{−1} in the rates predicted by the variational calculation.

Next, we consider a case where coherence is important by setting *J*=100 cm^{−1} and leaving all other parameters the same (again following Ishizaki *et al*. [69]). Since the energy transfer has coherent oscillations in this case, it is not as useful to think about a rate of energy transfer and we instead look at the dynamics as a function of time when the system is started in state 1. Figure 2 shows the population of site 1 as a function of time for four different values of the reorganization energy, *λ*.

For the two smallest values of the reorganization energy, the results obtained from the variational polaron transform are essentially identical to those obtained without performing a polaron transform. These results are in reasonable agreement with those obtained by Ishizaki *et al*. [69] using the hierarchical techniques developed by Tanimura & Kubo [45,46]. We point out that for these two cases, the full polaron result is qualitatively incorrect, yielding unphysical results for *λ*=2 cm^{−1} and failing to predict any coherence for *λ*=20 cm^{−1}.

For the two largest values of reorganization energy, the variational result is the same as the full polaron-transform result, for the reasons discussed earlier. These polaron-transform dynamics correctly show a decrease in the rate of energy transfer at large reorganization energy. In contrast, the results without a polaron transform show a qualitatively incorrect time scale for large reorganization energy, as was also seen in figure 1. We also note that the results from Ishizaki and Fleming's hierarchical technique show non-exponential kinetics at short times, as well as a few coherent oscillations at *λ*=100 cm^{−1} that are not reproduced by our variational technique. Nonetheless, the overall time scale of population transfer from our method agrees quite well with that of Ishizaki *et al*. [69].

Finally, we note that for the parameters in figure 2, the variational method always yields results that are identical to either performing a full polaron transform or not performing a polaron transform. We emphasize, however, that this is not always the case; in particular, our results at low reorganization in figure 1 show that the variational method can indeed give results distinct from the two limiting cases.

In general, we expect the variational method to be distinct from the two limiting cases when *J*_{p} is intermediate between 0 and *J*. When the reorganization energy is zero, *J*_{p}=*J*; as the reorganization energy increases, *J*_{p} decreases continuously until a critical reorganization energy is reached, at which point *J*_{p} jumps discontinuously to zero. For some parameters (such as those in figure 2), *J*_{p} has barely decreased from *J* before this discontinuous jump and thus *J*_{p} is never truly intermediate between 0 and *J*; for these parameters, the variational method always essentially reproduces one of the limiting cases. For other parameters (such as those in figure 1), *J*_{p} does decrease significantly from *J* before jumping to zero; for those parameters where *J*_{p} is intermediate, the variational calculation will yield results that are distinct from both limiting cases.

### (f) Conclusion

We have presented a method for computing the dynamics of excitonic energy transfer that combines a variational polaron transform with the second-order time-convolutionless master equation. By variationally minimizing the perturbation, we were able to apply a second-order master equation to obtain energy-transfer dynamics at small, intermediate and large reorganization energies from a single theory. Our theory has the advantage of being conceptually simple and easily computable while still being applicable over a wide range of parameter values.

## Acknowledgements

R.J.S. was supported by an Army Research Office grant (W911NF-09-0480). E.N.Z. acknowledges support from the Fonds québécois de la recherche sur la nature et les technologies (FQRNT). We would like to thank Dara McCutcheon and Ahsan Nazir for many useful discussions.

## Footnotes

One contribution of 14 to a Theo Murphy Meeting Issue ‘Quantum-coherent energy transfer: implications for biology and new energy technologies’.

- This journal is © 2012 The Royal Society