An efficient method to calculate excitation energy transfer in light harvesting systems. Application to the FMO complex

A master equation, derived from the non-Markovian quantum state diffusion (NMQSD), is used to calculate excitation energy transfer in the photosynthetic Fenna-Matthews-Olson (FMO) pigment-protein complex at various temperatures. This approach allows us to treat spectral densities that contain explicitly the coupling to internal vibrational modes of the chromophores. Moreover, the method is very efficient, with the result that the transfer dynamics can be calculated within about one minute on a standard PC, making systematic investigations w.r.t. parameter variations tractable. After demonstrating that our approach is able to reproduce the results of the numerically exact hierarchical equations of motion (HEOM) approach, we show how the inclusion of vibrational modes influences the transfer.

The FMO complex consists of three identical subunits, called "monomers", each containing eight [23,24] bacteriochlorophyll (BChl) molecules (see Fig. 1). It is assumed that BChls 1, 6, and 8 are located at the baseplate and receive electronic excitation captured by the chlorosomes [27]. BChls 3 and 4 are located in the vicinity of the reaction center. To understand the excitation transfer mechanisms, many experimental and theoretical investigations, that took advantage of the availability of high resolution crystal structure, have been performed. In the present work BChl 8 is not taken into account. The figure is created using PyMOL [26]. * Electronic address: eisfeld@mpipks-dresden.mpg.de Supported by theoretical modeling a good overall understanding of the local energies of the chlorophyll molecules and their dipole-dipole couplings has been obtained, although there is still a large variance between values found in the literature [28].
The present interest in the FMO complex has been greatly stimulated by the evidence for wavelike energy transfer through quantum coherence in photosynthetic systems [10]. Since the electronic coherences are strongly influenced by the coupling of electronic excitation to vibrational modes of the BChl molecules and by the coupling to the protein environment, a detailed modeling of this interaction is important. Information on this coupling has been obtained e.g. from fluorescence line narrowing [8,29] or from theoretical simulations [16,21,28], which indicate that the spectral density, describing the coupling of the electronic excitation to vibrational modes, is quite structured. Such a structured spectral density will lead to non-Markovian effects. The influence of these non-Markovian effects have been investigated by various theoretical methods [30][31][32][33][34][35][36]. In general these methods are difficult to apply to large systems and/or arbitrary spectral densities.
In the present paper we present an efficient approach, which allows us to treat structured spectral densities with minimal requirements on the computer capabilities. All the calculations presented below can be performed on a standard PC within a few minutes. The method is based on the non-Markovian quantum state diffusion (NMQSD) approach [37][38][39][40][41][42] and allows us to treat the whole range from coherent wavelike motion to incoherent hopping. We have applied this method, within the zeroth order functional expansion (ZOFE) approximation, to describe the optical and transfer properties of molecular aggregates [43][44][45]. Within the ZOFE approximation it is possible to derive a non-Markovian master equation [40], which we use in this paper (in the following we will call it ZOFE master equation). This approach enables us to calculate the excitation dynamics of fairly large systems, e.g. aggregates consisting of more than ∼ 30 chromophores, taking structured spectral densities for the exciton-phonon coupling into account.
Although our method contains an approximation, from a previous study on absorption spectra of smaller complexes we expect that for parameters relevant for the FMO complex the approximation works well [45]. This will be confirmed in the following, where we will compare excitation energy transfer, calculated using the ZOFE master equation, with results of the so-called hierarchical equations of motion (HEOM) approach. The HEOM approach has become popular [18,31,[46][47][48][49], since it allows in principle converged results to be obtained. However, the numerical effort grows rapidly with increasing system size and for spectral densities that deviate from the Drude-Lorentz (D-L) form. To overcome these difficulties extensive parallelization (e.g. by using graphics processing units [18]) is necessary.
Since in our method we are not restricted to a simple form of the spectral density, we can investigate how energy transfer depends on its details. In particular we show that the transfer dynamics is not sensitive to coupling to modes with energies above the largest energy difference between the purely electronic eigenstates of the FMO complex (which is roughly 500 cm −1 ). Furthermore, we demonstrate that the transfer dynamics depends on the details of the local structure of the spectral density in the region, where the FMO complex has electronic transition energies.
The paper is organized as follows: In Section II we introduce the model used to describe the FMO complex and we present the ZOFE master equation that we use for the numerical calculations. In Section III we first demonstrate that with the ZOFE master equation we are able to reproduce calculations obtained with the HEOM method. Then we study the influence of different parts of the spectral density. Finally, in Section IV we conclude with a summary and an outlook.

II. MODEL
In the following we will briefly introduce the model Hamiltonian that we use to describe the FMO complex. Further details can be found e.g. in Ref. [45,50].
In the equations below, we set = 1. For each chromophore we take two electronic states into account. Since we are interested in excitation transfer of a single electronic excitation, we restrict ourselves to the one-exciton subspace which is spanned by the states | π n = | g · · · | e · · · | g , where chromophore n is excited electronically and the other BChls are in their electronic ground state. In this basis the electronic "system" part of the FMO complex is ε n | π n π n | + N n,m=1 V nm | π n π m |, (2) where ε n are the electronic energies of the chromophores and V nm are the couplings between them. The "environment" of vibrational modes is taken in the form and the coupling of electronic excitation to these vibrations is expressed through Here the κ nλ , which describe the coupling of the electronic excitation on BChl n to the local vibrational mode λ at this chromophore, are related to the spectral density which will be taken as a continuous function in the following. For a thermal initial state ρ th of the environment the environment correlation function α n (t − t ′ ) is related to the spectral density by the standard expression Note that we assume that the local environments of different chromophores are uncorrelated, which is in accordance with recent calculations [17]. Finally the total Hamiltonian in the one-exciton space is given by In the following we consider electronic excitation transfer. Thus we will focus on the reduced density operator ρ(t) in the electronic subspace, which is the trace of the total density operator ρ tot (t) of system and environment over the vibrations Its diagonal elements in the | π n basis describe the timedependent populations of the chromophores.

A. Numerical method
We calculate the reduced density operator by applying a recently developed method, which is based on the non-Markovian quantum state diffusion approach [38,41,43,45]. Within the ZOFE approximation (described in detail in Ref. [39,45]) one can derive a convolutionless master equation for the reduced density operator [40] Here the operator is determined by the auxiliary equation 0 (s, s) = −| π n π n | (we use the same notation as in Ref. [45] where the subscript 0 indicates that the operators are obtained from the ZOFE approximation). For the numerical implementation, we approximate the environment correlation function as a sum of exponentials, i.e. α n (τ ) = M j p nj e iΩnj τ with in general complex prefactors p nj and complex frequencies Ω nj . Such a form allows us to derive a closed set of coupled differential equations for the operatorsŌ In the numerical implementation we expand the equations above in the | π n basis. To determine the Ω nj we proceed similarly as in Ref. [51,52]. However, in the present work we do not use the Matsubara decomposition for the coth appearing in the correlation function α n (τ ) (see Eq. (6)) but use the partial fraction decomposition proposed by Croy and Saalmann [53]. The latter has superior convergence properties and uses poles in the complex plane (not only on the imaginary axis, as is the case for the Matsubara decomposition). Since in the present approach one can use hundreds of exponentials, in principle it is no problem to approximate the environment correlation function α(τ ) for arbitrary spectral densities J(ω) and temperatures with high accuracy.

A. Comparison with HEOM calculations
First we will compare the energy transfer calculated with the ZOFE master equation approach with the HEOM results of Ref. [32]. For the comparison we restrict ourselves to the same model used in Ref. [32] where only a single monomer unit of the FMO trimer was considered and the recently discovered eighth BChl was not taken into account (see Hamiltonian in Table I).   For the comparison we have taken the transition energies and couplings from [16]. In Ref. [32] a Drude-Lorentz (D-L) spectral density was used to describe the coupling to the environment. For the decomposition of the environment correlation function as a sum of exponentials, described above, we represent the spectral density by a sum of anti-symmetrized Lorentzians. With 10 antisymmetrized Lorentzians we could fit the D-L spectral density of Ref. [32] perfectly in the relevant energy range. The resulting spectral density is shown as the solid curve in Fig. 2(a).
The calculated transfer resulting from this spectral density is shown in Fig. 3. The curves show the time-dependent excitation probabilities of the individual BChls -the dashed lines are the HEOM results taken from Ref. [32] and the solid lines show the results calcu- lated with the ZOFE master equation (see Eq. (9)). In the upper panels the transfer is calculated for a temperature of 77K and in the lower panels for 300K. Initially all excitation is either on BChl 1 (left panels) or on BChl 6 (right panels). In all these four cases there is quite good agreement between the HEOM and the ZOFE results, showing that the ZOFE approximation preserves all the detailed features of the transfer, e.g. the strong oscillations at 77K.

B. Influence of high-energy modes
Since the ZOFE method enables us to consider an almost arbitrary form of the spectral density, we can now investigate the influence of different contributions to the spectral density.
Often, the so-called reorganization energy is used as a global measure for the coupling strength to the (environmental) vibrational modes depending on the spectral density J(ω) [4]. Then strong coupling (compared to the electronic interactions V nm ) is identified as λ ≫ V nm and weak coupling as λ ≪ V nm . In the following we will investigate the influence of contributions to the spectral density in different energy regions and we will show that in many cases the reorganization energy is not a reasonable measure for the coupling strength.
First let us consider the influence of the high energy part of the spectral density. As an example we consider a spectral density (dotted curve in Fig. 2(a)), which coincides with the D-L spectral density of Ref. [32] for low energies, but has a large additional peak centered at about 1600 cm −1 , resulting in strong additional coupling to high energy modes. We have calculated the transfer through the FMO monomer for this modified spectral density using the ZOFE master equation. The resulting transfer is exactly the same as that calculated by the ZOFE method for the original spectral density of Ref. [32] shown in Fig. 3. This shows that despite contributing significantly to the reorganization energy, the high energy part above 500 cm −1 (up to 500 cm −1 the spectral densities are the same) is of no relevance for the system dynamics. The reason for that lies in the fact that the system energies, i.e. all the transition energies between the electronic eigenenergies of the FMO monomer (indicated by the vertical lines in Fig. 2), are in the energy region below 500 cm −1 . That means that the detuning between these system energies and the high energy part of the spectral density is so large (compared to the coupling strength) that the high energy modes do not couple to the system. However, the reorganization energy does not take this fact into account and thus gives in this case a wrong indication for the strength of the coupling to the modes.
Similarly, a faster drop-off, as given by the spectral density shown by the curve (triangles) in Fig. 2(a), does not affect the transfer dynamics, since in the region below 500 cm −1 it also agrees with the original spectral density.

C. Influence of low-energy modes
The spectral densities we considered so far have a smooth form without any prominent structure. However, a more realistic spectral density consists of several distinct peaks, that embody the strong coupling to the intramolecular modes of the BChls, and is taken e.g. from measured fluorescence line narrowing spectra [16,29] or from atomistic simulations (see e.g. Ref. [21,28]).
In the following we will investigate the influence of such a structured spectral density and the influence of how the coupling strength is shared among the peaks at different energies. To this end we consider a spectral density consisting of four distinct peaks. To obtain information about the importance of the individual peaks we vary their height in a range of parameters where from previous investigations [45] we are confident that the ZOFE approximation gives reliable results.
To distinguish the influence of the structure from that of the overall coupling strength, we want to roughly keep the total coupling strength in the relevant energy region (below ∼ 550 cm −1 ) constant. To this end we first intro- Transfer calculated with ZOFE master equation, as in Fig. 3, but for the three different spectral densities in Fig. 2(b)-(d) and the D-L spectral density in Fig. 2(a).
duce an effective reorganization energy where we take E min = 0 cm −1 and E max = 550 cm −1 as the relevant energy range, that contains all electronic transition energies of the FMO monomer (see the vertical lines in Fig. 2) [58]. For the D-L spectral density in Fig. 2(a) we get a value of λ eff = 31 cm −1 for the effective reorganization energy.
To investigate the influence of the structure of the spectral density, we consider three different spectral densities, that all have the same effective reorganization energy as the D-L spectral density. But in contrast to the D-L spectral density, these spectral densities consist of four prominent peaks, whose intensities are varied. These spectral densities are shown as the solid curves in Fig. 2(b)-(d). From Fig. 2(b) to (d), one of the peaks at a time is enhanced. For comparison again the D-L spectral density is shown in each figure as the dotted curve. Note, that due to the division by ω in the definition of the reorganization energy (see Eq. (14)), to keep the value of λ eff constant, the peak intensity of the high-energy peaks becomes much larger than that of the low-energy peaks, when they are enhanced.
To investigate the influence of these changes in the spectral density, we calculated the excitation transfer through the FMO monomer using the ZOFE master equation. The resulting time-dependent excitation probabilities are shown in Fig. 4 in the same way as in Fig. 3, but now for the three different spectral densities of Fig. 2(b)-(d) and the D-L spectral density.
One sees that even though the overall coupling strength in the relevant energy region is kept constant, the transfer changes clearly when the structure of the spectral density is varied. In Fig. 4(a) for instance, the population on BChl 3 after 1 ps changes almost by up to 50% depending on the structure of the spectral density.
Obviously, not only the overall coupling strength in the relevant energy region is important, but also the local distribution of the coupling strength.
When the temperature is increased (lower panels), the spreading of the curves decreases, showing that the structure of the spectral density becomes less important at higher temperatures.

IV. SUMMARY AND OUTLOOK
In the present paper we have discussed excitation energy transfer in the FMO complex, emphasizing the role of the coupling of the excitation to vibrational modes of the BChl molecules.
a. ZOFE master equation: To handle the coupling of the electronic excitation to vibrational modes the numerical calculations have been performed using the ZOFE master equation, utilizing the ZOFE approximation, which allows to calculate the transport for the FMO complex within a few minutes on a standard PC. We validated the applicability of the ZOFE approximation by comparing with the exact HEOM results of Ref. [32], where a Drude-Lorentz spectral density has been used.
b. Reorganization energy: Since we are not restricted to a certain spectral density (SD), we investigated the influence of the parts of the SD in different energy regions. As expected, the high energy part of the SD (above ∼ 500 cm −1 ) does not influence the transport properties. In this context it is worth noting that the so-called reorganization energy, which is often used as a measure for the strength of the coupling to the environment, is often not a reasonable measure, since it takes into account the SD at all energies. In particular SDs, that have the same reorganization energy but, in the relevant energy region, have different shapes, affect the system dynamics differently.
c. Structured spectral density: Being aware of this we used an "effective reorganization energy" (ERE), which is a local measure in the energy region of system transitions [59], to investigate the dependence of the transfer on the local structure of the SD. We find that even when the structure is changed only slightly and the ERE is kept constant, we observe a change in the transfer to BChl 3 up to 50%. This influence of the structure of the SD decreases when temperature is increased.
d. Electronic energies and couplings: Note that the transfer of course depends also strongly on the precise values of the electronic couplings V nm between the BChls and their transition energies. For example we found that the final population of BChl 3 increases by roughly 50% when the electronic couplings are increased by 25%. Thus care has to be taken not to overestimate theoretical results obtained for a certain set of parameters.
e. Outlook: Since measured spectral densities show many peaks in the energy region where electronic tran-sition energies are located [29], in further work the interplay between the electronic Hamiltonian and highly structured SDs should be studied in more detail.
For such systematic studies the NMQSD (and accordingly the ZOFE master equation derived from it) seems to be well suited, due to the short calculation times. Another big advantage of the computational efficiency of the approach is that larger systems are tractable now -the full trimer with 24 chromophores and a structured SD can be calculated on a standard PC within a few hours. One goal of the present study was to further investigate the rather abstract ZOFE approximation [39,45]. The quite good agreement with the exact HEOM results of Ref. [32] makes us confident that it is worth to further explore the NMQSD method.
The NMQSD approach is also suitable to describe other light harvesting complexes and molecular aggregates. While in the case of the FMO complex the coupling to high energy vibrational modes with energies larger than 1000 cm −1 does not influence the energy transfer markedly, in molecular aggregates of organic dyes these modes have a strong impact on the electronic excitation dynamics [54]. Furthermore they have energies comparable to the dipole-dipole interaction which influences strongly e.g. the aggregate absorption spectra [55]. In previous investigations we have found that typ-ical features of molecular aggregate spectra, such as the exchange narrowing of the J-band and the broad H-band, are reproduced by the NMQSD calculations [44].
The NMQSD method is also readily applicable to treat optical spectroscopy. In particular it should be well suited to investigate the influence of the coupling to vibrational modes on the signals observed in 2D experiments. Such a study might be revealing, since the "power spectrum" obtained from 2D measurements in Ref. [10] shows strong resemblance to the fluorescence line narrowing spectrum of Wendling et. al [29] .
To obtain further information how individual vibrational modes influence the transport it might be useful to simulate the transport using networks consisting of quantum [56] or classical oscillators [57].

V. ACKNOWLEDGMENT
Financial support from the DFG under Contract No. Ei 872/1-1 is acknowledged. AE thanks Alán Aspuru-Guzik for the hospitality. We thank John Briggs for helpful comments and Alexander Croy for many useful discussions and for providing the numerical code for the coth expansion of the environment correlation function.