Quantum effects in biology: master equation studies of exciton motion in photosynthetic systems

Quantum effects in biology: master equation studies of exciton motion in photosynthetic systems


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 chain-biomolecules, 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 non-Markovian 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 (system-system and system-bath) 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.

Key-words: Excitation energy transfer in photosynthesis; Quantum master equations; stochastic theories; line broadening; 2-D photon echo spectroscopy; Non-equilibrium statistical mechanics

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.

Figure 1: On the average, on earth, biomass worth twice the mass of the Great Pyramid of Giza ( million tons) is being produced every hour. (Right) FMO protein acts as a “wire” which connects an antenna to a reaction center.

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

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

  2. This excitation then migrates along the chain-biomolecules, 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:

  1. Inter BChl molecules (system-system) coupling (which is responsible for EET).

  2. The system-bath (BChl molecule-protein) coupling (which is responsible for decoherence).

Two important timescales:

  1. Excitation transfer time scale (usual in photosynthesis pigments).

  2. Decoherence time scale .

If the system-bath 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 system-bath 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).

Figure 2: Two extreme cases: (1) Upper: incoherent (classical diffusive motion), (2) Lower: ultra-quantum (delocalized excitation). The theoretical investigation of the intermediate case is a big challenge.

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 system-bath coupling case, the standard approach was the 2nd Born quantum master equation which is a perturbative quantum master equation (upto second order in system-bath interaction). Its Markovian and secular approximation is known as Redfield master equation(6). The 2nd Born quantum master equation can be obtained from the Nakajima-Zwanzig 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 non-equilibrium statistical mechanics. As the the light harvesting pigments fall in the intermediate regime (system-bath coupling is of the order of system-system 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 system-bath 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 degrees-of-freedom re-organize 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 Nakajima-Zwanzig projection operator technique by restricting the perturbation series upto second order in the system-bath interaction(6). In the following we will apply this master equation to a concrete model of a dimer (open two-state quantum system) which caricature the dynamics of decoherence in a typical photosynthetic system(9).

ii.1 Non-Markovian solution

Projection super-operators 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 Liouville-van Neumann equation (classical equivalent is the invariance of the ”extension” in phase space):


Here the interaction representation is used (see for details any standard refs(6)). is the system-bath interaction Hamiltonian (for the 2nd Born approximation . is the bath Hamiltonian. To construct the equation-of-motion for the relevant part of i.e., , one defines the super-operator:


called the projection super-operator, one also defines . By applying on the Liouville-van Neumann equation turn-by-turn, we get:


Solving the second equation formally for the irrelevant part , and inserting in the first, one obtains the required equation-of-motion 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 system-bath interaction(6)) and for a traditional dimer system(9):


master equation takes the form,


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 re-organization energy (the elastic energy related to the physical organization of the bath degrees-of-freedom). is the phonon Hamiltonian and is the system-bath 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 re-organization 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 an-harmonic terms are not important)) are:


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


Assuming the Drude-Lorentz model(9) for the bath spectral density, and assuming the high temperature approximation (), as appropriate for the FMO problem, we obtain


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 integro-differential 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:


with time evolution given as,


with and


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 integro-differential 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 integro-differential equations to a bigger set of ordinary differential equations:


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 integro-differential equations):


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 integro-differential equations(10). In the straightforward numerical method (traditional method) the integro-differential equation is first written as integral equation with double integration as converted to . The double integration is then done self-consistently with numerical integration(10).

Speed: On the Lenovo ThinkCentre-i7, 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).

Figure 3: Comparison of numerical results for a traditional method (blue-solid curve), to that introduced here (red-dotted curve). Parameters used are: . Abscissa is in fs (fs femto seconds).

ii.4 Markovian limit

As mentioned before in an important case of fast bath relaxation (when bath degrees-of-freedom re-organize 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 non-Markovian 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


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.

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:


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
Table 1: Regimes of validity of the Markov approximation
Figure 4: Sample verification of the analytic inequalities (Table I). Time evolution of population on site 1: blue (solid) curve is the non-Markovian solution and red (dotted) curve is the 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.

Figure 5: Time evolution of population on site 1 [] and the coherences [] [blue (solid) curve is the non-Markovian solution and red (dotted) curve is the Markovian approximation], for various values of (in cm). Other parameter are: = 100 cm, = 50 cm,   s.
Figure 6: Time evolution of population on site 1: blue (solid) curve is the non-Markovian and red (dotted) curve is the Markovian) for various values of the reorganization energy and inter-site coupling .
Figure 7: Time evolution of population on site 1 [blue (solid) curve is the non-Markovian and red (dotted) curve is the Markovian] for various values of the reorganization energy and phonon relaxation time .

First consider the non-Markovian 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 re-organization 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 re-organization 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 non-Markovian 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 photo-excitation, bath (nuclear co-ordinates) re-organize very fast representing an ”apt” bath (right hand picture) while in the non-Markovian 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 non-Markovian regime.

Figure 8: Approximate Physical Picture

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 photo-excitation. We have considered this problem in the second reference of the list(8) which considers photo-excitation of a dimer oscillator system with an ultra-short 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 damping-out 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 system-bath coupling () is of the same order of magnitude of the system-system 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 non-Markovian 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 non-perturbative 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 system-bath coupling Hamiltonian. But this re-normalize the system Hamiltonian. Then one re-partition 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 two-level electronic systems (Fig. 9). These two-level systems are electronically coupled with each other with coupling .

Figure 9: Two interacting molecules

Due to the electronic coupling between the molecules the exciton will transfer back-and-forth between the molecules. This will happen forever if the molecules are completely isolated—pure oscillatory quantum motion. Now consider that our two-molecular system is open i.e., interacting with the external bath degrees-of-freedom—-the phonons. It is well known that the dynamics remain quantum at short time scales and becomes classical at longer time scales(16). This quantum-to-classical crossover happens at a critical time scale which is inversely proportional to system-bath coupling energy (). Large (system-bath coupling) means fast quantum-to-classical crossover and vice versa. We denote system-bath 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


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


Here, as mentioned before, is the strength of system-bath coupling (also known as dynamical disorder) measured in the units of frequency.

We start with Liouville-von-Neumann equation for the total density matrix,


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:


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


Here is the functional derivative. Using the properties of stochastic noise and with some simplification(see appendix A), we get


Finally, one has a set of coupled differential equations:


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 2-D photon echo spectra which occurs at (Fig. 10).

Figure 10: (Left) Schematic line broadening information in 2-D photon echo spectrum shown without cross peaks. The linewidth due to homogeneous and in-homogeneous broadening are in orthogonal directions as shown. (Right) The energy levels of the model system

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


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: