Stochastic Liouville Equations for Coherent Multidimensional

Oct 16, 2008 - ... Diversity in Inorganic/Organic Chemistry. ACS Publications will host a free two-day forum featuring world-class researchers in the ...
0 downloads 0 Views 2MB Size
Download PDF
14212

J. Phys. Chem. B 2008, 112, 14212–14220

Stochastic Liouville Equations for Coherent Multidimensional Spectroscopy of Excitons Frantisˇek Sˇanda*,† and Shaul Mukamel‡ Faculty of Mathematics and Physics, Institute of Physics, Charles UniVersity, Ke KarloVu 5, Prague, 121 16 Czech Republic, and Department of Chemistry, UniVersity of California, IrVine, California 92697-2025 ReceiVed: February 18, 2008; ReVised Manuscript ReceiVed: August 13, 2008

Signatures of chemical exchange and spectral diffusion in 2D photon-echo line shapes 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 have different dynamics 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. (Un)jj′ ) exp[-iεj(tn)∆t]δjj′

1. Introduction Exciton models are widely used in the description of the coherent nonlinear optical response of complex systems.1 Vibrational spectra of proteins2-6 and liquid water7 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 toward 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 |φj(t)〉 and energies εj(t) (we set p ) 1)

H(t)|φj(t)〉 ) εj(t)|φj(t)〉 Direct evaluation of the time evolution operator, when the time t is divided into N f ∞ short segments ∆t ≡ t/N gives N

exp[-iH(tn)∆t] [ ∫0t dτH(τ)] ) T∏ n)1

U ≡ T exp -i

is the adiabatic propagator, and the nonadiabatic coupling matrix

(Sn)jj′ ) 〈φj(tn)|φj′(tn-1) 〉 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)15-17 allow the modeling 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

tn ) n∆t (1) where T is time ordering operator. By adopting the matrix ˜ jk ≡ 〈φj(t)|U|φk(0)〉, we get representation U

˜ )T U

N

UnSn ∏ n)1

(2)

Here * To whom correspondence should be addressed. E-mail: sanda@ karlov.mff.cuni.cz (F.Sˇ.); [email protected] (S.M.). † Charles University. ‡ University of California.

Figure 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 ) k1 - k2 + k3 (eq 7).

10.1021/jp801457c CCC: $40.75  2008 American Chemical Society Published on Web 10/17/2008

Coherent Multidimensional Spectroscopy of Excitons

J. Phys. Chem. B, Vol. 112, No. 45, 2008 14213 A discrete fluctuation model is introduced by an additional bath degree of freedom with two states u (up) and d (down), representing, for example, 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.

( ) ( )

dFu ) -kdFu + kuFd dt M dFd ) kdFu - kuFd dt M

Figure 2. Absorption line shapes (eq 6) of a dimer with two-state bath fluctuations (eqs 3 and 4) (solid line) and a monomer with fourstate bath fluctuations (eqs 18 and 19) (dotted line) in the static limit. The dimer peaks correspond to resonance with (from left to right) ε2u, ε2d, ε1d, ε1u levels. Parameters: Dimer: ε1u ) 30Γ, ε1d ) 10Γ, ε2u ) -30Γ, ε2d ) -10Γ, Ju ) Jd ) 0, ku ) kd ) k ) Γ. Monomer: ε1 ) 30Γ, ε2 ) 10Γ, ε3 ) -30Γ, ε4 ) -10Γ, k12k21 ) k13 ) k31 ) k14 ) k41 ) k23 ) k32 ) k24 ) k42 ) k34 ) k43 ) Γ/3, µ4sj ) √2 µdim. All the lineshapes shown at Figures 28 are calculated by solving the full SLE model, and not by using the various limiting approximations discussed in the main text.

calculating the density distributions in the joint system + bath space. The SLE generalize the Redfield equations, which are limited to fast fluctuations where the bath degrees of freedom can be eliminated. We have recently applied the SLE to calculate 2D line shapes 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 fluctuation 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 will account 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 nonexcitonic 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 the Molecular Dimer We consider a dimer with two one-exciton states, |e1〉 ) B†1|g〉 and |e2〉 ) B†2|g〉, and one doubly excited state, |f1〉 ≡ B†1B†2|g〉, described by the Frenkel-exciton Hamiltonian1,22-24

H ) ε1(t)B†1B1 + ε2(t)B†2B2 + J(t)(B†2B1 + B†1B2) + ∆B†1B†2B1B2 (3) Both the coupling J and energies j undergo stochastic fluctuations. We will assume that the doubly excited state is strongly shifted, ∆ f ∞, so that it does not contribute to the optical signals near the single-exciton frequency ω ≈ ε.

(4)

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

dF dF ) -i[H, F] + dt dt

( )

M

≡ LF

Here, F carries three indices representing bra and ket variables of the system Liouville space and one index denoting the state of the 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; u〉〉, |gg; d〉〉), with the block

Lgg,gg )

(

-kd ku kd -ku

)

The coherence (eg) evolution is described in the following basis set |e1g; u〉〉, |e1g; d〉〉, |e2g; u〉〉, and |e2g; d〉〉 by the Liouvillean

(

-kd - iε1u ku -iJu 0 kd -ku - iε1d -iJd 0 Leg,eg ) -iJu ku -kd - iε2u 0 k -k -iJ 0 d d u - iε2d

)

Finally, the excited-state block |e1e1; u〉〉, |e1e1; d〉〉, |e1e2; u〉〉, |e1e2; d〉〉, |e2e1; u〉〉, |e2e1; d〉〉, |e2e2; u〉〉, |e2e2; d〉〉 is represented by

(

Lee′,ee′ ) -kd ku iJu -iJu 0 0 0 0 kd -ku iJd -iJd 0 0 0 0 iJu 0 -kd - i∆εu ku -iJu 0 0 0 kd -ku - i∆εd 0 iJd 0 0 0 -iJd -iJu 0 ku iJu 0 -kd + i∆εu 0 0 -ku + i∆εd 0 iJd kd 0 -iJd 0 0 -iJu iJu -kd ku 0 0 0 0 kd -ku -iJd iJd 0 0 0 0

)

where we have denoted ∆εd ≡ ε1d - ε2d, ∆εu ≡ ε1u - ε2u. With Lge,ge ) L /eg,eg, the description is complete. The Green’s function solution of the Liouville equation is

G(t) ≡ θ(t)exp L t In the frequency domain, G(ω) ) ∫∞0 eiωtG(t)dt, we have

G(ω) ) -[iω + L ]-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

Sˇanda and Mukamel

14214 J. Phys. Chem. B, Vol. 112, No. 45, 2008

Ggg,gg(t) )

( )

(

ku ku 1 e-(kd+ku)t kd -ku + kd + ku kd kd kd + ku -kd ku

)

) S/(-ω)). We focus on the ω1 ≈ ε peaks. The absorption line shape is given by

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 of the block matrix inversions (eq 5) can be performed analytically in closed, though somewhat lengthy, form. Averaging over the bath variable is made with respect to the equilibrium bath density

|F(exp)〉 )

()

ku 1 k kd + ku d

and the summation over final states is represented by the scalar product with the vector 〈I| ) (1, 1). 3. The Two-Dimensional Signals

(eq) I(ω1) ) Im S(ω1) ) Re〈I|µgg,egGeg,eg(ω1)µ(-) 〉 (6) eg,gg|F

We next turn to the third-order resonant signals generated in the kI ) -k1 + k2 + k3 and kII ) k1 - k2 + k3 phase-matching directions. These will be displayed as ω1, ω3 frequency-frequency plots S(ω3,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 giving peaks at opposite frequencies S(ω3,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

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

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 four Feynman diagrams shown in Figure 1

Dint ) [µ1(B†1 + B1) + µ2(B†2 + B2)] is the dipole operator. We consider impulsive experiments performed with short laser pulses. The linear response function is

R1(ω3, t2, ω1) )

S(t1) ) i〈〈µ|G(t1)µ(-)|F(0) 〉〉

R2(ω3, t2, ω1) ) (eq) (-) (-) 〈I|µgg,egGeg,eg(ω3)µ(-) 〉 eg,eeGee,ee(t2)µee,geGge,ge(ω1)µge,gg|F

and the third-order response function26

S(t3, t2, t1) ) i3〈〈µ|G(t3)µ(-)G(t2)µ(-)G(t1)µ(-)|F(0)〉〉

R3(ω3, t2, ω1) )

where µ(-) ≡ [Dint,...] is the dipole superoperator and 〈〈µ|... ) 〈I|TrDint... provides the mean value of dipole moment. Initially, the system is in the ground state, and the bath is at equilibrium

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

(

µ1 0 0 µ1

)

µ ˆ2 )

(

µ2 0 0 µ2

(-) (-) (eq) 〈I|µgg,egGeg,eg(ω3)µ(-) 〉 eg,ggGgg,gg(t2)µgg,egGeg,eg(ω1)µeg,gg|F

R4(ω3, t2, ω1) ) (-) (-) (eq) 〈I|µgg,egGeg,eg(ω3)µ(-) 〉 (8) eg,eeGee,ee(t2)µee,egGeg,eg(ω1)µeg,gg|F

|F(0) 〉〉 ) |gg 〉〉 |F(eq) 〉

µ ˆ1 )

(-) (-) (eq) 〈I|µgg,egGeg,eg(ω3)µ(-) 〉 eg,ggGgg,gg(t2)µgg,geGge,ge(ω1)µge,gg|F

)

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

SA(ω3, t2, ω1) ≡ -Im SI(ω3, t2, -ω1) - Im SII(ω3, t2, ω1)

(9)

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

µ(-) eg,gg )

()

4. Slow Fluctuations

µ ˆ1 µ ˆ2

( ) ( )

The absorption spectrum (Figure 2) has four peaks, centered at two E1(2),u ) (ε1u + ε2u)/2 ( (∆ε2u + J2u)1/2 eigenvalues of Hu and two E1(2),d ) (ε1d + ε2d)/2 ( (∆ε2d + J2d)1/2 eigenvalues of Hd, in the static limiting case ku, kd , |E1(2),u - E1(2),d|

-µ ˆ1 0 -µ ˆ2 0 µ(-) ee,eg ) ˆ1 0 -µ ˆ2 0 -µ µ ˆ1 0 µ(-) ee,ge ) µ ˆ2 0

2

I(ω) )

0 ˆ1 µ 0 ˆ2 µ

(-) (-) (-) T (-) We further have µ(-) ge,gg ) -µeg,gg, µgg,eg ) [µeg,gg] , µgg,ge (-) T (-) (-) T (-) (-) T [µge,gg] , µeg,ee ) [µee,eg] , µee,ge ) [µge,ee] . We note that µR,β (-) (-) [µβ,R]T and µij,kl ) -µji,lk for ij(kl) * ee.

∑∑

j)1 v)u;d

2 Mj,v ΓvF(eq) ggv

Γ2v + (ω - Ej,v)2

(10)

where

Mj,v ) 〈Ej,v|µ|g 〉 ) )

The linear response is given by (eq) S(t1) ) i〈I|µgg,egGeg,eg(t1)µ(-) 〉 + c.c. eg,gg|F

In frequency space, S(ω1) ≡ ∫∞0 eiω1t1S(t1). These represent a symmetric peak structure around ε and -ε, respectively, (S(ω)

is the transition dipole and Γd ) ku and Γu ) kv are the line broadening parameters. The 2D SA line shape at short t2, displayed in Figure 3, top left panel, shows primarily four diagonal peaks and four selected cross-peaks. Three consecutive interactions with laser pulses -1 (applications of µ(-)) for times t2 , kd,u induce cross-peaks between peaks belonging to the same bath state but different excitons.

Coherent Multidimensional Spectroscopy of Excitons

J. Phys. Chem. B, Vol. 112, No. 45, 2008 14215

SA(ω3, t2 , k-1 d,u , ω1) ) 2

Re

[

2

1 ∑ ∑ ∑ M2kvMjv2F(eq) ggv Γv - i(ω3 - Ej,v)

k)1 j)1 v)u;d

2Γv 2 Γv + (ω1 - Ek;v)2

+ e-Γv,kjt2-i(Ej,v-Ek;v)t2 e

(

×

1 + Γv - i(ω1 - Ej,v)

1 Γv + i(ω1 - Ek,v)

)]

(11)

where exciton dephasing Γev,kj is negligible in the slow limit (compare fast fluctuation limit, section 5). Additional cross-peaks induced by jumps among the u and d peaks (top left panel) appear for t2 . k-1 d,u. The dynamics is quite rich on these time scales since the phase relation between u and d eigenstates 〈Ej,u|Ei,d〉 is important for ee evolution and may also induce fluctuations of the peaks magnitudes at the 2π(Ej,v - Ei;v)-1 period. We distinguish between two cases. If the u and d eigenstates are degenerate |Ej,u〉 ) |Ej,d〉, no exciton transport is allowed. This is shown in the top panels of Figure 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,u〉 * |Ej,d〉, exciton transfer is allowed, and asymptotically, we have

kv 1 GEj jEj kv,Ej lEj mw(t2 f ∞) ) δjkδlm 2 ku + kd The absorptive line shape can then be factorized as follows

SA(ω3, t2 f ∞, ω1) ) 3I(ω1)I(ω3)

(12)

This is similar to the asymptotic line shape of the monomer (eq 21).29,30 Equation 12 holds for an arbitrary time scale and for any type of stochastic modulation (once transport is allowed). Equation 12 can be generalized beyond the stochastic model31 to include finite temperature effects, where the emission line shape is red shifted with respect to the absorption (the Stokes shift). Line shapes of these 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 time scales (with differences in models and approximations used) have no phase matching but have 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 those used here, particularly adequate for interpreting 2D optical and infrared experiments. For a more detailed analysis of these line shapes, we display the four contributions Ri to SA in Figure 4. At t2 ) 0, R1 ) R2 (Figure 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. Figure 5 repeats these calculations in the long delay time 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 line shape may be recast as the product of the two line shapes, ∼2I(ω1)I(ω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 have 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 line shape. The features shown in Figure 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,v〉 eigenstates or when some levels are accidently degenerate. The 2D line shapes 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 fourstate jump bath model (FSJ) modulating a single chromophore,

Figure 3. The SA signal (eq 9) for a dimer with a two-state bath (top panels) and a monomer with a four-state bath (bottom panels) corresponding to Figure 2. For delay times Γt2 ) 0, 0.1, 0.2, 1.0, 100 increasing from left to right. Other parameters are the same as those in Figure 2.

Sˇanda and Mukamel

14216 J. Phys. Chem. B, Vol. 112, No. 45, 2008

Figure 4. The R1(ω3,t2,-ω1), R2(ω3,t2,-ω1), R3(ω3,t2,ω1), and R4(ω3,t2,ω1) (eq 8, real parts) contributions to the SA signal (from left to right, then top to bottom) for the dimer with a two-state bath for short delay times, t2 ) 0, corresponding to the most left top panel of Figure 2. Other parameters are the same as those in Figure 2.

Figure 6. (A) Top: The SA (eq 9) signal for the dimer with a twostate 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 splitting E∆ ≡ [(ε1u + ε1d)/ 4] - [(ε2u + ε2d)/4]. Parameters: Top: ε1u ) 1.5E∆, ε1d ) 0.5E∆, ε2u ) -1.5E∆, ε2d ) 0.5E∆, Ju ) -Jd ) 0.68E∆, ku ) kd ) k ) 0.5E∆, µ1 ) µ2 ) 1.189 µ. Bottom: ε1u ) 6.5E∆, ε1d ) -0.45E∆, ε2u ) -6.5E∆, ε2d ) 4.5E∆, Ju ) -Jd ) 0.75E∆, ku ) kd ) k ) 100E∆, µ1 ) µ2 ) µ. (B) Absorption line shapes of the model of (A). The dashed line corresponds to fast limit (bottom of A), and the solid line corresponds to the intermediate limit (top of A). Figure 5. The same as in Figure 4 but for long delay times Γt2 ) 100, corresponding to the most right top panel of Figure 3.

5. Fast and Intermediate Fluctuations Time Scale given in Appendix A. Comparison of eq 10 with eq 20 or the line shapes displayed in Figure 2 shows that the two models may have identical absorption line shapes in the static limit with the proper choice of the resonance frequencies and line width (as defined by eqs 10 and 20). 2D line shapes of the two models are markedly different at short times, as shown in Figure 3. The slow fluctuation FSJ limit (bottom left) SA line shape at short t2 jumps only reveals the diagonal peaks because no change of bath state can occur, in contrast to the 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 cannot be observed. Asymptotically, these two line shape agree (compare eqs 12, 21, and 22), except when exciton transfer is not allowed.

In the motional narrowing ku, kd . |Eju - Ekd| limit, fluctuations are averaged on the time scale of optical experiments. The averaged Hamiltonian

j ) ε1B†1B1 + ε2B†2B2 + jJ(B†2B1 + B†1B2) H

(13)

where

jJ )

kuJu + kdJd ku + kd

εj )

kuεju + kdεjd ku + 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) j )〉|E j 2〉|2〉/k and dephasing rates Γe ∼ ΓT + j 1|(H - H ΓT ∼ 〈|〈E 2 2 j j j ∑i)1 〈Ei|〈(H - H) 〉|Ei〉/k and dephasing between ground and j i|〈(H - H j )2〉|E j i〉/k + 〈g|〈(H excited states Γi ∼ ΓT/2 + 〈E

Coherent Multidimensional Spectroscopy of Excitons

J. Phys. Chem. B, Vol. 112, No. 45, 2008 14217

j )2〉|g〉/k. The averaged density matrix is approximately deH scribed by the Redfield master equations in the secular approximation in the high-temperature limit

dFEj 1Ej 1 dt dFEj 2Ej 2 dt dFEj 1Ej 2 dt

) -ΓTFEj 1Ej 1 + ΓTFEj 2Ej 2 ) -ΓTFEj 2Ej 2 + ΓTFEj 1Ej 1

j 2) - Γe]FEj Ej j1 - E ) [-i(E 1 2

dFEj jg dt

j j - Γj]FEj g ) [-iE j

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

GEj 1Ej 1,Ej 1Ej 1(t) ) GEj 2Ej 2,Ej 2Ej 2(t) )

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

-2ΓTt

1+e 2

1 - e-2Γ GEj 1Ej 1,Ej 2Ej 2(t) ) GEj 2Ej 2,Ej 1Ej 1(t) ) 2

Tt

j

j

GEj 1Ej 2,Ej 1Ej 2(t) ) GEj 2Ej 1,Ej 2Ej 1*(t) ) e[-i(E1-E2)-Γ ]t e

j

GEj jg,Ej jg(t) ) GgEj j,gEj j*(t) ) e[-iEj-Γj]t Note that the SLE provides an exact representation of stochastic Markovian fluctuations of Hamiltonian dynamics. Both the Liouville equation and the averaging over bath paths preserve the positive semidefiniteness 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 j 1 and E j 2 of the averaged dimer Hamiltonian eigenvalues E

Mi2Γi

2

I(ω) )

∑ (Γ )2 + (ω - Ej )2 i)1

i

i

j1 - E j 2, these can eventually overlap and merge If Γ1 + Γ2 > E into a single peak. The 2D line shape (Figure 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) ) 2

2

Γ

∑ ∑ 2Mi2Mj2 (Γ )2 + (ωi i)1 j)1

2

∑

i,j,k,l)1

j

1 - Ei)

i

Re

Γj 2

j j)2 (Γj) + (ω3 - E 2

MiMjMkMlGEj lEj k,Ej jEj i(t2) j l) Γl - i(ω3 - E

[

+

×

1 1 + j ji Γj - i(ω1 - Ej) Γi + i(ω1 - E

]

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

Figure 8. The SA line shape 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 ) Λ2 ) Λ ) 0.1J, ε1d ) ε1u ) ε1d ) ε1u ) 0.0.

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. Figure 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 the 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 line shapes. At longer delay times (right panels), the decoherence and exciton transport erase all memory, and line shapes 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 time scale shows additional delicate effects caused by altering the eigenvectors (rather then merely eigenvalues) of the Hamiltonian. These will be discussed in the next section using a Gaussian model for spectral diffusion. We next examine more closely the R3 and R4 contributions (Figure 7). At short times GEj lEj k,Ej jEj i ) δljδik, and R3 follows the four peak structure identical to R1, while R4 is diagonal, for the

Sˇanda and Mukamel

14218 J. Phys. Chem. B, Vol. 112, No. 45, 2008 same reason as that mentioned in Figure 4. At long times, t2 coherences are dumped. Asymptotically, GEj lEj k,Ej jEj i ) 1/2δlkδij, and the line shape (eq 12) is restored. This argument holds as long as exciton transfer is allowed, δJ * 0; otherwise, the picture of Figure 5 is repeated. and the cross-peaks never show up in R4. 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))B†1B1 + f2(Q2(t))B†2B2 to the Hamiltonian, together with a Fokker-Planck master equation for the density evolution of the two collective coordinates Q1 and Q2

(

∂F(Q1, Q2) ∂t

)

(

2

) Q

)

∑ Λw ∂Q∂ w Qw + σw2 ∂Q∂ w F(Q1, Q2) ) L QF

w)1

( (

[Q1]Rβ,γη ) Rσ1√2δR-1,γ +

dF dF ) -i[H + HSQ, F] + dt dt

( ) + ( dFdt ) M

Q

The eigenvectors of the Fokker-Planck operator eq 14 are direct products of two single Fokker-Planck variables Q1, Q2

φRβ ) φ1Rφβ2 The eigenvectors are given by

φwR )

exp[-(Qw/σw)2] R

2 √2πR ! σw

(15)

( )

HR

Qw

σw√2

where HR are Hermite polynomials

HR(x) ) (-1)Rex

2

dR -x2 e dxR

where {φRβ}, R ) 0, 1, 2,... and β ) 0, 1, 2,..., correspond to eigenvalues -RΛ1 - βΛ2 and form the suitable basis for a matrix representation of bath densities.36 The matrix elements of L Q are

[L Q]Rβ,γη ) (-RΛ1 - βΛ2)δR,γδβ,η The Q variable is represented by the following matrix

(16)

) )

δR,γ-1 δβη √2 σ2 [Q2]Rβ,γη ) βσ2√2δβ-1,η + δβ,η-1 δRγ √2

(17)

If HSQ is a polynomial in Q1 and Q2, it can be represented by the block-tridiagonal matrix, and the algorithm of refs 10, 20, and 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 an arbitrary function. The zero eigenvector FQ;(eq) ) φ00 (or in the matrix notation eq [F ]Rβ ) δR0δβ0) gives the equilibrium distribution. Summation over the bath ∫dQ1dQ2 corresponds to zero left eigenvector in the eigenbasis {φR}; thus, 〈IQ|Rβ ) δR0δβ0. The initial state is a direct product of equilibrium densities of all of the bath variables |F(eq)〉|FQ;(eq)〉 and, similarly, the final summation. The various Liouville space exciton manifold blocks of the Liouville superoperator L SQ ) -i[HSQ,...] are

(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 the averaged Hamiltonian, allowing an arbitrary time scale for Gaussian diagonal fluctuations (diagonal in the eigenbasis) but only fast (if any) fluctuations for the off-diagonal. This make the problem numerically tractable, however, it may not always hold. Our aim here is to treat arbitrary time scales for off-diagonal (in eigenbasis) fluctuations by using the SLE and to 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 the bath time scale. For smaller aggregates, the SLE can be directly used to calculate 2D spectra. The SLE assumes the form

σ1

(

-if1(Q1) 0 0 0 (Q ) -if 0 0 0 1 1 SQ Leg,eg ) (Q ) -if 0 0 0 2 2 -if2(Q2) 0 0 0

SQ Lee′,ee′ )

0 0 0 0 0 0 0 0

(

SQ Lgg,gg )0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 -i(f1 - f2) -i(f1 - f2) 0 0 0 0 i(f1 - f2) 0 0 0 0 i(f1 - f2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0

) 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

)

Here f1 ≡ f1(Q1) and f2 ≡ 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, and a Gaussian peak of an arbitrary relaxation time scale for peaks can be added.1 Static disorder in site frequencies can be easily incorporated. Slow fluctuations only account for peak shape, 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 |E1〉 ) |1〉 + |2〉 is radiatively coupled to the ground state. For fast spectral diffusion, the other peak |E2〉 ) |1〉 - |2〉 is dark. However, when the slow coordinate nonequilibrium eigenstates |E2(Q1, Q2)〉 couple to the dipole moment 〈E2(Q1, Q2)|µ|g〉 * 0, the other peak may be observed (Figure 8, left panel). The SLE can consistently describe the fluctuation time scales where the additional diagonal peak vanishes (Figure 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

Coherent Multidimensional Spectroscopy of Excitons

Figure 9. The relative intensities of peaks (against t2 value) at the I0(t2) ) SA(ω3, t2, ω1)/SA(ω3, 0, ω1) line shape of a dimer with spectral diffusion for varying delay times t2 at the diagonal peak ω1 ) ω3 ) J (solid line), the cross-peak at ω1 ) -ω3 ) J (dashed line), and the transitory diagonal ω1 ) ω3 ) -J (dotted line). Other parameters are as those in Figure 8.

the transport. Additional static disorder can induce some oscillator strength to the dark state but cannot describe its evolution. At zero delay time, the magnitude of the (-1,-1) peak is proportional to the probability p˜ that bath coordinates shift the site energies ε so that the dipole moment 〈E2(Q1, Q2)|µ|g〉 is substantial. Such configurations are rather rare and far from equilibrium; therefore, they disappear at relaxation time scale Λ-1. For t2 . Λ-1, the (-1,-1) peak can only be induced by two independent hits of the “p˜” area; therefore, its magnitude is quadratically proportional to p˜, in agreement with asymptotic formula eq 12. These effects are illustrated in Figure 9, which depicts the magnitudes (normalized with respect to t2 ) 0) of all peaks in Figure 8 as they vary with the delay time t2. We see significant decay of the (-1,-1) peak at the Λ-1 time scale, as discussed. The (1,1) diagonal and the cross-peak (1,-1) (the other crosspeak (-1,1) has the same magnitude) change less significantly because the cross-peak requires an instantaneous dipole moment for the dark state only in the t1 (peak (1,-1)) or t3 intervals (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 the cross-peak and not for the (1,1) peak. 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 the 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 an essential tool for describing the intermediate time scale of fluctations and may

J. Phys. Chem. B, Vol. 112, No. 45, 2008 14219 thus be used to test various approximate methods39,40 developed in order to avoid the numerically expensive explicit representation of the stochastic coordinate. The SLE assume that the stochastic process is independent of 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 p∆ε < kT limit usually holds for the infrared19,21 but not in the optical domain at room temperature. Finite-temperature corrections for linear coupling include system-dependent bath dynamics.16,41-44 Provided that 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 have to be added, increasing computation costs.16 ˇ R (Grant No. 202/ Acknowledgment. The support of the GAC 07/P245), the Ministry of Education, Youth and Sports of the Czech Republic (MSM 0021620835) (F.Sˇ.), the NSF (CHE0745892), and the NIH (GM59230) (S.M.) is gratefully acknowledged. Appendix A. The Four-State Jump (FSJ) Model We consider a single two-level chromophore described by the Hamiltonian

H ) ε(t)B†B

(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 (Figure 2). This model has 12 rates ksr for the transitions between r and s states. The time evolution is described by the master equation

dFs ksrFr ) - krsFs + dt r*s r*s

∑

∑

(19)

Equation 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

dFegs krsFegs + ksrFegr ) -iεsFegs dt r*s r*s

∑

∑

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

ΓVF(eq) ggV

4

I(ω) ) µ2

∑ [Γ2 + (ω

V)1

2 1 - εV) ]

V

(20)

with ΓV ) ∑r*VkrV. Consider the 2D line shape in the slow jumps limit krs/(εr εs) , 1. At short times, krst2 , 1 line shapes show only diagonal peaks at fundamental frequencies. Cross-peaks are absent (Figure 3 bottom left panel) 4

SA(ω3, t2 ) 0, ω1) ) µ

4

∑ [Γ2 + (ω

V)1

V

Γ2VF(eq) ggV

2 2 2 1 - εV) ][ΓV + (ω3 - εV) ]

Since there is no coherence during t2, the SA line shape will only change at the k-1 time scale, where the cross-peaks appear. Asymptotically, the line shape approaches the factorized form

Sˇanda and Mukamel

14220 J. Phys. Chem. B, Vol. 112, No. 45, 2008

SA(ω3, t2 f ∞, ω1) ) 4I(ω1)I(ω3)

(21)

This limit holds for all one-exciton stochastic line shapes with finite memory; see ref 30 for some exceptions when the memory persists for arbitrary time scales. A similar relation can be generalized for the first manifold contribution of any excitonic system, which allows the excitonic transfer toward the thermodynamic high-temperature uniform excitonic distribution

SA(ω3, t2 f ∞, ω1) ) 2

n+1 I(ω1)I(ω3) n

(22)

where n is the number of excitons.31 References and Notes (1) Abramavicius, D.; Mukamel, S. Chem. ReV. 2004, 104, 2073. (2) Woutersen, S.; Hamm, P. J. Phys.: Condens. Matter 2002, 14, R1035. (3) Meinhold, L.; Smith, J. C.; Kitao, A.; Zewail, A. H. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 17261. (4) Nilsson, L.; Halle, B. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13867. (5) Lim, M.; Hamm, P.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 15315. (6) Ganim, Z.; Tokmakoff, A. Biophys. J. 2006, 91, 2636. (7) Cowan, M. L.; Bruner, B. D.; Huse, N.; Dwyer, J. R.; Chugh, B.; Nibbering, E. T. J.; Elsaesser, T.; Miller, R. J. D. Nature 2005, 434, 199. (8) Brixner, T.; Stenger, J.; Vaswani, H. M.; Cho, M.; Blankenship, R. E.; Fleming, G. R. Nature 2005, 434, 625. (9) Jansen, T. L. C.; Knoester, J. J. Chem. Phys. 2007, 127, 234502. (10) Jansen, T. L. C.; Zhuang, W.; Mukamel, S. J. Chem. Phys. 2004, 121, 10577. (11) Abramavicius, D.; Zhuang, W.; Mukamel, S. J. Phys. Chem. B 2004, 108, 18034. (12) Redfield, A. G. AdV. Magn. Reson. 1965, 1, 1. (13) May V.; Kuhn, O. Charge and Energy Transfer Dynamics in Molecular Systems; Wiley-VCH: Weinheim, Germany, 1999. (14) Szo¨cz, V.; Palszegi, T.; Lukesˇ, V.; Sperling, J.; Milota, F.; Jakubetz, W.; Kauffmann, H. F. J. Chem. Phys. 2006, 124, 124511. (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: River Edge, NJ, 1995. (18) Jansen, T. L. C.; Hayashi, T.; Zhuang, W.; Mukamel, S. J. Chem. Phys. 2005, 123, 114504. (19) Kim, Y. S.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 1185. (20) Sˇanda, F.; Mukamel, S. J. Chem. Phys. 2006, 125, 014507.

(21) Zheng, J.; Kwak, K.; Asbury, J.; Chen, X.; Piletic, I. R.; Fayer, M. D. Science 2005, 309, 1338. (22) Davydov, A. S. Theory of Molecular Excitons; McGraw-Hill: New York, 1962. ˇ a´pek, V. Organic Molecular Crystals: Interaction, (23) Silinsh, E. A., C Localization and Transport Phenomena; AIP Press: New York, 1994. (24) van Amerogen, H.; Valkunas, L.; van Grondelle, R. Photosyntetic Excitons; World Scientific: Singapore, 2000. (25) van Kampen, N. G. Stochastic Processes in Physics and Chemistry; North Holland: Amsterdam, The Netherlands, 1992. (26) Mukamel, S. Principles of Nonlinear Optical Spectroscopy, Oxford University Press: New York, 1995. (27) Jonas, D. M. Annu. Phys. ReV. Chem. 2003, 54, 425. (28) Khalil, M.; Demirdo¨ven, N.; Tokmakoff, A. Phys. ReV. Lett. 2003, 90, 047401. (29) Sˇanda, F.; Mukamel, S. Phys. ReV. Lett. 2007, 98, 080603. (30) Sˇanda, F.; Mukamel, S. J. Chem. Phys. 2007, 127, 154107. (31) As noted in ref 30, for a single chromophore, eq 22 may be generalized beyond the stochastic model. For a multiexcitonic but excitonconserving Hamiltonian, the four diagrams of Figure 1 contribute pSA(ω3, t2 ) ∞, ω1) ) I(ω1)[I(ω3) + IE(ω3)], where IE is the emission line shape obtained from the excited-state equilibrium Fν, assuming the asymptotic Green’s function (of gg and ee manifold) approach [G]νν′(tf∞) ) FνTrν′ in some suitable representation, where ν is a 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, R. R.; Bodenhausen, G.; Wokaun, A. Principles of Nuclear Magnetic Resonance in One and Two Dimensions; Oxford University Press: New York, 1987. (34) Pisliakov, A. V.; Mancˇal, T.; Fleming, G. R. J. Chem. Phys. 2006, 124, 234505. (35) Alicki, R.; Lendi, K. Quantum Dynamical Semigroups and Applications; Springer: Berlin, Germany, 1987. (36) Risken H. The Fokker-Plank Equation; Springer: Berlin, Germany, 1989. (37) Kjellberg, P.; Bru¨ggemann, B.; Pullerits, T. Phys. ReV. B 2006, 74, 024303. (38) Engel, G. S.; Calhoun, T. R.; Read, E. L.; Ahn, T. K.; Mancˇal, T.; Cheng, Y.-C.; Blankenship, R. E.; Fleming, G. R. Nature 2007, 446, 782. (39) Kim, H. D.; 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, J. N.; Ambegaokar, V. J. Chem. Phys. 1985, 83, 4491. (42) Zusman, D. Chem. Phys. 1980, 49, 295. (43) Chernyak, V.; Mukamel, S. J. Chem. Phys. 1996, 105, 4565. (44) Ishizaki, A.; Tanimura, Y. J. Chem. Phys. 2006, 125, 084501.

JP801457C