PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Phys Chem B. Author manuscript; available in PMC 2010 July 11.
Published in final edited form as:
PMCID: PMC2901710
NIHMSID: NIHMS86524

Stochastic Liouville Equations for Coherent Multidimensional Spectroscopy of Excitons

Abstract

Signatures of chemical exchange and spectral diffusion in 2D photon-echo lineshapes of molecular aggregates are studied using model calculations for a dimer whose Hamiltonian parameters are stochastically modulated. Cross peaks induced by chemical exchange and by exciton transport has different dynamics, and and distinguish two models which have the same absorption spectrum (a two-state-jump bath modulation model of a dimer and a four state jump bath model of a single chromophore). Slow Gaussian-Markovian spectral diffusion of a symmetric dimer induces new peaks which are damped as the dipole moment is equilibrated. These effects require an explicit treatment of the bath and may not be described by lower-level theories such as the Redfield equations, which eliminate the bath.

1 Introduction

Exciton models are widely used in the description of the coherent nonlinear optical response of complex systems.1 Vibrational spectra of proteins,26 liquid water,7 and electronic spectra of photosynthetic antenae8 are a few examples. Bath-induced fluctuations of transition frequencies and couplings affect the electronic and transport properties of aggregates (e.g. photosynthetic antenae).

A broad arsenal of techniques has been developed towards the simulation of two dimensional coherent spectroscopy (2DCS) signals.1 Direct simulation of the optical response requires repeated diagonalizations of the system Hamiltonian at many time points. The response is then given by a path integral over the trajectories of the time dependent eigenvalues and the overlap factors of the instantaneous eigenstates at different times.9,10

The fluctuating time-dependent Hamiltonian H(t) defines the instantaneous (adiabatic) basis |[var phi]j(t)right angle bracket and energies εj(t) (we set ћ = 1)

H(t)|φj(t)=εj(t)|φj(t)

Direct evaluation of the time evolution operator, when the time t is divided into N → ∞ short segments Δt [equivalent] t/N gives

UTexp[i0tdτH(τ)]=Tn=1Nexp[iH(tn)Δt]    tn=nΔt
(1)

Adopting the matrix representation Ũjk [equivalent] left angle bracket[var phi]j(t)|U|[var phi]k(0)right angle bracket we get

Ũ=Tn=1NUnSn;
(2)

Here

(Un)jj=exp[iεj(tn)Δt]δjj

is the adiabatic propagator, and the nonadiabatic coupling matrix

(Sn)jj=ϕj(tn)|ϕj(tn1),

represents the overlap of the adiabatic states at two adjacent points.

Implementing Eq. (1) is numerically expensive. Moreover, one must typically account for various fluctuation scenarios and average over sufficiently long trajectories (or, equivalently by ergodic hypothesis, over possible paths) of the fluctuating Hamiltonian.

More affordable simulations are feasible when the fluctuations are small and fast. The response can then be calculated by sums over delocalized eigenstates of the average Hamiltonian, where dephasing rates are introduced phenomenologically.11 Exciton transport may be accounted for by using the Redfield equations1,12,13 for the reduced Liouville space dynamics of the chromophores. A quasiparticle approach based on the nonlinear exciton equations allows the interpretation of signals in terms of exciton scattering.1 Very slow fluctuations can be readily incorporated by static averaging of the signals over their realizations.14

The stochastic Liouville equations (SLE)1517 allow the modelling of strong fluctuations with arbitrary time scales provided they can be described by a few classical collective coordinates, which satisfy a Markovian master equation. The path integrals of Eq.(2) and the final averaging are avoided by calculating the density distributions in the joint system + bath space. The SLE generalize the Redfield equations, that are limited to fast fluctuations where the bath degrees of freedom can be eliminated. We have recently applied the SLE to calculate 2D lineshapes of a single chromophore whose transition frequency undergoes a stochastic modulation due to hydrogen bonding in water18,19 or organic solvents.20,21 Excitonic spectra of a multichromophoric system (trialanine) were described by the SLE approach for a limited class of models and delay times.10

In this paper we apply the SLE to a model dimer interacting with a classical bath. We predict the complete delay time evolution of 2D signals, and show new effects of interchromophore coupling, as well as spectral diffusion characteristics of the finite fluctuations time. Two types of stochastic models are considered: discrete chemical-exchange will be described by a two state jump model for the bath, and continuous spectral diffusion accounted for using a Gaussian-Markovian coordinate (high temperature limit of a Brownian overdamped oscillator). We show how the 2D signals may distinguish the model from a non excitonic system which has the same linear response. Some effects which are missed in the simpler Redfield treatment are pointed out.

2 Stochastic Liouville equations for molecular dimer

We consider a dimer with two one-exciton states |e1=B1|g, and |e2=B2|g, and one doubly excited state |f1B1B2|g, described by the Frenkel-exciton Hamiltonian.1,2224

H=ε1(t)B1B1+ε2(t)B2B2+J(t)(B2B1+B1B2)+ΔB1B2B1B2
(3)

Both the coupling J and energies εj undergo stochastic fluctuations. We will assume that the doubly-excited state is strongly shifted Δ → ∞ so that it does not contribute to the optical signals near the single exciton frequency ω ≈ ε.

A discrete fluctuation model is introduced by a additional bath degree of freedom with two states u (up) and d (down), representing e.g. two conformers of the molecule or solvent configurations. We assume that the Hamiltonian parameters jump stochastically between two values Ju, ε1u, ε2u and Jd, ε1d, ε2d corresponding to the two configurations. The transitions between these states are described by a Markovian Pauli master equation25 with down (kd) and up (ku) jump rates, respectively.

(dρudt)M=kdρu+kuρd(dρddt)M=kdρukuρd
(4)

The entire density matrix, defined in the joint system + bath space, is described by the stochastic Liouville equation.15,16

dρdt=i[H,ρ]+(dρdt)Mρ

ρ carries three indices representing bra, and ket variable of the system Liouville space and one index denoting the state of bath.

Since the Hamiltonian (Eq.(3)) conserves the number of excitons, the Liouville operator L is block-diagonal, with separate blocks corresponding to the gg, ee, ge, and eg manifolds. We define the basis set for ground state evolution (|gg; uright angle bracketright angle bracket, |gg; dright angle bracketright angle bracket), with the block

gg,gg=(kdkukdku)

The coherence (eg) evolution is described in the following basis set |e1g; uright angle bracketright angle bracket, |e1g; dright angle bracketright angle bracket, |e2g; uright angle bracketright angle bracket, and |e2g; dright angle bracketright angle bracket by the Liouvillean

eg,eg=(kdiε1ukuiJu0kdkuiε1d0iJdiJu0kdiε2uku0iJdkdkuiε2d)

Finally, the excited state block |e1e1; uright angle bracketright angle bracket, |e1e1; dright angle bracketright angle bracket, |e1e2; uright angle bracketright angle bracket, |e1e2; dright angle bracketright angle bracket, |e2e1; uright angle bracketright angle bracket, |e2e1; dright angle bracketright angle bracket, |e2e2; uright angle bracketright angle bracket, |e2e2; dright angle bracketright angle bracket is represented by

ee,ee=(kdkuiJu0iJu000kdku0iJd0iJd00iJu0kdiΔεuku00iJu00iJdkdkuiΔεd000iJdiJu000kd+iΔεukuiJu00iJd00kdku+iΔεd0iJd00iJu0iJu0kdku000iJd0iJdkdku)

where we have denoted Δεd [equivalent] ε1d−ε2d, Δεu [equivalent] ε1u−ε2u. With ge,ge=eg,eg* the description is complete.

The Green’s function solution of the Liouville equation is

𝒢(t)θ(t) exp t

In the frequency-domain G(ω) = ∫ eiωtG(t)dt we have

𝒢(ω)=[iω+]1
(5)

Like L, the Green’s function is block diagonal as well. For instance, the time domain Green’s function for the gg block reads

𝒢gg,gg(t)=1kd+ku(kukukdkd)+e(kd+ku)tkd+ku(kdkukdku)

In the other manifolds the time domain Green’s functions will be calculated numerically by spectral decomposition of the Liouville operator. In the frequency domain, all the block matrix inversions (Eq. 5) can be performed analytically in closed, though somewhat lengthy, form.

Averaging over bath variable is made with respect to the equilibrium bath density

|ρ(eq)=1kd+ku(kukd)

and the summation over final states is represented by the scalar product with the vector left angle bracketI| = (1, 1).

3 The two- dimensional signals

We shall calculate in the response of the system to three laser pulses. The interaction is given by Hint = −E(t)Dint where

Dint=[µ1(B1+B1)+µ2(B2+B2)]

is the dipole operator. We consider impulsive experiments performed with short laser pulses. The linear response function is

S(t1)=iµ|𝒢(t1)µ()|ρ(0)

and the third-order response function26

S(t3,t2,t1)=i3µ|𝒢(t3)µ()𝒢(t2)µ()𝒢(t1)µ()|ρ(0)

where µ(−) [equivalent] [Dint, …] is the dipole superoperator, and left angle bracketleft angle bracketµ|…= left angle bracketI|TrDint… provides the mean value of dipole moment. Initially the system is in the ground state and the bath is at equilibrium

|ρ(0)=|gg|ρ(eq)

The dipole moment is independent of the bath variables so that µ(−) is represented by the diagonal matrices in bath space

û1=(µ100µ1);µ^2=(µ200µ2).

In our basis set, µ(−) is represented by the following matrices

µeg,gg()=(µ^1µ^2)µee,gg()=(µ^10µ^200µ^10µ^2)µee,ge()=(µ^100µ^1µ^200µ^2)

We further have µge,gg()=µeg,gg(),µgg,eg()=[µeg,gg()]T,µgg,ge()=[µge,gg()]T,µeg,ee()=[µee,eg()]T,µee,ge()=[µge,ee()]T. We note that µα,β=[µβ,α]T, and µij,kl()=µji,lk()forij(kl)ee

The linear response is given by

S(t1)=iI|µgg,eg𝒢eg,eg(t1)µeg,gg()|ρ(eq)+c.c.

In frequency space S1) [equivalent]eiω1t1S(t1). These represent a symmetric peak structure around ε and −ε, respectively (S(ω) = S*(−ω)). We focus on the ω1 ≈ ε peaks. The absorption lineshape is given by

I(ω1)=ImS(ω1)=ReI|µgg,eg𝒢eg,eg(ω1)µeg,gg()|ρ(eq)
(6)

We next turn to the third order resonant signals generated in the kI = −k1 + k2 + k3 and kII = k1k2 + k3 phase-matching directions. These will be displayed as ω1, ω3 frequency-frequency plots S3, t2, ω1) for various time slices t2. As in the linear response (Eq. (6)), each contribution to the third order response has a complex conjugated counterpart guving peaks at opposite frequencies S3, t2, ω1) = S*(−ω3, t2, −ω1).27 We focus on the peak structure around ω3 = ε; The corresponding response function for the kI signal is given by:

SI(ω3,t2,ω1)=(i)3[R1(ω3,t2,ω1)+R2(ω3,t2,ω1),]

and for kII

SII(ω3,t2,ω1)=(i)3[R3(ω3,t2,ω1)+R4(ω3,t2,ω1)].
(7)

The various contributions R (Liouville space pathways) are represented by the 4 Feynman diagrams shown in Fig 1:

R1(ω3,t2,ω1)=I|µgg,eg𝒢eg,eg(ω3)µeg,gg()𝒢gg,gg(t2)µgg,ge()𝒢ge,ge(ω1)µge,gg()|ρ(eq)R2(ω3,t2,ω1)=I|µgg,eg𝒢eg,eg(ω3)µeg,ee()𝒢ee,ee(t2)µee,ge()𝒢ge,ge(ω1)µge,gg()|ρ(eq)R3(ω3,t2,ω1)=I|µgg,eg𝒢eg,eg(ω3)µeg,gg()𝒢gg,gg(t2)µgg,eg()𝒢eg,eg(ω1)µeg,gg()|ρ(eq)R4(ω3,t2,ω1)=I|µgg,eg𝒢eg,eg(ω3)µeg,ee()𝒢ee,ee(t2)µee,eg()𝒢eg,eg(ω1)µeg,gg()|ρ(eq)
(8)

Fig 1
Double sided Feynman diagrams for the Liouville space pathways contributing to the third order response in the phase matching direction kI = −k1 + k2 + k3 and kII = k1k2 + k3. (Eq.7)

We shall further display the following combination of signals which provide a clean absorptive signal28

SA(ω3,t2,ω1)ImSI(ω3,t2,ω1)ImSII(ω3,t2,ω1)
(9)

4 Slow fluctuations

The absorption spectrum (Fig 2) has four peaks, centered at two E1(2),u=(ε1u+ε2u)/2±Δεu2+Ju2 eigenvalues of Hu and two E1(2),d=(ε1d+ε2d)/2±Δεd2+Jd2 eigenvalues of Hd, in the static limiting case ku, kd << |E1(2),uE1(2),d|.

I(ω)=j=12υ=u;dMj,υ2Γυρggυ(eq)Γυ2+(ωEj,υ)2
(10)

where

Mj,υ=Ej,υ|µ|g

is the transition dipole, and Γd = ku, Γu = kv are the line broadening parameters.

Fig 2
Absorption lineshapes (Eq. (6)) of a dimer with two state bath fluctuations (Eqs. (3) and (4)) (solid line) and a monomer with four state bath fluctuations (Eqs. 18) and (19)) (dotted line) in static limit. The dimer peaks correspond to resonance with ...

The 2D SA lineshape at short t2, displayed in Fig 3, top left panel, shows primarily four diagonal peaks and four selected cross peaks. Three consecutive interactions with laser pulses (applications of µ(−)) for times t2kd,u1 induce cross-peaks between peaks belonging to the same bath state, but different excitons.

SA(ω3,t2kd,u1,ω1)=Rek=12j=12υ=u;dMkυ2Mjυ2ρggυ(eq)1Γυi(ω3Ej,υ)×[2ΓυΓυ2+(ω1Ek;υ)2+eΓυ,kjet2i(Ej,υEk;υ)t2(1Γυ+i(ω1Ej,υ)+1Γυ+i(ω1Ek,υ))]
(11)

where exciton dephasing Γυ,kje is negligible in the slow limit (compare fast limit).

Fig 3
The SA signal (Eq. 9) for a dimer with two state bath (top panels) and monomer with four state bath (bottom panels) corresponding to Fig 2. For delay times Γt2 = 0, 0.1, 0.2, 1.0, 100 increasing from left to right panels. Other parameters same ...

Additional cross-peaks induced by jumps among the u and d peaks (top left panel) appear for t2kd,u1. The dynamics is quite rich on these timescales, since the phase relation between u and d eigenstates left angle bracketEj,u|Ei,dright angle bracket is important for ee evolution and may also induce fluctuations of the peaks magnitudes at 2π(Ej,vEi;v)−1 period. We distinguish between two cases: If the u and d eigenstates are degenerate |Ej,uright angle bracket = |Ej,dright angle bracket no exciton transport is allowed. This is shown in the top panels of Fig 3 where the delay time t2 increases from the left to the right panel. The peaks gradually develop at all possible combinations of frequencies. However, peaks belonging to the same exciton are stronger even for very long times. When |Ej,uright angle bracket ≠ |Ej,dright angle bracket, exciton transfer is allowed and asymptotically we have

𝒢ĒjĒkυ,ĒlĒmw(t2)=12δjkδlmkυku+kd

The absorptive lineshape can then be factorized as follows

ћSA(ω3,t2,ω1)=3I(ω1)I(ω3).
(12)

This is similar to the asymptotic lineshape of the monomer (Eq.(21)).29,30 Eq. (12) holds for arbitrary timescale and for any type of stochastic modulation (once transport is allowed). Eq. (12) can be generalized beyond the stochastic model31 to include finite temperature effects, where the emission lineshape is red shifted with respect to the absorption (the Stokes shift).

Lineshapes of this simple spectra are qualitatively reminiscent of 2D NMR spectra33 measured on systems obeying similar dynamical rules. The general relation between 2D NMR methods and their optical counterpart has been discussed in Ref.32 2D NMR works on different timescales (with differences in models and approximations used), has no phase matching, but has a better control of the pulses. These differences and the rather weak fields applied in most optical experiments make the nonlinear response functions,26 such as used here, particularly adequate for interpreting 2D optical and infrared experiments.

For a more detailed analysis of these lineshapes we display the four contributions Ri to SA in Figure 4. At t2 = 0 R1 = R2 (Fig 4, top panels) since these diagrams only differ by the second interval evolution in the ee or gg manifolds. The ket and the bra action of the dipole moment connects different exciton states in t1 and t3 intervals resulting in their cross peaks. A similar peak-pattern is seen in R3, but this nonrephasing contribution is elongated in the antidiagonal direction compared to the diagonal elongation of the rephasing contribution. R4 only shows diagonal peaks, because the eg state in the third interval is caused by interaction with the first pulse and represents the same exciton state.

Fig 4
The R13, t2, −ω1), R23, t2, −ω1), R33, t2, ω1), R43, t2, ω1) (Eq.(8), real parts) contributions to SA signal (from left to right, then top to bottom) for the dimer ...

Fig 5 repeats these calculations in the long delay tuime regime kut2, kdt2 >> 1. The ground state t2 evolution for R1 and R3 mixes the u and d bath states, and their combination shows a complete peak structure. When the rephasing R1 and nonrephasing R3 contributions, characterized by the ground state (gg) evolution during t2, are summed over, the lineshape may be recast as the the product of the two lineshapes ~ 2I1)I3).

Fig 5
The same as Fig 4., but for long delay times Γt2 = 100, corresponding to the most right top panel of Fig 3.

The ee evolution in the R2 and R4 diagrams is more complex. During t2 coherences decay creating a statistical mixture of diagonal Ejv states. Since our parametrization does not allow exciton transport, we see decoherence. All contributions with different acting dipole moments at the first and second pulse decay. The only resulting cross peaks are connected with the u and d jumps. This explains why the peaks are stronger in the SA lineshape. The features shown in Fig 5 are not universal. When exciton transport is allowed the other cross-peaks appear, approaching Eq. (12).

The complete peak structure described above is simplified when the dipole moment does not couple to all of the |Ej,vright angle bracket eigenstates, or when some levels are accidently degenerate. The 2D lineshapes can distinguish between various dynamical models with a similar absorption (1D) spectrum. To illustrate this we compare the dimer with a two state bath with the four state jump bath model (FSJ) modulating single chromophore, given in Appendix B. Comparison of Eq. (10) with Eq. (20) or the lineshapes displayed in Fig 2 shows that the two models may have identical absorption lineshape in the static limit with the proper choice of the resonance frequencies and linewidth (as defined by Eqs. (10) and (20)).

2D lineshapes of the two models are markedly different at short times as shown in Fig 3. The slow fluctuation FSJ limit (bottom left) SA lineshape at short t2 jumps only reveals the diagonal peaks, because no change of bath state can occur in contrast to dimer. With increasing delay time the complete peak structure is developed, in contrast to the dimer, peak oscillations due to coherence dynamics (Eq. (11))34 can not be observed. Asymptotically these two lineshape agree (compare Eqs.(12), (21), and (22)), except when exciton transfer is not allowed.

5 Fast and intermediate fluctuations timescale

In the motional narrowing ku, kd >> |EjuEkd| limit fluctuations are averaged on the timescale of optical experiments. The averaged Hamiltonian

=ε̄1B1B1+ε̄2B2B2+(B2B1+B1B2)
(13)

where

=kuJu+kdJdku+kd,ε̄j=kuεju+kdεjdku+kd,

may be used.

The dynamics is however different from the isolated dimer, which follows the Hamiltonian Eq. (13) since the fluctuations lead to an additional population relaxation rate (the exciton transfer) ΓT ~ left angle bracket|left angle bracketĒ1|(H[H with macron])right angle bracket|Ē2right angle bracket|2right angle bracket/k and dephasing rates ΓeΓT+i=12Ēi|(H)2|Ēi/k and dephasing between ground and excited states Γi ~ ΓT/2 + left angle bracketĒi|left angle bracket(H[H with macron])2right angle bracket|Ēiright angle bracket/k + left angle bracketg|left angle bracket(H[H with macron])2right angle bracket|gright angle bracket/k. The averaged density matrix is approximately described by the Redfield master equations in the secular approximation in the high temperature limit

dρĒ1Ē1dt=ΓTρĒ1Ē1+ΓTρĒ2Ē2dρĒ2Ē2dt=ΓTρĒ2Ē2+ΓTρĒ1Ē1dρĒ1Ē2dt=[i(E1E2)Γe]ρĒ1Ē2dρĒjgdt=[iĒjΓj]ρĒjg

The dynamics may be described in the reduced system Liouville space. Signatures of the fluctuations are present in the relaxation and dephasing rates, which may not be described in the Hilbert space. The corresponding Green’s functions in system Liouville space are

𝒢Ē1Ē1,Ē1Ē1(t)=𝒢Ē2Ē2,Ē2Ē2(t)=1+e2ΓTt2𝒢Ē1Ē1,Ē2Ē2(t)=𝒢Ē2Ē2,Ē1Ē1(t)=1e2ΓTt2𝒢Ē1Ē2,Ē1Ē2(t)=𝒢Ē2Ē1,Ē2Ē1*(t)=e[i(Ē1Ē2)Γe]t𝒢Ējg,Ējg(t)=𝒢gĒj,gĒj*(t)=e[iĒjΓj]t

Note that the SLE provides an exact representation of stochastic Markovian fluctuations of Hamiltonian dynamics. Both Liouville equation, and the averaging over bath paths preserve the positive semidefinitness of the density matrix, so does the SLE. The positivity issue of master equations,35 which is a major problem for the Redfield equations when the secular approximation is not invoked, is avoided.

The absorption spectrum has two peaks centered at the eigenvalues Ē1 and Ē2 of the averaged dimer Hamiltonian

I(ω)=i=12Mi2Γi(Γi)2+(ωĒi)2

If Γ1 + Γ2 > Ē1Ē2 these can eventually overlap and merge into a single peak. The 2D lineshape (Fig 6, bottom line) has two cross peaks induced by the annihilation of one excitation by the second pulse and the creation of different excitation by the third pulse also for short t2.

SA(ω3,t2,ω1)=i=12j=122Mi2Mj2Γi(Γi)2+(ω1Ēi)2Γj(Γj)2+(ω3Ēj)2+Rei,j,k,l=12MiMjMkMlGĒlĒk,ĒjĒi(t2)Γli(ω3Ēl)×[1Γji(ω1Ēj)+1Γi+i(ω1Ēi)]

Fig 6Fig 6
Fig 6A Top: The SA (Eq.9) signal for dimer with two state bath (Eqs. (3) and (4)) for intermediate (top) and fast (bottom) fluctuations at short t2 = 0 (left panels) and long t2=2000EΔ1 (right panels) delay times. Measured in the peak ...

This is different from two state jump model lineshapes in the slow modulation limit.20 The cross peaks have intensities comparable with the diagonal peaks at long times.

The limiting cases of slow or fast fluctuations can be treated by a simpler methods, the stochastic approach is however essential for the intermediate regime.

Fig 6 compares the intermediate (top) and fast (bottom) fluctuation limit. The separated peaks of the static case overlap and merge into a single peak. In the intermediate regime, however, the memory of two bath states can still be tracked by the diagonal elongation of peaks at short times. In the fast limit the peaks become narrower and all signatures of memory are erased. Panel B shows the corresponding absorption lineshapes.

At longer delay times (right panels) the decoherence and exciton transport erase all memory and lineshapes approach Eq. (12). The t2 dependence follows from the ee t2 evolution of R2 and R4, while the R1 and R3 contributions give half of the product (Eq.21). (not shown)

The intermediate timescale shows additional delicate effects caused by altering the eigen-vectors (rather then merely eigenvalues) of Hamiltonian. These will be discussed in the next section using a Gaussian model for spectral diffusion.

We next examine more closely the R3, R4 contributions (Fig 7). At short times GĒlĒk,ĒjĒi = δljδik, and R3 follows the four peak structure identical to R1, while R4 is diagonal, from the same reason mentioned in Fig 4, At long times t2 coherences are dumped. Asymptotically 𝒢ĒlĒk,ĒjĒi=12δlkδij and the lineshape (Eq.12) is restored. This argument holds as long as exciton transfer is allowed δJ ≠ 0, otherwise the picture of Fig 5 is repeated and the cross peaks never show up in R4.

Fig 7
The R3(top panel) and R4 (bottom panel) signal (Eq.(8, real parts)) for the dimer in fast modulation limit for short t2 = 0 (left panels) and long t2=2000EΔ1 (right panels) delay times. Other parameters as in Fig 6A, bottom.

6 Spectral diffusion

We now turn to a different model representing the coupling to uncorrelated continuous diagonal frequency fluctuations. This situation, very common in molecular aggregates, is described by adding the following coupling term

HSQ=f1(Q1(t))B1B1+f2(Q2(t))B2B2

to the Hamiltonian, together with a Fokker-Planck master equation for the density evolution of the two collective coordinates Q1 and Q2.

(dρ(Q1,Q2)dt)Q=w=12ΛwQw(Qw+σw2Qw)ρ(Q1,Q2)=Qρ
(14)

This model assumes, that bath motions primarily alter (and in simple way) the local site frequencies. The resulting delocalized eigenenergy changes are secondary and more complicated. Standard simulations of aggregates rely on diagonalization of averaged Hamiltonian allowing arbitrary timescale for Gaussian diagonal fluctuations (diagonal in the eigenbasis), but only fast (if any) fluctuations off-diagonal. This make the problem numerically tractable, however it may not always hold. Our aim here is to treat arbitrary timescales for off-diagonal (in eigenbasis) fluctuations by using the SLE and show significant differences from the more approximate treatments. This may lead to the development of better approximations, that retain the reduced bath treatment (necessary for large aggregates), but properly includes effects of bath timescale. For smaller aggregates the SLE can be directly used to calculate 2D spectra.

The SLE assumes the form:

dρdt=i[H+HSQ,ρ]+(dρdt)M+(dρdt)Q

The eigenvectors of the Fokker-Planck operator Eq. (14) are direct products of two single Fokker-Planck variables Q1, Q2

φαβ=φα1φβ2
(15)

The eigenvectors are given by

φαw=exp[(Qw/σw)2]2α2πα!σwα(Qwσw2)

where Hα are Hermite polynomials

α(x)=(1)αex2dαdxαex2.

αβ}, α = 0, 1, 2, …; β = 0, 1, 2, …; correspond to eigenvalues −αΛ1 − βΛ2, and form the suitable basis for a matrix representation of bath densities.36 The matrix elements of LQ are

[Q]αβ,γη=(αΛ1βΛ2)δα,γδβ,η
(16)

The Q variable is represented by the following matrix

[Q1]αβ,γη=(ασ12δα1,γ+σ12δα,γ1)δβη[Q2]αβ,γη=(βσ22δβ1,η+σ22δβ,η1)δαγ
(17)

If HSQ is a polynomial in Q1, Q2 it can be represented by the block-tridiagonal matrix, and the algorithm of Ref.10,20,36 may be employed to obtain the Green’s function by inverting the block-tridiagonal matrices (Eq. (5)). For instance, the linear coupling c1Q1+c2Q2+c12Q1Q2 becomes tridiagonal in blocks {ϕ00},{ϕ10, ϕ11, ϕ01}, {ϕ20, ϕ21, ϕ22, ϕ12, ϕ02}, …. However, the SLE approach is not limited to any particular form of the coupling; f can be arbitrary function.

The zero eigenvector ρQ;(eq) = ϕ00 (or in the matrix notation [ρeq]αβ = δα0δβ0) gives the equilibrium distribution. Summation over the bath ∫ dQ1dQ2 corresponds to zero left eigenvector in the eigenbasis {ϕα} thus left angle bracketIQ|αβ = δα0δβ0. The initial state is direct product of equilibrium densities of all the bath variables |ρ(eq)right angle bracketQ;(eq)right angle bracket and similarly the final summation.

The various Liouville space exciton manifold blocks of the Liouville superoperator LSQ = −i[HSQ, …] are

gg,ggSQ=0;eg,egSQ=(if1(Q1)0000if1(Q1)0000if2(Q2)0000if2(Q2))ee,eeSQ=(000000000000000000i(f1f2)00000000i(f1f2)00000000i(f1f2)00000000i(f1f2)000000000000000000)

Here f1 [equivalent] f1(Q1), f2 [equivalent] f2(Q2). The response can still be calculated using (Eq.(8)) but now in the higher space which combines Q with the u,d space.

Several methods have been employed to describe the signatures of spectral diffusion in the nonlinear response of aggregates. Standard simulations use the Redfield equations for the evolution of the exciton density matrix during t2, Gaussian peak of arbitrary relaxation timescale for peaks can be added.1 Static disorder in site frequencies can be easily incorporated. Slow fluctuations only account for peakshape, interlevel coupling fluctuation are fast, and dynamical effects on dipole moments are neglected.

The SLE framework is more general. Consider a symmetric dimer Ju = Jd, ε1 = ε2, µ1 = µ2 where only the symmetric eigenvector |E1right angle bracket = |1right angle bracket + |2right angle bracket is radiatively coupled to the ground state. For fast spectral diffusion the other peak |E2right angle bracket = |1right angle bracket−|2right angle bracket is dark. However when the slow coordinate nonequilibrium eigenstates |E2(Q1,Q2)right angle bracket couple to the dipole moment left angle bracketE2(Q1, Q2)|µ|gright angle bracket ≠ 0, the other peak may be observed (Fig 8, left panel). The SLE can consistently describe the fluctuation timescales where the additional diagonal peak vanishes (Fig 8, right panel). This is not properly described by simulations at the Redfield level, which do not account for the necessary correlations between Liouville space and bath coordinates during the transport. Additional static disorder can induce some oscillator strength to the dark state, but cannot describe its evolution.

Fig 8
The SA lineshape of a dimer with spectral diffusion for short Λt2 = 0 (left panel) and long Λt2 = 2 (right panel) times. Other parameters µ1 = µ2 = 1, Ju = Jd = J, f(Q) = Q , σ1 = σ2 = 0.3J, Λ1 = ...

At zero delay time the magnitude of the (−1,−1) peak is proportional to the probability [p with tilde] that bath coordinates shifts the site energies ε so that dipole moment left angle bracketE2(Q1,Q2)|µ|gright angle bracket is substantial. Such configurations are rather rare and far from equilibrium, so they disappear at relaxation timescale Λ−1. For t2 >> Λ−1 the (−1,−1) peak can only be induced by two independent hits of the ”[p with tilde]” area, so its magnitude is quadratically proportional to [p with tilde], in agreement with asymptotic formula Eq. (12).

These effects are illustrated in Fig 9, which depicts the magnitudes (normalized with respect to t2 = 0) of all peaks at Fig 8 as they vary with the delay time t2. We see significant decay of the (−1,−1) peak at the Λ−1 timescale as discussed. The (1,1) diagonal and the cross peak (1,−1) (the other cross-peaks (−1,1) has the same magnitude) change less significantly, because the cross peak requires instantaneous dipole moment for dark state, only in the t1 (peak (1,−1)) or t3 interval (peak (−1,1)). Its magnitude thus remains linear in p (for small p) for any delay times. In addition, we see oscillation of the peaks, caused by the coherent dynamics, as predicted by the Redfield level of theory34,37 and observed experimentally.38 Such oscillations require the interplay of both exciton states. The effect is thus most pronounced for (−1,−1) and cross peak, and not for the (1,1) peak.

Fig 9
The relative intensities of peaks (against t2 value) at I0(t2) = SA3, t2, ω1)/SA3, 0, ω1) lineshape of a dimer with spectral diffusion for varying delay times t2 at diagonal peak ω1 = ω3 = J (solid line), ...

7 Conclusions

It is straightforward to extend the present model to other types of excitonic systems as long as the fluctuations are Markovian. In fact, the code developed for the present calculations can treat any number of ground, excited and doubly excited states, and multistate or Gaussian -Markovian fluctuations (coordinate linearly or quadratically coupled to Q coordinate). Since the computational time grows exponentially with the number of stochastic coordinates, applications to extended excitonic systems will require the identification of a few relevant coordinates.

We showed that the SLE approach is essential tool for describing intermediate timescale of fluctations and may thus be used to test various approximate methods39,40 developed in order to avoid the numerically expensive explicit representation of stochastic coordinate.

The SLE assume that the stochastic process is independent on the state of the system. This is justified at high temperatures. Finite-temperature effects such as the Stokes-shift or nonuniform exciton distributions are missed. The Redfield equations in contrast, contain finite temperature effects. The high temperature ћΔε < kT limit usually holds for infrared,19,21 but not in optical domain at room temperatures. Finite temperature corrections for linear coupling include system dependent bath dynamics.16,4144 Provided the temperature is not too low, the equations of motion have the same complexity as the SLE and may be readily implemented, if necessary. This level was widely used to describe electron transfer processes (adiabatic and nonadiabatic). At very low temperatures, when the quantization of bath degrees of freedom is necessary, some additional oscillator modes with Matsubara frequencies has to be added increasing computation costs.16

Acknowledgements

The support of the GAČR (Grant No. 202/07/P245), the Ministry of Education, Youth and Sports of the Czech Republic (MSM 0021620835) (F.Š.), NSF (CHE-0745892) and NIH (GM59230) (S.M.) is gratefully acknowledged.

Appendix B

A The four-state-jump (FSJ) model

We consider a single two level chromophore described by the Hamiltonian

H=ε(t)BB
(18)

Where εs assumes four possible values s = 1, 2, 3, 4. The model has a four-peak linear spectrum indistinguishable from the TSJ dimer (Fig 2). This model has 12 rates ksr for the transitions between r to s states. The time evolution is described by the master equation

dρsdt=rskrsρs+rsksrρr
(19)

Eq. (19) further describes the evolution in the gg and ee manifolds because no coherence is connected with the ground (excited) state evolution.

The SLE for the coherence evolution is

dρegsdt=iεjρegsrskrsρegs+rsksrρegr

The rates were chosen to yield the same linear spectrum as the TSJ dimer in the slow limit

I(ω)=µ2υ=14Γυρggυ(eq)[Γυ2+(ω1ευ)2]
(20)

with Γv = ∑rv krv.

Consider 2D lineshape in the slow jumps limit krs/(εr − εs) << 1. At short times krst2 << 1 lineshapes show only diagonal peaks at fundamental frequencies. Cross-peaks are absent (Fig 3 bottom left panel)

SA(ω3,t2=0,ω1)=µ4υ=14Γυ2ρggυ(eq)[Γυ2+(ω1ευ)2][Γυ2+(ω3ευ)2]

Since there is no coherence during t2, the SA lineshape will only change at k−1 timescale, where the cross peaks appear. Asymptotically the lineshape approaches the factorized form

ћSA(ω3,t2,ω1)=4I(ω1)I(ω3).
(21)

This limit holds for all one exciton stochastic lineshapes with finite memory; see30 for some exceptions when the memory persists for arbitrary timescales. A similar relation can be generalized for the first manifold contribution of any excitonic system, that allows the excitonic transfer toward the thermodynamic high temperature uniform excitonic distribution

ћSA(ω3,t2,ω1)=2n+1nI(ω1)I(ω3).
(22)

where n is the number of excitons.31

References

1. Abramavicius D, Mukamel S. Chem. Rev. 2004;104:2073. [PubMed]
2. Woutersen S, Hamm P. J. Phys.: Condens. Matter. 2002;14:R1035.
3. Meinhold L, Smith JC, Kitao A, Zewail AH. PNAS. 2007;104:17261. [PubMed]
4. Nilsson L, Halle B. PNAS. 2005;102:13867. [PubMed]
5. Lim M, Hamm P, Hochstrasser RM. PNAS. 1998;95:15315. [PubMed]
6. Ganim Z, Tokmakoff A. Biophys J. 2006;91:2636. [PubMed]
7. Cowan ML, Bruner BD, Huse N, Dwyer JR, Chugh B, Nibbering ETJ, Elsaesser T, Miller RJD. Nature. 2005;434:199. [PubMed]
8. Brixner T, Stenger J, Vaswani HM, Cho M, Blankenship RE, Fleming GR. Nature. 2005;434:625. [PubMed]
9. Jansen TLC, Knoester J. J. Chem. Phys. 2007;127:234502. [PubMed]
10. Jansen TLC, Zhuang W, Mukamel S. J. Chem. Phys. 2004;121:10577. [PubMed]
11. Abramavicius D, Zhuang W, Mukamel S. J. Phys. Chem. B. 2004;108:18034.
12. Redfield AG. Adv. Magn. Reson. 1965;1:1.
13. May V, Kuhn O. Charge and Energy Transfer Dynamics in Molecular Systems. Weinheim: Wiley-VCH; 1999.
14. Szöcz V, Palszegi T, Lukeš V, Sperling J, Milota F, Jakubetz W, Kauffmann HF. J. Chem. Phys. 2006;124:124511. [PubMed]
15. Kubo R. J. Math. Phys. 1963;4:174.
16. Tanimura Y. J. Phys. Soc. Jpn. 2006;75 082001.
17. Gamliel D, Levanon H. Stochastic processes in Magnetic Resonance. World Scientific; 1995.
18. Jansen TLC, Hayashi T, Zhuang W, Mukamel S. J. Chem. Phys. 2005;123:114504. [PubMed]
19. Kim YS, Hochstrasser RM. PNAS. 2005;102:1185. 2005.
20. Šanda F, Mukamel S. J. Chem. Phys. 2006;125 014507. [PubMed]
21. Zheng J, Kwak K, Asbury J, Chen X, Piletic IR, Fayer MD. Science. 2005;309:1338. [PubMed]
22. Davydov AS. Theory of Molecular Excitons. New York: McGraw-Hill; 1962.
23. Silinsh EA, Čápek V. Organic molecular crystals:Interaction, Localization and Transport phenomena. New York: AIP Press; 1994.
24. van Amerogen H, Valkunas L, van Grondelle R. Photosyntetic excitons. Singapore: World Scientific; 2000.
25. van Kampen NG. Stochastic processes in Physics and Chemistry. North Holland: Amsterdam; 1992.
26. Mukamel S. Principles of Nonlinear Optical Spectroscopy. New York: Oxford Univ. Press; 1995.
27. Jonas DM. Annu. Phys. Rev. Chem. 2003;54:425. [PubMed]
28. Khalil M, Demirdöven N, Tokmakoff A. Phys. Rev. Lett. 2003;90 047401. [PubMed]
29. Šanda F, Mukamel S. Phys. Rev. Lett. 2007;98 080603. [PubMed]
30. Šanda F, Mukamel S. J. Chem. Phys. 2007;127:154107. [PubMed]
31. As noted in30 for single chromophore, Eq. (22) may be generalized beyond the stochastic model. For multiexcitonic, but exciton conserving Hamiltonian, the four diagrams of Fig 1 contributes ћSA3, t2 = ∞, ω1) = I1)[I3)+IE3)] where IE is emission line-shape obtained from excited state equilibrium ρν, assuming that the asymptotic Green’s function (of gg and ee manifold) approach [G]νν′(t → ∞) = ρνTrν′ in some suitable representation, where ν is multiindex (quantum or classical) completely describing the excited (ground) state of the system. Additional contribution to SA could, however, come from the excited state absorption (involving the second exciton manifold).
32. Scheurer Ch, Mukamel S. J. Chem. Phys. 2001;115:4989.
33. Ernst RR, Bodenhausen G, Wokaun A. Principles of Nuclear Magnetic Resonance in One nad Two Dimensions. New York: Oxford University Press; 1987.
34. Pisliakov AV, Mančal T, Fleming GR. J. Chem. Phys. 2006;124:234505. [PubMed]
35. Alicki R, Lendi K. Quantum Dynamical Semigroups and Applications. Berlin: Springer; 1987.
36. Risken H. The Fokker-Plank equation. Berlin: Springer; 1989.
37. Kjellberg P, Brüggemann B, Pullerits T. Phys. Rev. B. 2006;74 024303.
38. Engel GS, Calhoun TR, Read EL, Ahn TK, Mančal T, Cheng Y-C, Blankenship RE, Fleming GR. Nature. 2007;446:782. [PubMed]
39. Kim HD, Tanimura Y, Cho M. J. Chem. Phys. 2007;127 075101.
40. Abramavicius D, Valkunas L, Mukamel S. Europhys. Lett. 2007;80:17005.
41. Garg A, Onuchic JN, Ambegaokar V. J. Chem. Phys. 1985;83:4491.
42. Zusman D. Chem Phys. 1980;49:295.
43. Chernyak V, Mukamel SJ. Chem. Phys. 1996;105:4565.
44. Ishizaki A, Tanimura YJ. Chem. Phys. 2006;125 084501.