Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Chem Rev. Author manuscript; available in PMC 2010 November 8.
Published in final edited form as:
PMCID: PMC2975548

Coherent Multidimensional Optical Spectroscopy of Excitons in Molecular Aggregates; Quasiparticle versus Supermolecule Perspectives

1. Introduction

Molecular aggregates are abundant in nature; they form spontaneously in concentrated solutions and on surfaces and can be synthesized by supramolecular chemistry techniques.13 Assemblies of chromophores play important roles in many biological processes such as light-harvesting and primary charge-separation in photosynthesis.4 Aggregates come in various types of structures: one-dimensional strands (H and J aggregates),2,5 two-dimensional layers,68 circular complexes,9 cylindrical nanotubes, multiwall cylinders and supercylinders,1012 branched fractal structures, dendrimers, and disordered globular complexes.13 They also can be fabricated on substrates by “dip-pen” technology.14 Figure 1 presents some typical structures of pigment-protein complexes found in natural light-harvesting membranes. These are made of chlorophyll and carotenoid chromophores held together by proteins.4,9,15,16

Figure 1
Some biophysical photosynthetic systems: FMO = Fenna-Matthews-Olson protein, a trimer, each unit has 7 chromophores; LH2 = light-harvesting complex 2, double-ring structure of chromophores, 27 chromophores; PS1 = photosystem 1, a complex of 96 chromophores; ...

We consider aggregates made out of chromophores with nonoverlapping charge distributions where intermolecular couplings are purely electrostatic. The optical excitations of such complexes are known as Frenkel excitons.1720 Applications of molecular exciton theory to dimers4,2124 show some key features shared by larger complexes. Their absorption spectra have two bands, whose intensities and Davydov splitting are related to the intermolecular interaction strength and the relative orientation.

Excitons in large aggregates can be delocalized across many chromophores and may show coherent or incoherent energy transfer. The optical absorption, time-resolved fluorescence, and pump-probe spectra in strongly coupled linear J-aggregates have been extensively studied.2,5,18,2529 These spectra contain signatures of cooperative optical response: the absorption splits into several Davydov sub-bands and is shifted compared to the monomer. Additionally, the absorption and fluorescence bands are narrower in strongly coupled J aggregates because of motional and exchange narrowing that dynamically average these properties over the inhomogeneous distribution of energies.30 Aggregates may show cooperative spontaneous emission, superradiance, which results in shorter radiative lifetime than the monomer.3136 Large molecular aggregates may show several optical absorption bands12 whose shapes and linewidths depend on their size and geometry. They undergo elaborate multistep energy-relaxation pathways that can be monitored by time-resolved techniques. These pathways are optimized in natural photosynthetic antennas to harvest light with high speed and efficiency.4,3741 Exciton energy and charge transport, as well as relaxation pathways, are of fundamental interest. Understanding their mechanisms may be used toward the development of efficient, inexpensive substitutes for semiconductor devices.

Nonlinear optical four-wave mixing (FWM) techniques have long been used for probing inter- and intramolecular interactions, excitation energies, vibrational relaxation, and charge-transferpathwaysinmolecularaggregates.Pump-probe, transient gratings, photon-echo, and time-resolved fluorescence were applied to study excited states, their interactions with the environment, and relaxation pathways.4,4252

Optical pulse sequences can be designed for disentangling the spectral features of the coherent nonlinear optical signals. Multidimensional correlation spectroscopy has been widely used in NMR to study the structure and microsecond dynamics of complex molecules.53 It has been proposed to extend these techniques to the femtosecond regime by using optical (Raman) or infrared (IR) lasers tuned in resonance with vibrations.5456 The connection with NMR was established.5759 Multidimensional infrared spectroscopy monitors hydrogen-bonding network and dynamics, protein structures, and their fluctuations.6067 Development of these IR techniques was reviewed recently67 and will not be repeated here.

Multidimensional techniques allow one to study molecular excitons in the visible region and reveal couplings and relaxation pathways. These applications were proposed in refs 6870 and demonstrated experimentally in molecules40,41,65,71,72 and semiconductor quantum wells.7375 Laser phase-locking during excitation and detection is required in these experiments, which measure the signal electric field (both amplitude and phase), not just its intensity. All excitation pulses as well as the detected signal must have a well-defined phase. This review covers multidimensional techniques carried out by applying four femtosecond pulses, as shown in Figure 2, and controlling the three time intervals, t1, t2, and t3, between them. In practice the t3 information is usually obtained interferometrically in one shot rather than by scanning t3. Spectral dispersion of the signal with the local oscillator gives the Fourier transform with respect to t3. Two-dimensional (2D) spectra are displayed by a Fourier transform of these signals with respect to a pair of these time variables. NMR signals do not depend on the wavevectors since the sample is much smaller than the wavelengths of the transitions. Different Liouville space pathways are then separated by combining measurements carried out with different phases of the pulses, rather than by detecting signals in different directions. This phase-cycling technique, which provides the same information as heterodyne detection,5759 has been first used in the optical regime using a collinear pulse configuration.76 Direct heterodyne detection by wave-vector selection was reported in ref 77. By controlling the time-ordering of incoming pulses and the signal wavevector (ks), we obtain a wealth of spectroscopic information about the system. Multidimensional techniques can target a broad variety of physical phenomena. Visible pulses probe electronically excited-state dynamical events: ultrafast intramolecular and intermolecular dephasing and relaxation, energy transport, charge photogeneration, and recombination.76,77 2D techniques can also eliminate certain types of inhomogeneous broadenings, show two-exciton resonances, and reveal couplings and correlations between chromophores.24,40,41,45,65,70,72,7881 Extensions to the UV8284 may target backbone transitions of proteins and DNA. In the future, X-ray attosecond techniques may reveal electronic wavepackets with high temporal and spatial resolution.8587

Figure 2
Coherent third-order nonlinear optical experiment. The four laser pulses are ordered in time; the signal is generated in the phase-matching direction k4. Data processing of time-domain signals and their parametric dependence on the three delays t1, t ...

Multidimensional techniques could further utilize the vector nature of the optical field by selecting specific polarization configurations. These may lead to pulse sequences sensitive to structural chirality (hereafter denoted chirality induced, CI)79,8892 These are 2D extensions of the 1D circular dichroism spectroscopy (CD).93 CI optical signals can be employed to probe specific correlations and couplings of chromophores. These have been demonstrated for probing protein structure in the infrared.94 Numerical sensitivity analysis algorithms and pulse-shaping and coherent-control techniques can help dissect and analyze coherent spectra and simplify congested signals.9597

This review surveys the broad arsenal of theoretical techniques developed toward the description of nonlinear optical signals in molecular complexes. By treating the system-field coupling perturbatively, the signals are expressed in terms of response functions, which allow a systematic classification and interpretation of the various possible signals.

The optical properties of molecular aggregates may be described by the Frenkel exciton model. This model has been first applied to molecular crystals17,98 and subsequently extended to aggregates.4,5,2729,99 The system is partitioned into units (chromophores) with nonoverlapping charge distributions; electron exchange between these units is neglected. The direct product of eigenstates of isolated chromophores forms a convenient basis set for the global excited states. In the Heitler-London approximation, the aggregate ground state is given by the product of ground states of all chromophores.17,99 Single-exciton states are formed by promoting one chromophore to its excited state, keeping all others in their ground state. Their number is equal to the number N of chromophores. Double- and higher-excited states are created similarly. The response is formulated in Liouville space in terms of the system’s density matrix, which allows one to incorporate energy dissipation due to interaction with the environment.52

The methods used for computing the optical response of aggregates may be broadly classified into two types.100 In the supermolecule or Sum Over States (SOS) approach, the response function is expanded in the global eigenstates. Optical spectra are interpreted in terms of transitions between these states. Feynman diagrams, which describe the evolution of the molecular density matrix (Liouville space pathways, LSP), are the key tool in the analysis.52 These provide a direct look at the relevant dynamics at each time interval between interactions with the fields. The signals are interpreted in terms of the relevant density matrix elements of the super-molecule and the sequence of transitions between eigenstates.

The LSP can be divided into two groups, depending on whether the density matrix during the second interval t2 is in a diagonal state (populations) or off diagonal (coherences).81 These groups have symmetry properties associated with permutations of pulse-polarization configurations that lead to distinct signatures in the multidimensional spectrum.

The necessary calculations of excited states make the supermolecule approach computationally expensive. Calculating the global eigenstates is not always feasible. Moreover, sums over states do not offer a simple physical picture when many states are involved.

In the alternative quasiparticle (QP) description, the aggregates are viewed as coupled, localized electronic oscillators. Two-exciton eigenstates are never calculated explicitly; instead, two-exciton resonances are obtained via exciton scattering and calculated using the nonlinear exciton equations (NEE).68,69,101107 Two-exciton propagators are calculated using the QP scattering matrix by solving the Bethe Salpeter equation. The lower computational cost, stemming from the more favorable scaling with size, makes this approach particularly suitable to large aggregates. The dominant contributions to the scattering matrix can be identified a priori by examining the exciton overlaps in real space, providing an efficient truncation strategy.

The spectral lineshapes of 2D signals contain signatures of interactions with phonons, vibrations, and the solvent. Slow and static fluctuations may be treated by statistical averaging. Fast fluctuations can be easily incorporated through population relaxation and dephasing rates. The intermediate fluctuation-time regime requires more careful attention. Many techniques exist for the modeling of quantum dissipative dynamics.108111 These treat the coupling of excitons to a bath at different levels of sophistication. In the simplest scheme, bath fluctuations are assumed to be fast, and second-order perturbation theory for the coupling is used to derive a Pauli master equation in the Markovian limit. Dephasing and energy-transfer processes enter as simple exponential decays. When exciton transport is neglected, energy (diagonal) fluctuations coming from a bath with Gaussian statistics can be exactly incorporated in the response functions using the cumulant expansion for Gaussian fluctuations (CGF).42,43,100,112 Exciton transport may be approximately incorporated into the CGF.24 A more detailed (and computationally expensive) approach is based on the stochastic Liouville equations (SLE), which explicitly include collective bath coordinates in the description.113 The SLE have been used to describe the signatures of chemical-exchange kinetics in coherent 2D signals.

The quasi-particle description of the Frenkel-exciton model is formally similar to that of Wannier excitons in multiband tight-binding models of semiconductors.108,114116 The number of variables is, however, different: since electrons and holes in the Frenkel-exciton model are tightly bound, the number of single-exciton variables coincides with the number N of sites; in contrast, the holes and electrons of Wannier excitons are loosely bound, and the number of single-exciton variables scales as ~N2. This unfavorable scaling is more severe for double-exciton states (~N2 for Frenkel and ~N4 for Wannier). Periodic infinite systems of Frenkel and Wannier excitons may be treated analytically, making use of translational invariance.

Several approaches for computing response functions and multidimensional optical signals are presented in this review. Closed expressions for the signals are derived based on both the QP representation and the supermolecule approach. We use Markovian limit with respect to bath fluctuations in the QP representation. Higher-level CGF and SLE approaches are described within the supermolecule approach. The semiclassical treatment of the molecular coupling with the radiation field whereby the classical radiation field interacts with the quantum exciton system.117 Applications to the Fenna-Matthews-Olson (FMO) photosynthetic complex made of seven coupled bacteriochlorophyll molecules illustrate the various levels of theory.

Section 2 introduces a simple multilevel model system. The response functions are derived in section 3. Section 4 develops the density-matrix formalism and derives the response functions for dissipative quantum systems. The response function theory is connected with experimental heterodyne-detected signals in section 5. Section 6 develops a microscopic model for excitonic aggregates in the quasi-particle representation. In section 7, we derive response functions for 2D signals based on the quasi-particle approach. Different models of system-bath coupling are analyzed. Exciton relaxation and dephasing rates are calculated in section 8, and section 9 provides closed expressions for multidimensional signals that include exciton population transport. Section 10 presents various applications of the quasi-particle theory to the FMO complex. In Section 11, we revisit the supermolecule approach to include slow bath fluctuations and exciton transport. In section 12, we describe a stochastically modulated multilevel system by including bath coordinates explicitly. Spectral line-shape parameters are obtained by solving the stochastic Liouville equations (SLE). Section 13 describes the response functions of an isotropic ensemble of molecular complexes (solutions, liquids). The response functions of the previous sections are extended by orientational averaging going beyond the dipole approximation to include first-order contributions in the optical wavevector. Chirality-induced and nonchiral signals are compared. Section 14 shows how coherent control pulse-shaping algorithms can be used to simplify 2D signals. General discussion and future directions are outlined in section 15.

Most mathematical background and technical details are given in the appendices. Linear optical signals are described in Appendix A. Appendix B describes different modes of signal detection of nonlinear signals. Appendix C provides the relation between a system of coupled two-level oscillators (hard-core bosons) and boson quasi-particles. In Appendix D, we present complete expressions for the response function for the model introduced in section 6. These complement the expressions of section 9. In Appendix E, we describe single- and double-exciton coherent propagation within the QP representation and introduce the scattering matrix. Appendix F derives exciton dephasing and transport rates in the real space representation. The final expressions for the system relaxation rates in terms of the bath spectral density are given in Appendix G. Exciton scattering in nonbosonic systems is described in Appendix H. Appendix I generalizes the quasi-particle description to infinite periodic systems. In Appendix J, we give expressions for the doorway and window functions using cumulant expansion technique for the supermolecule approach. Some quantum correlation functions used in section 11.3 are defined in Appendix K. Appendix L presents closed expressions for orientationally averaged signals.

2. Supermolecule Approach; Coherent Optical Response of Multilevel Systems

Conceptually, the simplest approach for describing and analyzing the optical response of a molecular aggregate is to view it as a supermolecule, expand the response in its global eigenstates, and add phenomenological relaxation rates. This level of theory will be described in this section. In the coming sections, we shall rederive and generalize these results using microscopic models for the bath. The alternative, quasi-particle approach developed in section 6 offers many computational advantages and often provides a simpler physical picture for the response. The supermolecule and quasiparticle approaches will be compared.

When viewed as a supermolecule, the aggregate is a multilevel quantum system described by the Hamiltonian,52




Here, ħεa is the energy of state a and Ĥ′ represents the dipole interaction with the external optical electric field E,


where r is the position of the molecule and


is the dipole operator. μab = Σαqαleft angle bracketa|rα|bright angle bracket is the transition dipole between states a and b where α runs over molecular charges qα (electrons and nuclei) with coordinates rα, relative to the center of charge.

The quantum state of an ensemble of optically driven aggregates is described by the density matrix:


The diagonal elements represent populations of various states, while off-diagonal elements, coherences, carry phase information.

The density matrix satisfies the Liouville-von Neumann equation,


Here An external file that holds a picture, illustration, etc.
Object name is nihms212366ig1.jpg, the Liouville superoperator, is given by the commutator, with the Hamiltonian An external file that holds a picture, illustration, etc.
Object name is nihms212366ig1.jpg[rho with circumflex]= [Ĥ, [rho with circumflex]]. Similar to eq 1, we partition the Liouville operator as


where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig2.jpg represents the isolated system and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig45.jpg[rho with circumflex] = [Ĥ′, [rho with circumflex]] is the interaction with the external field.

The retarded Green’s function (forward propagator) describes the free evolution of the molecular density matrix between the interaction events. Setting An external file that holds a picture, illustration, etc.
Object name is nihms212366ig45.jpg = 0, we get


where θ(t) is the Heavyside step-function, defined by θ(t)=tdτδ(τ). The density matrix of the driven system (eq 6) can be calculated perturbatively in An external file that holds a picture, illustration, etc.
Object name is nihms212366ig45.jpg by iterative integration of eqs 6 and 7. This yields52


Equation 9 provides the order-by-order expansion of the density matrix in the field and can be recast as [rho with circumflex](t) = [rho with circumflex](0)(t) +[rho with circumflex](1)(t) + [rho with circumflex](2)(t) +... . [rho with circumflex](n)(t) is the density matrix to nth order in the field. The n’th-order induced polarization is the quantity of interest in spectroscopy since it is a source of the signal field. It is given by the expectation value of the dipole operator


The induced polarization is given by P(t) = P(1)(t) + P(2)(t) +.... This polarization is a source of a new optical field. The generated field is calculated for a simple sample geometry in Appendix B. The signal optical field obtained by mixing of incoming laser fields characterized by wavevectors k1, k2, ... may be detected at certain directions ks = ± k1 ± k2 ± ... .

The linear response function connects the linear polarization with the field:52


Since P(1) and E are both vectors, R(1) is a second-rank tensor. For clarity, we avoid tensor notation through most of the review. It will be introduced only when necessary, in section 13. Expressions for the signals related to the linear response function are summarized in Appendix A.

The linear response function is calculated by substituting eq 9 into eq 10:52


where [rho with circumflex]0 is the equilibrium density matrix and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg [rho with circumflex] [equivalent] [[mu], [rho with circumflex]] is the dipole superoperator.

3. Response Functions of a Molecular Aggregate

Hereafter, we consider a molecular aggregate made of three-level molecules, whose exciton level scheme is shown in Figure 3..44 This model will be introduced microscopically in section 6. The eigenstates of this model form distinct manifolds (bands). The three lowest manifolds are g, the ground state; e, the single-exciton states; and f, the double-exciton states. The dipole operator only couples g to e (μeg) and e to f (μfe). The number of chromophores (and single-exciton states) is denoted by N. The linear response depends only on the e manifold, whereas the third-order response involves both the e and f bands.

Figure 3
Excitonic aggregate made of coupled three-level chromophores. εn is the excitation energy of each chromophore and 2(εn + Unnnn) is its double-excitation energy. The chromophores are quadratically coupled by Jmn (and quartically by Umnkl ...
Figure 4
Feynman diagrams for the signal generated in the direction kI = −k1 + k2 + k3; ESE = excited-state emission, ESA = excited-state absorption, and GSB = ground-state bleaching. Transport is not included in these diagrams (the coherent limit). Population ...

The linear response function can be calculated by expanding eq 12 in eigenstates,


where the equilibrium density matrix is [rho with circumflex]0 = |gright angle bracketleft angle bracketg|, ωeg = εeεg, and we have added phenomenological dephasing rate γeg to represent the decay of coherence. Dephasing will be treated microscopically for a model of fast bath fluctuations in section 8. Note that μge=μeg.

Higher-order response functions may be defined and calculated in the same way. Since the second-order response vanishes for our model because of the dipole selection rules, we shall focus on the third-order signals.

The third-order polarization induced in the system by the optical field is related to the incoming electric fields by52




is the third-order response function.

This expression represents a sequence of three interactions with the incoming fields described by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg, three free propagation periods between interactions (An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg), and signal generation described by the last [mu]. Since An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg is a commutator and [mu] can act either from the left or the right, eq 15 has 23 = 8 contributions known as Liouville space pathways.

We shall consider time-domain experiments, where the optical electric field consists of several pulses:


Here, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig5.jpg (tτj) is the complex envelope of the jth pulse, centered at time τj with carrier frequency ωj and wavevector kj. The indices uj = ±1 represent the positive (uj = +1) and the negative (uj = −1) frequency components of the field. Note that Ej=[Ej+]. The uj variables allow a compact notation for the response functions. We assume that the pulses are temporally well-separated: the pulse with wavevector k1 comes first, followed by k2 and finally k3.

The field is generated in the directions ±k1 ±k2 ±k3. The detector selects a signal in one direction. To calculate the optical signals, we introduce the wavevector-dependent polarization Pks


where ks = u1k1 + u2k2 + u3k3 is the signal wavevector, and the sum runs over the possible choices of uj = ± 1. We assume resonant excitation and make the rotating-wave approximation (RWA), where only low-frequency terms in eq 14 (where the field and molecular frequency subtract ωjωeg) are retained. High-frequency ωj + ωeg terms make a very small contribution to the polarization and are neglected. This is an excellent approximation for resonant techniques. By combining eq 17 with eqs 14 and 16, we get


with the signal frequency given by ωs = u1ω1 + u2ω2 + u3ω3. The wavevector-dependent response function, Rks(3), represents the signal generated in the ks direction. For short pulses, the integrand is finite when t3tτ3, t2τ3τ2, and t1τ2τ1 Therefore, the integral needs to be performed only in the narrow time window specified by pulse envelopes. When the pulses are much shorter than their delays, the pulse delays control the free propagation intervals between interactions described by the Green’s functions.

There are four independent third-order techniques, defined by the various choices of uj:





For our excitonic level scheme and its dipole selection rules, the kIV signal vanishes within the RWA and will not be considered any further.

The resonant third-order response may be interpreted using the Feynman diagrams shown in Figures 46. These represent the sequences of interactions with the various fields as well as the state of the excitonic density matrix during the intervals between interactions.52 The two vertical lines represent the ket (left) and the bra (right) of the density matrix. Time flows from bottom to top; the labels on the graph indicate the density-matrix elements during the free evolution periods between the radiative interactions. An interaction with the optical field is marked by a dot. Each interaction is accompanied by a transition dipole μ factor. The density matrix elements can only change by an interaction with the field. Wavy arrows pointing to the left indicate an interaction with the negative frequency field component An external file that holds a picture, illustration, etc.
Object name is nihms212366ig6.jpg exp[−ikjr], and those to the right indicate interaction with the positive frequency An external file that holds a picture, illustration, etc.
Object name is nihms212366ig7.jpg exp[ikjr]. Arrows coming into the diagram represent absorption of a photon accompanied by a molecular transition to a higher-energy state (g to e or e to f); outgoing arrows represent photon emission and a transition to a lower-energy state (e to g or f to e). The number of interactions, p (red dots), on the right side (right line) of the diagram determines its overall sign, which is (−1)p.

Figure 6
Feynman diagrams as in Figure 5 but for the signal generated in the direction kIII = +k1 + k2k3. ESE and GSB do not contribute in this case, and the technique does not show transport. The diagrams are labeled vii and viii for future reference. ...

The kI signal is depicted in Figure 4. During t1, the aggregate density matrix is in an optical coherence with a characteristic frequency given by the energy difference between states e and g, ωge. During the second interval (t2), it is either in the ground state (gg) (ii) or in the singly excited-state manifold (e′e) (i or iii). During the third interval (t3), it again oscillates with optical frequencies ωe′g (i and ii) or ωfe′ (iii).

The various contributions to the signal are denoted ground-state bleaching (GSB), excited-state stimulated emission (ESE), and excited-state absorption (ESA),52 as marked in Figure 4. This nomenclature reflects the physical processes during the course of interactions: in the GSB pathway, the system returns to the ground-state during the second interval t2; thus, the third interaction reflects the reduction of the ground-state population, which reduces (bleaches) the subsequent absorption. In the ESE pathway, during t2, the system is in the single-excited manifold and the third interaction brings it back to the ground state by stimulated emission. The ESA pathway has the same t1 and t2 dynamics as the ESE; however, the third interaction creates a doubly excited state f.

During the intervals between interactions (t1, t2, and t3), the density matrix evolves freely (no field). For instance, in the ESE (Figure 4) diagram, after the first interaction on the right (bra), the system density matrix in the t1 interval is in the state ρge, and its evolution gives an exp(−get1) factor, where ωge = εgεe is the interband frequency. The second interaction from the left with the ket creates the density matrix element ρe′e, which then evolves as exp(−e′et2). After the third interaction (with the bra), the density-matrix element ρe′g is created and its evolution is exp(−e′gt3). The overall contribution of this pathway to the response function is


where the Heavyside step function, θ(t), ensures causality (the response only depends on the fields at earlier times), and we have added the phenomenological dephasing rate γab for the ab coherence.

Similar to the three contributions to the kI technique shown in Figure 4, kII has three contributions (Figure 5) and kIII has two contributions (Figure 6). In the next section, we provide closed expressions for all techniques by extending this model to include exciton transport. The expressions corresponding to the diagrams shown in Figures 46 may be obtained from eqs 3239 by setting the population relaxation matrix to unity, Ge2e2,e1e1(N)(t2)=θ(t2)δe2e1.

Figure 5
Feynman diagrams as in Figure 4 but for the signal generated in the direction kII = + k1k2 + k3. Population diagrams are labeled iva, v, and via, and the coherence diagrams are labeled ivb and vib.

4. Optical Response of Excitons with Transport; The Lindblad Master Equation

The time evolution of closed quantum systems can be expanded in eigenstates that evolve independently by simply acquiring exp(−jt) phase factors. Open systems (i.e., systems coupled to a bath) undergo various types of relaxation and transport processes. These are commonly described by the evolution of the reduced density matrix of the system, where the bath variables have been projected out.111,118,119

The dynamics of the reduced density matrix may be described by a quantum master equation (QME). Its general form may be obtained directly from the quantum mechanical analogue of the classical Langevin equations of motion for open systems.111 The derivation starts by phenomenologically adding fluctuation and dissipation terms to the Schrödinger equation. This yields the SchrödingerLangevin equation,


where ψ is the wave function of the open system, V is an arbitrary set of system operators (independent of ξ) that describe its coupling with the bath, and ξα(t) are white noise variables with zero mean, left angle bracketξα(t) right angle bracket= 0, and short correlation time, ξα(t)ξβ(t)=δα,βδ(tt). The form of eq 24 guarantees conservation of the norm of the wave function when averaged over the noise.111

The QME that corresponds to this Schrödinger–Langevin equation is known as the Lindblad equation.120122


It is derived from the Liouville equation for [rho with circumflex](t) [equivalent] | ψ(t) right angle bracket left angle bracket (t)| using eq 24 followed by averaging over the noise variables. This equation may be also derived from firmer foundations that show how it preserves all essential properties of the density matrix: it is hermitian with time-independent trace (Tr([rho with circumflex]) [equivalent] 1) and positive definite.111 It will be shown in section 8 that the microscopic Redfield theory in the secular approximation can be recovered by a specific choice of Vα. The secular Redfield equations decouple the populations (diagonal elements of [rho with circumflex] in the eigenstate basis) from the coherences (off-diagonal elements). The populations then satisfy the Pauli master equation:109,111


Here, Kee, ee is the population transfer rate from state e′ into e. This is an N × N matrix with two pairs of identical indices (ee), (e′e′). In eq 26, the diagonal elements, e = e′, Kee, ee are positive, whereas the off-diagonal elements, ee′, Kee, e′e′ are negative. The rate matrix satisfies ΣeKee, e′e′ = 0 (probability conservation) and detailed-balance Ke2e2, e1e1/Ke1e1, e2e2 = exp(−[variant Planck's over 2pi]ωe2e1/(kBT)); here, kB is the Boltzmann constant and T is the temperature. These conditions are sufficient for the system to reach thermal equilibrium at long times.111

The evolution of diagonal density-matrix elements is described by the population Green’s function, ρee(t)=eGee,ee(N)(t)ρee(0), which is the probability of the transition from state e into e′ during time t. The formal solution of the Pauli master equation, eq 26, is given by


where λp is the pth eigenvalue of the left- and right-eigen equations eKee,eeχep(R)=λpχep(R) and eKee,eeχpe(L)=λpχpe(L). χ(L) χ (R) are matrices made out of the left (right) eigenvectors, and D = χ(L) χ (R) is a diagonal matrix. The eigenvectors are organized as columns of χ (R) and rows of χ (L). These two sets of eigenvectors are different since Ke′e′, ee is not Hermitian, but they share the same eigenvalues λp. χ(L) and χ(R) can be normalized as eχp,e(L)χe,p(R)=1 for each p; then χ (R) = χ (L)−1 and D is the unitary matrix.

In the secular Redfield equation, the off-diagonal elements satisfy


Its solution is ρee(t)=exp(iωeetγee(N)t), where


The first two terms represent population relaxation, and [gamma with circumflex]e,e′ is the pure-dephasing rate. For a compact representation of the response functions, we shall combine the diagonal (eq 26) and off-diagonal (eq 28) blocks into a single tetradic Green’s function representing both coherence and population relaxation:


Comparison of eq 30 with eqs 27 and 28 gives


This equation will be derived microscopically and further extended in section 8.

For the kI and kII techniques, populations are generated during the second time interval t2 in the ESA and ESE diagrams of Figures 4 and and55 when e = e′ (GSB includes ground-state population during t2; since our model contains a single ground state, no population relaxation occurs in the ground state). Population created in state e′ can relax to e during t2 according to eq 26. Coherences decay according to eq 28.

By including this relaxation model, we obtain the following sum over states (SOS) expression for the response function for the kI technique (Figures 4 and and77),

Figure 7
Feynman diagrams for the kI and kII signals, which extend diagrams ia, iiia and iva, via of Figures 4 and and55 by including population transport.




where we have introduced the complex frequencies ξab [equivalent] ωabab. Note that we do not include population relaxation between the excited-state manifold and the ground state; this would change An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg (t2) and would influence the GSB t2 dependence.52

The oscillation frequencies during t1 and t3 have opposite sense. For a two-level system (g and e), the frequency factors cancel out at t1 = t3. This eliminates inhomogeneous broadening and gives rise to the photon-echo signal,70,123 which has been used for studying dephasing of coherent dynamics in molecules and molecular complexes.124,125 For multilevel systems (with a band of e states), this cancelation will occur only for the diagonal e = e′ contributions, in the ESE and GSB pathways during t2 when transport is neglected.

The diagrams shown in Figure 7 extend Figures 4 and and55 to include population transfer during t2, as marked by the pairs of dashed arrows. For e = e′, the diagrams in Figure 7 coincide with Figures 4 and and55.

For the kII technique, we similarly have (Figures 5 and and77)


The two contributions to the kIII (Figure 6) technique are of the ESA type:



Note that exciton transport only enters kI and kII but not kIII where exciton populations are never created.

The present phenomenological model allows an easy calculation of the response functions. Exciton relaxation is described by Markovian rate equations. More elaborate decay patterns that show oscillations and nonexponential correlated dynamics often observed in molecular aggregates41,48,126 require higher levels of theory that will be presented in the coming sections.

5. Multidimensional Spectroscopy of a Three-Band Model; Heterodyne-Detected Signals

The response functions introduced earlier may be used to calculate the polarization induced by a resonant excitation. The relation between the signal and the induced polarization depends on the detection mode, as described in Appendix B. Heterodyne detection is the most advanced detection method that gives the signal field itself (both amplitude and phase).

The third-order heterodyne-detected signal is given by


where the signal depends parametrically on the delays between pulses T3 = τsτ3, T2 = τ3τ2, and T1 = τ2τ1. The notation τ [equivalent] T1, T [equivalent] T2, and t [equivalent] T3 is commonly used. Experimentally the signal field is often dispersed in a spectrometer. The entire T3 dependence is then measured in a single shot rather by scanning the delay between the third pulse and the signal. An external file that holds a picture, illustration, etc.
Object name is nihms212366ig42.jpg(t – τs) is the local oscillator field envelope used for heterodyne detection. The first three pulses are represented by eq 16, and the third-order polarization P(3)(t) is given by eq 18.

Multidimensional signals are displayed in the frequency domain by Fourier transforming Sks(3)(T3,T2,T1) with respect to the time intervals between the pulses.


Often, one uses a mixed time-frequency representation by performing double-Fourier transform. This gives, e.g., S3, T2, Ω1), etc.

Assuming that all four pulses are temporally well-separated, the integrations can be carried out using the response function eqs 3239.127 Neglecting population transport (Kee, ee [equivalent] 0), this gives the following three contributions to the kI signal (see Figures 46),




where e and e′ run over the single-exciton manifold, f runs over the two-exciton states, and η→ +0. ω1, ω2, and ω3 are the carrier frequencies of the first three pulses, and ωab and ξab were defined in sections 2 and 4. The spectral pulse envelopes


are centered around ω = 0.

For the kII signal, we similarly obtain


Finally, the kIII signal is given by




Note that resonances occur at both positive and negative Ωj.

Population relaxation may be added, as done in eqs 3239. However, the Ω2 dependence is then more complicated, and it is usually preferable to display the signal S3, T2, Ω1) in the time domain as a function of T2.

Equations 4250 show how the pulse envelopes select the possible transitions allowed by their bandwidths. The pulse envelopes serve as frequency filters, removing all transitions falling outside the pulse bandwidth. If the pulse bandwidth is larger than the exciton bandwidth, we can set &An external file that holds a picture, illustration, etc.
Object name is nihms212366ig43.jpg(ω) = 1. This impulsive signal then coincides with the response function.

This concludes the phenomenological supermolecule description of multidimensional signals in terms of the global eigenstates. In the next section, we present the alternative, quasiparticle, picture. Numerical applications of both approaches will be presented in the following sections.

6. Quasiparticle Representation of the Optical Response; The Nonlinear-Exciton Equations (NEE)

The calculation of double-exciton states requires the computationally expensive diagonalization of an An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg × An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg matrix, where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg = N (N ± 1)/2 is the number of double-exciton states, and N is the number of molecules. The quasiparticle picture avoids the explicit calculation of the two-exciton states making it more tractable in large aggregates.68,79,102,104,106108,128 This picture naturally emerges out of equations of motion for exciton variables, the NEE, that were derived and gradually developed at different levels.79,107,108 Spano and Mukamel had shown how theories of optical susceptibilities in the frequency domain based on the local mean-field approximation (MFA) can be extended by adding two-exciton variables to properly account for two-exciton resonances.102,128 The nonlinear response is then attributed to exciton–exciton scattering.104,105 The formalism was subsequently extended to molecular aggregates made of three-level molecules and to semiconductors.68,114 Additional relevant variables have been introduced to account for exciton relaxation due to coupling with phonons.68 The phonon degrees of freedom are formally eliminated; they only enter through the relaxation rates. This results in the Redfield equations for the reduced exciton density matrix.

By solving the NEE, we obtain closed Green’s function expressions for the optical response that maintain the complete bookkeeping of time ordering. Applications of this approach were made to J-aggregates,103 pump–probe, photonechoes, and other four-wave mixing techniques of light-harvesting antenna complexes.78,130,131

We first introduce the microscopic exciton model for an aggregate made out of N three-level chromophores. The special case of two-level chromophores is treated in Appendix C. Electronic excitations are expressed using the basis set of the electronic eigenstates of each chromophore: these are the ground state φm(g), the single-excited state φm(e), and the double-excited state φm(f). In the Heitler–London approximation,17 the aggregate ground state is given by a direct product of the ground states of all chromophores,


A single-excited-state basis is constructed by moving one of the chromophores into its excited state,


Double-excitations are obtained either when one of the chromophores is doubly excited


or when two chromophores are singly excited:


Overall, our model has N singly excited states, and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg = N( N + 1)/2 doubly-excited states.

We next introduce an exciton-creation operator on molecule m:131


These operators have the following properties: single-exciton states are created from the ground state Φem=B^mΦg, and double-excitons are obtained either by creating two excitations on different molecules, Φfmn;mn=B^mB^nΦg, or on the same molecule, Φfmm=21/2B^m2Φg (the √2 factor is introduced to resemble bosonic exciton properties ( B^n1=nn) within the single- and double-exciton space).

The Hermitian-conjugate, annihilation, operator will be denoted [B with circumflex]m. Operators corresponding to different chromophores commute. The commutation relations of these operators are (see Appendix C)


Only these operators are required to describe the linear and the third-order optical signals, which involve up to double-exciton states. Higher, e.g., triple, excitations, [B with circumflex][B with circumflex][B with circumflex]Φg, will not be considered.

We assume the following model Hamiltonian written in a normally ordered form (i.e., all [B with circumflex] are to the left of [B with circumflex])68


where [variant Planck's over 2pi]εn is the single-excitation energy of molecule n and quadratic coupling [variant Planck's over 2pi]Jmn is responsible for resonant exciton hopping.132134 The quartic couplings [variant Planck's over 2pi]Umn,kl represent various types of anharmonicities that only affect the two-exciton (and higher) manifolds. Terms that do not conserve the number of excitons, such as B^m and B^m2, have been neglected in this Hamiltonian (a more general Hamiltonian is discussed in ref 135).

A simplified form of eq 57,


captures the essential physical processes. The diagonal term Δm = 2Umm,mm modifies the energy of two excitations on the same site: left angle bracketΦfmm|Ĥ0fmmright angle bracket/[variant Planck's over 2pi] = 2εm + Δm. Two-level chromophores can be described by setting Δm → ∞. This prevents two excitations from residing on the same site. The couplings An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg/4 = Umn, mn [equivalent]Umn, nm shift the two-exciton energies with respect to the noninteracting excitons left angle bracketΦfmn|Ĥ0fmnright angle bracket/[variant Planck's over 2pi] = εm + εn + An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg. They can be either repulsive, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg > 0, or attractive, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg < 0.

We further assume the following form for the dipole operator where each interaction with the field can create or annihilate a single exciton:


Here, h.c denotes the Hermitian conjugate. The transition dipole left angle bracketΦg|[mu]emright angle bracket = μm gives the transition between the ground and a single-excited state and Φemμ^Φfmm=2(μm+μm(2)) is the transition between a single- and a double-excited state. Within the RWA, the coupling with the field (eq 3) is given by






are the time-dependent transition amplitudes with μ [equivalent] (μ+)*: μ annihilates the photon and, thus, is conjugated to the creation of exciton, while μ+ creates the photon accompanied by annihilation of the exciton. The total Hamiltonian is given by eq 1 together with eqs 57 and 60.

The exciton dynamics will be calculated starting with the Heisenberg equations of motion


Taking  = [B with circumflex]m and  = [B with circumflex]m[B with circumflex]n gives






and Vmn,kl = Umn,kl + Unm,kl. These are the first two members of an infinite hierarchy of coupled equations.

Different-order contributions in the field can be easily sorted out in these equations since interactions with the field change the number of excitons one at a time. Thus, μ is first order, [B with circumflex] is first order and higher, [B with circumflex][B with circumflex][B with circumflex] and μ(2)BmB^m are third order and higher, and [B with circumflex][B with circumflex][B with circumflex][B with circumflex] is at least fourth order, etc.

When the system is in a pure state, it can be described by a wave function (the expansion Φ=mcmB^m0+CmnB^mB^n0 is sufficient to represent the third-order response), and we have B^mB^n=B^mB^n. We can further factorize B^mB^nB^k=BmYnk where we used B = left angle bracketBright angle bracket and Ymn = left angle bracket[B with circumflex]m[B with circumflex]nright angle bracket. Taking expectation value of both sides of eqs 64 and 65, we obtain the NEE to third order in the field107


Quartic products [B with circumflex]†2[B with circumflex]2 contribute only to fourth and higher orders in the fields and were neglected. Within this space of relevant states, the exciton commutation relation, eq 56, can be, therefore, replaced by the boson commutation relation,


which we use in the following derivations. Equation 70 allows the creation of triple- and higher-exciton states on each chromophore; however, these lie outside of the physical space of states relevant for third-order spectroscopy. Exact bosonization procedures can be performed for arbitrary models of nonbosonic truncated oscillators99,136,137 as well as fermions.107,138,139 We review the bosonization of two-level molecules in Appendix C. The same ideas can be applied to multilevel molecules (truncated oscillators).

Different levels of approximation will be introduced for truncating the hierarchy of eqs 64 and 65. The simplest, mean-field approximation (MFA), which is equivalent to the Hartree–Fock approximation in many-electron problems114 uses the full factorization of all normally ordered products left angle bracket[B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket and left angle bracket[B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket, left angle bracket[B with circumflex][B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket.102,128 The excitons are treated as noninteracting quasiparticles, each moving in the mean field created by the others. Equation 64 is then closed, and the other equations become redundant.103 A higher-level approximation factorizes all daggered and undaggered variables left angle bracket[B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracketand left angle bracket[B with circumflex][B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex][B with circumflex]right angle bracket, while retaining the left angle bracket[B with circumflex][B with circumflex]right angle bracket variables.79,104,105 This coherent exciton dynamics (CED) takes into account two-particle statistics but neglects transport. Equations 68 and 69 then form a complete closed set. Another possible factorization left angle bracket[B with circumflex][B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex][B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket neglects double-exciton statistical properties, while maintaining transport. This level of theory is equivalent to the Semiconductor Bloch equations (SBE) with dephasing, used for Wannier excitons.114 The response functions predicted by these various levels of approximation will be discussed in the coming sections.

7. Coherent Nonlinear Optical Response of Quasiparticles in the Molecular Basis

The response functions, eqs 14 and 18, connect the induced polarization to the incoming laser electric fields. The third-order signals can be calculated using the Green’s function solution of the equation hierarchy, eqs 64 and 65 (eqs 68 and 69 as a special case). The induced polarization is defined as the expectation value of the dipole operator (eq 59).

7.1. Linear Response

In the linear regime, we get


Single-exciton dynamics when the field is switched off is described by


obtained from eq 68 with Y [equivalent]0. Propagation of excitons is described by a single-exciton Green’s function,


whose solution is G(t) = θ(t) exp(−iht), where h is the matrix with elements hmn. The solution of eq 72 is given by Bm(1)(t)=Gmn(t)Bn(1)(0).

The first-order terms in the NEE equations are obtained from


Substituting the Green’s function solution of eq 74,


into eq 71, and using eq 11, we get the linear response function:


Linear response techniques are surveyed in Appendix A.

7.2. Third-Order Response

The third-order induced polarization is


Below, we present two levels of approximation for the kI, kII, and kIII techniques: the MFA and the coherent exciton dynamics (CED) limit.

7.2.1. Mean-Field Approximation

At this level of theory, we make the factorizations left angle bracket[B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket. Equation 68 then becomes


and the polarization (eq 77) reduces to


The first-order variable is given by eq 75. The second-order variables vanish, and to third order from eq 78 we get


Here and later, summations over repeating indices are implied. Below we present simplified expressions obtained by setting μm(2)=0. The complete third-order polarization obtained by substituting eqs 75 and 80 into eq 79 is given in Appendix D.

The kI response function is generated by the VB(1)*B(1)B(1) term in eq 80 as follows: pulse 1 creates B(1)*, while the second and third pulses generate B(1) (in any order). This gives


For the kII technique, we similarly have


and for kIII,


In these expressions, the kj pulse interacts with chromophore nj, and j = 4 denotes the signal generated at chromophore n4. We thus keep track of all interaction wavevectors and polarizations.

Corresponding frequency-domain expressions of the signals are given as a special case of eqs 127129 and will be discussed in section 9.

7.2.2. The Coherent Exciton Dynamics (CED) Limit

When the system–bath coupling is neglected, the system remains in a pure state and can be described by a wave function at all times. Equations 68 and 69 then adequately describe the coherent short time dynamics.

A Green’s function for the Y variable is now required for solution of the NEE. The Y variables describe two-exciton dynamics. The corresponding Green’s function defined by Ymn(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig46.jpg(t)Ykl(0) satisfies


The solution, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(t) = θ(t) exp(−ih(Y)t), requires the diagonalization of a tetradic matrix h(Y). To avoid this diagonalization, we calculate this Green’s function using the Bethe–Salpeter equation, as described in Appendix E. In that appendix, we further introduce the exciton scattering matrix Γ, which allows one to interpret double-exciton resonances in terms of the exciton scattering.

By expanding the variables in powers of the fields, the first-order variable is given by eq 75 and the second-order variable


where we have further used the relation An external file that holds a picture, illustration, etc.
Object name is nihms212366ig47.jpg [equivalent] An external file that holds a picture, illustration, etc.
Object name is nihms212366ig48.jpg, which results from boson symmetry. Equation 68 then gives


The polarization (eq 77) is given by


For μm(2)=0, we get for the kI response function


We recast this result using the exciton scattering matrix (see Appendix E). Substituting eq 324 in eq 88 with Gmn,kl(0)(t)Gmk(t)Gnl(t), we get


Similarly, we obtain for kII


and for kIII:


Γn4n3, n2n1 is the exciton scattering matrix in two-exciton space. It is a tetradic matrix, whose indices n1 and n2 represent incoming excitons and n3 and n4 are outgoing excitons. The scattering is induced by the coupling V. The scattering matrix is given by




Closed frequency-domain expressions of the signals corresponding to eqs 8991 are obtained from eqs 127129 by neglecting transport and setting Ge4e3,e2e1(N)(ω)=δe4e2δe3e1[ω(εe2εe1)+i(γe2+γe1)]1, as discussed in section 9.

The physical picture for exciton dynamics emerging from this level of theory is very different from the MFA: excitons scatter by their mutual interactions as demonstrated by the diagrams in Figure 8 during the time interval ττ′ > 0. This scattering changes the resonant frequencies; thus, the correct double-exciton resonances are recovered.

Figure 8
Feynman diagrams representing the quasiparticle dynamics in the various techniques for μ(2) = 0 in the coherent limit (eqs 8991). Each interaction with the field is displayed by a solid dot. Time propagation is from the bottom up. The ...

8. Coupling with Phonons; Exciton-Dephasing and Transport

So far, we have expressed the third-order response in terms of exciton Green’s functions. These must be calculated in the presence of a bath to include dephasing and population relaxation. In section 4, we have introduced the bath phenomenologically using the Lindblad master equation. We now derive and extend these results using a microscopic model of the phonon bath.107

8.1. Reduced Dynamics of Excitons Coupled to a Bath

As in section 6, we consider an aggregate made out of N three-level molecules and described by the Frenkel exciton Hamiltonian:


Ĥ0 (eq 58) and Ĥ′ (eq 60) represent the isolated system and system-field Hamiltonians, respectively.

We shall use a harmonic model for the bath,


where α runs over the bath oscillators and a^α(a^α) are Bose creation (annihilation) operators that satisfy [a^α,a^β]=δαβ. Note that



where Qα and Pα are the coordinate and momentum of bath oscillator α with frequency wα. The system–bath coupling is assumed to be linear in bath coordinates


where the coupling parameters dmn represent bath-induced fluctuations of the transition energies (m = n) and the intermolecular couplings (mn). For strongly coupled systems (J [dbl greater-than sign] d), both contribute to energy-level and coupling fluctuations when transforming to the eigenstate basis.

Because of the coupling with the bath, the factorization left angle bracket [B with circumflex] [B with circumflex]right angle bracket into left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket in eqs 64 and 65 now no longer holds, and exciton dynamics for the first-to-third orders in the field depends on two additional variables, Nmn=B^mB^n and Zkmn=B^mB^mB^n. The bath-induced terms are calculated in Appendix F and result in the relaxation rates for the exciton variables. The full set of equations up to third order in the field finally reads





where we have used An external file that holds a picture, illustration, etc.
Object name is nihms212366ig49.jpg = hnnδmmhmmδnn, which is derived from the single-exciton Liouville operator An external file that holds a picture, illustration, etc.
Object name is nihms212366ig50.jpg; “[multiply sign in circle]” denotes the direct product of two vectors: [A [multiply sign in circle]B]mn = AmBn. The relaxation rate matrices are calculated using eqs 355 and 359. The rate matrix K is given by eq 362 (since ρnm = Nmn, the corresponding relaxation matrices satisfy Kmn,kl = Knm,lk).

Some useful approximations may be used for the relaxation matrices. γmn,kl(Y) will be represented approximately as γmn,kl(Y)γmkδnl+γnlδmk. For certain parameter regimes, the full Redfield relaxation operator K may lead to unphysical (negative and larger than 1) probabilities.111,140 In the next section, we convert the operator in the eigenstate basis into the Lindblad’s form, which guarantees a physically acceptable solution.

8.2. Redfield Equations in the Secular Approximation: The Lindblad Form

The secular approximation described in section 4 is widely used since it guarantees a physical behavior of the propagated density matrix.111,120,141145 We shall transform eq 101 into the single-exciton basis, i.e., the eigenstates of the h matrix, e = εeψe. This gives (the field is turned off)


Assuming that ωe2e1ωe2e1 and both frequencies are larger than values of K, the couplings between different coherences can be neglected within the RWA, and only the dephasing terms may be retained in the reduced relaxation operator. However, since ωee = 0, populations do couple. We thus obtain Ke4e3,e2e1Ke4e4,e2e2δe4e3δe2e1+γe4e3(N)δe4e2δe3e1 with γee(N)=0. For the Redfield equation in the eigenstate basis in the secular approximation, populations are decoupled from coherences. Furthermore, each coherence satisfies its own equation and is decoupled from the rest. Populations satisfy the Pauli master equation (eq 26) with population transport rates, Kee,ee. Closed expressions for the relaxation rates and for γee(N)=Kee,ee in the eigenstate representation in terms of the bath spectral density are given in Appendix G. This approximation for the relaxation operator, also known as the Davies procedure,146 guarantees that the density matrix remains positive-definite. Some examples of nonphysical behavior of the full Redfield equation are given in refs 111, 122, 140, and 142. Similar difficulties have been found for the Bloch equations, which describe the nuclear induction of a spin that interacts with a magnetic field.147 The problems associated with the Redfield equations are a consequence of the second-order perturbation theory and the Markovian approximation.140 Various schemes have been suggested to remedy this, by introducing short-time corrections to these equations.140,148,149 These were shown to give physical results for specific model systems and some ranges of parameters, where the Markovian Redfield equation breaks down. Unlike the secular approximation, these schemes do not guarantee the positivity of the density matrix. The secular approximation, when applied to the Redfield equation, brings it into the Lindblad form.111,120 The stochastic Liouville equation described in section 12 allows a more general description of the dynamics that does not rely on the secular approximation and guarantees a physically acceptable behavior.

In the eigenstate basis, the Lindblad equation (eq 25) can be written as


where we represent α by a pair of indices c and d. This will be convenient for connecting with the Redfield equation. Ve3e4(c,d)=Ψe3V^(c,d)Ψe4 where V(c,d) has been projected into the one-exciton eigenstate, Ψe = ψemΦem with Φem and ψlm defined in sections 6 and 9.

We next present the set of V’s that reproduce the Redfield relaxation rates of Appendix G. We define two types of V matrices: the first contains only a single nonzero off-diagonal element


with cd. Since the rank of the Ne2e1 matrix is N, there are N (N – 1) such V(c, d) matrices (these correspond to all off-diagonal elements of Ne2e1). The remaining N matrices are associated with the diagonal elements (c = d):


The elements a and b represent off-diagonal and diagonal fluctuations of the single-exciton Hamiltonian. They can be determined by writing equations for the population and the coherences using eqs 104, 105, and 106 and comparing with the Redfield equation within the secular approximation, eq 103. The a’s govern the population transport and are given by


To obtain the b’s, one needs to solve the equation


where [gamma with circumflex]e1,e2 is the pure-dephasing rate (as introduced in eq 29),


Multiple solutions for the b’s exist since eq 108 contains N( N – 1) equations with N2 unknowns (the N equations for e1 = e2 in eq 108 are satisfied for any choice of b’s). One particular solution is obtained from the rate expressions of Appendix G:


Finally, the frequency, ωe1e2 in eq 104 is related to ω e1e2 in eq 103 by a level-shift,


This demonstrates that the Redfield equations in the secular approximation can be recast in the Lindblad form. In general, however, the Lindblad equation (eq 104) can go beyond the secular Redfield equation and can couple populations and coherences. They always guarantee to yield a physically acceptable result.

8.3. Transport of Localized Excitons; Förster Resonant Energy Transfer (FRET)

The full Redfield Green’s function derived by a second-order expansion of the rates in the system–bath couplings, couples populations and coherences. It is invariant to the basis set used, i.e., a unitary transformation between the localized and delocalized basis sets does not affect the dynamics.144 However, because of the different way populations and coherences are treated, the secular approximation is basis-set dependent. The delocalized eigenstate basis represents strongly coupled chromophores. For weakly coupled chromophores, the aggregate eigenstates are essentially localized on individual chromophores and the real-space representation then constitutes the natural basis set, where the Redfield equations become


Making the secular approximation in this basis, the Redfield relaxation operator becomes Kmn,kl=δmnδklKmm,kk(F)+δmkδnlγmn(N) where K(F) is known as the Förster exciton-transfer rate. Note that γnn(N)=0 (see discussion following eq 103). When the intermolecular distance is large compared to molecular sizes, we can further invoke the dipole–dipole approximation for intermolecular interactions. The exciton transfer rate between molecules is then known as the Förster rate for the energy transfer between the “donor”, D, and the “acceptor”, A, molecules:4,150


where FD(ω) is the fluorescence spectrum of the donor molecule, εA(ω) is the absorption spectrum of the acceptor molecule, and VDA is the dipole–dipole coupling between the molecules:


Here, μ are the transition dipoles and RDA is the vector connecting the donor and acceptor molecules. Note that the exciton transfer rate includes the diagonal energy fluctuations. In practical applications, one uses the experimental absorption εA(ω) and emission FD(ω) spectra in the Förster formula. The strong ~R−6 dependence of K(F) is used in fluorescence resonant energy transfer (FRET) studies, to probe molecular structural fluctuations and energy transport.151,152

For small RDA, the molecular electron densities begin to overlap and energy transfer can proceed via a different, Dexter, mechanism involving the double exchange of an electron and a hole between the donor and acceptor molecules. This mechanism, which results in an exponential, ~e−αR; α > 0, dependence of the rate on the distance is particularly important for triplet excitons where the optical transition is forbidden and the Förster rate vanishes.153,154

9. Optical Response of Quasiparticles with Relaxation

We now derive closed QP expressions for the response functions in the exciton eigenstate basis. Exciton transport takes place during t2. During the delay periods t1 and t3, we only include dephasing and neglect transport. This is usually justified since population relaxation times (ps to ns) are typically longer than the interband dephasing times (100 fs).

The procedure for calculating the response function for the full set of eqs 99 and 100 is the same as in section 7.2.2: the dynamical variables are expanded to third order in the field using the necessary Green’s functions. These include the Green’s functions for the N and the Z variables: N(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig51.jpg(t)N(0) and Z(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg(t)Z(0) where N = left angle bracket [B with circumflex] [B with circumflex]right angle bracket and Z =left angle bracket [B with circumflex] [B with circumflex][B with circumflex]right angle bracket. An external file that holds a picture, illustration, etc.
Object name is nihms212366ig51.jpg is a tetradic matrix, while An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg is sextic. Note that the N variable corresponds to the transpose of the density matrix in the single-exciton space, i.e., Nmn [equivalent] ρnm. The same holds for the Green’s functions in the real-space representation (compare eq 31), so that Gmn,kl(N)=Gnm,lk(N), and so we will use only An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg. The response functions obtained by solving eqs 99 and 100 with the polarization equation (eq 77) are given by


The exciton creation and annihilation events and the propagation of exciton variables for our model with μ(2) = 0 are sketched diagrammatically in Figure 9. The excitons are created/annihilated at each solid dot, and time propagations are represented by solid arrows. Arrow-up represents G, while arrow-down stands for the conjugate propagation. Two lines within the t2 interval represent An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg (when both point up) and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg (when pointing in opposite directions). Three lines (triple propagation) within the t3 interval represent An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg. A horizontal dashed line stands for the V matrix.

Figure 9
Same as Figure 8 but including exciton transport (eqs 115117). Single solid line represents G(t) propagation, double line is either An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(when both lines point to the same direction) or An external file that holds a picture, illustration, etc.
Object name is nihms212366ig56.jpg (when the lines point to opposite directions), and triple line ...

RkI is depicted in the left diagram. The interval between the first and second interactions is described by the single-exciton variable left angle bracket [B with circumflex]right angle bracket, whose propagation oscillates with interband frequencies according to G(t1). During the second interval, the system is characterized by left angle bracket [B with circumflex] [B with circumflex]right angle bracket and its propagation given by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg(t2). This oscillates at low, intraband frequencies associated with differences between single-exciton energies on different sites and their couplings. The [B with circumflex] [B with circumflex][B with circumflex]right angle bracket variables generated during the third interval propagate according to An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg(t3t′), which again oscillate with high interband frequencies. The evolution during t2 is strongly influenced by the bath since the low, intraband oscillation frequencies are close to the bath frequencies. Therefore, the interplay of incoherent transport and dynamics will be most important during this interval as described by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg(t2). The high-frequency evolution during t1 and t3 is only weakly influenced by the bath (unless it contains resonances at optical frequencies).

Approximate factorizations will be used next to simplify the calculation of the Green’s functions An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg: by neglecting exciton population relaxation in the third interval, we can set Gkmn,kmn(Z)GkkGmn,mn.68 An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg is expressed in terms of the exciton-scattering matrix Γ (Appendix E). The secular approximation, eq 31, will be used for An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg so that exciton eigenstate populations will be decoupled from their coherences.

The Green’s function expression (eqs 115117) assumes a simpler form in the one-exciton eigenstate basis, e. The one-exciton Green’s functions become


where the single-exciton dephasing, γe, is derived in Appendix G. In the frequency domain, we have


Exciton scattering (Appendix E) is best described in the eigenstate basis.79,91,108 This leads to efficient truncation schemes of the scattering matrix, based on transition amplitudes and on the exciton-overlap integrals. The scattering picture has been applied to infinite symmetric systems and to large disordered aggregates with localized excitons.79,82,92,94,155,156

The transformation of the scattering matrix from the local to the exciton basis reads


Equation 325 expresses the scattering matrix in terms of the tetradic noninteracting-exciton Green’s function. In the eigenstate basis, the latter is diagonal,


and is related to the single-exciton Green’s functions by a convolution:


To describe the dynamics of the N variables, we transform their Green’s function to the exciton basis:


Since the population and coherence dynamics of eigenstates are decoupled in the secular approximation, this Green’s function assumes the form of eq 31.

The response functions assume a particularly compact form in the frequency domain (eq 41). For kI, we obtain from eq 115


The time-domain response functions RkI(QP)(t3,t2,t1) can be calculated by using the inverse transform (eq 317). We note that the Fourier transforms with respect to t1 and t2 only involve single Green’s functions and are trivial. We can thus derive simple expressions for mixed representation such as RkI(QP)(Ω3,t2,Ω1). When incoherent exciton transport is neglected, we can set Ge4e3,e2e1(N)(t)=δe4e2δe3e1Ge2(t)Ge1(t) and recover the coherent dynamics result (eq 89).

Starting from eqs 116 and 117, we similarly obtain the response functions for the other two techniques




Combining eqs 40 and 41 with the response functions (eqs 124126) and using the results of section 5, we obtain our final QP expressions for the multidimensional signals:





These expressions contain fewer terms than their supermolecule counterparts (eqs 4250) and allow one to make the approximations discussed above. The CED limit is obtained by setting Ge4e3,e2e1(N)(ω)=δe4e2δe3e1[ω(εe2εe1)+i(γe2+γe1)]1, which includes dephasing but no transport. The MFA is recovered by neglecting transport and using the MFA scattering matrix, Γ(MFA)( ω) = V, as discussed in Appendix E. Equations 127129 can represent nonbosonic systems by using the scattering matrix given in Appendix H. These expressions can be readily applied for infinite periodic structures; see Appendix I.

The numerical integration of eqs 40 and 18 with the response functions of eqs 124126 will be required when the optical pulses do overlap temporally or when the density matrix dynamics is nonexponential (as in section 11). Simulated finite-pulse signals for the FMO light-harvesting complex are presented in the next section.

10. Applications to the FMO Complex

Photosynthetic light-harvesting complexes found in biological membranes (see Figure 1) constitute an important class of chromophore aggregates. Photosynthesis starts with the absorption of a solar photon by one of the light-harvesting pigments, followed by transfer of the excitation energy to the reaction center where charge separation is initiated.4,150,157,158 This triggers a series of electron-transfer reactions, the redox chain, where ADP is eventually converted into ATP, which stores chemical free energy.150 The primary step, the absorption of a photon, takes place in light-harvesting antennae containing assemblies of pigment molecules that absorb light in a broad window of the solar spectrum. The transfer of this excitation toward the reaction center occurs with a remarkably high (98%) quantum yield.150 The active pigments are various types of chlorophylls (Chl) and bacterio-chlorophylls (BChl, BChla, and BChlb), which absorb light in the 600–700 nm and at longer wavelengths,159,160 and carotenoids, which absorb light at higher frequencies, in the 450–500 nm range.159,161 Most energy transferred to the reaction center is absorbed by the Chl and BChl molecules, but there is clear evidence for energy transfer from the carotenoids to the Chl molecules.162164 The carotenoids also act as photoprotective agents by thermally dissipating excess energy that could otherwise generate harmful photooxidative intermediates.165,166

Optical properties of several complexes, light-harvesting complex 1 (LH1) and 2 (LH2) and photosystem 1 (PSI) and 2 (PSII), have been extensively studied using linear and nonlinear optical techniques, revealing ultrafast exciton dynamics and transport pathways.4,37,158,167 Understanding the interaction of these complexes with light is crucial for unraveling their function and developing new generations of artificial light-harvesting and photonic devices and poses major theoretical and experimental challenges. Multidimensional spectroscopy provides an invaluable tool for unraveling and quantifying the photophysics of the photosynthetic apparatus.40,41,81 These techniques can pinpoint couplings between pigment molecules70 and can directly probe excitation energy transfer timescales.

The FMO complex168,169 in green sulfur bacteria is the first light-harvesting system whose X-ray structure has been determined (Figure 1). The FMO protein is a trimer made of identical subunits, each containing 7 BChl pigments and no carotenoids. This system has been extensively studied by 1D techniques such as absorption and linear and circular dichroism.170,171 The spectra were simulated using the exciton model (eq 58) whose parameters (site energies, εm, and interactions, Jmn) were fitted to experiment.40,170,171 Subsequent 2D optical spectroscopy40,41,81 revealed that the excitation energy is transferred toward the reaction center using specific pathways; the energy does not simply cascade stepwise down the energy ladder.40 By recording 77 K 2D spectra vs t2, Engel et al.41 observed 600 fs quantum beating in both the shape and intensity of the peaks, indicating long-lived coherences. The role of exciton delocalization and coherence in energy transport in photosynthesis is an intriguing issue: if they survive at room temperature, they could allow for fast transport and improve the efficiency of energy funneling to the reaction center. Pulse sequences that can distinguish between coherent quantum oscillations and incoherent energy dissipation during t2 will be presented in section 13.81

Multidimensional spectroscopy has also been applied to other light-harvesting systems: Lee et al. measured coherence dynamics in the bacterial reaction center,172 and Zigmantas et al. generated a 2D signal of photosynthetic antenna LH3 (27 chromophores) in purple bacteria.173 The reaction center of purple bacteria absorbs light in two well-separated bands at 800 and 850 nm. It is surrounded by a complex of BChl molecules forming a ringlike core light-harvesting antennae LH1.174,175 The membrane LH1 is surrounded by peripheral light-harvesting antennae LH2 and LH3. LH2 exists in different sizes. In Rps. acidophila, it is made of two concentric rings containing 18 and 9 BChl molecules, respectively, while in Rps. molischianum, the two rings have 16 and 8 BChl molecules, respectively.174 The two rings are connected by carotenoids.174,176 More details about the structure and the interactions of these light-harvesting complexes are given in refs 4, 174, and 177. Many theoretical studies of 1D spectroscopy of this system have been carried out. Schröder et al.178,179 tested various methods for calculating the linear absorption lineshapes, Jang and Silbey180,181 developed line-shape theories for single-molecule absorption, and Dahlbom et al.182 calculated transition-absorption of the B850 ring of LH2. The pump–probe and photon echo signals as well as superradiance were covered in our group.34,130,183185 Pullerits et al. have estimated the average delocalization length of exciton in LH2 from the simulated and experimental pump–probe spectra.186,187 They showed that the LH2 excitons are delocalized on 4 chromophores on average.

Photosystem I (PSI) is the largest photosynthetic system with known structure containing 96 BChl and 22 β-carotenoid molecules.16 Six of the BChl molecules form the electron transfer chain and constitute the reaction center. Ten “red” BChls absorb light at lower energies than the others, including the reaction center. Their precise role in energy transport is under debate.188 Photosystem II (PSII) is made of a light-harvesting antenna, LHCII, which contains 8 BChla, 6 BChlb, and 4 β-carotenoids,189 and its reaction center.190 Numerous experimental and theoretical 1D spectroscopic studies have been carried out on PSI37,38,188,191 and PSII.192195 These include absorption, linear and circular dichroism, and pump–probe.

The peak linewidths for the kI technique at t2 = 0 reveal the homogeneous and inhomonegeous dephasing rates.70,196 Variation of the cross-peak pattern with t2 allows one to monitor pigment interactions and excitation energy transfer timescales.40,81 During t2, the system’s density matrix evolves either in the first excited-state manifold, e, where population transport and coherent evolution take place, or in the ground state, as seen from the Feynman diagrams (Figure 4). Various quantum master equations for the evolution of the reduced density matrix have been applied to study energy transport in photosynthetic systems.120,178 These include the Redfield theory,120 a modified Redfield theory,185,197 and some non-Markovian variants.178,198 Other studies199 have demonstrated the differences between the Bloch model and the Redfield theory for the description of relaxation processes in LH2. Short-time, low-temperature, dynamical features of the pump–probe spectra of LH2 were explained using a simplified polaron model.200 Exciton annihilation has been shown to strongly affect nonlinear optical signals.201203

Various levels of theory developed in the earlier sections will be applied below to simulate 2D signals of the FMO complex of Chlorobium tepidum. Its chromophore configuration is sketched in Figure 10. The BChl molecules of the FMO monomer are modeled as 7 coupled chromophores. The Hamiltonian of eq 58 is used with An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg = 0 and Δ → ∞ (two-level chromophores). εm, Jmn, and μn are taken from ref 81. The overdamped Brownian oscillator spectral density (eq 372) with timescale Λ −1 = 100 fs and interaction strength λ = 35 cm−1 is used for the bath in the Markovian limit. The linear absorption, calculated using eq 269 shown in Figure 11 has five peaks labeled ae. The d and e bands contain unresolved doublets.

Figure 10
FMO complex structure of Chlorobium tepidum (left) and its simulated linear absorption (right), taken from ref 81 (eq 158).
Figure 11
Linear absorption spectrum for FMO of Chlorobium tepidum. The 5 identifiable peaks a–e, where a is the lowest-energy absorption peak and e is the highest. The power spectra of the three excitation pulses (Wb, Wd, and We) used to calculate the ...

10.1. Simulated Signals at the CED and MFA Levels

We start by neglecting transport and using the quasiparticle approach (eqs 127129) with the scattering matrix eq 387 (MFA) or eq 92 (CED). We first assume broad pulse bandwidths, setting the pulse spectral envelopes to 1.

The MFA and CED SkI signals at t2 = 0 are compared in Figure 12. Since the two approaches are identical in the single-exciton space, the diagonal section, which shows the single excitons, is similar and reflects absorption peaks. However, their amplitudes are different.

Figure 12
Im SkI(QP)(Ω3,t2,Ω1) signal (imaginary part) in the MFA and CED for FMO at t2 = 0 ps.

The large differences appear mostly in off-diagonal regions, where double-exciton properties show up. In the MFA, the two excitons evolve independently, and their correlations are neglected. Peak positions, which reflect the double-exciton energies, are then given by the combinations of single-exciton energies εe + εe′. The CED scattering matrix, in contrast, induces additional shifts of the resonant positions coming from exciton–exciton interactions.

A similar comparison is shown in Figures 13 and and1414 for the double quantum coherence SkIII signal.83,205,206 It projects the double excitons directly on the Ω2 axis, making it very sensitive to the fine features of the double-exciton manifold. The MFA provides an approximation to the double-exciton manifold by neglecting exciton–exciton correlation effects. The two peak patterns are mostly different in the higher-energy region (Ω2 >24 800 cm−1), where the MFA signal is featureless but the CED shows well-resolved peaks.

Figure 13
SkIII(QP)(t3,Ω2,Ω1) signal in the MFA and CED for FMO at t3 = 1 fs and t3 = 200 fs.
Figure 14
SkIII(QP)(Ω3,Ω2,t1) signal in the MFA and CED for FMO at t1 = 0 fs.

The remainder of this section shows simulated kI signals that go beyond the CED to include transport during t2 at the secular Redfield theory level.

10.2. Dissecting Energy-Transfer Pathways by Using Finite-Pulse Bandwidths

The 2D signals were calculated using eq 127 and include transport during t2. We used Gaussian pulse envelopes,


In Figure 15, we show three finite-pulse simulations that demonstrate how a proper choice of the pulse envelopes, using information from the absorption spectrum, can be used to reveal the energy-transfer pathways, without any prior knowledge of the system’s Hamiltonian and parameters.

Figure 15
2D signal SkI(QP)(Ω3,t2,Ω1) for the three experiments A, B, and C, as described in the text. (Top row, A) broad pulses; (Middle row, B) the first two pulses select the e band; (Bottom row, C) the first two pulses select ...

Experiment A (top row) uses four broad pulses, whose power spectra are given by Wb in Figure 11 and cover the entire exciton band. At t2 = 0 (left column of Figure 15), the kI signal is dominated by two diagonal peaks corresponding to bands b and d. The cross-peaks are weak. As t2 is increased, relaxation toward lower energies is observed; a higher-energy photon creates an excited state that decays toward lower-energy excited states, and the system emits at lower energies. At t2 = 1 ps, cross-peaks start to appear below the diagonal. The cross-peak at (−Ω1 = d, Ω3 = b) shows that some excitation energy of the d band has been transported to b. The right panel (t2 = 5 ps) shows additional cross-peaks. The (d, b) and (d, a) peaks indicate that energy from the d band has migrated to both the b and a bands. The absence of a diagonal (d, d) peak indicates that most of the initial exciton energy of the d band has decayed to lower-energy bands. The (b, a) peak is induced by energy transfer from b to a. Energy transport to the lowest excited state is completed in a few picoseconds. All energy-transfer pathways contribute to the signal; however, some are clearly dominant.

In Figure 15B, the first two pulses are narrow (We in Figure 11; σ1/2 = 50 cm−1 bandwidth) and selectively excite band e, and the last two pulses are broad (Wb). At t2 = 0 ps, the pulses select peaks such that −Ω1 lies in band e. The dominant diagonal peaks at (b, b) and (d, d) fall out of the pulse bandwidth, making weak cross-peaks more visible. At t2 = 1 ps, the population created in band e has been transferred to b, c, and d (as seen from the larger (e, d), (e, c), and (e, b) cross-peak intensities) but not to the lower a band (no (e, a) cross-peak). After 5 ps, we see four cross-peaks indicating that all lower bands are populated, but the population of the higher energy states, especially d, is very small (weaker (e, d) and (e, c) cross-peaks). It is difficult to observe this pathway with broadband pulses (Figure 15A), since it is an order of magnitude weaker than the dominant one.

In Figure 15C, the first two pulses are narrow (Wd in Figure 11) and selectively excite the d band, but the last two pulses are broad. At t2 = 0 ps, only the diagonal (d, d) peak is observed. As in Figure 15B, no significant excitation has been transferred to the lowest-energy band in 1 ps (no (d, a) cross-peak). What is different here is that only two bands are populated during the entire energy-transfer process (shown as the cross-peaks (d, b) and (d, a), the cross-peak at (d, c) never appears). The c band is not populated at t2 = 5 ps even though its energy is lower than the original excitation, suggesting that the c and d bands correspond to localized excitations residing on different (and weakly coupled) sites. Similar conclusions were reached by Brixner et al.40 who argued, by analyzing the exciton eigenvectors of a model system similar to ours, that the energy does not simply relax stepwise down the energy ladder.

10.3. Dissecting the Response Functions in Real Space

The third-order response function, whether calculated using the supermolecule or the quasiparticle approach, contains a product of four dipole moments associated with the various interactions with the fields. In order to analyze the NEE RkI(QP) response function, eq 124, we shall transform the dipoles back to real space using μe = Σnψe,n μn with the single-exciton eigenstates ψe,n defined in section 9. We further single out the n4 sum and write




RkI,n4(QP) gives the signal generated at site n4. By varying t2, we can further find the sites into which excitation energy flows, thus characterizing the population transport pathways and timescales in real space. Monitoring the energy transport in real-space is of particular interest in photosynthetic antennae, where the energy must eventually be funneled into specific sites (the reaction center).

A similar dissection procedure could be applied by singling out the third, second, or first interactions with the field (i.e., restricting the sums in eq 131 over n3, n2, and n1, respectively) and calculating the contributions of each site to the various interactions. The linear absorption signal can be similarly dissected by singling out the second interaction (see eq 269):


It gives the contribution of site n2. The dissected absorption spectrum of FMO displayed in Figure 16 shows that the lowest-energy excited state (band a) is localized on chromophore 3, band b is on chromophores 4 and 7, and band d is delocalized over chromophores 5 and 6.

Figure 16
Absorption spectrum for FMO dissected into its seven chromophores by the last interaction with the field (eq 133). The contribution from each chromophore is given by a different color.

The RkI,n4(QP)(Ω3,t2=0,Ω1) spectrum displayed in Figure 17 is dominated by chromophores 1, 4, 5, 6, and 7. As can also be seen from the absorption dissection, the (b, b) peak is delocalized over chromophores 4 and 7. The (d, d) peak is delocalized over sites 5 and 6, while (c, c) is localized on chromophore 1. For t2 = 10 ps (Figure 18), the signal is generated at chromophores 3, 4, and 7. This is to be expected since excitation energy flows in a downhill path. Band a is localized on chromophore 3. Hence, the cross-peaks (d, a) and (b, a) are dominated by RkI,3. The second lowest energy band, b, is still populated at t2 = 10 ps. Hence, RkI,4 and RkI,7 still contribute to the response function since the b band is delocalized across sites 4 and 7. These contributions contain the two peaks at (b, b) (the b band is still populated) and (d, b) (exciton transfer from d to b).

Figure 17
RkI,n4(QP)(Ω3,t2,Ω1) signal (eq 132) for the FMO light-harvesting complex at t2 = 0 ps. (Left column) The imaginary part of the response function; (Right column) the absolute value. The signal is dominated by chromophores 4, 5, 6, and ...
Figure 18
Same as Figure 17 but at t2 = 10 ps. The signal is dominated by chromophores 3 and 4 (and a small contribution of 7), onto which the lowest excited states (3) and the second lowest excited states (4 and 7) are delocalized.

This dissection analysis provides similar information to the finite pulse envelope techniques discussed above. The difference is that here we select the chromophores, whereas the finite envelopes select the exciton eigenstates. The two are close when the eigenstates are localized. One can achieve a similar dissection using a broad band for the first three pulses and a narrow band for the last pulse.

This dissection analysis provides information that goes beyond the eigenvalues and the eigenvectors: it can be used to highlight specific energy-transfer pathways and timescales (as was done in section 10.2 using finite pulses) and find the dephasing rates of coherences between chromophores. Exciton delocalization can be characterized by the participation ratio,182,207,208


Le varies between 1 and N for an N-chromophore aggregate and is a measure of the effective number of chromophores contributing to the exciton e. For FMO with the Hamiltonian parameters of ref 81, the participation ratios for the seven exciton states are (in increasing order of energy): L1 = 1.3, L2 = 2.3, L3 = 1.6, L4 = 1.8, L5 = 2.4, L6 = 2.8, and L7 = 1.6. As expected, the lowest-energy exciton is localized on a single site (chromophore number 3). The other excitons are more delocalized, especially excitons 2, 5, and 6, which are delocalized on chromophores 4 and 7 (exciton 2), 5 and 6 (exciton 5), and 5 and 6 (exciton 6). Another useful measure is the coherence length,70


which reveals the number of chromophores on which the exciton density matrix is delocalized along the antidiagonal direction.

The dissection analysis can be used to extract some additional measures of exciton-state delocalization. Let us consider the following quantity


where RkI,n4(QP)(Ω3,t2,Ω1) is the contribution of the n4th chromophore to the response function. Ld will vanish for regions dominated by a single chromophore and will grow with the degree of delocalization. In Figure 19, we display Ld at t2 = 0 and t2 = 10 ps. Obviously the diagonal peaks are localized on very few chromophores. This can be anticipated from our analysis of the linear absorption shown in Figure 16. The cross-peaks in Figure 19 for t2 = 10 ps are more delocalized because they result from population transfer facilitated by exciton delocalization (c.f., eq 370).

Figure 19
(Left column) Absolute value of the response function, RkI(QP)(Ω3,t2,Ω1) (eq 124). (Middle column) The sum of the absolute value of the dissected response function, n4RkI,n4(QP)(Ω3,t2,Ω ...

10.4. Signatures of Quartic Exciton Couplings

The quadratic excitonic couplings Jmn in eq 58 correspond to the “strong J coupling” regime in homonuclear NMR. Chromophores coupled through negative (positive) couplings form J (H) aggregates. The quartic couplings An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg are analogous to the Darling–Dennison coupling209 in infrared spectroscopy and are known as the “weak (heteronuclear) coupling” of heteronuclear NMR. A detailed comparison between NMR and multidimensional infrared and optical spectroscopy was given in refs 57 and 58. The FMO simulations presented so far assumed that each chlorophyll is a two-level chromophore, setting An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg = 0 and Δm → ∞. This anharmonic penalty prevents two excitations from residing on the same site.

Here we study the possible signatures of An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg.210 The energy of two-exciton states for excitations on sites m and n (i.e., Em,n = left angle bracketΦfmn|ĤSfmnright angle bracket, with Φ fmn defined in section 6) is


The An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg coupling implies that double excitations on neighboring sites have lower (negative An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg; exciton attraction) or higher (positive An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg; exciton repulsion) energy compared to the sum of energies of the two single excitations. In the NEE, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg enters solely through the scattering matrix. For finite An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg, the coherent time-evolution of the Y variables will show resonances shifted by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg couplings.

To demonstrate the possible effects of these off-diagonal quartic couplings, we have switched on the quartic coupling between the fourth and fifth chromophores in FMO. These sites are chosen because they contribute significantly to the strong b and d absorption bands (Figure 11), as can be easily seen from the dissected linear absorption spectrum (Figure 16). This is a simple illustration of possible signature of the An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg coupling, which we use as a free parameter. Figure 20 depicts the kI signals calculated with the NEE for An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig13.jpg = ±300 cm−1 and t2 = 0.

Figure 20
Effects of the off-diagonal quartic coupling, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg, on the kI spectra for t2 = 0. (Left column) The imaginary part of the response function; (Right column) absolute value. From top to bottom, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig13.jpg = 0, +300, −300 cm−1. Color code is the same ...

Effects of the quartic couplings are most easily explained by looking at the SOS expressions for the response function of a dimer,





where μ± are the transition dipoles between the ground and one of the single-excited states (“+” and “−”). They are related to the dipoles in the local basis by the following transformation: μ± = Ψ±,1μ1 + Ψ±,2μ2. The dimer model has a single doubly excited state, f, with energy ε+ + ε + An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg, and the transition dipole between the singly-excited and the double-excited state is μf = Ψ±,1μ2±,2μ1. The two terms of eq 140 are ESA contributions. The quartic coupling enters only in these terms.

As can be seen from eq 139, the GSB and ESE contributions to the cross-peaks interfere with the ESA. Hence, the cross-peaks at (ε+, ε) and (ε, ε+) interfere destructively. As An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg is increased, the ESA peaks are shifted along Ω3 by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg and the cross-peaks no longer cancel out, creating a doublet.

A similar splitting of the cross-peaks is expected in larger aggregates. This explains the trends seen in Figure 20. When An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = 0, the diagonal peaks (b, b) and (d, d) dominate. As An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg is turned on, the GSB and ESE contributions to the cross-peaks shift spectrally and no longer interfere with the ESA and the cross-peaks at (b, d) and (d, b) appear. Furthermore, the ESA cross-peak can be observed as weaker peaks shifted along Ω3 (these regions are circled in Figure 20). These cross-peaks appear at (b, d + An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg) and (d, b + An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg).

Signatures of the off-diagonal quartic couplings are also manifest for t2 > 0. In Figure 21, we show the kI response function for different An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg as in Figure 20 but for t2 = 10 ps. The b band is still populated even after long t2. This is evident from the large diagonal peak in the b region for An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = 0, shown in the top panel of Figure 21. Hence, the fourth site is populated since band b corresponds to an exciton localized mainly on this site. For An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = −300 cm−1, the doubly excited Φf45 state is energetically favored. The excitation energy, absorbed in the b or d band, which relaxes to b during t2, can re-emit from the d band, leading to a high-energy diagonal peak and corresponding cross-peaks. Since the doubly excited state Φf45 is blue-shifted when An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = +300 cm−1, these new features are weaker in the second row of Figure 21. In this case, the anharmonic coupling shifts the ESA contributions for the high-energy diagonal peak (d, d) and for the (b, d) cross-peak along Ω3.

Figure 21
Same as in Figure 20, but at t2 = 10 ps.

11. Interplay of Transport and Slow Fluctuations in the Third-Order Response

In section 7, we have derived quasiparticle expressions for the exciton-response functions by assuming fast fluctuations of the bath, which justifies the Markovian approximation for relaxation. This leads to a Pauli rate equation for populations and dephasing of exciton coherences. The supermolecule approach, in contrast, can describe the coupling to a bath with a broad range of fluctuation timescales. This higher-level treatment requires the explicit calculation of two-exciton eigenstates, which is only feasible in small aggregates. Fortunately, the detailed bath dynamics can be typically observed only in small systems where the electronic spectra are not too congested.

It was shown in section 3 that the resonant Liouville space pathways can be classified into ESE, ESA, and GSB types. The secular approximation for the Green’s function An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg additionally decouples coherences and populations. The response function may then be partitioned into contributions of coherences C and populations P during t2.

The resonant contributions to each technique are given by the Feynman diagrams of Figures 47. In the following expressions, we will use the same labeling for the various terms used in the diagrams: ia, ii, iiia, iva, v, and via are purely population contributions, and ib, iiib, ivb, vib, vii, and viii are purely coherence contributions.

The kI response function will be recast in the form:






As in section 4, we define the complex transition frequencies


where τa is the lifetime of state a. ϕ(t3, t2, t1) are multipoint line-shape functions induced by bath fluctuations. These are responsible for all other dephasing mechanisms beyond the lifetimes τ.

Similarly, we recast the kII response in the form:






Finally, kIII is composed of ESA1 and ESA2 (both are coherent contributions since populations are not created in kIII):



In the coming sections, we shall derive microscopically and extend all quantities appearing in the phenomenological expressions of section 4. The results of three approximations will be recast in the form of eqs 141153 and will only differ by the phase functions ϕ that will be given by eqs 166, 177, and 179. The response function expressions of section 4 are recovered by setting all ϕ = 0 and assuming that the dephasing is caused by finite lifetimes + pure dephasing in the Markovian approximation.

11.1. Coherent Exciton Dynamics; Diagonal Gaussian Fluctuations with Arbitrary Timescale

We assume that the system is described by the Hamiltonian equation (eq 94) with the system part in the form of eq 2 and the harmonic bath model given by eq 95. We partition Ĥ = Συ|υright angle bracketĤυleft angle bracketυ|, where the bath operator Ĥυ characterizes the bath-induced perturbation of the energy of system state υ,


qυυ is a collective bath coordinate:


qυυ induces the frequency fluctuations of state υ with respect to the reference Hamiltonian of state a:


Off-diagonal fluctuations responsible for population relaxation are neglected. The arguments presented here hold in the eigenstate-basis, where all off-diagonal elements of the Hamiltonian vanish. However, for coupled molecular aggregates, both diagonal and off-diagonal fluctuations in real-space representation lead to transport, since they result in nonzero off-diagonal Hamiltonian fluctuations in an arbitrary reference basis set. This model was studied in refs 100 and 112. Here, we present the final expressions for the exciton model of Figure 3. In this subsection, unlike the rest of the review, we do not make the RWA for the response function. This allows for more compact expressions. The RWA is only invoked for the final results (eq 166).

The linear and third-order response functions can be calculated exactly from eqs 12 and 15 using the second-order cumulant expansion.52,112 This is known as the CGF: cumulant expansion of Gaussian fluctuations.

The linear response function is given by


where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig14.jpg(τ2τ1) = left angle bracket[mu] (τ2)[mu] (τ1)right angle bracket is the two-point dipole correlation function. The angular brackets imply equilibrium statistical averaging over the bath. An external file that holds a picture, illustration, etc.
Object name is nihms212366ig14.jpg is calculated exactly for this model:


Here, ρa is the equilibrium population of state a, and we have introduced the line-shape function


Cvv(τ2τ1)q^vv(g)(τ2)q^vv(g)(τ1) is the correlation function of bath fluctuations whose time evolution is given by the ground-state Hamiltonian. The line-shape function can be alternatively recast in terms of the spectral density,


see eq 342 and 353. Its symmetries and other properties are summarized in Appendix F. We then have


The cumulant expansion similarly gives for the third-order response function,




is the four-point correlation function of the dipole operator in the Heisenberg picture (because we did not select pathways by invoking the RWA, R(3) is now independent of the wavevector ks). Expansion in the exciton eigenstates yields


where [var phi]dcba(τ4,τ3,τ2,τ1) is a four-point line-shape function. For our harmonic bath model, we get


where τij = τiτj. The expressions given in ref 100 were derived in the semiclassical limit, left angle bracketωab(tab(0)right angle bracket = left angle bracketωab(0)ωab(t) right angle bracket; here, we do not invoke that approximation.24 It is important to note that [var phi] does depend on the ground-state index g, even though it does not appear explicitly in the right-hand side (rhs), since all time evolutions are calculated with the ground-state reference Hamiltonian (eq 160).

By invoking the RWA, we can now select the resonant contributions in the response function, as warranted for each technique. The resulting response functions then assume the form of eqs 141153 with Gee,ee(N)(t)=θ(t)δee and no lifetime broadening τ−1 = 0. The phase functions are given by


The line-shape function can be calculated analytically for certain models of the bath spectral density. Commonly used examples are Ohmic, white-noise, and Brownian oscillator spectral densities. The overdamped Brownian oscillator spectral density, defined by


where Λ −1 is the fluctuation time scale and λ is the system–bath coupling strength, provides a simple expression for the line-shape function for t > 0 (for more details, see Appendix G),


where β = (kBT)−1 and νn = (2πn)/([variant Planck's over 2pi]β) are known as the Matsubara frequencies. For negative times, we have gvv(t)=gvv(t). The second term is a low-temperature correction, given by the additional Matsubara oscillators. The third term is responsible for bath-induced transition frequency shift. In the high-temperature limit, kBT [dbl greater-than sign] [variant Planck's over 2pi]Λ, we have


11.2. Response with Population-Transport and Intermediate-Timescale Bath Fluctuations

The results of the previous section hold for purely diagonal (frequency) fluctuations with arbitrary timescales; population relaxation caused by off-diagonal fluctuations has been neglected. The Redfield equations, on the other hand, include transport but are limited to fast fluctuations. Here, we present approximate expressions, which include exciton transport and intermediate timescale fluctuations.

Dephasing occurs during the delay times t1 and t3, whereas transport takes place during t2. Typically t1,t3 [double less-than sign]t2 (for photosynthetic complexes, population transport occurs in 500 fs-10 ps, while the coherence time scale is 100–300 fs). We assume an intermediate bath fluctuation time: fast compared to the transport but slow or comparable to the coherence-dephasing time. In this case, we can use the doorway-window representation of the response function and project the doorway and window functions onto the population/coherence blocks of the density matrix.184,185,211,212

Populations are created during t2 only in the kI and kII techniques. We therefore focus on these techniques. We first recast the response function, eq 15, in the form


where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig15.jpg(t) [equivalent] An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg (t) An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg[rho with circumflex]0 is the doorway exciton wavepacket prepared by the first two pulses, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(t) describes its propagation and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig16.jpg (t) [equivalent] P |An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(t)An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg is the window wavepacket, which represents the detection.

Both kI and kII techniques are described by three Feynman diagrams (ESA, ESE, and GSB) as shown in Figures 4, ,5,5, and and7,7, respectively. We shall separate the response function accordingly:


R(ESE) is given by diagrams i and iv, RGSB is given by ii and v, and RESA is given by iii and vi in Figures 4 and and55.

Expanding eq 170 in the system eigenstates, we get and





The superscripts of the doorway and window functions indicate whether interactions occur on the left (l, ket) or right (r, bra) side of the double-sided Feynman diagram, whereas the subscripts denote the density matrix elements. The angular brackets imply statistical averaging over bath fluctuations.

Making the secular approximation for An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg, we can separate the response function (eqs 172174), into two terms:


The coherent contribution, RC(3), only includes coherences (e1e2, e1 = e3, e2 = e4) during t2, and the entire optical process is completed before a relaxed exciton population is created. The second term ( RP(3)) only includes populations (e1 = e2, e3 = e4) and represents sequential exciton transport contributions. This representation thus keeps track on how populations are created and propagated during t2. The corresponding response functions have been formally derived using projection operator techniques;185 here, we present the final results.

We shall evaluate eq 175 using the secular Redfield relaxation operator. Rc(3) may then be calculated using the formalism of section 11.1. Rp(3) is calculated using projection operator techniques185 by factorizing the bath averages in eqs 172174 into An external file that holds a picture, illustration, etc.
Object name is nihms212366ig16.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg, and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig15.jpg factors. This gives


In the first, hopping, term, the doorway function De represents the population of the eth exciton created after two interactions with the radiation field. The hopping term includes the ESA diagrams (iiia for kI and via for kII) and ESE (ia for kI and iva for kII). The Green’s function Gee,ee(N) is the conditional probability for the exciton to hop from e to e′ during t2 obtained by solving the Pauli equation (eq 26 with the rate matrix in eq 362). The window function We represents the contribution of the eth exciton to the detected signal. The second term in eq 176 is the GSB contribution where, during t2, the system is in the ground state (diagram ii for kI and v for kII). Formal expressions for the doorway and the window functions are given in Appendix J.

The final expressions for the kI and kII response functions are given by eqs 141151. The relevant phase functions ϕ(ib), ϕ (iiib), ϕ (ivb), and ϕ (vib) (coherence pathways) are given in eq 166. The remaining phase functions ϕ (ia), ϕ (ii), ϕ (iiia), ϕ (iva), ϕ (v), and ϕ (via) (sequential pathways) are


The lifetimes (eq 146) in this case are given by τa1=bKbb,aa, and


is the bath reorganization energy.

The results in section 4 are based on the Markovian approximation for the entire response function. The present expressions treat the various time intervals and pathways differently. For the pathways that contain coherences during t2, we maintain the correlated bath dynamics by treating diagonal fluctuations exactly. Second-order perturbation theory with the Markovian approximation is used for off-diagonal fluctuations to get lifetime broadenings. The populations containing pathways during t2 ignore bath-fluctuation correlations between different propagation intervals t1, t2, and t3. During t1 and t3, the diagonal fluctuations are treated exactly (by the cumulant expansion). Again, second-order Markovian perturbation theory is made for off-diagonal fluctuations. During t2, second-order Markovian perturbation theory is performed for the coupling with all bath coordinates, which results in the secular Redfield equation.

The Redfield equations assume fast bath fluctuations and yield Lorentzian line shapes. The present expressions are more general since they allow for finite bath timescale and arbitrary line shapes. However, we assumed the complete factorization of the response function into the doorway and window functions, i.e., memory is erased during the t2 period; correlations between fluctuations during intervals t1 and t3 have been neglected. In the next subsection, we present a higher-level approximation that does retain such correlations.

11.3. Combining Slow Diagonal Fluctuations with Exciton Transport

We now derive approximate expressions for the kI and kII techniques that include both slow bath fluctuations and transport.24 kIII does not involve transport and needs not be considered here. In subsection 11.2, we assumed fluctuation timescale shorter than the time intervals t1, t2, and t3. Here, the bath fluctuations are slow compared to t1, t2, and t3.

The Feynman diagrams for kI (Figure 4) represent the time evolution of the density matrix in the coherent regime where population relaxation is neglected. Figure 7 contains additional terms that represent population transfer (ee′) during t2. The population contributions, therefore, have two terms: a coherent, population-conserving term when population during t2 does not change (Figure 4) and the population-transfer term, where populations hop from state e into e′ (Figure 7). The same partitioning also applies to kII.

We assume a bath with both slow Q(S) and fast Q(F) modes. The Markovian approximation is applied to the fast modes. During t1 and t3, these cause homogeneous line broadening (diagonal fluctuations), while during t2, they induce population relaxation and dephasing (off-diagonal fluctuations) described by the Pauli master equation (eq 26). The slow bath modes are responsible for spectral diffusion during all three intervals t1, t2, and t3, which cause correlations of density matrix between all these intervals. The diagonal parts, Qee(S), shift the system energies, while the off-diagonal parts, Qee(S), can be eliminated by diagonalizing the system Hamiltonian with Q(S) included explicitly (the adiabatic approximation). This is justified provided the fluctuation amplitude is smaller than the intraband transitions within the one-exciton band and the Q(S) timescale is longer than dephasing times γ−1. Within this model, the Pauli master equation (eq 26) still holds during t2, but the population rate matrix K is modulated by the fluctuations Q(S). We assume a smooth bath spectral density. Q(S) fluctuations then induce weak modulations of the transport rates, which will be neglected. The population Green’s function An external file that holds a picture, illustration, etc.
Object name is nihms212366ig17.jpg(t) does not depend on Q(S) and is taken out of the bath averaging in eqs 172 and 173.24 The contribution of the population-transfer pathways to the response functions presented below is derived in Appendix K.

The complete response functions for the kI and kII techniques (eq 175) contain both coherent and sequential contributions. As in section 11.2, eqs 141151 present the final expressions for kl and kll response functions. The phase functions ϕ(ib), ϕ (iiib), ϕ (ivb), and ϕ (vib) (coherence pathways) are given in eq 166, while ϕ (ia), ϕ (ii), ϕ (iiia), ϕ (iva), ϕ (v), and ϕ (via) (sequential pathways) are


where we define ζab = 1 − δab. As in eq 146, the lifetimes are given by τa1=bKbb,aa. In contrast to eq 177, these expressions include correlations between the time intervals t1, t2, and t3 in both population-conserving and population-transfer diagrams. This level of theory can reproduce some signatures of population transport during t2 such as the time-dependent fluorescence Stokes shift.24

This approach, therefore, extends the DW picture of the previous subsection. The full factorization of doorway and window functions and population propagator in the previous section implies that the density matrix phase is completely lost when the population is created, even if the population does not hop to the other exciton state. The approach presented in this section retains these phases despite population hopping.

12. Explicit Treatment of Bath Dynamics; The Stochastic Liouville Equations and Beyond

In the previous sections, bath-induced fluctuations in the exciton transport were described perturbatively. Bath variables were formally eliminated using projection operator techniques. A higher level of description is possible by explicitly including some collective bath coordinates in the supermolecule approach. This makes it possible to address a broader class of models with arbitrary bath timescales.

12.1. Stochastic Models in Liouville Space

The stochastic Liouville equations (SLE) approach provides a consistent accounting for fluctuations of arbitrary timescale by including some collective bath modes explicitly.109,213,214 The method extends the Redfield equations for exciton transport, which are limited to fast fluctuations. It further extends the CGF, which only applies to linear diagonal coupling of energies to a Gaussian bath. At high temperatures, both the CGF and Redfield approaches are obtained as limiting cases of the SLE. However, the SLE has several limitations: it neglects the system-bath entanglement, assumes an infinite bath temperature, and is computationally more expensive than the other methods. Finite temperature corrections will be discussed in section 12.4.

The stochastic model assumes that the exciton Hamiltonian (eq 57) is coupled to some classical stochastic variables σ(t), which represent the bath. The response function is calculated by averaging left angle bracket right angle bracket over the ensemble of stochastic paths


Here, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig18.jpg = [Ĥσ(t), · · · ] is the Liouville superoperator and subscript + denotes forward time ordering. It is important to note that Ĥ0;σ now depends on time only implicitly through the variable σ.

Equation 180 may be recast in the form


where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig19.jpg(ta,tb) is the Green’s function solution of the Liouville equation for free evolution (no field) between times ta and tb:


Direct implementation requires the generation of random paths σ(t), by, e.g., molecular dynamics simulations, constructing the Hamiltonian at each step, and evaluating eq 180. The quantum evolution for each path σ(t) may be described by a wave function in Hilbert space, because relaxation and decoherence effects only appear at the final averaging stage. Formally, this is obtained from eq 180 by setting


where the rhs has been recast in Hilbert space. A more economical approach that avoids the generation of ensembles of paths is possible when the stochastic bath fluctuations, σ(t), can be described by a few collective coordinates, satisfying simple Markovian dynamical rules: the probability density P(σ) of the random variables σ satisfies the master equation


where L(σ) is a linear operator. The dynamics is then fully described by the stochastic Liouville equations. The contribution to the density matrix ρ may be represented at a given time by a vector in the joint system + bath space. This density matrix has three indices: two represent the ket and the bra in Liouville space, while the third denotes the state of the stochastic variable. In the joint space, the Hamiltonian and the Liouville superoperator are time-independent (except for the external field). The evolution of the joint density matrix is described by the SLE:


An external file that holds a picture, illustration, etc.
Object name is nihms212366ig20.jpg ρ [equivalent] [Ĥσ,ρ] thus describes both the system evolution Hamiltonian and the system–bath coupling ( H^0σ=H^0+H^SB in eq 94). The system density matrix (and the response function) is finally obtained by the standard prescription of statistical mechanics: averaging over the initial and summing over the final states of the stochastic variables.

12.2. Stochastic Models for Nonlinear Spectroscopy

Since the bath dynamics is Markovian, the third-order response function may be factorized into the product of three Green’s function matrices in the joint system + bath space, each representing a single time interval, with




Here, the initial joint density matrix |ρ(0) right angle bracketright angle bracket = |ggright angle bracketright angle bracket|0right angle bracket is a direct product of the equilibrium distribution of the stochastic variables (zero right eigenvector of L(σ)) and the ground state of the system. The final summation left angle bracket0|TrS includes both tracing over the system and summation over the final bath states realized as a scalar product with the zero left eigenvector left angle bracket0| of L(σ).111

Equation 186 can be readily used to calculate the nonlinear response function; all correlations are built into the bath variables. The frequency-domain response function is obtained by a Fourier transform of the Green’s functions


The Green’s function matrices of the exciton-conserving Hamiltonian (eqs 57 and 58) are block diagonal, and each block (such as gg, eg, ee′) can be calculated separately. For each Liouville space pathway, eq 186 may be recast in terms of the matrix elements of the various block submatrices:


Equation 190 accounts for coherence transfer (e.g., gejgek) as well as exciton transport. Consequently, the Feynman diagrams (Figures 4 and and5)5) can be used, but e and f now denote the entire single-exciton and two-exciton blocks. The SLE is not limited to any particular basis.

Two types of stochastic models are commonly used in theories of spectral line shapes.215 The first, M state jump, assumes that the bath hops between M discrete states. The Hamiltonian Ĥα depends on the state α. The master equation for Pα [equivalent] P(σ = α) is


and the SLE reads


The zero left eigenvector of An external file that holds a picture, illustration, etc.
Object name is nihms212366ig20.jpg, which represents the final summation over discrete bath, is given by right angle bracket0| = (1, 1, ..., 1).

In the second model, σ is a continuous dimensionless stochastic variable Q undergoing an Ornstein–Uhlenbeck process and described by a Fokker–Planck equation216


where Λ is the relaxation rate (inverse autocorrelation time). Equation 193 describes an overdamped Brownian oscillator in the high-temperature limit.52

The Green’s function solution of eq 193 with the initial condition G(Q,Q′;0) = δ(QQ′) is


The bath density approaches the equilibrium distribution Peq(Q) = G(Q,Q′; ∞) = (2π)− 1/2 eQ2/2 at long times t [dbl greater-than sign] Λ−1. Q is, thus, a Gaussian–Markovian variable. The final summation over a continuous bath variable is represented by ∫dQ. A Gaussian process is uniquely characterized by its two-point correlation function.


All higher multipoint correlations may be calculated by using the Gaussian factorization rule:


It will be instructive to show how the CGF expressions can be recovered from the SLE. To that end, we assume linear coupling of the system to the Q coordinate


where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig2.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig21.jpg operate on the system Liouville space. We further assume noninteracting chromophores with purely diagonal bath fluctuations


The solution of the stochastic Liouville equations without the field




This integral can be recast using the Gaussian factorization rule (eq 196). To that end, we calculate the following functional that holds for arbitrary Gaussian, not necessarily Markovian, fluctuations.


Equation 199 can be solved by the second-order cumulant expansion


with the line-shape function


Equation 203 may be also obtained from eq 161 by taking the infinite temperature limit, and the following spectral density Cij(ω)=(β)diidjjωΛ/(ω2+Λ2). Comparing eqs 169 and 203, we find 2λνν′ [equivalent] (β[variant Planck's over 2pi])dννdν′ν′. At infinite temperature, the imaginary part of the line-shape function vanishes. Qualitatively new effects of the bath timescale, described by the SLE, will only appear when we go beyond eq 198 and include the coupling between excitons, add nonlinear coupling between the system and Q, or assume a non-Gaussian bath.

The third-order response functions can be obtained by using the cumulant expansion and recast in the form of eqs 162165 together with eq 203. To see that, we calculate, for instance, the ESE contribution


Making use of eqs 201 and 203, the second-order cumulant gives


This agrees with eqs 162165. The other contributions to the third-order response can be recovered in the same manner.

For interacting chromophores, the second-order cumulant expression is an approximation. However, the exact solution can be calculated from the SLE (eq 193) by expanding in the eigenbasis set of L(σ) (eq 193),


where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig22.jpg are Hermite polynomials and α = 0, 1, 2, · · · . The initial (equilibrium) bath state is represented by |0 right angle bracket= (1, 0, 0, ..., 0). The final summation assumes right angle bracket0| = (1, 0, 0, ..., 0).

The bath Liouvillian is thus diagonal [L(Q)]αα= −αΛ δα α, and the Q variable is given by a tridiagonal matrix.


Assuming linear system–bath coupling (eq 197), the stochastic Liouville equations become


By switching to Hilbert space notation, we have


Green’s function solution (eq 189) to eq 209 thus requires inversion of block tri-diagonal matrices.216 Since the Q blocks of the Green’s function are connected by raising Sα+ (lowering Sα) operators Gα±1,β=Sα±Gα,β, these can be calculated iteratively in the form of a continued fraction




The Green’s function is computed by starting with the diagonal element


For fast fluctuations |An external file that holds a picture, illustration, etc.
Object name is nihms212366ig21.jpg| [double less-than sign] [variant Planck's over 2pi]Λ, the continued fraction can be truncated at the lowest level and the stochastic variables may be eliminated; the dynamics can then be expressed in the reduced system phase space


As the fluctuations slow down, the continued fraction must be calculated to higher orders for a proper convergence.

A single Gaussian coordinate (eq 198) can describe only a special class of bath fluctuations. Multiexponential decays of arbitrarily correlated bath spectral functions Cij(t) can be reproduced by including several independent Gaussian–Markovian coordinates Qk (each has its relaxation rate Λk and coupling d(k)). An arbitrary C″ij is then decomposed into Lorentzians


The full SLE treatment (eq 208) requires no additional effort when relaxing the model of purely diagonal fluctuations, and a general four-index correlation matrix left angle bracketqij(t)qkl(0) right angle bracket can be decomposed similarly to eq 214 by allowing finite di≠j.

12.3. Application to Excitons

We now apply the SLE to the exciton model eq 57. The Redfield equations and other quantum master equations in the reduced system Liouville space hold in the limit of short bath correlation time (Λ → ∞ in eq 193) where the bath degrees of freedom may be projected out of the explicit treatment of the SLE. To compare with the Redfield equations, we focus on the single-exciton manifold (An external file that holds a picture, illustration, etc.
Object name is nihms212366ig53.jpg block). Assuming fast stochastic dynamics, we can employ eq 213. For a general coupling of the form Ĥ1 = [variant Planck's over 2pi]Σlm dlm|elright angle bracketleft angle bracketem|, we get for the relaxation term in the eigenbasis of Ĥ0


By making the secular approximation (see discussion below eq 103), we obtain the following expressions for exciton-transfer rates (ij) due to off-diagonal fluctuations (compare eqs 104 and 112):



In addition, we have dephasing rates originating from exciton transport (lifetime) and pure dephasing due to site diagonal fluctuations


The stochastic model thus agrees with the Redfield equations for fast fluctuations and infinite temperature (see Appendix G). Finite temperature corrections will be discussed below.

Note that, because the SLE is an exact representation of a physical stochastic model, it is guaranteed to yield a physically acceptable, positive definite density matrix for arbitrary values of all parameters. For the Redfield equation, in contrast, this is only guaranteed in the secular approximation.

We next demonstrate qualitatively new features predicted by stochastic models of exciton dynamics. We restrict the discussion to slow fluctuations, which show the most pronounced effects. We demonstrate how three assumptions made in the Redfield equations and the cumulant expansion may be relaxed by the SLE: Gaussian fluctuations, Markovian exciton transport, and linear system–bath coupling.

We first consider non-Gaussian fluctuations. The Kubo–Anderson two-state-jump model217219 (eq 192 with M = 2) has been applied220 to hydrogen-bonding fluctuations in the OD stretching mode of phenol in benzene,221 for which 2D line shapes show the timescale of creation/dissociation of the hydrogen bond. Below, we examine slow two-state-jump fluctuations in a dimer. We consider a model system with a few well-resolved absorption peaks. These may arise either from various exciton states or from different slowly inter-converting states of the bath. The absorption and the 2D SkI line shape of two noninteracting chromophores, undergoing two-state-jump spectral diffusion, are shown in Figure 22. In Figure 22 [top], the “up”, u (“down”, d) state is responsible for the two outer (inner) peaks. The 2D line shapes at short delay times t2 = 0 are plotted for the i, ii, iv, and v diagrams (Figure 22 [bottom]). In diagrams i, ii, and v, we see cross-peaks between different excitons, because one exciton can be annihilated and another can be created by two interactions with the laser field, and the aggregate can thus oscillate at a different frequency during t1 and t3 without any exciton transport during t2. However, no cross-peaks are observed at short delay times between various bath states before they can change. (For fast bath, motional narrowing will collapse the spectrum into a single peak.) In diagram iv, the exciton state in t1 and t3 interval must be the same and no cross-peaks are observed.

Figure 22
Linear (top) and 2D signal (bottom) for two-exciton with a two-state-jump (u, d) bath model. The signals corresponding to diagrams i, ii, iv, and v are displayed for t2 = 0, ε1u = 30Γ, ε1d = 10Γ, ε2u = −30Γ, ...

Second, the SLE can account for the finite bath fluctuation timescale. For instance, in ref 221, fluctuations originating from the complexation timescale and the 2D signals revealed the kinetic rate constants. The fluctuation timescale is crucial for some dynamical effects. For excitonic systems that possess some symmetry, the ground state may be dipole coupled to only a few eigenstates, and many levels are dark. Assuming fast spectral fluctuations (Redfield theory), these levels remain dark. A slow coordinate may, however, break these selections rules, and the dark eigenstates |ψi(Q)right angle bracket may temporarily acquire some oscillator strength, resulting in additional peaks. Similar peaks can be also reproduced by introducing static disorder on the top of the Redfield equations. However, the SLE account for the fluctuation timescales and can properly describe the dynamics of the additional peaks.222 Figure 23 illustrates this for a slow Gaussian spectral diffusion in a homodimer where the absorptive line shape (the sum of kI and kII signals)

Figure 23
1D and 2D absorptive signal (eq 219) of a model dimer with spectral diffusion described by two Gaussian–Markovian oscillators (one coupled to one site) for short t2 = 0 (left panel) and long t2 [dbl greater-than sign] Λ−1 (right panel) times. ...


is plotted for two t2 delay times. All oscillator strength is carried by the symmetric state (upper diagonal peak in Figure 23), and the only way to observe the other (lower) level is through the coupling to slow bath coordinates (lower diagonal peak in Figure 23), which disappears on the relaxation timescale.

Third, the SLE can describe nonlinear coupling of system parameters to Gaussian fluctuations. This effect is less dramatic than the other two, yet clearly observable. For instance, the quadratic expansion of interchromophore (J) coupling with stochastic motions of dihedral angles in trialanine is hardly affected by the 1D infrared absorption line shape but was necessary for the simulation of its 2D spectra.223 Additional effects, such as relaxing the Condon approximation and allowing for fluctuations of the dipole amplitude, may be readily described by the explicit representation of the bath.

The SLE show qualitative effects that are missed by the Redfield theory and can provide a test for other approximate methods, but at a higher computational cost. The dimensionality of the joint system + bath space grows exponentially with the number of collective coordinates. It is, therefore, not practical to assign an independent bath oscillator for each chromophore in large aggregates, as required for the most general description of fluctuations with an arbitrary degree of correlations. Few collective coordinates should be identified in order to keep the problem tractable.

12.4. Microscopic Bath Dynamics beyond the Stochastic Liouville Equations

The SLE, when introduced phenomenologically, can readily describe non-Gaussian fluctuations such as multistate jumps. However, they do not provide any guidance as to how to include finite temperature effects. These may be incorporated by adopting a more microscopic level of modeling of the collective coordinates.

We consider a bath consisting of a single harmonic oscillator (P,Q) coupled to many harmonic modes and described by the Hamiltonian183,224,225


This Hamiltonian is identical to eqs 95 and 98 but with a different choice of bath variables. This microscopic model for linearly coupled Gaussian fluctuations is exactly solvable at all temperatures and, therefore, allows one to discuss finite temperature corrections to the SLE. It can describe, e.g., the Stokes shift, finite temperature excitonic densities, and other temperature effects. Quartic coupling can be solved in a closed form as well, but the expressions are much more complex.226229

We introduce the bath spectral density, which describes the correlations of the force F = Σjcjqj, acting on the Q coordinates by the other bath modes.


Comparing with eq 352, we get


The noise acting on the system is described by a different spectral density C(Q)(ω), which represents the correlations of Q coordinate in the Hilbert picture:


Assuming that the microscopic bath modes q are faster than the collective coordinate Q, taking J″ (ω) = [variant Planck's over 2pi]MΩ2ω/Λ, and approximating Ω2ω2 ≈ Ω2, the bath spectral density reduces to the overdamped Brownian oscillator form:


At high temperatures [variant Planck's over 2pi] coth([variant Planck's over 2pi]βω/2) ≈ 2(βω)−1, we recover the Gaussian–Markovian stochastic process (eq 195) with (½)left angle bracket{Q (t), Q (0)}right angle bracket = (βMΩ2)−1 exp(− Λ|τ|), where {A, B}= AB + BA is the anticommutator. The Uhlenbeck–Ornstein process is a high-temperature limiting case of the model (eq 220).224 By rescaling (βM)1/2ΩQQ, to obtain a dimensionless unit equilibrium width, and setting MΩ2dmnQQAn external file that holds a picture, illustration, etc.
Object name is nihms212366ig21.jpg, we recover eq 197.

Finite temperature corrections to the SLE with the linear coupling model are proportional to the anticommutator {Ĥ1, ρ},109,183


This correction, which represents the back-action of the system on the bath, is obviously missed by stochastic models. However, its inclusion does not increase the level of complexity of the present SLE, because it retains the simple tridiagonal form of eq 208 in Q space


The resulting equations of motion


describe a noise at finite temperature, which generalizes the infinite temperature with finite fluctuation timescale of the SLE, and the finite-temperature fast fluctuations of the Redfield equations.

The finite-temperature correction (eq 226) represents bath reorganization. Consider the spectral diffusion model of purely diagonal fluctuations dmn = δmndmm (no off-diagonal coupling fluctuations). Expanding T in exciton eigenstates |emright angle bracket, it may be added to the Smoluchowski equation,230,231 yielding


For this model, the free-energy surfaces for the various states are linearly displaced by β(dmmdnn). In both states, the Smoluchowski equation represents diffusion of the bath in a potential whose equilibrium is displaced by dmm. The displacement depends on the state of the system. For the coherence block, we have


and the particle moves on a surface whose equilibrium displacement lies midway between those of levels |emright angle bracket and |enright angle bracket in eq 229.

Equation 228 holds at moderately high temperatures. At low temperatures, β[variant Planck's over 2pi]Ω > 1, the quantization of bath oscillators is necessary and higher-order terms in β obtained by expanding the coth ([variant Planck's over 2pi]βω/2) term in eq 161 must be included. This requires the introduction of additional bath oscillators with relaxation rates 2πn/β (Matsubara frequencies) and a corresponding width.183,109 These corrections have a similar structure as eq 228 but increase the computational cost.

The line-shape function for noninteracting chromophores linearly coupled to an overdamped Brownian oscillator (eqs 161 and 225) may be calculated exactly for arbitrary temperature using the second-order cumulant expression and is given by eq 168 with 2λνν′[equivalent] (β[variant Planck's over 2pi])dννdν′ν′. The imaginary part of this line-shape function represents the correction eq 228 and the Stokes shift.

For interacting chromophores, finite-temperature corrections establish thermal distribution of exciton densities, which makes this level of theory important for realistic modeling of many experimental systems in optical domain (such as FMO complex at 77 K).

The phenomenological viewpoint of the SLE model as a spectral random walk suggests another class of generalizations. The theory of random walks has many techniques and models, which go beyond ordinary Markovian diffusion.232 One possible extension is the continuous-time random walk (CTRW) model of spectral diffusion.233 This non-Markovian model assumes a waiting-time distribution function [An external file that holds a picture, illustration, etc.
Object name is nihms212366ig23.jpg(t)]lk for the random walk jump of the bath from state k to l. The model is tractable thanks to its renewal property; all memory is erased at the time of the jump. Rather than looking at the joint system + bath density matrix at a given time, we must follow the contributions to the density matrix from random walks that made a step to a given state at the specified time. This contains the information necessary to predict the future without reference to the past. In contrast, the full (unrestricted) density matrix sums over contributions where different time had elapsed from the last step. This is not sufficient for predicting the future jumps. The random walk paths for the third-order response may be sorted into eight groups, based on whether or not some step is made in any of the intervals t1, t2, and t3, and the contribution of each group to the response function must be calculated separately. In each group, we first focus on subintervals between the first and last steps in each applicable interval and calculate the propagator Z between these two times by solving the integral equation


which follows from the renewal property.

It remains to account for the evolution between the boundary steps in different intervals, during which the random walk state is fixed. The coherence evolution between the last jump at tm in the mth interval and the first jump at tl in some subsequent (lth) interval is determined by the fixed state of bath and contributes by factor


Additional factors for the first jump (that may have special waiting time distribution functions (WTDF), depending on sample preparation (initial conditions)) and for the final interval between the very last jump and the final time are constructed in the same way as in eq 232 but with special function An external file that holds a picture, illustration, etc.
Object name is nihms212366ig23.jpg. All of these contributions must be properly multiplied, averaged over initial and summed over final states, and convoluted over the possible tl and tm.

In Figure 24, we show the 2D absorptive line shape (eq 219) of a Kubo–Anderson two-state CTRW modulation of transition frequency of a single chromophore with long-tailed WTDFs An external file that holds a picture, illustration, etc.
Object name is nihms212366ig23.jpg(t) ≈ 1/tα+1 (1 < α < 2) resulting in the significant divergences δΩα −3 at fundamental frequencies of the two bath states and algebraic cross-peak relaxation proportional to 1/t2α1.234 The CTRW model further allows one to treat the aging of random walks with diverging average waiting time (0 < α < 1). Aging implies that the average mobility vanishes with time, because more and more paths are trapped in the long tails, resulting, e.g., in the anomalous relation between squared displacement and time left angle bracketx2 right angle bracket [proportional, variant]tα of the free diffusion. In addition, the individual paths σ do not uniformly track the space and the random walks are nonergodic.235 Clear signatures of aging can be seen in the third-order response function (with respect to the time interval from the start of the random walk to the first pulse).129 In Figure 25, the 2D absorptive line shape (eq 219) is shown to depend on the initial time t0 elapsed between the start of the random walk and the first interaction with laser pulse. The parameters correspond to fast bath limit at t0 = 0, with a single motionally narrowed peak. Diagonal peaks at the fundamental frequencies that grow with t0 correspond to immobilized particles and provide clear signatures of aging. Most notable is the coexistence of static diagonal peaks at fundamental bath frequencies and the central motional narrowing peak. This reflects the strongly inhomogenous (nonergodic) form of aging as described by the CTRW model. Measurements of the aging response would require the preparation of fresh samples before the application of each pulse sequence, making sure that the system has fully relaxed so that the age of the sample can be properly defined. 2D spectroscopy can thus provide a probe for glassy systems characterized by broad distributions of tunneling times.

Figure 24
Variation of the 2D absorptive signal (eq 219) for a single chromophore two-state-jump model with CTRW jump dynamics An external file that holds a picture, illustration, etc.
Object name is nihms212366ig32.jpg(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig33.jpg(t) = (θ(t))/(2π) ∫ dω eiωt[(1)/(1 + 1ω/[1 + (κα ...
Figure 25
Aging effects in the 2D absorptive signal (eq 219). The SA signal for the nonstationary random walk model An external file that holds a picture, illustration, etc.
Object name is nihms212366ig32.jpg(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig33.jpg(t) = (θ(t))/(2π) ∫ dω eiωt[(1)/(1 + (κiω)α)] at t2 = 0, κ ...

13. Nonlocal Nonlinear Response of Chiral Aggregates

Chiral molecules and aggregates can be probed by pulse sequences that yield chirality-induced higher-resolution signals. The term “chirality” has been coined more than 100 years ago by Lord Kelvin. According to his definition, “any geometrical figure, or group of points is chiral if its image in a plane mirror, ideally realized, cannot be brought to coincide with itself”. Chiral systems, therefore, come in two equally probable mirror-image configurations.236 Mirror reflection with respect to the xy plane can be described by two successive symmetry operations: (i) parity P (space inversion P (x, y, z) = (−x, −y, −z)) and (ii) coordinate system rotation by π around the z-axis. Isotropic ensembles of randomly oriented molecules are invariant to an overall rotation. A racemic mixture of chiral molecules with opposite sense of chirality is invariant to parity, whereas an unequal mixture forms an isotropic chiral ensemble. The parity operation thus transforms between isotropic ensembles with an opposite sense of chirality. Aggregate chirality can arise from both the chemical structure of individual chromophores and from a chiral arrangement (even when the chromophores are nonchiral).

Since the successive application of two parity operations restores the system back into its original state, optical signals (like all system properties and physical observables) can be classified as either parity-odd (PFo =Fo) or even (PFe = Fe).237,238 The former change sign when the chirality is reversed and, therefore, must vanish for nonchiral systems and racemates (equal mixtures of chiral molecules with opposite chirality). Parity-even signals, in contrast, are not sensitive to the sense of chirality. We shall distinguish between two types of signals: chirality-induced (CI) and nonchiral (NC).

Circularly-polarized light is chiral since it has either left (L) or right (R) screw symmetry with respect to the propagation direction. The simplest chirality-induced linear optical technique, circular dichroism (CD), measures the difference between absorption of L and R circularly polarized light.93,239,240 Thanks to its high structural sensitivity, the technique has been extensively applied to structure determination using electronic transitions in the visible241,242 and the UV243245 or vibrational transitions in the IR.246248 Raman optical activity (ROA) is a closely related technique that measures the difference between resonant scattering of L and R circularly polarized light.237,249,250

CI techniques have been extended to the nonlinear regime. The chirality of liquid and solid interfaces where the surface layer of molecules has some degree of order is commonly studied using second-order techniques251 such as sum-frequency generation (SFG) and second-harmonic generation (SHG).252,253 CI tensor elements of the second-order susceptibility have been determined for specific molecular geometries.254256 Second-order techniques may also be used to study chirality in bulk samples; however, such signals are very weak due to phase-mismatch.89,238,257,258

In the bulk, third-order signal sensitivity to system chirality requires going beyond the dipole approximation. Chirality induced (CI) two-dimensional pump–probe, photon-echo, double-quatum correlation signals, have been simulated for vibrations in chiral polypeptides and electronic spectra of chromophore aggregates.79,81,82,9092,94,97,156,259263

Geometry enters optical signals through the transition dipoles and higher multipoles. All our applications so far were based on the dipole approximation for the coupling with the field. The molecule was assumed to be small compared to the wavelength of the light and can be treated as a point particle, neglecting the variation of the phase of the field across the molecule. This may not be generally justified for large aggregates. To take these phase shifts into account, we must go beyond the electric dipole approximation and include the system’s interaction with both electric and magnetic optical fields. This requires the introduction of a large number of electric, magnetic, and mixed-response functions. A more compact but equivalent description can be developed by describing the system–field interaction using the minimal coupling (pA) Hamiltonian,


where A is the vector potential of the electromagnetic field.


is the electric current operator, and


is the charge density (rα and pα are the coordinate and momentum operators of electron α). The pA Hamiltonian has two coupling terms, which are linear and quadratic, respectively, in the field. For a time-domain impulsive experiment where each of the incoming laser fields interacts only once with the system and the fields do not overlap temporally, the AA term can be neglected. We thus adopt the simplified Hamiltonian of eq 1 with the system–field interaction:


Similar to eq 16, we shall expand the vector potential in modes


where Aj(uj) is the vector potential envelope corresponding to pulse j, centered at τj, with wavevector kj and carrier frequency ωj. Both electric and magnetic optical fields are represented by the vector potential since B = [nabla] × A and the Maxwell equations give E = − Å. We then have B = iΣj kj × Aj and E =iΣj ωjAj, where Aj is the vector potential of the jth pulse.

We now generalize eq 14 by introducing the nonlocal response function for the induced current,


where tj are the delays between successive interactions with the fields and the third-order response function R(3) now relates the current density to the vector potentials. From the Maxwell equations and the charge continuity equation, it follows that the current is related to the polarization P and magnetization M vectors:


We consider an isotropic ensemble of aggregates. We define the center of charge of the mth molecule rm and the displacement ρm [equivalent] rrm relative to the rm. The total current operator is given by a sum of operators localized around the chromophore charge centers, J(r) = Σm δ(rρmrm)Jm(ρm). Each |ρm| < l0, where l0 is the molecular size because Jm(ρm) = 0 for |ρm| ≥ l0. Performing the Fourier transform Jm(k) = ∫ dr eik(ρ), we can recast eq 238 for current amplitudes in a form similar to eq 18:


Here, ks = u1k1 + u2k2 + u3k3 and ωs = u1ω1 + u2ω2 + u3ω3. We shall make the RWA for all interactions with the k1, k2, and k3 fields. Equation 240 is the spatially nonlocal generalization of eq 18.

The current-density operator for the m’th molecule will be written in the form


Fourier transform with respect to r, Jm(k) = ∫ dr exp(−ikr)Jm(r), provides a direct connection with the optical signals. When the aggregate size is small (but not negligible) compared to the optical wavelength, we can expand the transition current densities to first order in the wavevectors k, thereby connecting them to the transition dipole, μ, the transition quadrupole q, and the magnetic transition dipole m for chromophore m:


All molecular properties have been calculated with respect to the origin Rm. The kRm contribution is required for molecular aggregates, where each chromophore has its own reference point. Because of exciton delocalization, this enters as an additional contribution of the single-exciton eigenstate quadrupole (a collection of dipoles has higher multipoles: quadrupole, etc.). The total quadrupole contribution can then be represented as a sum of two terms: the transition intramolecular quadrupole amplitude of each chromophore nψneqn(i) and a nonlocal interchromophore contribution originating from spatial distribution of chromophores of the aggregate Σn ψneRnμn, where Rn is the location (origin) of dipole μn. Thus, the total contribution of the e exciton to the quadrupole operator is


We note that, while in the local molecular basis, we have qn(i)ν2ν1=qn(i)ν1ν2; this is no longer the case for the delocalized excitons basis, where qeν2ν1qeν1ν2, due to the different origin in the calculation of the quadrupole moment of different molecule. A common model for molecular aggregates assumes that each chromophore is nonchiral and is described by its local dipole, neglecting the intramolecular quadrupole and magnetic contributions.4,241 Only the interchromophore quadrupole contribution is then taken into account.

The functional form of the transition current operator (eq 241) resemble our previous definition of the polarization operator (eq 59 with μ(2) = 0). All calculations of the response function, therefore, proceed as before, except that the transition dipole expectation value must be replaced with the induced current density. The linear response function (eq 76), for example, is now given by


The QP third-order response function in the coherent regime for the kI =−k1 + k2 + k3 (eq 89) technique similarly reads






is the orientationally averaged product: Jj(k) is the transition current density of transition j, kj is the jth-field wavevector, and νj is its polarization. The summation over molecules describes an isotropic orientationally averaged ensemble.

The induced currents in isotropic ensembles (eq 247) can be expressed in terms of orientationally averaged products of the various multipoles:


13.1. The Dipole Approximation

The first term in eq 248 gives the orientationally-averaged response functions in the dipole approximation (k = 0)


This product of four transition dipoles is invariant to the parity transformation and is thus independent of system chirality. The resulting response functions and all signals derived from them are nonchiral (NC). The following relation holds between current-density/vector potential response function and the induced-polarization/electric field response function: −i4 ωsu3ω3u2ω2u1ω1R(3) = R(3).

Orientational averaging of the transition dipoles performed using eq 424 in Appendix L requires the calculation of the vectors F{ν}(4) in eq 429. The elements of this vector define the three linearly independent configurations of transition dipoles, F{ν}(4)=(1,0,0)T, or F{ν}(4)=(0,1,0)T, or F{ν}(4)=(0,0,1)T. Assuming that all laser beams propagate along z, the independent polarization configurations, which survive the orientational averaging, are listed in Table 1; the labels in (ν4ν3ν2ν1) are ordered chronologically: ν1 is the polarization vector of the first pulse, ν2 is the second, etc.

Table 1
Linearly Independent NC Third-Order Experiments for kI, kII, and kIII Techniques (There Is No Restriction on the Wavevector Configuration in This Case)

We next go beyond the dipole approximation, where the response functions become wavevector-dependent.

13.2. Beyond the Dipole Approximation; Collinear Chiral Techniques

We assume that all laser beams propagate along z and have the same carrier frequency. This configuration has been used experimentally.76 The various phase-matching signals can be separated using phase-cycling (combining various experiments defined with different phases of the fields) as is commonly done in NMR.58,264 Note that (u1, u2, u3) = (−1, 1, 1) for kI, (1, −1, 1) for kII, and (1, 1, −1) for kIII. From eq 248, we get


where k = |k| is the wavevector amplitude and εγβα is the antisymmetric Levi–Civita tensor (εxyz = εyzx = εzxy = 1, any other permutation of indices changes sign, and εγβα = 0 when at least two indices coincide). The field polarizations are now restricted to νj = x, y, and we assumed m = im″ to be pure imaginary, which is the case for a real basis set.

Upon examination of the vector F(5) in eq 430, as was done for the dipole approximation, we find the three independent techniques listed in Table 2. These involve products of three transition dipoles and either one quadrupole or magnetic transition dipole. Since these products change sign upon parity transformation, the signals must be induced by system chirality and we denote them chirality-induced (CI).

Table 2
Linearly Independent CI Collinear Third-Order Experiments for kI, kII, and kIII Techniques (We Assume Collinear Configuration Where All Wavevectors Are Along z)

13.3. Noncollinear Chiral Techniques

The resonant response functions in noncollinear configurations depend on both pulse polarizations and wavevectors, R (k4,k3,k2,k1; ν4,ν3,ν2,ν1). Using eq 240, we decompose all vectors into elementary components and obtain the tensor expression of the response function for an arbitrary wavevector configuration:


Here, ν4ν3ν2ν1(k4k3k2k1) defines an elementary field polarization and wavevector configuration. Using the formalism of ref 265, we find that only few of these configurations are linearly independent and should be calculated. Within the dipole approximation, we only had three linearly independent configurations. To first order in the wavevector, the column vector F{ν}(5) in eq 430 defines six linearly independent techniques. These are obtained by specifying F{ν}(5)=(0,0,0,0,0,1)T, (0, 0, 0, 0, 1, 0)T, etc. The six noncollinear wavevector and polarization configurations that lead to this F{ν}(5)156 are listed in Table 3.

Table 3
Linearly Independent CI Noncollinear Third-Order Experiments for kI, kII, and kIII Techniques with Pulse Configurations Denoted as ν4ν3ν2ν1(k4k3k2k1) and [alpha] [equivalent] −α

13.4. Chirality-Induced Signals in the FMO Complex

The simulations presented in section 10 used the quasi-particle expressions of the response function within the dipole approximation. Orientational averaging was described in section 13.1. FMO is a small complex with only 7 single-exciton states and 21 double-exciton states; the higher-level line-shape expressions of section 11 may be used in this case. Orientational averaging can be performed beyond the dipole approximation as described in sections 13.2 and 13.3.

In Figure 26, we compare NC (xxxx) and CI (xxxy) signals calculated using the doorway–window expressions of section 11.2 for different delay times t2. The xxxx (dipole approximation) is displayed in the top row. At short time (t2 ≤ 500 fs), it shows strong blue (negative) features along the diagonal (GSB and ESE contributions). Positive (green) peaks above the diagonal show ESA. Slow redistribution of peaks observed at longer picosecond timescales reflects population transport to lowest-energy state. The second row depicts the corresponding xxxy collinear CI signal. At short delay times, the fine structure of the spectrum is much more detailed than in the NC signal. The peaks along the diagonal now resemble the CD spectrum, while off-diagonal features provide information about couplings. The variation with t2 can now be followed more closely since separate well-resolved peaks appear in addition to peak shoulders seen in (xxxx). These simulations use the realistic line width.

Figure 26
kI signal of the FMO photosynthetic complex calculated using the CGF approach of section 11.2, with population transport for various t2 (broad band pulses with uniform An external file that holds a picture, illustration, etc.
Object name is nihms212366ig43.jpg(ω) spectral shape were assumed; the signal coincides with the response function). ...

Signals specifically designed to probe coherent and dissipative dynamics81 are shown in Figure 27. The combination B1=Sxyzx(zxx¯z¯)(3)Sxyxx(zxzx)(3) selectively probes exciton-intraband coherences and their dynamics, since the contributions from populations exactly cancel out. On the other hand, the signal C1=Sxxyz(zzxx)(3)Sxyzx(zxx¯z¯)(3) highlights population dynamics: since exciton-coherence contributions from the ESE diagrams are canceled by quantum interference. Both are CI signals that capture structural information with high sensitivity. Since population contributions have been canceled, the off-diagonal peak signatures are now much better resolved and show the coherent exciton evolution. This is connected to coherent quantum effects in the system such as exciton delocalization, excitonic nature, and relaxation mechanisms.

Figure 27
Optimized kI signals B1 and C1 (for details, see the main text) of the FMO complex calculated using the same approach as in Figure 26. B1 shows only contributions from excitonic coherences, and C ...

The C1 signal shows a much slower dynamics. At zero delay, the signal is simply the inverse of B1, but the differences rapidly develop at longer delays. The C1 peak structure is richer than B1, and specific peaks can be attributed to coherences in ESE and ESA processes by comparing both B1 and C1: the peaks that coincide in B1 and C1 are from ESA, and the ones missing in C1 are from ESE. Peaks in C1, which are absent in B1, come from populations. Note that coherences decay within 150 fs. Thus, the longer time dynamics is now solely due to populations. The C1 signal vanishes in the absence of population transfer in the ESE diagram. Thus, the peaks in C1 are particularly sensitive to the population-transport pathways (note that ground-state populations in GSB do not evolve in our model and only give static contributions).

14. Manipulating 2D Signals by Coherent-Control Pulse-Shaping Algorithms

In complex biological systems such as photosynthetic aggregates, coherent and dissipative dynamics occur on similar ultrafast timescales and line broadening leads to congested, weak cross-peaks. Global data analysis may be used to extract dynamical information out of time-resolved signals.50 The technique allows one to unravel spectra of excited species covered under the background of other signals. This method has been used for data analysis of ultrafast pump–probe (transient absorption) spectroscopy to characterize energy-transfer dynamics in light-harvesting antenna,47,266 to identify a photoprotective energy dissipation mechanism in higher plants,51 and for other applications. The interpretation is, however, model-dependent.

The active control of ultrafast processes can provide detailed information on the underlying vibrational and electronic dynamics.267270 Coherent-control pulse-shaping techniques have been used to drive quantum systems into a desired state,271275 to manipulate excitons in multidimensional spectroscopy,96,97,262,263,276282 and to control the function of biological systems.283286

A number of NMR control techniques involving simultaneous and shaped pulses, composite pulses, refocusing schemes, and effective Hamiltonians have been proposed.287 These developments are expected to find future applications in quantum information and computation. The same ideas may be used to design and simplify 2D signals of excitons.

2D spectroscopy of photosynthetic excitons provides valuable insights into coherent and dissipative excited-state dynamics.40,41,80,173,288,289 Specific energy-transfer pathways may be directly observed through the temporal evolution of cross-peaks, and their oscillations provide information about the long-lived electronic quantum coherences. These provide signatures of quantum effects in primary biological events.290 Analysis of these 2D spectra is complicated by spectral overlaps and weak cross-peaks. Specific pulse polarization sequences designed to highlight the off-diagonal spectral features were described in section 13.81,288 In spite of this progress, resolving weak congested cross-peaks and providing general methods for simplifying 2D optical spectra as in 2D NMR remain an open challenge.

Below, we demonstrate how pulse-shaping and coherent-control algorithms may be used to simplify the multidimensional spectra of photosynthetic excitons, revealing information about dynamical correlations between exciton states. In the weak-field (third-order) regime considered in this review, the shaped laser pulses do not modify the response functions of the system but allow one to selectively amplify and resolve the targeted spectral features by manipulating the constructive and destructive interferences among Liouville space pathways.

We shall represent the optical electric field in the form


Compared to eq 16, here we have added the temporal phase function ϕjv (tτj)68,291 and introduced the pulse polarizations ν explicitly.

The electric field can be alternatively recast in the frequency domain


where Ejν+(ω)exp(iϕjν(ω)) is the Fourier transform of the complex envelope function Ejν+(t)exp(iϕjν(t)) and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig24.jpg = [An external file that holds a picture, illustration, etc.
Object name is nihms212366ig25.jpg]*. Note that An external file that holds a picture, illustration, etc.
Object name is nihms212366ig26.jpg (ω) and [phi with tilde](ω) are not the Fourier transforms of An external file that holds a picture, illustration, etc.
Object name is nihms212366ig43.jpg(t) and ϕ (t).

The signal (eq 40) will now be recast in the form


The frequency-domain signals may be obtained using eq 41. The nonlinear response contains a 4-fold summation over molecular eigenstates. These terms interfere. Other types of interferences arise from the multiple integrations over pulse envelopes An external file that holds a picture, illustration, etc.
Object name is nihms212366ig44.jpg(t) and the superposition of various tensor components. Later in this section we demonstrate how the pulse phases and envelopes may be manipulated in order to control the interferences among pathways and simplify the 2D signals.

14.1. Phase, Amplitude, and Polarization Pulse Shaping

The field of coherent control with laser pulse shaping has grown in many directions since the original proposals in 1980s292295 and the first experimental realizations in the early 1990s.296299 Shaped femtosecond laser pulses provide numerous “control knobs” for manipulating the outcome of light-matter interactions. Experimentally this is achieved with pulse shapers, which can independently modify the phase, amplitude, and polarization of each frequency component of a laser pulse. Various pulse-shaper designs have been introduced to create arbitrary pulse profiles over broad spectral ranges.300304

It is possible to perform single-parameter coherent-control protocols, where a pulse width, chirp, subpulse separation, or other single control knobs are varied. One can also perform free optimization control by modulating a large number of parameters simultaneously. Pulse-shaping may be performed in the time or frequency domain, by controlling either the spectral (eq 253) or temporal (eqs 252) phase and amplitude profiles. Shaping the polarization state of the electric field (degree of ellipticity and orientation angle) provides a different class of control parameters.275,305308 This technique, pioneered by Brixner, Gerber, and co-workers, opens up a new avenue of control based on manipulating vector properties of light propagation with promising applications to time-dependent chirality.309 In this method, the spectral phases of two orthogonal polarization components are modulated, and their interference leads to complex polarization-shaped field profiles. Various pulse-shaper designs provide full control of the polarization of ultrashort laser pulses by means of phase and amplitude shaping.310312

Early work had focused on the a priori design of optical pulse shapes and was limited to simple model systems and a few control parameters. In 1992, Judson and Rabitz proposed an iterative-control scheme, where the optimal pulse is found in a self-learning, adaptive loop that does not require any prior knowledge of the system.313 This paved the way for practical applications to complex systems with multiple control parameters. The approach involves the following steps: (i) generating an input trial laser pulse, (ii) applying the pulse to the sample and observing the signal, and (iii) using a genetic learning algorithm with a feedback loop to generate new pulse shapes based on the outcome of the previous trials.314 These steps are repeated iteratively until the desired target is met. Genetic algorithms perform a parallel search on an entire “population” of pulses, using payoff (cost function) information, rather than derivatives or auxiliary knowledge, and employ probabilistic, rather than deterministic, rules.

Adaptive control was demonstrated in 1997.315317 Various pulse-shaping schemes with genetic optimization were studied by Motzkus and co-workers.318 They investigated the stability and effectiveness of evolutionary algorithms, emphasizing the influence of steering parameters, number of configurations in search space, and noise. One of the most attractive applications of these closed-loop optimizations is the adaptive control of chemical reactions in gas and liquid phases, a subject that has been reviewed recently.269,319 Polarization pulse shaping has been applied to control the ionization of iodine320 and potassium dimers321 and was found to be more efficient than the shaping of linearly polarized pulses. Brixner and co-workers have also used polarization pulse shaping for the adaptive control of nano-optical fields.204 By controlling the interferences between near fields in nanostructures, they achieved subwavelength localization of electromagnetic energy, which can be used for nanoscopic ultrafast space-time-resolved spectroscopy.322 Silberberg and co-workers had used polarization pulse shaping to eliminate the nonresonant background and selectively excite closely lying Raman modes in coherent anti-Stokes Raman spectroscopy (CARS).307 The same technique was then applied to control angular-momentum states of atoms, demonstrating that it can access the three-dimensional vectorial properties of the system and excite states not accessible by linearly polarized control.323

Polarization pulse shaping has been applied to simulate the control of exciton localization in an ensemble of energetically disordered and randomly oriented FMO complexes278 and to control the exciton dynamics in FMO276 and in larger photosynthetic complexes.277 The excitation energy was localized on one of the chromophores by varying the spectral and temporal profiles of linearly polarized laser pulses. Linearly polarized pulses and simple forms of the spectral phase function have been applied toward the control of the energy flow in the light-harvesting complex LH2.286 Prokhorenko et al.283 varied the isomerization quantum yield of retinal in Bacteriorhodopsin by ± 20% using a genetic algorithm with a feedback loop and low-intensity shaped laser pulses.

Femtosecond pulse shaping has been applied to perform phase-coherent multidimensional spectroscopy.324327 Nelson et al. constructed a 2D femtosecond pulse-shaping apparatus325327 that provides control over the delay times, phases, and temporal shapes of all four interacting pulses and can be used for the coherent control of multidimensional spectroscopy. Zanni’s group demonstrated the use of a pulse shaper to generate a phase-stable collinear pulse pair that can be used to obtain 2D electronic spectra in a partially collinear configuration.324 These developments make it possible to combine coherent control with pulse shaping with 2D optical spectroscopy.

Here, we demonstrate how adaptive control can be used to simplify 2D optical spectra of complex systems and enhance desired spectral features by shaping the spectral, temporal, and polarization pulse profiles. Various pulse-shaping algorithms may be applied to improve the resolution of nonlinear optical spectra.96,97,279 Simulations demonstrate how the pump–probe spectrum of a model helical pentamer may be better resolved by using pure-phase polarization pulse shaping with genetic and iterative Fourier transform algorithms.279 The state of light was manipulated by varying the phases of two perpendicular polarization components of the pump, holding its total spectral and temporal intensity profiles fixed. Genetic and iterative Fourier transform algorithms were used to search for pulse phase profiles that optimize the ratio of the signal at two frequencies and allow for the observation of new features.

We adopt the description of polarization pulse shaping developed by Brixner et al.305,328 The time-evolution of the electric field vector within a single optical cycle around time t is represented by an ellipse (Figure 28B). A complete characterization of polarization-shaped laser pulses is provided by the quasi-3D representation obtained by specifying the temporal intensity ij(t), total phase ϕj(t), orientation of the ellipse θj(t), and ellipticity εj(t).

Figure 28
Phase, amplitude, and polarization pulse shaping. (A) Laser pulse sequence and time variables used in eqs 254 and 261 to describe a heterodyne-detected two-pulse 2D photon-echo experiment. (B) Polarization ellipse representation of the electric field ...

We consider a pulse propagating along z (eq 252) with two polarization components ν = x and y. The field will be described using the following variables. An auxiliary angle χj(t) ε [0, π/2] represents the ratio of the amplitudes An external file that holds a picture, illustration, etc.
Object name is nihms212366ig27.jpg(t):


A second angle δj(t) ε [-π, π] is given by the difference between the temporal phase modulations:


[theta w/ tilde]j(t) ε [-π/4, π/4] is a third angle


The orientation angle of the ellipse θ j(t) ε [-π/2, π/2] is then given by


The ellipticity angle εj(t) ε [-π/4, π/4] is


Finally, the total phase ϕj(t) of the laser-field oscillation with respect to the perihelion of the momentary light ellipse is defined as


The instantaneous pulse frequency ωj(t) is given by the time derivative of ϕj(t)


Various porphyrin structures, such as dimers, rings, tapes, wheels, and boxes, have been proposed as potential artificial light-harvesting systems.4,96,97,279,329,330 The Soret band was described by the Frenkel-exciton model using known parameters.263,331 Each porphyrin has two orthogonal electronic transitions (Figure 29A). The electronic level scheme consists of three manifolds: a ground state (g), four singly excited states (e), and six doubly excited states (f) (Figure 29B).

Figure 29
Pulse shaping control of excitons in the Soret band of a chiral porphyrin dimer. (A) Structure of the complex. (B) The exciton level scheme: the ground state, four single-exciton states, and six two-exciton states. (C) and (D) are the linear absorption ...

We shall compare four laser pulse-shaping protocols. The first (ST||) involves spectral and temporal shaping of both x- and y-components, creating parallel linearly polarized pulses. The real and imaginary components of the electric field are independently varied. The second form (ST[perpendicular]) involves spectral and temporal shaping with perpendicular linearly polarized pulses. The third (P) involves pure-phase polarization pulse shaping constructed using the iterative Fourier transform (IFT) algorithm.279,332 The temporal envelopes An external file that holds a picture, illustration, etc.
Object name is nihms212366ig28.jpg(t) and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig29.jpg (t) are varied holding the total temporal intensity profile Ij(t)=Ejx2(t)+Ejy2(t) as well as the spectral amplitudes of both electric field components An external file that holds a picture, illustration, etc.
Object name is nihms212366ig30.jpg(ω) and Ij(ω)=Ejx2(ω)+Ejy2(ω) fixed to the initial values. The temporal (ϕjv(t)) and spectral ([phi with tilde]jv(ω);) phases were calculated using the IFT algorithm. To obtain An external file that holds a picture, illustration, etc.
Object name is nihms212366ig31.jpg(t), we varied the auxiliary angle χj(t). The fourth (STP) is the most general and uses the entire parameter space involving spectral, temporal, and polarization pulse shaping.

We apply these pulse-shaping algorithms to control the amplification of weak spectral features in heterodyne-detected 2D two-pulse photon-echo (PE) signals with a Gaussian linearly polarized first pulse with wavevector k1, local oscillator ks pulse, and a second pulse k2 = k3 (Figure 28A) tuned with adaptive pulse shaping. The signal is generated in the direction ks = −k1 + 2k2 by varying t1 and t3 and holding t2 = 0. The signals were calculated by performing multiple time integrations and taking Fourier transforms with respect to the time delays between the first and second laser pulses, τ, and between the second and the local oscillator pulse, t,52


where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig31.jpg(t) are the complex laser pulse envelopes.

The simulated spectra are displayed in Figure 29E. The reference 2D PE spectrum obtained with two linearly polarized Gaussian pulses (top row) has two major diagonal peaks, D1 and D2, and weak off-diagonal peaks, C1 and C2. The major contributions to D1 and D2 are the one-exciton states at ΔΩ1 and ΔΩ3, respectively (Figure 29C). The chosen control target was to amplify the off-diagonal (C1 and C2) relative to the diagonal (D1 and D2) peaks. To specify the target, we define the integrated intensity of the jth peak.


The ratios of the integrated diagonal and cross-peaks T1 = (ID1 + ID2)/IC1 and T2 = (ID1 + ID2)/IC2 were used as cost functions and minimized using a genetic algorithm (GA). The optimal 2D PE spectra P and STP are displayed in rows 2 and 3 of Figure 29E. Peaks D1, D2, C1, and C2 are marked by circles. Arrows mark the amplified target peaks.

Complete resolution and amplification of the target peaks was only achieved using the STP (target T1) and STP′ (target T2). The resulting polarization-shaped pulses are displayed using the quasi-3D electric-field representation based on specifying the temporal intensity Ij(t), total phase ϕj(t), orientation of the ellipse θj(t), and ellipticity εj(t).305,328 In Figure 29F, we show these parameters as a set of instantaneous light ellipses at different times along the electric-field propagation axis. The ellipse sizes represent the intensity, and their shapes provide instantaneous snapshots of the polarization state. The projections represent the intensities of individual electric-field components. The variation of the total phase (chirp) is shown by color. Weak cross-peaks were enhanced compared to the intense diagonal peaks. The optimized electric fields with different polarizations select different Liouville space pathways of the tensor components of the response functions, highlighting several tensor elements with distinct features. It is clear by comparison of the four pulse-shaping strategies that modulation of the polarization profiles was necessary in order to achieve the complete resolution and amplification of the cross-peaks.97

14.2. Adaptive Optimization of Pulse-Polarization Configurations

The multiparameter coherent-control schemes described above involve elaborate pulse-shaping and characterization techniques. We now describe a different approach for simplifying 2D spectra by taking linear combinations of various tensor components of the response function.

It was shown in section 13 that the third-order 2D signals depend on various tensor components of the response function Rks4ν3ν2ν1.79 The nonchiral (NC) tensor components (dipole approximation) include ν4ν3ν2ν1 [equivalent] xxyy, xyxy, and xyyx (note that xxxx = xxyy + xyxy + xyyx). We label them T1(NC)T3(NC). The nine linearly independent chirality-induced (CI) kI = −k1 + k2 + k3 response function tensor components are labeled T1(CI)T9(CI). These were defined in Tables 13 in section 13.

We have constructed the following superposition of linearly independent response function tensor components:


where n = 3 for the NC and n = 9 the CI techniques. The complex coefficients ci=ci+ici were varied using a genetic algorithm in order to optimize the ratios of the amplitudes of selected peaks in the porphyrin dimer described above (Figure 29).263 The interference of CI tensor components was used to suppress the strong diagonal peaks and amplify weak cross-peaks (Figure 30). The initial nonchiral signals at delay time t2 = 0 were simulated with a homogeneous broadening γ = 500 cm−1. Only two diagonal peaks, P1 and P5, and two cross-peaks, P2 and P4, were resolved. In the chirality-induced T1T9 spectra, the remaining four cross-peaks P3, P6, P7, and P8 were only partially resolved. Using a genetic algorithm, we found optimal linear combinations of CI tensor components that amplify these peaks and reveal information about couplings of the x- and y-polarized transitions. This information was not available from the absorption spectra or the initial 2D NC spectra. The Soret-band in porphyrin aggregates is the strongest and is connected with interesting photophysical processes such as S2-exciton emission and dynamic intensity borrowing.333 The optimized cross-peaks provide additional information that is not revealed by superpositions of the nonchiral tensor components.

Figure 30
Coherent control of 2D signals of a porphyrin dimer by a linear combination of 2D CI tensor components. Shown are absolute magnitude at delay time t2 = 0: four nonchiral tensor components xxxx, xxyy = xyxy, xyyx; nine chirality-induced tensor components ...

We next demonstrate how coherent control may be applied to manipulate exciton-transport pathways262 in the FMO complex as revealed by 2D experiments. The simulated xxxx spectrum shown in Figure 31 (left column) fits the experimental data.40,80 Thin lines on both axes in the 2D spectra mark the exciton 1–7 positions in order of increasing energy. The cross-peaks (2,4) and (2,5) corresponding to exciton states 4 and 5 on the Ω1-axis, and state 2 on the Ω3-axis, overlap with the diagonal peaks of the same sign, forming an L-shape signal.

Figure 31
Coherent control of energy transfer in FMO by adaptive optimization of a linear combination of 2D CI tensor components. Absolute magnitude of 2D signals: nonchiral xxxx and optimal 2D CI (fast) and (slow) linear combinations corresponding to the optimizations ...

Using coherent control combined with a genetic learning algorithm, different energy-transfer channels were separated by optimizing a linear superposition of CI tensor components. The control targets were aimed at manipulating various peaks in these 2D spectra by selecting the ratios of the integrated cross-peaks in the absolute values of the 2D signals at long delay time t2 = 5 ps. Two cross-peaks (1,5) and (1,7) contributing to predicted fast and slow pathways40,80 are marked in Figure 31 by red and green circles, respectively. These cross-peaks represent energy transfer from higher to lower energy states, and grow as a function of t2. The (1,7) cross-peak was not resolved in the xxxx spectrum.

In the optimization of the fast energy-transfer pathway, the ratio of (1,5) to (1,7) is maximized. The red circle and arrow mark the cross-peaks (1,5) and (1,2), respectively. The optimal coefficients obtained from the optimizations at t2 = 5 ps were then used to obtain the 2D spectra for other t2 delay times. The exciton states 1, 2, and 5 are involved in the fast pathway, and the corresponding cross-peaks are enhanced. For the slow pathway, the ratio of (1,7) to (1,5) is maximized. Green circle and arrow indicate cross-peaks (1,7) and (1,3), respectively. As a result, both pathways are separately optimized and well-resolved. In addition, the fast pathway shows a node (a region where the signal vanishes) between the arrow and the red circle, corresponding to the elimination of the (1,3) cross-peak. Similarly, in the optimization of the slow pathway, we see nodes corresponding to the elimination of the (1,2), (1,4), and (1,5) cross-peaks of the fast pathway. Optimal laser-pulse polarization configurations obtained using coherent control and genetic algorithms enhance the chirality-induced spectral features, revealing a specific energy-transfer pathway that was not resolved in the xxxx spectra. These methods may be applied to other 2D techniques. They have also been applied, for example, in vibrational 2D spectroscopy to resolve the structure of amyloid fibrils.334

In conclusion, pulse shaping combined with coherent control can increase the sensitivity of multidimensional spectroscopy of molecular aggregates by unraveling weak cross-peaks from otherwise congested spectra.

15. Discussion and Conclusions

Femtosecond nonlinear optical spectroscopy had gradually evolved from simple pump–probe (transient absorption) applications in the early 1980s to the elaborate coherent multidimensional techniques surveyed here. The field has been driven by both ultrafast pulse-laser source technology (down to attoseconds)268 and progress in the theoretical modeling of complex molecular assemblies with dissipation.

The response-function formalism provides a unified picture of nonlinear optical techniques in systems driven by moderately strong optical fields. The signals are given as convolutions of response functions to various orders in the incoming fields, with the field envelopes. The response functions carry all necessary molecular information for calculating the optical signals.

We have surveyed the key aspects of the theory, focusing on coherent third-order heterodyne-detected signals. This holographic phase-sensitive detection gives the signal field itself (both amplitude and phase). Early third-order optical measurements used homodyne detection, which only gives the signal intensity. For completeness, we briefly describe these signals using the response-function formalism.

The simplest third-order technique, pump–probe, is performed with two optical pulses with variable delay time τ. The first pulse, the pump, perturbs the system, and the second, the probe, detects the time evolution of this perturbation. The technique can be performed in various configurations: the pump and probe can have same/different carrier frequencies, have same/different bandwidths, and be in completely different optical regions ranging from THz to the UV. Instead of using the response function, this signal can be alternatively described as an induced absorption/emission events from a nonstationary state. Using the third-order response function when the pump and the probe are temporally well-separated, the pump is associated with the first two interactions, while the probe induces the third interaction and the self-heterodyned detection. Since the first two interactions occur within the pulse envelope and are indistinguishable, the signal is the sum of the kI and kII techniques, which only differ by the order of the first two interactions (Figures 4, ,5).5). The delay time t1 is controlled by the pump-field, while t3 is controlled by the probe. The delay between pump and probe pulses then coincides with the time variable t2. The pump–probe signal is thus given by


where ωpu and ωpr are the pump and probe carrier frequencies. This gives the absorptive lineshape (eq 219).

Transient grating is another closely related technique. Here, two pump-pulses with wavevectors k1 and k2 are applied simultaneously; the probe comes after a delay τ and the homodyne signal is measured in the k1k2 + kpr direction. Usually ωpu and ωpr are tuned to the same optical transition and the 1D time-domain signal is recorded versus τ. The pulse bandwidth selects the relevant spectral region. This homodyne signal is given by


Similar to the pump-probe, this signal monitors excited-state population and relaxation.

Finally, we note the three-pulse photon-echo peak-shifts (3PEPS) often used to study bath fluctuations. It uses three pulses with the same carrier frequency and a broad bandwidth. 3PEPS measures the homodyne kI signal as a function of two delay times t1 and t2.


This signal is usually displayed as a family of 1D plots vs t1, for various values of t2. For an ensemble of two-level systems, the signal along t1 grows for short t1, reaches a maximum at t1 = τM, and eventually decays. The value of τM and its dependence on t2 reflects the bath correlation function (eq 343). Thus, the dependence τM(t2) (the 3PEPS signal) reveals direct information on the system–bath couplings.4244,48,335,336 This interpretation is limited to a single two-level chromophore.

We have introduced two types of response functions: one (R) connecting the induced polarization with the incoming electric fields (eq 14), and the other (R) relates the induced current with the vector potential of the field (eq 238). The relation between them can be established using the Maxwell equations. In the dipole approximation where the signals calculated using both approaches must coin-cide, we find that R → (ω1ω2ω3ω4)−1R, where ωj is the carrier frequency of the jth laser field. We derived R in the dipole approximation, ignoring the magnetic dipole and the electric quadrupole contributions. This is adequate for molecular complexes of nonchiral molecules. The response function R may be extended to include the magnetic corrections. However, it is much more convenient to calculate R instead.

We have used the Frenkel-exciton model in the boson representation. Two types of parameters are responsible for the nonlinear optical signals. The quartic potential Umnkl controls the exciton scattering when the field is turned off, while μm(2) represents anharmonic system–field interaction. Hard-core bosons are a special case of this model, as shown in Appendix C. They have Ummmm = ∞. This model developed here for electronic excitations also applies to vibrational excitations. We then model the excitations as intrinsically boson-like with a finite Ummmm and μm(2)=0 (soft-core bosons). This approach has been applied to coupled vibrations in polypeptides.67,100,337

The supermolecule and quasiparticle scattering are two conceptually different approaches for calculating the response functions of aggregates. The sources of nonlinearities in these descriptions look very different. In the supermolecular representation, the double-exciton states are included explicitly, and interband transitions do not follow a harmonic oscillator model. In the quasiparticle representation, the final expressions involve only single-exciton states. Double-exciton resonances are recovered from the frequency-dependent exciton-scattering matrix. The relation and equivalence of the two representations can be established when supermolecule double-exciton states are expanded in the single-exciton product basis for an isolated system. Differences remain since different approximations are made for the system–bath interaction.

In this review, the quasiparticle representation has been used in conjunction with simple Lorentzian line shapes in the Markovian limit. More elaborate line-shape theories could be developed (equivalent to the CGF or SLE); however, this is usually not necessary in large aggregates where the lineshapes are often controlled by inhomogeneities and spectral congestions. Each peak is then composed of many overlapping transitions. The detailed microscopic description of bath dynamics is crucial for small systems (~10 or less coupled oscillators), where the line shapes are better resolved.

The marriage of coherent-control pulse-shaping techniques with 2D spectroscopy offers numerous future possibilities. The rich information obtained through snapshots of electronic dynamics can then be further refined by manipulating the interference of various Liouville space pathways. These techniques may also be used to shift and isolate specific resonances, facilitating the assignment of electronic transitions.

Chirality-induced signals may lead to useful applications of multidimensional coherent spectroscopy to structure determination. X-ray crystallography and NMR determine structures of molecular complexes with atomic resolution. Optical spectroscopy can be used as a complimentary technique for further structure refinement.

Electronic nonlinear spectroscopy is most valuable for tracking dynamical events in molecules and their complexes. Unlike traditional 1D transient-absorption and transient-grating techniques, in multidimensional coherent techniques, a large amount of information is gathered in a single shot. Numerical data processing and analysis techniques need to be developed to help extract system information (e.g., Hamiltonian, bath, structural, and relaxation parameters) from the measured signals.

The various methods and techniques reviewed here are based on a generic exciton and supermolecule approaches within the RWA (sections 2, 5, 6) and may be applied to other frequency regimes: from THz to X-ray. 2D techniques, originally developed in NMR,53 were extended during the past decade to IR laser spectroscopy.60,338343 Simulation protocols based on molecular dynamics (MD) strategies have been reviewed recently in ref 67. These applications use the same formalism; however, the vibrations are bosons. Electrostatic parametrization of vibrational Hamiltonian parameters is currently a useful technique for including fluctuations, simulated by MD trajectories. Developing electrostatic maps for electronic Hamiltonians is necessary for affordable modeling of solvation effects.

The theoretical apparatus presented in this review has been implemented in the computer software package “SPECTRON”, which has been used for simulating coherent nonlinear optical responses in IR,79,91,92,94,95,155,344,345 visible,81,97 and UV82,156 regions. The code includes the super-molecule and quasiparticle approaches as well as the spectral lineshape-broadening models described in sections 2, 4, 9, and 11.

Figure 33
Feynman diagrams of the quasiparticle dynamics in kII technique. The left-most diagram is for the μ(0) = 0 case, while others include μ(0) ≠ 0. The symbols are the same as in Figure 32.


This research was supported by the National Science Foundation (CHE-0745892) and National Institutes of Health (GM59230), and the Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U. S. Department of Energy. This support is gratefully acknowledged. F. Š. acknowledges support of GAČR (Grant No. 202/07/P245) and MŠMT ČR (MSM 0021620835). We wish to thank Dr. W. Zhuang and Dr. O. Roslyak for many useful discussions.


An external file that holds a picture, illustration, etc.
Object name is nihms212366b1.gif

Darius Abramavicius was born in 1974 in Alytus, Lithuania. He received his M.S. degree from Vilnius University in 1998 and his Ph.D. degree under the supervision of Prof. Leonas Valkunas from the Institute of Physics at Vytautas Magnus University, Lithuania, in 2002. He then joined the group of Professor Shaul Mukamel at University of Rochester (New York, U.S.A.). He moved to the University of California Irvine in 2003. During 2005–2006, he was a temporary lecturer in Theoretical Physics Department of Vilnius University in Lithuania. Currently he is a visiting scientist at the University of California, Irvine, in the group of Prof. S. Mukamel. His research interests are in the field of optical and dynamical properties of molecular complexes.

An external file that holds a picture, illustration, etc.
Object name is nihms212366b2.gif

Benoit Palmieri was born in 1978 in Montréal, Canada. He received his Ph.D. degree under the direction of Professor David Ronis from the chemistry department of McGill University in 2007. The main subject of his Ph.D. thesis was the theoretical study of diffusion of gas or liquid molecules inside channeled structures (zeolites). He joined the group of Professor Shaul Mukamel in 2007, where he worked on the nonlinear optical response of photosynthetic systems and on the importance of interband excitonic couplings. He is now working at the Canadian Space Agency under the direction of Dr. Marcus Dejmek, developing theoretical models of crystal growth.

An external file that holds a picture, illustration, etc.
Object name is nihms212366b3.gif

Dmitri V. Voronine was born in 1979 in Moscow, Russia, and is currently a Postdoctoral Researcher in the group of Prof. Tobias Brixner at the Institut für Physikalische Chemie, Universität Würzburg, Germany. He received his Ph.D. in 2004 from the Center for Photochemical Sciences, Bowling Green State University, under the supervision of Prof. Michael A. J. Rodgers, and did postdoctoral work with Prof. Shaul Mukamel at the Department of Chemistry, University of California, Irvine. He has been the recipient of the McMaster Fellowship (BGSU) and the RCCM Excellent Young Investigator award (U. Würzburg). His research interests in theoretical and experimental physical chemistry and nanobiophysics include applying the methods of coherent control with laser pulse shaping to femtosecond 2D spectroscopy of complex systems such as biomolecules, photosynthetic complexes, chiral and fractal molecular aggregates, porphyrins and carotenoids; and ultrafast nano-optics of metal nanostructures for space-time-resolved spectroscopy with nanometer spatial and femtosecond temporal resolution.

An external file that holds a picture, illustration, etc.
Object name is nihms212366b4.gif

František Šanda was born in 1976 in Domažlice, Czech Republic. He received his Ph.D. in Theoretical Physics (supervisor Vladislav Čápek) in 2003 from Charles University, Prague. From 2004 to 2005, he had a postdoctoral appointment at the University of California, Irvine, in the group of Prof. Shaul Mukamel. He is an Assistant Professor at Institute of Physics, Faculty of Mathematics and Physics, Charles University in Prague since 2006. His research interest focuses on stochastic modeling of dynamic processes in Liouville space, with emphasis on applications in spectroscopy.

An external file that holds a picture, illustration, etc.
Object name is nihms212366b5.gif

Shaul Mukamel, currently the Chancellor Professor of Chemistry at the University of California, Irvine, received his Ph.D. in 1976 from Tel Aviv University. Following postdoctoral appointments at MIT and the University of California, Berkeley, he has held faculty positions at Rice University, the Weizmann Institute, and the University of Rochester. He is the recipient of the Sloan, Dreyfus, Guggenheim, and Alexander von Humboldt Senior Scientist award and the Lippincort Award and is a fellow of the American Physical Society and the Optical Society of America. His interests focus on developing computational techniques for the design of novel ultrafast laser pulse sequences for probing electronic and vibrational dynamics in molecules. Biophysical applications include folding and dynamical fluctuations in proteins and DNA, hydrogen bonding, long-range electron and energy transfer in photosynthetic complexes, and signatures of chirality. Other areas are attosecond x-ray spectroscopy, excitons in semiconductor nanostructures, many-body effects in quantum optics, and photon and electron statistics in single-molecule spectroscopy and open molecular junctions. He is the author of over 600 publications in scientific journals and the textbook, Principles of Nonlinear Optical Spectroscopy (Oxford University Press), 1995.

17. Appendix A: Survey of Linear Optical Techniques

Linear optical signals are generated by a single interaction with the incident optical field. Three types of linear signals are most widely used: absorption, where the transmission of the incident beam is measured; linear dichroism, the difference between two polarized absorptions (longitudinal, L, and transverse, T) in oriented samples; and circular dichroism, the difference between the absorptions of two circularly polarized beams (left, L, and right, R).

Using the response function formalism, the linear absorption may be interpreted as follows. The incident beam interacts once with the molecule and induces the linear polarization. This polarization then creates the outgoing optical field, which interferes with the incoming beam. This is a self-heterodyne detection. The detector measures the net optical field, which is attenuated by the interference. The mathematical formulation is very nicely described in the book by Loudon.117

The absorption of the field is defined by following the field intensity in the medium, which follows


The absorption coefficient κa is given by


Here, n(ω) is the refractive index of the medium, and c is the speed of light. The macroscopic response function of an isotropic ensemble of identical absorbers is related to the microscopic single-molecule response function by a three-dimensional orientational averaging. Using the response function eq 13, we get for the absorption coefficient


Assemblies of oriented molecules are described by their tensor properties. The response function eq 13 is a second-rank tensor with respect to system symmetry axis. Equation 13 then involves the following orientationally averaged product:


Here, we introduced the transformation tensor from the molecular (μα) to the laboratory frame (μi), μi = ΣαTiαμα We assume that the molecules are oriented along z and take the molecular z-axis to coincide with the laboratory z-axis. Orientational averaging in the xy plane then amounts to rotating the molecular xy plane around the z-axis. This can be easily calculated since, in this case, only five elements of the transformation tensor are nonzero: Txx = cos [var phi], Txy = −sin [var phi], Tyx = -sin [var phi], Tyy = cos [var phi], and Tzz = 1, and the orientational averaging amounts to an integration over the rotation angle [var phi]:


Thus, when μ2 = μ1 =μ, we obtain the following possible field configurations. For ν2 = ν1= z:


and for ν2 = ν1 = x (the same holds for y):


The linear dichroism (LD) technique measures the difference between the absorption of two perpendicularly polarized beams in an oriented sample. The LD signal vanishes in unoriented samples. For an oriented sample, it is given by


where κL is the longitudinal absorption given by eq 269 with the transition amplitude eq 272:


whereas κT is the transverse absorption (eq 269) together with eq 273:


Circular dichroism (CD) is a chirality-induced linear technique, involving two circularly polarized beams. A circularly polarized beam propagating along z is described by rotating unit electric vector:


where “− (+)” sign corresponds to R(L) circular polarization.

Using the response function, we get the induced polarization generated by the R/L beam:


The circular dichroism signal is commonly defined in terms of the ellipticity angle,


where ER/L are the transmitted field amplitudes of the R and L beams given by eq 282. For very small angles, we use tan(θ) ≈ θ and get


Instead of ellipticity angle, the molar ellipticity angle and the rotational strength are often used.93

18. Appendix B: Homodyne-Detected versus Heterodyne-Detected Coherent Signals

The induced polarization serves as a source for the signal field, which is eventually detected. Several detection schemes are possible. Neglecting the magnetic response, the electric field and polarization satisfy the Maxwell equation:


For transverse field and for pencil-shape sample geometry, it is possible to derive a simple relation between the emitted field envelope and the envelope of the induced polarization:


where n(ωs) is the dielectric constant at frequency ωs, l is the interaction length, and Δk = ksu1k1u2k2u3k3 is the wavevector mismatch between the polarization and the optical field; sinc(x) = sin(x)/x. Thus, the strongest amplitude occurs at the phase-matching condition Δk =0. The other important conclusion is that the signal field is directly proportional to the polarization amplitude, but is π/2 out of phase because of the i factor.

Two detection schemes for the signal are most commonly employed. These are known as homodyne and heterodyne detection. In homodyne detection (HOD), the detector is placed at the phase-matching direction and measures the integrated intensity of the signal field, or counts photons. The field intensity is proportional to the absolute value of the signal field. Thus, the homodyne signal is given by


An additional time-gating device can be introduced that filters the incoming signal in a given time window before it is detected by the homodyne technique. This time-gated ho-modyne (TGH) detected signal is given by


In heterodyne detection (HET), an additional laser pulse (local oscillator, EH) is applied in the direction of the signal field, ES. It interferes with the signal, and the resulting field is then detected by homodyne detector. We thus get


The local oscillator (second term) can be readily subtracted. The |Es(t)|2 term is typically weak since Es [double less-than sign] EH, and thus, the HET approach detects the signal amplitude


This measurement can be repeated with varying phase of EH to reveal the induced polarization Ps(t) using eq 282.

19. Appendix C: Bosonization Schemes for Frenkel Excitons

Transformations between fermion and boson representations widely used in many-body problems may be employed for describing excitations in semiconductors346353 and molecular aggregates.92,136,137,354,355

Consider molecular aggregates made of multilevel molecules. Molecule m has a set of transitions a from its ground state. Each transition acts as a two-level system, where the molecule can be either in the excited state a or the ground state. The molecule is now viewed as a collection of two-level systems. Within the space of a single excitation and double excitations, an aggregate made of such molecules is described by the Frenkel-exciton Hamiltonian


where m, n label different molecules and b^mb is a bosonic ( [b^ma,b^nb]=δmnδab) exciton creation operator on molecule m of transition a (this can be brought into the form of eq 57 by relabeling pairs ma, nb, kc and ld as single indices). Vma,nb is the exciton resonant coupling term, and Wma,nbkc,ld is the two-boson interaction potential. When W [double less-than sign] V (weak scattering regime), the excitations are soft-core bosons, where two particles have a small penalty (anharmonicity), when occupying the same state in the molecule. The requirement that two excitations cannot reside on the same molecule is satisfied by shifting the unphysical state energies out of the physical window of interest. In the Hamiltonian equation (eq 287), this is done by setting Wma,mbmc,md. The electronic excitations now become hard-core bosons. Note that, by varying Wma,mbmc,md, we can continuously turn soft-core into hard-core bosons.

As a simple example, consider two-level molecules. The bosons do not scatter when they are on different molecules. To describe this model, we set Wm,mm,m for all m (the indices a, b, c, and d are redundant for this model), while it is finite for all other index combinations and describes double-exciton binding energies. The exciton-scattering matrix for coupled two-level systems can be calculated by taking this limit in eq 325,


where all other elements vanish. Note that it involves the inversion of N × N matrix Gmm,nn<