Iterative path-integral algorithm versus cumulant time-nonlocal master equation approach for dissipative biomolecular exciton transport

We determine the real-time quantum dynamics of a biomolecular donor–acceptor system in order to describe excitonic energy transfer in the presence of slow environmental Gaussian fluctuations. For this, we compare two different approaches. On the one hand, we use the numerically exact iterative quasi-adiabatic propagator path-integral scheme that incorporates all non-Markovian contributions. On the other, we apply the second-order cumulant time-nonlocal quantum master equation that includes non-Markovian effects. We show that both approaches yield coinciding results in the relevant crossover regime from weak to strong electronic couplings, displaying coherent as well as incoherent transitions.


Introduction
The photosynthetic conversion of the physical energy of sunlight into its chemical form suitable for cellular processes involves a complex series of physical and chemical mechanisms [1,2]. Photosynthesis starts with the absorption of a photon by a light-harvesting pigment and the formation of an exciton, followed by the transfer of the exciton to the reaction center, where charge separation is initiated. It is the nature of this transfer in the form of a cascade that has recently become a subject of intensive research. The experimental progress [3][4][5][6][7] in terms of ultrafast time-resolved optical spectroscopy allows reconsideration of the longstanding question [8] of whether the transfer of excitonic energy is coherent or incoherent. Put differently, does the excitation show signatures of quantum coherent wave-like propagation, or does it move via a hopping process, which may be described by a (classical) incoherent master equation? Recent experimental results [3][4][5][6][7] have provided an indication that the simple model of incoherent exciton hopping through the network of chromophore sites embedded in proteins is not sufficient. Instead, quantum coherence is necessary to describe long-lasting beating signals in two-dimensional electronic spectra [9] recorded from various photosynthetic systems.
Theoretically, photosynthetic excitonic energy transfer (EET) processes in light-harvesting complexes are often discussed in terms of simplified low-dimensional models [8] describing a few individual chromophore sites which mutually interact by dipolar couplings and which are exposed to the fluctuations of the solvent molecules and the protein host [2]. The picture that arises illustrates that long-lived electronic coherence at room temperature can be understood in terms of the constructive role of the low reorganization energies λ and environmental modes slower than the inverse of electronic couplings [10,11]. In the Fenna-Matthews-Olson (FMO) complex [1,2] isolated from green sulfur bacteria such as Chlorobaculum tepidum, for example, electronic couplings are J da 1-100 cm −1 , whereas the frequencies of the environmental fluctuations span a similar range [12][13][14][15]. Furthermore, Lee et al [4] have demonstrated that strongly spatially correlated fluctuations between two excitonic sites preserve the quantum coherence in a bacterial reaction center by applying a two-color electronic coherence photon echo technique. Since then, spatial correlations of environmental fluctuations have been under 3 intensive investigation [16][17][18][19][20][21][22]. Furthermore, recent investigations have addressed quantum entanglement in photosynthetic pigment-protein complexes [11,19], [23][24][25][26][27].
Decoherence and energy relaxation in the exciton transfer in photosynthetic antenna complexes have been studied theoretically using various models. Many studies use one of the various weak-coupling Markovian approaches. It has been shown recently that the weakcoupling Markovian methods such as the Redfield equation [28] fail qualitatively as well as quantitatively [10,29,30], since they assume that the environmental fluctuations readjust to their respective equilibrium much faster than the speed with which the coherent changes in the donor-acceptor systems occur. For the bioexcitonic energy transfer that typically occurs in a protein host, this assumption, however, does not hold. Traditional Förster theory [31] that is a good treatment for reorganization energies larger than the electronic coupling between pigments is also insufficient, since both energy scales are of the same order of magnitude. In addition, each pigment has its own protein environment, which yields to site-dependent reorganization energies and dynamics. It was pointed out that the site-dependent reorganization dynamics also plays a significant role as a physical origin of the long-lived quantum coherence [10].
Hence, techniques beyond the weak-coupling Markovian approaches are required. Related to exciton transfer in molecular complexes, extensions to non-Markovian fluctuations have been formulated in terms of generalized master equations [32], based on the stochastic Liouville equation [8], in the Lindblad approaches [33] and by small polaron approaches [34,35]. More recently, Nalbach and Thorwart [11,21,30] adopted the numerically exact quasi-adiabatic propagator path integral [36][37][38][39][40] to the particular situation of slow polarization fluctuations. Moreover, Ishizaki and Fleming applied the second-order cumulant time-nonlocal quantum master equation (2CTNL) approach [10] to the case of photosynthetic excitation energy transfer [14,19,23]. The purpose of this paper is to directly compare both methods in the parameter range relevant to the photosynthetic energy transfer, namely where the reorganization energy is of the order of the electronic coupling, site-dependent reorganization dynamics and long environmental relaxation time, i.e. of the order of the electronic transfer times ω −1 c = J −1 da . This allows us to combine the conceptual beauty of 2CTNL with the numerical exactness of the 'ab initio' path-integral approach. Hence, this paper is organized as follows. First, we introduce the model considered and briefly describe both methods. Then, in section 4, we compare the results of both methods for typical parameter sets of biomolecular exciton transfer dynamics, before we conclude with a short summary.

Model
We employ the simplest model of Frenkel excitons that transfer their excitation energy between a donor and an acceptor site, leading to the standard Hamiltonian The state |d (|a ) denotes the exciton to be at the donor (acceptor), J da is the respective electronic coupling between the donor and the acceptor and d ( a ) is the excited electronic energy of the donor (acceptor). The second line in equation (1) describes the coupling of the exciton to environmental Gaussian fluctuations generated by the harmonic bath H B j = 1 2 κ p 2 j,κ + ω 2 j,κ q 2 j,κ , with momenta p j,κ , displacement q j,κ , frequency ω j,κ and coupling ν ( j) κ . Additionally, equation (1) indicates that each site is coupled to its local bath, as can be 4 derived from molecular Hamiltonians of pigment-protein complexes [41,42]. In such a system, electronic de-excitation of a donor site and excitation of an acceptor site take place via nonequilibrium bath states in accordance with the Franck-Condon principle. The bath modes associated with each site then relax to the respective equilibrium states. This process is the socalled site-dependent reorganization dynamics. It should be noted that the standard Hamiltonian describing photosynthetic exciton transfer is slightly different from a simple spin-boson model, where a two-state system is coupled to a single harmonic bath; for a detailed discussion, see [21]. This spin-boson model is usually employed to describe electron transfer reactions in condensed phases, [43,44] not the exciton transfer.
In the following, we will assume separate environmental fluctuations for the donor and the acceptor but both with identical spectra. Several forms of the spectral density are employed in the literature, either based on model assumptions or based on analyses of molecular simulations. In this paper, we use the Debye spectrum with reorganization energy λ and timescale of environmental relaxation ω −1 c . The Debye spectral density has been successfully used for theoretical analyses of experimental results [14], [45][46][47]. In the field of open quantum dynamics, the environmental fluctuation spectrum typically is categorized by a cut-off frequency ω c and a dimensionless coupling strength α = λ/(hω c ). Often, exponentially cut-off spectra are used [36,48] and thus we will discuss in parallel which gives rise to the dimensionless coupling strength α = πλ/(2hω c ). Finally, we are interested in the quantum mechanical dynamics of the donor and acceptor influenced by the environmental fluctuations, but not in the actual environmental dynamics itself. The full dynamics is determined by the von Neumann equation for the statistical operator W (t) of the total system. Its solution permits the definition of a time evolution operator . Next, we average over the environmental degrees of freedom and need to find the effective time evolution for the reduced statistical operator of the donor-acceptor system, for all operators of the Hilbert space of the donor-acceptor system alone, where A = A S ⊗ 1 B . In passing, we note that more refined models have been studied in order to describe realistic exciton transfer, including models with spatially correlated environments [49], non-Ohmic spectral densities and structured environments with strongly localized vibrational modes [13,50].

Methods
The evaluation of the dynamics becomes particularly involved when none of the timescales can be assumed as small. In this, more refined methods that go beyond a Markov approximation are to be used. Here, we briefly summarize the two methods that we employ and that are at present rather well established. For further details, see the literature.

Quasi-adiabatic propagator path integral (QUAPI)
In order to calculate the effective dynamics of the donor-acceptor system, ρ(t), two of the present authors have employed the numerically exact QUAPI [36][37][38][39][40] scheme. It has been adopted to include multiple environments [21] specifically addressing site-dependent reorganization of protein environments, which play a significant role in the EET [10,29,42]. The algorithm, in brief, is based on a symmetric Trotter splitting of the short-time propagator R(t k+1 , t k ) for the full Hamiltonian into a part depending on the donor-acceptor Hamiltonian and a part involving the environment and the coupling term. The short-time propagator describes time evolution over a Trotter time slice δt. This splitting is by construction exact in the limit δt → 0, but introduces a finite Trotter error for a finite time increment, which has to be eliminated by choosing δt small enough such that convergence is achieved. On the other hand, the environmental degrees of freedom generate correlations that are nonlocal in time. For any finite temperature, these correlations decay exponentially fast at asymptotic times, thereby setting the associated memory timescale. QUAPI now defines an object called the reduced density tensor, which lives on this memory time window and establishes a temporal iteration scheme in order to determine the time evolution of this object. Within the memory time window, all correlations are included exactly over the finite memory time τ mem = K δt, but are neglected for times beyond τ mem . Then, the memory parameter K has to be increased, until convergence is found, meaning that all memory effects are taken into account up to a desired accuracy for the quantity of interest. Please note that the parameter K does not allow for an interpretation in terms of number for individual phonon creation processes since they are all summed over and thus included exactly within the memory time window. The phonon propagator is included via the Feynman-Vernon influence functional and already appears in the exponent, i.e. all phonon processes are captured at any instant of time within this window. What can be stated is that processes, which are due to correlations over larger times than the memory time, are not captured due to the cut-off in the memory time. Their weight, however, is exponentially suppressed and they are thus negligible.
Since the summation over all possible paths within the memory time window is exact and deterministic, the method does not suffer from any sign problem that sometimes renders quantum Monte Carlo methods troublesome. Likewise, its iterative nature allows in principle arbitrary long times to be reached. In practice, the QUAPI code runs on a standard Intel Xeon Processor (X5650 6C 2.66 GHz), typically calculating a typical result of a time-dependent occupation within a few hours (the resulting error bars from an extrapolation to vanishing Trotter time step are smaller than the size of the symbols in the figures shown below). The CPU time scales linearly with the simulation time. The scaling behavior of the approach with the system size is less advantageous. The size of the computational resources, i.e. the involved matrices, grows at most as M 2K +2 , where M is the number of basis states of the system Hilbert space. For M = 2 as in our case here, standard hardware architectures permit us to choose K 12-14. The required CPU time scales correspondingly, while the simulation time still grows only linearly. The exponential scaling property certainly limits the applicability of the scheme to not too large model systems.
In addition, the general idea of using a finite memory time of the bath-induced correlations is limited to finite-temperature environments. Only then, their exponential decay ensures the 6 existence of a rather short memory time and thus tractable computational objects. Conversely, this excludes strictly zero-temperature environments whose correlation functions decay only algebraically and induce long-time correlations in the system.

Second-order cumulant time-nonlocal quantum master equation (2CTNL)
For the model, expressed in equation (1), a reliable theoretical framework has been presented to describe EET in pigment-protein complexes specifically addressing site-dependent reorganization dynamics of protein environments [10]. The theory reduces to the Redfield [28] and Förster [31] theories in their respective limits of validity. This capability of interpolating between the two is crucial in order for the equation to be reliable in describing photosynthetic EET, which can occur between these two regimes.
Owing to the Gaussian fluctuations caused by the harmonic bath, the environmental effects can be characterized fully by two-point correlation functions of the collective bath coordinate u j ≡ κ ν ( j) κ q j,κ , e.g. the symmetrized correlation function, D j j (t), and the response function, j j (t). Note that the imaginary part of Fourier-Laplace transform of the response function is identical to the spectral density, G(ω). Then, the formally exact equation is expressed as [10] ∂ ∂tρ where a tilde indicates an operator in the interaction picture, and the integration kernel is given by The important point to note is that equation (5) is a time-nonlocal equation because the chronological time-ordering operator T + resequences and mixes the operators V j (t) × and V j (t) • comprised in K j (t, s) andρ(t). This time nonlocality is crucial to a correct description of non-Markovian interplay between electronic excitations and protein environments, i.e. site-dependent reorganization dynamics [10,42]. In [10], the classical fluctuation-dissipation relation, ∂ t D j j (t) −k B T j j (t), and the Debye spectral density were assumed, and the recurrence or hierarchy expansion technique [51][52][53] was employed for numerical calculations. For low-temperature systems where the classical fluctuation-dissipation relation breaks down, the quantum fluctuation-dissipation relation should be employed, D j j (t) = (h/π ) ∞ 0 dω G(ω) coth(βhω/2)cos(ωt), and thus the expansion needs to be modified by including quantum corrections [54]. Furthermore, it is possible to extend the treatment to cases of arbitrary spectral density with the help of the Meier-Tannor numerical decomposition scheme [55] which allows us to express the symmetrized correlation and response functions as the sums of complex exponential functions. For example, the Ohmic spectral density with an exponential cutoff, equation (3), can be very accurately parameterized by only three terms, N = 3, and thus the correlation and response functions can be approximated as the 7 sums of six complex exponential functions [55]. Specific values for the parameterization can be found in table I in [55]. Although the extended Meier-Tannor scheme was also discussed [56], we employ the original scheme in equation (7) for the sake of simplicity. In general, if the integration kernel in equation (5) can be expressed as K j (t, s) where A j and B j,α j are complex c-numbers and j,α j is an operator, then the following auxiliary operator can be introduced for the hierarchy expansion: j,α j (s)] n j,α j . This auxiliary operator yields the same form as equation (6) in [14]. In practice, calculations (performed by a Macbook Pro, 2.53 GHz Intel Core Duo) presented in the following section take from merely sub-seconds (figure 2) to about 2 h (figure 4). More details and specifically the scaling behavior with increasing system size for multichromophoric systems can be found in [14], which treats a larger system, the FMO complex comprising seven pigments.

Weak electronic coupling
We discuss first the regime of weak electronic donor-acceptor coupling. We calculate the population difference P(t) between the donor and the acceptor. Parameters are fixed as typical for photosynthetic exciton transfer, with a − d = 100 cm −1 , J da = 20 cm −1 , ω c = 53 cm −1 and T = 300 K for a Debye spectrum. The dynamics turns out to be overdamped (not explicitly shown), except for the smallest reorganization energies. Thus, we can safely fit the occupation difference P(t) using with the relaxation rate γ = k a←d + k d←a with the forward (backward) transfer rate k a←d (k d←a ) and the equilibrium occupation difference P eq = (k d←a − k a←d )/γ according to detailed balance. In figure 1, we plot the intersite transfer rate k a←d versus the reorganization energy λ. The plot shows the results from 2CTNL in comparison with a full Redfield and a Förster treatment as reported earlier [10]. The plot now includes results obtained using the numerically exact QUAPI scheme and shows excellent agreement with the 2CTNL results. This highlights the applicability of both approaches, QUAPI and 2CTNL, for reorganization energies of the order of the electronic couplings, λ J da , where neither weak-coupling approaches nor the orthodox Förster treatment can describe the dynamics satisfactorily. This range of reorganization energies is typically encountered in photosynthetic energy transfer.

Strong electronic coupling
In the case of strong electronic coupling the dynamics shows coherent oscillations, and it is no longer possible to characterize it using a single exponential and a single transition rate as before. We therefore plot in figure 2 the occupation ρ dd (t) of the donor versus time for the parameters a − d = 100 cm −1 and ω c = 53 cm −1 as before, but for a stronger electronic coupling J da = 100 cm −1 . As above, we have used a Debye spectrum. The temperature was chosen to be room temperature, T = 300 K. Furthermore, the coupling strength was chosen as α = 0.24 (or, respectively, the reorganization energy λ = J da /5) in figure 2(a) and α = 1.2 (λ = J da ) in figure 2(b). Again, almost perfect agreement between the two results was found. The off-diagonal element (coherence) ρ da (t) of the statistical operator also shows very good agreement between the two methods. In figure 2, black squares (magenta crosses) present the real (imaginary) part of ρ da (t) obtained by QUAPI. The full red lines are the corresponding 2CTNL results. Figure 3 shows the corresponding plots for a temperature T = 77 K at which the first experiments for finding coherent exciton dynamics have been performed [3,4]. In all cases we find coherent oscillations for several hundreds of femtoseconds where increasing damping and temperature suppress coherence more strongly, as expected. The red lines show the 2CTNL results and the blue circles, black squares and magenta crosses the QUAPI results for the donor population and the real and imaginary parts of the coherence. The results of both approaches coincide with each other, which allows us to conclude that both approaches describe this non-Markovian coherent dynamical regime correctly.

Exponential cut-off spectrum
There are not many details of the high-frequency tail of the environmental fluctuation spectrum known experimentally. Often, Debye spectra are used to analyze experiments [14], [45][46][47].
In the theoretical literature on open quantum dynamics, exponential cut-off spectra are often used [36,48]. The different cut-off functions hardly result in any measurable difference for the quantum dynamics of the donor-acceptor model for fast environmental fluctuations with ω c J da . However, this no longer holds true for slow environments ω c J da [30], as they occur in the exciton transfer dynamics. Naturally, the details of the cut-off function change quantitatively the decoherence of the exciton transport, as the cut-off frequency is a relevant timescale.
For a description of the dynamics with our methods, QUAPI and 2CTNL, the cut-off function enters the calculation of the environmental correlations. Whereas in QUAPI these are calculated only once and thus a modified cut-off function does not modify the computational efforts, in 2CTNL integrals involving the environmental spectra are calculated within every time step. Nevertheless, exponential cut-off spectra can be integrated into the formalism and we show the corresponding results in figure 4. Here, the occupation ρ dd (t) of the donor versus time is shown for an exponential cut-off environmental fluctuation spectrum at a temperature T = 300 K with ω c = 53 cm −1 using the same coupling strength as before. Note that this results in slightly different reorganization energies. The donor-acceptor parameters are set to a − d = 100 cm −1 and J da = 100 cm −1 . In figure 4(a), the coupling strength is α = 0.24, resulting in λ = (2/π )J da /5 and α = 1.2 (λ = (2/π)J da ) in figure 4(b). The results agree qualitatively with the previous results from a Debye spectrum, but, on a quantitative level, the coherent oscillations persist longer. The red lines show the 2CTNL results and the blue circles, black squares and magenta crosses the QUAPI results for the donor population and the real and imaginary parts of the coherence. Again, the results from both approaches coincide almost perfectly.

Conclusions and discussion
In conclusion, we have studied the real-time transfer dynamics of a biomolecular donoracceptor system as the basic building block for a description of photosynthetic exciton energy transport. For this purpose, we have compared two numerical methods, QUAPI and 2CTNL. For the relevant parameter range, namely reorganization energy of the order of the electronic coupling, site-dependent reorganization dynamics and long environmental relaxation time, i.e. of the order of the electronic transfer times ω −1 c = J −1 da , both methods yield perfectly coinciding results (within the numerical accuracy) and thus allow an accurate description of the dynamics of the energy transfer. As the competing timescales of system and fluctuations are comparable, it is a highly nontrivial task to accurately include the ensuing non-Markovian correlations of the fluctuations and it is significant to have reliable dynamical schemes to hand.
Studying the exciton dynamics in the FMO complex (a seven-chromophoric system) when disturbed by fluctuations following a Debye spectrum, QUAPI and 2CTNL also yield coinciding results [14,57]. Both methods, however, are numerically expensive, and although the accuracy of both methods is independent of the system size once convergent results are achieved, with current computer powers they are not (yet) capable of describing larger systems, such as excitonic transfer in Photosystems I and II of cyanobacteria, algae and green plants [1,2]. Clearly, numerical approaches can hardly cover the entire parameter space. What we can state is that in all the considered cases, QUAPI delivers converged results, and then the two approaches yield coinciding results. Our methods could serve as benchmarks for testing simpler methods in the physically relevant parameter ranges in order to solve physiologically relevant models, i.e. large systems with many pigments.
Finally, we briefly mention other recent methods introduced to describe non-Markovian environmental features. An efficient simulation method for the nonperturbative regime was recently developed by Prior et al [50], based on the combination of the time-dependent density matrix renormalization group with the theory of orthogonal polynomials. A related approach has been formulated in terms of hierarchically coupled spectral densities [58]. As formulated up to the present, the former treats the limiting case of zero temperature and becomes increasingly costly for larger temperatures. In turn, it is general enough to treat an arbitrary form of the spectral density. An alternative scheme has been put forward by Vacchini and Breuer [59], which formulates exact master equations for the decay of a two-state system into a structured reservoir. Therefore, a perturbative expansion of the generator of the dynamical equation is performed. The approach treats a localized environmental mode but becomes increasingly complicated when the non-Markovian corrections become more important.
A full understanding of physiologically relevant models and the importance of the experimentally observed long-lived electronic quantum coherence opens up the possibility of understanding the design principle of photosynthetic complexes and of exploiting the nearunity efficiency of energy transfer for technological devices. It may also allow one to clarify whether the high quantum efficiency is a result of the constructive interplay between quantum coherence and slow, spatially correlated environmental fluctuations. This could open a path to efficient future artificial light-harvesting complexes finding applications, for instance, in optimized organic solar cells.