Quantum effects in biology: master equation studies of exciton motion in photosynthetic systems
Abstract
The present review is devoted to our recent studies on the excitonic motion in photosynthetic systems. In photosynthesis, the light photon is absorbed to create an exciton in the antenna complex of the photosynthetic pigments. This exciton then migrates along the chainbiomolecules, like FMO complex, to the reaction centre where it initiates the chemical reactions leading to biomass generation. Recently, it has been experimentally observed that the exciton motion is highly quantum mechanical in nature i.e., it involve long time ( femto sec) quantum coherence effects. Traditional semiclassical theories like Forrester’s and second Born master equations cannot be applied. We point out why the 2nd Born nonMarkovian master equation and its Markovian limit (also called the Redfield master equation) cannot be used to explain the observed long coherences. Briefly, the reason is that these approaches are perturbative in nature and in real light harvesting systems various couplings (systemsystem and systembath) are of the similar order of magnitude. Various new approaches are being developed to go beyond the above two limiting theories. The present review is not a review in the usual sense of the word as we summarize our own approaches and only refer to the literature for the other ones. A brief introduction to the sophisticated 2D photon echo spectroscopy is also given at the end with an emphasis on the underlying physics of the multidimensional echo spectroscopies.
Keywords: Excitation energy transfer in photosynthesis; Quantum master equations; stochastic theories; line broadening; 2D photon echo spectroscopy; Nonequilibrium statistical mechanics
Contents:
I Introduction to Excitation Energy Transfer (EET) and the statement of the Problem
Photosynthesis provides chemical energy for almost all life on earth. Understanding of the natural photosynthesis can enable us to construct artificial photosynthesis devices and a solution to the future energy problems. We know that the safe (green) energy resources is a big challenge of the future and we know that our present energy technology is seriously disturbing our environment(1). This prompts us to investigate “ Green” resources of energy and photosynthesis is one of them. Photosynthesis is an interesting phenomenon. On the average, on earth, biomass worth twice the mass of the Great Pyramid of Giza ( million tons) is being produced every hour or million kg/sec.
In photosynthetic systems a central role is played by the energy transport ”wire”–the FMO protein which is a trimer made of identical subunits containing seven bacteriochlorophyll (BChl) molecules each (Fig. 1). The photosynthesis process can be divided into the following steps(2):

A light photon is absorbed to create an excitation in the antenna complex of the photosynthetic pigments.

This excitation then migrates along the chainbiomolecules, like FMO complex, to the reaction center where it initiates the chemical reactions leading to biomass generation.
Quantum dynamics of EET (excitation energy transfer) can be very easily analyzed in the following two limiting cases. We identify first the couplings: Two important couplings:

Inter BChl molecules (systemsystem) coupling (which is responsible for EET).

The systembath (BChl moleculeprotein) coupling (which is responsible for decoherence).
Two important timescales:

Excitation transfer time scale (usual in photosynthesis pigments).

Decoherence time scale .
If the systembath coupling is very weak and , the system is almost closed and dynamics is quantum mechanical in nature, i.e., one can use the Shroedinger’s equation to analyze it in the extreme case. But in the opposite case (strong systembath coupling), the system is almost open, decoherence rate is very fast and the dynamics is almost incoherent. One can analyze the process with simple Pauli type master equation with rate of transfer of the excitation from one molecule to other calculated with Fermi’s golden rule (Forester’s theory)(3).
The important problem arises in the intermediate regime as real light harvesting systems do not fall in either of the extreme cases. Recently, it has been experimentally observed that the exciton motion is highly quantum mechanical in nature i.e., it involve long time ( fs) quantum coherence effects(4). These discoveries caused a lot of excitement in field(5) as the traditional view was of incoherent ”hopping” of the exciton(Fig. 2). In the weak systembath coupling case, the standard approach was the 2nd Born quantum master equation which is a perturbative quantum master equation (upto second order in systembath interaction). Its Markovian and secular approximation is known as Redfield master equation(6). The 2nd Born quantum master equation can be obtained from the NakajimaZwanzig projection operator technique by restricting the perturbation series upto second order(6). These simple and powerful projection operator techniques were introduced by Zawnzig in 1960’s in then active field of nonequilibrium statistical mechanics. As the the light harvesting pigments fall in the intermediate regime (systembath coupling is of the order of systemsystem coupling) one clearly cannot use the 2nd order perturbative quantum master equation for its study in its original form. But in its modified form its scope becomes wider(7).
It is of interest to quantitatively know upto what value of systembath coupling strength and other important couplings in the problem, one can use the 2nd Born master equation. Section II deals with this.
In an important case of fast bath relaxation (when bath degreesoffreedom reorganize very fast as compared to the transfer time scale of the exciton) a very useful approximation can be made. The details of which are given below. This is called the Markovian approximation. In the following sections we summarize our study(8) of the quantitative determination of the regime of validity of the second order approximation and the Markovian approximation.
Ii Microscopic approach: 2nd Born master equation
As is well known that the 2nd Born quantum master equation can be obtained from the NakajimaZwanzig projection operator technique by restricting the perturbation series upto second order in the systembath interaction(6). In the following we will apply this master equation to a concrete model of a dimer (open twostate quantum system) which caricature the dynamics of decoherence in a typical photosynthetic system(9).
ii.1 NonMarkovian solution
Projection superoperators and dynamics of the relevant system: The total system (relevant system (electronic part) + Bath (phonons)) dynamics is pure quantum in nature. The partial time derivative of the total density matrix is given by Liouvillevan Neumann equation (classical equivalent is the invariance of the ”extension” in phase space):
(1) 
Here the interaction representation is used (see for details any standard refs(6)). is the systembath interaction Hamiltonian (for the 2nd Born approximation . is the bath Hamiltonian. To construct the equationofmotion for the relevant part of i.e., , one defines the superoperator:
(2) 
called the projection superoperator, one also defines . By applying on the Liouvillevan Neumann equation turnbyturn, we get:
(3) 
Solving the second equation formally for the irrelevant part , and inserting in the first, one obtains the required equationofmotion for the relevant part (see for details(6)):
With 2nd Born approximation (i.e., by expanding the time evolution operator upto the first power in the systembath interaction(6)) and for a traditional dimer system(9):
(4) 
master equation takes the form,
(5) 
Here in the Hamiltonian, represents the state in which ONLY th site is excited and all others are in the ground state i.e., . is the system Hamiltonian which consists of the electronic Hamiltonian for the two level system, and the Hamiltonian for the reorganization energy (the elastic energy related to the physical organization of the bath degreesoffreedom). is the phonon Hamiltonian and is the systembath coupling Hamiltonian. In the absence of phonons, is the excited electronic energy of site and is the electronic coupling between the sites which is responsible for excitation transfer. The ground state energies of both donor and acceptor are set equal to zero and is the reorganization energy of the site (Dissipated energy in the bath after the electronic transition occurs). are the dimensionless displacement of the equilibrium configuration of the phonon mode, dimensionless coordinates, momenta of the phonon mode respectively.
In the master equation, the bath correlation functions (bath is assumed to be a continuum of harmonic oscillators (valid when anharmonic terms are not important)) are:
(6) 
We consider a case where the characteristics of the bath as seen by both the sites are the same, and there is no systematic bath correlations between the sites. Thus the bath correlation function takes the form: :
(7) 
(8) 
Assuming the DrudeLorentz model(9) for the bath spectral density, and assuming the high temperature approximation (), as appropriate for the FMO problem, we obtain
(9) 
ii.2 Representations
The equation (5) is an operator equation and this can be expressed in site or in energy representation. In site representation, with definitions (site), , and , and with lengthy but straightforward calculations (see for details(8)), the equation (5) can be written explicitly as a set of coupled integrodifferential delay equations:
with .
For energy representation we need the eigensystem of the Hamiltonian. Let the kets be the eigenstates of the Hamiltonian . The reduced density matrix in energy representation can be expressed as:
(10) 
with time evolution given as,
(11)  
with and
(12) 
Eigensystem of the Hamiltonian: Assuming the eigenvalues and eigenvectors of the system Hamiltonian
can be easily obtained as
Here the column vectors denote components in the site basis. The eigenkets are normalized and are orthogonal .
ii.3 Numerical Approach
It is well known that the numerical propagation of integrodifferential equations is an involved task and time consuming. In the following we construct a simple method to the solution of numerical integration(8). Specifically we utilize the exponential nature of the bath correlation function which helps to convert the set of coupled integrodifferential equations to a bigger set of ordinary differential equations:
(13) 
With . Here, and we also define .
We obtain a set of coupled ordinary differential equations (note that these are much simpler to solve as compared to coupled integrodifferential equations):
(14) 
Our aim is to use this to establish the parameter range over which the Markovian approximation is valid. Before doing so we compare this method with the traditional method of solution of integrodifferential equations(10). In the straightforward numerical method (traditional method) the integrodifferential equation is first written as integral equation with double integration as converted to . The double integration is then done selfconsistently with numerical integration(10).
Speed: On the Lenovo ThinkCentrei7, the traditional method took minutes to obtain the result, whereas the present method took only about a few milliseconds. A sample comparison is given in Fig. (3).
ii.4 Markovian limit
As mentioned before in an important case of fast bath relaxation (when bath degreesoffreedom reorganize very fast as compared to the transfer time scale of the exciton) a very useful approximation can be made. This is called the Markovian approximation(21). To do this approximation, note that it is particularly simple to invoke it in the energy representation (as one can “average out” the fast oscillations in the system’s density matrix and can compare the temporal envelope of the system’s density matrix with time decay of the bath correlation function). Hence, below we first utilize the energy basis and then convert the result back to the site representation.
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 (6). One can then introduce the following approximation:
To do this energy representation equations are first converted to dimensionless form with (see for details(8)). Putting in the resulting equations 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. The result is the set of Markovian equations:
Here . The results can then be transformed back to the site representation using .
ii.5 Limiting Cases: Analytical Results
Markovian and nonMarkovian results were obtained computationally and compared for various regimes. Before presenting them we show some interesting analytic results in two extreme cases on the parameter dependence of the region of validity of the Markov approximation:
Strong Coupling Case:
For , we have and for and for and for all . One can then analytically solve the coupled equations to obtain the simple expression
(15) 
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 autocorrelation function. Hence, must hold for the Markov approximation to be valid in the domain.
Weak Coupling Case:
For this case, we have . Hence, in this domain , and . This leads to . , and .
where . See for details the second paper in(8). By separating real and imaginary parts as and writing , we have:
(16) 
with initial conditions . Here and . Thus, for the Markov approximation to hold requires . These are summarized in Table 1. A numerical verification of these analytical inequalities is given in Fig. (4).
Case  Markovian approximation 

ii.6 Computational Results
To investigated the validity regime of Markovian approximation beyond the above two limiting cases we have to rely on the numerical computation (as the analytic approach becomes very cumbersome). We here display an extensive list of graphs (Figs. 5,6, and 7) which explore how Markovian approximation behaves for various values of parameters.
First consider the nonMarkovian results (solid curves). Coherent oscillatory dynamics up to substantial time scales are evident. Oscillatory behavior in the populations are seen (Fig. 5) to be accompanied by oscillations in the off diagonal elements representing coherences. The oscillations fall off faster with increasing reorganization energies, as expected. Another point to be noted is that the relaxation to equilibrium populations occurs on a longer time scale than does relaxation of the coherences to zero. This difference is more evident at smaller values of reorganization energy . The dependence of population relaxation on the bath correlation decay time is shown in Fig. 7. Clearly, the larger the (i.e., fast bath relaxation), the better is the Markov approximation.
Figures 5  7 also contain a comparison of the Markovian limit to the nonMarkovian solution for domains other than those in Table I. For the standard electronic coupling parameter values in photosynthetic EET: K, figure 5 shows that the Markovian approx is very good for , fair for , and invalid for reorganization energies . As typical values of in photosynthetic EET systems are considerably larger than cm, these results support the conclusions of Ref. (9), although with a different approach.
Figs. 6 and 7 show the validity of the Markov approximation obtained by varying and respectively, but keeping cm. The results show that the Markovian approximation is poor for large and small , and for large and small . Other parameter values can be easily examined using this approach.
From these qualitative conclusions a rough physical picture can be drawn as depicted in Fig. (8). In the Markovian regime, after photoexcitation, bath (nuclear coordinates) reorganize very fast representing an ”apt” bath (right hand picture) while in the nonMarkovian regime bath correlations exist for longer time scale (of the order of exciton transfer time scale). Thus the combined (system+bath) dynamics is much more complicated in the nonMarkovian regime.
Here, results are given 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 literature(9) but are somewhat unphysical(11), because they lack initial coherences which become important in preparatory photoexcitation. We have considered this problem in the second reference of the list(8) which considers photoexcitation of a dimer oscillator system with an ultrashort laser pulse. It appeared that the presence of initial coherence ( at ) effected the time scale on which the populations reach equilibrium value but had little effect on the overall dampingout of the coherences (see for details(8)).
Iii Phenomenological approach: a stochastic model
In the previous sections we saw that 2nd Born quantum master equation cannot be applied to the real light harvesting systems as in these systems the systembath coupling () is of the same order of magnitude of the systemsystem coupling (). To make the situation more intractable, it is not possible to justify the Markovian approximation when (given that Markovian master equations are much easier to solve than the nonMarkovian ones). Thus the use of Markovian Redfield theory to these systems is questionable as pointed out in(9); (8). This open up a difficult problem. One should formulate some nonperturbative theories. Recently Ishizaki and Fleming(9) have developed a formalism which goes beyond the limitations of the 2nd Born master equation. They use the reduced hierarchy equation approach previously developed by Tanimura and Kubo(12). There is an other route to the problem pioneered by people like Silbey(13). In this approach one uses a unitary transformation (called polaron transformation) to completely eliminate the systembath coupling Hamiltonian. But this renormalize the system Hamiltonian. Then one repartition the resulting system Hamiltonian to identify a weaker term which can be used as a perturbation. The remaining problem is done in line with 2nd Born master equation(14).
We have developed an alternative stochastic approach which is phenomenological in nature(15). This approach, as its input, takes the homogeneous line width from the experiment and uses Kubo’s stochastic theory of motional narrowing to get phenomenological decoherence rate.
iii.1 The model and its solution
We again consider the dynamics of exciton transfer between two molecules modeled as twolevel electronic systems (Fig. 9). These twolevel systems are electronically coupled with each other with coupling .
Due to the electronic coupling between the molecules the exciton will transfer backandforth between the molecules. This will happen forever if the molecules are completely isolated—pure oscillatory quantum motion. Now consider that our twomolecular system is open i.e., interacting with the external bath degreesoffreedom—the phonons. It is well known that the dynamics remain quantum at short time scales and becomes classical at longer time scales(16). This quantumtoclassical crossover happens at a critical time scale which is inversely proportional to systembath coupling energy (). Large (systembath coupling) means fast quantumtoclassical crossover and vice versa. We denote systembath coupling strength with in the subsequent subsections (before we used ).
In a recent contribution(15) we extracted from experimental information using Kubo’s stochastic theory of line shapes and observed upto what timescale one could see quantum effects.
The important point is that we modeled the dynamical effect of phonons as a stochastic noise. The total Hamiltonian takes the form
(17) 
Here is the energy of the upper electronic level of the first molecule and is that of the second molecule. The ground state energies of both the molecules are taken to be zero. The energy separation has a random component (due to phonons) which we denote with . is a stochastic process taken here as Gaussian White Noise (GWN):
(18) 
Here, as mentioned before, is the strength of systembath coupling (also known as dynamical disorder) measured in the units of frequency.
We start with LiouvillevonNeumann equation for the total density matrix,
(19) 
As is a stochastic operator, thus is also a stochastic operator. Therefore, we need to evaluate the averaged density matrix. So we need to do an averaging over the dynamical disorder which is denoted by . We define .
Averaging over dynamical disorder:
(20) 
In term I and II above we have terms like . As is a functional of (a stochastic quantity) the will also be a stochastic function. To decouple these we will use the famous theorem of Novikov(17):
(21) 
Here is the functional derivative. Using the properties of stochastic noise and with some simplification(see appendix A), we get
(22) 
Finally, one has a set of coupled differential equations:
(23) 
The above system of ODEs can be solved analytically, however, the exact expression is very cumbersome. We give the analytic solution only in the long time limit (see Appendix B).
To simulate the dynamics of decoherence in this dimer model with physical parameters of the FMO problem (), we need to find out our phenomenological parameter , for this we use Kubo’s stochastic theory of lineshapes.
iii.2 Kubo’s stochastic theory of lineshapes and estimation of from motionally narrowed lineshape
We now determine the phenomenological parameter . Let us focus on exciton 1 in 2D photon echo spectra which occurs at (Fig. 10).
We use Kubo’s randomly modulated oscillator model for the exciton under question. The energy levels of the model system are given in Fig. 10. The levels 1 and 2 are separated on the average by .
Let the level 2 be randomly modulated with a random process such that . The random frequency of the level 2 is then given as . We use the classical oscillator model for the molecule with resonance frequency . The equation of motion is , with solution
(24) 
The absorption spectrum measures a large number of possible realizations of the random process ( The radiation field acts on a macroscopic number of molecules present with different realizations of in the sample). Thus one need to consider an ensemble average to compare with absorption spectrum: