Non-Markovian Second-Order Quantum Master Equation and Its Markovian Limit: Electronic Energy Transfer in Model Photosynthetic Systems
A direct numerical algorithm for solving the time-nonlocal non-Markovian master equation in the second Born approximation is introduced and the range of utility of this approximation, and of the Markov approximation, is analyzed for the traditional dimer system that models excitation energy transfer in photosynthesis. Specifically, the coupled integro-differential equations for the reduced density matrix are solved by an efficient auxiliary function method in both the energy and site representations. In addition to giving exact results to this order, the approach allows us to computationally assess the range of the reorganization energy and decay rates of the phonon auto-correlation function for which the Markovian Redfield theory and the second order approximation is valid. For example, the use of Redfield theory for in systems like Fenna-Mathews-Olson (FMO) type systems is shown to be in error. In addition, analytic inequalities are obtained for the regime of validity of the Markov approximation in cases of weak and strong resonance coupling, allowing for a quick determination of the utility of the Markovian dynamics in parameter regions. Finally, results for the evolution of states in a dimer system, with and without initial coherence, are compared in order to assess the role of initial coherences.
The quantum mechanics of open systems, i.e. systems interacting with an external environment, is currently the focus of widespread attention. Of particular recent interest is the issue of the extent to which quantum coherence of the system is maintained in the presence of the environment. Two significant modern examples may be noted: (a) the need to maintain coherence in order to implement methods for quantum mechanically controlling molecular processes, and (b) issues of the role of such quantum coherent processes in natural environments, such as the observed long-lived coherent electronic energy transfer (EET) in photosynthesis.
In either of these cases, and in many other examples as well, dynamical evolution of the open system provide a significant computational challenge. As such, a variety of approximations are often invoked to propagate the system, such as the second Born approximation to master equations and the Markov approximation, both of which are the focus of this paper. In particular, in this paper we introduce a simple method to solve the second Born quantum master equation without doing the Markov approximation on the slowly decaying envelope of the density matrix. In addition to being straightforward, this approach allows, by comparing to results using the Markov approximation, a reliable determination of the range of coupling strengths and decay rates of the bath auto-correlation function, over which one can use the Markovian theory and the second order approximation. In addition, by examining the size of the fourth order term, this approach affords an estimate of the range of validity of the widely used second-order approximation for model photosynthetic systems.
Although the method developed here is applicable to general systems with exponential bath correlation functions, for computational simplicity we study, as do others, the dimer system, generally regarded as a simple photosynthesis EET model. In this case our Markovian analysis contrasts with, for example, that in Ref.  in which Markovian Redfield theory is used apriori and its consequences analyzed, as opposed to comparison with exact results.
The particular challenge arising from these EET types of systems relates to the parameter range in which the dynamics occurs. Specifically, quantum dynamics can be readily analyzed in two limiting cases defined by the relative contributions of the inter- system coupling responsible for excitation energy transfer (EET), and system-bath coupling constant , responsible for decoherence. These parameters define two important time scales: the excitation transfer time scale , and the decoherence time scale . If the system-bath coupling is very weak and , the system is almost closed and the Schrödinger equation can be used to study the dynamics. In the opposite case (strong system-bath coupling), the system is open, the decoherence rate is very fast, the dynamics is almost incoherent and a simple Pauli type master equation description suffices. These limiting regimes are well understood. Many real systems, such as a number of harvesting systems , however, fall between these extremes. Recent observations  of the long-lived EET has reactivated interest in these systems.
The standard approach used to treat this intermediate regime is to use the second Born quantum master equation , a perturbative master equation up to second order in system-bath interaction with weak system-bath coupling, plus its Markovian approximation (e.g., as in the Redfield master equation). Recently, two approaches have been studied for arbitrary coupling regimes. One is based on weakening the system-bath coupling removal of system-bath interaction and repartitioning the Hamiltonian term using a polaron transformation, followed by the standard second Born master equation . The second approach is based on a reduced hierarchy equation of Kubo and Tanimura, starting from the path integral approach for quantum dissipative systems . Additional methods are also being developed by the community. Here, as noted above, we introduce and utilize a particularly direct approach.
In the following section (Section II) we outline the basic model for a dimer. In Section III we introduce the second Born quantum master equation and phonon correlation function, diagonalize the Hamiltonian and cast the master equation into both the site and energy representations. Section IV gives a new auxiliary function method of solving these equations, and provides an analysis of results in the Markovian approximation. Discussion of the results and the underlying physical picture is given at the end of Section V. In Section VI we discuss the regime of validity of the second order approximation in the master equation by estimating the order of magnitude of the fourth order term.
The vast majority of treatments in the literature utilize coherence-free initial conditions and study the subsequent dynamics. Results for these initial conditions are compared to that obtained with model excitation with weak light in Section VII. The last section provides a brief conclusion.
2The Model: Dimer System
Consider a model dimer system given by the following standard Frenkel exciton Hamiltonian :
Here represents the state in which only the site is excited and all others are in the ground state. The quantity is the excited electronic energy of the site in the absence of phonons, and is the electronic coupling between the sites which is responsible for EET. The ground state energies of the donor and acceptor are set equal to zero and is the reorganization energy of the site that is dissipated in the bath after the electronic transition occurs. The quantity is the dimensionless displacement of the equilibrium configuration of the phonon mode between the ground and the excited electronic state of the site, and are the dimensionless coordinates and momenta of the phonon mode of frequency .
3The second-Born quantum master equation
The method of projection operators used to obtain open system master equations is well known . With the help of projection operators one can obtain the following quantum master equation for the reduced density matrix of the system in the second Born approximation, which is valid when system-bath coupling is weak as compared to the characteristic energy scale of the system [see, e.g., Ref. ].
Here, the interaction representation has been used, which is defined for system operators as,
where is the time evolution operator, and is the system Hamiltonian. Here, the bath is assumed to be a continuum of harmonic oscillators, and the bath correlation functions are defined as
Below, the canonical average of the bath operators, , which involve the averaging over the product of displacement and bath position co-ordinates is taken to be zero. The above master equation [Eq. (6)] is also termed the time convolution equation and can be obtained from the Nakajima-Zwanzig equation with a zeroth order approximation to the time evolution operator in the kernel .
Converting this master equation [Eq. (6)] back to the Schrödinger representation gives
We consider the case where the characteristics of the bath as seen by both the sites are the same, and there is no bath correlation between the sites. The bath correlation function is then of the form , where
where is the Bose-Einstein distribution function. Assuming the Drude-Lorentz model for the spectral density where is the reorganization energy, and assuming the high temperature approximation (), as is appropriate for the systems like the FMO model [see Refs.  and ], we obtain the correlation function as,
3.1Explicit Site Representation of the Non-Markovian Master Equation
For the explicit site representation we need the eigensystem of the Hamiltonian. Assuming the eigenvalues and eigenvectors for the system Hamiltonian
can be easily obtained as
Here the column vectors denote components in the site basis, the eigenkets are normalized and, since , they are orthogonal. With lengthy but straightforward calculations, Eq. (Equation 2) for the reduced density operator can be written explicitly in the site representation, using Eqs. (Equation 3) and (Equation 4), and as a set of coupled integro-differential delay equations,
with . Here (site), , and , with subscripts denoting the sites.
3.2Energy Representation of the Non-Markovian Master Equation
The kets in Eq. (Equation 7) are the eigenstates of the Hamiltonian . The equation for a general element of the reduced density matrix in energy representation
Results in the energy representation, using Eq. (Equation 5) are discussed below.
4Method of solution: Non-Markovian
To obtain a solution for the non-Markovian case, we first convert the coupled integro-differential equations in the site representation [Eq. (Equation 8)] to a larger number of coupled ordinary differential equations, a transformation made possible by the exponential form of the correlation function. The resultant coupled ordinary differential equations can be numerically solved easily. This transformation is performed as follows. First, for computational simplicity we put in Eq. (Equation 8) and then in the resulting equations, and define three auxiliary functions :
Here, and we also define . We then obtain six coupled ordinary differential equations, three from Eq. (Equation 8) and three from differentiating the three auxiliary functions, giving:
where overdots denote derivatives with respect to . These equations can be efficiently solved numerically.
For comparison with other studies, results are given below for the particular initial conditions: . These initial conditions (corresponding to all the population being on site 1, and no coherences), are those which have been used extensively in previous investigations [see Ref.] but are somewhat unphysical, because they lack initial coherences which become important in photo-excitation. We treat this problem of initial conditions and state preparation with a more plausible model in Section VII.
5Energy representation and Markovian limit
To consider the Markov approximation, we note that it is particularly simple to invoke in the energy representation. Hence, below we first utilize the energy basis and then convert the result back to the site representation for comparison with the non-Markovian solution.
The Markov approximation can be performed when the time scale on which the envelope of the density matrix decays is much longer than the decay time of the phonon correlation function . Then one can introduce the following approximation:
As discussed in Ref. , the non-Markovian regime is marked by slow dissipation of the reorganization energy (i.e., the slow decay of the phonon correlation function as compared to relaxation dynamics time scale, the decay of the envelope part of the density matrix). Transitions occur in accord with the vertical Franck-Condon principle. In the Markovian regime phonon relaxation is very fast (e.g., large ) as compared to the decay of the envelope of the density matrix. Thus, phonons remain effectively in equilibrium during the EET process in the Markovian regime .
To obtain the equations in the Markov approximation, Eq. (17) is first converted to dimensionless form with and . Putting in the resulting equation in the energy representation and then implementing the above approximation on the density matrix elements allows the time integration to be performed easily for the case of exponential phonon correlation function [Eq. (11)]. The result is the set of Markovian equations;
Here . Equation (Equation 13) constitutes a system of coupled ordinary differential equations that can be solved with given initial conditions.
The results can then be transformed back to the site representation using the transformation
where is in site representation, is in energy representation, and . Equation (Equation 14) constitutes four linear equations that provides the relationship between the representations.
5.2Limiting Cases: Analytical Results
Strong Coupling Case:
for the traditional initial conditions . The Markov approximation can be performed when the time scale on which the envelope of the density matrix decays is much longer than the decay time of the phonon auto-correlation function. Hence, must hold for the Markov approximation to be valid in the domain. (Note that the decay time constant for the phonon auto-correlation function is unity, since we defined .)
Weak Coupling Case:
The reduction from a large number (33) of terms in Eq. (Equation 13) to the smaller number of terms [in Eq. (Equation 16)] is made possible by neglecting small terms of the order of . However, even in this approximation the above equations do not admit a simple analytic solution, and no simple analytic expression can be given for the range of validity of the Markov approximation. Hence, we invoke a further approximation, neglecting second order terms as compared to first order . By separating real and imaginary parts as and writing , Eqs. (Equation 15) and (Equation 16) become
These coupled ordinary differential equations have the straightforward solution:
with initial conditions . Here and . Interestingly, for the density matrix elements do not change with time. This is due to our approximations of only retaining terms first order in . Thus, for the condition of Markov approximation to hold requires (again noting that the decay time constant for the phonon auto-correlation function is unity since ).
These analytic results are summarized in Table I, and these inequalities have been numerically verified (e.g., see Fig. ?). Note that the result goes over to the result as gets smaller.
In other parameter regimes, the validity of the Markov approximation [Eq. (Equation 13)] must be determined by numerical comparisons with the exact result [Eq. (Equation 11)]. Figure ? compares the solution for the Markovian master equation to the non-Markovian results for the standard electronic coupling parameter values in photosynthetic EET: K , a regime in which the estimates in Table I do not apply. The initial excitation is assumed to be on site one. The Markovian approximation is seen to be very good for , fair for and invalid for reorganization energies .
To explore regimes of validity of the Markovian approximation for other values of the physical constants, we present sample results in Figs. ? and ?, obtained by varying and , keeping cm. The results show for these cases that the Markovian approximation is poor for large and “small’ and for large and small . Other parameter values can be readily examined computationally using this approach.
6Regime of validity of the second-order approximation
We have considered the master equation up to the second order in system-bath interaction. The aim of this section is to determine the system-bath interaction energy range (as represented by the reorganization energy ) over which the second order master equation [or coupled system of equations (Eq. (20)] can be used. To do so we compare estimates of the second order and the fourth order terms. This can be done analytically for the parameter regime where . To do so we note that up to the second order, with the master equation written in dimensionless time form [Eq. (20)], the magnitude of the second order term is of the order of for the case . This arises by noting that , and the integral , for the standard set of parameters (K).
Similarly, we can estimate the parameter dependence of the fourth order term. To estimate this we recall the Nakajima-Zwanzig master equation (valid to all orders)
where the time evolution operator is
Here, is the well know projection operator and is the system-bath Liouvillian (). The time ordering operator orders time dependent operators from left to right with decreasing time arguments, to take into account the non-commutation of operators at different times.
The zeroth order approximation to the time evolution operator gives the second order quantum master equation, the first order approximation to the time evolution operator gives the third order contribution which vanishes as the bath average of odd bath operators vanish (see Ref. ), and the second order approximation to gives the fourth order contribution. In order to estimate the magnitude of the latter term we write the fourth order term A from Nakajima-Zwanzig equation as
We start from the interior commutator and recall that and (where the summation convection is used). The bath average of single bath operators vanish  so that times the interior commutator gives . Similarly, writing for the second from the interior commutator, and simplifying, the bath trace operation gives us the two time bath correlation functions . Repeating the same operations for the remaining commutators, noting that the bath averaging for the odd bath operators vanish and using the Wick theorem , we obtain that the fourth order term includes the product of two time bath correlation functions i.e., . On converting the equation into dimensionless time form as described in Section IV, and using the high temperature approximation for the correlation function, we conclude that the order-of-magnitude of the fourth order term is;
The ratio of the fourth-order term to the second-order term is therefore , which is of the order of for , and for . This suggests that the second order approximation for the master equation is good for for the standard set of parameters (K). However, for large the fourth order term cannot be neglected. Interestingly, the domain of applicability of the second-order approximation in this regime has a dependence on the same collection of parameters, small), as does the Markov approximation in the .
7Initial state preparation by an ultra-short laser pulse
To understand the effect of initial coherences on the subsequent quantum dynamics, we consider a model of two 1-D harmonic oscillators separated by a distance that are excited by an ultra short laser pulse (Figure 1). The results of this excitation are used below as sample coherent initial conditions for dimer propagation. Note that we restrict attention, as do most treatments of excitation of light-harvesting systems, to the “one-exciton manifold”, i.e. single excitations on each site. This model is useful to examine the dynamics, but should be augmented by excitation of states with bi-excitons of examining issues like entanglement, where the contributions from higher exciton states, no matter how small in magnitude, affect the entanglement measure.
In the system (Figure 1) the wave vector of the laser pulse along the propagation direction makes an angle with the line perpendicular that joining the oscillators. The laser field is treated semiclassically with the field sufficiently weak to allow first order perturbation theory for the light-oscillator interaction . The total Hamiltonian in the coordinate representation is then
with the vector potential . For a coherent laser pulse, . If , where is the speed of light, then changing the summation to integration assuming continuous distribution of modes, we have
For a Gaussian pulse