MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\$','\$']]}}); Long-Lived Coherence Originating from Electronic-Vibrational Couplings in Light-Harvesting Complexes
 Chinese Journal of Chemical Physics  2017, Vol. 30 Issue (2): 186-192

#### The article information

Xian-ting Liang

Long-Lived Coherence Originating from Electronic-Vibrational Couplings in Light-Harvesting Complexes

Chinese Journal of Chemical Physics, 2017, 30(2): 186-192

http://dx.doi.org/10.1063/1674-0068/30/cjcp1609188

### Article history

Accepted on: November 14, 2016
Long-Lived Coherence Originating from Electronic-Vibrational Couplings in Light-Harvesting Complexes
Xian-ting Liang
Dated: Received on September 25, 2016; Accepted on November 14, 2016
Department of Physics and Institute of Optics, Ningbo University, Ningbo 315211, China
Author: Xian-ting Liang, E-mail:liangxianting@nbu.edu.cn
Abstract: We theoretically investigate the evolutions of two-dimensional, third-order, nonlinear pho-ton echo rephasing spectra with population time by using an exact numerical path integral method. It is shown that for the same system, the coherence time and relaxation time of excitonic states are short, however, if the couplings of electronic and intra-pigment vibra-tional modes are considered, the coherence time and relaxation time of this vibronic states are greatly extended. It means that the couplings between electronic and vibrational modes play important roles in keeping long-lived coherence in light-harvesting complexes. Particularly, by using the method we can fix the transition path of the energy transfer in bio-molecular systems.
Key words: Light-harvesting complex     Coherence     Exciton model     Vibronic model
Ⅰ. INTRODUCTION

The significance of photosynthesis is indisputable because it is the fundamental biological process through which solar energy is converted into biomass [1], and it provides chemical energy for plants, algae, and bacteria, while heterotrophic organisms rely on these species as their ultimate food source. The initial step in photosynthesis involves the capture of energy from sunlight. It is followed by the excitation energy transfer (EET) to the reaction center, where charge separation is initiated and physical energy is converted into chemical energy. At low light intensities, the quantum efficiency of the transfer of energy is very high, with a near unity quantum yield.

There is a hypothesis that the high efficiency of near unity EET attributes to coherent superposition between excitation states that belongs to different chromophores of some photosynthesis system. It is considered that quantum coherence between the electronic states involved in energy or electron transfer introduces correlations between the wave functions of these states enabling the excitations to move rapidly (by virtue of wave packet collapse) and to sample multiple pathways (due to different coherent superposition states) in space. Thus, the quantum coherence effects may render the process of effective energy and electron transfer [2, 3].

Indeed, several experimental results show that there is coherence with long enough time in light-harvesting complexes which is inconceivable in inorganic materials at physiological temperature. Fleming et al . [4-6] found that in the Fenna-Matthew-Olson (FMO) complex the coherence between the two excitation states lasts for longer time which is similar to the time scale of the EET, and specially, in a photosynthetic purple bacterium Rb . Sphaeroides , the coherence of the excitations created from laser pulses on adjacent bacteriochlorophyll (Bchls) and bacteriopheophytin (BPhys) can prolong for longer time which is long enough for the coherent transfer of energies between the chromophores Bchls and BPhys [6]. These kinds of experimental results imply that electronic excitations may move coherently through these photosynthesis systems rather than by incoherent hopping motion as it has usually been assumed [7, 8].

However, how to probe the coherence of these coherent superposition states? Theoretically, one can measures the degree of the coherence through quantifying the off-diagonal elements of the density matrix of the states. And experimentally, the coherence can be detected through measuring oscillations of the photon echo rephasing and non-rephasing signals [9-11]. In Ref. [12] can help us to understand this problem. This kind of experimental results can be verified theoretically by using the simulating method of 2D spectroscopy [13, 14] based on the quantum dissipative dynamics [15-18]. Initially it was conjectured that long-lived oscillations are a manifestation of purely electronic coherence, however, this appeared to require unrealistically low reorganization energies of the vibrational environment or the existence of correlated fluctuations, both of which appear unlikely to be correct. Recent studies suggested that vibrational motion may modulate both diagonal and off-diagonal features of the 2D spectra.

The electronic, vibrational, and electronic-vibrational [19-23] coherence has been extensively investigated by using the nonlinear 2D spectroscopy [23-27]. Recently, Chin et al . [28] proposed a mechanism, which shows that the resonant electronic-vibrational configuration can sustain, regenerate, or even create coherence between electronic states during the time scale of energy and electron transfer. Lim et al . [12] theoretically verify the experimental long-lasting coherence of an artificial light harvester by using the exciton-vibrational (vibronic) model.

These efforts, however, still cannot ascertain: whether the couplings between the electronic and vibrational modes must lead to coherence time expanding? In this work, we design a theoretical scheme to accurately answer this question. We plot the evolutions of the two-dimensional (2D), third-order, nonlinear photon echo spectra of two models for the same system. One of which is the exciton model, while the other is the vibronic model, namely the electronic-vibrational coupling model. It will be shown that the diagonal peaks and the cross-peaks of the 2D, third-order, nonlinear photon echo spectra calculated from the vibronic model can survive longer time than that obtained from the exciton model as other conditions are kept unchanged.

Thus, our investigations support that the long-lived coherence of the light-harvesting complex may origin from the couplings of the electronic-vibrational modes. In this study protocol the decay of the cross-peaks shows the decoherence, and the decay of the diagonal peaks reflects the relaxation of the system. In this investigation, we have not deliberately use the unrealistically low reorganization energies and correlated fluctuations of the environment. The evolution of the 2D, third-order, nonlinear photon echo rephasing spectra are calculated at lower temperature 77 K. Lower temperature is also applicable to ambient temperature.

Ⅱ. MODEL AND 2D, THIRD-ORDER, NONLINEAR PHOTON ECHO SPECTRA

According to Sato's report [29], if we assign the zero-mode of the nuclear vibrations to an intra-pigment vibration and exempt it from reorganization dynamics, and assign the rest modes to the continuum of the protein bath. The vibronic model can be written as

 (1)
 (2)
 (3)
 (4)

Here, H Sve describes the intrinsic dynamics of an exciton coupled to the vibrational modes and S is system. H B is the bath Hamiltonian, and the H SBve describes the coupling between a pigment and its bath and SB is system-bath. m =1, 2 indexes the exciton two modes and α , β , α ', β ' index the vibrational modes coupled to the electronic modes. ϵ m is the exciton site energy, and J 12 describes the coupling between exciton states. ωξ is the frequency of the ξ -th bath mode; C ξ , m and C ξ , m are creation and annihilation operators of the ξ -th mode of bath coupled to the m-th mode of the system; and g ξ , m is the coupling strength between the degree of freedom m of system and the bath mode ξ . (FC)α , β is the Franck-Condon factor [30], which have been obtained in Ref. [29], as:

 (5)

This is the vibronic model when ω 0≠0, and Huang-Rhys factor R ≠0, while it is the exciton model when ω 0=0, and R =0. Thus, the Hamiltonian H Sve for the exciton model can be written as

 (6)

and the Hamiltonian H Sve for the vibronic model truncated at α or β =1 can be described in Table Ⅰ.

Table Ⅰ H Sve for the vibronic model. Here and in the text we set |g 0〉=|1, 0, 0〉, |g 1〉=|1, 0, 1〉, |g 2〉=|1, 1, 0〉, |e 0〉=|2, 0, 0〉, |e 1〉=|2, 0, 1〉, |e 2〉=|2, 1, 0〉.

It is known that all information about the interaction of the system with bath is contained in the spectral density function which is defined as

 (7)

By fixing the environmental spectrum J m(ω ), described with some math functions or discrete data, we can investigate not only the dynamics but also the 2D, third-order, nonlinear photon echo spectroscopy by using some quantum dissipative dynamical methods. In this work, we are interested in the evolution features of the 2D, third-order, nonlinear photon echo spectra for the exciton and vibronic models in principle, so we simply choose the single excitation dipole operator. For the exciton model, it is

 (8)
 (9)
 (10)

And for the vibronic model, it is

 (11)
 (12)

Simply, but without loss of generality, here we set µ 12=µ 21=1, and use the Debye spectral density function

 (13)

Similarly to Ref. [29], here we set the parameters with more practical values, namely, λ =35 cm-1, γ -1=100 fs, ω 0=191 cm-1, J 12=100 cm-1, ϵ 1=191 cm-1, ϵ 2=382 cm-1, and the Huang-Rhys factor S =0.025. And throughout the paper we set the Planck constant h =1.

Ⅲ. RESPONSE FUNCTIONS AND 2D SPECTROSCOPY

In the 2D, third-order photon-echo experiments, three laser pulses with wave-vector k1, k2, k3 are projected to the sample, and the third-order nonlinear response signals are detected along the phase-matched direction ks, see Fig.S1 (a) and (b) in supplementary materials. In the impulsive limit, the rephasing and nonrephasing spectra are calculated as the real part of the following double Fourier-Laplace transforms of R r (t 3, t 2, t 1) and R nr (t 3, t 2, t 1) with respect to t 1 and t 3 as Ref. [31-35]

 (14)

Here, R r (t 3, t 2, t 1) and R nr (t 3, t 2, t 1) are the rephasing (ks=−k1+k2+k3) and nonrephasing (ks=k1-k2+k3) signals, and they can be calculated through

 (15)
 (16)

Here, G (t ) is the retarded propagator of the total system in Liouville space, µ ±×x =[µ ±, x ], and we suppose the total system in equilibrium state, initially, represented with ρ toteq(0). t 1=τ 2-τ 1, t 2=τ 3-τ 2, and t 3=τ -τ 3, where τ 1, τ 2, and τ 3are the time points of laser projecting to a sample, and τ is the time point of receiving echo signals. It is straightforward to calculate the response function of photon echo by using the usual quantum dissipative dynamics. The calculation of R r (t 3, t 2, t 1) can be performed as follows. At first, calculating ρ (t 1)≡G (t 1)(µ ×ρ toteq(0)) for t 1from 0to τ 1; next, ρ (t 2)≡G (t 2)(µ +× ρ (τ 1)) for t 2 from τ 1 to τ 2; then ρ (t 3)≡G (t 3)(µ +×ρ (τ 2)) for t 3from τ 2to τ 3; and finally, calculating R r (t 3, t 2, t 1) Tr [µρ 3(t 3, t 2, t 1)]. The only difference of calculations for dynamics and for responce functions is that the initial state ρ (0) in dynamics calculations has been changed into µ ±×ρ toteq (0), or µ ±×ρ 1(τ 1), µ ±×ρ 2(τ 2) in the calculations of response functions. Therefore, all schemes to calculate the reduced density matrix can be used to calculate the response functions. In the following, at first, we plot the 2D, third-order, nonlinear photon echo spectra, at fixed population time, for the exciton and vibronic models. The results are plotted in Fig. 1.

 FIG. 1 2D, third-order, nonlinear photon echo rephasing spectra at 77 K, for the exciton model (Ⅰ) at population time t 2=720 fs; and for the vibronic model (Ⅱ), at the population time t 2=324 fs. Here, the environment is described with Debye spectral density function, and the values of the environmental spectral parameters. The evolutions of cross-peaks are marked by arrow line (a) in Fig. 1(Ⅰ) for the exciton model, and arrow lines (a), (b), (c) in Fig. 1(Ⅱ) for the vibronic model. The evolutions of diagonal peaks are marked by double-arrow lines (b) in Fig. 1(Ⅰ), and (d) in Fig. 1(Ⅱ).

In the culculations of data for plotting Fig. 1 we choose the environmental temperature T =77 K. The calculations are performed with an exact numerical dynamics approach-the quasi-adiabatic propagator path integral (QUAPI) method developed by Makris group [36-40]. The theoretical framework for simulating the 2D, third-order, nonlinear photon echo spectra with the QUAPI method may refer to Refs. [41-44]. In this work, we are interested in the evolutions of the 2D, third-order, nonlinear photon echo rephasing spectra with the population time t 2. The figures in Fig. 1 are obtained from exciton (Ⅰ) and vibronic (Ⅱ) models. In Fig. 2 we plot the evolutions of the amplitude of cross-peak marked with arrow line (a) in Fig. 1(Ⅰ), and diagonal peaks marked with double-headed arrow line (b) in Fig. 1(Ⅰ) for the exciton model. The evolutions of cross-peaks marked by (a), (b), (c) in Fig. 1() for the vibronic model have been plotted in Fig. 3(a), (b), (c), and the evolutions of diagonal peaks marked by (d) in Fig. 1() have been plotted in Fig. 3(d). We obtain the following results: (ⅰ) for not only the exciton but also the vibronic models, the amplitudes of the cross peaks of the 2D, third-order, nonlinear photon echo rephasing spectra are oscillating with population time t 2; (ⅱ) in our choice parameters, the oscillations of the diagonal peaks of the 2D spectrum are not clear in the exciton model, however, they are clearly shown in the vibronic model; (ⅲ) the most important is that the oscillating times of both the diagonal and cross peaks of the 2D spectra obtained from the vibronic model are much longer than that from the exciton model.

 FIG. 2 The evolutions of the 2D, third-order, nonlinear photon echo rephasing spectra at 77 K for the exciton model. (a) The evolution of cross peak corresponding to (a) in Fig. 1(Ⅰ) , and (b) the evolution of the diagonal peak corresponding to double-arrow line (b) in Fig. 1(Ⅰ) .
 FIG. 3 The evolutions of the 2D, third-order, nonlinear photon echo rephasing spectra at 77 K for the vibronic model. (a), (b), and (c) are the evolutions of cross peaks corresponding to (a), (b), and (c) in Fig. 1(Ⅱ) , and (d) is the evolution of the diagonal peaks corresponding to (d) in Fig. 1(Ⅱ) .
Ⅳ. DISCUSSION AND CONCLUSION

We theoretically discuss the numerical results. The eigenvalues of H Sve for the exciton model, see Eq.(6), are approximatively equal to E 1=-42.78 cm-1, and E 2=233.78 cm-1, from which we can obtain the frequency ω =(E 1-E 2)/h ≈1737 cm-1 (here h =1/2π ), and from Fig. 2 we can read out the oscillating period τ ≈120 fs, then ω '=2π /τ ≈1745 cm-1. It shows that they are in good agreement with each other. From Table Ⅰ we can obtain the eigenvalues of the vibronic model, they are E 1=-39.25 cm-1, E 2=417.68 cm-1, E 3=230.21 cm-1, E 4=155.28 cm-1, E 5=423.00 cm-1, and E 6=150.00 cm-1. Similarly, from Fig. 3 we see that the oscillation frequency is about 513 cm-1, and it corresponds to the transitions |e 0〉〈g 2|, see Table S2, Fig.S1(g) in supplementary materials.

We also analyse the process according to the Feynman diagrams shown in Fig. 4, other two kinds of cases are shown in Fig.S1 (c) and (e) in supplementary materials. The 2D, third-order, nonlinear photon echo rephasing spectra are the result of Fourier transform with time t 1, and inverse Fourier transform with time t 3. For the diagonal peaks, the Feynman diagram of one kind of pathways is shown in Fig. 4(a). We set the initial state of the system be |g 0〉〈g 0|. At t =0, the first laser pulse with ω β is projected, the system becomes the coherent superposition state, for example |e 1〉〈g 0|, and then the system evolves for time t 1. At t =τ 1, the second pulse is projected, the system turns into the population state |e 1〉〈e 1|, and subsequently relaxes for time t 2. At t =τ 2, the third laser pulse, with frequency ω β is projected, and the system again becomes a coherent superposition state |e 1〉〈g 0|, and then it evolves for time t 3. Fig. 2(b) and Fig. 3(d) are the evolutions of the diagonal peaks for the 2D, third-order, nonlinear photon echo rephasing spectra with time t 2for the exciton and vibronic models.

For the cross peaks, the Feynman diagram of one kind of pathways is shown in Fig. 4(b), other two kinds of cases are shown in Fig.S1 (d) and (f) in supplementary materials. The initial state of the system is still set as |g 0〉〈g 0|. At t =0, the first laser pulse with ω α is projected, the system becomes the coherent superposition state, for example |e 0〉〈g 0|, and then it evolves till t =τ 1, the second pulse is projected, the system turns into the state |e 0〉〈e 1|, or |e 0〉〈e 2| till t =τ 2, and then at the third laser pulse, frequency ω α is projected. Subsequently, the system again becomes coherent superposition state, |e 0〉〈g 0|, and evolves for time t 3. Fig. 2(a), and Fig. 3(a), (b), (c) are the evolutions of the cross-peaks of the 2D, third-order, nonlinear photon echo rephasing spectra with time t 2for the exciton and vibronic models. From Fig. 2(a) we see that the oscillation frequency is about 1745 cm-1, it just corresponds to the transition between state |g 〉and |e 〉. From Fig. 3 we see that the oscillation frequency finally becomes at about 513 cm-1, it means that the main transitions are |e 2〉↔|g 2〉(509 cm−1), see Table S2 in supplementary materials. Similarly, for any practical system we can read out the main pathways of energy transfer by using the evolutions of the 2D, third-order, nonlinear photon echo spectra. In above investigations we set the vibrational transition ω 0=191 cm-1, which is not quite high compared to the thermal vibration at temperature 77 and 300 K, but it is large enough compared to the thermal vibration at 30 K in Fig.S2 in supplementary materials. If we select the vibrational transition ω 0=300 cm-1, as Ref. [45], the cnculsions obtained are all still held.

 FIG. 4 Double-side Feynman diagrams representing the pathways contributing to a diagonal peak (a) and to a cross peak (b).

In summary, in this work, we have theoretically investigated the role of the couplings of the vibrational-electronic modes. We use the comparative data of the evolutions of 2D, third-order, nonlinear photon echo rephasing spectra, calculated from the exciton and vibronic models. We carefully analyze the single excitation exciton and vibronic models at temperature 77 K. In supplementary materials, we give the results of 300 and 30 K. It is shown that the coherence exists indeed at ambient temperature, and the conclusions obtained at 77 K are kept at ambient temperature, but the coherent time becomes shorter at 300 K than at 77 K. We clearly show that when the vibrational modes coupled to the electronic modes, the electronic coherence time becomes considerable long compared to the cases without couplings. This result is in agreement with the corresponding experiments and supports the theory about long-lived coherence for some light-harvesting systems [46-48]. However, other mechanisms of the long-lived coherence [6, 49, 50], such as correlated fluctuations of the environment, are not ruled out yet. We suggest that it may be coexistence of a variety of mechanisms for long-lived coherence in light-harvesting complexes.

Supplementary materials: Tables S1 and S2 show the eigenvalues, and transition scheme and frequencies for the vibronic model. Table S3 shows the oscillation periods of the evolution of the 2D, third-order, non-linear photon echo spectra for the vibronic model. Figure S1 shows the laser excitation pulses, non-colinear geometrical configuration, two kinds of Feynman diagrams for diagonal peaks and and cross-peaks, the energy levels, and their transition diagrams for the vibronic model. Figure S2 and S3 show the evolution figures of the 2D, third-order, non-linear photon echo rephasing spectra at 30 and 300 K for the exciton and vibronic models.

Ⅴ. ACKNOWLEDGMENTS

This work was supported by the Zhejiang Provincial Natural Science Foundation of China (No.LY13A040006), and the K. C. Wong Magna Foundation in Ningbo University.

 [1] F. Müh, M.E. A. Madjet, J. Adolphs, A. Abdurahman, B. Rabenstein, H. Ishikita, E. W. Knapp, and T. Renger, Proc. Natl. Acad. Sci. 104 , 16862 (2007). DOI:10.1073/pnas.0708222104 [2] T. Renger, V. May, and O. Kuhn, Phys. Rep. 343 , 137 (2001). DOI:10.1016/S0370-1573(00)00078-8 [3] S. Huelga, and M. Plenio, Contemp. Phys. 54 , 181 (2013). DOI:10.1080/00405000.2013.829687 [4] G. S. Engel, T. R. Calhoun, E. L. Read, T. K. Ahn, T. Mančal, Y. C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446 , 782 (2007). DOI:10.1038/nature05678 [5] T. Brixner, J. Stenger, H. M. Vaswani, M. Cho, R. E. Blankenship, and G. R. Fleming, Nature 434 , 625 (2005). DOI:10.1038/nature03429 [6] H. Lee, Y. C. Cheng, and G. R. Fleming, Science 316 , 1462 (2007). DOI:10.1126/science.1142188 [7] M. Mohseni, P. Rebentrost, S. Lloyd, and A. Aspuru-Guzik, J. Chem. Phys. 129 , 174106 (2008). DOI:10.1063/1.3002335 [8] P. Rebentrost, and A. Aspuru-Guzik, J. Chem. Phys. 134 , 101103 (2011). DOI:10.1063/1.3563617 [9] E. E. Fenn, D. B. Wong, and M. D. Fayer, J. Chem. Phys. 134 , 054512 (2011). DOI:10.1063/1.3532542 [10] Y.C. Cheng, and G. R. Fleming, J. Phys. Chem A112 , 4254 (2008). [11] M. Cho, H. M. Vaswani, T. Brixner, J. Stenger, and G. R. Fleming, J. Phys. Chem B109 , 10542 (2005). [12] J. Lim, D. Paleček, F. Caycedo-Soler, C. N. Lincoln, J. Prior, H. von Berlepsch, S. F. Huelga, M. B. Plenio, D. Zigmantas, and J. Hauer, Nature Comm. 6 , 8755 (2015). DOI:10.1038/ncomms9755 [13] S. Mukamel, Principles of Nonlinear Optical Spectroscopy . New York: Oxford University Press (1995). [14] P. Hamm, and M. Zanni, Concepts and Methods of 2D Infrared Spectroscopy . Cambridge: Cambridge University Press (2011). [15] U. Weiss, Quantum Dissipative Systems , Fourth edition. Singapore: World Scientific Publishing Co. Pte. Ltd. (2012). [16] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Mod. Phys. 59 , 1 (1987). DOI:10.1103/RevModPhys.59.1 [17] A.O. Caldeira, and A. J. Leggett, Phys A121 , 587 (1983). [18] C.P. Sun, and L. H. Yu, Phys. Rev A51 , 1845 (1995). [19] S. Mukamel, Annu. Rev. Phys. Chem. 51 , 691 (2000). DOI:10.1146/annurev.physchem.51.1.691 [20] Y. Jing, R. Zheng, H. X. Li, and Q. Shi, J. Phys. Chem B116 , 1164 (2012). [21] D. B. Turner, K. E. Wilk, P.M. G. Curmi, and G. D. Scholes, J. Phys. Chem. Lett. 2 , 1904 (2011). DOI:10.1021/jz200811p [22] D. Hayes, and G. S. Engel, Philos. Trans. 370 , 3692 (2012). DOI:10.1098/rsta.2011.0201 [23] M. T. Zanni, N. H. Ge, Y. S. Kim, and R. M. Hochstrasser, Proc. Natl. Acad. Sci. U. S. A. 98 , 11265 (2001). DOI:10.1073/pnas.201412998 [24] J. Xu, H. D. Zhang, R. X. Xu, and Y. Yan, J. Chem. Phys. 138 , 244 (2013). [25] Y. Tanimura, and S. Mukamel, J. Chem. Phys. 99 (1993). [26] L. Chen, R. Zheng, Y. Jing, and Q. Shi, J. Chem. Phys. 134 , 194508 (2011). DOI:10.1063/1.3589982 [27] W. G. Noid, G. S. Ezra, and R. F. Loring, J. Chem. Phys. 119 , 1003 (2003). DOI:10.1063/1.1577319 [28] A. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. Huelga, and M. Plenio, Nature Physics 9 , 113 (2013). DOI:10.1038/nphys2515 [29] Y. Sato, and B. Doolittle, J. Chem. Phys. 141 , 185102 (2014). DOI:10.1063/1.4901056 [30] V. May, and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems . 2nd Edn New York: John Wiley & Sons (2008). [31] A. Ishizaki, and Y. Tanimura, Chem. Phys. 347 , 185 (2008). DOI:10.1016/j.chemphys.2007.10.037 [32] A. Ishizaki, and Y. Tanimura, J. Chem. Phys. 125 , 084501 (2006). DOI:10.1063/1.2244558 [33] L. Chen, R. Zheng, Q. Shi, and Y. Yan, J. Chem. Phys. 132 , 024505 (2010). DOI:10.1063/1.3293039 [34] Y. Tanimura, arXiv preprint arXiv , 1211–1469 (2012). [35] A. Ishizaki, and Y. Tanimura, J. Phys. Chem A111 , 9269 (2007). [36] D.E. Makarov, and N. Makri, Chem. Phys. Lett. 221 , 482 (1994). DOI:10.1016/0009-2614(94)00275-4 [37] N. Makri, J. Mathe. Phys. 36 , 2430 (1995). DOI:10.1063/1.531046 [38] E. Sim, and N. Makri, Comput. Phys. Commun. 99 , 335 (1997). DOI:10.1016/S0010-4655(96)00130-0 [39] E. Sim, and N. Makri, J. Phys. Chem B101 , 5446 (1997). [40] N. Makri, J. Chem. Phys. 141 , 134117 (2014). DOI:10.1063/1.4896736 [41] M.M. Sahrapour, and N. Makri, J. Chem. Phys. 132 , 134506 (2010). DOI:10.1063/1.3336463 [42] X. T. Liang, J. Chem. Phys. 141 , 044116 (2014). DOI:10.1063/1.4890533 [43] X. T. Liang, Chin. Phys. Lett. 29 , 084211 (2012). DOI:10.1088/0256-307X/29/8/084211 [44] H. G. Duan, A. G. Dijkstra, P. Nalbach, and M. Thorwart, Phys. Rev E92 , 042708 (2015). [45] J. Xu, R. Xu, D. Abramavicius, H. Zhang, and Y. Yan, Chin. J. Chem. Phys. 24 , 497 (2011). DOI:10.1088/1674-0068/24/05/497-506 [46] K. Avinash, E.J. O'Reilly, G. D. Scholes, and O. C. Alexandra, J. Chem. Phys. 137 , 174109 (2012). DOI:10.1063/1.4764100 [47] A. W. Chin, S. F. Huelga, and M. B. Plenio, Philos. Trans. A Math. Phys. Eng. Sci. 370 , 3638 (2012). DOI:10.1098/rsta.2011.0224 [48] T. Vivek, W. K. Peters, and D. M. Jonas, Proc. Natl. Acad. Sci. USA. 110 , 1203 (2012). [49] I. Akihito, T. R. Calhoun, G.S. Schlau-Cohen, and G. R. Fleming, Phys. Chem. Chem. Phys. 12 , 7319 (2010). DOI:10.1039/c003389h [50] H. Dugan, G. B. Griffin, and G. S. Engel, Science 340 , 1431 (2013). DOI:10.1126/science.1233828