Efficient energy transfer in light-harvesting systems, I: optimal temperature, reorganization energy and spatial–temporal correlations

Understanding the mechanisms of efficient and robust energy transfer in light-harvesting systems provides new insights for the optimal design of artificial systems. In this paper, we use the Fenna–Matthews–Olson (FMO) protein complex and phycocyanin 645 (PC 645) to explore the general dependence on physical parameters that help maximize the efficiency and maintain its stability. With the Haken–Strobl model, the maximal energy transfer efficiency (ETE) is achieved under an intermediate optimal value of dephasing rate. To avoid the infinite temperature assumption in the Haken–Strobl model and the failure of the Redfield equation in predicting the Forster rate behavior, we use the generalized Bloch–Redfield (GBR) equation approach to correctly describe dissipative exciton dynamics, and we find that maximal ETE can be achieved under various physical conditions, including temperature, reorganization energy and spatial–temporal correlations in noise. We also identify regimes of reorganization energy where the ETE changes monotonically with temperature or spatial correlation and therefore cannot be optimized with respect to these two variables.


3
This paper is the first of a two-part series, with paper I on optimization of ETE and paper II on mechanisms and network kinetics. Although intrinsically quantum mechanical, exciton dynamics are often described by site-to-site hopping as in random walk [21], and quantum coherence is often related to delocalized tunneling [22]- [24]. The kinetic expansion unifies coherent tunneling and incoherent hopping by systematically reducing exciton dynamics to network kinetics [16]. In paper II, this approach will be applied to analyze network structures of light-harvesting systems and to predict the contribution of quantum coherence [25].
The paper is organized as follows. In section 1, we introduce the exciton dynamics model and define the ETE for light-harvesting energy transfer. In section 2, we apply the Haken-Strobl model to investigate the optimal ETE in FMO as a function of the dephasing rate. In contrast to earlier studies, our study emphasizes the initial condition dependence, the approximation of the efficiency with the average trapping time and the secular approximation. In section 3, we apply the generalized Bloch-Redfield (GBR) equation to explore the detailed energy transfer optimization conditions in FMO with reorganization energy, temperature and spatial-temporal correlations. In section 4, we calculate the optimal ETE in another light-harvesting system, PC 645, with the GBR approach. The general optimization feature of the energy transfer process is illustrated by our studies. In section 5, we discuss our conclusions.

Exciton dynamics in light-harvesting systems and energy transfer efficiency (ETE)
When sunlight falls on light-harvesting pigments, an absorbed photon can excite the photosynthetic system. The excitation energy is then transported through an energy transfer network to the reaction center for subsequent charge separation, resulting in energy trapping. During transport, the energy can decay when excitation energy is lost as heat via irreversible electron-hole recombination and can be redistributed through interaction with protein environments. Hence, the exciton dynamics for the light-harvesting system follows the Liouville equation [16], [26]- [28], ∂ t ρ = −Lρ = −[L sys + L trap + L decay + L dissip ]ρ, (1) where ρ is the reduced density matrix of the exciton system, and each term of the Liouville superoperator L describes a distinct dynamic process.
The quantum coherent evolution, L sys ρ = i[H, ρ]/h, is controlled by the system Hamiltonian, which is given by H nm = δ n,m ε n + (1 − δ n,m )J nm from the tight-binding model in the local site basis set representation [27]. The diagonal elements, ε n , define the site energies, whereas the off-diagonal elements, J nm = J * nm , define the dipole-dipole interaction coupling strength between two distinct sites.
The two irreversible energy decay channels are exciton decay and trapping. The electron-hole recombination process is described by the decay term [L decay ] nm = (k d;m + k d;n )/2, where k d;m is the decay rate at site m, and L nm represents the diagonal element of the Liouville operator, L nm = L nm,nm . Similarly, localization at the charge separation state is described by the exciton trapping term, [L trap ] nm = (k t;m + k t;n )/2, where k t;m is the trapping rate at site m.
The system-bath interaction, H SB = m Q m B m , is applied to describe exciton dissipative dynamics, where Q m = |m m| and B m are system and bath quantum operators, respectively [26]- [28]. The bath influence is characterized by the time-correlation function, C mn (t) = B m (t)B n , which is related to the spectral density J mn (ω) by C mn (t) = ∞ 0 [coth(hβω/2) cos(ωt) − i sin(ωt)]J mn (ω)dω. Alternatively, we can introduce a timedependent site energy fluctuation δε m (t). In the Haken-Strobl model, a spatially uncorrelated 4 classical white noise follows δε m (t) = 0 and δε m (t)δε n = * δ(t)δ m.n , which is valid in the infinite temperature limit [29]. In the local site basis set representation, dissipation becomes simply decoherence, and exciton dynamics can be solved exactly by the second-order expansion form [22,29]. At finite temperatures, a theoretical framework of quantum dissipative dynamics differing from the Haken-Strobl model is required to include detailed balance and the memory effect in slow bath relaxation. In the weak dissipation regime, the Redfield equation (with or without the secular approximation) [30] has been widely used in modeling exciton dynamics, e.g. several previous studies exploring ETE [31,32]. However, this approach can lead to unphysical predictions in the intermediate to strong dissipation regime. With the introduction of auxiliary fields, exciton dyanamics can be studied in the GBR equation approach, which is more reliable over a broad regime of dissipation and predicts the correct Forster rate limit [33].
In this paper, we first use the Haken-Strobl model to study the optimization of the exciton energy transfer (EET), obtaining a simple physical picture. We will then use the GBR approach with the Debye spectral density to explore optimal EET conditions with various control parameters, including temperature, reorganization energy and spatial-temporal correlations of the bath.
The energy trapped at the reaction center, i.e. the rate process described by L trap , represents effective energy transfer, whereas the rate process described by L decay represents ineffective energy loss. The ETE is defined by the branching ratio of the energy trapping process [4,13,34], i.e. q = Tr ∞ 0 L trap ρ(t) dt = n k l,n τ n n k l,n τ n + n k d,n τ n , where the trace denotes the sum of diagonal population elements of the matrix. This definition can be simplified to q = n k t,n τ n because the total depletion probability in the denominator is normalized to one. The mean residence time at each site of the exciton system is τ n = ∞ 0 dtρ n (t), where the population ρ n is represented as the diagonal elements on the density matrix ρ n = ρ n,n . To simplify, we assume that k d is identical at each site together with a necessary condition k t k d for the high-efficient EET network. Then, the k d dependence of the residence time becomes negligible, giving Here, the residence time is approximately τ (k d ) ≈ τ (0), with the normalization condition n k t,n τ n (k d = 0) = 1, and t = n τ n (k d = 0) is the mean first passage time to the trap state in the absence of decay (i.e. the average trapping time). As shown later in this paper, the optimal ETE is determined reliably by the minimal trapping time. The above definitions follow closely those introduced in [16]. Using the stationary solution to equation (1), L ∞ 0 ρ(t)dt = Lτ = ρ(0), we obtain the average trapping time as where ρ(0) is the initial condition.
The experimental Hamiltonian has the following matrix elements (in cm −1 ) [37,38]: The site closest to the reaction center, BChl 3, is designated as the only site with a nonzero trapping rate k t =1ps −1 ([13]- [15]andseefootnote2).Tocomparewithrealexperiments,the initial population is considered to be at either BChl 1 or BChl 6, i.e. ρ m=1 (0) = 1 or ρ m=6 (0) = 1. For optimal ETE, the exciton dynamics of the FMO complex uses environmental noise toincreasethetransferefficiency([13]- [15]andseefootnote2),previouslystudiedusingthe Haken-Strobl model. In this section, we will revisit this simple model but emphasize three new aspects of the EET process in FMO: (i) the effect of different initial conditions, (ii) verification of approximating the ETE with the average trapping time and (iii) limitation of the secular approximation in the Redfield equation.
We calculate the average transfer time t as a function of the pure dephasing rate * with the two different initial conditions. As shown in figure 1(a), a characteristic optimization phenomenon of t is observed: at dephasing rates less than 50 cm −1 , the average trapping time has a steep increase as * decreases; at dephasing rates more than 200 cm −1 , the average time marginally increases with * . The interplay between coherent population oscillation and dissipative population redistribution leads to an optimal transfer efficiency at an intermediate dephasing rate. The explicit optimal values are t min = 10.1 ps with * opt = 195 cm −1 and the initial condition defined at BChl 1, and t min = 8.65 ps with * opt = 175 cm −1 and the initial condition defined at BChl 6. Different initial conditions result in similar but distinguishable results, which can be interpreted by a simple physical argument. In FMO, there are two dominating EET pathways: 1 → 2 → 3 and 6 → (5, 7) → 4 → 3, where each number represents a specific BChl site [39]. With fewer sites involved, the first pathway exhibits stronger quantum coherence than the second pathway. To compensate the inefficient energy transfer in the fully coherent limit, a larger optimal dephasing rate is needed for the first pathway than the second pathway, i.e. * opt (Bchl 1) > * opt (Bchl 6). The physical quantity of interest is quantum efficiency, whereas the theoretical quantity computed in this and our preceding papers is the average trapping time. The approximate expression in (4) establishes the connection between these two quantities. We examine the validity of this approximation by evaluating the ETE in FMO with the exciton recombination rate k d =1ns −1 ([13]- [15]andseefootnote2).Asshowninfigure1(b),theaveragetrapping time provides a reliable measurement for the ETE as * changes by five orders of magnitude (0.1 cm −1 < * < 10 4 cm −1 ), except for a slight underestimation in the weak dephasing limit ( * < 0.1 cm −1 ). Thus, we will focus on the calculation of t in the rest of the paper.  (4). Here, the trapping rate is k t = 1 ps −1 and the decay rate is k d = 1 ns −1 .
We now examine the failure of the secular approximation in the intermediate to strong damping regime. As demonstrated in our forthcoming paper, a physically reasonable description of dissipation in energy transfer networks will always lead to optimal noise for the maximal efficiency [17]. Here the optimal noise refers to a finite value of noise where the ETE reaches the maximal. We present here a counter-example, where the secular approximations adopted in the Redfield equation lead to an unphysical prediction in the average trapping time. The Haken-Strobl model introduces classical Gaussian fluctuations in the site energy and therefore is rigorously described by the dephasing operator in (2) in the local basis set [22,29]. A unitary basis set transformation of the Liouville equation in (1) does not change the exciton dynamics and will give exactly the same prediction, which is confirmed by our calculation using the full Redfied equation in the eigen-exciton basis set [30]. In contrast, as shown in figure 1(a), an additional secular approximation in the Redfield equation leads to a plateau of the average trapping time as the dephasing rate increases. In the secular Redfield equation, the relaxation rate is proportional to * , suggesting instantaneous population transfer to the trap state (e.g. BChl 3 in FMO) for large * . Therefore, the plateau in figure 1 is given by the trapping rate, i.e. t ∝ 1/ * + 1/k t . However, strong dephasing should always lead to classical behavior, i.e. hopping kinetics, which results in a slow-down in the diffusion to the trap state. As in the diffusioncontrolled reaction, diffusion becomes the rate-limiting step in the intermediate to strong dissipative regime, giving t ∝ * . The failure of the secular approximation has been discussed in several contexts before [41,42], and is now used to explain the observed discrepancy between the Haken-Strobl model calculation [10], [13]- [15] and secular Redfield calculation 7 (i.e. the quantum jump method [32]). Although our discussion here is limited to the Haken-Strobl model defined for infinite temperature, we show in the next section that the conclusion about the secular approximation remains valid for exciton dynamics at finite temperatures.
Noise-assisted energy transfer using the Haken-Strobl model has been reported recently byseveralgroups( [10], [13]- [15]andseefootnote2),butabasicquestionremainstobe answered: is the optimal dephasing rate observed in FMO a general rule for realistic lightharvesting systems or a special case obtained from model exciton systems? Without including temperature and memory effects, predictions based on the Haken-Strobl model may not be relevant under realistic conditions. Although previous studies have attempted to include temperature dependence and realistic spectral densities, the theoretical methods involved can lead to qualitatively unphysical predictions in the intermediate and strong dissipation regimes [31,32]. In the next section, our study at finite temperatures will demonstrate the general optimal conditions in FMO, including temperature, reorganization energy and spatial-temporal correlation.

GBR equation
In the previous section, we used the Haken-Strobl model, defined in the infinite temperature limit, to illustrate the possibility of optimization in the EET process. At finite temperatures, forward and backward energy transfer rates need to satisfy the detailed balance condition. In the weak dissipation regime, the Redfield approach is able to reliably describe exciton dissipative dynamics. With an increase in system-bath coupling and the slow-down of bath relaxation, more expensive methods such as the hierarchic expansion for the Gaussian bath are required for describing exciton dynamics accurately [43]- [46]. Although the hierarchic approach has been applied to FMO, a rigorous investigation of quantum dissipative dynamics in large-scale exciton systems can be numerically difficult. To capture the relevant optimization feature of EET dynamics in a qualitatively reliable way, we will apply a GBR equation approach, derived from the second-order cumulant expansion. For many bath spectral densities such as Debye or Ohmic, the corresponding time correlation functions can be (numerically) expanded using exponential functions [45], where c mn is the spatial correlation coefficient between sites m and n, v i is the relaxation rate of the ith bath mode and f r i ( f i i ) is the real (imaginary) part of the expansion coefficient. In the frequency domain, the expansion is applied to both spectral density J (ω) and temperature dependence coth(hω/2k B T ). With the facilitation of auxiliary fields g m;i (t) at site m, the GBR equation for exciton dynamics is written aṡ where Q m = |m m| is the system operator together with m = n c mn Q n , and [A, B] + = AB + B A is the anti-commutator. The initial value of g m;i (t) is zero. Equation (8) is a generalization of its original form [33]. Here, the memory effect in the dissipative dynamics due to the interaction with the bath is represented by the auxiliary field, which can be considered as an element of the projection operator Q(I -P) in the Nakajima-Zwanzig projection operator technique.
In this paper, we will use the spatially correlated Debye spectral density, J mn (ω) = c mn J (ω) and J (ω) = (2h/π )λωD/(ω 2 + D 2 ), where the reorganization energy λ represents the system-bath coupling strength and the Debye frequency D is the bath relaxation rate (the inversion of bath temporal correlation). The choice of spectral density will not affect our general conclusion. Following the Matsubara expansion [26], the bath time correlation function can be written explicitly as the exponential expansion in (7), where v i = 2πik B T /h is the ith (i > 0) Matsbura frequency and v 0 = D. Dissipation is thus characterized by four parameters: reorganization energy λ, temperature T , bath relaxation rate D and bath spatial correlation c mn .
Due to the equivalence of the GBR equation and the second-order time-nonlocal expansion form, the trapping time can be alternatively calculated by the Laplace transform in the eigenexciton basis set using spectral density. The results from two approaches are exactly the same, which demonstrates the basis set invariance of the GBR approach for dissipative dynamics in EET systems.

Optimization with reorganization energy
We begin with studying the influence of reorganization energy on the EET process. As an example, we assume zero spatial correlation c mn = δ m,n , an experimentally reasonable temporal correlation D −1 = 50 fs, and room temperature T = 300 K. The trapping time calculated from the GBR approach is plotted in figure 2(a). As λ increases, the trapping time quickly decreases from the fully coherent value by two orders of magnitude to the minimal value, corresponding to the optimal ETE, and then slowly increases again. The explicit optimal values are t min = 3.39 ps with λ opt = 92 cm −1 for the initial population at BChl 1 and t min = 3.27 ps with λ opt = 73 cm −1 for the initial population at BChl 6, respectively. These results, including the initial condition dependence, are similar to those derived in the Haken-Strobl model, except that the reorganization energy λ is used instead of the dephasing rate * (see figure 1(a)). At high temperatures, the dephasing rate can be effectively written as * which is the zero-frequency limit of J (ω)coth(hω/2k B T ). The linear relation between * and λ indicates that these two parameters behave qualitatively in the same way for the EET process. As we will find later, equation (9) can also be used to interpret the optimization behavior with other control parameters.
In figure 2(b), we compare the predictions of t from several master equation approaches. Both the full Redfield equation and its secular form lead to an unphysical plateau of t for large λ. This plateau was observed in recent calculations of FMO [31,32] and, as discussed in the previous section, is due to the time-local approximation of the Redfield equation. In the third approach, we only keep the real part of Redfield tensors, which is equivalent to the Haken-Strobl model but with temporally correlated noise. This extended Haken-Strobl model recovers the optimization behavior but overestimates noise enhancement: the predicted λ opt is 50% smaller, whereas t min is 50% larger than those from the GBR approach. These observations are consistent with an early study by Ishizaki and Fleming in the two-site model system [44,45]. On the other hand, with a time-nonlocal form for dissipative dynamics, our GBR approach is qualitatively accurate in calculating the variation of t with λ. A comparison with the sophisticated hierarchic approach will be necessary to test the quantitative accuracy of results in this paper. However, from the practical point of view, the GBR serves as a simple but unifying approach to investigate exciton dynamics in a broad range of parameter spaces with low computational cost.

Optimization with temperature
Next we study the dependence of the trapping time on temperature. With the same spatial-temporal correlation (c mn = δ m,n and D −1 = 50 fs) used in the previous subsection, we fix the reorganization energy at λ = 35 cm −1 , which is a physically reasonable estimation from experiments [46]. The results under the two initial conditions are reported in figure 3. With the initial condition at BChl 1, the minimal trapping time is t min = 4.2 ps at optimal temperature T opt = 162 K, but the increase in t with decreasing T at low temperatures is weak. With the initial condition at Bchl 6, the trapping time monotonically increases with temperature, and the maximum ETE is at zero temperature. Since the trapping time is kept in the order of ps, FMO shows a strong robustness against the change of temperature.
In dissipative dynamics, the exact temperature dependence is determined by coth(hω/2k B T ), which changes from a linear function, 2k B T /hω, at high temperatures to a constant at low temperatures. Nonvanishing dissipation is needed to maintain the equilibrium energy distribution even at T → 0. Because dissipation always increases with temperature for a given spectral density, the zero-temperature limit of dissipation determines the general behavior of T dependence of the trapping time and the ETE. When the zero-temperature limit of dissipation is weaker than the optimal value, increasing the temperature can increase the dissipation strength and therefore help the system to achieve the maximum ETE at a finite temperature. When the zero-temperature limit of dissipation is stronger than the optimal value, the trapping time increases monotonically with temperature and the system will be further away from the optimal condition. Hence, depending on the reorganization energy, temperature can behave in two different ways in the EET process. To verify this prediction, we use a small reorganization energy of λ = 20 cm −1 , and start the initial population at BChl 6. The inset of figure 3 presents an optimization curve where T opt = 46 K is observed.

Optimization with temporal correlation
In the Markovian approximation for dissipative dynamics, bath relaxation is assumed to be much faster than system dynamics. This assumption may not be reliable in light-harvesting systems due to the complexity of the protein scaffold. The time scale of bath relaxation (D −1 ) in FMO is usually of the order of 10-100 fs [46], smaller but comparable to that of exciton dynamics. The non-Markovian memory effect can enhance quantum coherence up to 500 fs as observed experimentally [5,7]. As a time-integrated effect of exciton dynamics, ETE can depend strongly on the temporal correlation of the bath, which is shown by the D dependence of the effective dephasing rate at high temperatures in equation (9).
Here we apply zero spatial correlation c mn = δ m,n , reorganization energy λ = 35 cm −1 and room temperature T = 300 K. Under the two initial conditions, the trapping time t is calculated as a function of the Debye frequency D in the GBR approach. As shown in figure 4, the temporal correlation of the bath leads to optimal efficiency at an intermediate level of D. In our calculation, the optimal conditions are t min = 4.48 ps with D −1 opt = 32 fs and the initial condition at BChl 1, and t min = 3.41 ps with D −1 opt = 41 fs and the initial condition at BChl 6. Qualitatively, the influence of the Debye frequency on the EET process can be interpreted by the inverse relation between the effective dephasing rate * and D in equation (9). Our results further suggest that a more quantitative characterization of bath temporal correlation, rather than being labeled as Markovian or non-Markovian, is necessary for studying the EET process.

Optimization with spatial correlation
Recent experimental and theoretical studies have suggested that strong spatial correlation can be relevant in exciton dynamics in light-harvesting systems [6]. Spatial correlation can be quantified through different approaches. Here, we use an exponentially decaying function [32,40,47], to define the spatial correlation with the site-site distance R mn and the correlation length R 0 [6]. Other functions, such as Bessel functions, have also been used to describe the spatial correlation [48]. With the assumption that R 0 in FMO is the same as that in the reaction center, we rewrite equation (10) as where c = 0.9 and R RC = 11.93 Å are the coefficient of the spatial correlation and the site-site distance in the reaction center, respectively. Instead of the reported value for c, we use c as the free parameter, ranging from 0 to 1, to characterize spatial correlation.
To compare with experiments, we use realistic parameters, reorganization energy λ = 35 cm −1 , temporal correlation D −1 = 50 fs and room temperature T = 300 K, to calculate the dependence of the trapping time t on the spatial correlation coefficient c. For simplicity, we present the result with the initial population localized at BChl 1 (the result with BChl 6 as the initial condition is similar). As shown in figure 5(a), the average trapping time monotonically increases with spatial correlation. The suppression of transfer efficiency due to spatial correlation has been shown in an early numerical simulation [47]. However, when we use a much larger reorganization energy, λ = 150 cm −1 , the EET process reaches maximal efficiency with an intermediate level of spatial correlation, c opt = 0.66. To understand these two different behaviors, we note that a positive spatial correlation always reduces the dissipation strength. For example, for a quantum two-level system, the energy fluctuations between the two levels can be characterized by an effective reorganization energy defined as λ eff = λ(1 − c). In the completely correlated case (c = 1), quantum coherence between these two levels can persist under an arbitrarily strong energy fluctuation (λ eff = 0). Therefore, when the reorganization energy is smaller than the optimal value, spatial correlation further reduces the effective reorganization energy and there is no optimal ETE as a function of c, i.e. the ETE is maximized at c = 0. In contrast, when the reorganization energy is larger than the optimal value, the ETE can be maximized at a nonzero c by the reduced dissipation at an optimal value of spatial correlation. Indeed, as shown in figure 5(b), the nonzero optimal spatial correlation starts appearing for λ > 60 cm −1 with the initial condition at BChl 1.

Optimization in PC 645
In this section, we consider another light-harvesting system, PC 645, in cryptophyte algae. Exciton dynamics of PC 645 can be studied using an effective eight-site Hamiltonian model, where the excitation energy is transferred from a center dimer dihyrobiliverdin (DBV) β50/β61 with the highest energy to two phycocyanins (PCB) β82 with the lowest energy during the initial step [8,49,50]. To compare with experiments where the initial state deviates from an eigenstate [8], we apply an incoherent initial condition at two DBV bilins with equal population, and the trapping rates k t = 1 ps −1 at two PCB 82 bilins. Hence, we limit our investigation to the initial step from DBV bilins to PCB bilins within the time scale of ps. As in the study in the previous section, the GBR approach is applied to explore the optimal ETE and its dependence on various control parameters. Experimentally, the bath spectral density is approximated by a two-frequency Debye model with D −1 1 = 50 fs and D −1 2 = 1.5 fs [50]. The reorganization energies for these two frequencies are chosen to be the same, giving J mn (ω) = c mn [J 1 (ω) + J 2 (ω)] with J i=1,2 (ω) = (2h/π)λωD i /(ω 2 + D 2 i ). Here, we adopt the Gaussian bath model with these two bath time scales and investigate the influence of varying reorganization energy, temperature and spatial correlation on the EET process.
Without spatial correlation, we calculate the trapping time t as a function of reorganization energy λ at room temperature T = 300 K; the result is plotted as a solid line in figure 6(a). Similar to the case in FMO, the trapping time is large in the regimes of small and large reorganization energies, but the reorganization energy dependence of ETE at intermediate λ becomes complicated. Two local minima appear in the trapping time due to the fact that the two Debye frequencies are highly separated in time scales. The explicit values of these two local minima are t min = 4.53 ps with λ opt = 83 cm −1 and t min = 4.84 ps with λ opt = 193 cm −1 . Thus, the optimization behavior is also dependent on the detailed structure of the bath. In comparison, we also present the result with consideration of spatial correlation. Following the original exponentially decaying function in equation (10), the spatial correlation is characterized by the correlation length R 0 . The result of the trapping time at R 0 = 20 Å is plotted as a dashed line in figure 6(a), which also exhibits a range of intermediate λ for minimized trapping time and optimized ETE, but the values of λ opt become larger than those without spatial correlation. The increase in λ opt arises from the reduction of dissipation strength due to spatial correlation.
Using the reorganization energy, λ = 135 cm −1 [50], we next calculate the trapping time as a function of temperature, with and without spatial correlation. As shown in figure 6(b), t can be minimized at an intermediate temperature, t min = 2.95 ps with T = 167 K and zero spatial correlation, and t min = 3.78 ps with T = 199 K and R 0 = 20 Å. The GBR approach in PC 645 becomes less reliable at low temperatures due to much larger reorganization energy than that in FMO. Our study is thus limited to a physically acceptable regime of the GBR approach with T > 150 K, below which the inversion of the Liouville superoperator may become unreliable. The value of the cutoff temperature, 150 K, is chosen for numerical convenience.
In our final calculation, we study the dependence of the trapping time on spatial correlation (using the correlation length). With temperature T = 300 K and reorganization energy λ = 135 cm −1 , the result is plotted in figure 7. The maximum ETE is observed with the minimum trapping time at R 0;opt = 27 Å. Since the interaction between the chromophores of PC 645 and the protein backbone is covalent, the reorganization energy is much larger than that of FMO and the optimization on spatial correlation becomes more pronounced, as discussed in the previous section. In addition, the dependence of the trapping time on spatial correlation becomes complicated with a local maximum trapping time appearing around R 0 = 10.4 Å.

Conclusions and discussions
This paper is based on recent efforts at exploring quantum effects in light-harvesting systems for the purpose of optimizing ETE. Employing the analytical methods and concepts introduced in a recent review article [16], our comprehensive study of FMO and PC 645 with both the Haken-Strobl model and the GBR equation reveals a general feature of the EET process: the interplay of coherent dynamics and environmental noise leads to optimal ETE at an intermediate level for various variables. Explicitly, the reorganization energy and bath relaxation rate, i.e. Debye frequency, yield nonmonotonic dependence and thus lead to optimal ETE. On the other hand, a nonzero value of temperature can optimize the ETE only when the dissipation is weak (less than its optimal value); a nonzero value of spatial correlation can optimize the ETE only when the dissipation is strong (large reorganization energy). If the dissipation strength does not fall in the chosen regime, temperature and spatial correlation can lead to monotonic changes in ETE and we will not observe optimal ETE. This observation clearly demonstrates a difference between spatial correlation and temporal correlation, although both correlations can enhance quantum coherence. Although the findings in this paper are reported for FMO and PC 645, the underlying principles should be general. In a forthcoming paper, we will identify the exciton states orthogonal to the trap state as the sufficient condition for nonzero optimal noise [17], which extends the invariant subspace proposed by Plenio et al [15] extracted from a specific fully connected network model. The orthogonal exciton subspace provides a unifying framework to explain optimal noise, initial preparation, coherent phase modulation and spatial-temporal correlations.
For light-harvesting systems, parameters under realistic situations are often not the same as their optimal values shown in our calculation. However, light-harvesting systems exhibit robustness against the change of variables since the resulting trapping time varies slowly and the time scales of energy trapping and decay are highly separated. For example, the variation in the trapping time at room temperature in FMO is within 1 ps as the reorganization energy changes from 50 to 150 cm −1 . A physical interpretation will be found based on the kinetic expansion discussed in our forthcoming paper [25].
An interesting observation is that the optimal conditions of different control parameters have roughly the same order of magnitude. For example, at room temperature (k B T ≈ 200 cm −1 ) and the experimental relaxation time scale D −1 = 50 fs (D = 106 cm −1 ), the optimal reorganization energy in FMO with zero spatial correlation is λ opt = 92 cm −1 (λ opt = 73 cm −1 ) for the initial condition at BChl 1 (BChl 6). All these parameters are approximately in the same energy scale, 100 cm −1 , which is also close to values of energy difference and electronic coupling between nearest-neighbor sites. In our forthcoming paper, based on two asymptotic behaviors of the trapping time in the regimes of weak and strong dissipations, we find that a crude estimation of the optimal dissipation strength is of the same order as site-site coupling [17].
Reliable predictions of the EET process in the light-harvesting systems require an adequate description of dissipative effects. For example, the Forster theory is applicable to incoherent hopping, whereas the Redfield equation is applicable to the Markovian limit of exciton dynamics [44,45]. Due to the nature of classical noise, the Haken-Strobl model is limited to the regime of high temperatures [22]. As derived from the second-order cumulant expansion, the GBR equation is shown in this paper to provide a reliable qualitative description of the trapping time in a broad range of dissipation and successfully predicts the optimization of the ETE with respect to different variables. However, higher-order cumulants become relevant for strong/slow dissipation and the corrections to the GBR approach are then necessary. The hierarchical approach with Gaussian noises, applied to FMO [46] and LH2 [51], can be used to improve the quantitative accuracy of our predictions. Alternatively, we can begin with incoherent exicton dynamics as a reference and then include quantum coherence systematically, which can be accomplished by the extension of the noninteracting-blip approximation (NIBA) derived in the two-site systems [52]- [55].