# Quantum master equations for non-linear optical response of molecular systems

## Abstract

Generalized master equations valid for the third order response of an optically driven multi-level electronic system are derived within Zwanzig projection formalism. Each of three time intervals of the response function is found to require specific master equation and projection operator. Exact cumulant response functions for the harmonic profiles of the potential energy surfaces leading to Gaussian spectral diffusion are reproduced. The proposed method accounts for the nonequilibrium state of bath at the border of intervals, and can be used to improve calculations of ultrafast non-linear spectra of energy transferring systems.

## I Introduction

Non-linear response theory forms the theoretical basis for description of many important physical phenomena and experimental techniques, in particular the multi-dimensional techniques of non-linear optical spectroscopy (1). In recent years, two-dimensional (2D) coherent spectroscopy (2), developed first in NMR (3), has been brought into the near infra-red and optical (4); (5) domains. The 2D Fourier transformed spectrum completely characterizes the third order non-linear response of a molecular ensemble in amplitude and phase (2) providing thus, the maximal information accessible in a three pulse optical experiment. Since its first experimental realization, it has yielded new insights into photo-induced dynamics of electronic excited states of small molecules (6), polymers (7), large photosynthetic aggregates (8) and even solid state systems (9). For instance, long living electronic coherence has been detected in energy transfer dynamics of photosynthetic proteins (10), a find that stimulated a renewed interest in the properties of energy transfer in biological systems and its arguably quantum nature.

Most 2D experiments are well described by the third order perturbation semiclassical light-matter interaction response function theory. Simulations must include dephasing of electronic coherence and electronic energy transfer between levels which are both induced by interaction of the electronic degrees of freedom (DOF) with environment (bath). Response functions of model few level systems with pure dephasing (no energy transfer between the levels) can be calculated using second cumulant in Magnus expansion (1), and expressed in terms of the energy gap correlation function (EGCF), . For Gaussian stochastic bath, this result is exact, and the EGCF of the electronic transitions thus determines fully both the linear and the non-linear (transient) absorption spectra. For energy transferring systems, such as Frenkel excitons in photosynthetic complexes (11), the exact response functions cannot be constructed by second cumulant. Photosynthetic complexes are relatively large, and the proper methods to simulate finite timescale stochastic fluctuations at finite temperatures (12) carry a substantial numerical cost. Stochastic simulations can thus be readily implemented only for small systems (13); (14). An explicite inclusion of some important environmental modes into the Hamiltonian also increases the size of the system beyond feasibility. Practical calculations thus require some type of reduced dynamics where only electronic DOF are treated explicitly. A convenient form of the reduced dynamics, master equations, is realized by deriving kinetic equations for expectation values of relevant quantities averaged over the DOF of the bath (15); (16). If these quantities (such as transition dipole moment) do not depend on the bath DOF, all relevant information is carried by the reduced density matrix (RDM).

It was demonstrated that RDM master equations derived by projection operator technique reproduce the linear response exactly (17). For higher order response, however, the same approach neglects bath correlation between different time intervals of the molecular photo-induced evolution separated by interaction with laser pulses (see e.g. Section 3 in Ref. (18)). In other words, standard RDM master equations do not properly account for the non-equilibrium state of the phonon bath present after the second and the third laser pulse. This neglect can lead to a complete loss of the experimentally observed dynamics in simulated 2D spectra, such as in the case of the vibrational line shape modulation of a single electronic transition (19). This modulation interplays with simultaneous effect of electronic coherence. Accounting reliably for bath correlations among intervals is therefore of utmost importance for interpretation of the role of quantum coherence in natural light-harvesting. The difference between the 2D lineshape calculated with exact formula, and by RDM, i. e. neglecting the correlations, is demonstrated in Fig. 1. We employed the definitions 2D spectra from Ref. (19) and the response theory from Ref. (1).

In this Letter, we show that the exact third order response function
of a multilevel molecular system with pure dephasing can be calculated
by certain generalization of RDM master equations. We introduce special *parametric
projection operators* that enable us to derive distinct equation of
motion for each time interval of the response function. These parametric
projectors offer a systematic approach to improvement of current master equation method by
incorporating bath correlation effects into the calculation of non-linear response in energy transfer systems.

## Ii Model Hamiltonian

We model a molecular aggregate (e.g. photosynthetic antenna) as an assembly from electronic two-level molecules coupled by resonance interaction. Third order optical experiments on such systems address only the collective electronic ground state and the states with one or two excited molecules of the aggregate. Thus the -molecule aggregate has relevant excited states. The free evolution of the system between the times of interaction with external light pulses is described by Frenkel exciton Hamiltonian

(1) |

Here, represents kinetic energy of the bath DOF, and are the electronic energies of the ground- and the excited states, respectively, is the energy gap operator, and and are potential energy surfaces (PES) in the electronic ground- and excited states, respectively. Harmonic PES represents Gaussian fluctuations of excitons. The energy gap operator is defined so that it vanishes in the ground state equilibrium . Here, the ground state bath Hamiltonian reads . The operators describe a macroscopic number of DOF so that one cannot easily treat them explicitely in numerical calculations. Resonance coupling

(2) |

between the excited states causes electronic energy transfer between molecules. Diagonalization of average Hamiltonian is required to identify the positions of peaks in absorption spectrum, which correspond to transition frequencies between electronic eigenstates with eigenenergies , and the ground state. In the very same spirit, 2D spectra are conveniently interpreted in terms of energy transfer between these energetic eigenstates. When resonance coupling between molecules in the aggregate vanishes the eigenstates coincide with excited states , and so do the eigenenergies and energies .

## Iii Linear response

For harmonic PES and the calculation of response functions of all orders can be carried out exactly by second order cumulant. The linear response function

(3) |

is generated by time evolution of coherence elements of the RDM , which can be differentiated into the convolutionless master equation (17); (20)

(4) |

In Eq. (3), are the transition dipole moments between eigenstate and ground state , the line broadening function

(5) |

is related to the EGCF,

(6) |

in the ground state equilibrium, and . Further on in this work, we assume that the electronic energy gap operators are nonzero.

## Iv Master equations for non-linear response

In non-linear experiment, molecular system is probed by three short laser pulses separated by time intervals , , and its response is monitored after time . Non-linear response function is conveniently dissected into Liouville space pathways (1), each representing one combination of indices of RDM during the periods , , and . For our discussion we select a pathway usually denoted as (see Fig. 2 for the corresponding double sided Feynman diagram), i.e. the pathway which is responsible for photon echo signal, and which contains excited state evolution during the second interval (waiting time) (1). The response function involving excited states and reads

(7) |

Here, are abbreviated matrix elements of the evolution superoperator which is the Green’s function solution of the Liouville-von Neumann equation with Hamiltonian, Eq. (1). Between each two consecutive Green’s functions in the trace of Eq. (7), the action of the external light causes transition between electronic levels , and . For a harmonic PES the second order cumulant is exact, and it yields (21)

(8) | |||||

By differentiating Eq. (8) with respect to we find the response functions to be solutions of evolution (master) equations

(9) |

with the relaxation tensor

(10) |

and initial condition . Master equation can be similarly found in the second interval by differentiation of

(11) |

with relaxation tensor

(12) |

and initial condition . In the first interval we find an analogical master equation for which is (up to a complex conjugation) identical to Eq. (4). We have thus demonstrated that evolution equations can be constructed with the structure of the convolutionless master equation, and that they yield exact third order response function. Their forms, however, have now to be specific to each interval and to each Liouville space pathway.

## V Constructing master equations by parametric projectors

In Ref. (17) it was pointed out that Eq. (4)
agrees with *convolutionless master equation* (CLME) of Hashitsume
et al. (22), when it is evaluated to second order in .
Similar microscopic derivation of Eqs. (9) to
(12) by using projection techniques would open a way
to constructing improved master equations for the electronic energy transfer
in the cases where no exact cumulant solution is known. The Zwanzig
projection technique requires to define projection superoperator by
its action on a arbitrary operator . In general form ,
where the normalized bath matrix ()
can be chosen to achieve specific goals. Master equation
for the linear response, Eq. (4), is homogenous and it corresponds
to the choice , when
(), and thus the inhomogeneous part of the
CLME vanishes (17) for uncorrelated initial state .
In the second and third interval the choice
will not yield Eqs. (9) and (11), since the inhomogeneous
part of the CLME does not vanish. In fact, the inhomogeneous terms
cannot be eliminated in all three intervals simultaneously by any
single projector. This is a microscopic counterpart of the fact that
non-linear optical response function cannot be exactly evaluated by
simulating the reduced dynamics of the system by a single master equation.
Indeed, the homogeneous part of the CLME corresponds to only one of
the four (nonzero) contributions of the third order response function

(13) |

The master equations (4) to (12)
can be derived when we allow different projection operator to be defined
for each time interval of the response function. Let the projection
operators , and in the respective intervals be chosen so that they
eliminate the inhomogeneous terms in the CLME master equation, i.e.
so that ; , etc (see Fig. 2).
The corresponding projectors read as follows. The
*first coherence projection superoperator*

(14) |

is the widely used Argyres-Kelly projector (23).
The second interval can be treated with the *population projection
superoperator*

(15) |

which yields a single master equation for all RDM elements in the
second interval. Projector, Eq. (15), to be applied to
a single RDM element as in Eq. (11), reads .
In the third interval of the response a set of *second coherence
projection superoperators* reflecting different bras ()
during the second interval have to be defined as

(16) |

The factors and are present to ensure that the superoperators have the projection property and . From these conditions it follows e.g. . The factor evaluates in the second cumulant to .

An alternative to this formalism would be to account directly for the terms in Eq. (13), and to use the time convolutionless operator identity (22) for the derivation of the non-linear response as in Ref. (26). The two approaches lead to the same exact results for a system of non-interacting molecules. For a coupled system, the quality of the equations of motion derived by the two approaches has to be compared in a separate study. However, the central idea of the projector approach, i.e. its choice in a form which satisfied , is general, and it does not limit the present method to the three projectors suggested above. Depending on physical situation, we expect that specialized projectors can be introduced for interacting systems, possibly based on the comparison with the results of Ref. (26).

## Vi Recovering the cumulant expressions

Let us now demonstrate how these projection operators can be used for calculation of the non-linear response function, Eq. (7), in terms of the reduced dynamics of the system. We can write

(17) |

where the three reduced evolution superoperators denoted by are solutions of three different master equations following from the projectors, Eqs. (14) to (16). To simplify the matters, we will concentrate on the response function where , i.e. on . Since we deal with a multi-level electronic system with adiabatically separated states (), we can assume that and thus

(18) |

Note that all -dependence in Eq. (18) is due to the second coherence parametric projector.

Now we construct the master equations for the first and
second coherence intervals (i.e. in intervals characterized by time and , respectively), and demonstrate that they lead to correct
response function.
First, we split the total Hamiltonian into a sum of the pure bath
(denoted by *B*), pure electronic (*S*) and system–bath
interaction (*S-B*) parts. Bath Hamiltonian has already
been defined above, the electronic Hamiltonian reads ,
and the system–bath coupling Hamiltonian is
We define the Liouville superoperator using the
corresponding Hamilton operator as
and switch to the interaction picture induced by .
In the second order in ,
and assuming that for , and we obtain the following quantum master equation

(19) |

The time-dependence of the Liouvillian and the upper index of the density operator denote interaction picture with respect to the Hamiltonian

In the first coherence interval we use the projection superoperator . The construction of the operator ensures that , and so the first term in Eq. (19) contributes with zero. The second term is only non-zero when no occurs on the right hand side of the density operator, and thus only the term contributes. The evaluation leads to Eq. (4) and consequently yields A straightforward application of the second coherence projection superoperator leads to

(20) |

where

(21) |

and

(22) |

The coefficients and , Eqs. (21) and (22), can be recast into the second cumulant by a technique described in Ref. (24). A tedious but straightforward calculation gives

(23) |

and

(24) |

Combining the results of the second order cumulant expansion of Eqs. (21) and (22) and switching back from the interaction picture we arrive at the master equation Eq. (10). The solution of the master equation can be easily obtained in form of the time evolution superoperator . Finally, inserting and into Eq. (18) we obtain Eq. (8) which is the desired result. Calculation can be repeated for any Liouville space pathway reproducing the results of direct cumulant expansion.

Note that standard choice of ground state equilibrium for projector would yield response function

(25) | |||||

Comparison between 2D lineshapes simulated using Eq. (8) and Eq. (25) is made in Figure 1.

The projector technique developed above can be extended to finite resonance coupling by straightforward abandoning of the limit. We will show elsewhere (25) that provided the optical coherences in the first interval can be assumed independent of each other (secular approximation), the projectors, Eqs. (14) and (15) can be used to derive corresponding master equations. The complete response function involves a single master equation for the first coherence interval and two master equations for the waiting interval. In the second coherence interval of the response function there are master equations in each Liouville pathways, where is the number of single exciton states. When energy transfer between system’s eigenstates is induced by non-zero couplings, the application of the parametric projectors represents an approximation whose quality depends on the previous energy transfer history of the system. Nevertheless, unlike the usual master equations based on standard projectors, the parametric projectors lead to correct limit in . The prescription for construction of the parametric projection operators is general, and can be used to construct a proper projector in the cases where would fail.

## Vii Conclusions

In this Letter, we have outlined a general method for calculation of non-linear response functions by quantum master equations. We have tailored parametric projection superoperators to derive master equations for each time interval of the response function, and we demonstrated the validity of the method by reproducing a well-know result exact for a multi-level system. The results presented here can be extended to a master equation method for calculation of a general multi-point correlation functions. This opens a new research avenue with consequences in the theory of open systems, non-linear spectroscopy and potentially other branches of physics.

###### Acknowledgements.

This work was supported by the Czech Science Foundation (GACR) grant 205/10/0898 and by the Ministry of Education, Youth, and Sports of the Czech Republic through grant ME899 and the research plan MSM0021620835.### References

- S. Mukamel, Principles of Nonlinear Optical Spectroscopy (Oxford University Press, Oxford, 1995).
- D. M. Jonas, Annu. Rev. Phys. Chem. 54, 425 (2003).
- R. R. Ernst, G. Bodenhausen, and A. Wokaun, Principles of Nuclear Magnetic Resonance in One and Two Dimensions (Clarendon Press, Oxford, 2004).
- M. Cho, Two-Dimensional Optical Spectroscopy (CRC Press, Boca Ranton, 2009) .
- P. Hamm and M. Zanni, Concepts and Methods of 2D Infrared Spectroscopy, (Cambridge University Press, Cambridge, 2011).
- F. Milota, J. Sperling, A. Nemeth, T. Mančal, and H. F. Kauffmann, Acc. Chem. Res. 42, 1364 (2009).
- E. Collini and G. D. Scholes, Science 323, 369 (2009).
- N. S. Ginsberg, Y.-C. Cheng, and G. R. Fleming, Acc. Chem. Res. 42, 1352 (2009).
- K. W. Stone, K. Gundogdu, D. B. Turner, X. Li, S. T. Cundiff, and K. A. Nelson, Science 324, 1169 (2009).
- G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Mančal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446, 782 (2007).
- H. van Amerongen, L. Valkunas, and R. van Grondelle, Photosynthetic Excitons (World Scientific, Singapore, 2000).
- Y. Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006).
- A. Ishizaki and G. R. Fleming, Proc. Nat. Acad. Sci. USA 106, 17255 (2009).
- F. Šanda and S. Mukamel, J. Phys. Chem. B 112, 14212 (2008).
- F. Rossi and T. Kuhn, Rev. Mod. Phys. 74, 895 (2002).
- V. M. Axt and T. Kuhn, Rep. Prog. Phys. 67, 433 (2004).
- R. Doll, D. Zueco, M. Wubs, S. Kohler, and P. Hanggi, Chem. Phys. 347, 243 (2008).
- A. Ishizaki and Y. Tanimura, Chem. Phys. 347, 185 (2008).
- A. Nemeth, F. Milota, T. Mančal, V. Lukeš, H. F. Kauffmann, and J. Sperling, Chem. Phys. Lett. 459, 94 (2008).
- V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems (Wiley-VCH, Berlin, 2001).
- D. Abramavičius, B. Palmieri, D. V. Voronine, F. Šanda, and S. Mukamel, Chemical Reviews 109, 2350 (2009).
- N. Hashitsume, F. Shibata, and M. Shingu, J. Stat. Phys. 17, 155 (1977).
- P. N. Argyres and P. L. Kelley, Phys. Rev. A 134, 98 (1964).
- W. M. Zhang, T. Meier, V. Chernyak, and S. Mukamel, J. Chem. Phys. 108, 7763 (1998).
- J. Olšina and T. Mančal, in preparation (2011).
- M. Richter, and A. Knorr, Ann. Phys. 325, 711 (2010).