Using non-Markovian measures to evaluate quantum master equations for photosynthesis

When dealing with system-reservoir interactions in an open quantum system, such as a photosynthetic light-harvesting complex, approximations are usually made to obtain the dynamics of the system. One question immediately arises: how good are these approximations, and in what ways can we evaluate them? Here, we propose to use entanglement and a measure of non-Markovianity as benchmarks for the deviation of approximate methods from exact results. We apply two frequently-used perturbative but non-Markovian approximations to a photosynthetic dimer model and compare their results with that of the numerically-exact hierarchy equation of motion (HEOM). This enables us to explore both entanglement and non-Markovianity measures as means to reveal how the approximations either overestimate or underestimate memory effects and quantum coherence. In addition, we show that both the approximate and exact results suggest that non-Markonivity can, counter-intuitively, increase with temperature, and with the coupling to the environment.

or time-convolutionless equation. Some works suggest that TL2 shows better performance than TC2 at numerically approximating exact results 28 . Nevertheless, their respective domains of applicability have not been thoroughly investigated yet.
In each QME model (TC2, TL2), certain approximations and simplifications are introduced to obtain solvable equations. To investigate the deviation of each approximate QME model from the exact results, we first compare the explicit dynamics of these two approximative QMEs with that of the HEOM. The HEOM approach is considered to be numerically exact for the models with the Drude-Lorentz spectral density function investigated here. For more general bath models and dynamics at low temperatures, the need to truncate at a certain level of the hierarchy equations can lead to errors, and thus the exactness of the HEOM approach requires further scrutiny in such cases [29][30][31] . We focus on the intermediate system-environment coupling regime, which has proven to be the most challenging and relevant to the dynamics in realistic systems such as the photosynthetic Fenna-Matthews-Olson complex. Notably, the intermediate regime is also the one at which the region of validity of most approximations breaks down. Both approximate methods are perturbative in the system-bath coupling, but can in principle harbor memory effects of the environment.
Recently, much effort has been devoted to the quantification of memory effects [32][33][34][35] which has subsequently been studied in the context of various physical systems [36][37][38] . To investigate how well the models we study here capture the memory effect, we utilize the concept of the Choi-Jamiołkowski isomorphism 39,40 to encode complete information on the dynamics of the system into the entanglement with an ancilla. By comparing the time evolution of the entanglement between system and ancilla, and an associated measure of non-Markovianity 33 , one can find out to what extent the memory effects and coherence predicted by each approximate QME deviates from being numerically exact. Our results suggest that entanglement and non-Markovianity provide a useful benchmark for the performance of such approximative treatments, providing a more fine-grained insight into the deviation from exact results than quantities like the fidelity alone.
In performing this analysis we also discuss several interesting physical trends, including a counter-intuitive increase of non-Markovianity with both temperature and with the coupling strength to the environment. We attribute this increase to an enhancement of system-environment correlations when both the coupling and temperature are increased. Additionally, evidence from other studies [41][42][43] suggests that non-Markovian environments are capable of sustaining quantum coherence. The interplay of these factors finally results in the increase of non-Markovianity with both temperature and coupling strength that we see in our results.

Results
The spin-boson model. The spin-boson model 1 is one of the most extensively studied models of open quantum systems, and is the one we employ here. It describes a spinor-like two-state system interacting with a bosonic environment. First, let us consider this standard model, which can be divided into three components here J x σ is the coherent-coupling term, which enables the tunneling between the two system quantum states, labeled as |1〉 and |− 1〉 , with the energy level spacing ħω 0 . Usually, one adopts the delocalized basis |χ + 〉 and |χ − 〉 (exciton), which is defined by the following eigenvalue problem The environment, H env , is usually modelled as a large collection of harmonic oscillators where a k † (a k ) is the creation (annihilation) operator of the environment mode k with angular frequency ω k . For simplicity, a linear system-environment coupling, H int , is adopted throughout this work: where g k is the coupling constant between the environment mode k and the system. In most physical problems, the details of the microscopic description of g k are not clear, and one usually employs a spectral density function, J(ω) = ∑ k |g k | 2 δ(ω − ω k ), to characterize the coupling strength via the reorganization energy . The physical meaning of the spectral density function can be understood as the density of states of the environment, weighted by the coupling strengths. Moreover the way in which the environment modulates the dynamics of the system is described by the correlation function The real part is related to the dissipation process, while the imaginary part corresponds to the response function.
The statistical properties of the entire system can be described by the total density matrix ρ tot , which contains all the degrees of freedom of the system and environment. If the correlation between the system and environment is negligible, the Born approximation can be used and the total density matrix can be factorized into where ρ sys (t) describes the dynamics of the system and ρ = − is the environment density matrix in thermal equilibrium at temperature T. Here, k B is the Boltzmann constant and is the partition function. One notes that when ω 0 , J, and λ are comparable, this makes the conventional perturbative treatment unreliable. In the following, we will adopt the two frequently-used perturbative but non-Markovian QME formalisms discussed in the introduction and compare their results with the exact one in the intermediate-coupling regime, as they both begin to break down, and investigate ways in which to evaluate their accuracy.

Second-order time-convolution equation (TC2).
For the Hamiltonian defined above, the time evolution of the system density matrix ρ tot (t) under the TC2 approximation is expressed as ,ˆ, and μ, ν = χ + , χ − . Substituting Eq. (9) into (8) with the explicit expansion leads to a set of simultaneous integrodifferential equations of the density matrix elements ρ μ,ν (t) One notes that the memory effects are taken into account in terms of the convolution of the memory kernel f μ,ν (t − τ). A detailed expression for this kernel is given in the Appendix.
To solve the simultaneous integrodifferential components of Eq. (10), we invoke the Laplace transformation f f t e dt { } : , and transform them into a set of algebraic equations. After carefully analyzing the properties of the poles, the conventional residual theorem enables one to accomplish the inverse Laplace transformation and move back from Laplace space into the time-domain.

Second-order time-local equation (TL2).
In the TL2 formalism, the system is considered to be sluggish, hence the bath feedback on the system dynamics can be neglected by approximating t i H t iH exp e xp sys sys sys sys This assumption is reasonable because it is impossible for a system to change its configuration instantaneously. Consequently the system density matrix should be pulled out from the integral to obtain the following QME The detailed expression of the memory kernel h μ,ν (t − τ) is given in the Appendix. It should be emphasized that although ρ μ,ν (t) is pulled out from the integral, Eq. (12) is capable of predicting a non-Markovian dynamics because the time integral of h μ,ν (t) results in time-varying coefficients in front of ρ μ,ν (t). Whether or not such differential equations behave non-Markovianly crucially depends on these time-varying coefficients.
Comparisons with exact results. To illustrate the differences of the approximations explicitly, we apply these two QMEs to a photosynthetic dimer model, which has attracted considerable interest recently 8,9,[44][45][46][47][48][49][50] . We employ the Drude-Lorentz spectral density function (the over-damped Brownian oscillator model) 15,51 , J(ω) = (2λγ/π)[ω/(ω 2 + γ 2 )], which has been widely used for a range of theoretical studies of this type of system [46][47][48][49][50] . We use it here because it is convenient for the comparison with the HEOM. However, in reality, the spectral densities found in real photosynthetic systems tend to be much more complex 46 , and while the HEOM can be extended to model such environments it typically involves a substantial additional numerical overhead 52 . As mentioned in the previous section, within the Drude-Lorentz spectral density the reorganization energy, λ, characterizes the coupling strength to the environment, while the quantity γ determines the width of the spectral density. These two parameters have considerable influence on the dynamics of the system. In Fig. 1, we show the system dynamics given by (a) TC2, (b) HEOM, and (c) TL2 with varying λ and temperature T. The other parameters are fixed at ω 0 = 70 cm −1 , J = 100 cm −1 , and γ = 50 cm −1 (γ −1 = 106 fs). These parameters are typical in photosynthetic systems. The solid curves in each panel denote the populations of the |χ + 〉 state with temperatures T = 300 K (black), 250 K (red), and 200 K (blue), respectively. It can be seen that, at higher temperatures, the population of the |χ + 〉 state transfers to the |χ − 〉 state faster than at lower temperatures, but there is always a crossing so that the thermal equilibrium population of the |χ + 〉 state is larger at higher temperatures.
For small values of λ, the results of the two QME models show excellent agreement with that of the HEOM, indicating that both TC2 and TL2 perform well in the weak system-environment coupling regime and that the bath memory effect is insignificant at small λ. Moreover, the result of TC2 completely coincides with that of the HEOM for very small couplings. We show the comparison between TC2 (solid curve) and HEOM (dot-dashed curve) methods in the inset of Fig. 1(a) for λ = 5 cm −1 , and T = 250 K. This is in line with a recent comparative work in Ref. 53. When λ is increased, the TC2 population results exhibit vigorous beating and produce oscillatory curves up to 800 fs, which is absent in the HEOM result. We attribute these oscillations to the over-estimation of the coherence by TC2. Apart from these beatings, the overall magnitude of the population of HEOM is quantitatively better approximated by TC2 than TL2. The TL2 model yields monotonically-decaying population dynamics that tends to reach thermal equilibrium too rapidly. This leads to a significant over-estimation of the population relaxation rate by TL2, especially at large λ. This over-estimation of the population relaxation rate in Redfield theory has been reported previously 54 , and here we gain further insight into its origin by comparing to the TC2 results.
The dashed curves in Fig. 1 denote the absolute value of the off-diagonal elements of the system density matrix, i.e., the coherence between the |χ + 〉 and |χ − 〉 states. The results from the TC2 method manifestly show the over-estimation of the coherence even if λ is small. When λ is increased, the over-estimation of the coherence becomes quite pronounced. On the other hand, the coherence in the TL2 model decays more rapidly, leading to the sluggish dynamics discussed above. In summary, the coherence dynamics is better approximated by TL2, and the TC2 model may fail in approximating the true coherence for large λ. However, the overall population decay rate predicted by the TC2 is generally more correct than that of TL2. It is interesting to note that the TL2 model yields an exact QME for a pure dephasing spin-boson model (i.e. J = 0) 28 while the TC2 model underestimates the pure dephasing rate, which is in line with our findings here.
Benchmark of approximative QMEs. In the previous section, we analyzed how the coherence terms of the two approximations are qualitatively different from the HEOM exact results. However, those comparisons fail in providing an overall intuitive picture about which model performs better as they are basis-dependent. In other words, it is possible that one model may perform better or worse than another depending on the bases used. In this section, we apply a measure of the non-Markovianity to develop a bases-free benchmark which can quantitatively describe the performance of the approximate methods.
Entanglement and non-Markovianity. Let us consider an isolated ancilla possessing the same degrees of freedom of the system and with which the system forms a maximally entangled initial state j j Fig. 2). If the system evolves according to a process t 0  , , then the Choi-Jamiołkowski isomorphism 39,40 guarantees that the extended density matrix The other parameters are ω 0 = 70 cm −1 , J = 100 cm −1 , and γ = 50 cm −1 (γ −1 = 106 fs). For small λ, both QMEs yield excellent results, as expected. The inset in (a) shows the results given by TC2 (solid curve) and HEOM(dot-dashed curve) for λ = 5 cm −1 , and T = 250 K, illustrating how they almost overlap. However, due to over-estimation of the coherence, the result calculated from the TC2 method shows a slightly higher beating behavior in the population dynamics. In contrast, for large λ the population dynamics predicted by the TC2 method is in better agreement with those of the HEOM, whereas the populations given by the TL2 method are somewhat sluggish and tend to approach thermal equilibrium a bit faster. contains all the necessary information on the dynamics of the system, where anc  is the identity process acting on the ancilla. The entanglement, E(ρ sys,anc ), between the system and the ancilla is a physical quantity which is typically very sensitive to environmental effects.
Another related quantity is the degree of non-Markovianity, NM. Recently, many efforts have been devoted to construct a proper measure of the non-Markovianity 32,33 . Rivas et al. 33 combine the concept of the divisibility of a quantum process 55,56 and the fact that no local completely positive (CP) operation 40 can increase the entanglement E between a system and its corresponding ancilla sys anc sys a nc sys anc Consequently, Rivas et al. 33 proposed that the degree of non-Markovianity within a given time interval [0, t] can be estimated by The non-Markovianity of open-system quantum dynamics can be evaluated at many different theoretical levels [32][33][34][35][36] , and the quantity NM is an extremely strict indicator of non-Markovianity that measures the information exchange in time between the system and its environment. For NM to have a non-zero value, explicit environmental memory effects must be present.
Here we compare the time evolution of the entanglement, E t , and the corresponding degree of non-Markovianity, NM, for the two approximate system-bath models and show how they can provide an integrated picture as to what extent their dynamics deviate from the exact results.
Evaluating non-Markovianity. To analyse the behavior of the non-Markovianity in each method, in this section we will show how the concurrence, a well-known measure for bipartite entanglement 57 , between system and ancilla evolves in time and how the corresponding non-Markovianity [Eq. (15)] depends on the physical parameters of the original spin-boson model.
As an explicit visualization of the integrand in Eq. (15), in Fig. 3, we apply the measure to (a) TC2, (b) HEOM, and (c) TL2 and show the time evolution of the concurrence for different values of λ at temperatures T = 300 K (black), 250 K (red), and 200 K (blue), respectively. The other parameters are ω 0 = 70 cm −1 , J = 100 cm −1 , and γ = 50 cm −1 (γ −1 = 106 fs). It can be seen that, when increasing the temperature and λ, the decoherence becomes more pronounced. Hence, the concurrence will die out earlier for larger λ and higher temperature. As shown in Fig. 3(a), except for λ = 5 cm −1 , which produces monotonically-decreasing concurrence, the TC2 model produces oscillatory curves, in which Environment Figure 2. Schematic illustration of the entanglement measure. We consider a system and a copy of it acting as a well-isolated ancilla possessing the same degrees of freedom of the system. Initially, they form a maximally-entangled state 0 sys anc ρ ( ) = Ψ Ψ , . Then the system starts to feel contact with its environment (denoted by the gray shadow) and evolves according to t 0  , , whereas the ancilla is kept isolated.
a concurrence revival is exhibited around 100 fs and results in a finite degree of non-Markovianity (shown later). A similar entanglement revival can also be seen in biomolecular systems 58 . While in Fig. 3(b,c), HEOM and TL2 produce monotonically-decreasing concurrence and generate no visible non-Markovianity with this measure. In Fig. 4, we show the corresponding measure of the non-Markovianity, NM, calculated using the time evolution of the concurrence shown in Fig. 3(a). Only TC2, for larger λ values, leads to non-zero non-Markovianity, while TC2 at λ = 5 cm −1 , HEOM, and TL2 generate null results due to the monotonically-decreasing concurrence. This comparison not only shows that the TL2 yields a better approximation to the HEOM dynamics, but also explicitly demonstrates the degree to which TC2 deviates from HEOM. We again attribute this deviation to the over-estimation of coherence shown in Fig. 1. In addition, it can be seen in Fig. 4 that NM tends to increase with increasing λ and temperature. We will investigate this below in a regime where the HEOM results exhibit similar behavior.
Increase of non-Markovianity with λ and temperature. The other two important parameters in our spin-boson model are the level spacing ω 0 and the bath relaxation time γ. The former affects to what extent the state |χ + 〉 is delocalized, while the latter is related to the correlation time of the environment and is directly connected to the non-Markovianity of the system. In general, the concurrence will die out faster for larger λ and higher temperatures. The coherence over-estimation of the TC2 method is manifested by a concurrence revival around 100 fs for larger values of λ, whereas HEOM and TL2 produce a monotonically-decreasing concurrence.
Scientific RepoRts | 5:12753 | DOi: 10.1038/srep12753 In Fig. 5(a), we reduce ω 0 to 40 cm −1 and fix the other parameters at λ = 5 cm −1 , γ = 50 cm −1 , and T = 200 K. The reduction of ω 0 leads to a manifest concurrence revival around 100 fs in the TC2 concurrence dynamics, a result of stronger delocalization and significant enhancement of the coherence effect. An analogous result can be seen in Ref. 37. In the mean time, the concurrence of the HEOM result is still monotonically decreasing. The TC2 model further over-estimates this enhancement and ends up with finite non-Markovianity within all range of temperatures shown in Fig. 5(b). The TL2 model predicts almost-Markovian results, besides the very small non-Markovianity at low temperatures, again showing a better agreement with the HEOM exact results.
In Fig. 6(a), γ is further reduced to 20 cm −1 (γ −1 = 265 fs) to investigate the effect of slow environments. As the spectral density function is narrower, the correlation time of the environment becomes long compared with the characteristic time of the system dynamics. Hence the information on the system dynamics is more likely to be retained in the environment and flow back into the system. This back-flow of information in turn affects the behavior of the system and results in beating in the concurrence curves for all methods. As shown in Fig. 6(b), the TC2 model predicts a non-Markovianity much larger than the exact results. On the other hand, the TL2 model predicts a non-Markovianity in excellent agreement with the HEOM results, with only a small under-estimation of the non-Markovianity in this set of parameters.
The above comparisons exhibit an interesting tendency for NM to increase with λ and temperature. Several relevant theoretical and experimental works have reported 41-43 that strong system-environment correlations are helpful for maintaining quantum coherence even at high temperatures. As a result, higher temperature may in turn activate more phonon modes in the environment without destroying the quantum coherence significantly. This provides more channels via which the system can interact with the environment. In the language of quantum information science, smaller γ and strong system-environment correlation may help to preserve the dynamical information; while larger λ and higher temperature may increase the possibility that this information can flow back into the system from the environment. Consequently, this increase of NM with larger temperature and λ is a result of the competition between the back-flow of information and thermal fluctuations. Meanwhile, the magnitude of the concurrence is reduced by the stronger random fluctuations in the environment.

Discussion
In summary, we first investigate the dynamics of two perturbative second-order QME methods, TC2 and TL2, and compare their results with the numerically-exact results calculated by HEOM. We find that TC2 can approximate the HEOM population better than TL2. However, a drawback of the TC2 model is its over-estimation of the coherence. This drawback results in the TC2 model predicting too much beating behavior in the population dynamics and limits the accuracy of TC2. In constrast, the TL2 model predicts sluggish dynamics and loss of coherence faster than that of the exact HEOM. As a result, the population tends to reach thermal equilibrium too rapidly.
To further investigate the dynamics and establish a benchmark for the performance of perturbative QMEs, we combine the concept of Choi-Jamiokowski isomorphism 39,40 , entanglement with an ancilla 57 , and a measure of non-Markovianity 33 to provide a quantitative way to determine how much the coherence dynamics and memory effects are deviating from the exact result. This provides a deep physical insight on the effects of each parameter and a single quantity to determine how much the QME dynamics deviates from the exact results. Here we find that the non-Markovian measure indicates that the TL2 approximates HEOM better than TC2 in terms of the coherence dynamics and memory effects for the The parameters are the same as used in Fig. 3. Among the three methods investigated in this work, only TC2 at higher λ generates non-zero non-Markovianity for these parameters. At λ = 5 cm −1 TC2 correctly produces the expected Markovian dynamics for this regime.
dimer system studied here. In addition, while it is well understood that the reorganization energy λ and temperature enhance the effect of thermal fluctuations in the environment on the system, increasing these parameters can have surprising results. In particular, our results show that higher temperature increases information back-flow from the environment, thus increasing the non-Markovianity of the system dynamics, even though the concurrence itself undergoes faster decay. Note that photosynthetic systems and other molecular light-harvesting networks are in general far more complex than the models studied here 4,5 , and more general models should be considered for realistic systems [59][60][61][62] . Nevertheless, the focus of this work is the physics revealed in the comparison of the theoretical methods and the application of the non-Markovianity measure for revealing new physical insights. The theoretical methods examined here have often been applied to model real photosynthetic systems, and the quantitative measures we employ are themselves model independent. The non-Markovianity analysis proposed here could be easily used to investigate coherence dynamics in more complex systems and more general models. Therefore, these results could have important implications in the theoretical modeling of electronic coherence in photosynthetic systems 8,9,47 .

Methods
Full expressions for the TC2 and TL2 quantum master equations. The detailed expression of the TC2 integrodifferential QMEs Eq. (10) is given by