An efficient quantum jump method for coherent energy transfer dynamics in photosynthetic systems under the influence of laser fields

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 (2D) electronic spectroscopy (2DES) [1–3], coherent EET dynamics revealed by 2DES [4–6], and theoretical studies to elucidate mechanisms of EET [7–10]. It is intriguing to consider the possibility of using laser pulses 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 donor–acceptor 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 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–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–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 photo-excitation 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–29]. As a result, the interplay between the photoexcitation process and ultrafast EET dynamics may be highly non-trivial, 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–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 [3, 36, 37]. 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 [31, 32, 38–41]. In particular, Piilo et al have developed a 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 Fenna–Matthews–Olson (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 the 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.

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 the scope of this paper and will be published elsewhere [30].

2.1. 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:

$$\begin{eqnarray}&&H_{S}^{(1)}=\mathop{\mathop{\sum }}\limits_{n=1}^{2}{{E}_{n}}\left| n \right\rangle \left\langle n \right|+J\left( \left| 1 \right\rangle \left\langle 2 \right|+\left| 2 \right\rangle \left\langle 1 \right| \right)+{{E}_{0}}\left| G \right\rangle \left\langle G \right|,\end{eqnarray} \tag{ 1 }$$

where $|n\rangle$ is the product state of the excited state on the nth site and the ground state on the other site, i.e., $|1\rangle =|e{{\rangle }_{1}}|g{{\rangle }_{2}}$, $|2\rangle =|g{{\rangle }_{1}}|e{{\rangle }_{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\rangle =|g{{\rangle }_{1}}|g{{\rangle }_{2}}$ with energy ${{E}_{0}}$ is decoupled to the single-excitation states $|n\rangle$. 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.

Zoom In Zoom Out Reset image size

Furthermore, the environment is modeled as a collection of harmonic oscillators described by the Hamiltonian

$$\begin{eqnarray}&&{{H}_{B}}=\mathop{\mathop{\sum }}\limits_{n=1}^{2}\mathop{\mathop{\sum }}\limits_{q}\left[ \frac{p_{nq}^{2}}{2{{m}_{q}}}+\frac{1}{2}{{m}_{q}}\omega _{q}^{2}x_{nq}^{2} \right],\end{eqnarray} \tag{ 2 }$$

where ${{p}_{nq}}$ and ${{x}_{nq}}$ are respectively the momentum and position for the qth harmonic oscillator on the nth site with mass ${{m}_{q}}$ and frequency ${{\omega }_{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 [21, 44, 45].

To describe exciton–bath interactions, we assume linear system–bath couplings described by

$$\begin{eqnarray}&&{{H}_{SB}}=\mathop{\mathop{\sum }}\limits_{n=1}^{2}{{u}_{n}}\left| n \right\rangle \left\langle n \right|,\end{eqnarray} \tag{ 3 }$$

where the collective bath-dependent coupling is

$$\begin{eqnarray}&&{{u}_{n}}=\mathop{\mathop{\sum }}\limits_{q}{{m}_{q}}\omega _{q}^{2}x_{nq}^{\left( 0 \right)}{{x}_{nq}},\end{eqnarray} \tag{ 4 }$$

where $x_{nq}^{\left( 0 \right)}$ 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]

$$\begin{eqnarray}&&H_{S}^{(1)}=\mathop{\mathop{\sum }}\limits_{k=0}^{2}\epsilon _{k}^{\prime }\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{k}^{(1)} \right|,\end{eqnarray} \tag{ 5 }$$

where each exciton basis state is a superposition of site-localized excitations

$$\begin{eqnarray}&&\left| \epsilon _{k}^{(1)} \right\rangle =\mathop{\mathop{\sum }}\limits_{n=0}^{2}C_{kn}^{(1)}\left| n \right\rangle ,\end{eqnarray} \tag{ 6 }$$

with exciton energy $\epsilon _{k}^{\prime }$. Since there is no interaction between the ground state and the single-excitation states, we have

$$\begin{eqnarray}&&\left| \epsilon _{0}^{(1)} \right\rangle =\left| 0 \right\rangle =\left| G \right\rangle \end{eqnarray} \tag{ 7 }$$

with

$$\begin{eqnarray}&&\epsilon _{0}^{\prime }={{E}_{0}}.\end{eqnarray} \tag{ 8 }$$

In other words,

$$\begin{eqnarray}&&C_{0n}^{(1)}={{\delta }_{n0}}.\end{eqnarray} \tag{ 9 }$$

Hereafter, for the sake of simplicity, we shall label $|G\rangle$ as $|0\rangle$.

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 non-perturbatively, and we only treat the off-diagonal term of ${{H}_{SB}}$ in the exciton basis as the perturbation. In this case, the unperturbed and perturbation Hamiltonians read respectively

$$\begin{eqnarray}&&H_{0}^{(1)}=\mathop{\mathop{\sum }}\limits_{k=0}^{2}\left[ \epsilon _{k}^{\prime }+{{H}_{B}}+\mathop{\mathop{\sum }}\limits_{n=1}^{2}{{u}_{n}}a_{kk}^{(1)}\left( n \right) \right]\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{k}^{(1)} \right|,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{V}^{(1)}}=\mathop{\mathop{\sum }}\limits_{k=0}^{2}\mathop{\mathop{\sum }}\limits^{2}_{{{k}^{\prime }}=0\left( \ne {{k}^{\prime }} \right)}\mathop{\mathop{\sum }}\limits_{n=1}^{2}{{u}_{n}}a_{k{{k}^{\prime }}}^{(1)}\left( n \right)\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{(1)} \right|,\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&a_{k{{k}^{\prime }}}^{(1)}\left( n \right)=C_{kn}^{(1)*}C_{{{k}^{\prime }}n}^{(1)}.\end{eqnarray} \tag{ 12 }$$

Notice that $a_{k{{k}^{\prime }}}^{(1)}(n)=0$ as long as k = 0 or ${{k}^{\prime }}=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 [30, 34, 46]:

$$\begin{eqnarray} {{\partial }_{t}}\rho & = & -{\rm{i}}\mathop{\mathop{\sum }}\limits_{{\rm{k}},{{{\rm{k}}}^{\prime }}}\left( \epsilon _{k}^{(1)}-\epsilon _{{{k}^{\prime }}}^{(1)} \right){{\rho }_{k{{k}^{\prime }}}}\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{(1)} \right| \\ {} & {} & -\frac{1}{2}\mathop{\mathop{\sum }}\limits_{k,{{k}^{\prime }}}\left[ R_{{{k}^{\prime }}k}^{(1){\rm{dis}}}\left( t \right){{\rho }_{kk}}\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{k}^{(1)} \right|-R_{k{{k}^{\prime }}}^{(1){\rm{dis}}}\left( t \right){{\rho }_{{{k}^{\prime }}{{k}^{\prime }}}}\left| \epsilon _{{{k}^{\prime }}}^{(1)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{(1)} \right| \right] \\ {} & {} & -\mathop{\mathop{\sum }}\limits_{k,{{k}^{\prime }}}\left\{ R_{k{{k}^{\prime }}}^{(1){\rm{pd}}}\left( t \right)+\frac{1}{2}\mathop{\mathop{\sum }}\limits_{{{k}^{\prime\prime }}}\left[ R_{{{k}^{\prime\prime }}k}^{(1){\rm{dis}}}\left( t \right)+R_{{{k}^{\prime\prime }}{{k}^{\prime }}}^{(1){\rm{dis}}}\left( t \right) \right] \right\}{{\rho }_{k{{k}^{\prime }}}}\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{(1)} \right|, \\ \end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&\epsilon _{k}^{(1)}=\epsilon _{k}^{\prime }-\Lambda _{k}^{(1)},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\Lambda _{k}^{(1)}=\mathop{\mathop{\sum }}\limits_{n=1}^{2}{{\left[ a_{kk}^{(1)}\left( n \right) \right]}^{2}}{{\lambda }_{n}},\end{eqnarray} \tag{ 15 }$$

and the reorganization energy is

$$\begin{eqnarray}&&{{\lambda }_{n}}=\mathop{\mathop{\sum }}\limits_{q}\frac{1}{2}{{m}_{q}}\omega _{q}^{2}x_{nq}^{\left( 0 \right)2}.\end{eqnarray} \tag{ 16 }$$

Here $\epsilon _{k}^{(1)}$ can be viewed as the eigenvalue 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 population-transfer induced dephasing. The dissipation rate $R_{k{{k}^{\prime }}}^{(1){\rm{dis}}}(t)$ is defined as

$$\begin{eqnarray} R_{k{{k}^{\prime }}}^{(1){\rm{dis}}}\left( t \right) & = & 2\operatorname{Re}\mathop{\int }_{0}^{t}{\rm{d}}\tau {{{\rm{e}}}^{{\rm{i}}\left( \epsilon _{{{k}^{\prime }}}^{(1)}-\epsilon _{k}^{(1)} \right)\tau }}{{{\rm{e}}}^{-{\rm{i}}\left( {\rm{ \Lambda }}_{k}^{(1)}+{\rm{ \Lambda }}_{{{k}^{\prime }}}^{(1)}-2\lambda _{kk,{{k}^{\prime }}{{k}^{\prime }}}^{(1)} \right)\tau }}{{{\rm{e}}}^{-g_{kkkk}^{(1)}\left( \tau \right)-g_{{{k}^{\prime }}{{k}^{\prime }}{{k}^{\prime }}{{k}^{\prime }}}^{(1)}\left( \tau \right)+2g_{kk,{{k}^{\prime }}{{k}^{\prime }}}^{(1)}\left( \tau \right)}} \\ {} & {} & \times \left\{ \mathop{g}\limits^{\ddot{\ }{}}_{{{k}^{\prime }}k,k{{k}^{\prime }}}^{(1)}\left( \tau \right)-\left[ \dot{g}_{{{k}^{\prime }}k,kk}^{(1)}\left( \tau \right)-\dot{g}_{{{k}^{\prime }}k,{{k}^{\prime }}{{k}^{\prime }}}^{(1)}\left( \tau \right)-2{\rm{i}}\lambda _{{{k}^{\prime }}k,{{k}^{\prime }}{{k}^{\prime }}}^{(1)} \right] \right. \\ {} & {} & \left. \times \left[ \dot{g}_{k{{k}^{\prime }},kk}^{(1)}\left( \tau \right)-\dot{g}_{k{{k}^{\prime }},{{k}^{\prime }}{{k}^{\prime }}}^{(1)}\left( \tau \right)-2{\rm{i}}\lambda _{k{{k}^{\prime }},{{k}^{\prime }}{{k}^{\prime }}}^{(1)} \right] \right\}, \\ \end{eqnarray} \tag{ 17 }$$

where $g_{{{k}_{1}}{{k}_{2}}{{k}_{3}}{{k}_{4}}}^{(1)}(t)$ is the line-broadening function evaluated from the spectral density of the system–bath couplings [44]:

$$\begin{eqnarray}&&g_{{{k}_{1}}{{k}_{2}}{{k}_{3}}{{k}_{4}}}^{(1)}\left( t \right)=\mathop{\mathop{\sum }}\limits_{n}a_{{{k}_{1}}{{k}_{2}}}^{(1)}\left( n \right)a_{{{k}_{3}}{{k}_{4}}}^{(1)}\left( n \right){{g}_{n}}\left( t \right),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\lambda _{{{k}_{1}}{{k}_{2}}{{k}_{3}}{{k}_{4}}}^{(1)}=\mathop{\mathop{\sum }}\limits_{n}a_{{{k}_{1}}{{k}_{2}}}^{(1)}\left( n \right)a_{{{k}_{3}}{{k}_{4}}}^{(1)}\left( n \right){{\lambda }_{n}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{g}_{n}}\left( t \right)=\mathop{\int }{\rm{d}}\omega \frac{{{J}_{n}}\left( \omega \right)}{{{\omega }^{2}}}\left[ \left( 1-{\rm cos} \omega t \right){\rm coth} \left( \frac{\beta \omega }{2} \right)+{\rm{i}}\left( {\rm sin} \omega t-\omega t \right) \right].\end{eqnarray} \tag{ 20 }$$

In addition, the pure-dephasing rate $R_{k{{k}^{\prime }}}^{(1){\rm{pd}}}(t)$ is given by

$$\begin{eqnarray}&&R_{k{{k}^{\prime }}}^{(1){\rm{pd}}}\left( t \right)=\mathop{\mathop{\sum }}\limits_{n=1}^{2}{{\left[ a_{kk}^{(1)}\left( n \right)-a_{{{k}^{\prime }}{{k}^{\prime }}}^{(1)}\left( n \right) \right]}^{2}}\operatorname{Re}\left[ {{{\dot{g}}}_{n}}\left( t \right) \right].\end{eqnarray} \tag{ 21 }$$

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 non-perturbatively, 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 [3, 33]. Despite this, our method still shows 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}}(\omega )$ 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.

$$\begin{eqnarray}&&I\left( \omega \right)\equiv \operatorname{Re}\left[ \mathop{\int }_{0}^{\infty }{\rm{d}}t\chi \left( t \right){{{\rm{e}}}^{{\rm{i}}\omega t}} \right],\end{eqnarray} \tag{ 22 }$$

where the dipole–dipole correlation function is

$$\begin{eqnarray}&&\chi \left( t \right)\equiv {\rm{Tr}}\left[ \mu \left( t \right)\mu \left( 0 \right) \right].\end{eqnarray} \tag{ 23 }$$

Within the framework of the CMRT under the same second-order approximation, the dipole–dipole correlation function reads [36, 48]

$$\begin{eqnarray}&&\chi \simeq \mathop{\mathop{\sum }}\limits_{k}{{\left| {{\mu }_{kG}} \right|}^{2}}{\rm exp} \left[ {\rm{i}}\left( \epsilon _{0}^{(1)}-\epsilon _{k}^{(1)} \right)t-g_{kkkk}^{(1)}\left( t \right)-\frac{1}{2}\mathop{\mathop{\sum }}\limits_{{{k}^{\prime }}\ne k}R_{{{k}^{\prime }}k}^{(1){\rm{dis}}}\left( \infty \right)t \right].\end{eqnarray} \tag{ 24 }$$

It corresponds to several peaks around the level spacings between the exciton states and the ground state, i.e. $\epsilon _{k}^{(1)}-\epsilon _{0}^{(1)}$. The second term in the exponential $g_{kkkk}^{(1)}(t)$ is the exchange-narrowed line-shape function [48], which is proportional to the participation ratio ${{\mathop{\sum }}_{n}}{{\left| C_{kn}^{(1)} \right|}^{4}}$. Generally speaking, for a delocalized state ${{\mathop{\sum }}_{n}}{{\left| C_{kn}^{(1)} \right|}^{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].

2.2. CMRT 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 $|\epsilon _{0}^{(1)}\rangle$ and delocalized exciton states $|\epsilon _{k}^{(1)}\rangle$ (k = 1, 2). In this case, the electronic Hamiltonian reads

$$\begin{eqnarray}&&H_{S}^{\left( 2 \right)}=\mathop{\mathop{\sum }}\limits_{k=0}^{2}\epsilon _{k}^{\prime }\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{k}^{(1)} \right|+2{\rm cos} \omega t\mathop{\mathop{\sum }}\limits_{k=1}^{2}{{g}_{k}}\left( \left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{0}^{(1)} \right|+{\rm{h}}.{\rm{c}}. \right),\end{eqnarray} \tag{ 25 }$$

where $2{{g}_{k}}$ is the laser-induced coupling strength between the ground state $|\epsilon _{0}^{(1)}\rangle$ and the kth delocalized exciton state $|\epsilon _{k}^{(1)}\rangle$. Transformed to a rotating frame with

$$\begin{eqnarray*}&&U={\rm exp} \left[ {\rm{i}}\omega t\left| \epsilon _{0}^{(1)} \right\rangle \left\langle \epsilon _{0}^{(1)} \right| \right],\end{eqnarray*}$$

the effective Hamiltonian $H_{{\rm{eff}}}^{\left( 2 \right)}={{U}^{\dagger }}H_{S}^{\left( 2 \right)}U+{\rm{i}}{{\dot{U}}^{\dagger }}U$ of the electronic part becomes

$$\begin{eqnarray}&&{\boldsymbol{H}} _{{\rm{eff}}}^{\left( 2 \right)}\simeq \mathop{\mathop{\sum }}\limits_{k=1}^{2}\epsilon _{k}^{\prime }\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{k}^{(1)} \right|+\left( \epsilon _{0}^{\prime }+\omega \right)\left| \epsilon _{0}^{(1)} \right\rangle \left\langle \epsilon _{0}^{(1)} \right|+\mathop{\mathop{\sum }}\limits_{k=1}^{2}{{g}_{k}}\left( \left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{0}^{(1)} \right|+{\rm{h}}.{\rm{c}}. \right),\end{eqnarray} \tag{ 26 }$$

where we have applied the rotating-wave approximation and dropped the fast-oscillating terms with factors ${\rm exp} (\pm \;{\rm{i}}2\omega t)$. The above Hamiltonian can be diagonalized as

$$\begin{eqnarray}&&H_{{\rm{eff}}}^{\left( 2 \right)}=\mathop{\mathop{\sum }}\limits_{k=0}^{2}\epsilon _{k}^{\prime\prime }\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{k}^{\left( 2 \right)} \right|,\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle =\mathop{\mathop{\sum }}\limits_{{{k}^{\prime }}^{2}=0}C_{k{{k}^{\prime }}}^{\left( 2 \right)}\left| \epsilon _{{{k}^{\prime }}}^{(1)} \right\rangle =\mathop{\mathop{\sum }}\limits_{n=0}^{2}\mathop{\mathop{\sum }}\limits_{{{k}^{\prime }}^{2}=0}C_{k{{k}^{\prime }}}^{\left( 2 \right)}C_{{{k}^{\prime }}n}^{(1)}\left| n \right\rangle \end{eqnarray} \tag{ 28 }$$

is the electronic eigenstate with eigen energy $\epsilon _{k}^{\prime\prime }$.

Zoom In Zoom Out Reset image size

In the basis of $\{|\epsilon _{k}^{\left( 2 \right)}\rangle \}$ in the rotating frame, the unperturbed and perturbation Hamiltonians read respectively

$$\begin{eqnarray} &&H_{0}^{\left( 2 \right)} = \mathop{\sum }\limits_{k=0}^{2}\left[ \epsilon _{k}^{\prime\prime }+{{H}_{B}}+\mathop{\sum }\limits_{n=1}^{2}{{u}_{n}}a_{kk}^{\left( 2 \right)}\left( n \right) \right]\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{k}^{\left( 2 \right)} \right|, \\ &&{{V}^{\left( 2 \right)}} = \mathop{\sum }\limits_{k=0}^{2}\mathop{\sum }\limits_{{{k}}^{\prime }=0\left( \ne k \right)}^{2}\mathop{\sum }\limits_{n=1}^{2}{{u}_{n}}a_{k{{k}^{\prime }}}^{\left( 2 \right)}\left( n \right)\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{\left( 2 \right)} \right|, \\ \end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&a_{k{{k}^{\prime }}}^{\left( 2 \right)}\left( n \right)=\mathop{\mathop{\sum }}\limits_{{{k}_{1}}=0}^{2}C_{k{{k}_{1}}}^{\left( 2 \right)*}C_{{{k}_{1}}n}^{(1)*}\mathop{\mathop{\sum }}\limits_{{{k}_{2}}=0}^{2}C_{{{k}^{\prime }}{{k}_{2}}}^{\left( 2 \right)}C_{{{k}_{2}}n}^{(1)}.\end{eqnarray} \tag{ 30 }$$

Following the CMRT approach, we could obtain the master equation

$$\begin{eqnarray} {{\partial }_{t}}\;{{\rho }^{{\rm{T}}}} & = & -{\rm{i}}\mathop{\mathop{\sum }}\limits_{k,{{k}^{\prime }}}\left( \epsilon _{k}^{\left( 2 \right)}-\epsilon _{{{k}^{\prime }}}^{\left( 2 \right)} \right)\rho _{k{{k}^{\prime }}}^{{\rm{T}}}\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{\left( 2 \right)} \right| \\ {} & {} & -\frac{1}{2}\mathop{\mathop{\sum }}\limits_{k,{{k}^{\prime }}}\left[ R_{{{k}^{\prime }}k}^{\left( 2 \right){\rm{dis}}}\left( t \right)\rho _{kk}^{{\rm{T}}}\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{k}^{\left( 2 \right)} \right|-R_{k{{k}^{\prime }}}^{\left( 2 \right){\rm{dis}}}\left( t \right)\rho _{{{k}^{\prime }}{{k}^{\prime }}}^{{\rm{T}}}\left| \epsilon _{{{k}^{\prime }}}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{\left( 2 \right)} \right| \right] \\ {} & {} & -\mathop{\mathop{\sum }}\limits_{k,{{k}^{\prime }}}\left\{ R_{k{{k}^{\prime }}}^{\left( 2 \right){\rm{pd}}}\left( t \right)+\frac{1}{2}\mathop{\mathop{\sum }}\limits_{{{k}^{\prime\prime }}}\left[ R_{{{k}^{\prime\prime }}k}^{\left( 2 \right]{\rm{dis}}}\left( t \right)+R_{{{k}^{\prime\prime }}{{k}^{\prime }}}^{\left( 2 \right){\rm{dis}}}\left( t \right) \right) \right\}\rho _{k{{k}^{\prime }}}^{{\rm{T}}}\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{\left( 2 \right)} \right|, \\ \end{eqnarray} \tag{ 31 }$$

where ${{\rho }^{{\rm{T}}}}(t)={{U}^{\dagger }}\rho (t)U$ is the density matrix in the rotating frame,

$$\begin{eqnarray}&&\epsilon _{k}^{\left( 2 \right)}=\epsilon _{k}^{\prime\prime }-\Lambda _{k}^{\left( 2 \right)}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\Lambda _{k}^{\left( 2 \right)}=\mathop{\mathop{\sum }}\limits_{n=1}^{2}{{\left[ a_{kk}^{\left( 2 \right)}\left( n \right) \right]}^{2}}{{\lambda }_{n}}.\end{eqnarray} \tag{ 33 }$$

Here $\epsilon _{k}^{\left( 2 \right)}$ can be viewed as the eigenvalue of the electronic Hamiltonian

$$\begin{eqnarray}&&H_{e}^{\left( 2 \right)}=\mathop{\mathop{\sum }}\limits_{k}\epsilon _{k}^{\left( 2 \right)}\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{{{k}^{\prime }}}^{\left( 2 \right)} \right|.\end{eqnarray} \tag{ 34 }$$

In (31), the dissipation rate $R_{k{{k}^{\prime }}}^{\left( 2 \right){\rm{dis}}}(t)$ and pure-dephasing rate $R_{k{{k}^{\prime }}}^{\left( 2 \right){\rm{pd}}}(t)$ can be obtained following (17) and (21) by substituting $\epsilon _{k}^{(1)}$, $\Lambda _{k}^{(1)}$, $g_{{{k}_{1}}{{k}_{2}}{{k}_{3}}{{k}_{4}}}^{(1)}\left( t \right)$, $\lambda _{{{k}_{1}}{{k}_{2}}{{k}_{3}}{{k}_{4}}}^{(1)}$ respectively with $\epsilon _{k}^{\left( 2 \right)}$, $\Lambda _{k}^{\left( 2 \right)}$,

$$\begin{eqnarray}&&g_{{{k}_{1}}{{k}_{2}}{{k}_{3}}{{k}_{4}}}^{\left( 2 \right)}\left( t \right)=\mathop{\mathop{\sum }}\limits_{n}a_{{{k}_{1}}{{k}_{2}}}^{\left( 2 \right)}\left( n \right)a_{{{k}_{3}}{{k}_{4}}}^{\left( 2 \right)}\left( n \right){{g}_{n}}\left( t \right),\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&\lambda _{{{k}_{1}}{{k}_{2}}{{k}_{3}}{{k}_{4}}}^{\left( 2 \right)}=\mathop{\mathop{\sum }}\limits_{n}a_{{{k}_{1}}{{k}_{2}}}^{\left( 2 \right)}\left( n \right)a_{{{k}_{3}}{{k}_{4}}}^{\left( 2 \right)}\left( n \right){{\lambda }_{n}}.\end{eqnarray} \tag{ 36 }$$

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 single-excitation states in the rotating frame. Notice that due to the transformation, the density matrix in the original frame $\rho (t)$ is transformed to ${{\rho }^{{\rm{T}}}}(t)={{U}^{\dagger }}\rho (t)U$. Hence, after solving the master equation, the obtained density matrix should be transformed back to the original frame as $\rho (t)=U{{\rho }^{{\rm{T}}}}(t){{U}^{\dagger }}$. 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.

2.3. Lindblad-form CMRT

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]:

$$\begin{eqnarray}&&{{\partial }_{t}}\rho =-{\rm{i}}\left[ {{H}_{e}}\left( t \right),\rho \right]-\frac{1}{2}\mathop{\mathop{\sum }}\limits_{k,{{k}^{\prime }}}{{R}_{k{{k}^{\prime }}}}\left( t \right)\left[ \left\{ A_{k{{k}^{\prime }}}^{\dagger }\left( t \right){{A}_{k{{k}^{\prime }}}}\left( t \right),\rho \right\}-2{{A}_{k{{k}^{\prime }}}}\left( t \right)\rho A_{k{{k}^{\prime }}}^{\dagger }\left( t \right) \right].\end{eqnarray} \tag{ 37 }$$

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

$$\begin{eqnarray}&&{{A}_{k{{k}^{\prime }}}}\left( t \right)=\left| {{\epsilon }_{k}}\left( t \right) \right\rangle \left\langle {{\epsilon }_{{{k}^{\prime }}}}\left( t \right) \right|,\end{eqnarray} \tag{ 38 }$$

where $|{{\epsilon }_{k}}(t)\rangle$ is kth eigenstate of the system Hamiltonian ${{H}_{e}}(t)$, and the matrix elements of the population transfer and dephasing rates are defined respectively as

$$\begin{eqnarray}{{R}_{k{{k}^{\prime }}}}\equiv \left\{ \begin{array}{ccccccccccccccc} R_{k{{k}^{\prime }}}^{{\rm{dis}}}, & {\rm{for}}\;k\ne {{k}^{\prime }}, \\ {{\Gamma }_{k}}, & {\rm{for}}\;k={{k}^{\prime }}, \\ \end{array} \right.\end{eqnarray} \tag{ 39 }$$

where ${{\Gamma }_{k}}$s are determined by the pure-dephasing rates in the CMRT master equation:

$$\begin{eqnarray}&&\Gamma ={{B}^{-1}}A,\end{eqnarray} \tag{ 40 }$$

where the matrix elements of A and B are respectively

$$\begin{eqnarray}&&{{A}_{a}}=\mathop{\mathop{\sum }}\limits_{k=a+1}^{M}R_{ak}^{{\rm{pd}}}+\mathop{\mathop{\sum }}\limits_{k=1}^{M-1}R_{ka}^{{\rm{pd}}},\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}{{B}_{jk}}=\left\{ \begin{array}{ccccccccccccccc} \frac{1}{2}, & {\rm{for}}\;k<j, \\ \frac{1}{2}(2M-j), & {\rm{for}}\;k=j, \\ 1, & {\rm{for}}\;j<k<M, \\ \frac{1}{2}, & {\rm{otherwise}}. \\ \end{array} \right.\end{eqnarray} \tag{ 42 }$$

Note that since there are $M(M-1)/$2 independent pure-dephasing rates $R_{k{{k}^{\prime }}}^{{\rm{pd}}}$ in the CMRT master equation and only M Lindblad-form dephasing rates ${{\Gamma }_{k}}$, the problem is overdetermined. We found that a least-square fit of ${{\Gamma }_{k}}$ to $R_{k{{k}^{\prime }}}^{{\rm{pd}}}$ 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.

We apply the 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 $|{{\psi }_{M}}(0)\rangle =|G\rangle$ and thus the initial density matrix is

$$\begin{eqnarray}&&\rho \left( 0 \right)=\left| {{\psi }_{M}}\left( 0 \right) \right\rangle \left\langle {{\psi }_{M}}\left( 0 \right) \right|.\end{eqnarray} \tag{ 43 }$$

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. $|\epsilon _{k}^{\left( 2 \right)}\rangle$ for $k=0,1,\cdots M-1$, in addition to the deterministic evolution

$$\begin{eqnarray}&&\left| \psi _{M}^{{\rm{T}}}\left( t+\delta t \right) \right\rangle ={{{\rm{e}}}^{-{\rm{i}}H\left( t \right)\delta t}}\left| \psi _{M}^{{\rm{T}}}\left( t \right) \right\rangle \|{{{\rm{e}}}^{-{\rm{i}}H\left( t \right)\delta t}}{{\left| \psi _{M}^{{\rm{T}}}\left( t \right) \right\rangle }\| ^{-1}},\end{eqnarray} \tag{ 44 }$$

where the non-Hermitian Hamiltonian is

$$\begin{eqnarray}&&H\left( t \right)={{H}_{e}}\left( t \right)-\frac{{\rm{i}}}{2}\mathop{\mathop{\sum }}\limits_{k,{{k}^{\prime }}}{{R}_{k{{k}^{\prime }}}}\left( t \right)A_{k{{k}^{\prime }}}^{\dagger }\left( t \right){{A}_{k{{k}^{\prime }}}}\left( t \right).\end{eqnarray} \tag{ 45 }$$

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_{\alpha }^{({\rm{I}})}(t)$ ($\alpha =0,1,\cdots M$), with the initial value given as

$$\begin{eqnarray}N_{\alpha }^{\left( {\rm{I}} \right)}\left( 0 \right)=\left\{ \begin{array}{ccccccccccccccc} 0, & {\rm{for}}\ \alpha =0,1,\cdots M-1, \\ N, & {\rm{for}}\ \alpha =M. \\ \end{array} \right.\end{eqnarray} \tag{ 46 }$$

During the pulse duration, the Lindblad-form master equation (31) is solved by the NMQJ method to obtain the distribution of $N_{\alpha }^{\left( {\rm{I}} \right)}(t)$. Therefore, the density matrix in the rotating frame is straightforwardly given by

$$\begin{eqnarray}&&{{\rho }^{{\rm{T}}}}\left( t \right)=\frac{1}{N}\left[ \mathop{\mathop{\sum }}\limits_{k=0}^{M-1}N_{k}^{\left( {\rm{I}} \right)}\left( t \right)\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle \left\langle \epsilon _{k}^{\left( 2 \right)} \right|+N_{M}^{\left( {\rm{I}} \right)}\left( t \right)\left| \psi _{M}^{{\rm{T}}}\left( t \right) \right\rangle \left\langle \psi _{M}^{{\rm{T}}}\left( t \right) \right| \right].\end{eqnarray} \tag{ 47 }$$

At the end of the pulse, we transform the density matrix back to the original frame by $\rho (\tau )=U(\tau ){{\rho }^{{\rm{T}}}}(\tau ){{U}^{\dagger }}(\tau )$ to yield

$$\begin{eqnarray}&&\left| {{\psi }_{M}}\left( \tau \right) \right\rangle =U\left( \tau \right)\left| \psi _{M}^{{\rm{T}}}\left( \tau \right) \right\rangle ={{\beta }_{0}}{{e}^{{\rm{i}}\omega \tau }}\left| \epsilon _{0}^{(1)} \right\rangle +\mathop{\mathop{\sum }}\limits_{k=1}^{M-1}{{\beta }_{k}}\left| \epsilon _{k}^{(1)} \right\rangle ,\end{eqnarray} \tag{ 48 }$$

where ${{\beta }_{k}}{\rm{s}}$ are the expanding coefficients in the basis $\{|\epsilon _{k}^{(1)}\rangle \}$.

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 $|\epsilon _{k}^{(2)}\rangle$ are not the eigenstates of $H_{e}^{(1)}$, we must expand the space of jump states and initialize M new states, i.e. for $k=0,1,\cdots M-1$,

$$\begin{eqnarray}&&\left| {{\psi }_{k+M+1}}\left( \tau \right) \right\rangle =U\left( \tau \right)\left| \epsilon _{k}^{\left( 2 \right)} \right\rangle =C_{k0}^{\left( 2 \right)}{{{\rm{e}}}^{{\rm{i}}\omega \tau }}\left| \epsilon _{0}^{(1)} \right\rangle +\mathop{\mathop{\sum }}\limits_{{{k}^{\prime }}=1}^{M-1} C_{k{{k}^{\prime }}}^{\left( 2 \right)}\left| \epsilon _{{{k}^{\prime }}}^{(1)} \right\rangle ,\end{eqnarray} \tag{ 49 }$$

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. $|\epsilon _{k}^{(1)}\rangle$ for $k=0,1,\cdots M-1$, and $|{{\psi }_{j}}(t)\rangle$ for $j=M,M+1,\cdots 2M$. The corresponding numbers of ensemble members in the states $N_{\alpha }^{\left( {\rm{II}} \right)}$ ($\alpha =0,1,\ldots 2M$) are determined by the NMQJ method to solve the CMRT master equation (13) with the initial values given as

$$\begin{eqnarray}N_{\alpha }^{\left( {\rm{II}} \right)}\left( \tau \right)=\left\{ \begin{array}{ccccccccccccccc} 0, & {\rm{for}}\ \alpha =0,1,\cdots M-1, \\ N_{M}^{\left( {\rm{I}} \right)}\left( \tau \right), & {\rm{for}}\ \alpha =M, \\ N_{\alpha -M-1}^{\left( {\rm{I}} \right)}\left( \tau \right), & {\rm{for}}\ \alpha =M+1,M+2,\cdots 2\;M. \\ \end{array} \right.\end{eqnarray} \tag{ 50 }$$

Finally, the density matrix can be calculated from the time-dependent distribution of numbers of ensemble members in states $N_{\alpha }^{\left( {\rm{II}} \right)}\left( t \right)$ as

$$\begin{eqnarray}&&\rho \left( t \right)=\frac{1}{N}\left[ \mathop{\mathop{\sum }}\limits_{k=0}^{M-1}N_{k}^{\left( {\rm{II}} \right)}\left( t \right)\left| \epsilon _{k}^{(1)} \right\rangle \left\langle \epsilon _{k}^{(1)} \right|+\mathop{\mathop{\sum }}\limits_{j=M}^{2\,M}N_{j}^{\left( {\rm{II}} \right)}\left( t \right)\left| {{\psi }_{j}}\left( t \right) \right\rangle \left\langle {{\psi }_{j}}\left( t \right) \right| \right].\end{eqnarray} \tag{ 51 }$$

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 eigenstates in the rotating frame are always the superpositions of the ground state and single-excitation states. Therefore, once the instantaneous eigenstates are changed due to the change in system–field interaction, M new trajectories will be initialized as the previous eigenstates are replaced by the eigenstates 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 cannot be calculated without two-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 so-called 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 [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.

In the previous sections, by combining 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.

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 short-time regime but is also consistent with the HEOM in the long-time limit. In the short-time regime the oscillations are somewhat suppressed in the case of the 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] are applied to simulate the pure-dephasing process, we would expect 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.

Zoom In Zoom Out Reset image size

Note that because the NMQJ is a numerical propagation scheme, the accuracy of the relaxation dynamics calculated by this method actually depends on the Lindblad-form CMRT, which is effectively a generalization to the widely used modified-Redfield theory in terms of population dynamics. The 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.

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.

5.1. Absorption spectrum

Before we examine the laser-induced dynamics, we demonstrate that the absorption spectrum can be obtained by the CMRT method (see (24)). In figure 4, we show numerically calculated absorption lineshape for dimer system with ground state energy ${{E}_{0}}=-12800\;{\rm{c}}{{{\rm{m}}}^{-1}}$, site energies ${{E}_{1}}=120\;{\rm{c}}{{{\rm{m}}}^{-1}}$ and ${{E}_{2}}=100\;{\rm{c}}{{{\rm{m}}}^{-1}}$, electronic coupling J = 300 cm$^{-1}$, and the electronic transition dipole moments in the eigen basis are set to be ${{\mu }_{10}}=10\;{\rm{D}}$ and ${{\mu }_{20}}=5\;{\rm{D}}$, respectively. We further assume identical Ohmic baths described by the spectral density

$$\begin{eqnarray}&&J\left( \omega \right)=\lambda \frac{\omega }{{{\omega }_{c}}}{{{\rm{e}}}^{-\omega /{{\omega }_{c}}}},\end{eqnarray} \tag{ 52 }$$

with the reorganization energy $\lambda =35\;{\rm{c}}{{{\rm{m}}}^{-1}}$ and the cut-off frequency ${{\omega }_{c}}=50\;{\rm{c}}{{{\rm{m}}}^{-1}}$ and temperature T = 300 K. The simulated spectrum exhibits two peaks located at $\epsilon _{2}^{(1)}-\epsilon _{0}^{(1)}$ and $\epsilon _{1}^{(1)}-\epsilon _{0}^{(1)}$, corresponding to the transitions from the ground state $|\epsilon _{0}^{(1)}\rangle$ to two exciton states $|\epsilon _{1}^{(1)}\rangle$ and $|\epsilon _{2}^{(1)}\rangle$.

Zoom In Zoom Out Reset image size

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 $|{{\mu }_{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 $|\epsilon _{2}^{(1)}\rangle$ to the low-energy state $|\epsilon _{1}^{(1)}\rangle$ when the dimer system is excited. For other cases, the contribution from the relaxation-induced broadening might play a more important role [36].

5.2. 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 $\left[ 0,{{t}_{1}} \right]$. 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.

Zoom In Zoom Out Reset image size

In our numerical simulation, we adopt the following parameters to model photosynthetic systems: the ground state energy ${{E}_{0}}=-12800\;{\rm{c}}{{{\rm{m}}}^{-1}}$, the site energies ${{E}_{1}}=200\;{\rm{c}}{{{\rm{m}}}^{-1}}$ and ${{E}_{2}}=100\;{\rm{c}}{{{\rm{m}}}^{-1}}$. Moreover, we consider two specific cases J = 120 cm$^{-1}$and J = 20 cm$^{-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 $\lambda =35\;{\rm{c}}{{{\rm{m}}}^{-1}}$, cut-off frequency ${{\omega }_{c}}=50\;{\rm{c}}{{{\rm{m}}}^{-1}}$, and temperature T = 300 K. In all cases, the laser duration ${{t}_{1}}$ is set to 100 fs. In addition, in order to make the results converge within a reasonable computational time, we use a moderate time step ${\rm{d}}t=1{\rm{fs}}$ and average over a sufficiently-large number of trajectories, i.e. $N={{10}^{5}}$.

5.2.1. Laser-induced population dynamics

Figures 6(a) and (b) show the time-evolution of the site populations initialized by a laser pulse for strongly (J = 120 cm$^{-1}$) and weakly (J = 20 cm$^{-1}$) electronically coupled dimers, respectively. The laser carrier frequency is tuned at $\omega =13\;000\;{\rm{c}}{{{\rm{m}}}^{-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 = 0 fs to t = 100 fs), the populations on both sites oscillate with large amplitudes, for the laser field is nearly on resonance with the transitions between the ground state $|\epsilon _{0}^{(1)}\rangle$ and two eigenstates $|\epsilon _{1}^{(1)}\rangle$ and $|\epsilon _{2}^{(1)}\rangle$. After the laser is turned off at t = 100 fs, 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 (see (3)), however, dynamics in the single-exciton manifold shows strong J dependence. A comparison between figures 6(a) and (b) shows that although both systems evolve to reach a thermal equilibrium due to the system–bath couplings, the relaxation dynamics is 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 400 fs, and the system then quickly relaxes to equilibrium. Clearly, the coherent energy relaxation is extremely efficient in this case and the photoexcitation dynamics and the coherent dynamics are intertwined [42, 54]. 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.

Zoom In Zoom Out Reset image size

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 $\omega =13\;400\;{\rm{c}}{{{\rm{m}}}^{-1}}$. Because the laser is largely detuned from the transition between $|\epsilon _{0}^{(1)}\rangle$ and $|\epsilon _{1}^{(1)}\rangle$ but nearly on resonance with the transition between $|\epsilon _{0}^{(1)}\rangle$ and $|\epsilon _{2}^{(1)}\rangle$, i.e. ${{g}_{01}}\ll |\epsilon _{1}^{\prime }-(\epsilon _{0}^{\prime }+\omega )|$ and ${{g}_{02}}\sim |\epsilon _{2}^{\prime }-(\epsilon _{0}^{\prime }+\omega )|$, we expect strong Rabi oscillations for the $|\epsilon _{0}^{(1)}\rangle \rightleftarrows |\epsilon _{2}^{(1)}\rangle$ transition, and $|\epsilon _{1}^{(1)}\rangle$ cannot 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 a 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 $|\epsilon _{2}^{(1)}\rangle$.

Zoom In Zoom Out Reset image size

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 = 100 fs. 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 effort must be made 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.

5.2.2. Entanglement dynamics

In 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–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-qubit system with arbitrary density matrix ρ, the concurrence is defined as [63]

$$\begin{eqnarray}&&C\equiv {\rm max} \left( 0,{{\lambda }_{1}}-{{\lambda }_{2}}-{{\lambda }_{3}}-{{\lambda }_{4}} \right),\end{eqnarray} \tag{ 53 }$$

where ${{\lambda }_{j}}{\rm{s}}$ are the square roots of eigenvalues of the matrix $\rho \tilde{\rho }$ in descending order. Here,

$$\begin{eqnarray}&&\tilde{\rho }=\left( {{\sigma }_{y}}\otimes {{\sigma }_{y}} \right){{\rho }^{*}}\left( {{\sigma }_{y}}\otimes {{\sigma }_{y}} \right),\end{eqnarray} \tag{ 54 }$$

with ${{\rho }^{*}}$ 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 monotonically increases as the state is more quantum-mechanically entangled. An analytical expression that yields C(t) from the matrix elements of $\rho (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

$$\begin{eqnarray}\rho =\left( \begin{array}{ccccccccccccccc} 0 & 0 & 0 & 0 \\ 0 & {{\rho }_{22}} & {{\rho }_{21}} & {{\rho }_{20}} \\ 0 & {{\rho }_{12}} & {{\rho }_{11}} & {{\rho }_{10}} \\ 0 & {{\rho }_{02}} & {{\rho }_{01}} & {{\rho }_{00}} \\ \end{array} \right),\end{eqnarray} \tag{ 55 }$$

where the basis states are $|0\rangle =|g{{\rangle }_{1}}|g{{\rangle }_{2}}$, $|1\rangle =|e{{\rangle }_{1}}|g{{\rangle }_{2}}$, $|2\rangle =|g{{\rangle }_{1}}|e{{\rangle }_{2}}$ and $|3\rangle =|e{{\rangle }_{1}}|e{{\rangle }_{2}}$. Thus, the concurrence is simplified as

$$\begin{eqnarray} C\left( \rho \right) & = & {\rm max} \left\{ 0,\sqrt{{{\rho }_{11}}{{\rho }_{22}}}+\left| {{\rho }_{21}} \right|-\left( \sqrt{{{\rho }_{11}}{{\rho }_{22}}}-\left| {{\rho }_{21}} \right| \right) \right\} \\ {} & = & 2\left| {{\rho }_{21}} \right|. \\ \end{eqnarray} \tag{ 56 }$$

Because ${{\rho }_{11}}+{{\rho }_{22}}\geqslant 2|{{\rho }_{21}}|$, the more population is excited to the single-excitation states, the larger the concurrence could be. In addition, at some time, when ${{\rho }_{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 100 fs is 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 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Besides, we notice that for all cases there is an abrupt change at t = 100 fs due to the discontinuity of the system–field interaction, which is an artifact of the step-function 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 [14, 28].

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. $\lambda >J>\Delta \epsilon$. 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 adds 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 [32] 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}_{{\rm{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 [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 exploration in the microscopic world [66–68]. Its amazing power manifests in manipulating quantum dynamics by systematic control design methods [69–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 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 pulse-phase affect EET dynamics in photosynthetic systems.

We thank stimulating discussions with Jian Ma, Yao Shen and Yi-Hsien Liu. YCC thanks Akihito Ishizaki for kindly providing the HEOM results for FMO. QA thanks the National Science Council, Taiwan (Grant No. NSC 101-2811-M-002-148) for financial support. YCC thanks the National Science Council, Taiwan (Grant No. NSC 100-2113-M-002-008-MY3), National Taiwan University (Grant No. 103R891305) and Center for Quantum Science and Engineering (Subproject: 103R891401) for financial support. We are grateful to Computer and Information Networking Center, National Taiwan University and the National Center for High-performance Computing, Taiwan for the support of high-performance computing facilities.

In section 2.2, we presented the revised CMRT for the case with a step-function pulse. However, in a realistic case, one more often encounters a Gaussian pulse. Because the NMQJ approach propagates the dynamics in discrete time intervals, it is trivial to generalize the short-time square pulse propagation scheme given in section 2.2 to treat the influence of a laser pulse with arbitrary pulse profiles.

In figure A1, we schematically demonstrate how to decompose a Gaussian pulse with profile $g(t)={{g}_{0}}{\rm exp} [-2{{(t-{{t}_{0}})}^{2}}/{{T}^{2}}]$ into ${{N}_{p}}$ rectangular sub-pulses with equal width $\Delta T$. The height of each sub-pulse is set to the strength of the Gaussian pulse right at the middle of the sub-pulse. The width $\Delta T$ is determined in such a way that by increasing the number of sub-pulses and keeping the duration $[{{t}_{0}}-T,{{t}_{0}}+T]$ fixed, i.e. ${{N}_{p}}={{2}^{n}}+1$ ($n=1,2,\cdots$) and ${{N}_{p}}\Delta T=2\;T$, the area difference between the Gaussian pulse and the total summation of all sub-pulses does not exceed a pre-determined accuracy threshold. During the time propagation, when the Hamiltonian switches to the case with a different sub-pulse, M new trajectories are added to the propagation. At the end of the sub-pulse sequence, there will be ${{N}_{p}}M$ states to propagate in the simulation, which still maintains the linear scaling of the method. For instance, when the Gaussian pulse is decomposed into ${{N}_{p}}=513$ rectangular sub-pulses, the relative error between the areas of the discretized pulse and the original Gaussian pulse is less than ${{10}^{-6}}$, and the total number of states in the propagation only increases linearly to $514\;M+1$. We remark that the above decomposition procedure holds for arbitrary pulse shapes.

Zoom In Zoom Out Reset image size

In this appendix, we show that (13) can be recast into the generalized Lindblad form. First of all, inspired by the discussion in [31, 32], we shall specifically show that the pure-dephasing term in the CMRT can be rewritten in the Lindblad form:

$$\begin{eqnarray} -\mathop{\mathop{\sum }}\limits_{k}\frac{{{\Gamma }_{k}}}{2}\left[ \left\{ A_{kk}^{\dagger }{{A}_{kk}},\rho \right\}-2{{A}_{kk}}\rho A_{kk}^{\dagger } \right] & = & -\mathop{\mathop{\sum }}\limits_{k\ne {{k}^{\prime }}}\frac{{{\Gamma }_{k}}+{{\Gamma }_{{{k}^{\prime }}}}}{2}{{\rho }_{k{{k}^{\prime }}}}\left| {{\epsilon }_{k}} \right\rangle \left\langle {{\epsilon }_{{{k}^{\prime }}}} \right| \\ {} & = & -\mathop{\mathop{\sum }}\limits_{k\ne {{k}^{\prime }}}R_{k{{k}^{\prime }}}^{{\rm{pd}}}{{\rho }_{k{{k}^{\prime }}}}\left| {{\epsilon }_{k}} \right\rangle \left\langle {{\epsilon }_{{{k}^{\prime }}}} \right|, \\ \end{eqnarray} \tag{ B.1 }$$

where

$$\begin{eqnarray}&&R_{k{{k}^{\prime }}}^{{\rm{pd}}}\left( t \right)=\frac{1}{2}\left( {{\Gamma }_{k}}+{{\Gamma }_{{{k}^{\prime }}}} \right).\end{eqnarray} \tag{ B.2 }$$

Here, we have made use of the properties of the diagonal jump operators, i.e.

$$\begin{eqnarray}&&A_{kk}^{\dagger }{{A}_{kk}}=A_{kk}^{\dagger }={{A}_{kk}}.\end{eqnarray} \tag{ B.3 }$$

Clearly, the elements of the pure dephasing rate matrix, $R_{k{{k}^{\prime }}}^{{\rm{pd}}}(t)$, are real and symmetric with respect to the matrix diagonal. In addition, the diagonal terms vanish, i.e.

$$\begin{eqnarray}R_{k{{k}^{\prime }}}^{{\rm{pd}}}=\left\{ \begin{array}{ccccccccccccccc} {\rm{real}}, & {\rm{for}}\;k\ne {{k}^{\prime }}, \\ 0, & {\rm{for}}\;k={{k}^{\prime }}. \\ \end{array} \right.\end{eqnarray} \tag{ B.4 }$$

Since there are $M(M-1)/$2 independent pure-dephasing rates $R_{k{{k}^{\prime }}}^{{\rm{pd}}}$ but M Lindblad-form dephasing rates ${{\Gamma }_{k}}$, we need some additional constraints in order to fit M Lindblad-form dephasing rates ${{\Gamma }_{k}}$ from $M(M-1)/$2 independent pure-dephasing rates $R_{k{{k}^{\prime }}}^{{\rm{pd}}}$. In this work, we propose to use a least-square fit to obtain the Lindblad-form dephasing rates. Therefore, we require the mean-square displacement to be minimal with respect to all the Lindblad rates:

$$\begin{eqnarray}&&\frac{\partial }{\partial {{\Gamma }_{a}}}\mathop{\mathop{\sum }}\limits_{k=1}^{M-1}\mathop{\mathop{\sum }}\limits_{{{k}^{\prime }}=k+1}^{M}{{\left[ R_{k{{k}^{\prime }}}^{{\rm{pd}}}-\frac{1}{2}\left( {{\Gamma }_{k}}+{{\Gamma }_{{{k}^{\prime }}}} \right) \right]}^{2}}=0,\end{eqnarray} \tag{ B.5 }$$

where $a=1,2,\ldots M$. After simplification, we obtain a system of linear equations for the Lindblad-form dephasing rates, that is

$$\begin{eqnarray}&&\frac{1}{2}\mathop{\mathop{\sum }}\limits_{k=1}^{a-1}{{\Gamma }_{k}}+\frac{1}{2}\left( 2\;M-a \right){{\Gamma }_{a}}+\mathop{\mathop{\sum }}\limits_{k=a+1}^{M-1}{{\Gamma }_{k}}+\frac{1}{2}{{\Gamma }_{M}}={{A}_{a}},\end{eqnarray} \tag{ B.6 }$$

where the coefficient on the right hand side is

$$\begin{eqnarray}&&{{A}_{a}}=\mathop{\mathop{\sum }}\limits_{k=a+1}^{M}R_{ak}^{{\rm{pd}}}+\mathop{\mathop{\sum }}\limits_{k=1}^{M-1}R_{ka}^{{\rm{pd}}}.\end{eqnarray} \tag{ B.7 }$$

Or equivalently, (B.6) can be written in a matrix form as

$$\begin{eqnarray}&&B\Gamma =A,\end{eqnarray} \tag{ B.8 }$$

with the matrix elements of B given by

$$\begin{eqnarray}{{B}_{jk}}=\left\{ \begin{array}{ccccccccccccccc} \frac{1}{2}, & {\rm{for}}\;k<j, \\ \frac{1}{2}\;\left( 2\;M-j \right), & {\rm{for}}\;k=j, \\ 1, & {\rm{for}}\;j<k<M, \\ \frac{1}{2}, & {\rm{otherwise}}. \\ \end{array} \right.\end{eqnarray} \tag{ B.9 }$$

Therefore, by multiplying both sides with the inverse of B, we have

$$\begin{eqnarray}&&\Gamma ={{B}^{-1}}A.\end{eqnarray} \tag{ B.10 }$$

Therefore, the original CMRT master equation (cf (13)) can be written in the Lindblad form as

$$\begin{eqnarray} {{\partial }_{t}}\rho & = & -{\rm{i}}\left[ {{H}_{e}}\left( t \right),\rho \right]-\frac{1}{2}\mathop{\mathop{\sum }}\limits_{k\ne {{k}^{\prime }}}R_{k{{k}^{\prime }}}^{{\rm{dis}}}\left( t \right)\left[ \left\{ A_{k{{k}^{\prime }}}^{\dagger }{{A}_{k{{k}^{\prime }}}},\rho \right\}-2{{A}_{k{{k}^{\prime }}}}\rho A_{k{{k}^{\prime }}}^{\dagger } \right] \\ {} & {} & -\frac{1}{2}\mathop{\mathop{\sum }}\limits_{k}{{\Gamma }_{k}}\left( t \right)\left[ \left\{ A_{kk}^{\dagger }{{A}_{kk}},\rho \right\}-2{{A}_{kk}}\rho A_{kk}^{\dagger } \right], \\ \end{eqnarray} \tag{ B.11 }$$

as given in (37).