Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2975548

Formats

Article sections

- 1. Introduction
- 2. Supermolecule Approach; Coherent Optical Response of Multilevel Systems
- 3. Response Functions of a Molecular Aggregate
- 4. Optical Response of Excitons with Transport; The Lindblad Master Equation
- 5. Multidimensional Spectroscopy of a Three-Band Model; Heterodyne-Detected Signals
- 6. Quasiparticle Representation of the Optical Response; The Nonlinear-Exciton Equations (NEE)
- 7. Coherent Nonlinear Optical Response of Quasiparticles in the Molecular Basis
- 8. Coupling with Phonons; Exciton-Dephasing and Transport
- 9. Optical Response of Quasiparticles with Relaxation
- 10. Applications to the FMO Complex
- 11. Interplay of Transport and Slow Fluctuations in the Third-Order Response
- 12. Explicit Treatment of Bath Dynamics; The Stochastic Liouville Equations and Beyond
- 13. Nonlocal Nonlinear Response of Chiral Aggregates
- 14. Manipulating 2D Signals by Coherent-Control Pulse-Shaping Algorithms
- 15. Discussion and Conclusions
- 29. References

Authors

Related links

Chem Rev. Author manuscript; available in PMC 2010 November 8.

Published in final edited form as:

PMCID: PMC2975548

NIHMSID: NIHMS212366

Darius Abramavicius,^{†} Benoit Palmieri,^{†} Dmitri V. Voronine,^{‡} František Šanda,^{§} and Shaul Mukamel^{*}^{†}

University of California Irvine, California 92623; Institut für Physikalische Chemie, Universität Würzburg, Germany; and Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic

Molecular aggregates are abundant in nature; they form spontaneously in concentrated solutions and on surfaces and can be synthesized by supramolecular chemistry techniques.^{1}^{–}^{3} 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,^{6}^{–}^{8} circular complexes,^{9} cylindrical nanotubes, multiwall cylinders and supercylinders,^{10}^{–}^{12} 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}

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.^{17}^{–}^{20} Applications of molecular exciton theory to dimers^{4}^{,}^{21}^{–}^{24} 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}^{,}^{25}^{–}^{29} 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.^{31}^{–}^{36} Large molecular aggregates may show several optical absorption bands^{12} 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}^{,}^{37}^{–}^{41} 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}^{,}^{42}^{–}^{52}

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.^{54}^{–}^{56} The connection with NMR was established.^{57}^{–}^{59} Multidimensional infrared spectroscopy monitors hydrogen-bonding network and dynamics, protein structures, and their fluctuations.^{60}^{–}^{67} Development of these IR techniques was reviewed recently^{67} 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 ^{68}^{–}^{70} and demonstrated experimentally in molecules^{40}^{,}^{41}^{,}^{65}^{,}^{71}^{,}^{72} and semiconductor quantum wells.^{73}^{–}^{75} 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, *t*_{1}, *t*_{2}, and *t*_{3}, between them. In practice the *t*_{3} information is usually obtained interferometrically in one shot rather than by scanning *t*_{3}. Spectral dispersion of the signal with the local oscillator gives the Fourier transform with respect to *t*_{3}. 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,^{57}^{–}^{59} 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 (*k** _{s}*), 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.

Coherent third-order nonlinear optical experiment. The four laser pulses are ordered in time; the signal is generated in the phase-matching direction *k*_{4}. Data processing of time-domain signals and their parametric dependence on the three delays *t*_{1}, *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}^{,}^{88}^{–}^{92} 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.^{95}^{–}^{97}

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 crystals^{17}^{,}^{98} and subsequently extended to aggregates.^{4}^{,}^{5}^{,}^{27}^{–}^{29}^{,}^{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 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 *t*_{2} 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}^{,}^{101}^{–}^{107} 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.^{108}^{–}^{111} 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}^{,}^{114}^{–}^{116} 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 of sites; in contrast, the holes and electrons of Wannier excitons are loosely bound, and the number of single-exciton variables scales as ~^{2}. This unfavorable scaling is more severe for double-exciton states (~^{2} for Frenkel and ~^{4} 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.

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}

$$\widehat{H}={\widehat{H}}_{0}+{\widehat{H}}^{\prime}(t),$$

(1)

where

$${\widehat{H}}_{0}=\sum _{a}\hslash {\epsilon}_{a}\mid a\rangle \langle a\mid .$$

(2)

Here, *ħε _{a}* is the energy of state

$${\widehat{H}}^{\prime}(t)=-\widehat{\mathit{\mu}}\xb7\mathit{E}(\mathit{r},t),$$

(3)

where ** r** is the position of the molecule and

$$\widehat{\mathit{\mu}}=\sum _{ab}{\mathit{\mu}}_{ab}\mid a\rangle \langle b\mid $$

(4)

is the dipole operator. *μ** _{ab}* = Σ

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

$$\widehat{\rho}=\sum _{ab}{\rho}_{ab}\mid a\rangle \langle b\mid \phantom{\rule{0.16667em}{0ex}}.$$

(5)

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,

$$i\hslash \frac{\text{d}\widehat{\rho}}{\text{d}t}=\mathbb{L}\widehat{\rho}.$$

(6)

Here , the Liouville superoperator, is given by the commutator, with the Hamiltonian = [*Ĥ*, ]. Similar to eq 1, we partition the Liouville operator as

$$\mathbb{L}={\mathbb{L}}_{0}+{\mathbb{L}}^{\prime}(t),$$

(7)

where represents the isolated system and = [*Ĥ′*, ] 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 = 0, we get

$${\widehat{\rho}}^{(0)}(t)=\mathbb{G}(t)\widehat{\rho}(0)=\theta (t)exp\left(-\frac{i}{\hslash}{\widehat{H}}_{0}t\right)\widehat{\rho}(0)exp\left(+\frac{i}{\hslash}{\widehat{H}}_{0}t\right),$$

(8)

where θ(*t*) is the Heavyside step-function, defined by
$\theta (t)={\int}_{-\infty}^{t}\text{d}\tau \delta (\tau )$. The density matrix of the driven system (eq 6) can be calculated perturbatively in by iterative integration of eqs 6 and 7. This yields^{52}

$$\widehat{\rho}(t)=\mathbb{G}(t)\widehat{\rho}(0)+\sum _{n=1}^{\infty}{\left(-\frac{i}{\hslash}\right)}^{n}{\int}_{0}^{t}\text{d}{\tau}_{n}{\int}_{0}^{{\tau}_{n}}\text{d}{\tau}_{n-1}\cdots {\int}_{0}^{{\tau}_{2}}\text{d}{\tau}_{1}\times \mathbb{G}(t-{\tau}_{n}){\mathbb{L}}^{\prime}({\tau}_{n})\mathbb{G}({\tau}_{n}-{\tau}_{n-1}){\mathbb{L}}^{\prime}({\tau}_{n-1})\cdots \times \mathbb{G}({\tau}_{2}-{\tau}_{1}){\mathbb{L}}^{\prime}({\tau}_{1})\mathbb{G}({\tau}_{1})\widehat{\rho}(0).$$

(9)

Equation 9 provides the order-by-order expansion of the density matrix in the field and can be recast as (*t*) = ^{(0)}(*t*) +^{(1)}(*t*) + ^{(2)}(*t*) +... . ^{(}^{n}^{)}(*t*) is the density matrix to *n*th 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

$${\mathit{P}}^{(n)}=Tr[\widehat{\mathit{\mu}}{\widehat{\rho}}^{(n)}].$$

(10)

The induced polarization is given by ** P**(

The *linear response function* connects the linear polarization with the field:^{52}

$${\mathit{P}}^{(1)}(\mathit{r},t)={\int}_{0}^{\infty}\text{d}{t}_{1}{R}^{(1)}({t}_{1})\mathit{E}(\mathit{r},t-{t}_{1}).$$

(11)

Since *P*^{(1)} and ** E** are both vectors,

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

$${R}^{(1)}({t}_{1})=\left(\frac{i}{\hslash}\right)Tr[\widehat{\mathit{\mu}}\mathbb{G}({t}_{1})\mathbb{V}{\widehat{\rho}}_{0}],$$

(12)

where _{0} is the equilibrium density matrix and [**, ] is the dipole superoperator.**

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

Excitonic aggregate made of coupled three-level chromophores. *ε*_{n} is the excitation energy of each chromophore and 2(*ε*_{n} + *U*_{nnnn}) is its double-excitation energy. The chromophores are quadratically coupled by *J*_{mn} (and quartically by *U*_{mnkl} **...**

Feynman diagrams for the signal generated in the direction *k*_{I} = −*k*_{1} + *k*_{2} + *k*_{3}; 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,

$${R}^{(1)}(t)=\left(\frac{i}{\hslash}\right)\theta (t)\sum _{e}{\mathit{\mu}}_{ge}{\mathit{\mu}}_{eg}{\text{e}}^{-i{\omega}_{eg}t-{\gamma}_{eg}t}+c.c.,$$

(13)

where the equilibrium density matrix is _{0} *= |g**g*|, *ω _{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 by^{52}

$${\mathit{P}}^{(3)}(\mathit{r},t)={\int}_{0}^{\infty}\text{d}{t}_{3}{\int}_{0}^{\infty}\text{d}{t}_{2}{\int}_{0}^{\infty}\text{d}{t}_{1}{R}^{(3)}({t}_{3},{t}_{2},{t}_{1})\times \mathit{E}(\mathit{r},t-{t}_{3})\mathit{E}(\mathit{r},t-{t}_{3}-{t}_{2})\mathit{E}(\mathit{r},t-{t}_{3}-{t}_{2}-{t}_{1}),$$

(14)

where

$${R}^{(3)}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}Tr[\widehat{\mathit{\mu}}\mathbb{G}({t}_{3})\mathbb{V}\mathbb{G}({t}_{2})\mathbb{V}\mathbb{G}({t}_{1})\mathbb{V}{\rho}_{0}]$$

(15)

is the third-order response function.

This expression represents a sequence of three interactions with the incoming fields described by , three free propagation periods between interactions (), and signal generation described by the last **. Since is a commutator and **** can act either from the left or the right, eq 15 has 2**^{3} = 8 contributions known as *Liouville space pathways*.

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

$$\mathit{E}(t)=\sum _{j}\sum _{{u}_{j}=\pm 1}{\mathbb{E}}_{j}^{{u}_{j}}(t-{\tau}_{j})exp[{iu}_{j}({\mathit{k}}_{j}\mathit{r}-{\omega}_{j}(t-{\tau}_{j}))].$$

(16)

Here, (*t* − *τ _{j}*) is the complex envelope of the

The field is generated in the directions ±*k*_{1} ±*k*_{2} ±*k*_{3}. The detector selects a signal in one direction. To calculate the optical signals, we introduce the wavevector-dependent polarization *P*_{ks}

$${\mathit{P}}^{(3)}(\mathit{r},t)=\sum _{{\mathit{k}}_{s}}{\mathit{P}}_{{\mathit{k}}_{s}}(t)exp(i{\mathit{k}}_{s}\mathit{r}),$$

(17)

where *k** _{s}* =

$${\mathit{P}}_{{\mathit{k}}_{s}}(t)=exp[-i{\omega}_{s}(t-{\tau}_{3})-i({u}_{2}{\omega}_{2}+{u}_{1}{\omega}_{1})({\tau}_{3}-{\tau}_{2})-{iu}_{1}{\omega}_{1}({\tau}_{2}-{\tau}_{1})]\times \int \int {\int}_{0}^{\infty}\text{d}{t}_{3}\text{d}{t}_{2}\text{d}{t}_{1}{R}_{{\mathit{k}}_{s}}^{(3)}({t}_{3},{t}_{2},{t}_{1})\times exp[i{\omega}_{s}{t}_{3}+i({u}_{2}{\omega}_{2}+{u}_{1}{\omega}_{1}){t}_{2}+{iu}_{1}{\omega}_{1}{t}_{1}]{\mathbb{E}}_{3}^{{u}_{3}}(t-{t}_{3}-{\tau}_{3})\times {\mathbb{E}}_{2}^{{u}_{2}}(t-{t}_{3}-{t}_{2}-{\tau}_{2}){\mathbb{E}}_{1}^{{u}_{1}}(t-{t}_{3}-{t}_{2}-{t}_{1}-{\tau}_{1}),$$

(18)

with the signal frequency given by *ω _{s}* =

There are four independent third-order techniques, defined by the various choices of *u _{j}*:

$${\mathit{k}}_{\text{I}}=-{\mathit{k}}_{1}+{\mathit{k}}_{2}+{\mathit{k}}_{3}\phantom{\rule{0.38889em}{0ex}}({u}_{1},{u}_{2},{u}_{3})=(-1,1,1),$$

(19)

$${\mathit{k}}_{\text{II}}=+{\mathit{k}}_{1}-{\mathit{k}}_{2}+{\mathit{k}}_{3}\phantom{\rule{0.38889em}{0ex}}({u}_{1},{u}_{2},{u}_{3})=(1,-1,1),$$

(20)

$${\mathit{k}}_{\text{III}}=+{\mathit{k}}_{1}+{\mathit{k}}_{2}-{\mathit{k}}_{3}\phantom{\rule{0.38889em}{0ex}}({u}_{1},{u}_{2},{u}_{3})=(1,1,-1),$$

(21)

$${\mathit{k}}_{\text{IV}}=+{\mathit{k}}_{1}+{\mathit{k}}_{2}+{\mathit{k}}_{3}\phantom{\rule{0.38889em}{0ex}}({u}_{1},{u}_{2},{u}_{3})=(1,1,1).$$

(22)

For our excitonic level scheme and its dipole selection rules, the *k*_{IV} 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 4–6. 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 exp[−

The *k*_{I} signal is depicted in Figure 4. During *t*_{1}, 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 (

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 *t*_{2}; thus, the third interaction reflects the reduction of the ground-state population, which reduces (bleaches) the subsequent absorption. In the ESE pathway, during *t*_{2}, 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 *t*_{1} and *t*_{2} dynamics as the ESE; however, the third interaction creates a doubly excited state *f*.

During the intervals between interactions (*t*_{1}, *t*_{2}, and *t*_{3}), 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 *t*_{1} interval is in the state *ρ _{ge}*, and its evolution gives an exp(−

$${R}_{i}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{ee}^{\prime}}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{eg}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{ge}\times exp(-i{\omega}_{ge}{t}_{1}-i{\omega}_{{e}^{\prime}e}{t}_{2}-i{\omega}_{{e}^{\prime}g}{t}_{3}-{\gamma}_{ge}{t}_{1}-{\gamma}_{{e}^{\prime}e}{t}_{2}-{\gamma}_{{e}^{\prime}g}{t}_{3}),$$

(23)

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

Similar to the three contributions to the *k*_{I} technique shown in Figure 4, *k*_{II} has three contributions (Figure 5) and *k*_{III} 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 4–6 may be obtained from eqs 32–39 by setting the population relaxation matrix to unity,
${\mathbb{G}}_{{e}_{2}{e}_{2},{e}_{1}{e}_{1}}^{(N)}({t}_{2})=\theta ({t}_{2}){\delta}_{{e}_{2}{e}_{1}}$.

The time evolution of closed quantum systems can be expanded in eigenstates that evolve independently by simply acquiring exp(−*iε _{j}t*) 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

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ödinger*–*Langevin equation*,

$$\stackrel{.}{\psi}=-\frac{i}{\hslash}\widehat{H}\psi -\frac{1}{2}\sum _{\alpha}{\widehat{V}}_{\alpha}^{\u2020}{\widehat{V}}_{\alpha}\psi +\sum _{\alpha}{\xi}_{\alpha}(t){\widehat{V}}_{\alpha}\psi ,$$

(24)

where *ψ* is the wave function of the open system, is an arbitrary set of system operators (independent of *ξ*) that describe its coupling with the bath, and *ξ _{α}*(

The QME that corresponds to this Schrödinger–Langevin equation is known as the *Lindblad equation*.^{120}^{–}^{122}

$$\stackrel{.}{\widehat{\rho}}=-\frac{i}{\hslash}[\widehat{H},\widehat{\rho}]+\sum _{\alpha}\left({\widehat{V}}_{\alpha}\widehat{\rho}{\widehat{V}}_{\alpha}^{\u2020}-\frac{1}{2}{\widehat{V}}_{\alpha}^{\u2020}{\widehat{V}}_{\alpha}\widehat{\rho}-\frac{1}{2}\widehat{\rho}{\widehat{V}}_{\alpha}^{\u2020}{\widehat{V}}_{\alpha}\right).$$

(25)

It is derived from the Liouville equation for (*t*) | *ψ*(*t*) (*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*() 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 _{α}. The secular Redfield equations decouple the populations (diagonal elements of in the eigenstate basis) from the coherences (off-diagonal elements). The populations then satisfy the Pauli master equation:^{109}^{,}^{111}

$$\frac{\text{d}}{\text{d}t}{\rho}_{ee}(t)=-\sum _{{e}^{\prime}}{K}_{ee,{e}^{\prime}{e}^{\prime}}{\rho}_{{e}^{\prime}{e}^{\prime}}(t).$$

(26)

Here, *K _{ee}*

The evolution of diagonal density-matrix elements is described by the population Green’s function,
${\rho}_{{e}^{\prime}{e}^{\prime}}(t)={\sum}_{e}{\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}(t){\rho}_{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

$${\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}(t)={[exp(-Kt)]}_{{e}^{\prime}{e}^{\prime},ee}\equiv \sum _{p}{\mathit{\chi}}_{{e}^{\prime}p}^{(\text{R})}{D}_{pp}^{-1}exp(-{\lambda}_{p}t){\mathit{\chi}}_{pe}^{(\text{L})},$$

(27)

where *λ _{p}* is the

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

$$\frac{\text{d}}{\text{d}t}{\rho}_{{e}^{\prime}e}(t)=[-i{\omega}_{{e}^{\prime}e}-{\gamma}_{{e}^{\prime}e}^{(N)}]{\rho}_{{e}^{\prime}e}(t),\phantom{\rule{0.38889em}{0ex}}e\ne {e}^{\prime}.$$

(28)

Its solution is ${\rho}_{{e}^{\prime}e}(t)=exp(-i{\omega}_{{e}^{\prime}e}t-{\gamma}_{{e}^{\prime}e}^{(N)}t)$, where

$${\gamma}_{{e}^{\prime}e}^{(N)}=\frac{1}{2}({K}_{ee,ee}+{K}_{{e}^{\prime}{e}^{\prime},{e}^{\prime}{e}^{\prime}})+{\widehat{\gamma}}_{e,{e}^{\prime}},$$

(29)

The first two terms represent population relaxation, and _{e}_{,}* _{e′}* is the

$${\rho}_{{e}_{4}{e}_{3}}(t)=\sum _{{e}_{2}{e}_{1}}{\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}(t){\rho}_{{e}_{2}{e}_{1}}(0).$$

(30)

Comparison of eq 30 with eqs 27 and 28 gives

$${\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}(t)={\delta}_{{e}_{4}{e}_{3}}{\delta}_{{e}_{2}{e}_{1}}\theta (t){[exp(-Kt)]}_{{e}_{4}{e}_{4},{e}_{2}{e}_{2}}+(1-{\delta}_{{e}_{4}{e}_{3}}){\delta}_{{e}_{4}{e}_{2}}{\delta}_{{e}_{3}{e}_{1}}\theta (t)exp[-i{\omega}_{{e}_{4}{e}_{3}}t-{\gamma}_{{e}_{4}{e}_{3}}^{(N)}t].$$

(31)

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

For the *k*_{I} and *k*_{II} techniques, populations are generated during the second time interval *t*_{2} in the ESA and ESE diagrams of Figures 4 and and55 when *e* = *e′* (GSB includes ground-state population during *t*_{2}; 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 *t*_{2} 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 *k*_{I} technique (Figures 4 and and77),

$${R}_{{\mathit{k}}_{\text{I}},\text{i}}^{(\text{SOS})}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}{\mathit{\mu}}_{g{e}_{4}}{\mathit{\mu}}_{{e}_{3}g}{\mathit{\mu}}_{{e}_{2}g}{\mathit{\mu}}_{g{e}_{1}}\times exp(-i{\xi}_{{e}_{4}g}{t}_{3}-i{\xi}_{g{e}_{1}}{t}_{1}){\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}({t}_{2}),$$

(32)

$${R}_{{\mathit{k}}_{\text{I}},\text{ii}}^{(\text{SOS})}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{e}_{4},{e}_{1}}{\mathit{\mu}}_{g{e}_{4}}{\mathit{\mu}}_{{e}_{4}g}{\mathit{\mu}}_{{e}_{1}g}{\mathit{\mu}}_{g{e}_{1}}\times exp(-i{\xi}_{{e}_{4}g}{t}_{3}-i{\xi}_{g{e}_{1}}{t}_{1}),$$

(33)

$${R}_{{\mathit{k}}_{\text{I}},\text{iii}}^{(\text{SOS})}({t}_{3},{t}_{2},{t}_{1})=-{\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{f{e}_{4}{e}_{3}{e}_{2}{e}_{1}}{\mathit{\mu}}_{{e}_{3}f}{\mathit{\mu}}_{f{e}_{4}}{\mathit{\mu}}_{{e}_{2}g}{\mathit{\mu}}_{g{e}_{1}}\times exp(-i{\xi}_{f{e}_{3}}{t}_{3}-i{\xi}_{g{e}_{1}}{t}_{1}){\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}({t}_{2}),$$

(34)

where we have introduced the complex frequencies *ξ _{ab}*

The oscillation frequencies during *t*_{1} and *t*_{3} have opposite sense. For a two-level system (*g* and *e*), the frequency factors cancel out at *t*_{1} = *t*_{3}. 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 *t*_{2} when transport is neglected.

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

For the *k*_{II} technique, we similarly have (Figures 5 and and77)

$${R}_{{\mathit{k}}_{\text{II}},\text{iv}}^{(\text{SOS})}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}{\mathit{\mu}}_{g{e}_{4}}{\mathit{\mu}}_{{e}_{3}g}{\mathit{\mu}}_{g{e}_{1}}{\mathit{\mu}}_{{e}_{2}g}\times exp(-i{\xi}_{{e}_{4}g}{t}_{3}-i{\xi}_{{e}_{2}g}{t}_{1}){\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}({t}_{2}),$$

(35)

$${R}_{{\mathit{k}}_{\text{II}},\text{v}}^{(\text{SOS})}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{e}_{4}{e}_{2}}{\mathit{\mu}}_{g{e}_{4}}{\mathit{\mu}}_{{e}_{4}g}{\mathit{\mu}}_{g{e}_{2}}{\mathit{\mu}}_{{e}_{2}g}\times exp(-i{\xi}_{{e}_{4}g}{t}_{3}-i{\xi}_{{e}_{2}g}{t}_{1}),$$

(36)

$${R}_{{\mathit{k}}_{\text{II}},\text{vi}}^{(\text{SOS})}({t}_{3},{t}_{2},{t}_{1})=-{\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{f{e}_{4}{e}_{3}{e}_{2}{e}_{1}}{\mathit{\mu}}_{{e}_{3}f}{\mathit{\mu}}_{f{e}_{4}}{\mathit{\mu}}_{g{e}_{1}}{\mathit{\mu}}_{{e}_{2}g}\times exp(-i{\xi}_{f{e}_{3}}{t}_{3}-i{\xi}_{{e}_{2}g}{t}_{1}){\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}({t}_{2}).$$

(37)

The two contributions to the *k*_{III} (Figure 6) technique are of the ESA type:

$${R}_{{\mathit{k}}_{\text{III}},\text{vii}}^{(\text{SOS})}({t}_{3},{t}_{2},{t}_{1})=-{\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{ee}^{\prime}f}{\mathit{\mu}}_{{e}^{\prime}f}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{fe}{\mathit{\mu}}_{eg}\times exp(-i{\xi}_{f{e}^{\prime}}{t}_{3}-i{\xi}_{fg}{t}_{2}-i{\xi}_{eg}{t}_{1}),$$

(38)

$${R}_{{\mathit{k}}_{\text{III}},\text{viii}}^{(\text{SOS})}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{ee}^{\prime}f}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{{e}^{\prime}f}{\mathit{\mu}}_{fe}{\mathit{\mu}}_{eg}\times exp(-i{\xi}_{{e}^{\prime}g}{t}_{3}-i{\xi}_{fg}{t}_{2}-i{\xi}_{eg}{t}_{1}).$$

(39)

Note that exciton transport only enters *k*_{I} and *k*_{II} but not *k*_{III} 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 aggregates^{41}^{,}^{48}^{,}^{126} require higher levels of theory that will be presented in the coming sections.

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

$${S}_{{\mathit{k}}_{s}}^{(3)}({T}_{3},{T}_{2},{T}_{1})={\int}_{-\infty}^{\infty}\text{d}t\phantom{\rule{0.16667em}{0ex}}{\mathit{P}}_{{\mathit{k}}_{s}}(t){\mathbb{E}}_{s}^{\ast}(t-{\tau}_{s}){\text{e}}^{i{\omega}_{s}(t-{\tau}_{s})},$$

(40)

where the signal depends parametrically on the delays between pulses *T*_{3} = *τ _{s}* −

Multidimensional signals are displayed in the frequency domain by Fourier transforming ${S}_{{\mathit{k}}_{s}}^{(3)}({T}_{3},{T}_{2},{T}_{1})$ with respect to the time intervals between the pulses.

$${S}_{{\mathit{k}}_{s}}^{(3)}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})={\int}_{0}^{\infty}\text{d}{T}_{3}{\int}_{0}^{\infty}\text{d}{T}_{2}{\int}_{0}^{\infty}\text{d}{T}_{1}{\text{e}}^{i{\mathrm{\Omega}}_{3}{T}_{3}+i{\mathrm{\Omega}}_{2}{T}_{2}+i{\mathrm{\Omega}}_{1}{T}_{1}}{S}_{{\mathit{k}}_{s}}^{(3)}({T}_{3},{T}_{2},{T}_{1}).$$

(41)

Often, one uses a mixed time-frequency representation by performing double-Fourier transform. This gives, e.g., *S*(Ω_{3}, *T*_{2}, Ω_{1}), etc.

Assuming that all four pulses are temporally well-separated, the integrations can be carried out using the response function eqs 32–39.^{127} Neglecting population transport (*K _{ee}*

$${S}_{{\mathit{k}}_{\text{I}},\text{i}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-\sum _{{ee}^{\prime}}({\mathit{\mu}}_{{e}^{\prime}g}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{{e}^{\prime}g}-{\omega}_{s}))({\mathit{\mu}}_{eg}\xb7{\mathbb{E}}_{3}({\omega}_{eg}-{\omega}_{3}))\times \frac{({\mathit{\mu}}_{{e}^{\prime}g}\xb7{\mathbb{E}}_{2}({\omega}_{{e}^{\prime}g}-{\omega}_{2}))({\mathit{\mu}}_{eg}^{\ast}\xb7{\mathbb{E}}_{1}^{\ast}({\omega}_{eg}-{\omega}_{1}))}{{\hslash}^{3}({\mathrm{\Omega}}_{3}-{\xi}_{{e}^{\prime}g})({\mathrm{\Omega}}_{2}-{\xi}_{{e}^{\prime}e})({\mathrm{\Omega}}_{1}-{\xi}_{ge})},$$

(42)

$${S}_{{\mathit{k}}_{\text{I}},\text{ii}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-\sum _{{ee}^{\prime}}({\mathit{\mu}}_{{e}^{\prime}g}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{{e}^{\prime}g}-{\omega}_{s}))({\mathit{\mu}}_{{e}^{\prime}g}\xb7{\mathbb{E}}_{3}({\omega}_{{e}^{\prime}g}-{\omega}_{3}))\times \frac{({\mathit{\mu}}_{eg}\xb7{\mathbb{E}}_{2}({\omega}_{eg}-{\omega}_{2}))({\mathit{\mu}}_{eg}^{\ast}\xb7{\mathbb{E}}_{1}^{\ast}({\omega}_{eg}-{\omega}_{1}))}{{\hslash}^{3}({\mathrm{\Omega}}_{3}-{\xi}_{{e}^{\prime}g})({\mathrm{\Omega}}_{2}-i\eta )({\mathrm{\Omega}}_{1}-{\xi}_{ge})},$$

(43)

$${S}_{{\mathit{k}}_{\text{I}},\text{iii}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=\sum _{{\mathit{fee}}^{\prime}}({\mathit{\mu}}_{fe}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{fe}-{\omega}_{s}))({\mathit{\mu}}_{f{e}^{\prime}}\xb7{\mathbb{E}}_{3}({\omega}_{f{e}^{\prime}}-{\omega}_{3}))\times \frac{({\mathit{\mu}}_{{e}^{\prime}g}\xb7{\mathbb{E}}_{2}({\omega}_{{e}^{\prime}g}-{\omega}_{2}))({\mathit{\mu}}_{eg}^{\ast}\xb7{\mathbb{E}}_{1}^{\ast}({\omega}_{eg}-{\omega}_{1}))}{{\hslash}^{3}({\mathrm{\Omega}}_{3}-{\xi}_{fe})({\mathrm{\Omega}}_{2}-{\xi}_{{e}^{\prime}e})({\mathrm{\Omega}}_{1}-{\xi}_{ge})},$$

(44)

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

$$\mathbb{E}(\omega )=\int \text{d}texp(i\omega t)\mathbb{E}(t)$$

(45)

are centered around *ω* = 0.

For the *k*_{II} signal, we similarly obtain

$${S}_{{\mathit{k}}_{\text{II}},\text{iv}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-\sum _{{ee}^{\prime}}({\mathit{\mu}}_{eg}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{eg}-{\omega}_{s}))({\mathit{\mu}}_{{e}^{\prime}g}\xb7{\mathbb{E}}_{3}({\omega}_{{e}^{\prime}g}-{\omega}_{3}))\times \frac{({\mathit{\mu}}_{{e}^{\prime}g}^{\ast}\xb7{\mathbb{E}}_{2}^{\ast}({\omega}_{{e}^{\prime}g}-{\omega}_{2}))({\mathit{\mu}}_{eg}\xb7{\mathbb{E}}_{1}({\omega}_{eg}-{\omega}_{1}))}{{\hslash}^{3}({\mathrm{\Omega}}_{3}-{\xi}_{eg})({\mathrm{\Omega}}_{2}-{\xi}_{{ee}^{\prime}})({\mathrm{\Omega}}_{1}-{\xi}_{eg})},$$

(46)

$${S}_{{\mathit{k}}_{\text{II}},\text{v}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-\sum _{{ee}^{\prime}}({\mathit{\mu}}_{{e}^{\prime}g}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{{e}^{\prime}g}-{\omega}_{s}))({\mathit{\mu}}_{{e}^{\prime}g}\xb7{\mathbb{E}}_{3}({\omega}_{{e}^{\prime}g}-{\omega}_{3}))\times \frac{({\mathit{\mu}}_{eg}^{\ast}\xb7{\mathbb{E}}_{2}^{\ast}({\omega}_{eg}-{\omega}_{2}))({\mathit{\mu}}_{eg}\u2022{\mathbb{E}}_{1}({\omega}_{eg}-{\omega}_{1}))}{{\hslash}^{3}({\mathrm{\Omega}}_{3}-{\xi}_{{e}^{\prime}g})({\mathrm{\Omega}}_{2}+i\eta )({\mathrm{\Omega}}_{1}-{\xi}_{eg})},$$

(47)

$${S}_{{\mathit{k}}_{\text{II}},\text{vi}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=\sum _{{\mathit{fee}}^{\prime}}({\mathit{\mu}}_{f{e}^{\prime}}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{f{e}^{\prime}}-{\omega}_{s}))({\mathit{\mu}}_{fe}\xb7{\mathbb{E}}_{3}({\omega}_{fe}-{\omega}_{3}))\times \frac{({\mathit{\mu}}_{{e}^{\prime}g}^{\ast}\xb7{\mathbb{E}}_{2}^{\ast}({\omega}_{{e}^{\prime}g}-{\omega}_{2}))({\mathit{\mu}}_{eg}\xb7{\mathbb{E}}_{1}({\omega}_{eg}-{\omega}_{1}))}{{\hslash}^{3}({\mathrm{\Omega}}_{3}-{\xi}_{f{e}^{\prime}})({\mathrm{\Omega}}_{2}-{\xi}_{{ee}^{\prime}})({\mathrm{\Omega}}_{1}-{\xi}_{eg})},$$

(48)

Finally, the *k*_{III} signal is given by

$${S}_{{\mathit{k}}_{\text{III}},\text{vii}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=\sum _{{\mathit{fee}}^{\prime}}({\mathit{\mu}}_{f{e}^{\prime}}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{f{e}^{\prime}}-{\omega}_{s}))({\mathit{\mu}}_{{e}^{\prime}g}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{{e}^{\prime}g}-{\omega}_{3}))\times \frac{({\mathit{\mu}}_{fe}\xb7{\mathbb{E}}_{2}({\omega}_{fe}-{\omega}_{2}))({\mathit{\mu}}_{eg}\xb7{\mathbb{E}}_{1}({\omega}_{eg}-{\omega}_{1}))}{{\hslash}^{3}({\mathrm{\Omega}}_{3}-{\xi}_{f{e}^{\prime}})({\mathrm{\Omega}}_{2}-{\xi}_{fg})({\mathrm{\Omega}}_{1}-{\xi}_{eg})},$$

(49)

and

$${S}_{{\mathit{k}}_{\text{III}},\text{viii}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-\sum _{{e}^{\prime}fe}({\mathit{\mu}}_{{e}^{\prime}g}^{\ast}\xb7{\mathbb{E}}_{s}^{\ast}({\omega}_{{e}^{\prime}g}-{\omega}_{s}))({\mathit{\mu}}_{f{e}^{\prime}}^{\ast}\xb7{\mathbb{E}}_{3}^{\ast}({\omega}_{f{e}^{\prime}}-{\omega}_{3}))\times \frac{({\mathit{\mu}}_{fe}\xb7{\mathbb{E}}_{2}({\omega}_{fe}-{\omega}_{2}))({\mathit{\mu}}_{eg}\xb7{\mathbb{E}}_{1}({\omega}_{eg}-{\omega}_{1}))}{{\hslash}^{3}({\mathrm{\Omega}}_{3}-{\xi}_{{e}^{\prime}g})({\mathrm{\Omega}}_{2}-{\xi}_{fg})({\mathrm{\Omega}}_{1}-{\xi}_{eg})}.$$

(50)

Note that resonances occur at both positive and negative Ω* _{j}*.

Population relaxation may be added, as done in eqs 32–39. However, the Ω_{2} dependence is then more complicated, and it is usually preferable to display the signal *S*(Ω_{3}, *T*_{2}, Ω_{1}) in the time domain as a function of *T*_{2}.

Equations 42–50 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 &(*ω*) = 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.

The calculation of double-exciton states requires the computationally expensive diagonalization of an × matrix, where = ( ± 1)/2 is the number of double-exciton states, and 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}^{,}^{106}^{–}^{108}^{,}^{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 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
${\phi}_{m}^{(g)}$, the single-excited state
${\phi}_{m}^{(e)}$, and the double-excited state
${\phi}_{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,

$${\mathrm{\Phi}}_{g}={\prod}_{m}^{\mathbb{N}}{\phi}_{m}^{(g)}.$$

(51)

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

$${\mathrm{\Phi}}_{{e}_{m}}={\phi}_{m}^{(e)}{\prod}_{k}^{k\ne m}{\phi}_{k}^{(g)}.$$

(52)

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

$${\mathrm{\Phi}}_{{f}_{mm}}={\phi}_{m}^{(f)}{\prod}_{k}^{k\ne m}{\phi}_{k}^{(g)}$$

(53)

or when two chromophores are singly excited:

$${\mathrm{\Phi}}_{{f}_{mn}}={\phi}_{m}^{(e)}{\phi}_{n}^{(e)}{\prod}_{k}^{k\ne m,k\ne n}{\phi}_{k}^{(g)}.$$

(54)

Overall, our model has singly excited states, and = ( + 1)/2 doubly-excited states.

We next introduce an exciton-creation operator on molecule *m*:^{131}

$${\widehat{B}}_{m}^{\u2020}=\mid {\phi}_{m}^{(e)}\rangle \langle {\phi}_{m}^{(g)}\mid +\sqrt{2}\mid {\phi}_{m}^{(f)}\rangle \langle {\phi}_{m}^{(e)}\mid .$$

(55)

These operators have the following properties: single-exciton states are created from the ground state ${\mathrm{\Phi}}_{{e}_{m}}={\widehat{B}}_{m}^{\u2020}{\mathrm{\Phi}}_{g}$, and double-excitons are obtained either by creating two excitations on different molecules, ${\mathrm{\Phi}}_{{f}_{mn};m\ne n}={\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}^{\u2020}{\mathrm{\Phi}}_{g}$, or on the same molecule, ${\mathrm{\Phi}}_{{f}_{mm}}={2}^{-1/2}{\widehat{B}}_{m}^{\u20202}{\mathrm{\Phi}}_{g}$ (the √2 factor is introduced to resemble bosonic exciton properties ( ${\widehat{B}}^{\u2020}\mid n-1\rangle =\sqrt{n}\mid n\rangle $) within the single- and double-exciton space).

The Hermitian-conjugate, annihilation, operator will be denoted * _{m}*. Operators corresponding to different chromophores commute. The commutation relations of these operators are (see Appendix C)

$$[{\widehat{B}}_{n},{\widehat{B}}_{m}^{\u2020}]={\delta}_{mn}\left(1-\frac{3}{2}{\widehat{B}}_{m}^{\u20202}{\widehat{B}}_{m}^{2}\right).$$

(56)

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, ^{†}^{†}^{†}Φ_{g}, will not be considered.

We assume the following model Hamiltonian written in a normally ordered form (i.e., all ^{†} are to the left of )^{68}

$${\widehat{H}}_{0}=\hslash \sum _{m}{\epsilon}_{m}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}+\hslash \sum _{mn}^{m\ne n}{J}_{mn}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}+\hslash \sum _{\mathit{mnkl}}{U}_{mn,kl}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}^{\u2020}{\widehat{B}}_{k}{\widehat{B}}_{l},$$

(57)

where *ε _{n}* is the single-excitation energy of molecule

A simplified form of eq 57,

$${\widehat{H}}_{0}=\hslash \sum _{m}{\epsilon}_{m}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}+\hslash \sum _{mn}^{m\ne n}{J}_{mn}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}+\hslash \sum _{m}\frac{{\mathrm{\Delta}}_{m}}{2}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}{\widehat{B}}_{m}+\hslash \sum _{mn}^{m\ne n}\frac{{\mathbb{K}}_{mn}}{2}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}^{\u2020}{\widehat{B}}_{m}{\widehat{B}}_{n}$$

(58)

captures the essential physical processes. The diagonal term Δ* _{m}* = 2

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

$$\widehat{\mathit{\mu}}=\sum _{m}[{\mathit{\mu}}_{m}{\widehat{B}}_{m}^{\u2020}+{\mathit{\mu}}_{m}^{(2)}{\widehat{B}}_{m}^{\u20202}{\widehat{B}}_{m}+h.c].$$

(59)

Here, *h.c* denotes the Hermitian conjugate. The transition dipole Φ* _{g}*|

$${\widehat{H}}^{\prime}(t)=\hslash \sum _{m,j}({\mu}_{m,j}^{-}(t){\widehat{B}}_{m}^{\u2020}+{\mu}_{m,j}^{+}(t){\widehat{B}}_{m}+{\mu}_{m,j}^{(2)-}(t){\widehat{B}}_{m}^{\u20202}{\widehat{B}}_{m}+{\mu}_{m,j}^{(2)+}(t){\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}^{2}),$$

(60)

where

$${\mu}_{m,j}^{-}(t)=-{\hslash}^{-1}{\mathit{\mu}}_{m}\xb7{\mathbb{E}}_{j}(t-{\tau}_{j})exp(i{\mathit{k}}_{j}\mathit{r}-i{\omega}_{j}t)$$

(61)

and

$${\mu}_{m,j}^{(2)-}(t)=-{\hslash}^{-1}{\mathit{\mu}}_{m}^{(2)}\xb7{\mathbb{E}}_{j}(t-{\tau}_{j})exp(i{\mathit{k}}_{j}\mathit{r}-i{\omega}_{j}t)$$

(62)

are the time-dependent transition amplitudes with *μ ^{−}* (

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

$$i\hslash \frac{\text{d}\widehat{A}}{\text{d}t}=[\widehat{A},\widehat{H}].$$

(63)

Taking *Â* = * _{m}* and

$$i\frac{\text{d}}{\text{d}t}{\widehat{B}}_{m}=\sum _{{m}^{\prime}}{h}_{m{m}^{\prime}}{\widehat{B}}_{{m}^{\prime}}+\sum _{{m}^{\prime}kl}{V}_{m{m}^{\prime},kl}{\widehat{B}}_{{m}^{\prime}}^{\u2020}{\widehat{B}}_{k}{\widehat{B}}_{l}+\sum _{j}[{\mu}_{m,j}^{-}+2{\mu}_{m,j}^{(2)-}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}+{\mu}_{m,j}^{(2)+}{\widehat{B}}_{m}{\widehat{B}}_{m}],$$

(64)

$$i\frac{\text{d}}{\text{d}t}{\widehat{B}}_{m}{\widehat{B}}_{n}=\sum _{{m}^{\prime}{n}^{\prime}}{h}_{mn,{m}^{\prime}{n}^{\prime}}^{(Y)}{\widehat{B}}_{{m}^{\prime}}{\widehat{B}}_{{n}^{\prime}}+\sum _{{m}^{\prime}kl}[{V}_{m{m}^{\prime},kl}{\widehat{B}}_{{m}^{\prime}}^{\u2020}{\widehat{B}}_{k}{\widehat{B}}_{l}{\widehat{B}}_{n}+{V}_{n{m}^{\prime},kl}{\widehat{B}}_{{n}^{\prime}}^{\u2020}{\widehat{B}}_{k}{\widehat{B}}_{l}{\widehat{B}}_{m}]+\sum _{j}[{\mu}_{m,j}^{-}{\widehat{B}}_{n}+{\mu}_{n,j}^{-}{\widehat{B}}_{m}+2{\delta}_{mn}{\mu}_{m,j}^{(2)-}{\widehat{B}}_{m}+2{\mu}_{n,j}^{(2)-}{\widehat{B}}_{n}^{\u2020}{\widehat{B}}_{n}{\widehat{B}}_{m}+2{\mu}_{m,j}^{(2)-}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}{\widehat{B}}_{n}+2{\mu}_{n,j}^{(2)+}{\widehat{B}}_{n}{\widehat{B}}_{n}{\widehat{B}}_{m}],$$

(65)

with

$${h}_{mn}={J}_{mn}(1-{\delta}_{mn})+{\delta}_{mn}{\epsilon}_{m},$$

(66)

$${h}_{mn,{m}^{\prime}{n}^{\prime}}^{(Y)}={h}_{m{m}^{\prime}}{\delta}_{n{n}^{\prime}}+{h}_{n{n}^{\prime}}{\delta}_{m{m}^{\prime}}+{V}_{mn,{m}^{\prime}{n}^{\prime}}$$

(67)

and *V _{mn}*

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, is first order and higher,

When the system is in a pure state, it can be described by a wave function (the expansion
$\mathrm{\Phi}={\sum}_{m}{c}_{m}{\widehat{B}}_{m}^{\u2020}\mid 0\rangle +{C}_{mn}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}^{\u2020}\mid 0\rangle $ is sufficient to represent the third-order response), and we have
$\langle {\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}\rangle =\langle {\widehat{B}}_{m}^{\u2020}\rangle \langle {\widehat{B}}_{n}\rangle $. We can further factorize
$\langle {\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}{\widehat{B}}_{k}\rangle ={B}_{m}^{\ast}{Y}_{nk}$ where we used *B* = *B* and *Y _{mn}* =

$$i\frac{\text{d}}{\text{d}t}{B}_{m}=\sum _{{m}^{\prime}}{h}_{m{m}^{\prime}}{B}_{{m}^{\prime}}+\sum _{{m}^{\prime}kl}{V}_{m{m}^{\prime},kl}{B}_{{m}^{\prime}}^{\ast}{Y}_{kl}+\sum _{j}[{\mu}_{m,j}^{-}+2{\mu}_{m,j}^{(2)-}{B}_{m}^{\ast}{B}_{m}+{\mu}_{m,j}^{(2)+}{Y}_{mm}],$$

(68)

$$i\frac{\text{d}}{\text{d}t}{Y}_{mn}=\sum _{{m}^{\prime}{n}^{\prime}}{h}_{mn,{m}^{\prime}{n}^{\prime}}^{(Y)}{Y}_{{m}^{\prime}{n}^{\prime}}+\sum _{j}[{\mu}_{m,j}^{-}{B}_{n}+{\mu}_{n,j}^{-}{B}_{m}+2{\delta}_{mn}{\mu}_{m,j}^{(2)-}{B}_{m}].$$

(69)

Quartic products ^{†2}^{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,

$${[{\widehat{B}}_{m},{\widehat{B}}_{n}^{\u2020}]}_{\text{boson}}={\delta}_{mn},$$

(70)

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 oscillators^{99}^{,}^{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 problems^{114} uses the full factorization of all normally ordered products = and ^{†} = ^{†}, ^{†} = ^{†}.^{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 ^{†} = ^{†}and ^{†} = ^{†}, while retaining the 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 ^{†} = ^{†} 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.

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).

In the linear regime, we get

$${\mathit{P}}^{(1)}(t)=\sum _{m}{\mathit{\mu}}_{m}^{\ast}{B}_{m}^{(1)}(t)+c.c.$$

(71)

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

$$i\frac{\text{d}}{\text{d}t}{B}_{m}^{(1)}=\sum _{{m}^{\prime}}{h}_{m{m}^{\prime}}{B}_{{m}^{\prime}}^{(1)},$$

(72)

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

$$\frac{\text{d}}{\text{d}t}{G}_{mn}=-i\sum _{{m}^{\prime}}{h}_{m{m}^{\prime}}{G}_{{m}^{\prime}n}+\delta (t){\delta}_{mn},$$

(73)

whose solution is *G*(*t*) = θ(*t*) exp(−*iht*), where *h* is the matrix with elements *h _{mn}*. The solution of eq 72 is given by
${B}_{m}^{(1)}(t)={G}_{mn}(t){B}_{n}^{(1)}(0)$.

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

$$i\frac{\text{d}}{\text{d}t}{B}_{m}^{(1)}=\sum _{{m}^{\prime}}{h}_{m{m}^{\prime}}{B}_{{m}^{\prime}}^{(1)}+\sum _{j}{\mu}_{m,j}^{-}.$$

(74)

Substituting the Green’s function solution of eq 74,

$${B}_{m}^{(1)}(t)=-i{\int}_{-\infty}^{t}\text{d}\tau \sum _{n,j}{G}_{mn}(t-\tau ){\mu}_{n,j}^{-}(\tau )$$

(75)

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

$${R}^{(1)}(t)={\left(\frac{i}{\hslash}\right)}^{1}\sum _{mn}{\mathit{\mu}}_{m}^{\ast}{\mathit{\mu}}_{n}{G}_{mn}(t)+c.c.$$

(76)

Linear response techniques are surveyed in Appendix A.

The third-order induced polarization is

$${\mathit{P}}^{(3)}(t)=\sum _{m}{\mathit{\mu}}_{m}^{\ast}{\langle {\widehat{B}}_{m}\rangle}^{(3)}+{\mathit{\mu}}_{m}^{(2)\ast}{\langle {\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}{\widehat{B}}_{m}\rangle}^{(3)}+c.c.$$

(77)

Below, we present two levels of approximation for the *k*_{I}, *k*_{II}, and *k*_{III} techniques: the MFA and the coherent exciton dynamics (CED) limit.

At this level of theory, we make the factorizations = . Equation 68 then becomes

$$i\frac{\text{d}}{\text{d}t}{B}_{m}=\sum _{{m}^{\prime}}{h}_{m{m}^{\prime}}{B}_{{m}^{\prime}}+\sum _{{m}^{\prime}kl}{V}_{m{m}^{\prime},kl}{B}_{{m}^{\prime}}^{\ast}{B}_{k}{B}_{l}+\sum _{j}{\mu}_{m,j}^{-}+2{\mu}_{m,j}^{(2)-}{B}_{m}^{\ast}{B}_{m}+{\mu}_{m,j}^{(2)+}{B}_{m}{B}_{m}$$

(78)

and the polarization (eq 77) reduces to

$${\mathit{P}}^{(3)}(r,t)=\sum _{m}[{\mathit{\mu}}_{m}^{\ast}{B}_{m}^{(3)}(t)+{\mathit{\mu}}_{m}^{(2)\ast}{B}_{m}^{(1)\ast}(t){B}_{m}^{(1)}(t){B}_{m}^{(1)}(t)+c.c.].$$

(79)

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

$${B}_{m}^{(3)}(t)=-i{\int}_{-\infty}^{t}\text{d}\tau \phantom{\rule{0.16667em}{0ex}}{G}_{mn}(t-\tau )\times [{V}_{n{n}^{\prime},kl}{B}_{{n}^{\prime}}^{(1)\ast}(\tau ){B}_{k}^{(1)}(\tau ){B}_{l}^{(1)}(\tau )+2{\mu}_{n,j}^{(2)-}(\tau ){B}_{n}^{(1)\ast}(\tau ){B}_{n}^{(1)}(\tau )+{\mu}_{n,j}^{(2)+}(\tau ){B}_{n}^{(1)}(\tau ){B}_{n}^{(1)}(\tau )].$$

(80)

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

The *k*_{I} 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

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{MFA})}({t}_{3},{t}_{2},{t}_{1})=2i{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}{\mathit{\mu}}_{{n}_{2}}{\mathit{\mu}}_{{n}_{1}}^{\ast}{\int}_{0}^{{t}_{3}}\text{d}\tau {V}_{{n}_{4}^{\prime}{n}_{1}^{\prime},{n}_{3}^{\prime}{n}_{2}^{\prime}}\times {G}_{{n}_{4}{n}_{4}^{\prime}}({t}_{3}-\tau ){G}_{{n}_{3}^{\prime}{n}_{3}}(\tau ){G}_{{n}_{2}^{\prime}{n}_{2}}({t}_{2}+\tau ){G}_{{n}_{1}^{\prime}{n}_{1}}^{\ast}({t}_{1}+{t}_{2}+\tau ).$$

(81)

For the *k*_{II} technique, we similarly have

$${R}_{{\mathit{k}}_{\text{II}}}^{(\text{MFA})}({t}_{3},{t}_{2},{t}_{1})=2i{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}{\mathit{\mu}}_{{n}_{2}}^{\ast}{\mathit{\mu}}_{{n}_{1}}{\int}_{0}^{{t}_{3}}\text{d}\tau {V}_{{n}_{4}^{\prime}{n}_{2}^{\prime},{n}_{3}^{\prime}{n}_{1}^{\prime}}\times {G}_{{n}_{4}{n}_{4}^{\prime}}({t}_{3}-\tau ){G}_{{n}_{3}^{\prime}{n}_{3}}(\tau ){G}_{{n}_{2}^{\prime}{n}_{2}}^{\ast}({t}_{2}+\tau ){G}_{{n}_{1}^{\prime}{n}_{1}}({t}_{1}+{t}_{2}+\tau ),$$

(82)

and for *k*_{III},

$${R}_{{\mathit{k}}_{\text{III}}}^{(\text{MFA})}({t}_{3},{t}_{2},{t}_{1})=2i{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}^{\ast}{\mathit{\mu}}_{{n}_{2}}{\mathit{\mu}}_{{n}_{1}}{\int}_{0}^{{t}_{3}}\text{d}\tau {V}_{{n}_{4}^{\prime}{n}_{3}^{\prime},{n}_{2}^{\prime}{n}_{1}^{\prime}}\times {G}_{{n}_{4}{n}_{4}^{\prime}}({t}_{3}-\tau ){G}_{{n}_{3}^{\prime}{n}_{3}}^{\ast}(\tau ){G}_{{n}_{2}^{\prime}{n}_{2}}({t}_{2}+\tau ){G}_{{n}_{1}^{\prime}{n}_{1}}({t}_{1}+{t}_{2}+\tau ).$$

(83)

In these expressions, the *k** _{j}* pulse interacts with chromophore

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

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 *Y _{mn}*(

$$\frac{\text{d}}{\text{d}t}{\mathbb{G}}_{mn,kl}=-i\sum _{{m}^{\prime}{n}^{\prime}}{h}_{mn,{m}^{\prime}{n}^{\prime}}^{(Y)}{\mathbb{G}}_{{m}^{\prime}{n}^{\prime},kl}+\delta (t){\delta}_{mk}{\delta}_{nl}.$$

(84)

The solution, (*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

$${Y}_{mn}^{(2)}(t)=-2i{\int}_{-\infty}^{t}\text{d}\tau \phantom{\rule{0.16667em}{0ex}}{\mathbb{G}}_{mn,{m}^{\prime}{n}^{\prime}}(t-\tau )\times [{\mu}_{{m}^{\prime},j}^{-}(\tau ){B}_{{n}^{\prime}}^{(1)}(\tau )+{\delta}_{{m}^{\prime}{n}^{\prime}}{\mu}_{{m}^{\prime},j}^{(2)-}(\tau ){B}_{{m}^{\prime}}^{(1)}(\tau )],$$

(85)

where we have further used the relation , which results from boson symmetry. Equation 68 then gives

$${B}_{m}^{(3)}(t)=-i{\int}_{-\infty}^{t}\text{d}\tau \phantom{\rule{0.16667em}{0ex}}{G}_{mn}(t-\tau )[{V}_{n{n}^{\prime},kl}{B}_{n}^{(1)\ast}(\tau ){Y}_{kl}^{(2)}(\tau )+2{\mu}_{n,j}^{(2)-}(\tau ){B}_{n}^{(1)\ast}(\tau ){B}_{n}^{(1)}(\tau )+{\mu}_{n,j}^{(2)+}(\tau ){Y}_{nn}^{(2)}(\tau )].$$

(86)

The polarization (eq 77) is given by

$${\mathit{P}}^{(3)}(\mathit{r},t)=\sum _{m}[{\mathit{\mu}}_{m}^{\ast}{B}_{m}^{(3)}(t)+{\mathit{\mu}}_{m}^{(2)\ast}{B}_{m}^{(1)\ast}(t){Y}_{mm}^{(2)}(t)+c.c.].$$

(87)

For
${\mathit{\mu}}_{m}^{(2)}=0$, we get for the *k*_{I} response function

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{CED})}({t}_{3},{t}_{2},{t}_{1})=2i{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}{\mathit{\mu}}_{{n}_{2}}{\mathit{\mu}}_{{n}_{1}}^{\ast}{\int}_{0}^{{t}_{3}}\text{d}\tau {V}_{{n}_{4}^{\prime}{n}_{1}^{\prime},{n}_{3}^{\prime}{n}_{2}^{\prime}}\times {G}_{{n}_{4}{n}_{4}^{\prime}}({t}_{3}-\tau ){\mathbb{G}}_{{n}_{3}^{\prime}{n}_{2}^{\prime},{n}_{3}{n}_{2}^{\u2033}}(\tau ){G}_{{n}_{2}^{\u2033}{n}_{2}}({t}_{2}){G}_{{n}_{1}^{\prime}{n}_{1}}^{\ast}({t}_{1}+{t}_{2}+\tau ).$$

(88)

We recast this result using the exciton scattering matrix (see Appendix E). Substituting eq 324 in eq 88 with ${\mathbb{G}}_{mn,kl}^{(0)}(t)\equiv {G}_{mk}(t){G}_{nl}(t)$, we get

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{CED})}({t}_{3},{t}_{2},{t}_{1})=2{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}{\mathit{\mu}}_{{n}_{2}}{\mathit{\mu}}_{{n}_{1}}^{\ast}\times {\int}_{0}^{{t}_{3}}\text{d}\tau {\int}_{0}^{\tau}\text{d}{\tau}^{\prime}{G}_{{n}_{4}{n}_{4}^{\prime}}({t}_{3}-\tau ){\mathrm{\Gamma}}_{{n}_{4}^{\prime}{n}_{1}^{\prime},{n}_{3}^{\prime}{n}_{2}^{\prime}}(\tau -{\tau}^{\prime})\times {G}_{{n}_{3}^{\prime}{n}_{3}}({\tau}^{\prime}){G}_{{n}_{2}^{\prime}{n}_{2}}({t}_{2}+{\tau}^{\prime}){G}_{{n}_{1}^{\prime}{n}_{1}}^{\ast}({t}_{1}+{t}_{2}+\tau ).$$

(89)

Similarly, we obtain for *k*_{II}

$${R}_{{\mathit{k}}_{\text{II}}}^{(\text{CED})}({t}_{3},{t}_{2},{t}_{1})=2{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}{\mathit{\mu}}_{{n}_{2}}^{\ast}{\mathit{\mu}}_{{n}_{1}}\times {\int}_{0}^{{t}_{3}}\text{d}\tau {\int}_{0}^{\tau}\text{d}{\tau}^{\prime}{G}_{{n}_{4}{n}_{4}^{\prime}}({t}_{3}-\tau ){\mathrm{\Gamma}}_{{n}_{4}^{\prime}{n}_{2}^{\prime},{n}_{3}^{\prime}{n}_{1}^{\prime}}(\tau -{\tau}^{\prime})\times {G}_{{n}_{3}^{\prime}{n}_{3}}({\tau}^{\prime}){G}_{{n}_{2}^{\prime}{n}_{2}}^{\ast}({t}_{2}+\tau ){G}_{{n}_{1}^{\prime}{n}_{1}}({t}_{1}+{t}_{2}+{\tau}^{\prime})$$

(90)

and for *k*_{III}:

$${R}_{{\mathit{k}}_{\text{III}}}^{(\text{CED})}({t}_{3},{t}_{2},{t}_{1})=2{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}^{\ast}{\mathit{\mu}}_{{n}_{2}}{\mathit{\mu}}_{{n}_{1}}\times {\int}_{0}^{{t}_{3}+{t}_{2}}\text{d}\tau {\int}_{0}^{\tau}\text{d}{\tau}^{\prime}{G}_{{n}_{4}{n}_{4}^{\prime}}({t}_{3}-\tau ){\mathrm{\Gamma}}_{{n}_{4}^{\prime}{n}_{3}^{\prime},{n}_{2}^{\prime}{n}_{1}^{\prime}}(\tau -{\tau}^{\prime})\times {G}_{{n}_{3}^{\prime}{n}_{3}}^{\ast}(\tau ){G}_{{n}_{2}^{\prime}{n}_{2}}({t}_{2}+{\tau}^{\prime}){G}_{{n}_{1}^{\prime}{n}_{1}}({t}_{1}+{t}_{2}+{\tau}^{\prime}).$$

(91)

Γ_{n4n3, n2n1} is the exciton scattering matrix in two-exciton space. It is a tetradic matrix, whose indices *n*_{1} and *n*_{2} represent incoming excitons and *n*_{3} and *n*_{4} are outgoing excitons. The scattering is induced by the coupling *V*. The scattering matrix is given by

$$\mathrm{\Gamma}(\omega )={[1-V{\mathbb{G}}^{(0)}(\omega )]}^{-1}V,$$

(92)

where

$${\mathbb{G}}_{{n}_{1}{n}_{2},{n}_{3}{n}_{4}}^{(0)}(\omega )=i\int \frac{\text{d}{\omega}^{\prime}}{2\pi}{G}_{{n}_{1}{n}_{3}}({\omega}^{\prime}){G}_{{n}_{2}{n}_{4}}(\omega -{\omega}^{\prime}).$$

(93)

Closed frequency-domain expressions of the signals corresponding to eqs 89–91 are obtained from eqs 127–129 by neglecting transport and setting ${\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}(\omega )={\delta}_{{e}_{4}{e}_{2}}{\delta}_{{e}_{3}{e}_{1}}{[\omega -({\epsilon}_{{e}_{2}}-{\epsilon}_{{e}_{1}})+i({\gamma}_{{e}_{2}}+{\gamma}_{{e}_{1}})]}^{-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.

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}

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

$$\widehat{H}={\widehat{H}}_{0}+{\widehat{H}}_{B}+{\widehat{H}}_{SB}+{\widehat{H}}^{\prime}(t).$$

(94)

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

We shall use a harmonic model for the bath,

$${\widehat{H}}_{B}=\sum _{\alpha}\hslash {w}_{\alpha}({\widehat{a}}_{\alpha}^{\u2020}{\widehat{a}}_{\alpha}+1/2),$$

(95)

where α runs over the bath oscillators and ${\widehat{a}}_{\alpha}^{\u2020}({\widehat{a}}_{\alpha})$ are Bose creation (annihilation) operators that satisfy $[{\widehat{a}}_{\alpha},{\widehat{a}}_{\beta}^{\u2020}]={\delta}_{\alpha \beta}$. Note that

$${\widehat{a}}_{\alpha}^{\u2020}={(2{m}_{\alpha}\hslash {w}_{\alpha})}^{-1/2}({m}_{\alpha}{w}_{\alpha}{\widehat{Q}}_{\alpha}-i{\widehat{P}}_{\alpha}),$$

(96)

$${\widehat{a}}_{\alpha}={(2{m}_{\alpha}\hslash {w}_{\alpha})}^{-1/2}({m}_{\alpha}{w}_{\alpha}{\widehat{Q}}_{\alpha}+i{\widehat{P}}_{\alpha}),$$

(97)

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

$${\widehat{H}}_{SB}=\hslash \sum _{mn\alpha}{d}_{mn,\alpha}{\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}({\widehat{a}}_{\alpha}^{\u2020}+{\widehat{a}}_{\alpha}),$$

(98)

where the coupling parameters *d _{mn}*

Because of the coupling with the bath, the factorization ^{†} into ^{†} 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,
${N}_{mn}=\langle {\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{n}\rangle $ and
${Z}_{\mathit{kmn}}=\langle {\widehat{B}}_{m}^{\u2020}{\widehat{B}}_{m}{\widehat{B}}_{n}\rangle $. 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

$$i\frac{\text{d}}{\text{d}t}{B}_{m}={[hB]}_{m}-i{[\gamma B]}_{m}+\sum _{{m}^{\prime}kl}{V}_{m{m}^{\prime},kl}{Z}_{{m}^{\prime},kl}+\sum _{j}[{\mu}_{m,j}^{-}+2{\mu}_{m,j}^{(2)-}{N}_{mm}+2{\mu}_{m,j}^{(2)+}{Y}_{mm}],$$

(99)

$$i\frac{\text{d}}{\text{d}t}{Y}_{mn}={[{h}^{(Y)}Y]}_{mn}-i{[{\gamma}^{(Y)}Y]}_{mn}+\sum _{j}{[{\mu}_{j}^{-}\otimes B+B\otimes {\mu}_{j}^{-}]}_{mn}+2{\delta}_{mn}{\mu}_{m,j}^{(2)-}{B}_{m},$$

(100)

$$i\frac{\text{d}}{\text{d}t}{N}_{mn}={[\stackrel{\sim}{\mathbb{L}}N]}_{mn}-i{[\stackrel{\sim}{K}N]}_{mn}+\sum _{j}-{[{\mu}_{j}^{+}\otimes B]}_{mn}+{[{B}^{\ast}\otimes {\mu}_{j}^{-}]}_{mn},$$

(101)

$$i\frac{\text{d}}{\text{d}t}{Z}_{\mathit{kmn}}=\sum _{{m}^{\prime}{n}^{\prime}}{({h}^{(Y)}-i{\gamma}^{(Y)})}_{mn,{m}^{\prime}{n}^{\prime}}{Z}_{k{m}^{\prime}{n}^{\prime}}-\sum _{{k}^{\prime}}({h}_{{k}^{\prime}k}+i{\gamma}_{{k}^{\prime}k}^{\ast}){Z}_{{k}^{\prime}mn}+\sum _{j}[{N}_{kn}{\mu}_{m,j}^{-}+{N}_{km}{\mu}_{n,j}^{-}-{\mu}_{k,j}^{+}{Y}_{mn}+2{\delta}_{mn}{\mu}_{m,j}^{(2)-}{N}_{km}],$$

(102)

where we have used = *h _{nn}*

Some useful approximations may be used for the relaxation matrices.
${\gamma}_{mn,kl}^{(Y)}$ will be represented approximately as
${\gamma}_{mn,kl}^{(Y)}\approx {\gamma}_{mk}{\delta}_{nl}+{\gamma}_{nl}{\delta}_{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.

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

$$\frac{\text{d}}{\text{d}t}{N}_{{e}_{2}{e}_{1}}=i{\omega}_{{e}_{2}{e}_{1}}{N}_{{e}_{2}{e}_{1}}-\sum _{{e}_{2}^{\prime}{e}_{1}^{\prime}}{K}_{{e}_{1}{e}_{2},{e}_{1}^{\prime}{e}_{2}^{\prime}}{N}_{{e}_{2}^{\prime}{e}_{1}^{\prime}}.$$

(103)

Assuming that
${\omega}_{{e}_{2}{e}_{1}}\ne {\omega}_{{e}_{2}^{\prime}{e}_{1}^{\prime}}$ 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
${K}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}\equiv {K}_{{e}_{4}{e}_{4},{e}_{2}{e}_{2}}{\delta}_{{e}_{4}{e}_{3}}{\delta}_{{e}_{2}{e}_{1}}+{\gamma}_{{e}_{4}{e}_{3}}^{(N)}{\delta}_{{e}_{4}{e}_{2}}{\delta}_{{e}_{3}{e}_{1}}$ with
${\gamma}_{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,

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

$${\stackrel{.}{\rho}}_{{e}_{2}{e}_{1}}=-i{\omega}_{{e}_{2}{e}_{1}}^{\prime}{\rho}_{{e}_{2}{e}_{1}}+\sum _{c,d,{e}_{3},{e}_{4}}\left({V}_{{e}_{2}{e}_{3}}^{(c,d)}{\rho}_{{e}_{3}{e}_{4}}{({V}_{{e}_{1}{e}_{4}}^{(c,d)})}^{\ast}-\frac{1}{2}{\rho}_{{e}_{2}{e}_{3}}{({V}_{{e}_{4}{e}_{3}}^{(c,d)})}^{\ast}{V}_{{e}_{4}{e}_{1}}^{(c,d)}-\frac{1}{2}{({V}_{{e}_{3}{e}_{2}}^{(c,d)})}^{\ast}{V}_{{e}_{3}{e}_{4}}^{(c,d)}{\rho}_{{e}_{4}{e}_{1}}\right),$$

(104)

where we represent α by a pair of indices *c* and *d*. This will be convenient for connecting with the Redfield equation.
${V}_{{e}_{3}{e}_{4}}^{(c,d)}=\langle {\mathrm{\Psi}}_{{e}_{3}}\mid {\widehat{V}}^{(c,d)}\mid {\mathrm{\Psi}}_{{e}_{4}}\rangle $ where ^{(}^{c}^{,}^{d}^{)} has been projected into the one-exciton eigenstate, Ψ* _{e} = ψ_{em}*Φ

We next present the set of ’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

$${V}_{{e}_{1}{e}_{2}}^{(c,d)}={a}_{{e}_{1},{e}_{2}}{\delta}_{{e}_{1},c}{\delta}_{{e}_{2},d}(1-{\delta}_{c,d})$$

(105)

with *c* ≠ *d*. Since the rank of the *N*_{e2e1} matrix is , there are ( – 1) such *V*^{(}^{c}^{,} ^{d}^{)} matrices (these correspond to all off-diagonal elements of *N*_{e2e1}). The remaining matrices are associated with the diagonal elements (*c = d*):

$${V}_{{e}_{1}{e}_{2}}^{(c,d)}={b}_{{e}_{1}}^{(c)}{\delta}_{{e}_{1}{e}_{2}}{\delta}_{cd}.$$

(106)

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

$${\mid {a}_{c,d}\mid}^{2}={K}_{cc,dd}.$$

(107)

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

$${\widehat{\gamma}}_{{e}_{1},{e}_{2}}=-\frac{1}{2}\sum _{c}{({b}_{{e}_{1}}^{(c)}-{b}_{{e}_{2}}^{(c)})}^{2},$$

(108)

where _{e1,e2} is the pure-dephasing rate (as introduced in eq 29),

$$Re({K}_{{e}_{1}{e}_{2},{e}_{1}{e}_{2}})=\frac{{K}_{{e}_{1}{e}_{1},{e}_{1}{e}_{1}}+{K}_{{e}_{2}{e}_{2},{e}_{2}{e}_{2}}}{2}+{\widehat{\gamma}}_{{e}_{1},{e}_{2}}.$$

(109)

Multiple solutions for the *b*’s exist since eq 108 contains ( – 1) equations with ^{2} unknowns (the equations for *e*_{1} = *e*_{2} in eq 108 are satisfied for any choice of *b*’s). One particular solution is obtained from the rate expressions of Appendix G:

$${b}_{e}^{(c)}={\left(2\sum _{\beta}{\lambda}_{\beta}Re({\overline{M}}_{\beta}^{(+)}(0))\right)}^{1/2}{\psi}_{e,c}^{2}.$$

(110)

Finally, the frequency,
${\omega}_{{e}_{1}{e}_{2}}^{\prime}$ in eq 104 is related to ω _{e1e2} in eq 103 by a level-shift,

$${\omega}_{{e}_{1}{e}_{2}}^{\prime}={\omega}_{{e}_{1}{e}_{2}}-Im({K}_{{e}_{1}{e}_{2},{e}_{1}{e}_{2}}).$$

(111)

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.

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

$${\stackrel{.}{\rho}}_{mn}=-i{\omega}_{mn}{\rho}_{mn}-\sum _{kl}{K}_{mn,kl}{\rho}_{kl}.$$

(112)

Making the secular approximation in this basis, the Redfield relaxation operator becomes
${K}_{mn,kl}={\delta}_{mn}{\delta}_{kl}{K}_{mm,kk}^{(F)}+{\delta}_{mk}{\delta}_{nl}{\gamma}_{mn}^{(N)}$ where *K*^{(}^{F}^{)} is known as the Förster exciton-transfer rate. Note that
${\gamma}_{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}

$${K}_{AA,DD}^{(F)}=\frac{2\pi}{{\hslash}^{2}}{\mid {V}_{DA}\mid}^{2}\frac{\int \frac{\text{d}\omega}{{\omega}^{4}}{F}_{D}(\omega ){\epsilon}_{A}(\omega )}{\int \frac{\text{d}\omega}{{\omega}^{3}}{F}_{D}(\omega )\int \frac{\text{d}\omega}{\omega}{\epsilon}_{A}(\omega )},$$

(113)

where *F _{D}*(

$${V}_{DA}=\frac{1}{4\pi {\epsilon}_{0}\epsilon}\left(\frac{{\mathit{\mu}}_{D}\xb7{\mathit{\mu}}_{A}}{{\mid {\mathit{R}}_{DA}\mid}^{3}}-\frac{3({\mathit{\mu}}_{D}\xb7{\mathit{R}}_{DA})({\mathit{\mu}}_{A}\xb7{\mathit{R}}_{DA})}{{\mid {\mathit{R}}_{DA}\mid}^{5}}\right)$$

(114)

Here, ** μ** are the transition dipoles and

For small *R** _{DA}*, 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

We now derive closed QP expressions for the response functions in the exciton eigenstate basis. Exciton transport takes place during *t*_{2}. During the delay periods *t*_{1} and *t*_{3}, 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*) = (*t*)*N*(0) and *Z*(*t*) = (*t*)*Z*(0) where *N* = ^{†} and *Z* = ^{†} . is a tetradic matrix, while is sextic. Note that the *N* variable corresponds to the transpose of the density matrix in the single-exciton space, i.e., *N _{mn}* ρ

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({t}_{3},{t}_{2},{t}_{1})=2i{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}{\mathit{\mu}}_{{n}_{2}}{\mathit{\mu}}_{{n}_{1}}^{\ast}\times {\int}_{0}^{{t}_{3}}\text{d}{t}^{\prime}\phantom{\rule{0.16667em}{0ex}}{G}_{{n}_{4}{n}_{4}^{\prime}}({t}^{\prime}){V}_{{n}_{4}^{\prime}{n}_{1}^{\prime},{n}_{3}^{\prime}{n}_{2}^{\prime}}\times {\mathbb{G}}_{{n}_{1}^{\prime}{n}_{3}^{\prime}{n}_{2}^{\prime},{n}_{1}^{\u2033}{n}_{3}{n}_{2}^{\u2033}}^{(Z)}({t}_{3}-{t}^{\prime}){\mathbb{G}}_{{n}_{2}^{\u2033}{n}_{1}^{\u2033},{n}_{2}{n}_{1}^{\u2034}}^{(N)}({t}_{2}){G}_{{n}_{1}^{\u2034}{n}_{1}}^{\ast}({t}_{1}),$$

(115)

$${R}_{{\mathit{k}}_{\text{II}}}^{(\text{QP})}({t}_{3},{t}_{2},{t}_{1})=2i{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}{\mathit{\mu}}_{{n}_{2}}^{\ast}{\mathit{\mu}}_{{n}_{1}}\times {\int}_{0}^{{t}_{3}}\text{d}{t}^{\prime}{G}_{{n}_{4}{n}_{4}^{\prime}}({t}^{\prime}){V}_{{n}_{4}^{\prime}{n}_{2}^{\prime},{n}_{3}^{\prime}{n}_{1}^{\prime}}\times {\mathbb{G}}_{{n}_{2}^{\prime}{n}_{3}^{\prime}{n}_{1}^{\prime},{n}_{2}^{\u2033}{n}_{3}{n}_{1}^{\u2033}}^{(Z)}({t}_{3}-{t}^{\prime}){\mathbb{G}}_{{n}_{1}^{\u2033}{n}_{2}^{\u2033},{n}_{1}^{\u2034}{n}_{2}}^{(N)}({t}_{2}){G}_{{n}_{1}^{\u2034}{n}_{1}}({t}_{1}),$$

(116)

$${R}_{{\mathit{k}}_{\text{III}}}^{(\text{QP})}({t}_{3},{t}_{2},{t}_{1})=2i{\left(\frac{i}{\hslash}\right)}^{3}{\mathit{\mu}}_{{n}_{4}}^{\ast}{\mathit{\mu}}_{{n}_{3}}^{\ast}{\mathit{\mu}}_{{n}_{2}}{\mathit{\mu}}_{{n}_{1}}\times {\int}_{0}^{{t}_{3}}\text{d}{t}^{\prime}\phantom{\rule{0.16667em}{0ex}}{G}_{{n}_{4}{n}_{4}^{\prime}}({t}^{\prime}){V}_{{n}_{4}^{\prime}{n}_{3}^{\prime},{n}_{2}^{\prime}{n}_{1}^{\prime}}\times {\mathbb{G}}_{{n}_{3}^{\prime}{n}_{2}^{\prime}{n}_{1}^{\prime},{n}_{3}{n}_{2}^{\u2033}{n}_{1}^{\u2033}}^{(Z)}({t}_{3}-{t}^{\prime}){\mathbb{G}}_{{n}_{2}^{\u2033}{n}_{1}^{\u2033},{n}_{2}{n}_{1}^{\u2034}}({t}_{2}){G}_{{n}_{1}^{\u2034}{n}_{1}}({t}_{1}).$$

(117)

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 *t*_{2} interval represent (when both point up) and (when pointing in opposite directions). Three lines (triple propagation) within the *t*_{3} interval represent . A horizontal dashed line stands for the *V* matrix.

*R*_{kI} is depicted in the left diagram. The interval between the first and second interactions is described by the single-exciton variable ^{†}, whose propagation oscillates with interband frequencies according to *G*^{†}(*t*_{1}). During the second interval, the system is characterized by ^{†} and its propagation given by (*t*_{2}). This oscillates at low, intraband frequencies associated with differences between single-exciton energies on different sites and their couplings. The ^{†} variables generated during the third interval propagate according to (*t*_{3} – *t*′), which again oscillate with high interband frequencies. The evolution during *t*_{2} 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 (*t*_{2}). The high-frequency evolution during *t*_{1} and *t*_{3} 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 and : by neglecting exciton population relaxation in the third interval, we can set
${\mathbb{G}}_{\mathit{kmn},{k}^{\prime}{m}^{\prime}{n}^{\prime}}^{(Z)}\approx {G}_{k{k}^{\prime}}^{\ast}{\mathbb{G}}_{mn,{m}^{\prime}{n}^{\prime}}$.^{68} is expressed in terms of the exciton-scattering matrix Γ (Appendix E). The secular approximation, eq 31, will be used for so that exciton eigenstate populations will be decoupled from their coherences.

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

$${G}_{e}(t)\equiv \sum _{mn}{\psi}_{em}^{\ast}{\psi}_{en}{G}_{mn}(t)=\theta (t)exp(-i{\epsilon}_{e}t-{\gamma}_{e}t),$$

(118)

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

$${G}_{e}(\omega )\equiv \frac{1}{\omega -{\epsilon}_{e}+i{\gamma}_{e}}.$$

(119)

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

$${\mathrm{\Gamma}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}=\sum _{{n}_{4}{n}_{3}{n}_{2}{n}_{1}}{\psi}_{{e}_{4}{n}_{4}}^{\ast}{\psi}_{{e}_{3}{n}_{3}}^{\ast}{\psi}_{{e}_{2}{n}_{2}}{\psi}_{{e}_{1}{n}_{1}}{\mathrm{\Gamma}}_{{n}_{4}{n}_{3},{n}_{2}{n}_{1}}$$

(120)

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

$${\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(0)}\equiv {\mathbb{G}}_{{e}_{2}{e}_{1}}^{(0)}{\delta}_{{e}_{4}{e}_{2}}{\delta}_{{e}_{3}{e}_{1}},{\mathbb{G}}_{{e}_{2}{e}_{1}}^{(0)}(\mathrm{\Omega})\equiv \sum _{{n}_{4}{n}_{3}{n}_{2}{n}_{1}}{\psi}_{{e}_{2}{n}_{4}}^{\ast}{\psi}_{{e}_{1}{n}_{3}}^{\ast}{\psi}_{{e}_{2}{n}_{2}}{\psi}_{{e}_{1}{n}_{1}}{\mathbb{G}}_{{n}_{4}{n}_{3},{n}_{2}{n}_{1}}^{(0)}(\mathrm{\Omega})=\frac{1}{\mathrm{\Omega}-{\epsilon}_{{e}_{2}}-{\epsilon}_{{e}_{1}}+i{\gamma}_{{e}_{2}}+i{\gamma}_{{e}_{1}}}$$

(121)

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

$${\mathbb{G}}_{{e}^{\prime}e}^{(0)}(\omega )=i\int \frac{\text{d}{\omega}^{\prime}}{2\pi}{G}_{e}({\omega}^{\prime}){G}_{{e}^{\prime}}(\omega -{\omega}^{\prime}).$$

(122)

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

$${\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}(t)=\sum _{{n}_{4}{n}_{3}{n}_{2}{n}_{1}}{\psi}_{{e}_{4}{n}_{4}}^{\ast}{\psi}_{{e}_{3}{n}_{3}}{\psi}_{{e}_{2}{n}_{2}}{\psi}_{{e}_{1}{n}_{1}}^{\ast}{\mathbb{G}}_{{n}_{4}{n}_{3},{n}_{2}{n}_{1}}^{(N)}(t).$$

(123)

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 *k*_{I}, we obtain from eq 115

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=2{\left(\frac{i}{\hslash}\right)}^{3}\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}\sum _{{e}^{\prime}{e}^{\u2033}}{\mathit{\mu}}_{{e}_{4}}^{\ast}{\mathit{\mu}}_{{e}_{3}}{\mathit{\mu}}_{{e}_{2}}{\mathit{\mu}}_{{e}_{1}}^{\ast}\times {G}_{{e}_{4}}({\mathrm{\Omega}}_{3}){\mathrm{\Gamma}}_{{e}_{4}{e}^{\prime},{e}_{3}{e}^{\u2033}}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\prime}}+i{\gamma}_{{e}^{\prime}})\times {\mathbb{G}}_{{e}_{3}{e}^{\u2033}}^{(0)}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\prime}}+i{\gamma}_{{e}^{\prime}}){\mathbb{G}}_{{e}^{\u2033}{e}^{\prime},{e}_{2}{e}_{1}}^{(N)}({\mathrm{\Omega}}_{2}){G}_{{e}_{1}}^{\ast}(-{\mathrm{\Omega}}_{1}).$$

(124)

The time-domain response functions
${R}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({t}_{3},{t}_{2},{t}_{1})$ can be calculated by using the inverse transform (eq 317). We note that the Fourier transforms with respect to *t*_{1} and *t*_{2} only involve single Green’s functions and are trivial. We can thus derive simple expressions for mixed representation such as
${R}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})$. When incoherent exciton transport is neglected, we can set
${\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}(t)={\delta}_{{e}_{4}{e}_{2}}{\delta}_{{e}_{3}{e}_{1}}{G}_{{e}_{2}}(t){G}_{{e}_{1}}^{\ast}(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

$${R}_{{\mathit{k}}_{\text{II}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=2{\left(\frac{i}{\hslash}\right)}^{3}\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}\sum _{{e}^{\prime}{e}^{\u2033}}{\mathit{\mu}}_{{e}_{4}}^{\ast}{\mathit{\mu}}_{{e}_{3}}{\mathit{\mu}}_{{e}_{2}}^{\ast}{\mathit{\mu}}_{{e}_{1}}\times {G}_{{e}_{4}}({\mathrm{\Omega}}_{3}){\mathrm{\Gamma}}_{{e}_{4}{e}^{\u2033},{e}_{3}{e}^{\prime}}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\u2033}}+i{\gamma}_{{e}^{\u2033}})\times {\mathbb{G}}_{{e}_{3}{e}^{\prime}}^{(0)}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\u2033}}+i{\gamma}_{{e}^{\u2033}}){\mathbb{G}}_{{e}^{\prime}{e}^{\u2033},{e}_{1}{e}_{2}}^{(N)}({\mathrm{\Omega}}_{2}){G}_{{e}_{1}}({\mathrm{\Omega}}_{1})$$

(125)

and

$${R}_{{\mathit{k}}_{\text{III}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=2{\left(\frac{i}{\hslash}\right)}^{3}\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}{\mathit{\mu}}_{{e}_{4}}^{\ast}{\mathit{\mu}}_{{e}_{3}}^{\ast}{\mathit{\mu}}_{{e}_{2}}{\mathit{\mu}}_{{e}_{1}}\times [{\mathrm{\Gamma}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}({\mathrm{\Omega}}_{2}){\mathbb{G}}_{{e}_{2}{e}_{1}}^{(0)}({\mathrm{\Omega}}_{2})-{\mathrm{\Gamma}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}_{3}}+i{\gamma}_{{e}_{3}}){\mathbb{G}}_{{e}_{2}{e}_{1}}^{(0)}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}_{3}}+i{\gamma}_{{e}_{3}})]\times {G}_{{e}_{4}}({\mathrm{\Omega}}_{3}){G}_{{e}_{1}}({\mathrm{\Omega}}_{1}){G}_{{e}_{3}}^{\ast}({\mathrm{\Omega}}_{2}-{\mathrm{\Omega}}_{3}).$$

(126)

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

$${S}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-2{\hslash}^{-3}\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}\sum _{{e}^{\prime}{e}^{\u2033}}({\mathit{\mu}}_{{e}_{4}}^{\ast}\xb7{\mathbb{E}}_{4}^{\ast}({\epsilon}_{{e}_{4}}-{\omega}_{s}))({\mathit{\mu}}_{{e}_{3}}\xb7{\mathbb{E}}_{3}({\epsilon}_{{e}_{3}}-{\omega}_{3}))\times ({\mathit{\mu}}_{{e}_{2}}\xb7{\mathbb{E}}_{2}({\epsilon}_{{e}_{2}}-{\omega}_{2}))({\mathit{\mu}}_{{e}_{1}}^{\ast}\xb7{\mathbb{E}}_{1}^{\ast}({\epsilon}_{{e}_{1}}-{\omega}_{1}))\times {G}_{{e}_{4}}({\mathrm{\Omega}}_{3}){\mathrm{\Gamma}}_{{e}_{4}{e}^{\prime},{e}_{3}{e}^{\u2033}}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\prime}}+i{\gamma}_{{e}^{\prime}})\times {\mathbb{G}}_{{e}_{3}{e}^{\u2033}}^{(0)}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\prime}}+i{\gamma}_{{e}^{\prime}}){\mathbb{G}}_{{e}^{\u2033}{e}^{\prime},{e}_{2}{e}_{1}}^{(N)}({\mathrm{\Omega}}_{2}){G}_{{e}_{1}}^{\ast}(-{\mathrm{\Omega}}_{1}),$$

(127)

$${S}_{{\mathit{k}}_{\text{II}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-2{\hslash}^{-3}\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}\sum _{{e}^{\prime}{e}^{\u2033}}({\mathit{\mu}}_{{e}_{4}}^{\ast}\xb7{\mathbb{E}}_{4}^{\ast}({\epsilon}_{{e}_{4}}-{\omega}_{s}))({\mathit{\mu}}_{{e}_{3}}\xb7{\mathbb{E}}_{3}({\epsilon}_{{e}_{3}}-{\omega}_{3}))\times ({\mathit{\mu}}_{{e}_{2}}^{\ast}\xb7{\mathbb{E}}_{2}^{\ast}({\epsilon}_{{e}_{2}}-{\omega}_{2}))({\mathit{\mu}}_{{e}_{1}}\xb7{\mathbb{E}}_{1}({\epsilon}_{{e}_{1}}-{\omega}_{1}))\times {G}_{{e}_{4}}({\mathrm{\Omega}}_{3}){\mathrm{\Gamma}}_{{e}_{4}{e}^{\u2033},{e}_{3}{e}^{\prime}}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\u2033}}+i{\gamma}_{{e}^{\u2033}})\times {\mathbb{G}}_{{e}_{3}{e}^{\prime}}^{(0)}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\u2033}}+i{\gamma}_{{e}^{\u2033}}){\mathbb{G}}_{{e}^{\prime}{e}^{\u2033},{e}_{1}{e}_{2}}^{(N)}({\mathrm{\Omega}}_{2}){G}_{{e}_{1}}({\mathrm{\Omega}}_{1}),$$

(128)

and

$${S}_{{\mathit{k}}_{\text{III}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-2{\hslash}^{-3}\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}({\mathit{\mu}}_{{e}_{4}}^{\ast}\xb7{\mathbb{E}}_{4}^{\ast}({\epsilon}_{{e}_{4}}-{\omega}_{s}))({\mathit{\mu}}_{{e}_{3}}^{\ast}\xb7{\mathbb{E}}_{3}^{\ast}({\epsilon}_{{e}_{3}}-{\omega}_{3}))\times ({\mathit{\mu}}_{{e}_{2}}\xb7{\mathbb{E}}_{2}({\epsilon}_{{e}_{2}}-{\omega}_{2}))({\mathit{\mu}}_{{e}_{1}}\xb7{\mathbb{E}}_{1}({\epsilon}_{{e}_{1}}-{\omega}_{1}))\times [{\mathrm{\Gamma}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}({\mathrm{\Omega}}_{2}){\mathbb{G}}_{{e}_{2}{e}_{1}}^{(0)}({\mathrm{\Omega}}_{2})-{\mathrm{\Gamma}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}_{3}}+i{\gamma}_{{e}_{3}}){\mathbb{G}}_{{e}_{2}{e}_{1}}^{(0)}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}_{3}}+i{\gamma}_{{e}_{3}})]\times {G}_{{e}_{4}}({\mathrm{\Omega}}_{3}){G}_{{e}_{1}}({\mathrm{\Omega}}_{1}){G}_{{e}_{3}}^{\ast}({\mathrm{\Omega}}_{2}-{\mathrm{\Omega}}_{3})$$

(129)

These expressions contain fewer terms than their supermolecule counterparts (eqs 42–50) and allow one to make the approximations discussed above. The CED limit is obtained by setting
${\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}(\omega )={\delta}_{{e}_{4}{e}_{2}}{\delta}_{{e}_{3}{e}_{1}}{[\omega -({\epsilon}_{{e}_{2}}-{\epsilon}_{{e}_{1}})+i({\gamma}_{{e}_{2}}+{\gamma}_{{e}_{1}})]}^{-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 127–129 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 124–126 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.

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.^{162}^{–}^{164} 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 molecules^{70} and can directly probe excitation energy transfer timescales.

The FMO complex^{168}^{,}^{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,

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 Silbey^{180}^{,}^{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}^{,}^{183}^{–}^{185} 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 PSI^{37}^{,}^{38}^{,}^{188}^{,}^{191} and PSII.^{192}^{–}^{195} These include absorption, linear and circular dichroism, and pump–probe.

The peak linewidths for the *k*_{I} technique at *t*_{2} = 0 reveal the homogeneous and inhomonegeous dephasing rates.^{70}^{,}^{196} Variation of the cross-peak pattern with *t*_{2} allows one to monitor pigment interactions and excitation energy transfer timescales.^{40}^{,}^{81} During *t*_{2}, 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 studies^{199} 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.^{201}^{–}^{203}

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 = 0 and Δ → ∞ (two-level chromophores). *ε _{m}*,

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 (*W*_{b}, *W*_{d}, and *W*_{e}) used to calculate the **...**

We start by neglecting transport and using the quasiparticle approach (eqs 127–129) 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 *S*_{kI} signals at *t*_{2} = 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.

Im
${S}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})$ signal (imaginary part) in the MFA and CED for FMO at *t*_{2} = 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} +* ε

A similar comparison is shown in Figures 13 and and1414 for the double quantum coherence *S*_{kIII} 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.

$\mid {S}_{{\mathit{k}}_{\text{III}}}^{(\text{QP})}({t}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})\mid $ signal in the MFA and CED for FMO at *t*_{3} = 1 fs and *t*_{3} = 200 fs.

$\mid {S}_{{\mathit{k}}_{\text{III}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{t}_{1})\mid $ signal in the MFA and CED for FMO at *t*_{1} = 0 fs.

The remainder of this section shows simulated *k*_{I} signals that go beyond the CED to include transport during *t*_{2} at the secular Redfield theory level.

The 2D signals were calculated using eq 127 and include transport during *t*_{2}. We used Gaussian pulse envelopes,

$${\mathbb{E}}_{j}(\omega )={\mathbb{E}}_{j}{\text{e}}^{-{(\omega /{\sigma}_{j})}^{2}}.$$

(130)

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.

2D signal
$\mid {S}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})\mid $ 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 *W _{b}* in Figure 11 and cover the entire exciton band. At

In Figure 15B, the first two pulses are narrow (*W _{e}* in Figure 11;

In Figure 15C, the first two pulses are narrow (*W _{d}* in Figure 11) and selectively excite the

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
${R}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}$ response function, eq 124, we shall transform the dipoles back to real space using *μ** _{e}* = Σ

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=\sum _{{n}_{4}}{R}_{{\mathit{k}}_{\text{I}},{n}_{4}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1}),$$

(131)

with

$${R}_{{\mathit{k}}_{\text{I}},{n}_{4}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})=-2i{\mathit{\mu}}_{{n}_{4}}^{\ast}\sum _{{e}_{4}}{\psi}_{{e}_{4}{n}_{4}}\sum _{{e}_{3}{e}_{2}{e}_{1}}\sum _{{e}^{\prime}{e}^{\u2033}}{\mathit{\mu}}_{{e}_{3}}{\mathit{\mu}}_{{e}_{2}}{\mathit{\mu}}_{{e}_{1}}^{\ast}\times {G}_{{e}_{4}}({\mathrm{\Omega}}_{3}){\mathrm{\Gamma}}_{{e}_{4}{e}^{\prime},{e}_{3}{e}^{\u2033}}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\prime}}+i{\gamma}_{{e}^{\prime}})\times {\mathbb{G}}_{{e}_{3}{e}^{\u2033}}^{(0)}({\mathrm{\Omega}}_{3}+{\epsilon}_{{e}^{\prime}}+i{\gamma}_{{e}^{\prime}}){\mathbb{G}}_{{e}^{\u2033}{e}^{\prime},{e}_{2}{e}_{1}}^{(N)}({\mathrm{\Omega}}_{2}){G}_{{e}_{1}}^{\ast}(-{\mathrm{\Omega}}_{1}).$$

(132)

${R}_{{\mathit{k}}_{\text{I}},{n}_{4}}^{(\text{QP})}$ gives the signal generated at site *n*_{4}. By varying *t*_{2}, 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 *n*_{3}, *n*_{2}, and *n*_{1}, 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):

$${\kappa}_{a,{n}_{2}}(\omega )\propto \omega \sum _{e}{\psi}_{{n}_{2},e}{\mathit{\mu}}_{{n}_{2}}^{\ast}{\mathit{\mu}}_{e}\frac{{\gamma}_{eg}}{{(\omega -{\omega}_{eg})}^{2}+{\gamma}_{eg}^{2}}.$$

(133)

It gives the contribution of site *n*_{2}. 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.

The
${R}_{{\mathit{k}}_{\text{I}},{n}_{4}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2}=0,{\mathrm{\Omega}}_{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 *t*_{2} = 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 *R*_{kI,3}. The second lowest energy band, *b*, is still populated at *t*_{2} = 10 ps. Hence, *R*_{kI,4} and *R*_{kI,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*).

${R}_{{\mathit{k}}_{\text{I}},{n}_{4}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})$ signal (eq 132) for the FMO light-harvesting complex at *t*_{2} = 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 **...**

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}

$${L}_{e}=\frac{1}{{\displaystyle \sum _{n}}{\mid {\psi}_{en}\mid}^{4}},$$

(134)

*L _{e}* varies between 1 and for an -chromophore aggregate and is a measure of the effective number of chromophores contributing to the exciton

$${L}_{\rho}(t)=\frac{{({\displaystyle \sum _{nm}}\mid {\rho}_{mn}(t)\mid )}^{2}}{\mathbb{N}{\displaystyle \sum _{nm}}{\mid {\rho}_{mn}(t)\mid}^{2}},$$

(135)

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

$${L}_{d}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})=\sum _{{n}_{4}}\mid {R}_{{\mathit{k}}_{\text{I}},{n}_{4}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})\mid -\mid {R}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})\mid ,$$

(136)

where
${R}_{{\mathit{k}}_{\text{I}},{n}_{4}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})$ is the contribution of the *n*_{4}th chromophore to the response function. *L _{d}* will vanish for regions dominated by a single chromophore and will grow with the degree of delocalization. In Figure 19, we display

(Left column) Absolute value of the response function,
$\mid {R}_{{\mathit{k}}_{\text{I}}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})\mid $ (eq 124). (Middle column) The sum of the absolute value of the dissected response function,
${\sum}_{{n}_{4}}\mid {R}_{{\mathit{k}}_{\text{I}},{n}_{4}}^{(\text{QP})}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{}$ **...**

The quadratic excitonic couplings *J _{mn}* in eq 58 correspond to the “strong

Here we study the possible signatures of .^{210} The energy of two-exciton states for excitations on sites *m* and *n* (i.e., *E _{m}*

$${E}_{m,n}{\hslash}^{-1}={\epsilon}_{m}+{\epsilon}_{n}+{\mathbb{K}}_{mn}(1-{\delta}_{mn})+{\mathrm{\Delta}}_{m}{\delta}_{mn}.$$

(137)

The coupling implies that double excitations on neighboring sites have lower (negative ; exciton attraction) or higher (positive ; exciton repulsion) energy compared to the sum of energies of the two single excitations. In the NEE, enters solely through the scattering matrix. For finite , the coherent time-evolution of the *Y* variables will show resonances shifted by 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 coupling, which we use as a free parameter. Figure 20 depicts the *k*_{I} signals calculated with the NEE for = = ±300 cm^{−1} and *t*_{2} = 0.

Effects of the off-diagonal quartic coupling, , on the *k*_{I} spectra for *t*_{2} = 0. (Left column) The imaginary part of the response function; (Right column) absolute value. From top to bottom, = = 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,

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{SOS})}({\mathrm{\Omega}}_{3},{t}_{2}=0,{\mathrm{\Omega}}_{1})={R}_{{\mathit{k}}_{\text{I}}}^{(\text{GSB})}({\mathrm{\Omega}}_{3},{t}_{2}=0,{\mathrm{\Omega}}_{1})+{R}_{{\mathit{k}}_{\text{I}}}^{(\text{ESE})}({\mathrm{\Omega}}_{3},{t}_{2}=0,{\mathrm{\Omega}}_{1})+{R}_{{\mathit{k}}_{\text{I}}}^{(\text{ESA})}({\mathrm{\Omega}}_{3},{t}_{2}=0,{\mathrm{\Omega}}_{1}),$$

(138)

where

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{GSB})}({\mathrm{\Omega}}_{3},{t}_{2}=0,{\mathrm{\Omega}}_{1})={R}_{{\mathit{k}}_{\text{I}}}^{(\text{ESE})}({\mathrm{\Omega}}_{3},{t}_{2}=0,{\mathrm{\Omega}}_{1})=-{\left(\frac{i}{\hslash}\right)}^{3}\left[\frac{{\mid {\mathit{\mu}}_{+}\mid}^{4}}{({\mathrm{\Omega}}_{3}-{\epsilon}_{+}+i{\gamma}_{+})({\mathrm{\Omega}}_{1}+{\epsilon}_{+}+i{\gamma}_{+})}+\frac{{\mid {\mathit{\mu}}_{-}\mid}^{4}}{({\mathrm{\Omega}}_{3}-{\epsilon}_{-}+i{\gamma}_{-})({\mathrm{\Omega}}_{1}+{\epsilon}_{-}+i{\gamma}_{-})}+\frac{{\mid {\mathit{\mu}}_{-}\mid}^{2}{\mid {\mathit{\mu}}_{+}\mid}^{2}}{({\mathrm{\Omega}}_{3}-{\epsilon}_{+}+i{\gamma}_{+})({\mathrm{\Omega}}_{1}+{\epsilon}_{-}+i{\gamma}_{-})}+\frac{{\mid {\mathit{\mu}}_{-}\mid}^{2}{\mid {\mathit{\mu}}_{+}\mid}^{2}}{({\mathrm{\Omega}}_{3}-{\epsilon}_{-}+i{\gamma}_{-})({\mathrm{\Omega}}_{1}+{\epsilon}_{+}+i{\gamma}_{+})}\right],$$

(139)

$${R}_{{\mathit{k}}_{\text{I}}}^{(\text{ESA})}({\mathrm{\Omega}}_{3},{t}_{2}=0,{\mathrm{\Omega}}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\left[\frac{({\mathit{\mu}}_{+}^{2}{\mathit{\mu}}_{f,+}^{2}+{\mathit{\mu}}_{-}{\mathit{\mu}}_{f,-}{\mathit{\mu}}_{+}{\mathit{\mu}}_{f,+})}{({\mathrm{\Omega}}_{3}-{\epsilon}_{-}+\mathbb{K}+i{\gamma}_{f,+})({\mathrm{\Omega}}_{1}+{\epsilon}_{+}+i{\gamma}_{+})}+\frac{({\mathit{\mu}}_{-}^{2}{\mathit{\mu}}_{f,-}^{2}+{\mathit{\mu}}_{-}{\mathit{\mu}}_{f,-}{\mathit{\mu}}_{+}{\mathit{\mu}}_{f,+})}{({\mathrm{\Omega}}_{3}-{\epsilon}_{+}+\mathbb{K}+i{\gamma}_{f,-})({\mathrm{\Omega}}_{1}+{\epsilon}_{-}+i{\gamma}_{-})}\right],$$

(140)

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: ** μ_{±}** = Ψ

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 is increased, the ESA peaks are shifted along Ω_{3} by 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 = 0, the diagonal peaks (*b*, *b*) and (*d*, *d*) dominate. As 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* + ) and (*d*, *b* + ).

Signatures of the off-diagonal quartic couplings are also manifest for *t*_{2} > 0. In Figure 21, we show the *k*_{I} response function for different as in Figure 20 but for *t*_{2} = 10 ps. The *b* band is still populated even after long *t*_{2}. This is evident from the large diagonal peak in the *b* region for = 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 = −300 cm^{−1}, the doubly excited Φ_{f}_{45} state is energetically favored. The excitation energy, absorbed in the *b* or *d* band, which relaxes to *b* during *t*_{2}, can re-emit from the *d* band, leading to a high-energy diagonal peak and corresponding cross-peaks. Since the doubly excited state Φ_{f}_{45} is blue-shifted when = +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}.

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 additionally decouples coherences and populations. The response function may then be partitioned into contributions of coherences *C* and populations *P* during *t*_{2}.

The resonant contributions to each technique are given by the Feynman diagrams of Figures 4–7. 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 *k*_{I} response function will be recast in the form:

$${R}_{{\mathit{k}}_{\text{I}},\text{ia}}({t}_{3},{t}_{2},{t}_{1})=-i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime}}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{eg}{\mathit{\mu}}_{ge}\times {\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}({t}_{2})exp[-i{\xi}_{{e}^{\prime}g}{t}_{3}-i{\xi}_{ge}{t}_{1}+{\varphi}_{{e}^{\prime}e}^{(\text{ia})}({t}_{3},{t}_{2},{t}_{1})],$$

(141)

$${R}_{{\mathit{k}}_{\text{I}},\text{ib}}({t}_{3},{t}_{2},{t}_{1})=-i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime}}^{e\ne {e}^{\prime}}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{eg}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{ge}\times exp[-i{\xi}_{{e}^{\prime}g}{t}_{3}-i{\xi}_{{e}^{\prime}e}{t}_{2}-i{\xi}_{ge}{t}_{1}+{\varphi}_{{e}^{\prime}e}^{(\text{ib})}({t}_{3},{t}_{2},{t}_{1})],$$

(142)

$${R}_{{\mathit{k}}_{\text{I}},\text{ii}}({t}_{3},{t}_{2},{t}_{1})=-i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{e}^{\prime},e}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{eg}{\mathit{\mu}}_{ge}\times exp[-i{\xi}_{{e}^{\prime}g}{t}_{3}-i{\xi}_{ge}{t}_{1}+{\varphi}_{{e}^{\prime}e}^{(\text{ii})}({t}_{3},{t}_{2},{t}_{1})],$$

(143)

$${R}_{{\mathit{k}}_{\text{I}},\text{iiia}}({t}_{3},{t}_{2},{t}_{1})=+i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime},f}{\mathit{\mu}}_{{e}^{\prime}f}{\mathit{\mu}}_{f{e}^{\prime}}{\mathit{\mu}}_{eg}{\mathit{\mu}}_{ge}\times {\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}({t}_{2})exp[-i{\xi}_{f{e}^{\prime}}{t}_{3}-i{\xi}_{ge}{t}_{1}+{\varphi}_{f{e}^{\prime}e}^{(\text{iiia})}({t}_{3},{t}_{2},{t}_{1})],$$

(144)

$${R}_{{\mathit{k}}_{\text{I}},\text{iiib}}({t}_{3},{t}_{2},{t}_{1})=+i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime},f}^{e\ne {e}^{\prime}}{\mathit{\mu}}_{ef}{\mathit{\mu}}_{f{e}^{\prime}}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{ge}\times exp[-i{\xi}_{fe}{t}_{3}-i{\xi}_{{e}^{\prime}e}{t}_{2}-i{\xi}_{ge}{t}_{1}+{\varphi}_{f{e}^{\prime}e}^{(\text{iiib})}({t}_{3},{t}_{2},{t}_{1})].$$

(145)

As in section 4, we define the complex transition frequencies

$${\xi}_{ab}\equiv {\epsilon}_{a}-{\epsilon}_{b}-i({\tau}_{a}^{-1}+{\tau}_{b}^{-1})/2,$$

(146)

where *τ _{a}* is the lifetime of state

Similarly, we recast the *k*_{II} response in the form:

$${R}_{{\mathit{k}}_{\text{II}},\text{iva}}({t}_{3},{t}_{2},{t}_{1})=-i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime}}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{ge}{\mathit{\mu}}_{eg}\times {\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}({t}_{2})exp[-i{\xi}_{{e}^{\prime}g}{t}_{3}-i{\xi}_{eg}{t}_{1}+{\varphi}_{{e}^{\prime}e}^{(\text{iva})}({t}_{3},{t}_{2},{t}_{1})],$$

(147)

$${R}_{{\mathit{k}}_{\text{II}},\text{ivb}}({t}_{3},{t}_{2},{t}_{1})=-i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime}}^{e\ne {e}^{\prime}}{\mathit{\mu}}_{ge}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{eg}\times exp[-i{\xi}_{eg}{t}_{3}-i{\xi}_{{ee}^{\prime}}{t}_{2}-i{\xi}_{eg}{t}_{1}+{\varphi}_{{e}^{\prime}e}^{(\text{ivb})}({t}_{3},{t}_{2},{t}_{1})],$$

(148)

$${R}_{{\mathit{k}}_{\text{II}},\text{v}}({t}_{3},{t}_{2},{t}_{1})=-i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{{e}^{\prime},e}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{ge}{\mathit{\mu}}_{eg}\times exp[-i{\xi}_{{e}^{\prime}g}{t}_{3}-i{\xi}_{eg}{t}_{1}+{\varphi}_{{e}^{\prime}e}^{(\text{v})}({t}_{3},{t}_{2},{t}_{1})],$$

(149)

$${R}_{{\mathit{k}}_{\text{II}},\text{via}}({t}_{3},{t}_{2},{t}_{1})=+i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime}f}{\mathit{\mu}}_{{e}^{\prime}f}{\mathit{\mu}}_{f{e}^{\prime}}{\mathit{\mu}}_{ge}{\mathit{\mu}}_{eg}\times {\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}({t}_{2})exp[-i{\xi}_{f{e}^{\prime}}{t}_{3}-i{\xi}_{eg}{t}_{1}+{\varphi}_{f{e}^{\prime}e}^{(\text{via})}({t}_{3},{t}_{2},{t}_{1})],$$

(150)

$${R}_{{\mathit{k}}_{\text{II}},\text{vib}}({t}_{3},{t}_{2},{t}_{1})=+i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime},f}^{e\ne {e}^{\prime}}{\mathit{\mu}}_{{e}^{\prime}f}{\mathit{\mu}}_{fe}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{eg}\times exp[-i{\xi}_{f{e}^{\prime}}{t}_{3}-i{\xi}_{{ee}^{\prime}}{t}_{2}-i{\xi}_{eg}{t}_{1}+{\varphi}_{f{e}^{\prime}e}^{(\text{vib})}({t}_{3},{t}_{2},{t}_{1})].$$

(151)

Finally, *k*_{III} is composed of ESA_{1} and ESA_{2} (both are coherent contributions since populations are not created in *k*_{III}):

$${R}_{{\mathit{k}}_{\text{III}},\text{vii}}({t}_{3},{t}_{2},{t}_{1})=+i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime},f}{\mathit{\mu}}_{{e}^{\prime}f}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{fe}{\mathit{\mu}}_{eg}\times exp[-i{\xi}_{f{e}^{\prime}}{t}_{3}-i{\xi}_{fg}{t}_{2}-i{\xi}_{eg}{t}_{1}+{\varphi}_{f{e}^{\prime}e}^{(\text{vii})}({t}_{3},{t}_{2},{t}_{1})],$$

(152)

$${R}_{{\mathit{k}}_{\text{III}},\text{viii}}({t}_{3},{t}_{2},{t}_{1})=-i{\hslash}^{-3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\sum _{e,{e}^{\prime},f}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{{e}^{\prime}f}{\mathit{\mu}}_{fe}{\mathit{\mu}}_{eg}\times exp[-i{\xi}_{{e}^{\prime}g}{t}_{3}-i{\xi}_{fg}{t}_{2}-i{\xi}_{eg}{t}_{1}+{\varphi}_{f{e}^{\prime}e}^{(\text{viii})}({t}_{3},{t}_{2},{t}_{1})].$$

(153)

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 141–153 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.

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 *Ĥ = Σ _{υ}*|

$${\widehat{H}}_{v}={\hslash}^{-1}\langle v\mid \widehat{H}\mid v\rangle ={\epsilon}_{v}+{\widehat{q}}_{vv}+{\widehat{H}}_{B},$$

(154)

* _{υυ}* is a collective bath coordinate:

$${\widehat{q}}_{vv}\equiv {\hslash}^{-1}\langle v\mid {\widehat{H}}_{SB}\mid v\rangle =\sum _{\alpha}{d}_{vv,\alpha}({\widehat{a}}_{\alpha}^{\u2020}+{\widehat{a}}_{\alpha}).$$

(155)

* _{υυ}* induces the frequency fluctuations of state

$${\widehat{q}}_{vv}^{(a)}(t)\equiv {\text{e}}^{i{\hslash}^{-1}{\widehat{H}}_{a}t}{\widehat{q}}_{vv}{\text{e}}^{-i{\hslash}^{-1}{\widehat{H}}_{a}t}.$$

(156)

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

$${R}^{(1)}(t)=\frac{i}{\hslash}[\mathbb{J}(t)-{\mathbb{J}}^{\ast}(t)],$$

(157)

where (*τ*_{2} − *τ*_{1}) = ** (***τ*_{2})** (***τ*_{1}) is the two-point dipole correlation function. The angular brackets imply equilibrium statistical averaging over the bath. is calculated exactly for this model:

$$\mathbb{J}(t)=\sum _{ab}{\mathit{\mu}}_{ab}{\mathit{\mu}}_{ba}{\rho}_{a}exp[-i{\omega}_{ba}t-{g}_{bb}(t)-{g}_{aa}(t)+{g}_{ab}(t)+{g}_{ba}(t)].$$

(158)

Here, *ρ _{a}* is the equilibrium population of state

$${g}_{v{v}^{\prime}}(t)={\int}_{0}^{t}\text{d}{\tau}_{2}{\int}_{0}^{{\tau}_{2}}\text{d}{\tau}_{1}{C}_{v{v}^{\prime}}({\tau}_{2}-{\tau}_{1}).$$

(159)

${C}_{v{v}^{\prime}}({\tau}_{2}-{\tau}_{1})\equiv \langle {\widehat{q}}_{vv}^{(g)}({\tau}_{2}){\widehat{q}}_{{v}^{\prime}{v}^{\prime}}^{(g)}({\tau}_{1})\rangle $ 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,

$${C}_{v{v}^{\prime}}^{\u2033}(\omega )=\frac{1}{2}{\int}_{0}^{\infty}\text{d}texp(i\omega t)\langle [{\widehat{q}}_{vv}^{(g)}(t),{\widehat{q}}_{{v}^{\prime}{v}^{\prime}}^{(g)}(0)]\rangle ;$$

(160)

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

$${g}_{v{v}^{\prime}}(t)=\int \frac{\text{d}\omega}{2\pi}\frac{{C}_{v{v}^{\prime}}^{\u2033}(\omega )}{{\omega}^{2}}[coth\left(\frac{\beta \hslash \omega}{2}\right)(1-cos\omega t)+isin\omega t-i\omega t].$$

(161)

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

$${R}^{(3)}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}[F({t}_{1},{t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},0)+F(0,{t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1})+F(0,{t}_{1},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2})+F({t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},{t}_{1},0)]+c.c.,$$

(162)

where

$$F({\tau}_{4},{\tau}_{3},{\tau}_{2},{\tau}_{1})=\langle \widehat{\mathit{\mu}}({\tau}_{4})\widehat{\mathit{\mu}}({\tau}_{3})\widehat{\mathit{\mu}}({\tau}_{2})\widehat{\mathit{\mu}}({\tau}_{1})\rangle $$

(163)

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 *k** _{s}*). Expansion in the exciton eigenstates yields

$$F({\tau}_{4},{\tau}_{3},{\tau}_{2},{\tau}_{1})=\sum _{\mathit{dcba}}{\rho}_{a}{\mathit{\mu}}_{ad}{\mathit{\mu}}_{dc}{\mathit{\mu}}_{cb}{\mathit{\mu}}_{ba}exp[-i({\omega}_{da}{\tau}_{43}+{\omega}_{ca}{\tau}_{32}+{\omega}_{ba}{\tau}_{21})+{\phi}_{\mathit{dcba}}({\tau}_{4},{\tau}_{3},{\tau}_{2},{\tau}_{1})],$$

(164)

where * _{dcba}*(

$${\phi}_{\mathit{cbag}}({\tau}_{4},{\tau}_{3},{\tau}_{2},{\tau}_{1})=-{g}_{cc}({\tau}_{43})-{g}_{bb}({\tau}_{32})-{g}_{aa}({\tau}_{21})-{g}_{cb}({\tau}_{42})+{g}_{cb}({\tau}_{43})+{g}_{cb}({\tau}_{32})-{g}_{ca}({\tau}_{41})+{g}_{ca}({\tau}_{42})+{g}_{ca}({\tau}_{31})-{g}_{ca}({\tau}_{32})-{g}_{ba}({\tau}_{31})+{g}_{ba}({\tau}_{32})+{g}_{ba}({\tau}_{21}),$$

(165)

where *τ _{ij}* =

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 141–153 with
${\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}(t)=\theta (t){\delta}_{{e}^{\prime}e}$ and no lifetime broadening *τ*^{−1} = 0. The phase functions are given by

$$\begin{array}{l}{\varphi}_{ee}^{(\text{ia})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{\mathit{egeg}}^{\ast}({t}_{1},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},0),\\ {\varphi}_{{e}^{\prime}e}^{(\text{ii})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{geg}}^{\ast}({t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1},0),\\ {\varphi}_{\mathit{fee}}^{(\text{iiia})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{\mathit{efeg}}^{\ast}({t}_{1},{t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},0),\\ {\varphi}_{{e}^{\prime}e}^{(\text{ib})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{geg}}^{\ast}({t}_{1},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},0),\\ {\varphi}_{f{e}^{\prime}e}^{(\text{iiib})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{feg}}^{\ast}({t}_{1},{t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},0),\\ {\varphi}_{ee}^{(\text{iva})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{\mathit{egeg}}({t}_{1},{t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},0),\\ {\varphi}_{{e}^{\prime}e}^{(\text{v})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{geg}}({t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},{t}_{1},0),\\ {\varphi}_{\mathit{fee}}^{(\text{via})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{\mathit{efeg}}({t}_{1},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},0),\\ {\varphi}_{{e}^{\prime}e}^{(\text{ivb})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{geg}}({t}_{1},{t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},0),\\ {\varphi}_{f{e}^{\prime}e}^{(\text{vib})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{feg}}({t}_{1},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},0),\\ {\varphi}_{f{e}^{\prime}e}^{(\text{vii})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{feg}}({t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1},0),\\ {\varphi}_{f{e}^{\prime}e}^{(\text{viii})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{feg}}({t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},{t}_{1},0).\end{array}$$

(166)

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

$${C}_{v{v}^{\prime}}^{\u2033}(\omega )=2{\lambda}_{v{v}^{\prime}}\frac{\omega \mathrm{\Lambda}}{{\omega}^{2}+{\mathrm{\Lambda}}^{2}},$$

(167)

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),

$${g}_{v{v}^{\prime}}(t)=\frac{{\lambda}_{v{v}^{\prime}}}{\mathrm{\Lambda}}coth(\hslash \mathrm{\Lambda}\beta /2)[exp(-\mathrm{\Lambda}t)+\mathrm{\Lambda}t-1]+\frac{4{\lambda}_{vv\text{'}}\mathrm{\Lambda}}{\hslash \beta}\sum _{n=1}^{\infty}\frac{{\text{e}}^{-{\nu}_{n}t}+{\nu}_{n}t-1}{{\nu}_{n}({\nu}_{n}^{2}-{\mathrm{\Lambda}}^{2})}-i\frac{{\lambda}_{v{v}^{\prime}}}{\mathrm{\Lambda}}[exp(-\mathrm{\Lambda}t)+\mathrm{\Lambda}t-1],$$

(168)

where *β* = (*k*_{B}*T*)^{−1} and *ν _{n}* = (2

$${g}_{v{v}^{\prime}}(t)=\frac{{\lambda}_{v{v}^{\prime}}}{\mathrm{\Lambda}}\left(\frac{2{k}_{\text{B}}T}{\hslash \mathrm{\Lambda}}-i\right)(exp(-\mathrm{\Lambda}t)+\mathrm{\Lambda}t-1).$$

(169)

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 *t*_{1} and *t*_{3}, whereas transport takes place during *t*_{2}. Typically *t*_{1},*t*_{3} *t*_{2} (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 *t*_{2} only in the *k*_{I} and *k*_{II} techniques. We therefore focus on these techniques. We first recast the response function, eq 15, in the form

$${R}^{(3)}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}Tr[\mathbb{W}({t}_{3})\mathbb{G}({t}_{2})\mathbb{D}({t}_{1})],$$

(170)

where (*t*) (*t*) _{0} is the *doorway* exciton wavepacket prepared by the first two pulses, (*t*) describes its propagation and (*t*) |(*t*) is the *window* wavepacket, which represents the detection.

Both *k*_{I} and *k*_{II} 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}^{(3)}={R}^{(\text{ESE})}+{R}^{(\text{GSB})}+{R}^{(\text{ESA})}.$$

(171)

*R*^{(ESE)} is given by diagrams i and iv, *R*^{GSB} is given by ii and v, and *R*^{ESA} is given by iii and vi in Figures 4 and and55.

Expanding eq 170 in the system eigenstates, we get and

$${R}^{(\text{ESE})}={\left(\frac{i}{\hslash}\right)}^{3}\sum _{{e}_{4}{e}_{3}{e}_{2}{e}_{1}}\langle {\mathbb{W}}_{{e}_{4}g}^{\text{lr}}({t}_{3}){\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}({t}_{2})\times [{\mathbb{D}}_{g{e}_{1}}^{\text{lr}}({t}_{1})+{\mathbb{D}}_{{e}_{2}g}^{\text{rl}}({t}_{1})]\rangle ,$$

(172)

$${R}^{(\text{ESA})}={\left(\frac{i}{\hslash}\right)}^{3}\sum _{f{e}_{4}{e}_{3}{e}_{2}{e}_{1}}\langle {\mathbb{W}}_{f{e}_{3}}^{\text{ll}}({t}_{3}){\mathbb{G}}_{{e}_{4}{e}_{3},{e}_{2}{e}_{1}}^{(N)}({t}_{2})\times [{\mathbb{D}}_{g{e}_{1}}^{\text{lr}}({t}_{1})+{\mathbb{D}}_{{e}_{2}g}^{\text{rl}}({t}_{1})]\rangle ,$$

(173)

and

$${R}^{(\text{GSB})}={\left(\frac{i}{\hslash}\right)}^{3}\sum _{{e}_{2}{e}_{1}}\langle {\mathbb{W}}_{{e}_{2}g}^{\text{ll}}({t}_{3})[{\mathbb{D}}_{g{e}_{1}}^{\text{rr}}({t}_{1})+{\mathbb{D}}_{{e}_{1}g}^{\text{ll}}({t}_{1})]\rangle .$$

(174)

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 , we can separate the response function (eqs 172–174), into two terms:

$${R}^{(3)}({t}_{3},{t}_{2},{t}_{1})={R}_{\text{C}}^{(3)}({t}_{3},{t}_{2},{t}_{1})+{R}_{\text{P}}^{(3)}({t}_{3},{t}_{2},{t}_{1}).$$

(175)

The coherent contribution,
${R}_{\text{C}}^{(3)}$, only includes coherences (*e*_{1} ≠ *e*_{2}, *e*_{1} = *e*_{3}, *e*_{2} = *e*_{4}) during *t*_{2}, and the entire optical process is completed before a relaxed exciton population is created. The second term (
${R}_{\text{P}}^{(3)}$) only includes populations (*e*_{1} = *e*_{2}, *e*_{3} = *e*_{4}) and represents *sequential* exciton transport contributions. This representation thus keeps track on how populations are created and propagated during *t*_{2}. 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.
${R}_{c}^{(3)}$ may then be calculated using the formalism of section 11.1.
${R}_{p}^{(3)}$ is calculated using projection operator techniques^{185} by factorizing the bath averages in eqs 172–174 into , , and factors. This gives

$${R}_{\text{P}}^{(3)}({t}_{3}{t}_{2}{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}[{W}_{{e}^{\prime}}({t}_{3}){\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}({t}_{2}){D}_{e}({t}_{1})+{W}_{0}({t}_{3})+{D}_{0}({t}_{1})].$$

(176)

In the first, hopping, term, the doorway function *D*_{e} represents the population of the *e*th exciton created after two interactions with the radiation field. The hopping term includes the ESA diagrams (iiia for *k*_{I} and via for *k*_{II}) and ESE (ia for *k*_{I} and iva for *k*_{II}). The Green’s function
${\mathbb{G}}_{{e}^{\prime}{e}^{\prime},ee}^{(N)}$ is the conditional probability for the exciton to hop from *e* to *e′* during *t*_{2} obtained by solving the Pauli equation (eq 26 with the rate matrix in eq 362). The window function *W _{e}* represents the contribution of the

The final expressions for the *k*_{I} and *k*_{II} response functions are given by eqs 141–151. 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

$$\begin{array}{l}{\varphi}_{{e}^{\prime}e}^{(\text{ia})}({t}_{3},{t}_{2},{t}_{1})=-{g}_{ee}^{\ast}({t}_{1})-{g}_{{e}^{\prime}{e}^{\prime}}^{\ast}({t}_{3})+2i{\lambda}_{{e}^{\prime}{e}^{\prime}}{t}_{3},\\ {\varphi}_{ee}^{(\text{ii})}({t}_{3},{t}_{2},{t}_{1})=-{g}_{ee}^{\ast}({t}_{1})-{g}_{{e}^{\prime}{e}^{\prime}}({t}_{3}),\\ {\varphi}_{\mathit{fee}}^{(\text{iiia})}({t}_{3},{t}_{2},{t}_{1})=-{g}_{ee}^{\ast}({t}_{1})-{g}_{{e}^{\prime}{e}^{\prime}}^{\ast}({t}_{3})-{g}_{ff}({t}_{3})+{g}_{{e}^{\prime}f}({t}_{3})+2i({\lambda}_{ff}-{\lambda}_{{e}^{\prime}{e}^{\prime}}){t}_{3},\\ {\varphi}_{ee}^{(\text{iva})}({t}_{3},{t}_{2},{t}_{1})=-{g}_{ee}({t}_{1})-{g}_{{e}^{\prime}{e}^{\prime}}^{\ast}({t}_{3})+2i{\lambda}_{{e}^{\prime}{e}^{\prime}}{t}_{3},\\ {\varphi}_{{e}^{\prime}e}^{(\text{v})}({t}_{3},{t}_{2},{t}_{1})=-{g}_{ee}({t}_{1})-{g}_{{e}^{\prime}{e}^{\prime}}({t}_{3}),\\ {\varphi}_{f{e}^{\prime}e}^{(\text{via})}({t}_{3},{t}_{2},{t}_{1})=-{g}_{ee}({t}_{1})-{g}_{{e}^{\prime}e}^{\ast}({t}_{3})-{g}_{ff}({t}_{3})+{g}_{{e}^{\prime}f}({t}_{3})+2i({\lambda}_{ff}-{\lambda}_{{e}^{\prime}{e}^{\prime}}){t}_{3}.\end{array}$$

(177)

The lifetimes (eq 146) in this case are given by ${\tau}_{a}^{-1}={\mathrm{\sum}}_{b}{K}_{bb,aa}$, and

$${\lambda}_{v{v}^{\prime}}=Im\underset{t\to \infty}{lim}{\stackrel{.}{g}}_{v{v}^{\prime}}(t)$$

(178)

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 *t*_{2}, 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 *t*_{2} ignore bath-fluctuation correlations between different propagation intervals *t*_{1}, *t*_{2}, and *t*_{3}. During *t*_{1} and *t*_{3}, the diagonal fluctuations are treated exactly (by the cumulant expansion). Again, second-order Markovian perturbation theory is made for off-diagonal fluctuations. During *t*_{2}, 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 *t*_{2} period; correlations between fluctuations during intervals *t*_{1} and *t*_{3} have been neglected. In the next subsection, we present a higher-level approximation that does retain such correlations.

We now derive approximate expressions for the *k*_{I} and *k*_{II} techniques that include both slow bath fluctuations and transport.^{24} *k*_{III} does not involve transport and needs not be considered here. In subsection 11.2, we assumed fluctuation timescale shorter than the time intervals *t*_{1}, *t*_{2}, and *t*_{3}. Here, the bath fluctuations are slow compared to *t*_{1}, *t*_{2}, and *t*_{3}.

The Feynman diagrams for *k*_{I} (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 (*e* ≠ *e*′) during *t*_{2}. The population contributions, therefore, have two terms: a coherent, population-conserving term when population during *t*_{2} 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 *k*_{II}.

We assume a bath with both slow *Q*^{(S)} and fast *Q*^{(F)} modes. The Markovian approximation is applied to the fast modes. During *t*_{1} and *t*_{3}, these cause homogeneous line broadening (diagonal fluctuations), while during *t*_{2}, 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 *t*_{1}, *t*_{2}, and *t*_{3}, which cause correlations of density matrix between all these intervals. The diagonal parts,
${Q}_{ee}^{(\text{S})}$, shift the system energies, while the off-diagonal parts,
${Q}_{{ee}^{\prime}}^{(\text{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 *t*_{2}, 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 (*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 *k*_{I} and *k*_{II} techniques (eq 175) contain both coherent and sequential contributions. As in section 11.2, eqs 141–151 present the final expressions for *k*_{l} and *k*_{ll} 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

$$\begin{array}{l}{\varphi}_{{e}^{\prime}e}^{(\text{ia})}({t}_{3},{t}_{2},{t}_{1})={\delta}_{{e}^{\prime}e}{\phi}_{\mathit{egeg}}^{\ast}({t}_{1},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},0)+{\zeta}_{{e}^{\prime}e}{\overline{\phi}}_{{e}^{\prime}g{e}^{\prime}e}^{\ast}({t}_{3},{t}_{2},{t}_{1}),\\ {\varphi}_{{e}^{\prime}e}^{(\text{ii})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{geg}}^{\ast}({t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1},0),\\ {\varphi}_{{e}^{\prime}e}^{(\text{iiia})}({t}_{3},{t}_{2},{t}_{1})={\delta}_{{e}^{\prime}e}{\phi}_{\mathit{efeg}}^{\ast}({t}_{1},{t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},0)+{\zeta}_{{e}^{\prime}e}{\overline{\phi}}_{f{e}^{\prime}{e}^{\prime}e}^{\ast}({t}_{3},{t}_{2},{t}_{1}),\\ {\varphi}_{{e}^{\prime}e}^{(\text{iva})}({t}_{3},{t}_{2},{t}_{1})={\delta}_{{e}^{\prime}e}{\phi}_{\mathit{egeg}}({t}_{1},{t}_{1}+{t}_{2},{t}_{1}+{t}_{2}+{t}_{3},0)+{\zeta}_{{e}^{\prime}e}{\overline{\phi}}_{g{e}^{\prime}{e}^{\prime}e}({t}_{3},{t}_{2},{t}_{1}),\\ {\varphi}_{{e}^{\prime}e}^{(\text{v})}({t}_{3},{t}_{2},{t}_{1})={\phi}_{{e}^{\prime}\mathit{geg}}({t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},{t}_{1},0),\\ {\varphi}_{f{e}^{\prime}e}^{(\text{via})}({t}_{3},{t}_{2},{t}_{1})={\delta}_{{e}^{\prime}e}{\phi}_{\mathit{efeg}}({t}_{1},{t}_{1}+{t}_{2}+{t}_{3},{t}_{1}+{t}_{2},0)+{\zeta}_{{e}^{\prime}e}{\phi}_{{e}^{\prime}f{e}^{\prime}e}({t}_{3},{t}_{2},{t}_{1}),\end{array}$$

(179)

where we define ζ* _{ab}* = 1 −

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.

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.

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 over the ensemble of stochastic paths

$${R}^{(3)}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\times \langle Tr\left[\widehat{\mathit{\mu}}\mid {exp}_{+}\left[-\frac{i}{\hslash}{\int}_{{t}_{1}+{t}_{2}}^{{t}_{1}+{t}_{2}+{t}_{3}}{\mathbb{L}}_{0;\sigma ({t}^{\prime})}\text{d}{t}^{\prime}\right]\mathbb{V}\times {exp}_{+}\left[-\frac{i}{\hslash}{\int}_{{t}_{1}}^{{t}_{1}+{t}_{2}}{\mathbb{L}}_{0;\sigma ({t}^{\prime})}\text{d}{t}^{\prime}\right]\mathbb{V}{exp}_{+}\left[-\frac{i}{\hslash}{\int}_{0}^{{t}_{1}}{\mathbb{L}}_{0;\sigma ({t}^{\prime})}\text{d}{t}^{\prime}\right]\mathbb{V}{\rho}_{0}\right]\rangle .$$

(180)

Here, = [*Ĥ _{σ}*

Equation 180 may be recast in the form

$${R}^{(3)}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\theta ({t}_{1})\theta ({t}_{2})\theta ({t}_{3})\times \langle Tr[\widehat{\mathit{\mu}}\mid {\mathbb{G}}_{\sigma (t)}({t}_{3}+{t}_{2}+{t}_{1},{t}_{2}+{t}_{1})\times \mathbb{V}{\mathbb{G}}_{\sigma (t)}({t}_{2}+{t}_{1},{t}_{1})\mathbb{V}{\mathbb{G}}_{\sigma (t)}({t}_{1},0)\mathbb{V}{\rho}_{0}]\rangle ,$$

(181)

where (*t _{a}*,

$$\frac{\text{d}{\rho}_{\sigma (t)}}{\text{d}t}=-\frac{i}{\hslash}[{\widehat{H}}_{0;\sigma (t)},{\rho}_{\sigma (t)}].$$

(182)

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

$${exp}_{+}\left[-\frac{i}{\hslash}{\int}_{0}^{t}{\mathbb{L}}_{\sigma ({t}^{\prime})}\text{d}{t}^{\prime}\right]{\rho}_{\sigma}={exp}_{+}\left[-\frac{i}{\hslash}{\int}_{0}^{t}{\widehat{H}}_{\sigma ({t}^{\prime})}\text{d}{t}^{\prime}\right]{\rho}_{\sigma}{exp}_{-}\left[\frac{i}{\hslash}{\int}_{0}^{t}{\widehat{H}}_{\sigma ({t}^{\prime})}\text{d}{t}^{\prime}\right],$$

(183)

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

$$\frac{\text{d}P(\sigma )}{\text{d}t}={L}^{(\sigma )}P(\sigma ),$$

(184)

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:

$$\frac{\text{d}\rho}{\text{d}t}=-\frac{i}{\hslash}{\mathbb{L}}^{(\sigma )}\rho +{L}^{(\sigma )}\rho ;$$

(185)

*ρ* [*Ĥ _{σ}*,

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

$${R}^{(3)}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\langle 0\mid T{r}_{S}[\widehat{\mathit{\mu}}\mid \mathbb{G}({t}_{3})\mathbb{V}\mathbb{G}({t}_{2})\mathbb{V}\mathbb{G}({t}_{1})\mathbb{V}\mid \rho (0)\rangle \rangle ]$$

(186)

with

$$\mathbb{G}(t)\equiv \theta (t)exp\left(-\frac{i}{\hslash}{\mathbb{L}}_{0}^{(\sigma )}(t)+{L}^{(\sigma )}t\right).$$

(187)

Here, the initial joint density matrix |*ρ*(0) = |*gg*|0 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 0|*Tr _{S}* includes both tracing over the system and summation over the final bath states realized as a scalar product with the zero left eigenvector 0| of

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

$${R}^{(3)}({\mathrm{\Omega}}_{3},{\mathrm{\Omega}}_{2},{\mathrm{\Omega}}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\langle 0\mid T{r}_{S}[\widehat{\mathit{\mu}}\mid \mathbb{G}({\mathrm{\Omega}}_{3})\mathbb{V}\mathbb{G}({\mathrm{\Omega}}_{2})\mathbb{V}\mathbb{G}({\mathrm{\Omega}}_{1})\mathbb{V}\mid \rho (0)\rangle ],$$

(188)

$$\mathbb{G}(\mathrm{\Omega})=-{\left[i\mathrm{\Omega}-\frac{i}{\hslash}{\mathbb{L}}_{0}^{(\sigma )}\rho +{L}^{(\sigma )}\right]}^{-1}.$$

(189)

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:

$$\begin{array}{l}{R}_{\text{i}}({t}_{3},{t}_{2},{t}_{1})\equiv {\left(\frac{i}{\hslash}\right)}^{3}\times \sum _{\mathit{ijklmn}}\langle 0\mid {\mathit{\mu}}_{g{e}_{n}}{\mathbb{G}}_{{e}_{n}g,{e}_{m}g}({t}_{3}){\mathit{\mu}}_{{e}_{i}g}{\mathbb{G}}_{{e}_{m}{e}_{i},{e}_{k}{e}_{j}}({t}_{2}){\mathit{\mu}}_{{e}_{k}g}{\mathbb{G}}_{g{e}_{j},g{e}_{i}}({t}_{1}){\mathit{\mu}}_{g{e}_{i}}\mid 0\rangle ,\\ {R}_{\text{ii}}({t}_{3},{t}_{2},{t}_{1})\equiv {\left(\frac{i}{\hslash}\right)}^{3}\times \sum _{\mathit{ijkl}}\langle 0\mid {\mathit{\mu}}_{g{e}_{l}}{\mathbb{G}}_{{e}_{l}g,{e}_{k}g}({t}_{3}){\mathit{\mu}}_{g{e}_{k}}{\mathbb{G}}_{gg,gg}({t}_{2}){\mathit{\mu}}_{{e}_{j}g}{\mathbb{G}}_{g{e}_{j},g{e}_{i}}({t}_{1}){\mathit{\mu}}_{g{e}_{i}}\mid 0\rangle ,\\ {R}_{\text{iii}}({t}_{3},{t}_{2},{t}_{1})\equiv -{\left(\frac{i}{\hslash}\right)}^{3}\times \sum _{\mathit{ijklmnop}}\langle 0\mid {\mathit{\mu}}_{{f}_{p}{e}_{o}}{\mathbb{G}}_{{f}_{p}{e}_{o},{f}_{n}{e}_{m}}({t}_{3}){\mathit{\mu}}_{{f}_{n},{e}_{l}}{\mathbb{G}}_{{e}_{l}{e}_{m},{e}_{k}{e}_{j}}({t}_{2}){\mathit{\mu}}_{{e}_{k}g}{\mathbb{G}}_{g{e}_{j},g{e}_{i}}({t}_{1}){\mathit{\mu}}_{g{e}_{i}}\mid 0\rangle ,\\ {R}_{\text{iv}}({t}_{3},{t}_{2},{t}_{1})\equiv {\left(\frac{i}{\hslash}\right)}^{3}\times \sum _{\mathit{ijklmn}}\langle 0\mid {\mathit{\mu}}_{{e}_{n}g}{\mathbb{G}}_{{e}_{n}g,{e}_{i}g}({t}_{3}){\mathit{\mu}}_{{e}_{m}g}{\mathbb{G}}_{{e}_{l}{e}_{m},{e}_{j}{e}_{k}}({t}_{2}){\mathit{\mu}}_{g{e}_{k}}{\mathbb{G}}_{{e}_{j}g,{e}_{i}g}({t}_{1}){\mathit{\mu}}_{{e}_{i}g}\mid 0\rangle ,\\ {R}_{\text{v}}({t}_{3},{t}_{2},{t}_{1})\equiv {\left(\frac{i}{\hslash}\right)}^{3}\times \sum _{\mathit{ijkl}}\langle 0\mid {\mathit{\mu}}_{{e}_{i}g}{\mathbb{G}}_{{e}_{i}g,{e}_{k}g}({t}_{3}){\mathit{\mu}}_{{e}_{k}g}{\mathbb{G}}_{gg,gg}({t}_{2}){\mathit{\mu}}_{g{e}_{j}}{\mathbb{G}}_{{e}_{j}g,{e}_{i}g}({t}_{1}){\mathit{\mu}}_{{e}_{i}g}\mid 0\rangle ,\\ {R}_{\text{vi}}({t}_{3},{t}_{2},{t}_{1})\equiv -{\left(\frac{i}{\hslash}\right)}^{3}\times \sum _{\mathit{ijklmnop}}\langle 0\mid {\mathit{\mu}}_{{f}_{p}{e}_{o}}{\mathbb{G}}_{{f}_{p}{e}_{o},{f}_{n}{e}_{i}}({t}_{3}){\mathit{\mu}}_{{f}_{n}{e}_{m}}{\mathbb{G}}_{{e}_{m}{e}_{i},{e}_{j}{e}_{k}}({t}_{2}){\mathit{\mu}}_{g{e}_{k}}{\mathbb{G}}_{{e}_{j}g,{e}_{i}g}({t}_{1}){\mathit{\mu}}_{{e}_{i}g}\mid 0\rangle ,\\ {R}_{\text{vii}}({t}_{3},{t}_{2},{t}_{1})\equiv -{\left(\frac{i}{\hslash}\right)}^{3}\times \sum _{\mathit{ijklmnop}}\langle 0\mid {\mathit{\mu}}_{{f}_{p}{e}_{o}}{\mathbb{G}}_{{f}_{p}{e}_{o},{f}_{i}{e}_{m}}({t}_{3}){\mathit{\mu}}_{g{e}_{m}}{\mathbb{G}}_{{f}_{i}g,{f}_{k}g}({t}_{2}){\mathit{\mu}}_{{f}_{k}{e}_{j}}{\mathbb{G}}_{{e}_{j}g,{e}_{i}g}({t}_{1}){\mathit{\mu}}_{{e}_{i}g}\mid 0\rangle ,\\ {R}_{\text{viii}}({t}_{3},{t}_{2},{t}_{1})\equiv {\left(\frac{i}{\hslash}\right)}^{3}\times \sum _{\mathit{ijklmn}}\langle 0\mid {\mathit{\mu}}_{{e}_{n}g}{\mathbb{G}}_{{e}_{n}g,{e}_{m}g}({t}_{3}){\mathit{\mu}}_{{e}_{m}{f}_{i}}{\mathbb{G}}_{{f}_{i}g,{f}_{k}g}({t}_{2}){\mathit{\mu}}_{{f}_{k}{e}_{j}}{\mathbb{G}}_{{e}_{j}g,{e}_{i}g}({t}_{1}){\mathit{\mu}}_{{e}_{i}g}\mid 0\rangle .\end{array}$$

(190)

Equation 190 accounts for coherence transfer (e.g., *ge _{j}* →

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*_{α} *P*(*σ* = α) is

$$\frac{\text{d}{P}_{\alpha}}{\text{d}t}=\sum _{\beta =1}^{M}{L}_{\alpha \beta}^{(M)}{P}_{\beta}$$

(191)

and the SLE reads

$$\frac{\text{d}{\rho}_{ij,\alpha}}{\text{d}t}=\sum _{\beta =1}^{M}{L}_{\alpha \beta}^{(M)}{{\rho}_{ij,}}_{\beta}-\frac{i}{\hslash}\sum _{l}({[{H}_{\alpha}]}_{il}{\rho}_{lj,\alpha}-{[{H}_{\alpha}]}_{lj}{\rho}_{il,\alpha}),$$

(192)

The zero left eigenvector of , which represents the final summation over discrete bath, is given by 0| = (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 equation^{216}

$$\frac{\text{d}P(Q)}{\text{d}t}={L}^{(Q)}P(Q)=\mathrm{\Lambda}\frac{\partial}{\partial Q}\left(Q+\frac{\partial}{\partial Q}\right)P(Q),$$

(193)

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) = *δ*(*Q* − *Q′*) is

$$G(Q,{Q}^{\prime};t)=\sqrt{\frac{1}{2\pi (1-{\text{e}}^{-2\mathrm{\Lambda}t})}}exp\left[\frac{-{(Q-{\text{e}}^{-\mathrm{\Lambda}t}{Q}^{\prime})}^{2}}{2(1-{\text{e}}^{-2\mathrm{\Lambda}t})}\right].$$

(194)

The bath density approaches the equilibrium distribution *P _{eq}*(

$$\langle Q(t)Q(0)\rangle =\int \int Q{Q}^{\prime}G(Q,{Q}^{\prime};t){P}_{eq}(Q)\text{d}Q\phantom{\rule{0.16667em}{0ex}}\text{d}{Q}^{\prime}={\text{e}}^{-\mathrm{\Lambda}t}.$$

(195)

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

$$\langle Q(t)Q({t}_{2})\cdots Q({t}_{n})\rangle =\langle Q({t}_{1})Q({t}_{2})\rangle \langle Q({t}_{3})\cdots Q({t}_{n})\rangle +\langle Q({t}_{1})Q({t}_{3})\rangle \langle Q({t}_{2})\cdots Q({t}_{n})\rangle +\cdots +\langle Q({t}_{1})Q({t}_{n})\rangle \langle Q({t}_{2})\cdots Q({t}_{n-1})\rangle .$$

(196)

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

$${\mathbb{L}}_{0}^{(Q)}={\mathbb{L}}_{0}+Q{\mathbb{L}}_{1},$$

(197)

where and operate on the system Liouville space. We further assume noninteracting chromophores with purely diagonal bath fluctuations

$${\widehat{H}}_{0}+{\widehat{H}}_{SB}=\sum _{i}\hslash ({\epsilon}_{i}+{d}_{ii}Q){B}_{i}^{\u2020}{B}_{i}$$

(198)

The solution of the stochastic Liouville equations without the field

$$\frac{\text{d}{\rho}_{ij}}{\text{d}t}=-i[({\epsilon}_{i}-{\epsilon}_{j})+({d}_{ii}-{d}_{jj})Q]{\rho}_{ij}$$

(199)

is

$$\langle {\rho}_{ij}(t)\rangle ={\text{e}}^{-i{\omega}_{ij}t}\langle exp[-i{\int}_{0}^{t}({d}_{ii}-{d}_{jj})Q({t}^{\prime})\text{d}{t}^{\prime}]{\rho}_{ij}(0)\rangle .$$

(200)

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.

$$\begin{array}{l}S[\phi ]\equiv \langle exp[i\int \phi (t)Q(t)\phantom{\rule{0.16667em}{0ex}}\text{d}t]\rangle \\ =\sum _{n}\frac{{(-i)}^{n}}{n!}\int \cdots \int \phi ({t}_{1})\cdots \phi ({t}_{n})\langle Q({t}_{1})\cdots Q({t}_{n})\rangle \\ =\sum _{n}\frac{{(-i)}^{n}(n-1)!!}{n!}\times \int \int \phi ({t}^{\prime})\phi ({t}^{\u2033})\langle Q({t}^{\prime})Q({t}^{\u2033})\rangle \text{d}{t}^{\prime}\text{d}{t}^{\u2033}\times \cdots \times \int \int \phi ({t}^{(n-1)})\phi ({t}^{(n)})\langle Q({t}^{(n-1)})Q({t}^{(n)})\rangle \text{d}{t}^{(n-1)}\text{d}{t}^{(n)}\\ =exp[\frac{-1}{2}\int \int \phi ({t}^{\prime})\phi ({t}^{\u2033})\langle Q({t}^{\prime})Q({t}^{\u2033})\rangle \text{d}{t}^{\prime}\text{d}{t}^{\u2033}].\end{array}$$

(201)

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

$${\rho}_{ij}(t)={\text{e}}^{-i{\omega}_{ij}t-{g}_{ii}(t)-{g}_{jj}(t)+{g}_{ij}(t)+{g}_{ji}(t)}{\rho}_{ij}(0)$$

(202)

with the line-shape function

$${g}_{ij}(t)={\int}_{0}^{t}{\int}_{0}^{{t}^{\prime}}{d}_{ii}{d}_{jj}\langle Q({t}^{\prime})Q({t}^{\u2033})\rangle \text{d}{t}^{\u2033}\text{d}{t}^{\prime}=\frac{{d}_{ii}{d}_{jj}}{{\mathrm{\Lambda}}^{2}}({\text{e}}^{-\mathrm{\Lambda}t}+\mathrm{\Lambda}t-1).$$

(203)

Equation 203 may be also obtained from eq 161 by taking the infinite temperature limit, and the following spectral density
${C}_{ij}^{\u2033}(\omega )=(\beta \hslash ){d}_{ii}{d}_{jj}\omega \mathrm{\Lambda}/({\omega}^{2}+{\mathrm{\Lambda}}^{2})$. Comparing eqs 169 and 203, we find 2λ* _{νν′}* (

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

$${R}_{i}^{(3)}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\sum _{{ee}^{\prime}}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{eg}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{ge}{\text{e}}^{i{\epsilon}_{e}({t}_{1}+{t}_{2})-i{\epsilon}_{{e}^{\prime}}({t}_{2}+{t}_{3})}\times \langle exp(i{\int}_{0}^{{t}_{1}+{t}_{2}}{d}_{ee}Q({t}^{\prime})\text{d}{t}^{\prime}-i{\int}_{{t}_{1}}^{{t}_{1}+{t}_{2}+{t}_{3}}{d}_{{e}^{\prime}{e}^{\prime}}Q({t}^{\prime})\text{d}{t}^{\prime})\rangle .$$

(204)

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

$${R}_{i}^{(3)}({t}_{3},{t}_{2},{t}_{1})={\left(\frac{i}{\hslash}\right)}^{3}\sum _{{ee}^{\prime}}{\mathit{\mu}}_{g{e}^{\prime}}{\mathit{\mu}}_{eg}{\mathit{\mu}}_{{e}^{\prime}g}{\mathit{\mu}}_{ge}{\text{e}}^{i{\epsilon}_{e}({t}_{1}+{t}_{2})-i{\epsilon}_{{e}^{\prime}}({t}_{2}+{t}_{3})}\times exp[-{g}_{ee}({t}_{1}+{t}_{2})-{g}_{{e}^{\prime}{e}^{\prime}}({t}_{2}+{t}_{3})+{g}_{{ee}^{\prime}}({t}_{1}+{t}_{2}+{t}_{3})+{g}_{{ee}^{\prime}}({t}_{2})-{g}_{{ee}^{\prime}}({t}_{1})-{g}_{{ee}^{\prime}}({t}_{3})].$$

(205)

This agrees with eqs 162–165. 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),

$${\mid \alpha \rangle}_{Q}=\frac{exp[-({Q}^{2}/2)]}{{2}^{\alpha}\sqrt{2\pi}\alpha !}{\mathbb{H}}_{\alpha}\left(\frac{Q}{\sqrt{2}}\right),$$

(206)

where are Hermite polynomials and α = 0, 1, 2, · · · . The initial (equilibrium) bath state is represented by |0 = (1, 0, 0, ..., 0). The final summation assumes 0| = (1, 0, 0, ..., 0).

The bath Liouvillian is thus diagonal [*L*^{(}^{Q}^{)}]_{αα}*′*= −αΛ *δ*_{α α}* _{′}*, and the

$${[Q]}_{\alpha {\alpha}^{\prime}}={\alpha}^{\prime}\sqrt{2}{\delta}_{\alpha ,{\alpha}^{\prime}+1}+\frac{1}{\sqrt{2}}{\delta}_{\alpha ,{\alpha}^{\prime}-1}.$$

(207)

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

$$\frac{\text{d}{\rho}_{ij,\alpha}}{\text{d}t}=-\alpha \mathrm{\Lambda}{\rho}_{ij,\alpha}-\frac{i}{\hslash}\sum _{lm}{[{\mathbb{L}}_{0}]}_{ij,lm}{\rho}_{lm,\alpha}-\frac{i}{\hslash}\sum _{lm}{[{\mathbb{L}}_{1}]}_{ij,lm}\times \left[\alpha \sqrt{2}{\rho}_{lm,\alpha -1}+\frac{1}{\sqrt{2}}{\rho}_{lm,\alpha +1}\right].$$

(208)

By switching to Hilbert space notation, we have

$$\frac{\text{d}{\rho}_{ij,\alpha}}{\text{d}t}=-k\mathrm{\Lambda}{\rho}_{ij,\alpha}-\frac{i}{\hslash}\sum _{l}({[{\widehat{H}}_{0}]}_{il}{\rho}_{lj,\alpha}-{[{\widehat{H}}_{0}]}_{lj}{\rho}_{il,\alpha})-\frac{i}{\hslash}\sum _{l}{[{\widehat{H}}_{1}]}_{il}\left[\alpha \sqrt{2}{\rho}_{lj,\alpha -1}+\frac{1}{\sqrt{2}}{\rho}_{lj,\alpha +1}\right]-{[{\widehat{H}}_{1}]}_{lj}\left[\alpha \sqrt{2}{\rho}_{il,\alpha -1}+\frac{1}{\sqrt{2}}{\rho}_{il,\alpha +1}\right].$$

(209)

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
${\mathbb{S}}_{\alpha}^{+}$ (lowering
${\mathbb{S}}_{\alpha}^{-}$) operators
${\mathbb{G}}_{\alpha \pm 1,\beta}={\mathbb{S}}_{\alpha}^{\pm}{\mathbb{G}}_{\alpha ,\beta}$, these can be calculated iteratively in the form of a continued fraction

$${\mathbb{S}}_{\alpha}^{+}=-{\hslash}^{-1}\sqrt{2}(\alpha +1){\left[i\omega +{\hslash}^{-1}{\mathbb{L}}_{0}-(\alpha +1)\mathrm{\Lambda}+\frac{1}{\sqrt{2}}{\hslash}^{-1}{\mathbb{L}}_{1}{\mathbb{S}}_{\alpha +1}^{-}\right]}^{-1}{\mathbb{L}}_{1}.$$

(210)

and

$${\mathbb{S}}_{\alpha}^{-}=-{\hslash}^{-1}\frac{1}{\sqrt{2}}{[i\omega +{\hslash}^{-1}{\mathbb{L}}_{0}-(\alpha -1)\mathrm{\Lambda}+\sqrt{2}(\alpha -1){\hslash}^{-1}{\mathbb{L}}_{1}{\mathbb{S}}_{\alpha -1}^{-}]}^{-1}{\mathbb{L}}_{1}.$$

(211)

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

$${\mathbb{G}}_{\alpha ,\alpha}=-{\left[i\omega +{\hslash}^{-1}{\mathbb{L}}_{0}-\alpha \mathrm{\Lambda}+\sqrt{2}\alpha {\hslash}^{-1}{\mathbb{L}}_{1}{\mathbb{S}}_{\alpha}^{-}+\frac{1}{\sqrt{2}}{\hslash}^{-1}{\mathbb{L}}_{1}{\mathbb{S}}_{\alpha}^{+}\right]}^{-1}.$$

(212)

For fast fluctuations || Λ, 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

$$\frac{\text{d}\rho}{\text{d}t}=[-i{\hslash}^{-1}{\mathbb{L}}_{0}-{\mathrm{\Lambda}}^{-1}{\hslash}^{-2}{\mathbb{L}}_{1}{\mathbb{L}}_{1}]\rho .$$

(213)

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 *C _{ij}*(

$${C}_{ij}^{\u2033}(\omega )=\sum _{k}\beta \hslash {d}_{ii}^{(k)}{d}_{jj}^{(k)}\frac{\omega {\mathrm{\Lambda}}_{k}}{{\omega}^{2}+{\mathrm{\Lambda}}_{k}^{2}}.$$

(214)

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 *q _{ij}*(

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 ( block). Assuming fast stochastic dynamics, we can employ eq 213. For a general coupling of the form *Ĥ*_{1} = Σ* _{lm} d_{lm}*|

$${\hslash}^{-2}{[{\mathbb{L}}_{1}{\mathbb{L}}_{1}]}_{ij,kl}=-2{d}_{ik}{d}_{ij}+{\delta}_{jl}\sum _{m}{d}_{im}{d}_{mk}+{\delta}_{ik}\sum _{m}{d}_{lm}{d}_{mj}.$$

(215)

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

$${K}_{ii,jj}=-2{\mathrm{\Lambda}}^{-1}{d}_{ij}{d}_{ji},$$

(216)

$${K}_{ii,ii}=2{\mathrm{\Lambda}}^{-1}\sum _{j\ne i}{d}_{ij}{d}_{ji}.$$

(217)

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

$${K}_{ij,ij}={\mathrm{\Lambda}}^{-1}\left[{({d}_{ii}-{d}_{jj})}^{2}+\sum _{k\ne i}{d}_{ik}{d}_{kj}+\sum _{k\ne j}{d}_{jk}{d}_{kj}\right].$$

(218)

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 model^{217}^{–}^{219} (eq 192 with *M* = 2) has been applied^{220} 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 *S*_{kI} 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 *t*_{2} = 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 *t*_{1} and *t*_{3} without any exciton transport during *t*_{2}. 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 *t*_{1} and *t*_{3} interval must be the same and no cross-peaks are observed.

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 *t*_{2} = 0, *ε*_{1}_{u} = 30Γ, *ε*_{1}_{d} = 10Γ, *ε*_{2}_{u} = −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}*(

$${R}_{A}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})\equiv -Im[{R}_{{\mathit{k}}_{\text{I}}}({\mathrm{\Omega}}_{3},{t}_{2},-{\mathrm{\Omega}}_{1})+{R}_{{\mathit{k}}_{\text{II}}}({\mathrm{\Omega}}_{3},{t}_{2},{\mathrm{\Omega}}_{1})]$$

(219)

is plotted for two *t*_{2} 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.

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 (,) coupled to many harmonic modes and described by the Hamiltonian^{183}^{,}^{224}^{,}^{225}

$${\widehat{H}}_{B}=\frac{{\widehat{P}}^{2}}{2M}+\frac{M{\mathrm{\Omega}}^{2}{\widehat{Q}}^{2}}{2}+\sum _{j}\frac{{\widehat{p}}_{\text{j}}^{2}}{2{m}_{j}}+\frac{{m}_{j}{\omega}_{j}^{2}}{2}{\left({\widehat{q}}_{j}-\frac{{c}_{j}\widehat{Q}}{{m}_{j}{\omega}_{j}^{2}}\right)}^{2},$$

(220)

$${\widehat{H}}_{SB}=\sum _{mn}M{\mathrm{\Omega}}^{2}{\stackrel{\sim}{d}}_{mn}{B}_{m}^{+}{B}_{n}\widehat{Q}.$$

(221)

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.^{226}^{–}^{229}

We introduce the bath spectral density, which describes the correlations of the force = Σ* _{j}c_{j}_{j},* acting on the coordinates by the other bath modes.

$${J}^{\u2033}(\omega )=\frac{1}{2}\int \text{d}t\phantom{\rule{0.16667em}{0ex}}{\text{e}}^{i\omega t}\langle [\widehat{F}(t),\widehat{F}(0)]\rangle .$$

(222)

Comparing with eq 352, we get

$${J}^{\u2033}(\omega )\equiv \sum _{j}\frac{\hslash \pi {c}_{j}^{2}}{2{m}_{j}{\omega}_{j}}[\delta (\omega -{\omega}_{j})-\delta (\omega +{\omega}_{j})].$$

(223)

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

$${C}^{(Q)}(\omega )=\frac{1}{2}\int \text{d}t\phantom{\rule{0.16667em}{0ex}}{\text{e}}^{i\omega t}\langle [\widehat{Q}(t),\widehat{Q}(0)]\rangle =\frac{1}{{M}^{2}}\frac{{J}^{\u2033}(\omega )}{{({\mathrm{\Omega}}^{2}-{\omega}^{2})}^{2}+{[{M}^{-1}{\hslash}^{-1}{J}^{\u2033}(\omega )]}^{2}}.$$

(224)

Assuming that the microscopic bath modes are faster than the collective coordinate , taking *J″* (*ω*) = *M*Ω^{2}*ω/*Λ, and approximating Ω^{2} − *ω*^{2} ≈ Ω^{2}, the bath spectral density reduces to the overdamped Brownian oscillator form:

$${C}^{(Q)}(\omega )=\frac{\hslash}{M{\mathrm{\Omega}}^{2}}\frac{\omega \mathrm{\Lambda}}{{\omega}^{2}+{\mathrm{\Lambda}}^{2}}.$$

(225)

At high temperatures coth(*βω*/2) ≈ 2(*βω*)^{−1}, we recover the Gaussian–Markovian stochastic process (eq 195) with (½){ (*t*), (0)} = (*β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}*Ω* → *Q*, to obtain a dimensionless unit equilibrium width, and setting *M*Ω^{2}* _{mn}* →

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

$$T=-\frac{\mathrm{\Lambda}\beta}{2}\left(\frac{\partial}{\partial Q}\right)\{{H}_{1},\dots \}.$$

(226)

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

$${\left[\frac{\partial}{\partial Q}\right]}_{{\alpha}^{\prime}\alpha}=\sqrt{2\alpha}{\delta}_{{\alpha}^{\prime},\alpha +1}.$$

(227)

The resulting equations of motion

$$\frac{\text{d}\rho}{\text{d}t}=-\frac{i}{\hslash}{\mathbb{L}}^{(\sigma )}\rho +{L}^{(\sigma )}\rho +T\rho $$

(228)

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 *d _{mn} = δ_{mn}d_{mm}* (no off-diagonal coupling fluctuations). Expanding

$${[L+T]}_{\mid {e}_{m}{e}_{m}\rangle \rangle}^{(Q)}=\mathrm{\Lambda}\frac{\partial}{\partial Q}\left(Q-\beta {d}_{mn}+\frac{\partial}{\partial Q}\right).$$

(229)

For this model, the free-energy surfaces for the various states are linearly displaced by *β*(*d _{mm}* −

$${[L+T]}_{\mid {e}_{m}{e}_{n}\rangle \rangle}^{(Q)}=\mathrm{\Lambda}\frac{\partial}{\partial Q}\left(Q-\beta \frac{{d}_{mm}+{d}_{nn}}{2}\frac{\partial}{\partial Q}\right)$$

(230)

and the particle moves on a surface whose equilibrium displacement lies midway between those of levels |*e _{m}* and |

Equation 228 holds at moderately high temperatures. At low temperatures, *β*Ω > 1, the quantization of bath oscillators is necessary and higher-order terms in *β* obtained by expanding the coth (*βω*/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*λ _{νν′}* (

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 [(*t*)]* _{lk}* for the random walk jump of the bath from state

$$\mathbb{Z}(\tau )={\int}_{0}^{\tau}\mathbb{Y}(\tau -{\tau}^{\prime})exp\left[-\frac{i}{\hslash}\mathbb{L}(\tau -{\tau}^{\prime})\right]\mathbb{Z}({\tau}^{\prime})\phantom{\rule{0.16667em}{0ex}}\text{d}{\tau}^{\prime}+\delta (\tau )$$

(231)

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
${t}_{m}^{\prime}$ in the *m*th interval and the first jump at
${t}_{l}^{\prime}$ in some subsequent (*l*th) interval is determined by the fixed state of bath and contributes by factor

$$\mathbb{Y}({t}_{l}^{\prime}+{t}_{m}^{\prime}+\sum _{i=m+1}^{l-1}{t}_{i})exp\left(-\frac{i}{\hslash}\mathbb{L}{t}_{l}^{\prime}\right)\mathbb{V}\times \prod _{i=m+1}^{l-1}\left[exp\left(-\frac{i}{\hslash}\mathbb{L}{t}_{i}\right)\mathbb{V}\right]exp\left(-\frac{i}{\hslash}\mathbb{L}{t}_{m}^{\prime}\right).$$

(232)

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 . All of these contributions must be properly multiplied, averaged over initial and summed over final states, and convoluted over the possible ${t}_{l}^{\prime}$ and ${t}_{m}^{\prime}$.

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 (*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
$\sim 1/{t}_{2}^{\alpha -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 *x*^{2} *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 *t*_{0} elapsed between the start of the random walk and the first interaction with laser pulse. The parameters correspond to fast bath limit at *t*_{0} = 0, with a single motionally narrowed peak. Diagonal peaks at the fundamental frequencies that grow with *t*_{0} 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.

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 (space inversion (*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 (*F _{o} =* −

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 visible^{241}^{,}^{242} and the UV^{243}^{–}^{245} or vibrational transitions in the IR.^{246}^{–}^{248} 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 techniques^{251} 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.^{254}^{–}^{256} 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}^{,}^{90}^{–}^{92}^{,}^{94}^{,}^{97}^{,}^{156}^{,}^{259}^{–}^{263}

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 (*p* • *A*) Hamiltonian,

$${\widehat{H}}_{pA}=-\frac{1}{c}\int \text{d}\mathit{r}[\widehat{\mathit{J}}(\mathit{r})\xb7\mathit{A}(\mathit{r})+\widehat{\sigma}(\mathit{r})\mathit{A}(\mathit{r})\xb7\mathit{A}(\mathit{r})],$$

(233)

where ** A** is the vector potential of the electromagnetic field.

$$\widehat{\mathit{J}}=\sum _{\alpha}\delta (\mathit{r}-{\widehat{\mathit{r}}}_{\alpha})\frac{{q}_{\alpha}}{{m}_{\alpha}}{\widehat{\mathit{p}}}_{\alpha}$$

(234)

is the electric current operator, and

$$\widehat{\sigma}(\mathit{r})=-\sum _{\alpha}\delta (\mathit{r}-{\widehat{\mathit{r}}}_{\alpha})\frac{{q}_{\alpha}^{2}}{2{m}_{\alpha}c}$$

(235)

is the charge density (*r*_{α} and *p*_{α} are the coordinate and momentum operators of electron α). The *p* • *A* 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:

$${\widehat{H}}^{\prime}=-\frac{1}{c}\int \text{d}\mathit{r}\widehat{\mathit{J}}(\mathit{r})\mathit{A}(\mathit{r},t).$$

(236)

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

$$\mathit{A}(\mathit{r},t)=\sum _{j}\sum _{{u}_{j}=\pm 1}{\mathbb{A}}_{j}^{({u}_{j})}(t-{\tau}_{j}){\text{e}}^{{iu}_{j}{\mathit{k}}_{j}\mathit{r}-{iu}_{j}{\omega}_{j}(t-{\tau}_{j})},$$

(237)

where
${\mathbb{A}}_{j}^{({u}_{j})}$ is the vector potential envelope corresponding to pulse *j*, centered at *τ _{j}*, with wavevector

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

$${\mathit{J}}^{(3)}(\mathit{r},t)=\int \int \int \text{d}{\mathit{r}}_{3}\text{d}{\mathit{r}}_{2}\text{d}{\mathit{r}}_{1}{\int}_{0}^{\infty}\text{d}{t}_{3}{\int}_{0}^{\infty}\text{d}{t}_{2}{\int}_{0}^{\infty}\text{d}{t}_{1}{\mathbb{R}}^{(3)}(\mathit{r},{\mathit{r}}_{3},{\mathit{r}}_{2},{\mathit{r}}_{1};{t}_{3},{t}_{2},{t}_{1})\times \mathit{A}({\mathit{r}}_{3},t-{t}_{3})\mathit{A}({\mathit{r}}_{2},t-{t}_{3}-{t}_{2})\mathit{A}({\mathit{r}}_{1},t-{t}_{3}-{t}_{2}-{t}_{1}),$$

(238)

where *t _{j}* are the delays between successive interactions with the fields and the third-order response function

$$\mathit{J}(\mathit{r})=\stackrel{.}{\mathit{P}}(\mathit{r})+{\nabla}_{r}\times \mathit{M}(\mathit{r}).$$

(239)

We consider an isotropic ensemble of aggregates. We define the center of charge of the *m*th molecule *r** _{m}* and the displacement

$${\mathit{J}}_{{\mathit{k}}_{s}}(t)=exp[-i{\omega}_{s}(t-{\tau}_{3})-i({u}_{2}{\omega}_{2}+{u}_{1}{\omega}_{1})({\tau}_{3}-{\tau}_{2})-{iu}_{1}{\omega}_{1}({\tau}_{2}-{\tau}_{1})]\times \int \int {\int}_{0}^{\infty}\text{d}{t}_{3}\text{d}{t}_{2}\text{d}{t}_{1}{\mathbb{R}}_{{\mathit{k}}_{s},{\mathit{k}}_{3},{\mathit{k}}_{2},{\mathit{k}}_{1}}^{(3)}({t}_{3},{t}_{2},{t}_{1})\times exp[i{\omega}_{s}{t}_{3}+i({u}_{2}{\omega}_{2}+{u}_{1}{\omega}_{1}){t}_{2}+{iu}_{1}{\omega}_{1}{t}_{1}]\times {\mathbb{A}}_{3}^{{u}_{3}}(t-{t}_{3}){\mathbb{A}}_{2}^{{u}_{2}}(t-{t}_{3}-{t}_{2}){\mathbb{A}}_{1}^{{u}_{1}}(t-{t}_{3}-{t}_{2}-{t}_{1}).$$

(240)

Here, *k** _{s}* =

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

$${\widehat{\mathit{J}}}_{m}(\mathit{r})={\mathit{J}}_{m}^{\ast}(\mathit{r}){\widehat{B}}_{m}+{\mathit{J}}_{m}(\mathit{r}){\widehat{B}}_{m}^{\u2020}.$$

(241)

Fourier transform with respect to ** r**,

$${\mathit{J}}_{m}(u\mathit{k})\approx iu\omega \left({\mathit{\mu}}_{m}+iu\mathit{k}\xb7{\mathit{q}}_{m}+iu(\mathit{k}\xb7{\mathit{R}}_{m}){\mathit{\mu}}_{m}-\frac{\mathit{k}}{\omega}\times {\mathit{m}}_{m}\right).$$

(242)

All molecular properties have been calculated with respect to the origin *R** _{m}*. The

$${\mathit{q}}_{e}^{{\nu}_{2}{\nu}_{1}}=\sum _{n}{\psi}_{ne}{\mathit{q}}_{n}^{(i){\nu}_{2}{\nu}_{1}}+\sum _{n}{\psi}_{ne}{\mathit{R}}_{n}^{{\nu}_{2}}{\mathit{\mu}}_{n}^{{\nu}_{1}}.$$

(243)

We note that, while in the local molecular basis, we have
${\mathit{q}}_{n}^{(i){\nu}_{2}{\nu}_{1}}={\mathit{q}}_{n}^{(i){\nu}_{1}{\nu}_{2}}$; this is no longer the case for the delocalized excitons basis, where
${\mathit{q}}_{e}^{{\nu}_{2}{\nu}_{1}}\ne {\mathit{q}}_{e}^{{\nu}_{1}{\nu}_{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

$${\mathbb{R}}^{(1)}(t)={\left(\frac{i}{\hslash}\right)}^{1}\sum _{mn}\langle {\mathit{J}}_{m}^{{\nu}_{2}^{\ast}}(-\mathit{k}){\mathit{J}}_{n}^{{\nu}_{1}}(-\mathit{k})\rangle {G}_{mn}(t)+c.c.$$

(244)

The QP third-order response function in the coherent regime for the *k*_{I} =−*k*_{1} + *k*_{2} + *k*_{3} (eq 89) technique similarly reads

$${\mathbb{R}}_{{\mathit{k}}_{\text{I}}}^{(\text{CED})}({t}_{3},{t}_{2},{t}_{1})=2{\left(\frac{i}{\hslash}\right)}^{3}\langle {\mathit{J}}_{{n}_{4}}^{{\nu}_{4}^{\ast}}(-{\mathit{k}}_{s}){\mathit{J}}_{{n}_{3}}^{{\nu}_{3}^{\ast}}(-{\mathit{k}}_{3}){\mathit{J}}_{{n}_{2}}^{{\nu}_{2}^{\ast}}(-{\mathit{k}}_{2}){\mathit{J}}_{{n}_{1}}^{{\nu}_{1}^{\ast}}(-{\mathit{k}}_{1})\rangle \times {\int}_{0}^{{t}_{3}}\text{d}\tau {\int}_{0}^{\tau}\text{d}{\tau}^{\prime}{G}_{{n}_{4}{n}_{4}^{\prime}}({t}_{3}-\tau ){\mathrm{\Gamma}}_{{n}_{4}^{\prime}{n}_{1}^{\prime},{n}_{3}^{\prime}{n}_{2}^{\prime}}(\tau -{\tau}^{\prime})\times {G}_{{n}_{4}^{\prime}{n}_{3}}({\tau}^{\prime}){G}_{{n}_{2}^{\prime}{n}_{2}}({t}_{2}+{\tau}^{\prime}){G}_{{n}_{1}^{\prime}{n}_{1}}^{\ast}({t}_{1}+{t}_{2}+\tau ).$$

(245)

Here,

$$\langle {\mathit{J}}_{2}^{{\nu}_{2}}({\mathit{k}}_{2}){\mathit{J}}_{1}^{{\nu}_{1}}({\mathit{k}}_{1})\rangle =\langle ({\mathit{J}}_{2}({\mathit{k}}_{2})\xb7{\mathit{\nu}}_{2})({\mathit{J}}_{1}({\mathit{k}}_{1})\xb7{\mathit{\nu}}_{1})\rangle $$

(246)

and

$$\langle {\mathit{J}}_{4}^{{\nu}_{4}}({\mathit{k}}_{4}){\mathit{J}}_{3}^{{\nu}_{3}}({\mathit{k}}_{3}){\mathit{J}}_{2}^{{\nu}_{2}}({\mathit{k}}_{2}){\mathit{J}}_{1}^{{\nu}_{1}}({\mathit{k}}_{1})\rangle =\langle \prod _{j=1}^{4}({\mathit{J}}_{j}({\mathit{k}}_{j})\xb7{\mathit{\nu}}_{j})\rangle $$

(247)

is the orientationally averaged product: *J** _{j}*(

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

$$\langle {\mathit{J}}_{4}({\mathit{k}}_{S}){\mathit{J}}_{3}(-{u}_{3}{\mathit{k}}_{3}){\mathit{J}}_{2}(-{u}_{2}{\mathit{k}}_{2}){\mathit{J}}_{1}(-{u}_{1}{\mathit{k}}_{1})\rangle =-{u}_{3}{u}_{2}{u}_{1}{\omega}_{4}{\omega}_{3}{\omega}_{2}{\omega}_{1}\{\langle {\mathit{\mu}}_{4}{\mathit{\mu}}_{3}{\mathit{\mu}}_{2}{\mathit{\mu}}_{1}\rangle +i\langle ({\mathit{k}}_{S}\xb7{\mathit{q}}_{4}){\mathit{\mu}}_{3}{\mathit{\mu}}_{2}{\mathit{\mu}}_{1}\rangle -\frac{1}{{\omega}_{S}}\langle ({\mathit{k}}_{S}\times {\mathit{m}}_{4}){\mathit{\mu}}_{3}{\mathit{\mu}}_{2}{\mathit{\mu}}_{1}\rangle -{iu}_{3}\langle {\mathit{\mu}}_{4}({\mathit{k}}_{3}\xb7{\mathit{q}}_{3}){\mathit{\mu}}_{2}{\mathit{\mu}}_{1}\rangle -\frac{1}{{\omega}_{3}}\langle {\mathit{\mu}}_{4}({\mathit{k}}_{3}\times {\mathit{m}}_{3}){\mathit{\mu}}_{2}{\mathit{\mu}}_{1}\rangle -{iu}_{2}\langle {\mathit{\mu}}_{4}{\mathit{\mu}}_{3}({\mathit{k}}_{2}\xb7{\mathit{q}}_{2}){\mathit{\mu}}_{1}\rangle -\frac{1}{{\omega}_{2}}\langle {\mathit{\mu}}_{4}{\mathit{\mu}}_{3}({\mathit{k}}_{2}\times {\mathit{m}}_{2}){\mathit{\mu}}_{1}\rangle -{iu}_{1}\langle {\mathit{\mu}}_{4}{\mathit{\mu}}_{3}{\mathit{\mu}}_{2}({\mathit{k}}_{1}\xb7{\mathit{q}}_{1})\rangle -\frac{1}{{\omega}_{1}}\langle {\mathit{\mu}}_{4}{\mathit{\mu}}_{3}{\mathit{\mu}}_{2}({\mathit{k}}_{1}\times {\mathit{m}}_{1})\rangle .$$

(248)

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

$$\langle {\mathit{J}}_{4}^{{\nu}_{4}}{\mathit{J}}_{3}^{{\nu}_{3}}{\mathit{J}}_{2}^{{\nu}_{2}}{\mathit{J}}_{1}^{{\nu}_{1}}\rangle =-{i}^{4}({u}_{3}{u}_{2}{u}_{1})({\omega}_{4}{\omega}_{3}{\omega}_{2}{\omega}_{1})\langle {\mathit{\mu}}_{4}^{{\nu}_{4}}{\mathit{\mu}}_{3}^{{\nu}_{3}}{\mathit{\mu}}_{2}^{{\nu}_{2}}{\mathit{\mu}}_{1}^{{\nu}_{1}}\rangle .$$

(249)

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: −*i*^{4} *ω _{s}u*

Orientational averaging of the transition dipoles performed using eq 424 in Appendix L requires the calculation of the vectors
${\mathit{F}}_{\{\nu \}}^{(4)}$ in eq 429. The elements of this vector define the three linearly independent configurations of transition dipoles,
${\mathit{F}}_{\{\nu \}}^{(4)}={(1,0,0)}^{\text{T}}$, or
${\mathit{F}}_{\{\nu \}}^{(4)}={(0,1,0)}^{\text{T}}$, or
${\mathit{F}}_{\{\nu \}}^{(4)}={(0,0,1)}^{\text{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.

Linearly Independent NC Third-Order Experiments for *k*_{I}, *k*_{II}, and *k*_{III} 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.

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 (*u*_{1}, *u*_{2}, *u*_{3}) = (−1, 1, 1) for *k*_{I}, (1, −1, 1) for *k*_{II}, and (1, 1, −1) for *k*_{III}. From eq 248, we get

$$\langle {\mathit{J}}_{4}^{{\nu}_{4}}{\mathit{J}}_{3}^{{\nu}_{3}}({u}_{3}){\mathit{J}}_{2}^{{\nu}_{2}}({u}_{2}){\mathit{J}}_{1}^{{\nu}_{1}}({u}_{1})\rangle =-{i}^{4}{\omega}_{s}{u}_{3}{\omega}_{3}{u}_{2}{\omega}_{2}{u}_{1}{\omega}_{1}\times \{\langle {\mathit{\mu}}_{4}^{{\nu}_{4}}{\mathit{\mu}}_{3}^{{\nu}_{3}}{\mathit{\mu}}_{2}^{{\nu}_{2}}{\mathit{\mu}}_{1}^{{\nu}_{1}}\rangle +ik[\langle {\mathit{q}}_{4}^{z,{\nu}_{4}}{\mathit{\mu}}_{3}^{{\nu}_{3}}{\mathit{\mu}}_{2}^{{\nu}_{2}}{\mathit{\mu}}_{1}^{{\nu}_{1}}\rangle -{u}_{3}\langle {\mathit{\mu}}_{4}^{{\nu}_{4}}{\mathit{q}}_{3}^{z,{\nu}_{3}}{\mathit{\mu}}_{2}^{{\nu}_{2}}{\mathit{\mu}}_{1}^{{\nu}_{1}}\rangle -{u}_{2}\langle {\mathit{\mu}}_{4}^{{\nu}_{4}}{\mathit{\mu}}_{3}^{{\nu}_{3}}{\mathit{q}}_{2}^{z,{\nu}_{2}}{\mathit{\mu}}_{1}^{{\nu}_{1}}\rangle -{u}_{1}\langle {\mathit{\mu}}_{4}^{{\nu}_{4}}{\mathit{\mu}}_{3}^{{\nu}_{3}}{\mathit{\mu}}_{2}^{{\nu}_{2}}{\mathit{q}}_{1}^{z,{\nu}_{1}}\rangle ]+\frac{i}{c}\sum _{\alpha}[{\epsilon}_{{\nu}_{4}z\alpha}\langle {\mathit{m}}_{4}^{\u2033\alpha}{\mathit{\mu}}_{3}^{{\nu}_{3}}{\mathit{\mu}}_{2}^{{\nu}_{2}}{\mathit{\mu}}_{1}^{{\nu}_{1}}\rangle -{u}_{3}{\epsilon}_{{\nu}_{3}z\alpha}\langle {\mathit{\mu}}_{4}^{{\nu}_{4}}{\mathit{m}}_{3}^{\u2033\alpha}{\mathit{\mu}}_{2}^{{\nu}_{2}}{\mathit{\mu}}_{1}^{{\nu}_{1}}\rangle -{u}_{2}{\epsilon}_{{\nu}_{2}z\alpha}\langle {\mathit{\mu}}_{4}^{{\nu}_{4}}{\mathit{\mu}}_{3}^{{\nu}_{3}}{\mathit{m}}_{2}^{\u2033\alpha}{\mathit{\mu}}_{1}^{{\nu}_{1}}\rangle -{u}_{1}{\epsilon}_{{\nu}_{1}z\alpha}\langle {\mathit{\mu}}_{4}^{{\nu}_{4}}{\mathit{\mu}}_{3}^{{\nu}_{3}}{\mathit{\mu}}_{2}^{{\nu}_{2}}{\mathit{m}}_{1}^{\u2033\alpha}\rangle ],$$

(250)

where *k* = |** k**| is the wavevector amplitude and

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).

The resonant response functions in noncollinear configurations depend on both pulse polarizations and wavevectors, (*k*_{4},*k*_{3},*k*_{2},*k*_{1}; *ν*_{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:

$${\mathit{J}}_{{\mathit{k}}_{s}}^{{\nu}_{4}}\propto {\mathbb{R}}_{{\mathit{k}}_{s},{\mathit{k}}_{3},{\mathit{k}}_{2},{\mathit{k}}_{1}}^{(3)}({t}_{3},{t}_{2},{t}_{1})\times {\mathbb{A}}_{3}^{{u}_{3}}(t-{t}_{3}){\mathbb{A}}_{2}^{{u}_{2}}(t-{t}_{3}-{t}_{2}){\mathbb{A}}_{1}^{{u}_{1}}(t-{t}_{3}-{t}_{2}-{t}_{1})=\sum _{{\nu}_{3}=x,y,z}\cdots \sum _{{\nu}_{1}=x,y,z}{\mathbb{R}}_{{\mathit{k}}_{s};{\nu}_{4}{\nu}_{3}{\nu}_{2}{\nu}_{1}({\mathit{k}}_{s}{\mathit{k}}_{3}{\mathit{k}}_{2}{\mathit{k}}_{1})}^{(3)}({t}_{3},{t}_{2},{t}_{1})\times {\mathbb{A}}_{3}^{{u}_{3}{\nu}_{3}}(t-{t}_{3}){\mathbb{A}}_{2}^{{u}_{2}{\nu}_{2}}(t-{t}_{3}-{t}_{2}){\mathbb{A}}_{1}^{{u}_{1}{\nu}_{1}}(t-{t}_{3}-{t}_{2}-{t}_{1}).$$

(251)

Here, *ν*_{4}*ν*_{3}*ν*_{2}*ν*_{1}(*k*_{4}*k*_{3}*k*_{2}*k*_{1}) 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
${\mathit{F}}_{\{\nu \}}^{(5)}$ in eq 430 defines six linearly independent techniques. These are obtained by specifying
${\mathit{F}}_{\{\nu \}}^{(5)}={(0,0,0,0,0,1)}^{\text{T}}$, (0, 0, 0, 0, 1, 0)^{T}, etc. The six noncollinear wavevector and polarization configurations that lead to this
${\mathit{F}}_{\{\nu \}}^{(5)}$^{156} are listed in Table 3.

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 *t*_{2}. The *xxxx* (dipole approximation) is displayed in the top row. At short time (*t*_{2} ≤ 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 *t*_{2} 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.

Signals specifically designed to probe coherent and dissipative dynamics^{81} are shown in Figure 27. The combination
${\stackrel{\sim}{B}}_{1}={S}_{\mathit{xyzx}(zx\overline{x}\overline{z})}^{(3)}-{S}_{\mathit{xyxx}(\mathit{zxzx})}^{(3)}$ selectively probes exciton-intraband coherences and their dynamics, since the contributions from populations exactly cancel out. On the other hand, the signal
${\stackrel{\sim}{C}}_{1}={S}_{\mathit{xxyz}(\mathit{zzxx})}^{(3)}-{S}_{\mathit{xyzx}(zx\overline{x}\overline{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.

The _{1} signal shows a much slower dynamics. At zero delay, the signal is simply the inverse of _{1}, but the differences rapidly develop at longer delays. The _{1} peak structure is richer than _{1}, and specific peaks can be attributed to coherences in ESE and ESA processes by comparing both _{1} and _{1}: the peaks that coincide in _{1} and _{1} are from ESA, and the ones missing in _{1} are from ESE. Peaks in _{1}, which are absent in _{1}, come from populations. Note that coherences decay within 150 fs. Thus, the longer time dynamics is now solely due to populations. The _{1} signal vanishes in the absence of population transfer in the ESE diagram. Thus, the peaks in _{1} 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).

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.^{267}^{–}^{270} Coherent-control pulse-shaping techniques have been used to drive quantum systems into a desired state,^{271}^{–}^{275} to manipulate excitons in multidimensional spectroscopy,^{96}^{,}^{97}^{,}^{262}^{,}^{263}^{,}^{276}^{–}^{282} and to control the function of biological systems.^{283}^{–}^{286}

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

$$\mathit{E}(t)=\sum _{j=1}^{4}\sum _{{u}_{j}=\pm 1}\sum _{\nu =x,y,z}\mathit{\nu}{\mathbb{E}}_{j\nu}^{{u}_{j}}(t-{\tau}_{j})\times exp[{iu}_{j}({\mathit{k}}_{j}\mathit{r}-{\omega}_{j}(t-{\tau}_{j})-{\varphi}_{j\nu}(t-{\tau}_{j}))].$$

(252)

Compared to eq 16, here we have added the temporal phase function *ϕ _{jv}* (

The electric field can be alternatively recast in the frequency domain

$$\mathit{E}(\omega )=\sum _{j}\sum _{{u}_{j}}\sum _{\nu}\mathit{\nu}{\stackrel{\sim}{\mathbb{E}}}_{j\nu}^{{u}_{j}}({u}_{j}\omega -{\omega}_{j})\times exp[{iu}_{j}{\mathit{k}}_{j}\mathit{r}+i\omega {\tau}_{j}+{iu}_{j}{\stackrel{\sim}{\varphi}}_{j\nu}({u}_{j}\omega -{\omega}_{j})],$$

(253)

where
${\stackrel{\sim}{\mathbb{E}}}_{j\nu}^{+}(\omega )\phantom{\rule{0.16667em}{0ex}}exp(i{\stackrel{\sim}{\varphi}}_{j\nu}(\omega ))$ is the Fourier transform of the complex envelope function
${\stackrel{\sim}{\mathbb{E}}}_{j\nu}^{+}(t)\phantom{\rule{0.16667em}{0ex}}exp(-i{\varphi}_{j\nu}(t))$ and = []*. Note that (*ω*) and (*ω*) are *not* the Fourier transforms of (*t*) and *ϕ* (*t*).

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

$${S}_{{\mathit{k}}_{s},{\nu}_{s}}({\tau}_{s},{\tau}_{3},{\tau}_{2},{\tau}_{1})={\text{e}}^{-i{\omega}_{s}({\tau}_{s}-{\tau}_{3})-i({u}_{2}{\omega}_{2}+{u}_{1}{\omega}_{1})({\tau}_{3}-{\tau}_{2})-{iu}_{1}{\omega}_{1}({\tau}_{2}-{\tau}_{1})}\times {\int}_{-\infty}^{+\infty}\text{d}t\int \int {\int}_{0}^{\infty}\text{d}{t}_{3}\text{d}{t}_{2}\text{d}{t}_{1}\sum _{{\nu}_{3}{\nu}_{2}{\nu}_{1}}{R}_{{\mathit{k}}_{s};{\nu}_{s}{\nu}_{3}{\nu}_{2}{\nu}_{1}}^{(3)}({t}_{3},{t}_{2},{t}_{1})\times {\text{e}}^{i{\omega}_{s}{t}_{3}+i({u}_{2}{\omega}_{2}+{u}_{1}{\omega}_{1}){t}_{2}+{iu}_{1}{\omega}_{1}{t}_{1}}\times {\mathbb{E}}_{s,{\nu}_{s}}^{-}(t-{\tau}_{s}){\mathbb{E}}_{3,{\nu}_{3}}^{{u}_{3}}(t-{t}_{3}-{\tau}_{3}){\mathbb{E}}_{2,{\nu}_{2}}^{{u}_{2}}(t-{t}_{3}-{t}_{2}-{\tau}_{2})\times {\mathbb{E}}_{1,{\nu}_{1}}^{{u}_{1}}(t-{t}_{3}-{t}_{2}-{t}_{1}-{\tau}_{1}).$$

(254)

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 (*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.

The field of coherent control with laser pulse shaping has grown in many directions since the original proposals in 1980s^{292}^{–}^{295} and the first experimental realizations in the early 1990s.^{296}^{–}^{299} 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.^{300}^{–}^{304}

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}^{,}^{305}^{–}^{308} 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.^{310}^{–}^{312}

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.^{315}^{–}^{317} 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 iodine^{320} and potassium dimers^{321} 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 complexes^{278} and to control the exciton dynamics in FMO^{276} 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.^{324}^{–}^{327} Nelson et al. constructed a 2D femtosecond pulse-shaping apparatus^{325}^{–}^{327} 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 *i _{j}*(

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}*(

$${\chi}_{j}(t)=arctan\frac{{\mathbb{E}}_{jy}(t)}{{\mathbb{E}}_{jx}(t)}.$$

(255)

A second angle *δ _{j}*(

$${\delta}_{j}(t)={\varphi}_{jy}(t)-{\varphi}_{jx}(t),$$

(256)

* _{j}*(

$${\stackrel{\sim}{\theta}}_{j}(t)=\frac{1}{2}arctan[tan(2{\chi}_{j}(t))cos{\delta}_{j}(t)].$$

(257)

The orientation angle of the ellipse *θ _{j}*(

$${\theta}_{j}(t)=\{\begin{array}{ll}{\stackrel{\sim}{\theta}}_{j}(t)\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{\chi}_{j}(t)\le \pi /4,\hfill \\ {\stackrel{\sim}{\theta}}_{j}(t)+\pi /2\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{\chi}_{j}(t)>\pi /4\wedge {\stackrel{\sim}{\theta}}_{j}(t)<0,\hfill \\ {\stackrel{\sim}{\theta}}_{j}(t)-\pi /2\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{\chi}_{j}(t)>\pi /4\wedge {\stackrel{\sim}{\theta}}_{j}(t)\ge 0.\hfill \end{array}$$

The ellipticity angle *ε _{j}*(

$${\epsilon}_{j}(t)=\frac{1}{2}arcsin[sin(2{\chi}_{j}(t))sin{\delta}_{j}(t)].$$

(258)

Finally, the total phase *ϕ _{j}*(

$${\varphi}_{j}(t)={\varphi}_{jx}(t)+\text{sign}[{\theta}_{j}(t){\epsilon}_{j}(t)]\times arccos\left[\sqrt{\frac{{I}_{j}(t)}{{[{\mathbb{E}}_{jx}(t)]}^{2}}}cos{\theta}_{j}(t)cos{\epsilon}_{j}(t)\right].$$

(259)

The instantaneous pulse frequency *ω _{j}*(

$${\omega}_{j}(t)={\overline{\omega}}_{j}+\frac{\text{d}{\varphi}_{j}(t)}{\text{d}t}.$$

(260)

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).

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_{}) 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 (*t*) and (*t*) are varied holding the total temporal intensity profile
${I}_{j}(t)={\mathbb{E}}_{jx}^{2}(t)+{\mathbb{E}}_{jy}^{2}(t)$ as well as the spectral amplitudes of both electric field components (*ω*) and
${I}_{j}(\omega )={\stackrel{\sim}{\mathbb{E}}}_{jx}^{2}(\omega )+{\stackrel{\sim}{\mathbb{E}}}_{jy}^{2}(\omega )$ fixed to the initial values. The temporal (*ϕ _{jv}*(

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 *k*_{1}, local oscillator *k** _{s}* pulse, and a second pulse

$$\begin{array}{l}{S}_{{\mathit{k}}_{\text{I}}}({\mathrm{\Omega}}_{1},{\mathrm{\Omega}}_{2})=\sum _{\alpha ,\beta ,\gamma ,\delta}{\int}_{0}^{\infty}{\int}_{0}^{\infty}\text{d}t\phantom{\rule{0.16667em}{0ex}}\text{d}\tau \phantom{\rule{0.16667em}{0ex}}{\text{e}}^{i{\mathrm{\Omega}}_{1}\tau +i{\mathrm{\Omega}}_{2}t}{\int}_{0}^{\infty}\text{d}{t}_{3}{\int}_{0}^{\infty}\text{d}{t}_{2}{\int}_{0}^{\infty}\text{d}{t}_{1}\times \\ {R}_{\delta \gamma \beta \alpha}^{(3)}({t}_{3},{t}_{2},{t}_{1}){\text{e}}^{i(2{\omega}_{2}-{\omega}_{1}){t}_{3}+i({\omega}_{2}-{\omega}_{1}){t}_{2}-i{\omega}_{1}{t}_{1}}\times {\mathbb{E}}_{3\delta}^{\ast}(t){\mathbb{E}}_{2\gamma}(t-{t}_{3}){\mathbb{E}}_{2\beta}(t-{t}_{3}-{t}_{2}){\mathbb{E}}_{1\alpha}^{\ast}(t+\tau -{t}_{3}-{t}_{2}-{t}_{1}),\end{array}$$

(261)

where (*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
$\mathrm{\Delta}{\mathrm{\Omega}}_{1}^{\prime}$ and
$\mathrm{\Delta}{\mathrm{\Omega}}_{3}^{\prime}$, 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 *j*th peak.

$${I}_{j}\equiv {\int}_{-\delta}^{\delta}\text{d}{\mathrm{\Omega}}_{1}{\int}_{-\delta}^{\delta}\text{d}{\mathrm{\Omega}}_{2}{S}_{{\mathit{k}}_{\text{I}}}({\mathrm{\Omega}}_{1}+{\mathrm{\Omega}}_{1}^{j},{\mathrm{\Omega}}_{2}+{\mathrm{\Omega}}_{2}^{j}).$$

(262)

The ratios of the integrated diagonal and cross-peaks *T*_{1} = (*I*_{D1} + *I*_{D2})/*I*_{C1} and *T*_{2} = (*I*_{D1} + *I*_{D2})/*I*_{C2} 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 *T*_{1}) and STP′ (target *T*_{2}). The resulting polarization-shaped pulses are displayed using the quasi-3D electric-field representation based on specifying the temporal intensity *I _{j}*(

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 _{ks;ν4ν3ν2ν1}.^{79} The nonchiral (NC) tensor components (dipole approximation) include *ν*_{4}*ν*_{3}*ν*_{2}*ν*_{1} *xxyy*, *xyxy*, and *xyyx* (note that *xxxx* = *xxyy* + *xyxy* + *xyyx*). We label them
${T}_{1}^{(\text{NC})}-{T}_{3}^{(\text{NC})}$. The nine linearly independent chirality-induced (CI) *k*_{I} = −*k*_{1} + *k*_{2} + *k*_{3} response function tensor components are labeled
${T}_{1}^{(\text{CI})}-{T}_{9}^{(\text{CI})}$. These were defined in Tables 1–3 in section 13.

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

$${\mathbb{S}}^{(\text{NC}/\text{CI})}({\omega}_{1},{\omega}_{3})=\sum _{j=1}^{n}{c}_{j}{T}_{j}^{(\text{NC}/\text{CI})}$$

(263)

where *n* = 3 for the NC and *n* = 9 the CI techniques. The complex coefficients
${c}_{i}={c}_{i}^{\prime}+i{c}_{i}^{\u2033}$ 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 *t*_{2} = 0 were simulated with a homogeneous broadening *γ* = 500 cm^{−1}. Only two diagonal peaks, *P*_{1} and *P*_{5}, and two cross-peaks, *P*_{2} and *P*_{4}, were resolved. In the chirality-induced *T*_{1}–*T*_{9} spectra, the remaining four cross-peaks *P*_{3}, *P*_{6}, *P*_{7}, and *P*_{8} 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 *S*_{2}-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.

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 *t*_{2} = 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 pathways^{262} 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.

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 *t*_{2} = 5 ps. Two cross-peaks (1,5) and (1,7) contributing to predicted *fast* and *slow* pathways^{40}^{,}^{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 *t*_{2}. 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 *t*_{2} = 5 ps were then used to obtain the 2D spectra for other *t*_{2} 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.

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 *k*_{I} and *k*_{II} techniques, which only differ by the order of the first two interactions (Figures 4, ,5).5). The delay time *t*_{1} is controlled by the pump-field, while *t*_{3} is controlled by the probe. The delay between pump and probe pulses then coincides with the time variable *t*_{2}. The pump–probe signal is thus given by

$${I}_{\text{pp}}({\omega}_{\text{pr}},\tau ,{\omega}_{\text{pu}})=Im[{S}_{{\mathit{k}}_{\text{I}}}({\omega}_{\text{pr}},\tau ,-{\omega}_{\text{pu}})+{S}_{{\mathit{k}}_{\text{II}}}({\omega}_{\text{pr}},\tau ,{\omega}_{\text{pu}})],$$

(264)

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 *k*_{1} and *k*_{2} are applied simultaneously; the probe comes after a delay *τ* and the homodyne signal is measured in the *k*_{1} − *k*_{2} + *k*_{pr} 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

$${I}_{\text{TG}}(\tau ,{\omega}_{\text{pu}})={\int}_{0}^{\infty}\text{d}{t}_{3}{\mid {S}_{{\mathit{k}}_{\text{I}}}({t}_{3},\tau ,-{\omega}_{\text{pu}})+{S}_{{\mathit{k}}_{\text{II}}}({t}_{3},\tau ,{\omega}_{\text{pu}})\mid}^{2}.$$

(265)

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 *k*_{I} signal as a function of two delay times *t*_{1} and *t*_{2}.

$${I}_{3\text{PEPS}}({t}_{2},{t}_{1})={\int}_{0}^{\infty}\text{d}{t}_{3}{\mid {S}_{{\mathit{k}}_{\text{I}}}({t}_{3},{t}_{2},{t}_{1})\mid}^{2}.$$

(266)

This signal is usually displayed as a family of 1D plots vs *t*_{1}, for various values of *t*_{2}. For an ensemble of two-level systems, the signal along *t*_{1} grows for short *t*_{1}, reaches a maximum at *t*_{1} = *τ _{M}*, and eventually decays. The value of

We have introduced two types of response functions: one (*R*) connecting the induced polarization with the incoming electric fields (eq 14), and the other () 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}^{)−1}, where *ω _{j}* is the carrier frequency of the

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 *U _{mnkl}* controls the exciton scattering when the field is turned off, while
${\mathit{\mu}}_{m}^{(2)}$ represents anharmonic system–field interaction. Hard-core bosons are a special case of this model, as shown in Appendix C. They have

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}^{,}^{338}^{–}^{343} 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 UV^{82}^{,}^{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.

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.

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.

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.

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.

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.

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.

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

$$I(z)={I}_{0}(z)exp(-{\kappa}_{a}z).$$

(267)

The absorption coefficient *κ _{a}* is given by

$${\kappa}_{a}(\omega )=\frac{4\pi \omega}{n(\omega )c}Im[{R}^{(1)}(\omega )].$$

(268)

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

$${\kappa}_{a}(\omega )\propto \omega \sum _{e}{\mid {\mathit{\mu}}_{e}\mid}^{2}\frac{{\gamma}_{eg}}{{(\omega -{\omega}_{eg})}^{2}+{\gamma}_{eg}^{2}}.$$

(269)

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:

$$\langle ({\mathit{\mu}}_{2}\xb7{\mathit{\nu}}_{2})({\mathit{\mu}}_{1}\xb7{\mathit{\nu}}_{1})\rangle =\sum _{i,j,\alpha ,\beta =x,y,z}\langle {T}_{i\alpha}{T}_{j\beta}\rangle {\mathit{\nu}}_{2}^{\alpha}{\mathit{\nu}}_{1}^{\beta}{\mathit{\mu}}_{2}^{i}{\mathit{\mu}}_{1}^{j}.$$

(270)

Here, we introduced the transformation tensor from the molecular (*μ*^{α}) to the laboratory frame (*μ** ^{i}*),

$$\begin{array}{l}\langle {T}_{i\alpha}{T}_{j\beta}\rangle \equiv {(2\pi )}^{-1}{\int}_{0}^{2\pi}\text{d}\phi {T}_{i\alpha}(\phi ){T}_{j\beta}(\phi )\\ ={2}^{-1}[2{\delta}_{iz}{\delta}_{\alpha z}{\delta}_{jz}{\delta}_{\beta z}+{\delta}_{i\alpha}{\zeta}_{iz}{\delta}_{j\beta}{\zeta}_{jz}+({\delta}_{jy}{\delta}_{\beta x}-{\delta}_{jx}{\delta}_{\beta y})({\delta}_{iy}{\delta}_{\alpha x}-{\delta}_{ix}{\delta}_{\alpha y})].\end{array}$$

(271)

Thus, when *μ*_{2} = *μ*_{1} =** μ**, we obtain the following possible field configurations. For

$$\langle (\mathit{\mu}\xb7\mathit{z})(\mathit{\mu}\xb7\mathit{z})\rangle ={({\mathit{\mu}}^{z})}^{2}$$

(272)

and for *ν*_{2} = *ν*_{1} = ** x** (the same holds for

$$\langle (\mathit{\mu}\xb7\mathit{x})(\mathit{\mu}\xb7\mathit{x})\rangle ={2}^{-1}[{({\mathit{\mu}}^{x})}^{2}+{({\mathit{\mu}}^{y})}^{2}].$$

(273)

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

$${\kappa}_{\text{LD}}(\omega )={\kappa}_{\text{L}}(\omega )-{\kappa}_{\text{T}}(\omega ),$$

(274)

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

$${\kappa}_{\text{L}}=\omega \sum _{e}\frac{{({\mathit{\mu}}_{e}^{z})}^{2}{\gamma}_{eg}}{{(\omega -{\omega}_{eg})}^{2}+{\gamma}_{eg}^{2}},$$

(275)

whereas *κ*_{T} is the transverse absorption (eq 269) together with eq 273:

$${\kappa}_{\text{T}}=\omega \sum _{e}\frac{{\gamma}_{eg}}{2}\frac{{({\mathit{\mu}}_{e}^{x})}^{2}+{({\mathit{\mu}}_{e}^{y})}^{2}}{{(\omega -{\omega}_{eg})}^{2}+{\gamma}_{eg}^{2}}.$$

(276)

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:

$$\mathit{E}(t)=\frac{1}{2}(\mathit{x}\mp i\mathit{y})exp(\mathit{ikz}-i\omega t)+c.c.,$$

(277)

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

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

$${\mathit{P}}_{\text{R}/\text{L}}=(\mathit{x}\mp i\mathit{y})({R}_{xx}^{(1)}(\omega )\mp {iR}_{xy}^{(1)}(\omega ))exp(\mathit{ikz}-i\omega t)+c.c.$$

(278)

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

$$tan(\theta )=\frac{{E}_{\text{R}}-{E}_{\text{L}}}{{E}_{\text{R}}+{E}_{\text{L}}},$$

(279)

where *E*_{R/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

$$\theta (\omega )=-\frac{Re{R}_{xy}^{(1)}}{Im{R}_{xx}^{(1)}}.$$

(280)

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

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:

$$\nabla \times \nabla \times \mathit{E}+\frac{1}{{c}^{2}}\frac{{\partial}^{2}\mathit{E}}{\partial {t}^{2}}=-\frac{4\pi}{{c}^{2}}\frac{{\partial}^{2}\mathit{P}}{\partial {t}^{2}}.$$

(281)

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:

$${E}_{s}(t)=\frac{2\pi i}{n({\omega}_{s})}\frac{{\omega}_{s}}{c}l{P}_{s}(t)\phantom{\rule{0.16667em}{0ex}}\text{sinc(}\mathrm{\Delta}\mathit{k}l/2\text{)exp(}i\mathrm{\Delta}\mathit{k}l/2\text{)},$$

(282)

where *n*(*ω _{s}*) is the dielectric constant at frequency

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

$${S}^{(\text{HOD})}=\int \text{d}t{\mid {\mathit{P}}_{s}(t)\mid}^{2}.$$

(283)

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

$${S}^{(\text{TGH})}(t)={\mid {\mathit{P}}_{s}(t)\mid}^{2}.$$

(284)

In heterodyne detection (HET), an additional laser pulse (local oscillator, *E*_{H}) is applied in the direction of the signal field, *E*_{S}. It interferes with the signal, and the resulting field is then detected by homodyne detector. We thus get

$${S}^{(\text{HET})}=\int \text{d}t{\mid {\mathit{E}}_{s}(t)+{\mathit{E}}_{\text{H}}(t)\mid}^{2}\propto 2Re({\mathit{E}}_{\text{H}}^{\ast}(t)\xb7{\mathit{E}}_{s}(t))+{\mid {\mathit{E}}_{\text{H}}(t)\mid}^{2}+{\mid {\mathit{E}}_{s}(t)\mid}^{2}.$$

(285)

The local oscillator (second term) can be readily subtracted. The |*E** _{s}*(

$${S}^{(\text{HET})}(t)=2Re({\mathit{E}}_{\text{H}}^{\ast}(t)\xb7{\mathit{E}}_{s}(t)).$$

(286)

This measurement can be repeated with varying phase of *E*_{H} to reveal the induced polarization *P _{s}*(

Transformations between fermion and boson representations widely used in many-body problems may be employed for describing excitations in semiconductors^{346}^{–}^{353} 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

$${\widehat{H}}_{0}=\hslash \sum _{mn}\sum _{ab}{V}_{ma,nb}{\widehat{b}}_{ma}^{\u2020}{\widehat{b}}_{nb}+\hslash \sum _{\mathit{mnkl}}\sum _{\mathit{abcd}}{W}_{ma,nb}^{kc,ld}{\widehat{b}}_{ma}^{\u2020}{\widehat{b}}_{kc}^{\u2020}{\widehat{b}}_{nb}{\widehat{b}}_{ld},$$

(287)

where *m*, *n* label different molecules and
${\widehat{b}}_{mb}^{\u2020}$ is a bosonic (
$[{\widehat{b}}_{ma},{\widehat{b}}_{nb}^{\u2020}]={\delta}_{mn}{\delta}_{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). *V _{ma}*

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
${W}_{m,m}^{m,m}\to \infty $ 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,

$${\mathrm{\Gamma}}_{mm,nn}(\omega )=-{[{\mathbb{G}}^{(0)}(\omega )]}_{mm,nn}^{-1},$$

(288)

where all other elements vanish. Note that it involves the inversion of × matrix ${\mathbb{G}}_{mm,nn}^{<}$