Non-Markovian Quantum Jumps in Excitonic Energy Transfer

We utilize the novel non-Markovian quantum jump (NMQJ) approach to stochastically simulate exciton dynamics derived from a time-convolutionless master equation. For relevant parameters and time scales, the time-dependent, oscillatory decoherence rates can have negative regions, a signature of non-Markovian behavior and of the revival of coherences. This can lead to non-Markovian population beatings for a dimer system at room temperature. We show that strong exciton-phonon coupling to low frequency modes can considerably modify transport properties. We observe increased exciton transport, which can be seen as an extension of recent environment-assisted quantum transport (ENAQT) concepts to the non-Markovian regime. Within the NMQJ method, the Fenna-Matthew-Olson protein is investigated as a prototype for larger photosynthetic complexes.


I. INTRODUCTION
The initial step in photosynthesis is the excitonic transport of the energy captured from photons to a reaction center [1]. In this process, highly efficient transport occurs between interacting chlorophyll molecules embedded in a solvent and/or a protein environment [2]. The exciton transfer dynamics has been studied utilizing Förster theory in the limit of weak intermolecular coupling [3] or Redfield master equations in the limit of weak exciton-phonon coupling [4]. The latter approach describes the transport as dissipative dynamics for the reduced excitonic density matrix. Master equations are developed starting from projector operator techniques that separate relevant (system) from less relevant (phonon) degrees of freedom. Formally exact dynamics for the system is described by the Nakajima-Zwanzig equation with time-convolution [5,6] and the time-convolutionless (TCL) equation [7,8,9,10,11]. The first is equivalent to the chronological ordering prescription (COP) while the second corresponds to a partial ordering prescription (POP) of the time-ordering in a system-bath cumulant expansion [7,12,13]. The convolution kernel can be transformed into the TCL kernel by including the appropriate backward propagation for the density matrix [10,14]. The time-dependent Redfield equation is derived from a secondorder approximation in the system-bath interaction Hamiltonian [14]. Further imposing the Markov approximation leads to the standard time-independent Redfield equation. The dynamics of the populations of the density matrix and that of the coherences is separated when the secular approximation is employed, in which the master equation can be cast into Lindblad form. Recently, Palmieri et al. developed a prescription for reintroducing the coupling of populations and coherences with suitably defined Lindblad operators [15].
Non-Markovian (NM) effects can be taken into account by the time-convolutionless approach. Other frequently used methods explicitly include strongly coupled modes or environmental two-level systems into the system dynamics [16,17]. The Hilbert space size and thus numerical effort increases concomitantly. Kubo, Tanimura et al. developed a hierarchical treatment where auxiliary systems describe higher order system-bath interactions [18,19,20]. Xu et al. introduced an elegant filtering method for this approach [21,22]. While the hierarchical treatment is formally exact for Gaussian fluctuations, the infinite set of equations is truncated for numerical propagation. The method was recently applied to investigate long-lived coherence in the Fenna-Matthews-Olson (FMO) protein complex [23,24].
A Markovian master equation in Lindblad form can be simulated by means of the Monte-Carlo wavefunction method (MCWF) [25]. This numerical technique relies on the property that density matrix evolution is equivalent to an averaging of wavefunction trajectories, each of which is interrupted by stochastic, discontinuous quantum jumps. In this work, we employ the non-Markovian quantum jump approach (NMQJ), recently developed by Piilo et al. [26,27]. This method is a generalization of the MCWF to the case of non-Markovian dynamics derived from a TCL approach. The TCL approach can lead to time-dependent, oscillatory decoherence rates that have negative regions, a signature of NM behavior. These negative rates are shown to lead to a reversal of decoherence by defining appropriate quantum jumps and jump probabilities.
Compared to the explicit numerical integration of the master equation, the NMQJ approach has interesting features in the context of exciton transfer in chromophoric networks. First, the NMQJ approach is, similar to the MCWF, based on the propagation of wavefunctions and thus scales considerably better with the system size than approaches in Liouville space. It is therefore especially suitable for simulating large chromophoric networks of photosynthetic antenna systems. Second, positivity violation [10] can be efficiently detected during the simulation by a simple criterion for the negative jump probabilities [28]. Third, the trajectory picture allows for new insights into exciton dynamics. Quantum jumps related to negative transition rates can restore coherence and thus can provide an additional theoretical perspective on long-lived quantum coherence found in photosynthetic systems such as the Fenna-Matthews-Olson complex [23] and the reaction center of purple bacteria [29]. Here, we apply the NMQJ approach to dimer systems and the FMO complex. We observe population beatings arising from recurrence effects of the NM environment. We also find that in the non-Markovian regime transport can be enhanced compared to purely Markovian dynamics and thereby provide an extension to the recent ENAQT concept [30,31,32]. These effects are pronounced in situations when the main phonon-mode frequencies are much smaller than (i.e. "off-resonant" to) a particular system transition frequency.
In Sec. 2, we develop a time-convolutionless master equation for the dynamics of a single excitation, leading to timedependent rates. In Sec. 3, we discuss the spectral density and time-dependent rates in more detail and explain the physical situation where NM effects are considerable. In Sec. 4, we introduce the NMQJ method in the context of excitonic energy transfer. Finally, in Secs. 5-7, we analyze dimer systems and the Fenna-Matthews-Olson complex.

II. MASTER EQUATION
The transport dynamics of a single excitation is described by a master equation for the density matrix that includes coherent evolution, relaxation, and dephasing. In this work, we are mainly interested in the effect of slow phonon fluctuations and the memory of the bath on the excitonic energy transport dynamics. We utilize a time-convolutionless master equation to second order in the exciton-phonon coupling. We employ the secular approximation and focus on non-Markovian decoherence rates. The removal of the secular approximation requires modifications to the NMQJ approach and is left for future work. The validity and limitation of the Redfield approach with respect to parameters such as environmental coupling and temperature and the with repect to the neglect of fundamental processes such as the molecular reorganization has been discussed in detail in [4,33,34]. The complete Hamiltonian for an interacting N -chromophoric system in the single exciton manifold and including the phonon part is given by H = H S + H SB + H B . The system part is in tight-binding form: with the system part A m = |m m| and the bath part Each site is separately interacting with a set of phonon modes indicated by the subindex m of the bath part. The dimensionless coefficients λ i describe the coupling strength to each phonon mode. The phonon Hamil- , where the sum is over all phonon modes (at each site) described by the bosonic operators b i and frequencies ω i . For this work, we assume that the chromophores are coupled independently to their respective baths. Recent studies include spatial correlations into the exciton dynamics [35,36,37]. The second order TCL master equation for the reduced system density matrix in the interaction picture is given by [14]: (3) The characteristic double commutator arises from a secondorder perturbation treatment of the system-bath Hamiltonian. Note that it is assumed that the effect of the system on the bath is small such that the total system approximately factorizes for all times. Multiphonon processes, arising from higher order commutators, are not taken into account. An additional approximation in standard Redfield theory is that the phonon bath is always in equilibrium, which neglects molecular reorganization effects. Next, we introduce the operators A m (ω) = M − N = ω c * m (M )c m (N )|M N |, which describe the effect of the system-bath Hamiltonian in the eigenbasis of the system Hamiltonian, i.e. A m = ω A m (ω), where the sum is over all transitions in the single exciton manifold [14]. This leads to the master equation of the form: with the symmetrized correlation function, and the associated response function, The quantities in Eq. (5) and Eq. (6) are related to the real and imaginary parts of the bath correlator, (4) avoids the Markov approximation in the sense that the upper integration limit goes to t instead of ∞. One also observes the usual oscillating terms that mix population and coherences. Next, we perform the secular approximation, essentially averaging over these fast oscillating terms. Here, we would like to study the effect of the Markovian versus NM decoherence rates and formulate the master equation such that the NMQJ method can be straightforwardly applied. The secular approximation is justified in the slow decoherence regime when |ω − ω | −1 τ D for all transition frequency differences, where τ D is a general time scale of decoherence. Finally, we assume that every chromophore is embedded in an identical phonon environment, thus the m subscript for the correlator and the response function can be dropped [30]. One arrives at the master equation in the interaction picture: with the time-dependent Lamb shift, and the time-dependent rates, A transformation into the Schrödinger picture can be readily performed, resulting in the usual system Hamiltonian commutator term of the master equation [14]: leads to a Lamb-type renormalization of the energy levels. In the present work, we do not consider this term since we do not expect a qualitatively new contribution to the exciton dynamics [34,36].

III. SPECTRAL DENSITY AND TIME-DEPENDENT RATES
The main result for a Redfield master equation without Markov approximation is the time dependence of the rates. The decoherence rates depend on the phonon coupling strengths λ i . The spectral density (units of frequency) describes the coupling strength at a particular frequency: Assuming a continuous distribution of modes the spectral density can be modeled with various functional forms. In molecular energy transfer often an Ohmic spectral density with exponential or Drude-Lorentz cutoff is employed [24,30,38,39]. In Ref. [13] non-Markovian phonon sidebands in fluorescence spectra of the B777 complex were reproduced with a super-Ohmic spectral density. In this paper, we assume an Ohmic spectral density with exponential cutoff: The relevant quantities are the cutoff ω c and the reorganization energy λ. The cutoff determines the position of the peak of the spectral density and the reorganization energy is given by λ = dω J(ω) ω . In Fig. 1 (upper panel) we show the spectral density for a particular choice of parameters. We choose λ = 30cm −1 , which is typical for chromophores in photosynthetic systems [20], and ω c = 30cm −1 which corresponds to relatively slow phonon modes. Note that typical transition frequencies in the single exciton manifold such as ω ≈ 200cm −1 are located at the tail of the spectral density. The resulting Markovian relaxation rates are small. The strongly coupled, "off-resonant" modes at around 30cm −1 can lead to considerable non-Markovian effects of the decoherence rates.
For any spectral density and for the bosonic bath, the timedependent decoherence rate is derived from Eq. (9): Here, n(ω) is the bosonic distribution function. In the Markovian case (t → ∞), the spectral density is sampled only at the frequency ω, seen from the limiting behavior of the terms 1 ω±ω sin((ω ±ω)t). For the dephasing rate one obtains from Eq. (13) in the limit ω → 0: In the Markovian limit, the dephasing rate becomes linearly proportional to the temperature and the derivative of the spectral density at zero frequency [14]. In the non-Markovian case, a greater part of the spectral density is taken into account. This can lead to rich behavior of both relaxation and dephasing rates. In Fig. 1 (lower panel), we show the rates that follow from the spectral density (12) and the above choice of parameters. The relaxation rates in the NM case oscillate around the Markovian rates, have positive and negative regions, and finally converge to the Markovian rate on a time scale of 1ps. The NM dephasing rate converges from below to the Markovian limit on a similar time scale. This transient regime has been discussed in terms of slippage in the initial conditions e.g. in [40].

IV. NON-MARKOVIAN QUANTUM JUMPS
In this work, we perform a stochastic unraveling of the master equation with the non-Markovian quantum jump approach established in Ref. [26,27]. The master equation (10) is precisely in the form required for the NMQJ method. We give a brief summary of this technique. For every time t, one can separate the set of jump generators A m (ω) into A + m (ω) and A − m (ω) depending on the overall sign of the corresponding rate. That is for all A + m (ω) the rate is γ + (ω, t) > 0, while for all A − m (ω) the rate is γ − (ω, t) < 0. In the presence of only positive channels, the original MCWF method  13) and (14). The upper panel shows the Ohmic spectral density with exponential cutoff for the parameters λ = 30cm −1 and ωc = 30cm −1 . In the physical situation studied in the present work, the main strongly coupled modes are "off-resonant" to a transition frequency (200cm −1 in this figure). This leads to rich behavior of the corresponding rates at relevant time scales of around 1ps. The relaxation rates oscillate, turn negative, and converge to their Markovian limit (lower panel, left); the inset shows a magnification at times just before 1ps. The dephasing rate rises from zero and converges to the Markovian limit (lower panel, right).
can be employed [25,26]. The particular structure of the jump generators A m (ω) (see previous sections) leads to a relatively straightforward description of the quantum mechanical ensemble. The density matrix at all time can be written as: Here, |ψ i (t) is the initial state with statistical weight N 0 (t)/N. The exciton states |M are as defined above and have a statistical weight N M (t)/N . Initially, N 0 (0) = N and at all times the numbers N 0 (t) and N M (t) conserve probability, i.e.N 0 (t)+ M N M (t) = N. The time evolution consists of propagation of |ψ i (t) and stochastic changes of the weights N 0 (t) and N M (t). In general, one defines the effective Hamiltonian: where the sum is over positive and negative channels. The NMQJ method now describes the time evolution of the ensemble ρ(t) as a wavefunction evolution of the ensemble states with H eff interrupted by probabilistic, discontinuous jumps corresponding to the jump operators of all channels. Consider now a particular ensemble member |ψ(t) at time t evolving for a small time step δt. As in the MCWF, the no-jump evolution is: The positive jumps occur with probability P + mω (t) = δt γ + (ω, t) ψ(t)|A + † m (ω)A + m (ω)|ψ(t) and an ensemble member jumps according to: The negative jumps occur from the source state |ψ(t) to target states |ψ (t) if the source |ψ(t) has the property that it can be reached by a jump with A − m (ω) from the target state: Note that a negative jump can "undo" positive jumps that occurred earlier. The negative jump probability depends on the target state |ψ (t) and is P − where N (t) is the number of ensemble members in the target state and N (t) are the number of ensemble members in the source state. A Monte Carlo unraveling according to this prescription is shown to be equivalent to master Eq. (10), see Ref. [26]. Non-Markovian quantum jumps can explicitly lead to restored quantum coherence with this jump description that correctly handles negative rates in the master equation.
We end this section with a note regarding positivity of the density matrix. The master equation (10) with time-dependent rates is not guaranteed to ensure positivity of the density matrix. However, the NMQJ method yields a simple criterion for detecting when positivity is about to be violated based on the negative jump probability [28]. The P − mω (t) is inversely proportional to the number of ensemble members in the source state N (t). The case when N (t) becomes zero and the rate is negative at the same time is precisely when the master equation would violate positivity. The interpretation of this is that the environment tries to "undo" an event that never happened. Thus, based on the singularity of the negative jump probability one can easily detect unphysical time evolution in the algorithm. All results in this work originate from physical time evolution.
For the simulations we take the excitonic Hamiltonian to be of a particular form: V 12 = 87cm −1 and 2 = 120cm −1 . This form is equal to the Hamiltonian of the site 1 and 2 subsystem in the Fenna-Matthews-Olson complex given in [24]. Nevertheless, the effects presented here hold for a large variety of dimer Hamiltonians. The initial state is localized at site 1, i.e. ρ(0) = |1 1|. In Fig. 2, we show the time dependence of the decoherence rates (left panels) and the time evolution of the population and coherence elements of the density matrix in the site basis, i.e. ρ mn (t) = m|ρ(t)|n . We compare different temperatures in the Markovian (middle panels) and the non-Markovian description (right panels). The spectral density parameters are reorganization energy λ = 50cm −1 and cutoff ω c = 50cm −1 . The rates (left panels in Fig. 2) are similar as discussed for Fig. 1. The relaxation rates are generally smaller than the pure dephasing rate. The Markovian limit is reached on a time scale on the order of 1ps. Higher temperature leads to larger oscillations of the decoherence rates. In the Markovian case (middle panels in Fig. 2), we observe that the population and coherence oscillations die out very fast, especially at high temperatures. This is explained by the linear dependence of the dephasing rate on temperature.
In the non-Markovian case, the dynamics is considerably different. At low temperatures the beatings are similar to the Markovian case; decoherence rates are generally rather small and differences in NM versus Markovian are not pronounced. At increasing temperature, the quantum mechanical beatings live slightly longer in the NM case due to the smaller dephasing rate at short times. Quantum beatings can be recognized by an oscillating imaginary part of the coherence element of the density matrix. At large temperatures, another type of beating arises, which is due to the oscillatory relaxation rates. It leads to beatings of the populations matrix elements and the real part of the coherence matrix element. These beatings can be interpreted as recurrence of the NM environment; energy is emitted from the system into the environment during the positive regions of the decoherence rates and reabsorbed in the same decoherence channel during the negative regions.

VI. TRANSPORT IN THE NON-MARKOVIAN REGIME
In this section, we focus on transport properties in the non-Markovian regime. The general behavior of the exciton transport as a function of different parameters, contributions of various physical processes, and robustness of a chromophoric network can be investigated by theoretical measures such as the energy transfer efficiency and the transfer time [36,41]. Trapping sites model the reaction centers where charge separation and energy storage occurs in the photosynthetic system, neglecting further chemical detail. In this paper, we utilize a simpler measure to elucidate energy transport. We define the integrated probability of a particular excitonic state up to a certain time τ , the only free parameter. An explicit introduction of trapping sites and additional free parameters is not required. Formally, the measure is given by: where τ is the total integration time and |M is a particular exciton state given by the problem at hand. For this measure we choose an exciton state |M and focus on relaxation dynamics in the exciton basis. For example, in the Fenna-Matthews-Olson complex the exciton with the lowest energy, localized at site 3 and 4, would be an appropriate choice. A site-basis definition is straightforward within the NMQJ method. We analyze the transport properties of master equation (10) and the NMQJ unraveling and compare to the Markovian regime. We choose the Hamiltonian parameters in Eq. (20) as V 12 = 50cm −1 and 2 = 2V 12 . We assume that the system is initially in the energetically higher eigenstate |E 2 and investigate relaxation to the lower lying eigenstate |E 1 . We quantify the transport by the integrated probability of Eq. (21) using |E 1 , i.e.P 1 . In Fig. 3, we show the dependence of the transport on the essential parameters of the spectral density (reorganization energy λ and cutoff ω c ) and the temperature. If not explored as variables, the default parameters are λ = 30cm −1 , ω c = 30cm −1 , and room temperature (T = 300K). We investigate two integration times, τ = 1ps and τ = 4ps. The results are more pronounced for the shorter time. Shorter time scales are more relevant for smaller photosynthetic complexes such as the Fenna-Matthews-Olson complex [24]. As a result, we observe that transport can be enhanced in the NM situation where the rates are gives by Eq. (13). This study can be seen as an extension of ENAQT concepts to the non-Markovian regime.
In the upper panels of Fig. 3, the dependence ofP 1 as a function of the reorganization energy λ is investigated. Relaxation and dephasing rates are scaled linearly with λ. In general, the probabilityP 1 increases as a function of the reorganization energy. When relaxation rates are larger, thermal equilibration of the exciton populations is faster. At λ = 30cm −1 and τ = 1ps the non-Markovian probability isP 1 = 0.44 while the Markovian probability isP 1 = 0.27, a substantial difference for this rather small system. The improvement can be rationalized by the fact that the initially large NM relaxation rates lead to fast equilibration. The negative regions of the rate partially "undo" the positive region but the integrated population of the target exciton is overall larger than in the Markovian case. For reorganization energies beyond ≈90cm −1 we observe positivity violating time evolution identified with the criterion discussed earlier [28].
The middle panels of Fig. 3 show the dependence ofP 1 on the temperature. In the present case of energy transfer from a high to low exciton state, temperature can have an assisting effect for short times and for both Markovian and NM treatment [30]. For example, see the graphs for τ = 1ps. Increased thermal population of the phonon modes can lead to increased stimulated emission of exciton energy into the phonon bath and thus transport towards the lower exciton state. This effect becomes weaker for longer times, see the graphs for τ = 4ps. Absorption of energy from the phonon bath comes into play which transports the excitation from the lower exciton state back to the higher one. The temperature is more significant in the NM regime since the Bose functions in the rate integral are sampled at all frequencies instead of only at ω 21 .
In the lower panels of Fig. 3, the averaged probabilityP 1 is shown as a function of the cutoff ω c of the spectral density. stronger coupling of modes that are resonant with the transition frequency leads to increased transport. The differences in terms of transport of both cases vanish when ω c ≈ 60cm −1 ; resonant modes dominate the relaxation dynamics. However, in the presence of only slow modes, the NM treatment shows substantially larger transfer probabilities. For τ = 1ps and ω c = 20cm −1 , we obtainP 1 = 0.31 in the NM case and P 1 = 0.06 in the Markovian case. This improvement is due to sampling of the a broader range of the spectral density in Eq. (13). The physical interpretation is that the NM treatment allows the system to temporarily access energy nonconserving phonons for quantum jumps that would be inaccessible otherwise.

VII. FENNA-MATTHEWS-OLSON COMPLEX
We also performed simulations for the Fenna-Matthews-Olson (FMO) complex. The FMO complex acts as an energy transfer wire in green sulphur bacteria Chlorobium tepidum [2]. It is subject of recent experimental [23] and theoretical studies [24,30,35,42]. We derived the TCL master equa-tion (10) for the seven-site FMO subunit and performed simulations with NMQJ method. We used the Hamiltonian of ref. [24] and the spectral density (12) with λ = 35cm −1 and ω c = 150cm −1 [30]. The decoherence rates for the 42 relaxation channels (absorption + emission) and the 7 dephasing channels are evaluated. The rates oscillate and some have negative regions. The initial states are chosen to be localized at site 1 or 6, the sites that are close to the chlorophyll antenna complex. As a result, we find that the dynamics is not substantially affected by the time-dependent rates, see Fig. 4. Quantum beatings are slightly longer-lived for both temperatures 77K and 300K and both initial states. This is because the NM dephasing rates converge from below to the Markovian limit similar to Fig. 1 (bottom right panel). The main relaxation rates stay positive and oscillate around their Markovian values. The spectral density is rather broad, covering all transition frequencies, cf. Fig. 2 of [35], such that the effects described in the previous sections turn out to be not dominant. The Markovian approximation alone in the presence of the other approximations such as Born and secular does not have a substantial effect. Recently, Ishizaki and Fleming utilized the hierarchical equation of motion approach for explaining long- lived coherences in the FMO complex [24]. Since Born and secular approximations are avoided for Gaussian fluctuations, this approach has a larger range of validity than the Redfield model, especially at large temperatures, and correctly incorporates molecular reorganization effects.

VIII. CONCLUSION
In conclusion, we have applied the non-Markovian quantum jump method to excitonic energy transfer. The NM decoherence rates resulting from a time-convolutionless treatment of the master equation are oscillatory and negative for parameter regimes and time scales that are relevant to the problem. In the present work, NM effects are large when a system is strongly coupled to "off-resonant" phonon modes of the environment. These slow modes can lead to population beatings at room temperature, which are a signature of bath recurrence effects. Additionally, our computations show that Markovian versus non-Markovian dynamics can thus crucially affect transport dynamics. Quantum transport can be enhanced over strictly Markovian dynamics due to a sampling of broader regions of the spectral density. We thus have provided a non-Markovian extension to recent environment-assisted quantum transport (ENAQT) concepts. For example, a system with a transition of around 140cm −1 shows considerable NM improvement of transport in the presence of strong modes at around 30cm −1 .
Recently, Jang et al. [11] developed a novel theory of coherent resonance energy transfer. A small polaron transformation is applied before the second-order time-convolutionless expansion that leads to time-dependent decoherence rates. Nonequilibrium reorganization effects are taken into account by the exciton-phonon dressed state description. This treatment can lead to an increased range of validity with respect to the exciton-phonon couplings and temperatures compared to the standard Redfield approach. In this context, the NMQJ method in its present form or with suitable extensions could prove especially powerful to efficiently simulate larger donoracceptor systems and to correctly incorporate negative decoherence rates in a quantum jump description. Furthermore, the NMQJ method can be used for the stochastic computation of spectroscopic signals [43,44].