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}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**), 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.*_{s}^{76}^{,}^{77}2D techniques can also eliminate certain types of inhomogeneous broadenings, show two-exciton resonances, and reveal couplings and correlations between chromophores.^{24}^{,}^{40}^{,}^{41}^{,}^{45}^{,}^{65}^{,}^{70}^{,}^{72}^{,}^{78}^{–}^{81}Extensions to the UV^{82}^{–}^{84}may target backbone transitions of proteins and DNA. In the future, X-ray attosecond techniques may reveal electronic wavepackets with high temporal and spatial resolution.^{85}^{–}^{87}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.