Parametric projection operator technique for second order non-linear response

Chemical Physics 404 (2012) 103–115

Contents lists available at SciVerse ScienceDirect

Chemical Physics journal homepage: www.elsevier.com/locate/chemphys

Parametric projection operator technique for second order non-linear response Jan Olšina, Tomáš Mancˇal ⇑ Faculty of Mathematics and Physics, Charles University in Prague, Ke Karlovu 5, CZ-121 16 Prague 2, Czech Republic

a r t i c l e

i n f o

Article history: Available online 20 March 2012 Keywords: Non-linear optical response Open quantum systems Projection operator techniques

a b s t r a c t We demonstrate the application of the recently introduced parametric projector operator technique to a calculation of the second order non-linear optical response of a multilevel molecular system. We derive a parametric quantum master equation (QME) for the time evolution of the excited state of an excitonic system after excitation by the first two pulses in the usual spectroscopic four-wave-mixing scheme. This master equation differs from the usual QME by a correction term which depends on the delay s between the pulses. In the presence of environmental degrees of freedom with finite bath correlation time and in the presence of intramolecular vibrations we find distinct dynamics of both the excited state populations and the electronic coherence for different delays s. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Recent years have seen a rapid experimental development in multidimensional coherent spectroscopy. The two-dimensional (2D) Fourier transformed spectrum completely characterizes the third order non-linear response of a molecular ensemble in amplitude and phase [1,2] providing thus, the maximal information accessible in a three pulse experiment. In the three pulse degenerate four wave mixing (FWM) experiment [3], two independent adjustable parameters – the delays between the three pulses – allow experimentalists to address the time and delay dependent third order non-linear response of the molecular system Sð3Þ ðt; T; sÞ via complete characterization of the third order non-linear signal Eð3Þ s ðt; T; sÞ (see Fig. 1). This information is conveniently represented by a 2D plot of the Fourier transformed signal Eð3Þ s ðxt ; T; xs Þ which correlates the absorption frequency xs with the frequency xt of stimulated emission, ground state bleaching and excited state absorption contributions, separated from the absorption event by an adjustable waiting time T. Femtosecond photo-induced evolution of molecular systems is reflected in the amplitude, position and shape of 2D spectral features as functions of the waiting time revealing thus important information about the molecular structure, intramolecular dynamics, as well as the interaction of molecular electronic transitions with their environment [4–6]. By analyzing 2D spectra of photosynthetic Fenna–Mathews–Olson (FMO) protein, signatures of electronic quantum coherence in this system have been theoretically predicted [5] and experimentally verified [7]. The same measurements led to a surprising discovery of the long life time of these electronic coherence

⇑ Corresponding author. E-mail address: [email protected] (T. Mancˇal). 0301-0104/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.chemphys.2012.02.022

and to an identification of a host of other effects which contribute to the wave-like nature of the energy transfer in these systems. For the description of the most 2D experiments, the third order semi-classical light-matter interaction response function theory is well established. Response functions of model few-level systems with pure dephasing (i.e. with no energy transfer between the levels) can even be expressed analytically in terms of the so-called energy gap correlation function (EGCF), using the second order cumulant in Magnus expansion [3]. Some examples of small chromophores in solution fall in this category [8,9] when investigated on time scales shorter then radiative life time. For a Gaussian bath this analytical theory is exact, and thus knowing the EGCF of the electronic transitions enables us to determine linear (absorption) as well as non-linear spectra. However, the construction of exact response functions for realistic energy transferring systems, such as Frenkel excitons in photosynthetic aggregates [10], is no longer possible. Photosynthetic complexes are relatively large, and the proper methods to simulate finite timescale stochastic fluctuations at finite temperatures [11– 13] carry a substantial numerical cost. Practical calculations thus require some type of reduced dynamics where only electronic degrees of freedom (DOF) are treated explicitly. These approaches usually rely on a host of approximations that seem to work well for most spectroscopic techniques [14–17] and in systems where the details of system-bath coupling are apparently less important [18]. With increasing details of the excited state dynamics revealed by the 2D spectroscopy [19–21] and with increasing size of the studied systems, it becomes more and more important to keep the numerical cost of simulations low, while simultaneously account for experimentally observed quantum effects. One of the possible research directions is to extend on existing reduced density matrix (RDM) theories [22,23] or relax certain approximations [24].

104

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

Fig. 1. A scheme of a typical FWM experiment. Three laser pulses with delays s and T are incident on a sample. A non-linear signal generated into a selected direction, k1 þ k2 þ k3 here, is detected.

description of photosynthetic systems and their spectroscopy [10,4]. The second order response, which we have in mind here, corresponds to the interaction of the excitonic system with first two pulses in a standard FWM spectroscopy (see Fig. 1), and it should correspondingly describe the time evolution of the system in the excited or ground state. The part of the response associated with evolution in the excited state also corresponds to the material quantities that govern excitation of molecular systems by arbitrary quantum light [28]. For a second order response we need to consider only single exciton band, i.e. the collective excited states of the aggregate with only a single aggregate member excited. The Frenkel Hamiltonian including environmental contributions has the general form H ¼ HS þ HB þ HSB , where the purely electronic (system) Hamiltonian HS reads

HS  eg jgihgj þ

P n

It was demonstrated that RDM master equations derived by projection operator technique reproduce exactly the linear response [25]. In the case of higher order response functions, however, the same approach necessarily neglects bath correlations between the different periods of photo-induced system evolution (i.e. between so-called coherence time s and the waiting time T). The failure to account for this correlation leads sometimes to a complete loss of the experimentally observed dynamics in simulated 2D spectra, such as in the case of the vibrational modulation of electronic 2D spectra [8]. Because vibrational modulation leads to effects similar to those attributed to electronic coherence, developing methods that can account for its effect reliably in complex systems is of utmost importance. One possible approach to the problem is to derive equations of motion for the response as a whole, and to take the previous time evolution of the system into account explicitly [26]. In this paper, we take a different route in which the correlation effects are treated by a specific choice of the projection operator. Once a projector is specified, it can be used for any method of treating system bath interaction. We apply the previously suggested parametric projection operator technique [27] to a calculation of the second order response of a quantum system. The second order response operators can be used to determine the state of a molecular system subject to excitation by a weak light with arbitrary properties, and thus its importance goes beyond the semi-classical system light interaction theory [28]. Alternatively, calculation of the second order response can be viewed as a first step towards more involved calculation of the third order response functions that is required for the modeling of the third order non-linear spectra. The paper is organized as follows: in the next section we specify the model Hamiltonian and the description of the system’s interaction with the light and its environment (the thermodynamic bath). In Section 3 we discuss the third order response functions of a multilevel systems, and we point out general limitations of their evaluation by the reduced density matrix master equations. In Section 4 we write out the density operator describing the excited state of a molecular system in terms of the second order response operator, and we discuss the advantages of usage of the parametric projection operator over the standard projectors. The details of the parametric projection operator and the corresponding master equation are introduced in Section 5. Numerical results for the excited state dynamics are then discussed in Section 6. Some definitions and calculation details can be found in the appendices.

2. Model system and block formalism Our primary interest lies in photosynthetic aggregates of chlorophylls. Frenkel exciton model is well established for the



en þ hUen  Ug i jnihnj þ

P jmn jmihnj;

ð1Þ

mn

the bath Hamiltonian HB can be identified with

HB  T þ Ug ;

ð2Þ

and the system-bath interaction HSB reads as

HSB 

P

DUn jnihnj:

ð3Þ

n

Here, jgi is the electronic ground state of the aggregate of molecules with energy eg ; jmi are electronic states with the mth molecule of the aggregate of N molecules excited with energy en ; jnm represents resonance coupling between excited states jmi and jni; U represents potential energy surfaces (PES) of the bath DOF in corresponding electronic states, and T is the kinetic energy of the bath DOF. We also defined so-called energy gap operator

DUn  Uen  Ug  hUen  Ug i;

ð4Þ

where hi represents the equilibrium quantum mechanical average over the bath DOF. Symbol  denotes an arbitrary operator. The interaction of the system with light will be described by the usual semi-classical light-matter interaction term in dipole approximation

HSE ðtÞ ¼ l  EðtÞ:

ð5Þ

Here, l is the transition dipole moment operator of the aggregate system and EðtÞ ¼ eEðtÞ is the electric field vector of the light with polarization vector e. The states jmi with excitations localized on individual member molecules of the aggregate allow us to define the properties of the aggregate based on crystal structure and quantum chemical inputs. It is, however, often more convenient to work in a basis of the ~ of the electronic Hamiltonian HS . For later conveeigenstates jmi nience we will define projection operators K m ¼ jmihmj projecting ~ mj ~ projecting on on a state in a state in the local basis and K m~ ¼ jmih a state in the eigenstate (excitonic) basis. For the discussion of the non-linear response functions, it is advantageous to use a Liouville space superoperator formalism. We define the dipole moment superoperator as

V ¼

1 1 ½l;  ¼ ½l  e;  ; h  h

ð6Þ

and the Liouvillian as L ¼ 1h ½H;  . We denote UðtÞ the evolution operator and we use the evolution superoperator UðtÞ where convenient. In the second order system-bath coupling theories, the energy gap operators, Eq. (4), are present in the equations via the EGCFs

C mn ðtÞ  hU y ðtÞDUm UðtÞDUn i:

ð7Þ

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

Double integration of C mn ðtÞ over time yields so-called lineshape functions

g mn ðtÞ 

1 h

Z

2

t

ds

0

Z s

ds0 C mn ðs0 Þ;

ð8Þ

0

which determine absorption and emission line shapes. Later in the text, it will be convenient to define the energy gap operators in the excitonic basis

DUa~b~ ¼

P

~ DUk ha~jkihkjbi;

ð9Þ

k

the excitonic EGCFs

C a~b~ ðtÞ  hU y ðtÞDUa~a~ UðtÞDUb~b~ i ¼

P

~ 2 C ðtÞ ~jkij2 jhbjlij jha kl

ð10Þ

kl

and the excitonic lineshape functions

g a~b~ ðtÞ 

Z

1 h

2

t

ds

Z s 0

0

ds0 C a~b~ ðs0 Þ:

ð11Þ

We are interested in the influence of the bath DOF on the dynamics of the electronic DOF in the excited state. The bath can consists of several types of collective modes which can be described by various types of bath correlation function. Experimental situation is typically well described by one or several overdamped Brownian oscillator modes standing for an macroscopic number of harmonic DOF, and several underdamped oscillator modes standing for some important vibrational coordinates, e.g. normal modes of the chromophore molecules [8]. Both of these limits can be conveniently described by the general Brownian oscillator correlation function CðtÞ (see [3]) which is a function of temperature T, frequency XBath of the oscillator, damping coefficient c and the reorganization energy k. In the case of a strongly overdamped mode i.e. when c  2Xbath the formula simplifies significantly. At high temperatures we have

CðtÞ ¼ hkK cotðKhb=2Þ expðKtÞ  ihkK expðKtÞ;

ð12Þ

with K ¼ 1=c . In the opposite case of a non-damped oscillator c ! 0 the correlation function reads

CðtÞ ¼ kXBath h½cothðbhXBath =2Þ cosðXBath tÞ  i sinðXBath tÞ:

ð13Þ

The Hamiltonian, Eq. (1), has a block form. In case we would consider excited states with up to N excitations per aggregate a more general Hamiltonian would have to be written including N þ 1 blocks with one block containing just the ground state jgi, and the N blocks containing the excited states with one, two, etc, up to N excitations. The block structure is enabled by the fact that the Hamiltonian operators, Eqs. (1)–(3), do not contain any terms that enable de-excitation of excited states. This is a well justified assumption in studies of ultra-fast dynamics of the chlorophyll based photosynthetic aggregates [29]. The transitions between different blocks of the Hamiltonian are only enabled by the interaction with the light, Eq. (5). The second order response, which is of interest in this paper, requires on the ground- and the singleexciton blocks. We will therefore use upper indices g (ground state) and e (single exciton block) to denote different blocks of operators and superoperators whenever we are interested in operations on and evolution of their individual blocks. Thus, e.g. the time evolution of the system in excited state single exciton manifold is described by the RDM block qðeeÞ ðtÞ ¼ U ðeeeeÞ ðtÞqðeeÞ ð0Þ ¼ U ðeeÞ ðtÞqðeeÞ ð0ÞU ðeeÞy ðtÞ, and the action of the dipole operator (superoperator) promotes the ground state block into a coherence block as qðegÞ ðtÞ ¼ lðegÞ qðggÞ ðtÞ ¼ V ðegggÞ qðggÞ ðtÞ. When working in a particular basis of states (e.g. fjmig ¼ fj1i; j2i; . . .g) we can write out the sums over the states explicitly, e.g. as ðeeÞ qmn ðtÞ ¼

P kl

ðeeeeÞ

ðeeÞ

U mnkl ðtÞqkl ð0Þ ¼

P ðeeÞ ðeeÞ ðeeÞy U mk ðtÞqkl ð0ÞU ln ðtÞ:

105

non-linear experiments, the block structure has to be considered up to the two-exciton block. 3. Multi-point correlation functions in non-linear spectroscopy The problem of system-field interaction described by the Hamiltonian, Eq. (5), can be solved either by exact solution of equations of motion or perturbatively with respect to the electric field. In the context of non-linear spectroscopy, the perturbative approach is usually preferred since the order of the perturbation can be very well controlled as an experimental parameter. In the perturbative solution, the external field effect is taken into account by a response function which has the form of a multi-point correlation functions (MPCF). The field then does not enter the equations of motion. The MPCF of an nth order describes the response of the system to n delta pulses of unit intensity. The nth order response is then calculated as an n-fold convolution of the response function with the external field [3]. For a two band system (ground state and single excitons), the third order non-linear spectroscopy is completely described by four response functions [3], listed in Appendix A. As an example, we will consider the response usually denoted as R2 , Eq. (A2). We can notice that the block (upper) indices of the evolution operators in Eq. (A2) (read from left to right) follow the double-sided Feynman diagram in Fig. 2A (see [3] for details). During the so-called population interval T of the response, the system evolves in the excited state band (ee). If no resonance coupling is present (jmn ¼ 0 in Eq. (1)) each response function can be split into a sum of independent components. For R2 this means

R2 ðt; T; sÞ ¼

P R2;kl ðt; T; sÞ;

ð15Þ

kl

where the components

R2;kl ðt; T; sÞ ¼ jlkg j2 jllg j2 n o ðegegÞ ðeeeeÞ ðgegeÞ ðggÞ  trB U kgkg ðtÞU klkl ðTÞU glgl ðsÞW eq

ð16Þ

can be evaluated analytically in terms of the line shape functions, Eq. (8). Here, W ðggÞ eq ¼ weq jgihgj is the total equilibrium density operator (the electronic energy gap is assumed to be much higher than kB T) and weq represents the equilibrium density operator of the bath alone. For jmn > 0 no analytical result is available, and one needs to resort to some master equation simulations of the evolution superoperators U. However, the master equations can only deliver certain reduced (i.e. averaged over the state of the bath) version of this superoperator. The usual way of deriving master equations for the RDM is to apply a projection operator P which reduces the full density matrix WðtÞ of the system to its selected part [30], in our case, to the electronic DOF

ð14Þ

kl

The rules are rooted in the simple well-known fact that matrices can be multiplied by blocks. When discussing the most common

Fig. 2. Double-sided Feynman diagrams of the third and the second order response. Part A: The Feynman diagram of the Liouville pathway R2 . Part B: In full lines, the Feynman diagrams of the response operators RII (left) and RI (right). The dashed part completes the diagram into a corresponding third order response.

106

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

qðtÞ ¼ trB fWðtÞg:

ð17Þ

The most popular prescription, the so-called Argyres–Kelly (AK) projector [31], reads as

P 0 WðtÞ ¼ qðtÞweq :

ð18Þ

It is easy to see that inserting the identity 1 ¼ P þ Q, where Q ¼ 1  P into Eq. (17) leads to

qðtÞ ¼ trB fUðtÞPWð0Þ þ UðtÞQWð0Þg:

ð19Þ

Provided that QWð0Þ ¼ 0 (a condition satisfied in the first interval of the response function, Eq. (A2)) one can define the reduced evoe ðtÞ ¼ trB fUðtÞweq g such that qðtÞ ¼ U e ðtÞqð0Þ. lution superoperator U e ðtÞ can be calculated from a master The evolution superoperator U equation that can be derived by projection operator formalism [30]. Let us apply the same method to higher order response functions. Let use assume we have derived a master equation by applying the corresponding projector operator P and calculated e ðtÞ. We can then elements of the reduced evolution superoperator U assemble an approximate response function

e 2;kl ðt; T; sÞ ¼ U e ðegegÞ ðtÞ U e ðeeeeÞ ðTÞ U e ðgegeÞ ðsÞqðggÞ ð0Þ: R kgkg klkl glgl

ð20Þ

Here, we set all transition dipole moment elements to one for the e 2;kl can also be written as sake of brevity. It can be shown that R

n o e 2;kl ðt; T; sÞ ¼ trB U ðegegÞ ðtÞPU ðeeeeÞ ðTÞPU ðgegeÞ ðsÞqðggÞ ð0Þ : R kgkg klkl glgl

4. Second order response to light The second order optical response can yield an optical signal on the sum or difference frequencies of the exciting field. Here, we consider only the degenerate experiment when the excitation field have the same frequency, and we are interested in particular in the zero frequency part of the response which does not generate an optical signal. Unlike in spectroscopy where we study a quantum mechanical expectation value (polarization or field), here we would like to study the state of the system achieved by the excitation. Correspondingly, the second order response will not be expressed in terms of response function, but rather in terms of some response operators. By the ‘‘system’’ we mean electronic DOF, and thus the response operators correspond to the reduced density matrix of the system in the same way the response function corresponds to an expectation value of e.g. polarization. 4.1. Excited state of a molecular system

ð21Þ

However, the exact expression for R2;kl reads as n o ðegegÞ ðeeeeÞ ðgegeÞ R2;kl ðt; T; sÞ ¼ trB U kgkg ðtÞðP þ QÞU klkl ðTÞðP þ QÞU glgl ðsÞqðggÞ ð0Þ :

ð22Þ Eqs. (21) and (22) differ by Q-containing terms that cannot be in general eliminated, and consequently one cannot expect master equations based on a single projector operator P of any type to reproduce the third order response functions. This applies also to the validity of non-perturbative schemes of calculations of non-linear response, such as those derived in Refs. [32,33]. Note that we do e ðtÞ in not claim that master equations cannot lead to exact form of U a single interval. They sometimes do, even if they are based on finite order perturbation in system-bath coupling [25]. This exactness is enabled by the fact that the solution of the finite order equation contains all orders of perturbation

Rt @ qðtÞ yields dsKðsÞ ¼ KðtÞqðtÞ ! qðtÞ ¼ e 0 qð0Þ: @t

order convolutionless QME (CL-QME) as in Ref. [26]. While both approaches should lead to similar results, the projection operator technique is not limited to any specific way of expanding the equations of motion in terms of system-bath coupling, and we believe it is therefore somewhat more flexible.

ð23Þ

What we address here is their inability to compensate the Q-containing terms between intervals of the response function, i.e. to take into account the effect of the bath history on the change of electronic state after the delta pulse excitation. A general solution of this problem was proposed in Ref. [27]. It was argued that one cannot write down a single exact master equation for all three intervals of the response function. Rather, one has to write a different master equation for each interval. This is formally possible by introducing three different projectors P 0 (i.e. standard AK projector) for the first, P s for the second and PsþT for the third interval of the response [27]. The projectors P s and P sþT are constructed so as to cancel the Q-containing term exactly for jmn ¼ 0. In this limit, all response functions, Eqs. (A1)–(A4), can be exactly reproduced by the corresponding master equations. In the following sections, we will treat the case of jmn – 0, for the case of the second order response, i.e. we will derive equations of motion for the reduced density matrix using the projector P s . The application of this approach to the full response function will be treated elsewhere. An alternative to the projector approach is to attempt to derive specific equations of motion for each response function as a whole in a specific perturbation scheme, e.g. second

Assuming that we excited the system by some external field EðtÞ, the second order reduced density operator of the system reads as

q^ ð2Þ ðtÞ ¼ 

Z

1

dt 1

0

Z

1

  dt2 trB Uðt2 ÞVUðt 1 ÞVW eq jgihgj

0

 Eðt  t 2 ÞEðt  t2  t1 Þ:

ð24Þ

^ ð2Þ then The excited state part of the second order density operator q yields

qð2Þ ðtÞ ¼

Z

1

Z

1

 dt2 RI ðt2 ; t 1 ÞA ðt  t 2 ÞAðt  t2  t1 Þeixt1  þ RII ðt2 ; t 1 ÞAðt  t2 ÞA ðt  t2  t1 Þeixt1 ; 0

dt1

0

ð25Þ

where

  RI ðt2 ; t 1 Þ ¼ Tr U ðeeeeÞ ðt 2 ÞV ðeeegÞ U ðegegÞ ðt1 ÞV ðegggÞ weq jgihgj ;

ð26Þ

and

  RII ðt 2 ; t 1 Þ ¼ TrB U ðeeeeÞ ðt 2 ÞV ðeegeÞ U ðgegeÞ ðt 1 ÞV ðgeggÞ weq jgihgj ;

ð27Þ

are the second order response operators. In Eq. (25) we introduced the electric field envelopes AðtÞ and the carrier frequency x so that EðtÞ ¼ AðtÞeixt þ c:c. The interaction of the molecular system with arbitrary light can be expressed using the operators RI and RII as discussed in Ref. [28]. For a general state of the light jwrad i, the terms A ðt  t 2 ÞAðt  t 2  t1 Þeixt1 have to be replaced by the light correlation function Iðt  t 2 ; t  t2  t 1 Þ ¼ hwrad j b E y ðt  t 2 Þ b Eðt  t2  t 1 Þjwrad i, where b EðtÞ is the operator of the electric field in Heisenberg representation. This allows us to calculate the state of a system excited by light with arbitrary properties. In this paper, we will concentrate on the situation when the molecular system is excited by two ultra-short laser pulses traveling in two different directions k1 and k2 of which the second one arrives with a delay s. Thus we have A1 ðtÞ ¼ A0 eik1 r dðt þ sÞ; A2 ðtÞ ¼ A0 eik2 r dðtÞ. The time zero is set to the center of the second pulse. This corresponds to a typical situation in the coherent nonlinear spectroscopy, and the time t in which the excited state evolves corresponds to the population time T of the non-linear spectroscopy. Inserting the delta function envelopes into Eq. (25)

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

and assuming that the pulse with wave-vector k1 precedes the pulse with wave-vector k2 , we arrive at

q^ Ið2Þ ðt; sÞ ¼ RI ðt; sÞjA0 j2 :

ð28Þ

This operator corresponds to the excited state time evolution in the left hand side (l.h.s.) diagram of Fig. 2BB. The base of this diagram corresponds to the so-called rephasing pathways. When the pulse sequence is reverted (pulse k1 arrives second with delay s) we obtain so-called non-rephasing pathways and the time evolution corresponds to the operator

^ IIð2Þ ðt;

q

2

sÞ ¼ RII ðt; sÞjA0 j :

ð29Þ ^ ð2Þ I ðt;

^ ð2Þ I ðt;

It is easy to verify that q sÞ and q sÞ do not have to be Hermitian, and they alone cannot be said to represent a state of a molecular system. This is the result of them being only a portion of the perturbation expansion of the non-linear response operator. Only the sum of Eqs. (26) and (27) yields a Hermitian operator which can describe excited state of a molecular system. This situation has an analogy in the case of 2D coherent spectroscopy where the rephasing and non-rephasing signals alone cannot be interpreted as representing absorption and emission events, while the sum spectrum can be assigned this interpretation [2]. When the system is excited by a single finite length pulse, the two contributions to the excited state are equally weighted, guaranteeing thus the proper properties of the corresponding density matrix. In the following section we will discuss the equations of motion ^ ð2Þ ^ ð2Þ for q I ðt; sÞ and qII ðt; sÞ to asses the influence of the delay s on the dynamics of their diagonal (populations) and off-diagonal (coherence) elements. We will thus attempt to answer the question whether in the non-linear spectroscopic methods, such as 2D coherent spectroscopy, we observe the expected excited state dynamics, i.e. the one unaffected by the delay s. 4.2. Simulation of the response by reduced density matrix equations In order to calculate the second order non-linear response, we need to evaluate the component expressions of Eq. (27). We start with equations of motion for the perturbation of the total density matrix WðtÞ. We define the first order operator as ð1Þ

W II ðtÞ ¼ U ðgegeÞ ðtÞV ðgeggÞ weq jgihgj:

ð30Þ

This operator corresponds to the evolution after the first interaction of the system with light in response function, Eq. (27). The response ð1Þ function has to be read from left to right. The operator W II ðtÞ satisfies the following equation

@ ð1Þ ð1Þ W ðtÞ ¼ iLðgegeÞ W II ðtÞ; @t II

ð31Þ ð1Þ

with the initial condition W II ðt ¼ 0Þ ¼ V ðgeggÞ weq jgihgj ¼ weq jgihgjlðgeÞ . RDM equation of motion can be found applying stanð1Þ dard AK projection operator, Eq. (18), for which ð1  P 0 ÞW II ðt ¼ 0Þ is equal to zero identically. The procedure of the derivation of the RDM equation of motion leads to

@ ð1Þ q ðtÞ ¼ iL0ðgegeÞ qIIð1Þ ðtÞ  R0ðgegeÞ ðtÞqð1Þ II ðtÞ; @t II

ð32Þ

ðgegeÞ

where L0 is the coherence block of the electronic Liouville superðgegeÞ operator, and R0 ðtÞ is some second order dephasing tensor. The ðgegeÞ particular form of R0 depends on the approximations and the theory applied. For a pure dephasing model and a harmonic bath, the dephasing tensor which follows from a second order CL-QME ð1Þ (see e.g. [24]) can be shown to yield an exact dynamics for qII ðtÞ [25]. In secular approximation, Eq. (32) yields

@ ð1Þ ð1Þ q ~ ðtÞ ¼ ixkg ~ q ~ ðtÞ  II;kg @t II;kg

Z 0

t

107

which is a set of independent equations for optical coherences. Note ~ to denote the excitonic representation. In the that we use index k ~ and jki coincide.The limit of jmn ! 0 (pure dephasing) the states jki evolution in the second propagation interval can be expressed through the density operator ð2Þ

ð1Þ

W II ðt; sÞ ¼ U ðeeeeÞ ðtÞV ðeeegÞ W II ðsÞ;

ð34Þ

which satisfies

@ ð2Þ ð2Þ W ðt; sÞ ¼ iLðeeeeÞ W II ðt; sÞ; @t II

ð35Þ ð2Þ

ð1Þ

with the initial condition W II ðt ¼ 0; sÞ ¼ V ðeegeÞ W II ðsÞ ð1Þ ¼ lðegÞ W II ðsÞ. The application of the AK projector would lead to a loss of information in the construction of the response function, because the term Q0 lðegÞ W ð1Þ ðsÞ ¼ ð1  P 0 ÞW ð1Þ ðsÞ – 0. Consequently, even if we successfully derive an exact master equation for ð2Þ qð2Þ II ðtÞ ¼ TrB fW II ðtÞg, the response function constructed from this equation would be missing the Q-containing terms. The parametric projection operator recently suggested in Ref. [27] has the property of eliminating the initial term approximately, ð1  P s Þleg ð1Þ W II ðsÞ 0 . In the case of pure dephasing, it turns to zero exactly. ð2Þ We can therefore write an equation of motion for qII which is analogical to Eq. (32)

@ ð2Þ ð2Þ ðeeeeÞ q ðt; sÞ ¼ iLðeeeeÞ qð2Þ ðtÞqII ðt; sÞ: s II ðt; sÞ  Rs @t II

ð36Þ

ðeeeeÞ

Again, the Ls is the electronic Liouville superoperator and RðeeeeÞ ðtÞ is the superoperator describing the electronic energy relaxs ation and the electronic coherence dephasing in the single-exciton band. In the following section, we will derive the relaxation tensor RðeeeeÞ ðtÞ corresponding to the projection operator P s . s 5. Master equations with parametric projector In this section, we will apply the well-known Nakajima–Zwanzig identity in the second order in system bath interaction Hamiltonian together with the parametric projection operator proposed in Ref. [27]. We are interested in the second interval of the non-linear response, and in the evolution in the single exciton band in particular. 5.1. The choice of projector The projector P s is chosen in such way that it contains time evolution of the bath in the first interval of the response, where the relevant system dynamics is the one of an optical coherence. The projector corresponding to the pathway RII and the length of the first interval s reads according to Ref. [27]

P s;II  ¼ trB fgU g ðsÞweq

P ey U n~ ðsÞK n~ egn~n~ ðsÞ :

ð37Þ

~ n

This choice gives an exact description of the bath for jmn ¼ 0. For a non-zero coupling it corresponds to the secular approximation in the first interval equation of motion, Eq. (33). Further on in the text, all derivations will be done for the pathway RII , and we omit the index II. The treatment of the pathway RI is analogical. We define evolution operators describing the dynamics of the environmental DOF while the system is in its electronic ground state

 i U g ðsÞ ¼ exp  ðT þ Ug Þs h

ð38Þ

~i and in the excited eigenstate jn

 i U en~ ðsÞ ¼ exp  ðT þ Uen~n~  hUen~n~  Ug iÞs : h

ð39Þ

By the choice of the projector, Eq. (37), we prescribe an ansatz ð1Þ

dsC k~k~ ðsÞqII;kg ~ ðtÞ;

ð33Þ

W ð2Þ ð0Þ ¼ W ð1ÞðIÞ ðsÞ  qð1ÞðIÞ ðsÞwðIÞ s ;

ð40Þ

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

108

where ðIÞ

ws 

P

g

K n~ U ðs

~ n

Þweq U ney ~ ð

sÞe

g n~ n~ ðsÞ

ð41Þ

:

For zero resonance coupling jmn , this is an exact prescription for the bath. For non-zero resonance coupling, it is an approximation comparable to the secular approximation. Projector, Eq. (37), can be written in short as

P s  ¼ trB fgwðIÞ s :

ð42Þ ðIÞ

It is important to note that ws is not a purely bath operator, and it does not generally commute with q. Also the interaction picture with respect to the system Hamiltonian HS denoted by ðIÞ applies to it. It stands on the right hand side (r.h.s.) of q when evaluating RII , while in RI it would stand on the l.h.s. of q. This follows from the diagrams in Fig. 2B. The full form of W ð1Þ ðsÞ is

W ð1Þ ðsÞ ¼ U g ðsÞweq 

P ey e g U n~ ðsÞK n~ eiðen~ þhUn~ U iÞs=h qð1Þ ð0Þ:

ð43Þ

~ n

From now on, the upper index (2) will be omitted in text for the sake of brevity. We can verify that the projector property P 2s ¼ P s is fulfilled

  2g n~n~ ðsÞ ~ n ~jwÞ ¼ jmih ~ n ~ jtrB U g ðsÞweq U ey P s P s ðjmih ~ ðsÞ trB fwge n



~ n ~ jegn~n~ ðsÞ trB fwge2gn~n~ ðsÞ ¼ P s ðjmih ~ n ~ jwÞ: ¼ jmih

ð44Þ

Here, we used expression for the line shape function gðtÞ in the second cumulant approximation

  egðtÞ ¼ trB U gy ðtÞU e ðtÞweq :

ð45Þ

The action of the projector P s ð P s;II Þ on the electronic state is asymmetric, because the projector was derived for the Liouville pathway II (see Fig. 2B). 5.2. Parametric master equation Now, we apply the projector P s , Eq. (42), to the Nakajima– Zwanzig identity in the interaction picture @ P s W ðIÞ ðtÞ ¼ NZ1 þ NZ2 þ NZ3 @t  Z t ¼ P s LðIÞ ðtÞ exp i ds0 Qs LðIÞ ðs0 ÞQs Qs Wðt0 Þ t0 " ! Z 0 Z tt 0



0

only the CL-QME leads in the limit of jmn ¼ 0 to a result coinciding with the one obtained by the second cumulant treatment of the non-linear response functions (see e.g. Ref. [27] and the Appendix D) Let us point out the most important aspect of the application of the projector P s with the identity, Eq. (46). It is important to note that the projector P s itself contains the system-bath interaction to all orders in the form of the exponential of the line-shape function (see Eq. (37)). The success of the second order master equations (of the form q_ ¼ a2 q, where the dot denotes the time derivative, and a2 is some second order operator) lies in the fact that their solutions includes all orders of the perturbation ðq ¼ ea2 t q0 Þ. The solution corresponds to a partial summation of the perturbative series to infinity. For some types of bath, such as the bath consisting of harmonic oscillators, this may even lead to exact master equations [25]. When higher order terms are added to the right hand side of the equation motion by a procedure that does not respect the form of higher order terms dictated by the Nakajima–Zwanzig identity, the resulting equation of motion may lead to unphysical results. Therefore, one has to take care in application of the projector P s , not to allow higher than second order contributions to appear on the right hand side of Eq. (46). Since the difference between projectors P 0 and P s is only in dynamics of the system-bath coupling during time s, their difference is at least of the first order in DU. Since the NZ2 with AK projector is already of the second order in DU in all its terms, the difference caused by using projector P s will be of higher order in DU. In the second order master equation, the term NZ2 with projector P s has to be equivalent to the form obtained from with AK projector (see e.g. [30]). Applying the following approximation qðIÞ ðt  s0 Þ qðIÞ ðtÞ in the term NZ2 we obtain it in the form the second order relaxation term of CL-QME [34]. The details of the evaluation of the term NZ3 are presented in Appendix B. As expected it yields a s-dependent term. Putting all results together and tracing over bath DOF we obtain Z P t  @ ðIÞ q ðt; sÞ ¼ ds g€mn ðsÞK m ðtÞK n ðt  sÞqðIÞ ðtÞ @t mn 0 þ g€ mn ðsÞK m ðtÞqðIÞ ðtÞK n ðt  sÞ þ g€mn ðsÞK n ðt  sÞqðIÞ ðtÞK m ðtÞ   g€ nm ðsÞqðIÞ ðtÞK n ðt  sÞK m ðtÞ  P þ K m ðtÞqðIÞ ðtÞK k~  qðIÞ ðtÞK k~ K m ðtÞ ~ mk

  g_ ~ ðt þ sÞ  g_ ~ ðtÞ : mk

s

ds0 P s LðIÞ ðtÞ exp i i

ds00 Qs LðIÞ ðs00 ÞQs

t0

Qs LðIÞ ðt  s0 ÞP s W ðIÞ ðt  s0 Þ  iP s LðIÞ ðtÞP s W ðIÞ ðtÞ: ð46Þ

We will use Eq. (46) up to the second order in LðIÞ , and we set t0 ¼ 0. The term NZ1 corresponds to so-called initial term, which we made equal to zero by the choice of the projector. The last term NZ3 corresponds to an effective Liouvillian, and it is usually a purely electronic operator. Now, with the parametric projector it contains additional terms originating from the system bath interaction. Its ðeeeeÞ purely electronic part will stand for the effective Liouvillian Ls in Eq. (36), while its additional s-depending contribution we will add to the relaxation tensor RðeeeeÞ . Finally, the term NZ2 is a starting s point for the derivation of the second order relaxation term, which has a form of a convolution between the RDM and some memory kernel. In the second order expansion of Eq. (46) and with an approximation W ðIÞ ðt  s0 Þ W ðIÞ ðtÞ, the resulting equation for P s W ðIÞ ðtÞ coincides with the second order approximation of the equivalent time-convolutionless identity [34]. Thus the NZ2 will lead to a contribution to the tensor RðeeeeÞ in Eq. (36). However, s the parametric projection operator technique allows us to keep the convolution form of the equation if it is desired. While the convolution form may have some advantages over the CL-QME [24],

mk

ð47Þ

The first three lines of Eq. (47) corresponds to the standard CL-QME, while the last two lines represent the s-dependent contribution which the standard CL-QME does not predict. Let us investigate the s-dependent term only. First, we turn to Schrödinger picture by substituting qðIÞ ðtÞ ¼ U yS ðtÞqðtÞU S ðtÞ. We denote the new s-dependent term of the CL-QME by Dðt; sÞ

Dðt; sÞqðtÞ ¼

P ~ mk

 K m ðtÞqðtÞK k~  qðtÞK k~ K m ðtÞ

  g_ mk~ ðt þ sÞ  g_ mk~ ðtÞ :

ð48Þ

It can be easily verified that the new term preserves the trace of qðtÞ, because Tr K m ðtÞqðtÞK k~  qðtÞK k~ K m ðtÞ ¼ 0. From the inspection of the term g_ mk~ ðt þ sÞ  g_ mk~ ðtÞ we can conclude that the effect of the parameter s is transient. If EGCF tends to zero on a time scale given by some bath correlation time, this term also tends to zero. The dynamics at long times is therefore not affected by the delay between the two excitation interactions. In Appendix C we derived the difference term, Eq. (48), of a special case of molecular homodimer with bath fluctuations uncorrelated between the sides. We found the s-dependent term to be identically equal to zero in this case. In the next section we study the effect of the s-dependent term numerically for a heterodimer.

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

109

Fig. 3. The time evolution of the real parts of RI and RII for a heterodimer with energy gap D ¼ 100 cm1 , interacting with a bath represented by a general Brownian oscillator with parameters k ¼ 50 cm1 ; c ¼ 1 ps1 and XBath ¼ 100 ps1 . Red lines, the matrix element 12 of RI according to standard CL-QME (dashed) and the corresponding parametric CL-QME (full). Black lines, the matrix element 12 of RII according to standard CL-QME (dashed) and the corresponding parametric CLQME (full). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

6. Numerical results and discussion In this section, we study the dynamics of the elements of the response functions, Eqs. (26) and (27). The RDM qðt; sÞ for which we derived Eq. (47) corresponds the response function RII (see Eq. (29) for the case of A0 ¼ 1). We compare the dynamics in two cases. In the first case, the evolution operator U ðeeeeÞ ðTÞ appearing in Eqs. (26) and (27) is calculated by standard time dependent CL-QME derived using AK projector, Eq. (18). In the second case, it is calcu-

Fig. 4. The time evolution of the real parts of RI and RII for a heterodimer with energy gap D ¼ 100 cm1 , interacting with a bath represented by a general Brownian oscillator with parameters k ¼ 5 cm1 ; c ¼ 1 ps1 and XBath ¼ 19 ps1 . Red lines, the matrix element 12 of RI according to standard CL-QME (dashed) and the corresponding parametric CL-QME (full). Black lines, the matrix element 12 of RII according to standard CL-QME (dashed) and the corresponding parametric CL-QME (full). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

110

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

Fig. 5. The time evolution of the real parts of RI and RII and their dependency on the resonance coupling j for a heterodimer with energy gap D ¼ 100 cm1 and an anti-parallel arrangement of the transition dipole moments, interacting with a bath represented by an overdamped Brownian oscillator with reorganization energy k ¼ 120 cm1 and correlation time sc ¼ K1 ¼ 50 fs. The delay between the two pulses is fixed to s ¼ 60 fs. Black lines, the matrix elements 11; 22 (left column) and 12 (right column) of RII according to standard CL-QME (dashed) and the corresponding parametric CL-QME (full). Red lines, the matrix element 12 (right column) of RI according to standard CL-QME (dashed) and the corresponding parametric CL-QME (full). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

111

Fig. 6. The same as Fig. 5 but for a heterodimer with parallel transition dipole moment arrangement.

lated using Eq. (47). We use the relation qI;ij ðtÞ ¼ q II;ji ðtÞ between the RDMs from different Liouville pathways. The most simple system which exhibits s-dependent correction to the standard CL-QME is a molecular heterodimer. In general, it is characterized by orientation and magnitude of its transition dipole

moments, the excited state energies e1 ; e2 of the component molecules, their resonance coupling J and the properties of the bath. We denote the difference of the excited state energies by D  e1  e2 , and we set the magnitudes of the transition dipole moments to unity. The initial condition is assumed in a form qð1Þ ð0Þ ¼ jgihgj.

112

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

The evolution during the first interval of the response follows Eq. (33). For the calculations presented on Figs. 3–5, we choose the anti-parallel orientation of the transition dipole moments, while for the calculation shown on Figs. 6 and 7, we choose the parallel orientation. The temperature is set to T = 300 K in all calculations. First, let us investigate the sensitivity of the excited state dynamics to the interplay of the delay s and the phase of the bath vibrations. Fig. 3 shows the dynamics of the electronic coherence between the excited states je1 i and je2 i in presence of the bath represented by a single-mode general Brownian oscillator with parameters k ¼ 50 cm1 ; c ¼ 1 ps1 ; XBath ¼ 100 ps1 and D ¼ 100 cm1 . Both calculations show that the standard CL-QME calculation of U ðeeeeÞ ðTÞ is insensitive to the phase of bath vibration mode during the time evolution in the first interval. The parametric CLQME, Eq. (47), shows a distinct sensitivity to this phase. In both

cases, the resonance coupling J is chosen to be zero, and the dynamics can therefore be evaluated exactly by the cumulant expansion technique (see Appendix D). The numerical evaluation using Eq. (47) indeed matches the analytical result. In the calculation on the Fig. 3, the vibration of the bath is much faster than the period (333 fs) of the electronic coherence. The two theories give the same result for s ¼ 0 fs, then they start to deviate and after one period of bath oscillator, at s ¼ 60 fs, they coincide again. Initial phase of qI;12 and qII;12 is in general different at T ¼ 0 fs because of their different time evolution in the first interval (they are not simply complex conjugates of each other). In the Fig. 4, we calculated the same system, but we used bath with parameters k ¼ 5 cm1 , c ¼ 1 ps1 ; XBath ¼ 19 ps1 . The frequency is now resonant with the frequency of the electronic coherence, which makes the effect more significant. The two theories give the same result

Fig. 7. The time evolution of the real parts of RI and RII a heterodimer with energy gap D ¼ 100 cm1 , resonance coupling j ¼ 33:1 cm1 and a parallel arrangement of the transition dipole moments, interacting with a bath represented by one overdamped Brownian oscillator with parameters K1 ¼ 100 fs and koverdamped ¼ 12 cm1 and an underdamped vibration with reorganization energy kharmonic ¼ 30 cm1 and frequency XBath ¼ 60 cm1 . Black lines, the matrix elements 11, 22 (left column) and 12 (right column) of RII according to standard CL-QME (dashed) and the corresponding parametric CL-QME (full). Red lines, the matrix element 12 (right column) of RI according to standard CL-QME (dashed) and the corresponding parametric CL-QME (full). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

q_ ðtÞ ¼ aðqðtÞ  q0 Þ þ ðf ðt þ sÞ  f ðtÞÞqðtÞ;

ð49Þ

which describes exponential decay to a limiting value q0 with the rate constant a and a modulation by the function

f ðtÞ ¼

Z

t

dsA cosðxsÞ:

ð50Þ

0

Eq. (50) represents a first integral of the EGCF of an underdamped vibrational mode, Eq. (13). The solution of Eq. (50) with the parameters a ¼ 0:01; A ¼ 104 ; x1 ¼ 60 fs and q0 ¼ 0:6 is shown in Fig. 8. We can see that the oscillation phase and period of the oscillations is indeed not proportional to s, and the amplitude increases with s similarly to Fig. 7. The behavior is therefore a direct consequence of the parametric term in the parametric QME. The overall picture arising from the numerical simulations is the following: Except for special cases, such as the molecular homodimer, the excited state dynamics of an open quantum system, as it is observed by the non-linear spectroscopy, indeed depends on the delay s between the FWM scheme. The 2D spectroscopy therefore observes a certain averaged dynamics. The effects seem to be rather small in most cases, but they might be

τ=

1

100 fs 80 fs 60 fs

0.8

40 fs 20 fs

p(t)

for s ¼ 0 fs, and at s ¼ 360 fs, which is approximately one period of the bath vibration mode. Unlike in Fig. 3, the initial phase of qI;12 and qII;12 differs significantly in T ¼ 0 fs because of their different time evolution in s. In both the cases studied above, the time evolution of the offdiagonal elements of the second order response operator is slightly modulated by the time evolution of the vibrational DOF. The phase of the oscillations seems to be mostly unaffected. Figs. 5 and 6 demonstrate the influence of the resonance coupling on the population and electronic coherence dynamics in the homodimer. As above, we perform calculation according to standard CL-QME and the parametric CL-QME. This time, we choose the overdamped Brownian oscillator with fixed s ¼ 60 fs and D ¼ 100 cm1 , reorganization energy k ¼ 120 cm1 and correlation time sc ¼ K1 ¼ 50 fs to represent the bath, and we change the resonance coupling. We calculate both diagonal and off-diagonal elements (‘‘populations’’ and ‘‘coherences’’) of the operators qI=II . For J ¼ 0 cm1 , there is no population dynamics. By increasing the coupling, the difference between the theories in the diagonal elements increases. In the Fig. 5, the dipole moments of the molecules are anti-parallel, while in Fig. 6 they are parallel. Let us now investigate a molecular dimer coupled to overdamped harmonic bath and to a single harmonic mode with frequency XBath . The harmonic mode is assumed to continue oscillating even long after thermalization in the overdamped part of the bath has taken place. Therefore, the s-dependent term of Eq. (48) changes the system dynamics also at long times. The time dependence of the density operator elements for different s is shown in Fig. 7. Parameters of the overdamped bath are K1 ¼ 100 fs and koverdamped ¼ 12 cm1 and of the harmonic mode kharmonic ¼ 30 cm1 ; XBath ¼ 60 cm1 . The dimer is characterized by D ¼ 100 cm1 ; J ¼ 33:1 cm1 and the parallel electronic transition dipole moments of the molecules. We can notice that the interaction of the vibrational mode with the electronic DOF induces oscillations in both the diagonal and off-diagonal elements of the density operator. The amplitude of the oscillations increases with increasing s. Since the EGCF of the harmonic mode, Eq. (13), is periodic, we expect the relative amplitude of the oscillations to decrease again for sufficiently long s, and to become zero for s ¼ 2p=XBath . However, in such long s the second order response would be almost zero due to the decay of the optical coherence in the first interval. Fig. 7 suggests that the phase of the oscillation does not change linearly with s. To understand this behavior, we can study a simple equation of the form

113

0.6 0 fs 0.4 0

200

400

600

800

1000

T [fs] Fig. 8. The time evolution of ‘‘population’’ according to model Eq. (49). The curves are solutions for s increasing from 0 fs to 100 fs with step of 20 fs.

observable by an advanced implementation of 2D spectroscopy. They are especially pronounced in the case of the intramolecular vibrational modes, which have frequency similar to the electronic energy gap between excitonic levels. Both the dynamics of electronic level populations and electronic coherences are affected. In order to identify these effects in the experimental data, the theory has to be extended to include also the third interval of the third order non-linear response. The corresponding projection operator P Tþs which now depends on the duration of both the coherence and the population intervals s and T, respectively, has been already proposed an tested for jmn ¼ 0 in Ref. [27]. The formulation of the theory for the third interval of the response of a multilevel excitonic system in a similar manner as performed in this paper for the second interval will be the subject of our future work.

7. Conclusions In this paper, we have demonstrated that the third and the second order non-linear response of a multilevel system cannot be completely evaluated by propagating reduced density matrix by equations of motion derived using a single projection operator. Such treatment would inevitably neglect correlations between the time evolution of the bath during the neighboring intervals of the non-linear response. We have derived equation of motion, the parametric quantum master equation, which takes these correlations into account approximately, and showed that in the absence of resonance coupling the method yields an agreement with the result obtained by the second order cumulant method. We confirm by numerical simulations that for different delays s between the excitation pulses, distinct dynamics of both excite state populations and electronic coherence occurs, in the presence of environmental degrees of freedom with finite bath correlation time and in the presence of intramolecular vibrations. The twodimensional Fourier transformed spectroscopy sees in these cases some averaged dynamics.

Acknowledgments This work was supported by the Czech Science Foundation (GACR) grant no. 205/10/0989 and the Ministry of Education, Youth and Sports of the Czech Republic via grant KONTAKT ME899 and the research plan MSM0021620835. J. O. acknowledges the support by grant GAUK 416311.

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

114

Appendix A. Third-order non-linear response functions

NZ3 ¼

The four standard response functions (so-called Liouville pathways) of a two-band electronic system read in our block formalism as

Z s

1 h 

2

ds

0

" P

0

~ mk

# ðIÞ

ðIÞ

K m ðtÞq ðtÞK k~  q ðtÞK k~ K m ðtÞ

   trB DUk~ ðs0 ÞDUm ðtÞweq wðIÞ s  P ¼ K m ðtÞqðIÞ ðtÞK k~  qðIÞ ðtÞK k~ K m ðtÞ ~ mk

  g_ mk~ ðt þ sÞ  g_ mk~ ðtÞ wðIÞ s :

R1 ðt; T; sÞ ¼ trflðgeÞ U ðegegÞ ðtÞV ðegeeÞ ðggÞ  U ðeeeeÞ ðTÞV ðeeegÞ U ðegegÞ ðsÞV ðegggÞ W eq g;

R2 ðt; T; sÞ ¼ trfl U

ðgeÞ

U

ðeeeeÞ

ðegegÞ

ðTÞV

ðtÞV

ðeegeÞ

ðA1Þ

ðegeeÞ

ðggÞ U ðgegeÞ ðsÞV ðgeggÞ W eq g;

ðA2Þ

In order to obtain the last two lines of Eq. (47) we will trace Eq. (B7) over the bath DOF. Appendix C. The s-dependent term for a homodimer

R3 ðt; T; sÞ ¼ trflðgeÞ U ðegegÞ ðtÞV ðegggÞ

Let us assume a molecular homodimer with the Hamiltonian

 U ðggggÞ ðTÞV ðgggeÞ U ðgegeÞ ðsÞV ðgeggÞ W ðggÞ eq g; R4 ðt; T; sÞ ¼ trfl U

ðgeÞ

U

ðggggÞ

ðegegÞ

ðTÞV

ðtÞV

ðegggÞ

ðggegÞ

ðggÞ U ðegegÞ ðsÞV ðegggÞ W eq g:

ðA3Þ

ðA4Þ

Here, W ðggÞ eq ¼ weq jgihgj where weq is the bath equilibrium density operator.

HS ¼

Appendix B. Derivation of relaxation terms



P DUm ðtÞK m ðtÞ; qðIÞ ðtÞwðIÞ s m

 wðIÞ s

ðB1Þ

n

o

qðIÞ ðtÞ ¼ trB P s W ðIÞ ðtÞ :

ðB2Þ

o i P g ~ ðsÞ n ðIÞ NZ3 ¼  e k trB DUm ðtÞK m ðtÞqðIÞ ðtÞweq  U g ðsÞU ey ~ ws ~ ðsÞK k k h mk~ o i P g ~ ðsÞ n ðIÞ ðIÞ þ e k trB q ðtÞweq: U g ðsÞU key ~ DUm ðtÞK m ðtÞ ws : ~ ð sÞ  K k h mk~ ðB3Þ g ðsÞ

Now we have to set e k~ 1 since its contribution is of higher order in DU. Applying the first order expansion

i h

Z s

ds0 DUn~n~ ðs0 Þ;

ðB4Þ

0

yields    Z iP i s 0 trB DUm ðtÞK m ðtÞqðIÞ ðtÞweq  1 þ ds DUk~k~ ðs0 Þ K k~ wðIÞ s h mk~ h 0   Z s iP i þ trB 1þ ds0 DUk~k~ ðs0 Þ  K k~ DUm ðtÞK m ðtÞqðIÞ ðtÞweq wðIÞ s : h mk~ h 0

NZ3 ¼ 

ðB5Þ

Now, we define a line shape function

g ab~ ðtÞ ¼

1 h

2

Z

t

ds 0

Z s 0

j 0 DU1

ðC1Þ

;

0 : DU2

0

ðC2Þ



j 0

; 0 j 1 DU1 þ DU2 ¼ 2 DU2  DU1

ðC3Þ

DU2  DU1 : DU1 þ DU2

ðC4Þ

g_ mk~ ðtÞ ¼

1 2

Z

t

dt 0

0

2 P

hDUm ðt 0 ÞDUn i:

ðC5Þ

n¼1

We neglect the cross-terms hDUi ðtÞDUj i for i – j , because we assume only local correlations of the bath. In a homodimer _ DU1 ¼ DU2 , and we get g_ 1k~ ðtÞ ¼ g_ 2k~ ðtÞ ¼ gðtÞ from which it follows that

g_ mk~ ðt þ sÞ  g_ mk~ ðtÞ ¼ g_ ðt þ sÞ  g_ ðtÞ:

Using Eq. (41) yields

U g ðsÞU ey ~ ðsÞ 1 þ n

j

The derivative of the line shape function which enters the s-dependent correction Dðt; sÞ of the CL-QME then reads as



where

0

This yields in the exciton basis

HSB

In this section we will evaluate the term NZ3 of Eq. (46). Applying the definitions of its component operators from Section 5 we obtain



HSB ¼

HS ¼

i NZ3 ¼  trB h

ðB7Þ

s0 hDUa ðs0 ÞDUb~ i;

ðC6Þ

By substituting this equation into Eq. (48), the s-containing term Dðt; sÞ turns to zero, because

P ðK m qðtÞK k~  qðtÞK k~ K m Þ ¼ 0:

ðC7Þ

~ mk

Thus the s-dependent correction to the CL-QME is zero for the case of a homodimer. Appendix D. Pure dephasing of an electronic coherence In a pure dephasing model of the system-bath interaction, analytical solution for the dynamics of electronic coherences can be obtained with the cumulant expansion technique. The time evoluð2Þ tion of coherence qee0 ðt; sÞ reads as





ð2Þ ð2Þ y y qee 0 ðt; sÞ ¼ trB U e ðtÞU g ðsÞW eq U e0 ðsÞU e0 ðtÞ qee0 ð0Þ;

ðD1Þ

where we set all transition dipole moments to one. By using the second cumulant procedure [3] we can evaluate the expression into





ð2Þ g ee ðtÞg 0 0 ðtþsÞþg e0 e ðtÞ þg 0 ðtþsÞg 0 ðsÞ ð2Þ ee ee qee e ee qee0 ð0Þ: 0 ðt; sÞ ¼ e

ðD2Þ

Equivalently, we can rewrite Eq. (D2) in form of a master equation

ðB6Þ

where it is noteworthy that one of its indices represents the siteand the other the exciton-basis. The NZ3 term can then be expressed as

   d ð2Þ q 0 ðt; sÞ ¼ g_ e0 e0 ðtÞ  g_ ee ðtÞ þ g_ e0 e ðtÞ þ g_ ee0 ðtÞ  g_ e0 e0 ðt þ sÞ  g_ e0 e0 ðtÞ dt ee   ðD3Þ þ g_ ee0 ðt þ sÞ  g_ ee0 ðtÞ qð2Þ ee0 ðt; sÞ:

The first four terms on the right hand side of Eq. (D3) are described by the CL-QME. The remaining terms can be matched with the term

J. Olšina, T. Mancˇal / Chemical Physics 404 (2012) 103–115

Dðt; sÞ derived in this paper. Eq. (47) therefore gives the correct result for this simple exactly solvable case. References [1] S. Mukamel, Ann. Rev. Phys. Chem. 51 (2000) 691. [2] D.M. Jonas, Ann. Rev. Phys. Chem. 54 (2003) 425. [3] S. Mukamel, Principles of Nonlinear Spectroscopy, Oxford University Press, Oxford, 1995. [4] M.H. Cho, H.M. Vaswani, T. Brixner, J. Stenger, G.R. Fleming, J. Phys. Chem. B 109 (2005) 10542. [5] A.V. Pisliakov, T. Mancˇal, G.R. Fleming, J. Chem. Phys. 124 (2006) 234505. [6] P. Kjellberg, B. Bruggemann, T. Pullerits, Phys. Rev. B 74 (2006) 024303. [7] G.S. Engel, T.R. Calhoun, E.L. Read, T.-K. Ahn, T. Mancˇal, Y.-C. Cheng, R.E. Blankenship, G.R. Fleming, Nature 446 (2007) 782. [8] A. Nemeth, F. Milota, T. Mancˇal, V. Lukeš, H.F. Kauffmann, J. Sperling, Chem. Phys. Lett. 459 (2008) 94. [9] A. Nemeth, F. Milota, T. Mancˇal, V. Lukeš, J. Hauer, H.F. Kauffmann, J. Sperling, J. Chem. Phys. 132 (2010) 184514. [10] H. van Amerongen, L. Valkunas, R. van Grondelle, Photosynthetic Excitons, World Scientific, Singapore, 2000. [11] Y. Tanimura, J. Phys. Soc. Jpn. 75 (2006) 082001. [12] A. Ishizaki, G.R. Fleming, J. Chem. Phys. 130 (2009) 234111. [13] A. Ishizaki, G.R. Fleming, J. Chem. Phys. 130 (2009) 234110. [14] V. Novoderezhkin, M. Wendling, R. van Grondelle, J. Phys. Chem. B 107 (2003) 11534. [15] V.I. Novoderezhkin, M.A. Palacios, H. van Amerongen, J. Phys. Chem. B 109 (2005) 10493.

115

[16] R. van Grondelle, V.I. Novoderezhkin, Phys. Chem. Chem. Phys. 8 (2006) 793. [17] D. Abramavicˇius, B. Palmieri, D.V. Voronine, F. Šanda, S. Mukamel, Chem. Rev. 109 (2009) 2350. [18] D. Zigmantas, E.L. Read, T. Mancˇal, T. Brixner, A.T. Gardiner, R.J. Cogdell, G.R. Fleming, P. Natl. Acad. Sci. USA 103 (2006) 12672. [19] E.L. Read, G.S. Engel, T.R. Calhoun, T. Mancˇal, T.K. Ahn, R.E. Blankenship, G.R. Fleming, P. Natl. Acad. Sci. USA 104 (2007) 14203. [20] H. Rhee, Y.-G. June, J.-S. Lee, K.-K. Lee, J.-H. Ha, Z.H. Kim, S.-J. Jeon, M. Cho, Nature 548 (2009) 310. [21] D. Abramavicius, S. Mukamel, J. Chem. Phys. 133 (2010) 064510. [22] W.M. Zhang, T. Meier, V. Chernyak, S. Mukamel, J. Chem. Phys. 108 (1998) 7763. [23] B. Palmieri, D. Abramavicius, S. Mukamel, J. Mol. Model 130 (2010) 204512. [24] J. Olšina, T. Mancˇal, J. Mol. Model. 110 (2010) 23456. [25] R. Doll, D. Zueco, M. Wubs, S. Kohler, P. Hanggi, Chem. Phys. 347 (2008) 243. [26] M. Richter, A. Knorr, Ann. Phys. 325 (2010) 711. [27] T. Mancˇal, F. Šanda, Chem. Phys. Lett. 530 (2012) 140. [28] T. Mancˇal, L. Valkunas, New J. Phys. 12 (2010) 065044. [29] Chlorophylls and bacteriochlorophylls, B. Grimm, R.J. Porra, W. Rüdiger, H. Scheer (eds.), (Springer, Dordrecht, 2006). [30] V. May, O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, Wiley-VCH, Berlin, 2001. [31] P.N. Argyres, P.L. Kelley, Phys. Rev. A 134 (1964) 98. [32] M.F. Gelin, D. Egorova, W. Domcke, J. Chem. Phys. 123 (2005) 164112. [33] T. Mancˇal, A. Pisliakov, G. Fleming, J. Chem. Phys. 124 (2006) 234504. [34] B. Fain, Irreversibilities in Quantum Mechanics, Kluwer Academic Publishers, Dordrecht, 2000.