An efficient method to calculate excitation energy transfer in light-harvesting systems: application to the Fenna–Matthews–Olson complex

Since the pioneering work by Franck and Teller [1], exciton theory has been widely used to describe the optical and energy transfer properties of chromophore complexes such as J-aggregates [2, 3] or biological light-harvesting systems [4]. In recent years, particular focus has been put on the so-called Fenna–Matthews–Olson (FMO) complex (see e.g. [4–26]), which promotes energy transfer from the main light-harvesting complex to the reaction center in green sulfur bacteria [4].

The FMO complex consists of three identical subunits, called 'monomers', each containing eight [27, 28] bacteriochlorophyll (BChl) molecules (see figure 1). It is assumed that BChls 1, 6 and 8 are located at the baseplate and receive electronic excitation energy captured by the chlorosomes [30]. 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. Supported by theoretical modeling, a good overall understanding of the local energies of chlorophyll molecules and their dipole–dipole couplings has been obtained, even though there is still a large variance between the values found in the literature [31].

Zoom In Zoom Out Reset image size

The present interest in the FMO complex has been greatly stimulated by evidence for wavelike energy transfer through quantum coherence in photosynthetic systems [13]. 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 [11, 32] or from theoretical simulations [20, 25, 31], which indicate that the spectral density (SD), describing the coupling of the electronic excitation to vibrational modes, is quite structured. Such a structured SD will lead to non-Markovian effects. The influence of these non-Markovian effects has been investigated by various theoretical methods [33–39]. In general, these methods are difficult to apply to large systems and/or arbitrary spectral densities.

In this paper, we present an efficient approach that allows us to treat structured spectral densities with minimal requirements in terms of 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 [40–45] 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 [46–48]. Within the ZOFE approximation it is possible to derive a non-Markovian master equation [43], which we use in this paper (in the following we will call it the 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.

The ZOFE approximation is introduced and motivated in [42] and its application to the present problem is explained in [48]. In the present context, we would like to emphasize that the ZOFE approximation is not an expansion in a small parameter that can be interpreted easily (like the system bath coupling) but is the zeroth-order term of a functional expansion that leads to a hierarchy of coupled equations [42]. Therefore its range of validity is hard to estimate. The ZOFE approximation has been confirmed to be true in many cases: it contains, e.g., the Markov limit (Lindblad) and the weak coupling (Redfield) limit. In [48], a systematic study has been conducted for small oligomers. It was found that the method works well for small coupling between the chromophores, and also in the range from intermediate to strong coupling. It turned out that the quality of the approximation may decrease when for a peaked SD the strength of the coupling to the environment is increased (i.e. the magnitude of the environment correlation function becomes larger, which could also be caused by increasing the temperature). One aim of the present study is to further investigate the applicability of this approximation. From a previous study on the absorption spectra of smaller oligomers [48], we expect that for parameters relevant to the FMO complex the approximation works well. 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 [22, 34, 49–52], 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 [22]) is necessary.

Since in our method we do not restrict ourselves to a simple form of the SD, 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⁻¹). Furthermore, we demonstrate that the transfer dynamics depends on the details of the local structure of the SD in the region where the FMO complex has electronic transition energies.

This paper is organized as follows. In section 2, 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 3, 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 SD. Finally, in section 4 we conclude with a summary and an outlook.

In the following, we will briefly introduce the model Hamiltonian that we use to describe the FMO complex. Further details can be found in e.g. [48, 53].

In the equations below, we set ℏ = 1. For each chromophore we take into account two electronic states. Since we are interested in the excitation transfer of a single electronic excitation, we restrict ourselves to the one-exciton subspace that is spanned by the states

$$\begin{equation} |{\pi_n}\rangle=|{g}\rangle\cdots |{e}\rangle \cdots |{g}\rangle, \end{equation} \tag{ 1 }$$

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

$$\begin{equation} H_{\rm sys}=\sum_{n=1}^N\varepsilon_n|{\pi_n}\rangle\langle{\pi_n} |+\sum_{n,m=1}^N V_{nm}|{\pi_n}\rangle\langle{\pi_m}|, \end{equation} \tag{ 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

$$\begin{equation} H_{\rm env}=\sum_{n=1}^N\sum_{\lambda}\omega_{n\lambda}a^{\dagger}_{n\lambda}a_{n\lambda} \end{equation} \tag{ 3 }$$

and the coupling of electronic excitation to these vibrations is expressed through

$$\begin{equation} H_{\rm int}=-\sum_{n=1}^N|\,{\pi_n}\rangle\langle{\pi_n} |\sum_{\lambda}\kappa_{n\lambda}(a^{\dagger}_{n\lambda}+a_{n\lambda}). \end{equation} \tag{ 4 }$$

Here κ_nλ, which describe the coupling of the electronic excitation on BChl n to the local vibrational mode λ at this chromophore, are related to the SD J(ω) by [53]

$$\begin{equation} J_n(\omega)=\sum_{\lambda}|\kappa_{n\lambda}|^2\ \delta(\omega-\omega_{n\lambda}), \end{equation} \tag{ 5 }$$

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 SD by the standard expression

$$\begin{equation} \alpha_n(\tau) =\int_0^{\infty}\! \!\!\! {\rm d}\omega \,J_n(\omega)\left[\coth\left(\frac{\omega}{2kT}\right)\cos(\omega\tau)-\mathrm{i}\sin(\omega\tau)\right]. \end{equation} \tag{ 6 }$$

Note that we assume that the local environments of different chromophores are uncorrelated, which is in agreement with recent calculations [21]. Finally, the total Hamiltonian in the one-exciton space is given by

$$\begin{equation} H^e=H_{\rm sys}+H_{\rm int}+H_{\rm env}. \end{equation} \tag{ 7 }$$

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 the system and the environment over the vibrations

$$\begin{equation} \rho(t)={\rm Tr}_{\rm env}\,\rho_{\rm tot}(t). \end{equation} \tag{ 8 }$$

Its diagonal elements in the |π_n〉 basis describe the time-dependent populations of the chromophores.

2.1. The numerical method

We calculate the reduced density operator by applying a recently developed method, which is based on the NMQSD approach [41, 44, 46, 48]. In this approach, the dynamics of a quantum system coupled to an environment of harmonic oscillators is obtained from a non-Markovian stochastic Schrödinger equation. The ensemble average recovers the reduced density matrix without approximation. However, this equation is hard to solve in the general case. Therefore, we have recently applied the ZOFE approximation (described in detail in [42, 48]) to treat molecular aggregates and light-harvesting systems [46]. Within the ZOFE approximation, one can derive a convolutionless master equation for the reduced density operator [43]

$$\begin{eqnarray} &&\fl\partial_t\rho(t)=- \mathrm{i}[H_{\rm sys},\rho(t)]-\sum_n\big[|{\pi_{n}}\rangle\langle{\pi_{n}}|,\rho(t)\,\bar{O}_0^{(n)\,\dagger}(t)\big]-\sum_n\big[\bar{O}_0^{(n)}(t)\,\rho(t), |{\pi_{n}}\rangle\langle{\pi_{n}}|\big]. \end{eqnarray}\noindent \tag{ 9 }$$

Here the operator

$$\begin{equation} \bar{O}_0^{(n)}(t)=\int_0^t \mathrm{d}s\ \alpha_n(t-s)O^{(n)}_0(t,s) \end{equation}\noindent \tag{ 10 }$$

is determined by the auxiliary equation

$$\begin{eqnarray} &&\partial_tO_0^{(n)}(t,s)=-\mathrm{i}[H_{\rm sys},O^{(n)}_0(t,s)]+\sum_m[|{\pi_{m}}\rangle\langle{\pi_{m}}|\bar{O}^{(m)}_0(t),O^{(n)}_0(t,s)] \end{eqnarray}\noindent \tag{ 11 }$$

with the initial condition O⁽ⁿ⁾₀(s,s) = −|π_n〉〈π_n|. We use the same notation as in [48] where the subscript 0 indicates that the operators are obtained from the ZOFE approximation. Note that, while the density matrix obtained from the stochastic Schrödinger equation is always positive, positivity of the density operator of equation (9) is not guaranteed.

For the numerical implementation, we approximate the environment correlation function as the sum of exponentials, i.e. $\alpha _n(\tau )=\sum _j^M\! p_{nj}{\rm e}^{\mathrm {i} \Omega _{nj}\tau }$ 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

$$\begin{eqnarray} \partial_t\bar{O}_0^{(nj)}(t)& =& -\mathrm{i}[H_{\rm sys},\bar{O}_0^{(nj)}(t)]\nonumber\\ &&-p_{nj}|{\pi_{n}}\rangle\langle{\pi_{n}}|+\mathrm{i}\,\Omega_{nj}\bar{O}_0^{(nj)}(t)\nonumber\\ &&+\sum_m[|{\pi_{m}}\rangle\langle{\pi_{m}}|\bar{O}^{(m)}_0(t),\bar{O}_0^{(nj)}(t)] \end{eqnarray} \tag{ 12 }$$

for the operators $\bar {O}^{(n)}_0(t)=\sum _j\!\bar {O}^{(nj)}_0(t)$ with the initial condition $\bar {O}_0^{(nj)}(t=0)=0$. In the numerical implementation we expand the above equations in the |π_n〉 basis. To determine Ω_nj, we proceed similarly as in [54, 55]. However, in the present work we do not use the Matsubara decomposition for the coth appearing in the correlation function α_n(τ) (see equation (6)) but use the partial fraction decomposition proposed by Croy and Saalmann [56]. 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.

3.1. Comparison with hierarchical equations of motion calculations

First, we will compare the energy transfer calculated using the ZOFE master equation approach with the HEOM results of [35]. For the comparison, we restrict ourselves to the same model as that used in [35], where only a single monomer unit of the FMO trimer was considered and the recently discovered eighth BChl was not taken into account (see the Hamiltonian in table 1).

Table 1. Matrix elements of the Hamiltonian H_sys of the FMO monomer in cm⁻¹. The upper triangle (couplings) is the same as the lower triangle (not shown here). From the energies on the diagonal, an offset of 12 000 cm⁻¹ is subtracted (for the calculation of the transfer, only the energy differences are relevant). The Hamiltonian is that used in [35] with the energies and couplings from [20].

	1	2	3	4	5	6	7
1	410	−87.7	5.5	−5.9	6.7	−13.7	−9.9
2		530	30.8	8.2	0.7	11.8	4.3
3			210	−53.5	−2.2	−9.6	6.0
4				320	−70.7	−17.0	−63.3
5					480	81.1	−1.3
6						630	39.7
7							440

For the comparison we have taken the transition energies and couplings from [20]. In [35], a Drude–Lorentz (D–L) SD was used to describe the coupling to the environment. For the decomposition of the environment correlation function as the sum of exponentials, described above, we represent the SD by the sum of anti-symmetrized Lorentzians. With ten anti-symmetrized Lorentzians, we could fit the D–L SD of [35] perfectly in the relevant energy range. The resulting SD is shown as the solid curve in figure 2(a).

Zoom In Zoom Out Reset image size

The calculated transfer resulting from this SD is shown in figure 3. The curves show the time-dependent excitation probabilities of the individual BChls—the dashed lines are the HEOM results taken from [35] and the solid lines show the results calculated with the ZOFE master equation (see equation (9)). In the upper panels, the transfer is calculated for a temperature of 77 K, and in the lower panels, for 300 K. Initially all excitation is either on BChl 1 (left panels) or on BChl 6 (right panels). In all the 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 77 K.

Zoom In Zoom Out Reset image size

3.2. Influence of high-energy modes

Since the ZOFE method enables us to consider an almost arbitrary form of the SD, we can now investigate the influence of different contributions to the SD.

Often, the so-called reorganization energy

$$\begin{equation} \lambda=\int_0^{\infty}\mathrm{d}\omega\, J(\omega)/\omega \end{equation}\noindent \tag{ 13 }$$

is used as a global measure for the coupling strength to the (environmental) vibrational modes depending on the SD J(ω) [8]. 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 SD 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 SD. As an example we consider an SD (dotted curve in figure 2(a)), which coincides with the D–L SD of [35] for low energies, but has a large additional peak centered at about 1600 cm⁻¹, resulting in strong additional coupling to high-energy modes. We have calculated the transfer through the FMO monomer for this modified SD using the ZOFE master equation. The resulting transfer is exactly the same as that calculated by the ZOFE method for the original SD of [35] shown in figure 3. This shows that despite contributing significantly to the reorganization energy, the high-energy part above 500 cm⁻¹ (for up to 500 cm⁻¹ the spectral densities are the same) is of no relevance to the system dynamics. The reason 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 figure 2), are in the energy region below 500 cm⁻¹. This means that the detuning between these system energies and the high-energy part of the SD 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 into account this fact and thus gives, in this case, a wrong indication of the strength of the coupling to the modes.

Similarly, a faster drop-off, as given by the SD shown by the curve (triangles) in figure 2(a), does not affect the transfer dynamics, since in the region below 500 cm⁻¹ it also agrees with the original SD.

3.3. Influence of low-energy modes

The spectral densities we have considered so far have a smooth form without any prominent structure. However, a more realistic SD consists of several distinct peaks that embody the strong coupling to the intra-molecular modes of the BChls, and is taken, e.g., from measured fluorescence line narrowing spectra [20, 32] or from atomistic simulations (see, e.g., [25, 31]).

In the following, we will investigate how structures in the SD influence the excitation transfer dynamics. As an example, we consider an SD 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 we are confident from previous investigations [48] 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⁻¹) constant. To this end we first introduce an effective reorganization energy (ERE)

$$\begin{equation} \lambda_{\rm eff}=\int_{E_{\rm min}}^{E_{\rm max}}\mathrm{d}\omega\ J(\omega)/\omega, \end{equation}\noindent \tag{ 14 }$$

where we take E_min = 0 cm⁻¹ and E_max = 550 cm⁻¹ as the relevant energy range that contains all electronic transition energies of the FMO monomer (see the vertical lines in figure 2)⁵.

For the D–L SD in figure 2(a), we get a value of λ_eff = 31 cm⁻¹ for the ERE.

To investigate the influence of the structure of the SD, we consider three different spectral densities that all have the same ERE as the D–L SD. But in contrast to the D–L SD, these spectral densities consist of four prominent peaks whose intensities are varied. These spectral densities are shown as solid curves in figures 2(b)–(d). As shown in these figures, one of the peaks is enhanced at a time. For comparison again the D–L SD is shown in each figure as the dotted curve. Note that due to the division by ω in the definition of the reorganization energy (see equation (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 on the SD, we calculated the excitation transfer through the FMO monomer using the ZOFE master equation. The resulting time-dependent excitation probabilities are shown in figure 4 in the same way as in figure 3, but now for the three different spectral densities of figures 2(b)–(d) and the D–L SD.

Zoom In Zoom Out Reset image size

One sees that although the overall coupling strength in the relevant energy region is kept constant, the transfer changes clearly when the structure of the SD is varied. In figure 4(a) for instance, the population on BChl 3 after 1 ps changes by up to ∼40% depending on the structure of the SD.

Obviously, not only the overall coupling strength in the relevant energy region is important, but also the local distribution of the coupling strength.

As the temperature is increased (lower panels), the spreading of the curves decreases, showing that the structure of the SD becomes less important at higher temperatures.

In this paper, we have discussed the excitation energy transfer in the FMO complex, emphasizing the role of the coupling of the excitation to vibrational modes of BChl molecules.

The ZOFE master equation: To handle the coupling of the electronic excitation to vibrational modes, numerical calculations have been performed using the ZOFE master equation, utilizing the ZOFE approximation, which allows us to calculate the transport for the FMO complex within about 1 min on a standard PC. We validated the applicability of the ZOFE approximation by comparing with the exact HEOM results of [35], where a D–L SD has been used.

Reorganization energy: Since we do not restrict ourselves to a certain 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⁻¹) 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.

Structured SD: Being aware of this, we used an 'ERE', which is a local measure in the energy region of system transitions⁶, to investigate the dependence of the transfer on the local structure of the SD. We found 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 by up to ∼40%. This influence of the structure of the SD decreases when the temperature is increased.

Electronic energies and couplings: Note that the transfer, of course, depends strongly also 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 to not overestimate theoretical results obtained for a certain set of parameters.

Outlook: Since measured spectral densities show many peaks in the energy region where electronic transition energies are located [32], 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 [57]. One goal of the present study was to further investigate the rather abstract ZOFE approximation [42, 48]. The quite good agreement with the exact HEOM results of [35] makes us confident that it is worthwhile to further explore the NMQSD method. In a next step we plan to investigate the next order in the functional expansion. We further note that in the present work the ZOFE approximation was derived in the basis of localized states, but a formulation in the basis of delocalized electronic eigenstates is also possible. To compare with the HEOM results, we restricted ourselves to the coupling of the local excitations to a local environment. It is straightforward to include also the coupling of the off-diagonal elements to an environment.

The NMQSD approach is also suitable for describing 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⁻¹ does not influence the energy transfer markedly, in molecular aggregates of organic dyes these modes have a strong impact on the electronic excitation dynamics [58]. Furthermore, they have energies comparable to the dipole–dipole interaction which influences strongly e.g. the aggregate absorption spectra [59]. In previous investigations, we found that the typical 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 [47].

The NMQSD method is also readily applicable to treat optical spectroscopy. In particular, it should be well suited for the investigation of the influence of the coupling to vibrational modes on the signals observed in two-dimensional (2D) experiments. Such a study might be revealing, since the 'power spectrum' obtained from 2D measurements in [13] shows strong resemblance to the fluorescence line narrowing spectrum of Wendling et al [32].

To obtain further information on how individual vibrational modes influence the transport, it might be useful to simulate the transport using networks consisting of quantum [60] or classical oscillators [61].

Financial support from the DFG under contract no. Ei 872/1-1 is acknowledged. AE thanks Alán Aspuru-Guzik for 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.