An Efficient Quantum Jump Method for Coherent Energy Transfer Dynamics in Photosynthetic Systems under the Influence of Laser Fields

We present a non-Markovian quantum jump approach for simulating coherent energy transfer dynamics in molecular systems in the presence of laser fields. By combining a coherent modified Redfield theory (CMRT) and a non-Markovian quantum jump (NMQJ) method, this new approach inherits the broad-range validity from the CMRT and highly efficient propagation from the NMQJ. To implement NMQJ propagation of CMRT, we show that the CMRT master equation can be casted into a generalized Lindblad form. Moreover, we extend the NMQJ approach to treat time-dependent Hamiltonian, enabling the description of excitonic systems under coherent laser fields. As a benchmark of the validity of this new method, we show that the CMRT-NMQJ method accurately describes the energy transfer dynamics in a prototypical photosynthetic complex. Finally, we apply this new approach to simulate the quantum dynamics of a dimer system coherently excited to coupled single-excitation states under the influence of laser fields, which allows us to investigate the interplay between the photoexcitation process and ultrafast energy transfer dynamics in the system. We demonstrate that laser-field parameters significantly affect coherence dynamics of photoexcitations in excitonic systems, which indicates that the photoexcitation process must be explicitly considered in order to properly describe photon-induced dynamics in photosynthetic systems. This work should provide a valuable tool for efficient simulations of coherent control of energy flow in photosynthetic systems and artificial optoelectronic materials.


Introduction
During the past decade, much progress has been achieved in both experimental and theoretical explorations of photosynthetic excitation energy transfer (EET), e.g., EET pathways determined by pump-probe as well as two-dimensional electronic spectroscopy (2DES) [1,2,3], coherent EET dynamics revealed by 2DES [4,5,6], and theoretical studies to elucidate mechanisms of EET [7,8,9,10]. It is intriguing to consider the possibility of using laser pulse to coherently control energy flow in photosynthetic complexes. Notably, coherent light sources were adopted to control energy transfer pathways in LH2 [11]. The phases of the laser field can efficiently adjust the ratio of energy transfer between intra-and inter-molecule channels in the complex's donoracceptor system. In addition, theoretical control schemes were put forward to prepare specified initial states and probe the subsequent dynamics [12], and the optimal control theory has been adopted to optimize the effects of electromagnetic field's polarization and structural and energetic disorder on localizing excitation energy at a certain chromophore within a Fenna-Matthews-Olson complex of green bacteria [13].
However, although the breakthroughs in experimental techniques have provided much insight into the coherent dynamics in small photosynthetic systems [14,15], efficient theoretical methods that can describe quantum dynamics in realistic photosynthetic antenna are still lacking, because the typical size of an antenna (> 100 chromophores) and the warm and wet environment make the accurate description of coherent EET in these systems a formidable task. Methods for simulating coherent EET dynamics in photosynthetic systems have been put forward (see [9] for a comprehensive review). For example, Förster theory directly provides EET rates, however it neglects quantum coherence effects that have been shown to be critical in photosynthetic systems [7]. In contrast, Redfield theory describes coherent EET dynamics, but it is valid only when the system-bath interactions are sufficiently weak [16,17]. In order to bridge the gap between these two methods, recent advances in theoretical methods have focused on new approaches valid in the intermediate regime. For example, the small-polaron master equation approach [18,19,20,21] should provide a reasonable description of the EET dynamics in the intermediate regime, and the hierarchical equation of motion (HEOM) actually yields a numerically exact means to calculate EET dynamics [22,23,24].
Nevertheless, it is still difficult to apply these theoretical methods to a realistic photosynthetic antenna, which often exhibits > 100 chromophores and static disorder, due to the formidable computational resources required. Moreover, to understand the fundamentals of light-matter interactions in photosynthetic light harvesting and realize laser control of energy flow in photosynthetic systems, explicit treatment of photoexcitation processes would be necessary. Furthermore, it was suggested recently that the photoexcitation condition determined by the nature of the photon source could strongly affect the appearance of EET dynamics in photosynthetic complexes [25,26,27,28,29]. As a result, the interplay between the photoexcitation process and ultrafast EET dynamics may be highly nontrivial, and it would be important to explicitly consider system-field interactions in a theoretical description of light harvesting. Therefore, an accurate yet numerically efficient method for simulating quantum dynamics in photosynthetic systems in the presence of coherent laser fields is highly desirable.
In this work, we combine a coherent modified Redfield theory (CMRT) [30] and a non-Markovian quantum jump (NMQJ) method [31,32,33] to develop an efficient approach for coherent EET dynamics in photosynthetic systems under the influence of light-matter interactions. The modified-Redfield theory [34,35] yields reliable population dynamics of excitonic systems and has been successfully applied to describe spectra and population dynamics in many photosynthetic complexes. [36,37,3] As a generalization of the modified-Redfield theory, the CMRT approach provides a complete description of quantum dynamics of the system's density matrix. Furthermore, we adopt a quantum jump method to achieve efficient propagation of the CMRT equations of motion. Using quantum jumps or quantum trajectories to unravel a set of equations of motion for the density matrix into a stochastic Schrödinger equation has proven to be a powerful tool for the simulation of quantum dynamics in an open quantum system [38,39,40,41,31,32]. In particular, Piilo et al. have developed a non-Markovian quantum jump (NMQJ) method that can be implemented for simulating non-Markovian dynamics in multilevel systems [31,32]. To make use of the NMQJ method, which has been shown to be efficient for simulating quantum dynamics in EET [33,42], we recast the CMRT equation into a generalized Lindblad form [43]. In addition, by transforming to a rotating frame, we extend CMRT to describe a time-dependent Hamiltonian, which is necessary to describe systems in the presence of laser fields. To examine the validity of the combined CMRT-NMQJ approach, we simulate the EET dynamics in the FMO complex. Moreover, we investigate a simple dimer system under the influence of laser pulses. We aim to shed light on how the photoexcitation process affects the coherence of EET dynamics in light harvesting. This paper is organized as follows. In the next section, we briefly describe CMRT for EET in a dimer system with time-independent Hamiltonian. Then, we extend the CMRT to treat time-dependent case by including interactions with a laser pulse applied to the dimer. Furthermore, in order to make the CMRT suitable for time-propagation by the NMQJ method, we show that the CMRT quantum master equation can be rewritten in a generalized Lindblad form. In section 3, we detail the implementation of the NMQJ method to solve the CMRT master equation. To verify the applicability of the new method, in section 4, we apply it to simulate the EET dynamics in FMO and compare the results to those from numerically exact simulations. Furthermore, in section 5, we calculate the dynamics of site populations and entanglement in a dimer system after it has been excited by a laser pulse to investigate the interplay between laser excitation and ultrafast EET in the system. Finally, the main results of this work are summarized in section 6. Figure 1. Energy diagram of a dimer for a time-independent case. The site energy for the single-excitation state |n is E n and the energy of the ground state |G is E 0 . The electronic coupling J between the single-excitation states is also considered.

Coherent Modified Redfield Theory
For the sake of self-consistency, we briefly describe the CMRT method in this section. Note that the focus of this work is the combined CMRT-NMQJ approach. The details in the derivation and validity of the CMRT approach is outside of the scope of this paper and will be published elsewhere [30].

Time-Independent Hamiltonian
Before proceeding to the case with a time-dependent Hamiltonian, we briefly review the CMRT approach for a time-independent Hamiltonian. The CMRT approach is a generalization of the original modified Redfield theory [34,35] to treat coherent evolution of the full density matrix of an excitonic system. For simplicity, we investigate the quantum dynamics in a dimer system (figure 1) described by the following excitonic Hamiltonian: where |n is the product state of the excited state on the nth site and the ground state on the other site, i.e., |1 = |e 1 |g 2 , |2 = |g 1 |e 2 , E n is the corresponding site energy, J is the electronic coupling between these two sites. Here, when there is no pulse applied to the system, the ground state |G = |g 1 |g 2 with energy E 0 is decoupled to the single-excitation states |n . Note that although the Hamiltonian for a dimer system is considered in this work, the results can be straightforwardly generalized to treat multi-site systems, e.g. the FMO with seven sites as shown in section 4. Furthermore, the environment is modeled as a collection of harmonic oscillators described by the Hamiltonian where p nq and x nq are respectively the momentum and position for qth harmonic oscillator on the nth site with mass m q and frequency ω q . Without loss of generality, we have assumed an independent bath for each site. Nevertheless, the CMRT can deal with a general system with correlated baths, whose effects on EET dynamics have been investigated in detail previously [44,21,45].
To describe exciton-bath interactions, we assume linear system-bath couplings described by where the collective bath-dependent coupling is where nq denotes the displacement of the qth harmonic mode when the nth site is excited.For the electronic part, we diagonalize the system Hamiltonian in the exciton basis [7] H (1) where each exciton basis state is a superposition of site-localized excitations with exciton energy ε ′ k . Since there is no interaction between the ground state and the single-excitation states, we have with In other words, Hereafter, for the sake of simplicity, we shall label |G as |0 .
Transforming H SB to the exciton basis, the system-bath interactions contain both diagonal and off-diagonal terms. Following the essence of the modified Redfield theory [34], we include the diagonal term of H SB into the zeroth-order Hamiltonian in the exciton basis to be treated nonperturbatively, and we only treat the off-diagonal term of H SB in the exciton basis as the perturbation. In the case, the unperturbed and perturbation Hamiltonians read respectively where a (1) Notice that a (1) kk ′ (n) = 0 as long as k = 0 or k ′ = 0 as a result of (9). We then apply the second-order cumulant expansion approach using V (1) as the perturbation followed by secular approximation to derive a quantum master equation for the reduced density matrix of the electronic system [34,46,30]: and the reorganization energy is Here ε (1) k can be viewed as the eigen value of the electronic Hamiltonian shifted by bath reorganization.The CMRT master equation (13) describes the coherent dynamics driven by the exciton Hamiltonian, the population transfer dynamics, and dephasing dynamics including pure-dephasing induced by diagonal fluctuations and populationtransfer induced dephasing. The dissipation rate R where g (1) is the line-broadening function evaluated from the spectral density of the system-bath couplings [46]: In addition, the pure-dephasing rate R (1)pd kk ′ (t) is given by The CMRT equation (13) can be considered as a generalization of the original modified Redfield theory to treat full coherent dynamics of the excitonic system. By treating the diagonal system-bath coupling in the exciton basis nonperturbatively, the CMRT considers multi-phonon effects and generally interpolates between the strong system-bath coupling and the weak system-bath coupling limits, in contrast to the Redfield theory [35,3]. Despite this, our method still shows its promising validity as compared to the numerically exact HEOM approach shown in section 4.
Functions required to propagate the CMRT master equation can be readily calculated if the system Hamiltonian H S and the spectral densities J n (ω) are provided. Generally speaking, the spectral densities used in theoretical investigations of photosynthetic systems are obtained by fitting to available experimental data, e.g., linear absorption lineshape or nonlinear spectra [7,47]. For self-consistency, the same theoretical method should be used in the fitting process of the parameters and dynamical simulations. Therefore, it is important that a dynamical method can also describe spectroscopy within the same framework. For example, the absorption lineshape defined as the Fourier transform of the dipole-dipole correlation function can be calculated, i.e., where the dipole-dipole correlation function is Within the framework of the CMRT under the same second-order approximation, the dipole-dipole correlation function reads [48,36] It corresponds to several peaks around the level spacings between the exciton states and the ground state, i.e., ε 0 . The second term in the exponential g kkkk (t) is the exchange-narrowed line-shape function [48], which is proportional to the participation ratio n C (1) kn 4 . Generally speaking, for a delocalized state n C (1) kn 4 is much smaller than that for a completely-localized state. Therefore, the widths of the peaks in the absorption line shape are significantly narrowed due to exciton delocalization. The third term in the exponential is the relaxation-induced broadening. It originates from the population relaxation within the single-excitation subspace when the molecule is excited by light, and has been demonstrated to be important in the spectra of molecular assemblies [36].

Coherent Modified Redfield Theory in the Presence of Laser Fields
To treat system-field interactions in the framework of the CMRT, we need to extend the CMRT equations of motion to include influences of laser fields. As shown in figure 2, during the time-evolution, we consider a dimer system interacting with a monochromatic pulse whose frequency is ω, thus inducing couplings between the ground state ε (1) 0 and delocalized exciton states ε (1) k (k = 1, 2). In this case, the electronic Hamiltonian reads where 2g k is the laser-induced coupling strength between the ground state ε . Transformed to a rotating frame with where we have applied the rotating-wave approximation and dropped the fast-oscillating terms with factors exp (±i2ωt). The above Hamiltonian can be diagonalized as where is the electronic eigen state with eigen energy ε ′′ k . In the basis of ε (2) k in the rotating frame, the unperturbed and perturbation Hamiltonians read respectively where a (2) Following the CMRT approach, we could obtain the master equation where ρ T (t) = U † ρ(t)U is the density matrix in the rotating frame, Here ε (2) k can be viewed as the eigen value of the electronic Hamiltonian In (31), the dissipation rate R (2)dis kk ′ (t) and pure-dephasing rate R (2)pd kk ′ (t) can be obtained following (17) and (21) by substituting ε Thus, the influence of the laser field can be considered as modifying the electronic part of the Hamiltonian by pulse induced couplings between ground state and singleexcitation states in the rotating frame. Notice that due to the transformation, the density matrix in the original frame ρ(t) is transformed to ρ T (t) = U † ρ(t)U. Hence, after solving the master equation, the obtained density matrix should be transformed back to the original frame as ρ(t) = Uρ T (t)U † . We further remark that (31) enables us to simulate field-induced and dissipative dynamics for an excitonic system due to a step-function pulse. However, this is not a limitation because of the time-local nature of the equations of motion. Within any propagation scheme based on discrete time steps, this formalism can be easily generalized to arbitrary pulse shapes as discussed in Appendix A.

Lindblad-form Coherent Modified Redfield Theory
We have obtained the CMRT master equation that governs the dynamics of an excitonic system in the presence of laser fields. For a system with M chromophores, the CMRT master equation corresponds to a set of M 2 ordinary differential equations. As a consequence, it would be a formidable task to propagate the dynamics if the system has hundreds of sites, just as in a typical photosynthetic complex, e.g., 96 chlorophylls in PSI. We can apply the NMQJ method [31,32] to reduce the computational complexity to the order of M to achieve efficient simulation of the CMRT dynamics. The NMQJ method propagates a set of equations of motion in the generalized Lindblad form [49]: Thus, in order to solve the master equation (cf. (13) and (31)) by means of the NMQJ method, we should rewrite them in the generalized Lindblad form. To this end, we define the following jump operators where |ε k (t) is kth eigen state of the system Hamiltonian H e (t), and the matrix elements of the population transfer and dephasing rates are defined respectively as where Γ k s are determined by the pure-dephasing rates in the CMRT master equation: where the matrix elements of A and B are respectively Note that since there are M(M − 1)/2 independent pure-dephasing rates R pd kk ′ in the CMRT master equation and only M Lindblad-form dephasing rates Γ k , the problem is overdetermined. We found that a least-square fit of Γ k to R pd kk ′ preserves the CMRT dynamics extremely well and is also numerically straightforward to implement. The details of the numerical determination of Lindblad-form dephasing rates are presented in Appendix B.

Non-Markovian Quantum Jump Method
We apply the non-Markovian quantum jump (NMQJ) approach proposed by Piilo et al. to simulate the CMRT equation of motion for multilevel systems [31,32]. The original implementation of the NMQJ approach only considers a time-independent Hamiltonian. In order to apply the NMQJ method to simulate CMRT dynamics under the influence of laser fields, here we briefly describe the NMQJ method and show how to explicitly implement the method for the time-dependent case to propagate the quantum dynamics of an M-level system under the influence of external time-dependent fields.
In our theoretical investigation, we consider the system to be initially in the ground state |ψ M (0) = |G and thus the initial density matrix is Then, the laser is turned on and applied to the system for a duration τ starting at t = 0. According to the NMQJ method, for an M-level system, there are M + 1 possible states for the system to propagate in, i.e., ε (2) k for k = 0, 1, · · · M − 1, in addition to the deterministic evolution where the non-Hermitian Hamiltonian is If N ensemble members are considered in the simulation, the NMQJ method describes the number of ensemble members in each of the individual states at time t, N (I) α (t) (α = 0, 1, · · · M), with the initial value given as During the pulse duration, the Lindblad-form master equation (31) is solved by the NMQJ method to obtain the distribution of N (I) α (t). Therefore, the density matrix in the rotating frame is straightforwardly given by At the end of the pulse, we transform the density matrix back to the original frame by where β k s are the expanding coefficients in the basis ε (1) k . Following the end of the pulse, the laser field is switched off and the system then evolves free of laser pulse influence, yet still follows the CMRT master equation describing population relaxation and dephasing in the exciton basis. In this case, since ε Thus, to describe both laser-induced dynamics and dissipative dynamics of an excitonic system with M chromophores, we must consider 2M + 1 possible states for the system to propagate in, i.e., ε Finally, the density matrix can be calculated from the time-dependent distribution of numbers of ensemble members in states N (II) We remark that the above procedure can be generalized to simulate arbitrary pulse shape as long as the laser pulses only induce transitions between the ground state and single-excitation states. In this case, the eigen states in the rotating frame are always the superpositions of the ground state and single-excitation states. Therefore, once the instantaneous eigen states are changed due to the change in system-field interaction, M new trajectories will be initialized as the previous eigen states are replaced by the eigen states of the new Hamiltonian. In Appendix A we present the numerical details for simulating the interaction between an excitonic system and a laser pulse with arbitrary pulse profile.
The generalization of the combined CMRT-NMQJ method for a system in the presence of laser fields can be easily extended to describe multiple pulses to enable the simulation of N-wave-mixing optical experiments. For example, the method described in this work can be applied to calculate third-order response functions [46]. Additionally, since the combined CMRT-NMQJ approach describes the full laser-driven dynamics of a dissipative excitonic system, it can be used together with the density matrix based approach for the calculation of photon-echo signals [50,51] to evaluate nonlinear spectra, including 2D electronic spectra and photon-echo peakshift signals. Although 2D electronic spectra can not be calculated without 2-exciton manifold, it would be straightforward to generalize the results presented in this work to include dynamics in the two-exciton manifold [34,50]. As a result, nonlinear spectra for coupled systems, which are often difficult to calculate due to the need to average over static disorder, can be calculated by using the combined CMRT-NMQJ approach. We emphasize that in addition to the favorable system-size scaling due to dynamical propagation in the wavefunction space, this approach benefits from the fact that the average over static disorder can be performed simultaneously with the average over quantum trajectories. Therefore, the combined CMRT-NMQJ method should provide a much more numerically efficient means for the evaluation of nonlinear spectra in large condensed-phase systems.
Note that in addition to the NMQJ method, quantum trajectory methods such as the non-Markovian quantum state diffusion approach [39,40] offer alternatives to the unraveling of quantum dynamics in light-harvesting systems [33,52]. In that approach, an integration of memory kernel and a functional derivative of the wavefunction with respect to the noise are required to obtain the stochastic Schrödinger equation. Recently, de Vega has applied the quantum state diffusion approach within the socalled post-Markov approximation to describe non-Markovian dissipative dynamics of a network mimicking a photosynthetic system [52]. However, generally speaking it is not straightforward to unravel a non-Markovian quantum master equation using the quantum state diffusion approach. In this work we aim to focus on the CMRT method, which is a time-local quantum master equation with a generalized Lindblad form. As a result the quantum jump approach proposed by Piilo and coworkers seems to be a natural way to implement the propagation of the time-local CMRT equation. A detailed comparison of the numerical efficiency and domains of applicability of the various quantum jump and quantum trajectory based methods is given in Refs. [41] and [32]. We believe these methods could provide superior numerical performance as well as useful physical insights for EET dynamics in complex and large systems.

EET Dynamics of FMO
In the previous sections, by combing the CMRT with the NMQJ method, we have proposed an efficient approach to simulate the EET dynamics of photosynthetic systems in the presence of laser fields. In order to implement the NMQJ method, the master equation obtained by the CMRT is revised in the Lindblad form. Since several approximations are made in order to efficiently simulate the quantum dynamics of a large quantum open system for a broad range of parameters, it is necessary to verify the validity of this new approach for photosynthetic systems.  [24]) is used in both calculations. Solid curves are obtained by the CMRT-NMQJ approach, and dashed curves depict results from the HEOM method. Red curves are for site 6, green curves for site 3, blue curves for site 4, and black curves for site 5. In order to make the figure clear, population dynamics for other sites whose populations always stay less than 0.1 are not shown here.
As a numerical demonstration, in figure 3, we compare the population dynamics of EET in FMO simulated by two different approaches, i.e., the combined CMRT-NMQJ approach and the numerically exact HEOM method from Ishizaki and Fleming [24]. Clearly, the results from the two methods are in good agreement. The combined CMRT-NMQJ method does not only reproduce the coherent oscillations in the shorttime regime but is also consistent with the HEOM in the long-time limit. In the shorttime regime the oscillations are somewhat suppressed in the case of CMRT-NMQJ approach. This is likely due to the simplification of the pure-dephasing process by the direct projection operations. If more realistic dephasing operators as implemented in [53] is applied to simulate the pure-dephasing process, we would expect a better accuracy from the CMRT-NMQJ approach. However, this additional effort may require much more computational time with only marginal improvements for EET dynamics in photosynthetic systems. Therefore, it is a reasonable tradeoff for numerical efficiency that we adopt the Lindblad-form CMRT and make use of the NMQJ method to simulate quantum dynamics of realistic photosynthetic systems Note that because the NMQJ is a numerical propagation scheme, the accuracy of the relaxation dynamics calculated by this method actually depends on the Lindbladform CMRT, which is effectively a generalization to the widely used modified-Redfield theory in terms of population dynamics. Modified-Redfield approach has been proven to provide excellent results for biased-energy systems in photosynthesis [3], therefore we expect the combined CMRT-NMQJ method described in this work to be broadly applicable to EET in photosynthetic systems.

Dimer Model
To demonstrate the combined CMRT-NMQJ simulation of laser-driven dynamics, we investigate the behavior of a model dimer system under the influence of laser excitation in this section. Note that in previous simulations of EET dynamics in photosynthetic systems the photoexcitation process was often overlooked and the instantaneous preparation of a given initial condition at zero time is normally assumed. Given that the selective preparation of a specific excitonic state is a highly non-trivial task (in particular the site-localized states often assumed in coherent EET simulations), it is crucial that an experimentally realistic photoexcitation process can be considered in a complete dynamical simulation [27,28]. Here we show that the combined CMRT-NMQJ approach is capable of describing the excitation process. We further investigate the behavior of a model dimer system (figure 2) under the influence of laser excitation to study the interplay between the laser excitation and ultrafast EET dynamics.

Absorption Spectrum
Before we examine the laser-induced dynamics we demonstrate that the absorption spectrum can be obtained by the CMRT method (cf. (24)). In figure 4, we show numerically calculated absorption lineshape for dimer system with ground state energy E 0 = −12800cm −1 , site energies E 1 = 120cm −1 and E 2 = 100cm −1 , electronic coupling J = 300cm −1 , and the electronic transition dipole moments in the eigen basis are set to be µ 10 = 10D and µ 20 = 5D, respectively. We further assume identical Ohmic baths described by the spectral density with the reorganization energy λ = 35cm −1 and the cut-off frequency ω c = 50cm −1 , and temperature T = 300K. The simulated spectrum exhibits two peaks located at ε 2 . Figure 4 shows that the absorption lineshape is well captured by the CMRT approach, which is important if a comparison to experimental data is to be made [36,47]. The relative peak heights are given by the ratio of |µ k0 | 2 . Here, the widths of the two peaks are nearly the same as the effect of relaxation-induced broadening is not significant in this particular parameter set. Physically speaking, the pure-dephasing process is much faster than the energy relaxation from the high-energy state ε (1) 2 to the low-energy state ε (1) 1 when the dimer system is excited. For other cases, the contribution from the relaxation-induced broadening might play a more important role [36].

Laser-induced Dynamics
In this section, we apply the CMRT-NMQJ approach to investigate the quantum dynamics of a dimer system after laser excitation. As shown in figure 5, the dimer system initially in the ground state is applied with a square laser pulse during the time period [0, t 1 ]. Note that a square pulse is used here for the sake of simplicity (see Appendix A). Due to the coupling induced by the laser pulse, the dimer is coherently excited into a superposition between the ground state and the single-excitation states. After the laser is turned off at t = t 1 , the dimer is left alone to evolve under the influence of system-bath couplings. The energy diagrams for the case with and without a laser pulse are depicted in figures 2 and 1, respectively. In our numerical simulation, we adopt the following parameters to model photosynthetic systems: the ground state energy E 0 = −12800cm −1 , the site energies E 1 = 200cm −1 and E 2 = 100cm −1 . Moreover, we consider two specific cases J = 120cm −1 and J = 20cm −1 to explore the electronic coupling's influence on the dynamics. We further assume identical Ohmic spectral densities (cf. (52)) for each site with a reorganization energy of λ = 35cm −1 , cut-off frequency ω c = 50cm −1 , and temperature T = 300K. In all cases, the laser duration t 1 is set to 100fs. In addition, in order to make the results converge within a reasonable computational time, we use a moderate time step dt = 1fs and average over a sufficiently-large number of trajectories, i.e., N = 10 5 . Figures 6(a) and 6(b) show the timeevolution of the site populations initialized by a laser pulse for strongly (J = 120cm −1 ) and weakly (J = 20cm −1 ) electronically coupled dimers, respectively. The laser carrier frequency is tuned at ω = 13000cm −1 , close to both transition energies. Therefore we expect the laser pulse to induce significant coherence between the two exciton states. As a result, during the pulse duration (from t = 0fs to t = 100fs), the populations on both sites oscillate with large amplitudes, for the laser field is nearly on resonance with the transitions between the ground state |ε  After the laser is turned off at t = 100fs, the population on the ground state remains invariant since there is no relaxation from the single-excitation states to the ground state in our model (cf. (3)), however, dynamics in the single-exciton manifold show strong J dependence. A comparison between figures 6(a) and 6(b) shows that although both systems evolve to reach a thermal equilibrium due to the system-bath couplings, the relaxation dynamics are qualitatively different due to the different electronic coupling strengths. In the case with strong electronic coupling (figure 6(a)) coherent oscillations persist for up to 400fs, and the system then quickly relaxes to the equilibrium. Clearly, the coherent energy relaxation is extremely efficient in this case and the photoexcitation dynamics and the coherent dynamics are intertwined [54,42]. We conclude that in strongly coupled systems the nature of the photoexcitation process must be explicitly considered in order to reasonably describe the photon-induced dynamics, as a fictitiously assumed initial condition prepared instantaneously at t = 0 could not correctly capture the photoexcitation process. In contrast, for the weak electronic coupling case in figure  6(b), the coherent oscillations stop immediately after the pulse is turned off, and the system relaxes incoherently and also slowly. In this case a clear time-scale separation exists and a specific consideration of the photoexcitation process may not be necessary.

Laser-induced Population dynamics
To investigate the effects of laser excitation condition, we explore the effect of laser detuning on the quantum dynamics in figure 7. In this case, the laser carrier frequency is tuned at ω = 13400cm −1 . Because the laser is largely detuned from the transition between |ε  can not be effectively excited when the pulse is present [55]. Indeed, figure 7(b) shows that in this weak electronic coupling case, site 1 is selectively excited within the pulse duration, since it is the higher-energy site in our dimer model. Consequently, figure 7(b) indicates that specific initial state, e.g., a specific site being excited, may be prepared if a proper pulse is applied for an appropriate duration. In contrast, figure 7(a) shows that the oscillation amplitudes of site populations, P 1 and P 2 , are nearly the same in the strong electronic coupling case, as in this case both sites contribute significantly to the delocalized exciton state |ε (1) 2 . Noticeably, the excitation-relaxation dynamics presented in figures 6(a) and 7(a) are markedly different, i.e. no coherent evolution in figure 7(a) after t = 100fs. Clearly, the dynamical behavior of the strongly coupled dimer depends critically on the excitation conditions. Our results further emphasize the notion that the photoexcitation process must be explicitly considered [27,28]. Note that we aim to demonstrate the validity of the combined CMRT-NMQJ method developed in this work. It is clear that much more efforts must be paid to elucidate the effects of light source in natural light-harvesting processes. In this regard, the methodology developed in this work appears to be valuable for the study of coherent EET dynamics and quantum coherent control in photosynthetic systems.

Entanglement dynamics
In the population dynamics, coherent evolutions are present when the laser excitation or the electronic coupling is strong. To more effectively follow coherent EET dynamics, it has been shown that entanglement is a excellent probe for coherent dynamics in photosynthetic systems [56,57,58,59,60,61,62]. Here, we investigate the time-evolution of entanglement in the dimer system as a demonstration of the combined CMRT-NMQJ method for coherent EET dynamics. Because each pigment in the dimer is a natural qubit, we especially focus on the evolution of concurrence [63]. For a two-qubits system with arbitrary density matrix ρ, the concurrence is defined as [63]    where λ j s are the square roots of eigen values of the matrix ρρ in the descending order. Here,ρ = (σ y ⊗ σ y ) ρ * (σ y ⊗ σ y ) with ρ * being the complex conjugate of the density matrix ρ. For a purely classical system, its concurrence vanishes, while it would be unity for the maximum entangled states, e.g., Bell states. For other states, the concurrence monotonously increases as the state is more quantum-mechanically entangled. An analytical expression that yields C(t) from the matrix elements of ρ(t) can be derived. In our model, since we restrict the evolution of the system within the sub-space without the double-excitation state, the density matrix is of the block type where the basis states are |0 = |g 1 |g 2 , |1 = |e 1 |g 2 , |2 = |g 1 |e 2 , and |3 = |e 1 |e 2 . Thus, the concurrence is simplified as Because ρ 11 +ρ 22 ≥ 2 |ρ 21 |, the more population is excited to the single-excitation states, the larger the concurrence could be. In addition, at some time, when ρ 21 vanishes, we would expect the system evolves into the disentangled state as the concurrence disappears. In figure 8, we plot concurrence dynamics for the condition that the excitation pulse is in near resonance with both transitions. For both strong and weak electronic coupling cases, there exists significant entanglement induced by the optical excitation. However, the concurrence dynamics after the laser is switched off at 100fs are markedly different. Figure 8(a) shows that for the strongly coupled dimer the concurrence quickly evolves to reach a steady state with a considerable value as a result of the large electronic coupling. In contrast, for the weak electronic coupling case shown in figure 8(b), the concurrence rapidly decays to near zero and then follows by a continuous but slow rise. These results are consistent with the characteristics of population dynamics as shown in figure 6, indicating that the concurrence provides an effective tool to describe coherent energy transfer process.
Moreover, the evolution of concurrence with only one exciton state selectively excited is shown in figure 9. In the duration of pulse, the maximum entanglement generated by the field is significantly smaller than that by a resonance field. This is because less coherence between the two transitions can be created as a consequence of the laser de-tuning. In addition, because the laser is only in close resonance with one of the transitions, as shown in figure 9(a), the concurrence at the steady state is less than that achieved in figure 8(a), as less population can be excited to the single-excitation states by the weaker pulse. This discovery again emphasizes the importance of the state preparation by photoexcitation.
Besides, we notice that for all cases there is an abrupt change at t = 100fs due to the discontinuity of the system-field interaction, which is an artifact of the stepfunction pulse shape. We remark that this problem can be solved by using a smooth pulse profile. In a realistic experiment, a gaussian pulse should be applied to the system. The combined CMRT-NMQJ approach can be generalized to consider a gaussian pulse by decomposing it into a series of rectangular sub-pulses with equal width (see details in Appendix A). Again, these results indicate that the photoexcitation process plays an important role in the induction of initial coherence in excitonic systems, therefore system-field couplings must be explicitly considered in order to elucidate electronic quantum coherence effects in photosynthesis [28,14].

Conclusions
In this paper, we have combined the CMRT [30] approach and the NMQJ method [31,32] to develop a formalism that provides efficient simulations of coherent EET dynamics of photosynthetic complexes under the influence of laser fields. In order to implement the NMQJ propagation of the CMRT master equation, the original CMRT master equation is revised in the Lindblad form. In addition, the NMQJ approach is generalized to be suitable for the case with a time-dependent Hamiltonian due to interactions with laser fields. These new developments allow the efficient calculation of quantum dynamics of photosynthetic complexes in the presence of laser fields.
To demonstrate the effectiveness and efficiency of this new approach, we apply the CMRT-NMQJ approach to simulate the coherent energy transfer dynamics in the FMO complex, and compare the result to the one obtained by the numerically exact HEOM method. We show that both results are consistent in the long-time and short-time regimes. Furthermore, we investigate photon-induced dynamics in a dimer system. For strongly coupled dimers, coherent oscillations are observed in the population dynamics during the pulse duration as well as the free evolution as predicted by theory. In addition to the quantum dynamics in the excited states, we also consider coherent effects induced by the laser detuning in the stage of state preparation. We show that the dynamical behavior of a strongly coupled dimer depends critically on the excitation conditions, which further emphasizes the point that the photoexcitation process must be explicitly considered in order to properly describe photon-induced dynamics in photosynthetic systems. In addition, we investigate the evolution of concurrence of a dimer in the presence of laser fields. These results demonstrate the combined CMRT-NMQJ approach is capable of simulation of EET in photosynthetic complexes under the influence of laser control. Note that although the phases of the fields are fixed in our calculations, the combined CMRT-NMQJ method is fully capable of describing quantum coherent control phenomena by shaping phases of the laser fields.
Compared to other numerical simulation methods, this new approach has several key advantages. First of all, because it is based on the CMRT, its validity in a broad range of parameters has been well characterized [3,30], which is confirmed in simulating the energy transfer dynamics of the FMO complex. Note that the CMRT might fail in the regime where the exciton basis is not a suitable choice [3], e.g., λ > J > ∆ε. In that case, a combined Redfield-Förster picture could be used to provide accurate simulations of EET dynamics [64,65]. Because the generalized Förster dynamics only add diagonal population transfer terms between block-delocalized exciton states, which directly fit into a Lindblad form, an extension of the methods describe here to simulate EET dynamics in the Redfield-Förster picture will be straightforward. Second, as this approach incorporates the NMQJ method, it effectively reduces the calculation time [31] and therefore is capable of simulating the quantum dynamics in large photosynthetic systems. For an M-level system in the absence of laser fields, the effective ensemble size N eff = M + 1 in the CMRT-NMQJ approach can be sufficiently smaller than that needed in the non-Markovian quantum state diffusion, e.g. N = 10 4 for a two-level system in Ref. [40]. In addition, the problem of positivity violation, which may be encountered by master equation approaches, can be resolved by turning off the quantum jump once the population of the source state vanishes. Third, it can make use of parallel computation, because the sampling of trajectories can be efficiently parallelized, and it can take the average over static disorder at the same time as the average over trajectories. Finally, since the absorption spectrum [36] and other nonlinear optical signals should be efficiently simulated within the same theoretical framework, all the parameters used in the simulation can be self-consistently obtained by fitting the experimental data using this theoretical approach.
Our motivation to extend the CMRT theory to include time-dependent fields and NMQJ propagation is to provide a tool for the study of quantum control of photoexcitation and energy flow. In order to simulate EET dynamics in a coherent control experiment, it is necessary to consider a time-dependent Hamiltonian that includes the influence of control laser fields. Quantum coherent control lies at the heart of human's exploration in the microscopic world [66,67,68]. Its amazing power manifests in manipulating quantum dynamics by systematic control design methods [69,70,71,72]. The strong motivation to apply quantum coherent control to EET lies in not only the aspiration to knowledge but also the energy need to learn from the highly-efficient photosynthetic systems to improve artificial light harvesting. We believe the methodologies presented in this work would be useful for further investigation of using coherent quantum control methods to direct energy flow in photosynthetic as well as artificial light-harvesting systems. For example, we are currently using the theoretical framework described in this work to investigate how laser pulse-shape and -phase affect EET dynamics in photosynthetic systems. where a = 1, 2, . . . M. After simplification, we obtain a system of linear equations for the Lindblad-form dephasing rates, that is where the coefficient on the right hand side is Or equivalently, (B.6) can be written in a matrix form as with the matrix elements of B given by as given in (37).