PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Chem Rev. Author manuscript; available in PMC Nov 8, 2010.
Published in final edited form as:
PMCID: PMC2975548
NIHMSID: NIHMS212366
Coherent Multidimensional Optical Spectroscopy of Excitons in Molecular Aggregates; Quasiparticle versus Supermolecule Perspectives
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
* To whom correspondence should be addressed. smukamel/at/uci.edu
University of California Irvine.
Universität Würzburg.
§Charles University.
Molecular aggregates are abundant in nature; they form spontaneously in concentrated solutions and on surfaces and can be synthesized by supramolecular chemistry techniques.13 Assemblies of chromophores play important roles in many biological processes such as light-harvesting and primary charge-separation in photosynthesis.4 Aggregates come in various types of structures: one-dimensional strands (H and J aggregates),2,5 two-dimensional layers,68 circular complexes,9 cylindrical nanotubes, multiwall cylinders and supercylinders,1012 branched fractal structures, dendrimers, and disordered globular complexes.13 They also can be fabricated on substrates by “dip-pen” technology.14 Figure 1 presents some typical structures of pigment-protein complexes found in natural light-harvesting membranes. These are made of chlorophyll and carotenoid chromophores held together by proteins.4,9,15,16
Figure 1
Figure 1
Some biophysical photosynthetic systems: FMO = Fenna-Matthews-Olson protein, a trimer, each unit has 7 chromophores; LH2 = light-harvesting complex 2, double-ring structure of chromophores, 27 chromophores; PS1 = photosystem 1, a complex of 96 chromophores; (more ...)
We consider aggregates made out of chromophores with nonoverlapping charge distributions where intermolecular couplings are purely electrostatic. The optical excitations of such complexes are known as Frenkel excitons.1720 Applications of molecular exciton theory to dimers4,2124 show some key features shared by larger complexes. Their absorption spectra have two bands, whose intensities and Davydov splitting are related to the intermolecular interaction strength and the relative orientation.
Excitons in large aggregates can be delocalized across many chromophores and may show coherent or incoherent energy transfer. The optical absorption, time-resolved fluorescence, and pump-probe spectra in strongly coupled linear J-aggregates have been extensively studied.2,5,18,2529 These spectra contain signatures of cooperative optical response: the absorption splits into several Davydov sub-bands and is shifted compared to the monomer. Additionally, the absorption and fluorescence bands are narrower in strongly coupled J aggregates because of motional and exchange narrowing that dynamically average these properties over the inhomogeneous distribution of energies.30 Aggregates may show cooperative spontaneous emission, superradiance, which results in shorter radiative lifetime than the monomer.3136 Large molecular aggregates may show several optical absorption bands12 whose shapes and linewidths depend on their size and geometry. They undergo elaborate multistep energy-relaxation pathways that can be monitored by time-resolved techniques. These pathways are optimized in natural photosynthetic antennas to harvest light with high speed and efficiency.4,3741 Exciton energy and charge transport, as well as relaxation pathways, are of fundamental interest. Understanding their mechanisms may be used toward the development of efficient, inexpensive substitutes for semiconductor devices.
Nonlinear optical four-wave mixing (FWM) techniques have long been used for probing inter- and intramolecular interactions, excitation energies, vibrational relaxation, and charge-transferpathwaysinmolecularaggregates.Pump-probe, transient gratings, photon-echo, and time-resolved fluorescence were applied to study excited states, their interactions with the environment, and relaxation pathways.4,4252
Optical pulse sequences can be designed for disentangling the spectral features of the coherent nonlinear optical signals. Multidimensional correlation spectroscopy has been widely used in NMR to study the structure and microsecond dynamics of complex molecules.53 It has been proposed to extend these techniques to the femtosecond regime by using optical (Raman) or infrared (IR) lasers tuned in resonance with vibrations.5456 The connection with NMR was established.5759 Multidimensional infrared spectroscopy monitors hydrogen-bonding network and dynamics, protein structures, and their fluctuations.6067 Development of these IR techniques was reviewed recently67 and will not be repeated here.
Multidimensional techniques allow one to study molecular excitons in the visible region and reveal couplings and relaxation pathways. These applications were proposed in refs 6870 and demonstrated experimentally in molecules40,41,65,71,72 and semiconductor quantum wells.7375 Laser phase-locking during excitation and detection is required in these experiments, which measure the signal electric field (both amplitude and phase), not just its intensity. All excitation pulses as well as the detected signal must have a well-defined phase. This review covers multidimensional techniques carried out by applying four femtosecond pulses, as shown in Figure 2, and controlling the three time intervals, t1, t2, and t3, between them. In practice the t3 information is usually obtained interferometrically in one shot rather than by scanning t3. Spectral dispersion of the signal with the local oscillator gives the Fourier transform with respect to t3. Two-dimensional (2D) spectra are displayed by a Fourier transform of these signals with respect to a pair of these time variables. NMR signals do not depend on the wavevectors since the sample is much smaller than the wavelengths of the transitions. Different Liouville space pathways are then separated by combining measurements carried out with different phases of the pulses, rather than by detecting signals in different directions. This phase-cycling technique, which provides the same information as heterodyne detection,5759 has been first used in the optical regime using a collinear pulse configuration.76 Direct heterodyne detection by wave-vector selection was reported in ref 77. By controlling the time-ordering of incoming pulses and the signal wavevector (ks), we obtain a wealth of spectroscopic information about the system. Multidimensional techniques can target a broad variety of physical phenomena. Visible pulses probe electronically excited-state dynamical events: ultrafast intramolecular and intermolecular dephasing and relaxation, energy transport, charge photogeneration, and recombination.76,77 2D techniques can also eliminate certain types of inhomogeneous broadenings, show two-exciton resonances, and reveal couplings and correlations between chromophores.24,40,41,45,65,70,72,7881 Extensions to the UV8284 may target backbone transitions of proteins and DNA. In the future, X-ray attosecond techniques may reveal electronic wavepackets with high temporal and spatial resolution.8587
Figure 2
Figure 2
Coherent third-order nonlinear optical experiment. The four laser pulses are ordered in time; the signal is generated in the phase-matching direction k4. Data processing of time-domain signals and their parametric dependence on the three delays t1, t (more ...)
Multidimensional techniques could further utilize the vector nature of the optical field by selecting specific polarization configurations. These may lead to pulse sequences sensitive to structural chirality (hereafter denoted chirality induced, CI)79,8892 These are 2D extensions of the 1D circular dichroism spectroscopy (CD).93 CI optical signals can be employed to probe specific correlations and couplings of chromophores. These have been demonstrated for probing protein structure in the infrared.94 Numerical sensitivity analysis algorithms and pulse-shaping and coherent-control techniques can help dissect and analyze coherent spectra and simplify congested signals.9597
This review surveys the broad arsenal of theoretical techniques developed toward the description of nonlinear optical signals in molecular complexes. By treating the system-field coupling perturbatively, the signals are expressed in terms of response functions, which allow a systematic classification and interpretation of the various possible signals.
The optical properties of molecular aggregates may be described by the Frenkel exciton model. This model has been first applied to molecular crystals17,98 and subsequently extended to aggregates.4,5,2729,99 The system is partitioned into units (chromophores) with nonoverlapping charge distributions; electron exchange between these units is neglected. The direct product of eigenstates of isolated chromophores forms a convenient basis set for the global excited states. In the Heitler-London approximation, the aggregate ground state is given by the product of ground states of all chromophores.17,99 Single-exciton states are formed by promoting one chromophore to its excited state, keeping all others in their ground state. Their number is equal to the number N of chromophores. Double- and higher-excited states are created similarly. The response is formulated in Liouville space in terms of the system’s density matrix, which allows one to incorporate energy dissipation due to interaction with the environment.52
The methods used for computing the optical response of aggregates may be broadly classified into two types.100 In the supermolecule or Sum Over States (SOS) approach, the response function is expanded in the global eigenstates. Optical spectra are interpreted in terms of transitions between these states. Feynman diagrams, which describe the evolution of the molecular density matrix (Liouville space pathways, LSP), are the key tool in the analysis.52 These provide a direct look at the relevant dynamics at each time interval between interactions with the fields. The signals are interpreted in terms of the relevant density matrix elements of the super-molecule and the sequence of transitions between eigenstates.
The LSP can be divided into two groups, depending on whether the density matrix during the second interval t2 is in a diagonal state (populations) or off diagonal (coherences).81 These groups have symmetry properties associated with permutations of pulse-polarization configurations that lead to distinct signatures in the multidimensional spectrum.
The necessary calculations of excited states make the supermolecule approach computationally expensive. Calculating the global eigenstates is not always feasible. Moreover, sums over states do not offer a simple physical picture when many states are involved.
In the alternative quasiparticle (QP) description, the aggregates are viewed as coupled, localized electronic oscillators. Two-exciton eigenstates are never calculated explicitly; instead, two-exciton resonances are obtained via exciton scattering and calculated using the nonlinear exciton equations (NEE).68,69,101107 Two-exciton propagators are calculated using the QP scattering matrix by solving the Bethe Salpeter equation. The lower computational cost, stemming from the more favorable scaling with size, makes this approach particularly suitable to large aggregates. The dominant contributions to the scattering matrix can be identified a priori by examining the exciton overlaps in real space, providing an efficient truncation strategy.
The spectral lineshapes of 2D signals contain signatures of interactions with phonons, vibrations, and the solvent. Slow and static fluctuations may be treated by statistical averaging. Fast fluctuations can be easily incorporated through population relaxation and dephasing rates. The intermediate fluctuation-time regime requires more careful attention. Many techniques exist for the modeling of quantum dissipative dynamics.108111 These treat the coupling of excitons to a bath at different levels of sophistication. In the simplest scheme, bath fluctuations are assumed to be fast, and second-order perturbation theory for the coupling is used to derive a Pauli master equation in the Markovian limit. Dephasing and energy-transfer processes enter as simple exponential decays. When exciton transport is neglected, energy (diagonal) fluctuations coming from a bath with Gaussian statistics can be exactly incorporated in the response functions using the cumulant expansion for Gaussian fluctuations (CGF).42,43,100,112 Exciton transport may be approximately incorporated into the CGF.24 A more detailed (and computationally expensive) approach is based on the stochastic Liouville equations (SLE), which explicitly include collective bath coordinates in the description.113 The SLE have been used to describe the signatures of chemical-exchange kinetics in coherent 2D signals.
The quasi-particle description of the Frenkel-exciton model is formally similar to that of Wannier excitons in multiband tight-binding models of semiconductors.108,114116 The number of variables is, however, different: since electrons and holes in the Frenkel-exciton model are tightly bound, the number of single-exciton variables coincides with the number N of sites; in contrast, the holes and electrons of Wannier excitons are loosely bound, and the number of single-exciton variables scales as ~N2. This unfavorable scaling is more severe for double-exciton states (~N2 for Frenkel and ~N4 for Wannier). Periodic infinite systems of Frenkel and Wannier excitons may be treated analytically, making use of translational invariance.
Several approaches for computing response functions and multidimensional optical signals are presented in this review. Closed expressions for the signals are derived based on both the QP representation and the supermolecule approach. We use Markovian limit with respect to bath fluctuations in the QP representation. Higher-level CGF and SLE approaches are described within the supermolecule approach. The semiclassical treatment of the molecular coupling with the radiation field whereby the classical radiation field interacts with the quantum exciton system.117 Applications to the Fenna-Matthews-Olson (FMO) photosynthetic complex made of seven coupled bacteriochlorophyll molecules illustrate the various levels of theory.
Section 2 introduces a simple multilevel model system. The response functions are derived in section 3. Section 4 develops the density-matrix formalism and derives the response functions for dissipative quantum systems. The response function theory is connected with experimental heterodyne-detected signals in section 5. Section 6 develops a microscopic model for excitonic aggregates in the quasi-particle representation. In section 7, we derive response functions for 2D signals based on the quasi-particle approach. Different models of system-bath coupling are analyzed. Exciton relaxation and dephasing rates are calculated in section 8, and section 9 provides closed expressions for multidimensional signals that include exciton population transport. Section 10 presents various applications of the quasi-particle theory to the FMO complex. In Section 11, we revisit the supermolecule approach to include slow bath fluctuations and exciton transport. In section 12, we describe a stochastically modulated multilevel system by including bath coordinates explicitly. Spectral line-shape parameters are obtained by solving the stochastic Liouville equations (SLE). Section 13 describes the response functions of an isotropic ensemble of molecular complexes (solutions, liquids). The response functions of the previous sections are extended by orientational averaging going beyond the dipole approximation to include first-order contributions in the optical wavevector. Chirality-induced and nonchiral signals are compared. Section 14 shows how coherent control pulse-shaping algorithms can be used to simplify 2D signals. General discussion and future directions are outlined in section 15.
Most mathematical background and technical details are given in the appendices. Linear optical signals are described in Appendix A. Appendix B describes different modes of signal detection of nonlinear signals. Appendix C provides the relation between a system of coupled two-level oscillators (hard-core bosons) and boson quasi-particles. In Appendix D, we present complete expressions for the response function for the model introduced in section 6. These complement the expressions of section 9. In Appendix E, we describe single- and double-exciton coherent propagation within the QP representation and introduce the scattering matrix. Appendix F derives exciton dephasing and transport rates in the real space representation. The final expressions for the system relaxation rates in terms of the bath spectral density are given in Appendix G. Exciton scattering in nonbosonic systems is described in Appendix H. Appendix I generalizes the quasi-particle description to infinite periodic systems. In Appendix J, we give expressions for the doorway and window functions using cumulant expansion technique for the supermolecule approach. Some quantum correlation functions used in section 11.3 are defined in Appendix K. Appendix L presents closed expressions for orientationally averaged signals.
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
equation M1
(1)
where
equation M2
(2)
Here, ħεa is the energy of state a and Ĥ′ represents the dipole interaction with the external optical electric field E,
equation M3
(3)
where r is the position of the molecule and
equation M4
(4)
is the dipole operator. μab = Σαqαleft angle bracketa|rα|bright angle bracket is the transition dipole between states a and b where α runs over molecular charges qα (electrons and nuclei) with coordinates rα, relative to the center of charge.
The quantum state of an ensemble of optically driven aggregates is described by the density matrix:
equation M5
(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,
equation M6
(6)
Here An external file that holds a picture, illustration, etc.
Object name is nihms212366ig1.jpg, the Liouville superoperator, is given by the commutator, with the Hamiltonian An external file that holds a picture, illustration, etc.
Object name is nihms212366ig1.jpg[rho with circumflex]= [Ĥ, [rho with circumflex]]. Similar to eq 1, we partition the Liouville operator as
equation M7
(7)
where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig2.jpg represents the isolated system and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig45.jpg[rho with circumflex] = [Ĥ′, [rho with circumflex]] is the interaction with the external field.
The retarded Green’s function (forward propagator) describes the free evolution of the molecular density matrix between the interaction events. Setting An external file that holds a picture, illustration, etc.
Object name is nihms212366ig45.jpg = 0, we get
equation M8
(8)
where θ(t) is the Heavyside step-function, defined by equation M9. The density matrix of the driven system (eq 6) can be calculated perturbatively in An external file that holds a picture, illustration, etc.
Object name is nihms212366ig45.jpg by iterative integration of eqs 6 and 7. This yields52
equation M10
(9)
Equation 9 provides the order-by-order expansion of the density matrix in the field and can be recast as [rho with circumflex](t) = [rho with circumflex](0)(t) +[rho with circumflex](1)(t) + [rho with circumflex](2)(t) +... . [rho with circumflex](n)(t) is the density matrix to nth order in the field. The n’th-order induced polarization is the quantity of interest in spectroscopy since it is a source of the signal field. It is given by the expectation value of the dipole operator
equation M11
(10)
The induced polarization is given by P(t) = P(1)(t) + P(2)(t) +.... This polarization is a source of a new optical field. The generated field is calculated for a simple sample geometry in Appendix B. The signal optical field obtained by mixing of incoming laser fields characterized by wavevectors k1, k2, ... may be detected at certain directions ks = ± k1 ± k2 ± ... .
The linear response function connects the linear polarization with the field:52
equation M12
(11)
Since P(1) and E are both vectors, R(1) is a second-rank tensor. For clarity, we avoid tensor notation through most of the review. It will be introduced only when necessary, in section 13. Expressions for the signals related to the linear response function are summarized in Appendix A.
The linear response function is calculated by substituting eq 9 into eq 10:52
equation M13
(12)
where [rho with circumflex]0 is the equilibrium density matrix and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg [rho with circumflex] [equivalent] [[mu], [rho with circumflex]] is the dipole superoperator.
Hereafter, we consider a molecular aggregate made of three-level molecules, whose exciton level scheme is shown in Figure 3..44 This model will be introduced microscopically in section 6. The eigenstates of this model form distinct manifolds (bands). The three lowest manifolds are g, the ground state; e, the single-exciton states; and f, the double-exciton states. The dipole operator only couples g to e (μeg) and e to f (μfe). The number of chromophores (and single-exciton states) is denoted by N. The linear response depends only on the e manifold, whereas the third-order response involves both the e and f bands.
Figure 3
Figure 3
Excitonic aggregate made of coupled three-level chromophores. εn is the excitation energy of each chromophore and 2(εn + Unnnn) is its double-excitation energy. The chromophores are quadratically coupled by Jmn (and quartically by Umnkl (more ...)
Figure 4
Figure 4
Feynman diagrams for the signal generated in the direction kI = −k1 + k2 + k3; ESE = excited-state emission, ESA = excited-state absorption, and GSB = ground-state bleaching. Transport is not included in these diagrams (the coherent limit). Population (more ...)
The linear response function can be calculated by expanding eq 12 in eigenstates,
equation M14
(13)
where the equilibrium density matrix is [rho with circumflex]0 = |gright angle bracketleft angle bracketg|, ωeg = εeεg, and we have added phenomenological dephasing rate γeg to represent the decay of coherence. Dephasing will be treated microscopically for a model of fast bath fluctuations in section 8. Note that equation M15.
Higher-order response functions may be defined and calculated in the same way. Since the second-order response vanishes for our model because of the dipole selection rules, we shall focus on the third-order signals.
The third-order polarization induced in the system by the optical field is related to the incoming electric fields by52
equation M16
(14)
where
equation M17
(15)
is the third-order response function.
This expression represents a sequence of three interactions with the incoming fields described by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg, three free propagation periods between interactions (An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg), and signal generation described by the last [mu]. Since An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg is a commutator and [mu] can act either from the left or the right, eq 15 has 23 = 8 contributions known as Liouville space pathways.
We shall consider time-domain experiments, where the optical electric field consists of several pulses:
equation M18
(16)
Here, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig5.jpg (tτj) is the complex envelope of the jth pulse, centered at time τj with carrier frequency ωj and wavevector kj. The indices uj = ±1 represent the positive (uj = +1) and the negative (uj = −1) frequency components of the field. Note that equation M19. The uj variables allow a compact notation for the response functions. We assume that the pulses are temporally well-separated: the pulse with wavevector k1 comes first, followed by k2 and finally k3.
The field is generated in the directions ±k1 ±k2 ±k3. The detector selects a signal in one direction. To calculate the optical signals, we introduce the wavevector-dependent polarization Pks
equation M20
(17)
where ks = u1k1 + u2k2 + u3k3 is the signal wavevector, and the sum runs over the possible choices of uj = ± 1. We assume resonant excitation and make the rotating-wave approximation (RWA), where only low-frequency terms in eq 14 (where the field and molecular frequency subtract ωjωeg) are retained. High-frequency ωj + ωeg terms make a very small contribution to the polarization and are neglected. This is an excellent approximation for resonant techniques. By combining eq 17 with eqs 14 and 16, we get
equation M21
(18)
with the signal frequency given by ωs = u1ω1 + u2ω2 + u3ω3. The wavevector-dependent response function, equation M22, represents the signal generated in the ks direction. For short pulses, the integrand is finite when t3tτ3, t2τ3τ2, and t1τ2τ1 Therefore, the integral needs to be performed only in the narrow time window specified by pulse envelopes. When the pulses are much shorter than their delays, the pulse delays control the free propagation intervals between interactions described by the Green’s functions.
There are four independent third-order techniques, defined by the various choices of uj:
equation M23
(19)
equation M24
(20)
equation M25
(21)
equation M26
(22)
For our excitonic level scheme and its dipole selection rules, the kIV signal vanishes within the RWA and will not be considered any further.
The resonant third-order response may be interpreted using the Feynman diagrams shown in Figures 46. These represent the sequences of interactions with the various fields as well as the state of the excitonic density matrix during the intervals between interactions.52 The two vertical lines represent the ket (left) and the bra (right) of the density matrix. Time flows from bottom to top; the labels on the graph indicate the density-matrix elements during the free evolution periods between the radiative interactions. An interaction with the optical field is marked by a dot. Each interaction is accompanied by a transition dipole μ factor. The density matrix elements can only change by an interaction with the field. Wavy arrows pointing to the left indicate an interaction with the negative frequency field component An external file that holds a picture, illustration, etc.
Object name is nihms212366ig6.jpg exp[−ikjr], and those to the right indicate interaction with the positive frequency An external file that holds a picture, illustration, etc.
Object name is nihms212366ig7.jpg exp[ikjr]. Arrows coming into the diagram represent absorption of a photon accompanied by a molecular transition to a higher-energy state (g to e or e to f); outgoing arrows represent photon emission and a transition to a lower-energy state (e to g or f to e). The number of interactions, p (red dots), on the right side (right line) of the diagram determines its overall sign, which is (−1)p.
Figure 6
Figure 6
Feynman diagrams as in Figure 5 but for the signal generated in the direction kIII = +k1 + k2k3. ESE and GSB do not contribute in this case, and the technique does not show transport. The diagrams are labeled vii and viii for future reference. (more ...)
The kI signal is depicted in Figure 4. During t1, the aggregate density matrix is in an optical coherence with a characteristic frequency given by the energy difference between states e and g, ωge. During the second interval (t2), it is either in the ground state (gg) (ii) or in the singly excited-state manifold (e′e) (i or iii). During the third interval (t3), it again oscillates with optical frequencies ωe′g (i and ii) or ωfe′ (iii).
The various contributions to the signal are denoted ground-state bleaching (GSB), excited-state stimulated emission (ESE), and excited-state absorption (ESA),52 as marked in Figure 4. This nomenclature reflects the physical processes during the course of interactions: in the GSB pathway, the system returns to the ground-state during the second interval t2; thus, the third interaction reflects the reduction of the ground-state population, which reduces (bleaches) the subsequent absorption. In the ESE pathway, during t2, the system is in the single-excited manifold and the third interaction brings it back to the ground state by stimulated emission. The ESA pathway has the same t1 and t2 dynamics as the ESE; however, the third interaction creates a doubly excited state f.
During the intervals between interactions (t1, t2, and t3), the density matrix evolves freely (no field). For instance, in the ESE (Figure 4) diagram, after the first interaction on the right (bra), the system density matrix in the t1 interval is in the state ρge, and its evolution gives an exp(−get1) factor, where ωge = εgεe is the interband frequency. The second interaction from the left with the ket creates the density matrix element ρe′e, which then evolves as exp(−e′et2). After the third interaction (with the bra), the density-matrix element ρe′g is created and its evolution is exp(−e′gt3). The overall contribution of this pathway to the response function is
equation M27
(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 ab coherence.
Similar to the three contributions to the kI technique shown in Figure 4, kII has three contributions (Figure 5) and kIII has two contributions (Figure 6). In the next section, we provide closed expressions for all techniques by extending this model to include exciton transport. The expressions corresponding to the diagrams shown in Figures 46 may be obtained from eqs 3239 by setting the population relaxation matrix to unity, equation M28.
Figure 5
Figure 5
Feynman diagrams as in Figure 4 but for the signal generated in the direction kII = + k1k2 + k3. Population diagrams are labeled iva, v, and via, and the coherence diagrams are labeled ivb and vib.
The time evolution of closed quantum systems can be expanded in eigenstates that evolve independently by simply acquiring exp(−jt) phase factors. Open systems (i.e., systems coupled to a bath) undergo various types of relaxation and transport processes. These are commonly described by the evolution of the reduced density matrix of the system, where the bath variables have been projected out.111,118,119
The dynamics of the reduced density matrix may be described by a quantum master equation (QME). Its general form may be obtained directly from the quantum mechanical analogue of the classical Langevin equations of motion for open systems.111 The derivation starts by phenomenologically adding fluctuation and dissipation terms to the Schrödinger equation. This yields the SchrödingerLangevin equation,
equation M29
(24)
where ψ is the wave function of the open system, V is an arbitrary set of system operators (independent of ξ) that describe its coupling with the bath, and ξα(t) are white noise variables with zero mean, left angle bracketξα(t) right angle bracket= 0, and short correlation time, equation M30. The form of eq 24 guarantees conservation of the norm of the wave function when averaged over the noise.111
The QME that corresponds to this Schrödinger–Langevin equation is known as the Lindblad equation.120122
equation M31
(25)
It is derived from the Liouville equation for [rho with circumflex](t) [equivalent] | ψ(t) right angle bracket left angle bracket (t)| using eq 24 followed by averaging over the noise variables. This equation may be also derived from firmer foundations that show how it preserves all essential properties of the density matrix: it is hermitian with time-independent trace (Tr([rho with circumflex]) [equivalent] 1) and positive definite.111 It will be shown in section 8 that the microscopic Redfield theory in the secular approximation can be recovered by a specific choice of Vα. The secular Redfield equations decouple the populations (diagonal elements of [rho with circumflex] in the eigenstate basis) from the coherences (off-diagonal elements). The populations then satisfy the Pauli master equation:109,111
equation M32
(26)
Here, Kee, ee is the population transfer rate from state e′ into e. This is an N × N matrix with two pairs of identical indices (ee), (e′e′). In eq 26, the diagonal elements, e = e′, Kee, ee are positive, whereas the off-diagonal elements, ee′, Kee, e′e′ are negative. The rate matrix satisfies ΣeKee, e′e′ = 0 (probability conservation) and detailed-balance Ke2e2, e1e1/Ke1e1, e2e2 = exp(−[variant Planck's over 2pi]ωe2e1/(kBT)); here, kB is the Boltzmann constant and T is the temperature. These conditions are sufficient for the system to reach thermal equilibrium at long times.111
The evolution of diagonal density-matrix elements is described by the population Green’s function, equation M33, 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
equation M34
(27)
where λp is the pth eigenvalue of the left- and right-eigen equations equation M35 and equation M36. χ(L) χ (R) are matrices made out of the left (right) eigenvectors, and D = χ(L) χ (R) is a diagonal matrix. The eigenvectors are organized as columns of χ (R) and rows of χ (L). These two sets of eigenvectors are different since Ke′e′, ee is not Hermitian, but they share the same eigenvalues λp. χ(L) and χ(R) can be normalized as equation M37 for each p; then χ (R) = χ (L)−1 and D is the unitary matrix.
In the secular Redfield equation, the off-diagonal elements satisfy
equation M38
(28)
Its solution is equation M39, where
equation M40
(29)
The first two terms represent population relaxation, and [gamma with circumflex]e,e′ is the pure-dephasing rate. For a compact representation of the response functions, we shall combine the diagonal (eq 26) and off-diagonal (eq 28) blocks into a single tetradic Green’s function representing both coherence and population relaxation:
equation M41
(30)
Comparison of eq 30 with eqs 27 and 28 gives
equation M42
(31)
This equation will be derived microscopically and further extended in section 8.
For the kI and kII techniques, populations are generated during the second time interval t2 in the ESA and ESE diagrams of Figures 4 and and55 when e = e′ (GSB includes ground-state population during t2; since our model contains a single ground state, no population relaxation occurs in the ground state). Population created in state e′ can relax to e during t2 according to eq 26. Coherences decay according to eq 28.
By including this relaxation model, we obtain the following sum over states (SOS) expression for the response function for the kI technique (Figures 4 and and77),
Figure 7
Figure 7
Feynman diagrams for the kI and kII signals, which extend diagrams ia, iiia and iva, via of Figures 4 and and55 by including population transport.
equation M43
(32)
equation M44
(33)
equation M45
(34)
where we have introduced the complex frequencies ξab [equivalent] ωabab. Note that we do not include population relaxation between the excited-state manifold and the ground state; this would change An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg (t2) and would influence the GSB t2 dependence.52
The oscillation frequencies during t1 and t3 have opposite sense. For a two-level system (g and e), the frequency factors cancel out at t1 = t3. This eliminates inhomogeneous broadening and gives rise to the photon-echo signal,70,123 which has been used for studying dephasing of coherent dynamics in molecules and molecular complexes.124,125 For multilevel systems (with a band of e states), this cancelation will occur only for the diagonal e = e′ contributions, in the ESE and GSB pathways during t2 when transport is neglected.
The diagrams shown in Figure 7 extend Figures 4 and and55 to include population transfer during t2, as marked by the pairs of dashed arrows. For e = e′, the diagrams in Figure 7 coincide with Figures 4 and and55.
For the kII technique, we similarly have (Figures 5 and and77)
equation M46
(35)
equation M47
(36)
equation M48
(37)
The two contributions to the kIII (Figure 6) technique are of the ESA type:
equation M49
(38)
equation M50
(39)
Note that exciton transport only enters kI and kII but not kIII where exciton populations are never created.
The present phenomenological model allows an easy calculation of the response functions. Exciton relaxation is described by Markovian rate equations. More elaborate decay patterns that show oscillations and nonexponential correlated dynamics often observed in molecular aggregates41,48,126 require higher levels of theory that will be presented in the coming sections.
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
equation M51
(40)
where the signal depends parametrically on the delays between pulses T3 = τsτ3, T2 = τ3τ2, and T1 = τ2τ1. The notation τ [equivalent] T1, T [equivalent] T2, and t [equivalent] T3 is commonly used. Experimentally the signal field is often dispersed in a spectrometer. The entire T3 dependence is then measured in a single shot rather by scanning the delay between the third pulse and the signal. An external file that holds a picture, illustration, etc.
Object name is nihms212366ig42.jpg(t – τs) is the local oscillator field envelope used for heterodyne detection. The first three pulses are represented by eq 16, and the third-order polarization P(3)(t) is given by eq 18.
Multidimensional signals are displayed in the frequency domain by Fourier transforming equation M52 with respect to the time intervals between the pulses.
equation M53
(41)
Often, one uses a mixed time-frequency representation by performing double-Fourier transform. This gives, e.g., S3, T2, Ω1), etc.
Assuming that all four pulses are temporally well-separated, the integrations can be carried out using the response function eqs 3239.127 Neglecting population transport (Kee, ee [equivalent] 0), this gives the following three contributions to the kI signal (see Figures 46),
equation M54
(42)
equation M55
(43)
equation M56
(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 ξab were defined in sections 2 and 4. The spectral pulse envelopes
equation M57
(45)
are centered around ω = 0.
For the kII signal, we similarly obtain
equation M58
(46)
equation M59
(47)
equation M60
(48)
Finally, the kIII signal is given by
equation M61
(49)
and
equation M62
(50)
Note that resonances occur at both positive and negative Ωj.
Population relaxation may be added, as done in eqs 3239. However, the Ω2 dependence is then more complicated, and it is usually preferable to display the signal S3, T2, Ω1) in the time domain as a function of T2.
Equations 4250 show how the pulse envelopes select the possible transitions allowed by their bandwidths. The pulse envelopes serve as frequency filters, removing all transitions falling outside the pulse bandwidth. If the pulse bandwidth is larger than the exciton bandwidth, we can set &An external file that holds a picture, illustration, etc.
Object name is nihms212366ig43.jpg(ω) = 1. This impulsive signal then coincides with the response function.
This concludes the phenomenological supermolecule description of multidimensional signals in terms of the global eigenstates. In the next section, we present the alternative, quasiparticle, picture. Numerical applications of both approaches will be presented in the following sections.
The calculation of double-exciton states requires the computationally expensive diagonalization of an An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg × An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg matrix, where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg = N (N ± 1)/2 is the number of double-exciton states, and N is the number of molecules. The quasiparticle picture avoids the explicit calculation of the two-exciton states making it more tractable in large aggregates.68,79,102,104,106108,128 This picture naturally emerges out of equations of motion for exciton variables, the NEE, that were derived and gradually developed at different levels.79,107,108 Spano and Mukamel had shown how theories of optical susceptibilities in the frequency domain based on the local mean-field approximation (MFA) can be extended by adding two-exciton variables to properly account for two-exciton resonances.102,128 The nonlinear response is then attributed to exciton–exciton scattering.104,105 The formalism was subsequently extended to molecular aggregates made of three-level molecules and to semiconductors.68,114 Additional relevant variables have been introduced to account for exciton relaxation due to coupling with phonons.68 The phonon degrees of freedom are formally eliminated; they only enter through the relaxation rates. This results in the Redfield equations for the reduced exciton density matrix.
By solving the NEE, we obtain closed Green’s function expressions for the optical response that maintain the complete bookkeeping of time ordering. Applications of this approach were made to J-aggregates,103 pump–probe, photonechoes, and other four-wave mixing techniques of light-harvesting antenna complexes.78,130,131
We first introduce the microscopic exciton model for an aggregate made out of N three-level chromophores. The special case of two-level chromophores is treated in Appendix C. Electronic excitations are expressed using the basis set of the electronic eigenstates of each chromophore: these are the ground state equation M63, the single-excited state equation M64, and the double-excited state equation M65. In the Heitler–London approximation,17 the aggregate ground state is given by a direct product of the ground states of all chromophores,
equation M66
(51)
A single-excited-state basis is constructed by moving one of the chromophores into its excited state,
equation M67
(52)
Double-excitations are obtained either when one of the chromophores is doubly excited
equation M68
(53)
or when two chromophores are singly excited:
equation M69
(54)
Overall, our model has N singly excited states, and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg = N( N + 1)/2 doubly-excited states.
We next introduce an exciton-creation operator on molecule m:131
equation M70
(55)
These operators have the following properties: single-exciton states are created from the ground state equation M71, and double-excitons are obtained either by creating two excitations on different molecules, equation M72, or on the same molecule, equation M73 (the √2 factor is introduced to resemble bosonic exciton properties ( equation M74) within the single- and double-exciton space).
The Hermitian-conjugate, annihilation, operator will be denoted [B with circumflex]m. Operators corresponding to different chromophores commute. The commutation relations of these operators are (see Appendix C)
equation M75
(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, [B with circumflex][B with circumflex][B with circumflex]Φg, will not be considered.
We assume the following model Hamiltonian written in a normally ordered form (i.e., all [B with circumflex] are to the left of [B with circumflex])68
equation M76
(57)
where [variant Planck's over 2pi]εn is the single-excitation energy of molecule n and quadratic coupling [variant Planck's over 2pi]Jmn is responsible for resonant exciton hopping.132134 The quartic couplings [variant Planck's over 2pi]Umn,kl represent various types of anharmonicities that only affect the two-exciton (and higher) manifolds. Terms that do not conserve the number of excitons, such as equation M77 and equation M78, have been neglected in this Hamiltonian (a more general Hamiltonian is discussed in ref 135).
A simplified form of eq 57,
equation M79
(58)
captures the essential physical processes. The diagonal term Δm = 2Umm,mm modifies the energy of two excitations on the same site: left angle bracketΦfmm|Ĥ0fmmright angle bracket/[variant Planck's over 2pi] = 2εm + Δm. Two-level chromophores can be described by setting Δm → ∞. This prevents two excitations from residing on the same site. The couplings An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg/4 = Umn, mn [equivalent]Umn, nm shift the two-exciton energies with respect to the noninteracting excitons left angle bracketΦfmn|Ĥ0fmnright angle bracket/[variant Planck's over 2pi] = εm + εn + An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg. They can be either repulsive, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg > 0, or attractive, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg < 0.
We further assume the following form for the dipole operator where each interaction with the field can create or annihilate a single exciton:
equation M80
(59)
Here, h.c denotes the Hermitian conjugate. The transition dipole left angle bracketΦg|[mu]emright angle bracket = μm gives the transition between the ground and a single-excited state and equation M81 is the transition between a single- and a double-excited state. Within the RWA, the coupling with the field (eq 3) is given by
equation M82
(60)
where
equation M83
(61)
and
equation M84
(62)
are the time-dependent transition amplitudes with μ [equivalent] (μ+)*: μ annihilates the photon and, thus, is conjugated to the creation of exciton, while μ+ creates the photon accompanied by annihilation of the exciton. The total Hamiltonian is given by eq 1 together with eqs 57 and 60.
The exciton dynamics will be calculated starting with the Heisenberg equations of motion
equation M85
(63)
Taking  = [B with circumflex]m and  = [B with circumflex]m[B with circumflex]n gives
equation M86
(64)
equation M87
(65)
with
equation M88
(66)
equation M89
(67)
and Vmn,kl = Umn,kl + Unm,kl. These are the first two members of an infinite hierarchy of coupled equations.
Different-order contributions in the field can be easily sorted out in these equations since interactions with the field change the number of excitons one at a time. Thus, μ is first order, [B with circumflex] is first order and higher, [B with circumflex][B with circumflex][B with circumflex] and equation M90 are third order and higher, and [B with circumflex][B with circumflex][B with circumflex][B with circumflex] is at least fourth order, etc.
When the system is in a pure state, it can be described by a wave function (the expansion equation M91 is sufficient to represent the third-order response), and we have equation M92. We can further factorize equation M93 where we used B = left angle bracketBright angle bracket and Ymn = left angle bracket[B with circumflex]m[B with circumflex]nright angle bracket. Taking expectation value of both sides of eqs 64 and 65, we obtain the NEE to third order in the field107
equation M94
(68)
equation M95
(69)
Quartic products [B with circumflex]†2[B with circumflex]2 contribute only to fourth and higher orders in the fields and were neglected. Within this space of relevant states, the exciton commutation relation, eq 56, can be, therefore, replaced by the boson commutation relation,
equation M96
(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 oscillators99,136,137 as well as fermions.107,138,139 We review the bosonization of two-level molecules in Appendix C. The same ideas can be applied to multilevel molecules (truncated oscillators).
Different levels of approximation will be introduced for truncating the hierarchy of eqs 64 and 65. The simplest, mean-field approximation (MFA), which is equivalent to the Hartree–Fock approximation in many-electron problems114 uses the full factorization of all normally ordered products left angle bracket[B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket and left angle bracket[B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket, left angle bracket[B with circumflex][B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket.102,128 The excitons are treated as noninteracting quasiparticles, each moving in the mean field created by the others. Equation 64 is then closed, and the other equations become redundant.103 A higher-level approximation factorizes all daggered and undaggered variables left angle bracket[B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracketand left angle bracket[B with circumflex][B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex][B with circumflex]right angle bracket, while retaining the left angle bracket[B with circumflex][B with circumflex]right angle bracket variables.79,104,105 This coherent exciton dynamics (CED) takes into account two-particle statistics but neglects transport. Equations 68 and 69 then form a complete closed set. Another possible factorization left angle bracket[B with circumflex][B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex][B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket neglects double-exciton statistical properties, while maintaining transport. This level of theory is equivalent to the Semiconductor Bloch equations (SBE) with dephasing, used for Wannier excitons.114 The response functions predicted by these various levels of approximation will be discussed in the coming sections.
The response functions, eqs 14 and 18, connect the induced polarization to the incoming laser electric fields. The third-order signals can be calculated using the Green’s function solution of the equation hierarchy, eqs 64 and 65 (eqs 68 and 69 as a special case). The induced polarization is defined as the expectation value of the dipole operator (eq 59).
7.1. Linear Response
In the linear regime, we get
equation M97
(71)
Single-exciton dynamics when the field is switched off is described by
equation M98
(72)
obtained from eq 68 with Y [equivalent]0. Propagation of excitons is described by a single-exciton Green’s function,
equation M99
(73)
whose solution is G(t) = θ(t) exp(−iht), where h is the matrix with elements hmn. The solution of eq 72 is given by equation M100.
The first-order terms in the NEE equations are obtained from
equation M101
(74)
Substituting the Green’s function solution of eq 74,
equation M102
(75)
into eq 71, and using eq 11, we get the linear response function:
equation M103
(76)
Linear response techniques are surveyed in Appendix A.
7.2. Third-Order Response
The third-order induced polarization is
equation M104
(77)
Below, we present two levels of approximation for the kI, kII, and kIII techniques: the MFA and the coherent exciton dynamics (CED) limit.
7.2.1. Mean-Field Approximation
At this level of theory, we make the factorizations left angle bracket[B with circumflex][B with circumflex]right angle bracket = left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket. Equation 68 then becomes
equation M105
(78)
and the polarization (eq 77) reduces to
equation M106
(79)
The first-order variable is given by eq 75. The second-order variables vanish, and to third order from eq 78 we get
equation M107
(80)
Here and later, summations over repeating indices are implied. Below we present simplified expressions obtained by setting equation M108. The complete third-order polarization obtained by substituting eqs 75 and 80 into eq 79 is given in Appendix D.
The kI response function is generated by the VB(1)*B(1)B(1) term in eq 80 as follows: pulse 1 creates B(1)*, while the second and third pulses generate B(1) (in any order). This gives
equation M109
(81)
For the kII technique, we similarly have
equation M110
(82)
and for kIII,
equation M111
(83)
In these expressions, the kj pulse interacts with chromophore nj, and j = 4 denotes the signal generated at chromophore n4. We thus keep track of all interaction wavevectors and polarizations.
Corresponding frequency-domain expressions of the signals are given as a special case of eqs 127129 and will be discussed in section 9.
7.2.2. The Coherent Exciton Dynamics (CED) Limit
When the system–bath coupling is neglected, the system remains in a pure state and can be described by a wave function at all times. Equations 68 and 69 then adequately describe the coherent short time dynamics.
A Green’s function for the Y variable is now required for solution of the NEE. The Y variables describe two-exciton dynamics. The corresponding Green’s function defined by Ymn(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig46.jpg(t)Ykl(0) satisfies
equation M112
(84)
The solution, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(t) = θ(t) exp(−ih(Y)t), requires the diagonalization of a tetradic matrix h(Y). To avoid this diagonalization, we calculate this Green’s function using the Bethe–Salpeter equation, as described in Appendix E. In that appendix, we further introduce the exciton scattering matrix Γ, which allows one to interpret double-exciton resonances in terms of the exciton scattering.
By expanding the variables in powers of the fields, the first-order variable is given by eq 75 and the second-order variable
equation M113
(85)
where we have further used the relation An external file that holds a picture, illustration, etc.
Object name is nihms212366ig47.jpg [equivalent] An external file that holds a picture, illustration, etc.
Object name is nihms212366ig48.jpg, which results from boson symmetry. Equation 68 then gives
equation M114
(86)
The polarization (eq 77) is given by
equation M115
(87)
For equation M116, we get for the kI response function
equation M117
(88)
We recast this result using the exciton scattering matrix (see Appendix E). Substituting eq 324 in eq 88 with equation M118, we get
equation M119
(89)
Similarly, we obtain for kII
equation M120
(90)
and for kIII:
equation M121
(91)
Γn4n3, n2n1 is the exciton scattering matrix in two-exciton space. It is a tetradic matrix, whose indices n1 and n2 represent incoming excitons and n3 and n4 are outgoing excitons. The scattering is induced by the coupling V. The scattering matrix is given by
equation M122
(92)
where
equation M123
(93)
Closed frequency-domain expressions of the signals corresponding to eqs 8991 are obtained from eqs 127129 by neglecting transport and setting equation M124, as discussed in section 9.
The physical picture for exciton dynamics emerging from this level of theory is very different from the MFA: excitons scatter by their mutual interactions as demonstrated by the diagrams in Figure 8 during the time interval ττ′ > 0. This scattering changes the resonant frequencies; thus, the correct double-exciton resonances are recovered.
Figure 8
Figure 8
Feynman diagrams representing the quasiparticle dynamics in the various techniques for μ(2) = 0 in the coherent limit (eqs 8991). Each interaction with the field is displayed by a solid dot. Time propagation is from the bottom up. The (more ...)
So far, we have expressed the third-order response in terms of exciton Green’s functions. These must be calculated in the presence of a bath to include dephasing and population relaxation. In section 4, we have introduced the bath phenomenologically using the Lindblad master equation. We now derive and extend these results using a microscopic model of the phonon bath.107
8.1. Reduced Dynamics of Excitons Coupled to a Bath
As in section 6, we consider an aggregate made out of N three-level molecules and described by the Frenkel exciton Hamiltonian:
equation M125
(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,
equation M126
(95)
where α runs over the bath oscillators and equation M127 are Bose creation (annihilation) operators that satisfy equation M128. Note that
equation M129
(96)
equation M130
(97)
where Qα and Pα are the coordinate and momentum of bath oscillator α with frequency wα. The system–bath coupling is assumed to be linear in bath coordinates
equation M131
(98)
where the coupling parameters dmn represent bath-induced fluctuations of the transition energies (m = n) and the intermolecular couplings (mn). For strongly coupled systems (J [dbl greater-than sign] d), both contribute to energy-level and coupling fluctuations when transforming to the eigenstate basis.
Because of the coupling with the bath, the factorization left angle bracket [B with circumflex] [B with circumflex]right angle bracket into left angle bracket[B with circumflex]right angle bracketleft angle bracket[B with circumflex]right angle bracket in eqs 64 and 65 now no longer holds, and exciton dynamics for the first-to-third orders in the field depends on two additional variables, equation M132 and equation M133. 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
equation M134
(99)
equation M135
(100)
equation M136
(101)
equation M137
(102)
where we have used An external file that holds a picture, illustration, etc.
Object name is nihms212366ig49.jpg = hnnδmmhmmδnn, which is derived from the single-exciton Liouville operator An external file that holds a picture, illustration, etc.
Object name is nihms212366ig50.jpg; “[multiply sign in circle]” denotes the direct product of two vectors: [A [multiply sign in circle]B]mn = AmBn. The relaxation rate matrices are calculated using eqs 355 and 359. The rate matrix K is given by eq 362 (since ρnm = Nmn, the corresponding relaxation matrices satisfy Kmn,kl = Knm,lk).
Some useful approximations may be used for the relaxation matrices. equation M138 will be represented approximately as equation M139. For certain parameter regimes, the full Redfield relaxation operator K may lead to unphysical (negative and larger than 1) probabilities.111,140 In the next section, we convert the operator in the eigenstate basis into the Lindblad’s form, which guarantees a physically acceptable solution.
8.2. Redfield Equations in the Secular Approximation: The Lindblad Form
The secular approximation described in section 4 is widely used since it guarantees a physical behavior of the propagated density matrix.111,120,141145 We shall transform eq 101 into the single-exciton basis, i.e., the eigenstates of the h matrix, e = εeψe. This gives (the field is turned off)
equation M140
(103)
Assuming that equation M141 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 equation M142 with equation M143. For the Redfield equation in the eigenstate basis in the secular approximation, populations are decoupled from coherences. Furthermore, each coherence satisfies its own equation and is decoupled from the rest. Populations satisfy the Pauli master equation (eq 26) with population transport rates, Kee,ee. Closed expressions for the relaxation rates and for equation M144 in the eigenstate representation in terms of the bath spectral density are given in Appendix G. This approximation for the relaxation operator, also known as the Davies procedure,146 guarantees that the density matrix remains positive-definite. Some examples of nonphysical behavior of the full Redfield equation are given in refs 111, 122, 140, and 142. Similar difficulties have been found for the Bloch equations, which describe the nuclear induction of a spin that interacts with a magnetic field.147 The problems associated with the Redfield equations are a consequence of the second-order perturbation theory and the Markovian approximation.140 Various schemes have been suggested to remedy this, by introducing short-time corrections to these equations.140,148,149 These were shown to give physical results for specific model systems and some ranges of parameters, where the Markovian Redfield equation breaks down. Unlike the secular approximation, these schemes do not guarantee the positivity of the density matrix. The secular approximation, when applied to the Redfield equation, brings it into the Lindblad form.111,120 The stochastic Liouville equation described in section 12 allows a more general description of the dynamics that does not rely on the secular approximation and guarantees a physically acceptable behavior.
In the eigenstate basis, the Lindblad equation (eq 25) can be written as
equation M145
(104)
where we represent α by a pair of indices c and d. This will be convenient for connecting with the Redfield equation. equation M146 where V(c,d) has been projected into the one-exciton eigenstate, Ψe = ψemΦem with Φem and ψlm defined in sections 6 and 9.
We next present the set of V’s that reproduce the Redfield relaxation rates of Appendix G. We define two types of V matrices: the first contains only a single nonzero off-diagonal element
equation M147
(105)
with cd. Since the rank of the Ne2e1 matrix is N, there are N (N – 1) such V(c, d) matrices (these correspond to all off-diagonal elements of Ne2e1). The remaining N matrices are associated with the diagonal elements (c = d):
equation M148
(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
equation M149
(107)
To obtain the b’s, one needs to solve the equation
equation M150
(108)
where [gamma with circumflex]e1,e2 is the pure-dephasing rate (as introduced in eq 29),
equation M151
(109)
Multiple solutions for the b’s exist since eq 108 contains N( N – 1) equations with N2 unknowns (the N equations for e1 = e2 in eq 108 are satisfied for any choice of b’s). One particular solution is obtained from the rate expressions of Appendix G:
equation M152
(110)
Finally, the frequency, equation M153 in eq 104 is related to ω e1e2 in eq 103 by a level-shift,
equation M154
(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.
8.3. Transport of Localized Excitons; Förster Resonant Energy Transfer (FRET)
The full Redfield Green’s function derived by a second-order expansion of the rates in the system–bath couplings, couples populations and coherences. It is invariant to the basis set used, i.e., a unitary transformation between the localized and delocalized basis sets does not affect the dynamics.144 However, because of the different way populations and coherences are treated, the secular approximation is basis-set dependent. The delocalized eigenstate basis represents strongly coupled chromophores. For weakly coupled chromophores, the aggregate eigenstates are essentially localized on individual chromophores and the real-space representation then constitutes the natural basis set, where the Redfield equations become
equation M155
(112)
Making the secular approximation in this basis, the Redfield relaxation operator becomes equation M156 where K(F) is known as the Förster exciton-transfer rate. Note that equation M157 (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
equation M158
(113)
where FD(ω) is the fluorescence spectrum of the donor molecule, εA(ω) is the absorption spectrum of the acceptor molecule, and VDA is the dipole–dipole coupling between the molecules:
equation M159
(114)
Here, μ are the transition dipoles and RDA is the vector connecting the donor and acceptor molecules. Note that the exciton transfer rate includes the diagonal energy fluctuations. In practical applications, one uses the experimental absorption εA(ω) and emission FD(ω) spectra in the Förster formula. The strong ~R−6 dependence of K(F) is used in fluorescence resonant energy transfer (FRET) studies, to probe molecular structural fluctuations and energy transport.151,152
For small RDA, the molecular electron densities begin to overlap and energy transfer can proceed via a different, Dexter, mechanism involving the double exchange of an electron and a hole between the donor and acceptor molecules. This mechanism, which results in an exponential, ~e−αR; α > 0, dependence of the rate on the distance is particularly important for triplet excitons where the optical transition is forbidden and the Förster rate vanishes.153,154
We now derive closed QP expressions for the response functions in the exciton eigenstate basis. Exciton transport takes place during t2. During the delay periods t1 and t3, we only include dephasing and neglect transport. This is usually justified since population relaxation times (ps to ns) are typically longer than the interband dephasing times (100 fs).
The procedure for calculating the response function for the full set of eqs 99 and 100 is the same as in section 7.2.2: the dynamical variables are expanded to third order in the field using the necessary Green’s functions. These include the Green’s functions for the N and the Z variables: N(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig51.jpg(t)N(0) and Z(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg(t)Z(0) where N = left angle bracket [B with circumflex] [B with circumflex]right angle bracket and Z =left angle bracket [B with circumflex] [B with circumflex][B with circumflex]right angle bracket. An external file that holds a picture, illustration, etc.
Object name is nihms212366ig51.jpg is a tetradic matrix, while An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg is sextic. Note that the N variable corresponds to the transpose of the density matrix in the single-exciton space, i.e., Nmn [equivalent] ρnm. The same holds for the Green’s functions in the real-space representation (compare eq 31), so that equation M160, and so we will use only An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg. The response functions obtained by solving eqs 99 and 100 with the polarization equation (eq 77) are given by
equation M161
(115)
equation M162
(116)
equation M163
(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 t2 interval represent An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg (when both point up) and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg (when pointing in opposite directions). Three lines (triple propagation) within the t3 interval represent An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg. A horizontal dashed line stands for the V matrix.
Figure 9
Figure 9
Same as Figure 8 but including exciton transport (eqs 115117). Single solid line represents G(t) propagation, double line is either An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(when both lines point to the same direction) or An external file that holds a picture, illustration, etc.
Object name is nihms212366ig56.jpg (when the lines point to opposite directions), and triple line (more ...)
RkI is depicted in the left diagram. The interval between the first and second interactions is described by the single-exciton variable left angle bracket [B with circumflex]right angle bracket, whose propagation oscillates with interband frequencies according to G(t1). During the second interval, the system is characterized by left angle bracket [B with circumflex] [B with circumflex]right angle bracket and its propagation given by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg(t2). This oscillates at low, intraband frequencies associated with differences between single-exciton energies on different sites and their couplings. The [B with circumflex] [B with circumflex][B with circumflex]right angle bracket variables generated during the third interval propagate according to An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg(t3t′), which again oscillate with high interband frequencies. The evolution during t2 is strongly influenced by the bath since the low, intraband oscillation frequencies are close to the bath frequencies. Therefore, the interplay of incoherent transport and dynamics will be most important during this interval as described by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg(t2). The high-frequency evolution during t1 and t3 is only weakly influenced by the bath (unless it contains resonances at optical frequencies).
Approximate factorizations will be used next to simplify the calculation of the Green’s functions An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg: by neglecting exciton population relaxation in the third interval, we can set equation M164.68 An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg is expressed in terms of the exciton-scattering matrix Γ (Appendix E). The secular approximation, eq 31, will be used for An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg so that exciton eigenstate populations will be decoupled from their coherences.
The Green’s function expression (eqs 115117) assumes a simpler form in the one-exciton eigenstate basis, e. The one-exciton Green’s functions become
equation M165
(118)
where the single-exciton dephasing, γe, is derived in Appendix G. In the frequency domain, we have
equation M166
(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
equation M167
(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,
equation M168
(121)
and is related to the single-exciton Green’s functions by a convolution:
equation M169
(122)
To describe the dynamics of the N variables, we transform their Green’s function to the exciton basis:
equation M170
(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 kI, we obtain from eq 115
equation M171
(124)
The time-domain response functions equation M172 can be calculated by using the inverse transform (eq 317). We note that the Fourier transforms with respect to t1 and t2 only involve single Green’s functions and are trivial. We can thus derive simple expressions for mixed representation such as equation M173. When incoherent exciton transport is neglected, we can set equation M174 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
equation M175
(125)
and
equation M176
(126)
Combining eqs 40 and 41 with the response functions (eqs 124126) and using the results of section 5, we obtain our final QP expressions for the multidimensional signals:
equation M177
(127)
equation M178
(128)
and
equation M179
(129)
These expressions contain fewer terms than their supermolecule counterparts (eqs 4250) and allow one to make the approximations discussed above. The CED limit is obtained by setting equation M180, which includes dephasing but no transport. The MFA is recovered by neglecting transport and using the MFA scattering matrix, Γ(MFA)( ω) = V, as discussed in Appendix E. Equations 127129 can represent nonbosonic systems by using the scattering matrix given in Appendix H. These expressions can be readily applied for infinite periodic structures; see Appendix I.
The numerical integration of eqs 40 and 18 with the response functions of eqs 124126 will be required when the optical pulses do overlap temporally or when the density matrix dynamics is nonexponential (as in section 11). Simulated finite-pulse signals for the FMO light-harvesting complex are presented in the next section.
Photosynthetic light-harvesting complexes found in biological membranes (see Figure 1) constitute an important class of chromophore aggregates. Photosynthesis starts with the absorption of a solar photon by one of the light-harvesting pigments, followed by transfer of the excitation energy to the reaction center where charge separation is initiated.4,150,157,158 This triggers a series of electron-transfer reactions, the redox chain, where ADP is eventually converted into ATP, which stores chemical free energy.150 The primary step, the absorption of a photon, takes place in light-harvesting antennae containing assemblies of pigment molecules that absorb light in a broad window of the solar spectrum. The transfer of this excitation toward the reaction center occurs with a remarkably high (98%) quantum yield.150 The active pigments are various types of chlorophylls (Chl) and bacterio-chlorophylls (BChl, BChla, and BChlb), which absorb light in the 600–700 nm and at longer wavelengths,159,160 and carotenoids, which absorb light at higher frequencies, in the 450–500 nm range.159,161 Most energy transferred to the reaction center is absorbed by the Chl and BChl molecules, but there is clear evidence for energy transfer from the carotenoids to the Chl molecules.162164 The carotenoids also act as photoprotective agents by thermally dissipating excess energy that could otherwise generate harmful photooxidative intermediates.165,166
Optical properties of several complexes, light-harvesting complex 1 (LH1) and 2 (LH2) and photosystem 1 (PSI) and 2 (PSII), have been extensively studied using linear and nonlinear optical techniques, revealing ultrafast exciton dynamics and transport pathways.4,37,158,167 Understanding the interaction of these complexes with light is crucial for unraveling their function and developing new generations of artificial light-harvesting and photonic devices and poses major theoretical and experimental challenges. Multidimensional spectroscopy provides an invaluable tool for unraveling and quantifying the photophysics of the photosynthetic apparatus.40,41,81 These techniques can pinpoint couplings between pigment molecules70 and can directly probe excitation energy transfer timescales.
The FMO complex168,169 in green sulfur bacteria is the first light-harvesting system whose X-ray structure has been determined (Figure 1). The FMO protein is a trimer made of identical subunits, each containing 7 BChl pigments and no carotenoids. This system has been extensively studied by 1D techniques such as absorption and linear and circular dichroism.170,171 The spectra were simulated using the exciton model (eq 58) whose parameters (site energies, εm, and interactions, Jmn) were fitted to experiment.40,170,171 Subsequent 2D optical spectroscopy40,41,81 revealed that the excitation energy is transferred toward the reaction center using specific pathways; the energy does not simply cascade stepwise down the energy ladder.40 By recording 77 K 2D spectra vs t2, Engel et al.41 observed 600 fs quantum beating in both the shape and intensity of the peaks, indicating long-lived coherences. The role of exciton delocalization and coherence in energy transport in photosynthesis is an intriguing issue: if they survive at room temperature, they could allow for fast transport and improve the efficiency of energy funneling to the reaction center. Pulse sequences that can distinguish between coherent quantum oscillations and incoherent energy dissipation during t2 will be presented in section 13.81
Multidimensional spectroscopy has also been applied to other light-harvesting systems: Lee et al. measured coherence dynamics in the bacterial reaction center,172 and Zigmantas et al. generated a 2D signal of photosynthetic antenna LH3 (27 chromophores) in purple bacteria.173 The reaction center of purple bacteria absorbs light in two well-separated bands at 800 and 850 nm. It is surrounded by a complex of BChl molecules forming a ringlike core light-harvesting antennae LH1.174,175 The membrane LH1 is surrounded by peripheral light-harvesting antennae LH2 and LH3. LH2 exists in different sizes. In Rps. acidophila, it is made of two concentric rings containing 18 and 9 BChl molecules, respectively, while in Rps. molischianum, the two rings have 16 and 8 BChl molecules, respectively.174 The two rings are connected by carotenoids.174,176 More details about the structure and the interactions of these light-harvesting complexes are given in refs 4, 174, and 177. Many theoretical studies of 1D spectroscopy of this system have been carried out. Schröder et al.178,179 tested various methods for calculating the linear absorption lineshapes, Jang and Silbey180,181 developed line-shape theories for single-molecule absorption, and Dahlbom et al.182 calculated transition-absorption of the B850 ring of LH2. The pump–probe and photon echo signals as well as superradiance were covered in our group.34,130,183185 Pullerits et al. have estimated the average delocalization length of exciton in LH2 from the simulated and experimental pump–probe spectra.186,187 They showed that the LH2 excitons are delocalized on 4 chromophores on average.
Photosystem I (PSI) is the largest photosynthetic system with known structure containing 96 BChl and 22 β-carotenoid molecules.16 Six of the BChl molecules form the electron transfer chain and constitute the reaction center. Ten “red” BChls absorb light at lower energies than the others, including the reaction center. Their precise role in energy transport is under debate.188 Photosystem II (PSII) is made of a light-harvesting antenna, LHCII, which contains 8 BChla, 6 BChlb, and 4 β-carotenoids,189 and its reaction center.190 Numerous experimental and theoretical 1D spectroscopic studies have been carried out on PSI37,38,188,191 and PSII.192195 These include absorption, linear and circular dichroism, and pump–probe.
The peak linewidths for the kI technique at t2 = 0 reveal the homogeneous and inhomonegeous dephasing rates.70,196 Variation of the cross-peak pattern with t2 allows one to monitor pigment interactions and excitation energy transfer timescales.40,81 During t2, the system’s density matrix evolves either in the first excited-state manifold, e, where population transport and coherent evolution take place, or in the ground state, as seen from the Feynman diagrams (Figure 4). Various quantum master equations for the evolution of the reduced density matrix have been applied to study energy transport in photosynthetic systems.120,178 These include the Redfield theory,120 a modified Redfield theory,185,197 and some non-Markovian variants.178,198 Other studies199 have demonstrated the differences between the Bloch model and the Redfield theory for the description of relaxation processes in LH2. Short-time, low-temperature, dynamical features of the pump–probe spectra of LH2 were explained using a simplified polaron model.200 Exciton annihilation has been shown to strongly affect nonlinear optical signals.201203
Various levels of theory developed in the earlier sections will be applied below to simulate 2D signals of the FMO complex of Chlorobium tepidum. Its chromophore configuration is sketched in Figure 10. The BChl molecules of the FMO monomer are modeled as 7 coupled chromophores. The Hamiltonian of eq 58 is used with An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg = 0 and Δ → ∞ (two-level chromophores). εm, Jmn, and μn are taken from ref 81. The overdamped Brownian oscillator spectral density (eq 372) with timescale Λ −1 = 100 fs and interaction strength λ = 35 cm−1 is used for the bath in the Markovian limit. The linear absorption, calculated using eq 269 shown in Figure 11 has five peaks labeled ae. The d and e bands contain unresolved doublets.
Figure 10
Figure 10
FMO complex structure of Chlorobium tepidum (left) and its simulated linear absorption (right), taken from ref 81 (eq 158).
Figure 11
Figure 11
Linear absorption spectrum for FMO of Chlorobium tepidum. The 5 identifiable peaks a–e, where a is the lowest-energy absorption peak and e is the highest. The power spectra of the three excitation pulses (Wb, Wd, and We) used to calculate the (more ...)
10.1. Simulated Signals at the CED and MFA Levels
We start by neglecting transport and using the quasiparticle approach (eqs 127129) with the scattering matrix eq 387 (MFA) or eq 92 (CED). We first assume broad pulse bandwidths, setting the pulse spectral envelopes to 1.
The MFA and CED SkI signals at t2 = 0 are compared in Figure 12. Since the two approaches are identical in the single-exciton space, the diagonal section, which shows the single excitons, is similar and reflects absorption peaks. However, their amplitudes are different.
Figure 12
Figure 12
Im equation M369 signal (imaginary part) in the MFA and CED for FMO at t2 = 0 ps.
The large differences appear mostly in off-diagonal regions, where double-exciton properties show up. In the MFA, the two excitons evolve independently, and their correlations are neglected. Peak positions, which reflect the double-exciton energies, are then given by the combinations of single-exciton energies εe + εe′. The CED scattering matrix, in contrast, induces additional shifts of the resonant positions coming from exciton–exciton interactions.
A similar comparison is shown in Figures 13 and and1414 for the double quantum coherence SkIII signal.83,205,206 It projects the double excitons directly on the Ω2 axis, making it very sensitive to the fine features of the double-exciton manifold. The MFA provides an approximation to the double-exciton manifold by neglecting exciton–exciton correlation effects. The two peak patterns are mostly different in the higher-energy region (Ω2 >24 800 cm−1), where the MFA signal is featureless but the CED shows well-resolved peaks.
Figure 13
Figure 13
equation M370 signal in the MFA and CED for FMO at t3 = 1 fs and t3 = 200 fs.
Figure 14
Figure 14
equation M371 signal in the MFA and CED for FMO at t1 = 0 fs.
The remainder of this section shows simulated kI signals that go beyond the CED to include transport during t2 at the secular Redfield theory level.
10.2. Dissecting Energy-Transfer Pathways by Using Finite-Pulse Bandwidths
The 2D signals were calculated using eq 127 and include transport during t2. We used Gaussian pulse envelopes,
equation M181
(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.
Figure 15
Figure 15
2D signal equation M595 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 the d band. The black line marks the diagonal. (more ...)
Experiment A (top row) uses four broad pulses, whose power spectra are given by Wb in Figure 11 and cover the entire exciton band. At t2 = 0 (left column of Figure 15), the kI signal is dominated by two diagonal peaks corresponding to bands b and d. The cross-peaks are weak. As t2 is increased, relaxation toward lower energies is observed; a higher-energy photon creates an excited state that decays toward lower-energy excited states, and the system emits at lower energies. At t2 = 1 ps, cross-peaks start to appear below the diagonal. The cross-peak at (−Ω1 = d, Ω3 = b) shows that some excitation energy of the d band has been transported to b. The right panel (t2 = 5 ps) shows additional cross-peaks. The (d, b) and (d, a) peaks indicate that energy from the d band has migrated to both the b and a bands. The absence of a diagonal (d, d) peak indicates that most of the initial exciton energy of the d band has decayed to lower-energy bands. The (b, a) peak is induced by energy transfer from b to a. Energy transport to the lowest excited state is completed in a few picoseconds. All energy-transfer pathways contribute to the signal; however, some are clearly dominant.
In Figure 15B, the first two pulses are narrow (We in Figure 11; σ1/2 = 50 cm−1 bandwidth) and selectively excite band e, and the last two pulses are broad (Wb). At t2 = 0 ps, the pulses select peaks such that −Ω1 lies in band e. The dominant diagonal peaks at (b, b) and (d, d) fall out of the pulse bandwidth, making weak cross-peaks more visible. At t2 = 1 ps, the population created in band e has been transferred to b, c, and d (as seen from the larger (e, d), (e, c), and (e, b) cross-peak intensities) but not to the lower a band (no (e, a) cross-peak). After 5 ps, we see four cross-peaks indicating that all lower bands are populated, but the population of the higher energy states, especially d, is very small (weaker (e, d) and (e, c) cross-peaks). It is difficult to observe this pathway with broadband pulses (Figure 15A), since it is an order of magnitude weaker than the dominant one.
In Figure 15C, the first two pulses are narrow (Wd in Figure 11) and selectively excite the d band, but the last two pulses are broad. At t2 = 0 ps, only the diagonal (d, d) peak is observed. As in Figure 15B, no significant excitation has been transferred to the lowest-energy band in 1 ps (no (d, a) cross-peak). What is different here is that only two bands are populated during the entire energy-transfer process (shown as the cross-peaks (d, b) and (d, a), the cross-peak at (d, c) never appears). The c band is not populated at t2 = 5 ps even though its energy is lower than the original excitation, suggesting that the c and d bands correspond to localized excitations residing on different (and weakly coupled) sites. Similar conclusions were reached by Brixner et al.40 who argued, by analyzing the exciton eigenvectors of a model system similar to ours, that the energy does not simply relax stepwise down the energy ladder.
10.3. Dissecting the Response Functions in Real Space
The third-order response function, whether calculated using the supermolecule or the quasiparticle approach, contains a product of four dipole moments associated with the various interactions with the fields. In order to analyze the NEE equation M182 response function, eq 124, we shall transform the dipoles back to real space using μe = Σnψe,n μn with the single-exciton eigenstates ψe,n defined in section 9. We further single out the n4 sum and write
equation M183
(131)
with
equation M184
(132)
equation M185 gives the signal generated at site n4. By varying t2, we can further find the sites into which excitation energy flows, thus characterizing the population transport pathways and timescales in real space. Monitoring the energy transport in real-space is of particular interest in photosynthetic antennae, where the energy must eventually be funneled into specific sites (the reaction center).
A similar dissection procedure could be applied by singling out the third, second, or first interactions with the field (i.e., restricting the sums in eq 131 over n3, n2, and n1, respectively) and calculating the contributions of each site to the various interactions. The linear absorption signal can be similarly dissected by singling out the second interaction (see eq 269):
equation M186
(133)
It gives the contribution of site n2. The dissected absorption spectrum of FMO displayed in Figure 16 shows that the lowest-energy excited state (band a) is localized on chromophore 3, band b is on chromophores 4 and 7, and band d is delocalized over chromophores 5 and 6.
Figure 16
Figure 16
Absorption spectrum for FMO dissected into its seven chromophores by the last interaction with the field (eq 133). The contribution from each chromophore is given by a different color.
The equation M187 spectrum displayed in Figure 17 is dominated by chromophores 1, 4, 5, 6, and 7. As can also be seen from the absorption dissection, the (b, b) peak is delocalized over chromophores 4 and 7. The (d, d) peak is delocalized over sites 5 and 6, while (c, c) is localized on chromophore 1. For t2 = 10 ps (Figure 18), the signal is generated at chromophores 3, 4, and 7. This is to be expected since excitation energy flows in a downhill path. Band a is localized on chromophore 3. Hence, the cross-peaks (d, a) and (b, a) are dominated by RkI,3. The second lowest energy band, b, is still populated at t2 = 10 ps. Hence, RkI,4 and RkI,7 still contribute to the response function since the b band is delocalized across sites 4 and 7. These contributions contain the two peaks at (b, b) (the b band is still populated) and (d, b) (exciton transfer from d to b).
Figure 17
Figure 17
equation M372 signal (eq 132) for the FMO light-harvesting complex at t2 = 0 ps. (Left column) The imaginary part of the response function; (Right column) the absolute value. The signal is dominated by chromophores 4, 5, 6, and 7 (Figure 10), which also dominate (more ...)
Figure 18
Figure 18
Same as Figure 17 but at t2 = 10 ps. The signal is dominated by chromophores 3 and 4 (and a small contribution of 7), onto which the lowest excited states (3) and the second lowest excited states (4 and 7) are delocalized.
This dissection analysis provides similar information to the finite pulse envelope techniques discussed above. The difference is that here we select the chromophores, whereas the finite envelopes select the exciton eigenstates. The two are close when the eigenstates are localized. One can achieve a similar dissection using a broad band for the first three pulses and a narrow band for the last pulse.
This dissection analysis provides information that goes beyond the eigenvalues and the eigenvectors: it can be used to highlight specific energy-transfer pathways and timescales (as was done in section 10.2 using finite pulses) and find the dephasing rates of coherences between chromophores. Exciton delocalization can be characterized by the participation ratio,182,207,208
equation M188
(134)
Le varies between 1 and N for an N-chromophore aggregate and is a measure of the effective number of chromophores contributing to the exciton e. For FMO with the Hamiltonian parameters of ref 81, the participation ratios for the seven exciton states are (in increasing order of energy): L1 = 1.3, L2 = 2.3, L3 = 1.6, L4 = 1.8, L5 = 2.4, L6 = 2.8, and L7 = 1.6. As expected, the lowest-energy exciton is localized on a single site (chromophore number 3). The other excitons are more delocalized, especially excitons 2, 5, and 6, which are delocalized on chromophores 4 and 7 (exciton 2), 5 and 6 (exciton 5), and 5 and 6 (exciton 6). Another useful measure is the coherence length,70
equation M189
(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
equation M190
(136)
where equation M191 is the contribution of the n4th chromophore to the response function. Ld will vanish for regions dominated by a single chromophore and will grow with the degree of delocalization. In Figure 19, we display Ld at t2 = 0 and t2 = 10 ps. Obviously the diagonal peaks are localized on very few chromophores. This can be anticipated from our analysis of the linear absorption shown in Figure 16. The cross-peaks in Figure 19 for t2 = 10 ps are more delocalized because they result from population transfer facilitated by exciton delocalization (c.f., eq 370).
Figure 19
Figure 19
(Left column) Absolute value of the response function, equation M373 (eq 124). (Middle column) The sum of the absolute value of the dissected response function, equation M374 (eq 132). (Right column) The delocalization length, Ld3, t2, Ω1) (eq 136). The delay (more ...)
10.4. Signatures of Quartic Exciton Couplings
The quadratic excitonic couplings Jmn in eq 58 correspond to the “strong J coupling” regime in homonuclear NMR. Chromophores coupled through negative (positive) couplings form J (H) aggregates. The quartic couplings An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg are analogous to the Darling–Dennison coupling209 in infrared spectroscopy and are known as the “weak (heteronuclear) coupling” of heteronuclear NMR. A detailed comparison between NMR and multidimensional infrared and optical spectroscopy was given in refs 57 and 58. The FMO simulations presented so far assumed that each chlorophyll is a two-level chromophore, setting An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg = 0 and Δm → ∞. This anharmonic penalty prevents two excitations from residing on the same site.
Here we study the possible signatures of An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg.210 The energy of two-exciton states for excitations on sites m and n (i.e., Em,n = left angle bracketΦfmn|ĤSfmnright angle bracket, with Φ fmn defined in section 6) is
equation M192
(137)
The An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg coupling implies that double excitations on neighboring sites have lower (negative An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg; exciton attraction) or higher (positive An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg; exciton repulsion) energy compared to the sum of energies of the two single excitations. In the NEE, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg enters solely through the scattering matrix. For finite An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg, the coherent time-evolution of the Y variables will show resonances shifted by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig10.jpg couplings.
To demonstrate the possible effects of these off-diagonal quartic couplings, we have switched on the quartic coupling between the fourth and fifth chromophores in FMO. These sites are chosen because they contribute significantly to the strong b and d absorption bands (Figure 11), as can be easily seen from the dissected linear absorption spectrum (Figure 16). This is a simple illustration of possible signature of the An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg coupling, which we use as a free parameter. Figure 20 depicts the kI signals calculated with the NEE for An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig13.jpg = ±300 cm−1 and t2 = 0.
Figure 20
Figure 20
Effects of the off-diagonal quartic coupling, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg, on the kI spectra for t2 = 0. (Left column) The imaginary part of the response function; (Right column) absolute value. From top to bottom, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig13.jpg = 0, +300, −300 cm−1. Color code is the same (more ...)
Effects of the quartic couplings are most easily explained by looking at the SOS expressions for the response function of a dimer,
equation M193
(138)
where
equation M194
(139)
equation M195
(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: μ± = Ψ±,1μ1 + Ψ±,2μ2. The dimer model has a single doubly excited state, f, with energy ε+ + ε + An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg, and the transition dipole between the singly-excited and the double-excited state is μf = Ψ±,1μ2±,2μ1. The two terms of eq 140 are ESA contributions. The quartic coupling enters only in these terms.
As can be seen from eq 139, the GSB and ESE contributions to the cross-peaks interfere with the ESA. Hence, the cross-peaks at (ε+, ε) and (ε, ε+) interfere destructively. As An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg is increased, the ESA peaks are shifted along Ω3 by An external file that holds a picture, illustration, etc.
Object name is nihms212366ig11.jpg and the cross-peaks no longer cancel out, creating a doublet.
A similar splitting of the cross-peaks is expected in larger aggregates. This explains the trends seen in Figure 20. When An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = 0, the diagonal peaks (b, b) and (d, d) dominate. As An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg is turned on, the GSB and ESE contributions to the cross-peaks shift spectrally and no longer interfere with the ESA and the cross-peaks at (b, d) and (d, b) appear. Furthermore, the ESA cross-peak can be observed as weaker peaks shifted along Ω3 (these regions are circled in Figure 20). These cross-peaks appear at (b, d + An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg) and (d, b + An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg).
Signatures of the off-diagonal quartic couplings are also manifest for t2 > 0. In Figure 21, we show the kI response function for different An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg as in Figure 20 but for t2 = 10 ps. The b band is still populated even after long t2. This is evident from the large diagonal peak in the b region for An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = 0, shown in the top panel of Figure 21. Hence, the fourth site is populated since band b corresponds to an exciton localized mainly on this site. For An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = −300 cm−1, the doubly excited Φf45 state is energetically favored. The excitation energy, absorbed in the b or d band, which relaxes to b during t2, can re-emit from the d band, leading to a high-energy diagonal peak and corresponding cross-peaks. Since the doubly excited state Φf45 is blue-shifted when An external file that holds a picture, illustration, etc.
Object name is nihms212366ig12.jpg = +300 cm−1, these new features are weaker in the second row of Figure 21. In this case, the anharmonic coupling shifts the ESA contributions for the high-energy diagonal peak (d, d) and for the (b, d) cross-peak along Ω3.
Figure 21
Figure 21
Same as in Figure 20, but at t2 = 10 ps.
In section 7, we have derived quasiparticle expressions for the exciton-response functions by assuming fast fluctuations of the bath, which justifies the Markovian approximation for relaxation. This leads to a Pauli rate equation for populations and dephasing of exciton coherences. The supermolecule approach, in contrast, can describe the coupling to a bath with a broad range of fluctuation timescales. This higher-level treatment requires the explicit calculation of two-exciton eigenstates, which is only feasible in small aggregates. Fortunately, the detailed bath dynamics can be typically observed only in small systems where the electronic spectra are not too congested.
It was shown in section 3 that the resonant Liouville space pathways can be classified into ESE, ESA, and GSB types. The secular approximation for the Green’s function An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg additionally decouples coherences and populations. The response function may then be partitioned into contributions of coherences C and populations P during t2.
The resonant contributions to each technique are given by the Feynman diagrams of Figures 47. In the following expressions, we will use the same labeling for the various terms used in the diagrams: ia, ii, iiia, iva, v, and via are purely population contributions, and ib, iiib, ivb, vib, vii, and viii are purely coherence contributions.
The kI response function will be recast in the form:
equation M196
(141)
equation M197
(142)
equation M198
(143)
equation M199
(144)
equation M200
(145)
As in section 4, we define the complex transition frequencies
equation M201
(146)
where τa is the lifetime of state a. ϕ(t3, t2, t1) are multipoint line-shape functions induced by bath fluctuations. These are responsible for all other dephasing mechanisms beyond the lifetimes τ.
Similarly, we recast the kII response in the form:
equation M202
(147)
equation M203
(148)
equation M204
(149)
equation M205
(150)
equation M206
(151)
Finally, kIII is composed of ESA1 and ESA2 (both are coherent contributions since populations are not created in kIII):
equation M207
(152)
equation M208
(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 141153 and will only differ by the phase functions ϕ that will be given by eqs 166, 177, and 179. The response function expressions of section 4 are recovered by setting all ϕ = 0 and assuming that the dephasing is caused by finite lifetimes + pure dephasing in the Markovian approximation.
11.1. Coherent Exciton Dynamics; Diagonal Gaussian Fluctuations with Arbitrary Timescale
We assume that the system is described by the Hamiltonian equation (eq 94) with the system part in the form of eq 2 and the harmonic bath model given by eq 95. We partition Ĥ = Συ|υright angle bracketĤυleft angle bracketυ|, where the bath operator Ĥυ characterizes the bath-induced perturbation of the energy of system state υ,
equation M209
(154)
qυυ is a collective bath coordinate:
equation M210
(155)
qυυ induces the frequency fluctuations of state υ with respect to the reference Hamiltonian of state a:
equation M211
(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
equation M212
(157)
where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig14.jpg(τ2τ1) = left angle bracket[mu] (τ2)[mu] (τ1)right angle bracket is the two-point dipole correlation function. The angular brackets imply equilibrium statistical averaging over the bath. An external file that holds a picture, illustration, etc.
Object name is nihms212366ig14.jpg is calculated exactly for this model:
equation M213
(158)
Here, ρa is the equilibrium population of state a, and we have introduced the line-shape function
equation M214
(159)
equation M215 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,
equation M216
(160)
see eq 342 and 353. Its symmetries and other properties are summarized in Appendix F. We then have
equation M217
(161)
The cumulant expansion similarly gives for the third-order response function,
equation M218
(162)
where
equation M219
(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 ks). Expansion in the exciton eigenstates yields
equation M220
(164)
where [var phi]dcba(τ4,τ3,τ2,τ1) is a four-point line-shape function. For our harmonic bath model, we get
equation M221
(165)
where τij = τiτj. The expressions given in ref 100 were derived in the semiclassical limit, left angle bracketωab(tab(0)right angle bracket = left angle bracketωab(0)ωab(t) right angle bracket; here, we do not invoke that approximation.24 It is important to note that [var phi] does depend on the ground-state index g, even though it does not appear explicitly in the right-hand side (rhs), since all time evolutions are calculated with the ground-state reference Hamiltonian (eq 160).
By invoking the RWA, we can now select the resonant contributions in the response function, as warranted for each technique. The resulting response functions then assume the form of eqs 141153 with equation M222 and no lifetime broadening τ−1 = 0. The phase functions are given by
equation M223
(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
equation M224
(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),
equation M225
(168)
where β = (kBT)−1 and νn = (2πn)/([variant Planck's over 2pi]β) are known as the Matsubara frequencies. For negative times, we have equation M226. The second term is a low-temperature correction, given by the additional Matsubara oscillators. The third term is responsible for bath-induced transition frequency shift. In the high-temperature limit, kBT [dbl greater-than sign] [variant Planck's over 2pi]Λ, we have
equation M227
(169)
11.2. Response with Population-Transport and Intermediate-Timescale Bath Fluctuations
The results of the previous section hold for purely diagonal (frequency) fluctuations with arbitrary timescales; population relaxation caused by off-diagonal fluctuations has been neglected. The Redfield equations, on the other hand, include transport but are limited to fast fluctuations. Here, we present approximate expressions, which include exciton transport and intermediate timescale fluctuations.
Dephasing occurs during the delay times t1 and t3, whereas transport takes place during t2. Typically t1,t3 [double less-than sign]t2 (for photosynthetic complexes, population transport occurs in 500 fs-10 ps, while the coherence time scale is 100–300 fs). We assume an intermediate bath fluctuation time: fast compared to the transport but slow or comparable to the coherence-dephasing time. In this case, we can use the doorway-window representation of the response function and project the doorway and window functions onto the population/coherence blocks of the density matrix.184,185,211,212
Populations are created during t2 only in the kI and kII techniques. We therefore focus on these techniques. We first recast the response function, eq 15, in the form
equation M228
(170)
where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig15.jpg(t) [equivalent] An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg (t) An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg[rho with circumflex]0 is the doorway exciton wavepacket prepared by the first two pulses, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(t) describes its propagation and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig16.jpg (t) [equivalent] P |An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(t)An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg is the window wavepacket, which represents the detection.
Both kI and kII techniques are described by three Feynman diagrams (ESA, ESE, and GSB) as shown in Figures 4, ,5,5, and and7,7, respectively. We shall separate the response function accordingly:
equation M229
(171)
R(ESE) is given by diagrams i and iv, RGSB is given by ii and v, and RESA is given by iii and vi in Figures 4 and and55.
Expanding eq 170 in the system eigenstates, we get and
equation M230
(172)
equation M231
(173)
and
equation M232
(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 An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg, we can separate the response function (eqs 172174), into two terms:
equation M233
(175)
The coherent contribution, equation M234, only includes coherences (e1e2, e1 = e3, e2 = e4) during t2, and the entire optical process is completed before a relaxed exciton population is created. The second term ( equation M235) only includes populations (e1 = e2, e3 = e4) and represents sequential exciton transport contributions. This representation thus keeps track on how populations are created and propagated during t2. The corresponding response functions have been formally derived using projection operator techniques;185 here, we present the final results.
We shall evaluate eq 175 using the secular Redfield relaxation operator. equation M236 may then be calculated using the formalism of section 11.1. equation M237 is calculated using projection operator techniques185 by factorizing the bath averages in eqs 172174 into An external file that holds a picture, illustration, etc.
Object name is nihms212366ig16.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg, and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig15.jpg factors. This gives
equation M238
(176)
In the first, hopping, term, the doorway function De represents the population of the eth exciton created after two interactions with the radiation field. The hopping term includes the ESA diagrams (iiia for kI and via for kII) and ESE (ia for kI and iva for kII). The Green’s function equation M239 is the conditional probability for the exciton to hop from e to e′ during t2 obtained by solving the Pauli equation (eq 26 with the rate matrix in eq 362). The window function We represents the contribution of the eth exciton to the detected signal. The second term in eq 176 is the GSB contribution where, during t2, the system is in the ground state (diagram ii for kI and v for kII). Formal expressions for the doorway and the window functions are given in Appendix J.
The final expressions for the kI and kII response functions are given by eqs 141151. The relevant phase functions ϕ(ib), ϕ (iiib), ϕ (ivb), and ϕ (vib) (coherence pathways) are given in eq 166. The remaining phase functions ϕ (ia), ϕ (ii), ϕ (iiia), ϕ (iva), ϕ (v), and ϕ (via) (sequential pathways) are
equation M240
(177)
The lifetimes (eq 146) in this case are given by equation M241, and
equation M242
(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 t2, we maintain the correlated bath dynamics by treating diagonal fluctuations exactly. Second-order perturbation theory with the Markovian approximation is used for off-diagonal fluctuations to get lifetime broadenings. The populations containing pathways during t2 ignore bath-fluctuation correlations between different propagation intervals t1, t2, and t3. During t1 and t3, the diagonal fluctuations are treated exactly (by the cumulant expansion). Again, second-order Markovian perturbation theory is made for off-diagonal fluctuations. During t2, second-order Markovian perturbation theory is performed for the coupling with all bath coordinates, which results in the secular Redfield equation.
The Redfield equations assume fast bath fluctuations and yield Lorentzian line shapes. The present expressions are more general since they allow for finite bath timescale and arbitrary line shapes. However, we assumed the complete factorization of the response function into the doorway and window functions, i.e., memory is erased during the t2 period; correlations between fluctuations during intervals t1 and t3 have been neglected. In the next subsection, we present a higher-level approximation that does retain such correlations.
11.3. Combining Slow Diagonal Fluctuations with Exciton Transport
We now derive approximate expressions for the kI and kII techniques that include both slow bath fluctuations and transport.24 kIII does not involve transport and needs not be considered here. In subsection 11.2, we assumed fluctuation timescale shorter than the time intervals t1, t2, and t3. Here, the bath fluctuations are slow compared to t1, t2, and t3.
The Feynman diagrams for kI (Figure 4) represent the time evolution of the density matrix in the coherent regime where population relaxation is neglected. Figure 7 contains additional terms that represent population transfer (ee′) during t2. The population contributions, therefore, have two terms: a coherent, population-conserving term when population during t2 does not change (Figure 4) and the population-transfer term, where populations hop from state e into e′ (Figure 7). The same partitioning also applies to kII.
We assume a bath with both slow Q(S) and fast Q(F) modes. The Markovian approximation is applied to the fast modes. During t1 and t3, these cause homogeneous line broadening (diagonal fluctuations), while during t2, they induce population relaxation and dephasing (off-diagonal fluctuations) described by the Pauli master equation (eq 26). The slow bath modes are responsible for spectral diffusion during all three intervals t1, t2, and t3, which cause correlations of density matrix between all these intervals. The diagonal parts, equation M243, shift the system energies, while the off-diagonal parts, equation M244, can be eliminated by diagonalizing the system Hamiltonian with Q(S) included explicitly (the adiabatic approximation). This is justified provided the fluctuation amplitude is smaller than the intraband transitions within the one-exciton band and the Q(S) timescale is longer than dephasing times γ−1. Within this model, the Pauli master equation (eq 26) still holds during t2, but the population rate matrix K is modulated by the fluctuations Q(S). We assume a smooth bath spectral density. Q(S) fluctuations then induce weak modulations of the transport rates, which will be neglected. The population Green’s function An external file that holds a picture, illustration, etc.
Object name is nihms212366ig17.jpg(t) does not depend on Q(S) and is taken out of the bath averaging in eqs 172 and 173.24 The contribution of the population-transfer pathways to the response functions presented below is derived in Appendix K.
The complete response functions for the kI and kII techniques (eq 175) contain both coherent and sequential contributions. As in section 11.2, eqs 141151 present the final expressions for kl and kll response functions. The phase functions ϕ(ib), ϕ (iiib), ϕ (ivb), and ϕ (vib) (coherence pathways) are given in eq 166, while ϕ (ia), ϕ (ii), ϕ (iiia), ϕ (iva), ϕ (v), and ϕ (via) (sequential pathways) are
equation M245
(179)
where we define ζab = 1 − δab. As in eq 146, the lifetimes are given by equation M246. In contrast to eq 177, these expressions include correlations between the time intervals t1, t2, and t3 in both population-conserving and population-transfer diagrams. This level of theory can reproduce some signatures of population transport during t2 such as the time-dependent fluorescence Stokes shift.24
This approach, therefore, extends the DW picture of the previous subsection. The full factorization of doorway and window functions and population propagator in the previous section implies that the density matrix phase is completely lost when the population is created, even if the population does not hop to the other exciton state. The approach presented in this section retains these phases despite population hopping.
In the previous sections, bath-induced fluctuations in the exciton transport were described perturbatively. Bath variables were formally eliminated using projection operator techniques. A higher level of description is possible by explicitly including some collective bath coordinates in the supermolecule approach. This makes it possible to address a broader class of models with arbitrary bath timescales.
12.1. Stochastic Models in Liouville Space
The stochastic Liouville equations (SLE) approach provides a consistent accounting for fluctuations of arbitrary timescale by including some collective bath modes explicitly.109,213,214 The method extends the Redfield equations for exciton transport, which are limited to fast fluctuations. It further extends the CGF, which only applies to linear diagonal coupling of energies to a Gaussian bath. At high temperatures, both the CGF and Redfield approaches are obtained as limiting cases of the SLE. However, the SLE has several limitations: it neglects the system-bath entanglement, assumes an infinite bath temperature, and is computationally more expensive than the other methods. Finite temperature corrections will be discussed in section 12.4.
The stochastic model assumes that the exciton Hamiltonian (eq 57) is coupled to some classical stochastic variables σ(t), which represent the bath. The response function is calculated by averaging left angle bracket right angle bracket over the ensemble of stochastic paths
equation M247
(180)
Here, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig18.jpg = [Ĥσ(t), · · · ] is the Liouville superoperator and subscript + denotes forward time ordering. It is important to note that Ĥ0;σ now depends on time only implicitly through the variable σ.
Equation 180 may be recast in the form
equation M248
(181)
where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig19.jpg(ta,tb) is the Green’s function solution of the Liouville equation for free evolution (no field) between times ta and tb:
equation M249
(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
equation M250
(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
equation M251
(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:
equation M252
(185)
An external file that holds a picture, illustration, etc.
Object name is nihms212366ig20.jpg ρ [equivalent] [Ĥσ,ρ] thus describes both the system evolution Hamiltonian and the system–bath coupling ( equation M253 in eq 94). The system density matrix (and the response function) is finally obtained by the standard prescription of statistical mechanics: averaging over the initial and summing over the final states of the stochastic variables.
12.2. Stochastic Models for Nonlinear Spectroscopy
Since the bath dynamics is Markovian, the third-order response function may be factorized into the product of three Green’s function matrices in the joint system + bath space, each representing a single time interval, with
equation M254
(186)
with
equation M255
(187)
Here, the initial joint density matrix |ρ(0) right angle bracketright angle bracket = |ggright angle bracketright angle bracket|0right angle bracket is a direct product of the equilibrium distribution of the stochastic variables (zero right eigenvector of L(σ)) and the ground state of the system. The final summation left angle bracket0|TrS includes both tracing over the system and summation over the final bath states realized as a scalar product with the zero left eigenvector left angle bracket0| of L(σ).111
Equation 186 can be readily used to calculate the nonlinear response function; all correlations are built into the bath variables. The frequency-domain response function is obtained by a Fourier transform of the Green’s functions
equation M256
(188)
equation M257
(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:
equation M258
(190)
Equation 190 accounts for coherence transfer (e.g., gejgek) as well as exciton transport. Consequently, the Feynman diagrams (Figures 4 and and5)5) can be used, but e and f now denote the entire single-exciton and two-exciton blocks. The SLE is not limited to any particular basis.
Two types of stochastic models are commonly used in theories of spectral line shapes.215 The first, M state jump, assumes that the bath hops between M discrete states. The Hamiltonian Ĥα depends on the state α. The master equation for Pα [equivalent] P(σ = α) is
equation M259
(191)
and the SLE reads
equation M260
(192)
The zero left eigenvector of An external file that holds a picture, illustration, etc.
Object name is nihms212366ig20.jpg, which represents the final summation over discrete bath, is given by right angle bracket0| = (1, 1, ..., 1).
In the second model, σ is a continuous dimensionless stochastic variable Q undergoing an Ornstein–Uhlenbeck process and described by a Fokker–Planck equation216
equation M261
(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) = δ(QQ′) is
equation M262
(194)
The bath density approaches the equilibrium distribution Peq(Q) = G(Q,Q′; ∞) = (2π)− 1/2 eQ2/2 at long times t [dbl greater-than sign] Λ−1. Q is, thus, a Gaussian–Markovian variable. The final summation over a continuous bath variable is represented by ∫dQ. A Gaussian process is uniquely characterized by its two-point correlation function.
equation M263
(195)
All higher multipoint correlations may be calculated by using the Gaussian factorization rule:
equation M264
(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
equation M265
(197)
where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig2.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig21.jpg operate on the system Liouville space. We further assume noninteracting chromophores with purely diagonal bath fluctuations
equation M266
(198)
The solution of the stochastic Liouville equations without the field
equation M267
(199)
is
equation M268
(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.
equation M269
(201)
Equation 199 can be solved by the second-order cumulant expansion
equation M270
(202)
with the line-shape function
equation M271
(203)
Equation 203 may be also obtained from eq 161 by taking the infinite temperature limit, and the following spectral density equation M272. Comparing eqs 169 and 203, we find 2λνν′ [equivalent] (β[variant Planck's over 2pi])dννdν′ν′. At infinite temperature, the imaginary part of the line-shape function vanishes. Qualitatively new effects of the bath timescale, described by the SLE, will only appear when we go beyond eq 198 and include the coupling between excitons, add nonlinear coupling between the system and Q, or assume a non-Gaussian bath.
The third-order response functions can be obtained by using the cumulant expansion and recast in the form of eqs 162165 together with eq 203. To see that, we calculate, for instance, the ESE contribution
equation M273
(204)
Making use of eqs 201 and 203, the second-order cumulant gives
equation M274
(205)
This agrees with eqs 162165. The other contributions to the third-order response can be recovered in the same manner.
For interacting chromophores, the second-order cumulant expression is an approximation. However, the exact solution can be calculated from the SLE (eq 193) by expanding in the eigenbasis set of L(σ) (eq 193),
equation M275
(206)
where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig22.jpg are Hermite polynomials and α = 0, 1, 2, · · · . The initial (equilibrium) bath state is represented by |0 right angle bracket= (1, 0, 0, ..., 0). The final summation assumes right angle bracket0| = (1, 0, 0, ..., 0).
The bath Liouvillian is thus diagonal [L(Q)]αα= −αΛ δα α, and the Q variable is given by a tridiagonal matrix.
equation M276
(207)
Assuming linear system–bath coupling (eq 197), the stochastic Liouville equations become
equation M277
(208)
By switching to Hilbert space notation, we have
equation M278
(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 equation M279 (lowering equation M280) operators equation M281, these can be calculated iteratively in the form of a continued fraction
equation M282
(210)
and
equation M283
(211)
The Green’s function is computed by starting with the diagonal element
equation M284
(212)
For fast fluctuations |An external file that holds a picture, illustration, etc.
Object name is nihms212366ig21.jpg| [double less-than sign] [variant Planck's over 2pi]Λ, the continued fraction can be truncated at the lowest level and the stochastic variables may be eliminated; the dynamics can then be expressed in the reduced system phase space
equation M285
(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 Cij(t) can be reproduced by including several independent Gaussian–Markovian coordinates Qk (each has its relaxation rate Λk and coupling d(k)). An arbitrary C″ij is then decomposed into Lorentzians
equation M286
(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 left angle bracketqij(t)qkl(0) right angle bracket can be decomposed similarly to eq 214 by allowing finite di≠j.
12.3. Application to Excitons
We now apply the SLE to the exciton model eq 57. The Redfield equations and other quantum master equations in the reduced system Liouville space hold in the limit of short bath correlation time (Λ → ∞ in eq 193) where the bath degrees of freedom may be projected out of the explicit treatment of the SLE. To compare with the Redfield equations, we focus on the single-exciton manifold (An external file that holds a picture, illustration, etc.
Object name is nihms212366ig53.jpg block). Assuming fast stochastic dynamics, we can employ eq 213. For a general coupling of the form Ĥ1 = [variant Planck's over 2pi]Σlm dlm|elright angle bracketleft angle bracketem|, we get for the relaxation term in the eigenbasis of Ĥ0
equation M287
(215)
By making the secular approximation (see discussion below eq 103), we obtain the following expressions for exciton-transfer rates (ij) due to off-diagonal fluctuations (compare eqs 104 and 112):
equation M288
(216)
equation M289
(217)
In addition, we have dephasing rates originating from exciton transport (lifetime) and pure dephasing due to site diagonal fluctuations
equation M290
(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 model217219 (eq 192 with M = 2) has been applied220 to hydrogen-bonding fluctuations in the OD stretching mode of phenol in benzene,221 for which 2D line shapes show the timescale of creation/dissociation of the hydrogen bond. Below, we examine slow two-state-jump fluctuations in a dimer. We consider a model system with a few well-resolved absorption peaks. These may arise either from various exciton states or from different slowly inter-converting states of the bath. The absorption and the 2D SkI line shape of two noninteracting chromophores, undergoing two-state-jump spectral diffusion, are shown in Figure 22. In Figure 22 [top], the “up”, u (“down”, d) state is responsible for the two outer (inner) peaks. The 2D line shapes at short delay times t2 = 0 are plotted for the i, ii, iv, and v diagrams (Figure 22 [bottom]). In diagrams i, ii, and v, we see cross-peaks between different excitons, because one exciton can be annihilated and another can be created by two interactions with the laser field, and the aggregate can thus oscillate at a different frequency during t1 and t3 without any exciton transport during t2. However, no cross-peaks are observed at short delay times between various bath states before they can change. (For fast bath, motional narrowing will collapse the spectrum into a single peak.) In diagram iv, the exciton state in t1 and t3 interval must be the same and no cross-peaks are observed.
Figure 22
Figure 22
Linear (top) and 2D signal (bottom) for two-exciton with a two-state-jump (u, d) bath model. The signals corresponding to diagrams i, ii, iv, and v are displayed for t2 = 0, ε1u = 30Γ, ε1d = 10Γ, ε2u = −30Γ, (more ...)
Second, the SLE can account for the finite bath fluctuation timescale. For instance, in ref 221, fluctuations originating from the complexation timescale and the 2D signals revealed the kinetic rate constants. The fluctuation timescale is crucial for some dynamical effects. For excitonic systems that possess some symmetry, the ground state may be dipole coupled to only a few eigenstates, and many levels are dark. Assuming fast spectral fluctuations (Redfield theory), these levels remain dark. A slow coordinate may, however, break these selections rules, and the dark eigenstates |ψi(Q)right angle bracket may temporarily acquire some oscillator strength, resulting in additional peaks. Similar peaks can be also reproduced by introducing static disorder on the top of the Redfield equations. However, the SLE account for the fluctuation timescales and can properly describe the dynamics of the additional peaks.222 Figure 23 illustrates this for a slow Gaussian spectral diffusion in a homodimer where the absorptive line shape (the sum of kI and kII signals)
Figure 23
Figure 23
1D and 2D absorptive signal (eq 219) of a model dimer with spectral diffusion described by two Gaussian–Markovian oscillators (one coupled to one site) for short t2 = 0 (left panel) and long t2 [dbl greater-than sign] Λ−1 (right panel) times. (more ...)
equation M291
(219)
is plotted for two t2 delay times. All oscillator strength is carried by the symmetric state (upper diagonal peak in Figure 23), and the only way to observe the other (lower) level is through the coupling to slow bath coordinates (lower diagonal peak in Figure 23), which disappears on the relaxation timescale.
Third, the SLE can describe nonlinear coupling of system parameters to Gaussian fluctuations. This effect is less dramatic than the other two, yet clearly observable. For instance, the quadratic expansion of interchromophore (J) coupling with stochastic motions of dihedral angles in trialanine is hardly affected by the 1D infrared absorption line shape but was necessary for the simulation of its 2D spectra.223 Additional effects, such as relaxing the Condon approximation and allowing for fluctuations of the dipole amplitude, may be readily described by the explicit representation of the bath.
The SLE show qualitative effects that are missed by the Redfield theory and can provide a test for other approximate methods, but at a higher computational cost. The dimensionality of the joint system + bath space grows exponentially with the number of collective coordinates. It is, therefore, not practical to assign an independent bath oscillator for each chromophore in large aggregates, as required for the most general description of fluctuations with an arbitrary degree of correlations. Few collective coordinates should be identified in order to keep the problem tractable.
12.4. Microscopic Bath Dynamics beyond the Stochastic Liouville Equations
The SLE, when introduced phenomenologically, can readily describe non-Gaussian fluctuations such as multistate jumps. However, they do not provide any guidance as to how to include finite temperature effects. These may be incorporated by adopting a more microscopic level of modeling of the collective coordinates.
We consider a bath consisting of a single harmonic oscillator (P,Q) coupled to many harmonic modes and described by the Hamiltonian183,224,225
equation M292
(220)
equation M293
(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.226229
We introduce the bath spectral density, which describes the correlations of the force F = Σjcjqj, acting on the Q coordinates by the other bath modes.
equation M294
(222)
Comparing with eq 352, we get
equation M295
(223)
The noise acting on the system is described by a different spectral density C(Q)(ω), which represents the correlations of Q coordinate in the Hilbert picture:
equation M296
(224)
Assuming that the microscopic bath modes q are faster than the collective coordinate Q, taking J″ (ω) = [variant Planck's over 2pi]MΩ2ω/Λ, and approximating Ω2ω2 ≈ Ω2, the bath spectral density reduces to the overdamped Brownian oscillator form:
equation M297
(225)
At high temperatures [variant Planck's over 2pi] coth([variant Planck's over 2pi]βω/2) ≈ 2(βω)−1, we recover the Gaussian–Markovian stochastic process (eq 195) with (½)left angle bracket{Q (t), Q (0)}right angle bracket = (βMΩ2)−1 exp(− Λ|τ|), where {A, B}= AB + BA is the anticommutator. The Uhlenbeck–Ornstein process is a high-temperature limiting case of the model (eq 220).224 By rescaling (βM)1/2ΩQQ, to obtain a dimensionless unit equilibrium width, and setting MΩ2dmnQQAn external file that holds a picture, illustration, etc.
Object name is nihms212366ig21.jpg, we recover eq 197.
Finite temperature corrections to the SLE with the linear coupling model are proportional to the anticommutator {Ĥ1, ρ},109,183
equation M298
(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
equation M299
(227)
The resulting equations of motion
equation M300
(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 dmn = δmndmm (no off-diagonal coupling fluctuations). Expanding T in exciton eigenstates |emright angle bracket, it may be added to the Smoluchowski equation,230,231 yielding
equation M301
(229)
For this model, the free-energy surfaces for the various states are linearly displaced by β(dmmdnn). In both states, the Smoluchowski equation represents diffusion of the bath in a potential whose equilibrium is displaced by dmm. The displacement depends on the state of the system. For the coherence block, we have
equation M302
(230)
and the particle moves on a surface whose equilibrium displacement lies midway between those of levels |emright angle bracket and |enright angle bracket in eq 229.
Equation 228 holds at moderately high temperatures. At low temperatures, β[variant Planck's over 2pi]Ω > 1, the quantization of bath oscillators is necessary and higher-order terms in β obtained by expanding the coth ([variant Planck's over 2pi]βω/2) term in eq 161 must be included. This requires the introduction of additional bath oscillators with relaxation rates 2πn/β (Matsubara frequencies) and a corresponding width.183,109 These corrections have a similar structure as eq 228 but increase the computational cost.
The line-shape function for noninteracting chromophores linearly coupled to an overdamped Brownian oscillator (eqs 161 and 225) may be calculated exactly for arbitrary temperature using the second-order cumulant expression and is given by eq 168 with 2λνν′[equivalent] (β[variant Planck's over 2pi])dννdν′ν′. The imaginary part of this line-shape function represents the correction eq 228 and the Stokes shift.
For interacting chromophores, finite-temperature corrections establish thermal distribution of exciton densities, which makes this level of theory important for realistic modeling of many experimental systems in optical domain (such as FMO complex at 77 K).
The phenomenological viewpoint of the SLE model as a spectral random walk suggests another class of generalizations. The theory of random walks has many techniques and models, which go beyond ordinary Markovian diffusion.232 One possible extension is the continuous-time random walk (CTRW) model of spectral diffusion.233 This non-Markovian model assumes a waiting-time distribution function [An external file that holds a picture, illustration, etc.
Object name is nihms212366ig23.jpg(t)]lk for the random walk jump of the bath from state k to l. The model is tractable thanks to its renewal property; all memory is erased at the time of the jump. Rather than looking at the joint system + bath density matrix at a given time, we must follow the contributions to the density matrix from random walks that made a step to a given state at the specified time. This contains the information necessary to predict the future without reference to the past. In contrast, the full (unrestricted) density matrix sums over contributions where different time had elapsed from the last step. This is not sufficient for predicting the future jumps. The random walk paths for the third-order response may be sorted into eight groups, based on whether or not some step is made in any of the intervals t1, t2, and t3, and the contribution of each group to the response function must be calculated separately. In each group, we first focus on subintervals between the first and last steps in each applicable interval and calculate the propagator Z between these two times by solving the integral equation
equation M303
(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 equation M304 in the mth interval and the first jump at equation M305 in some subsequent (lth) interval is determined by the fixed state of bath and contributes by factor
equation M306
(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 An external file that holds a picture, illustration, etc.
Object name is nihms212366ig23.jpg. All of these contributions must be properly multiplied, averaged over initial and summed over final states, and convoluted over the possible equation M307 and equation M308.
In Figure 24, we show the 2D absorptive line shape (eq 219) of a Kubo–Anderson two-state CTRW modulation of transition frequency of a single chromophore with long-tailed WTDFs An external file that holds a picture, illustration, etc.
Object name is nihms212366ig23.jpg(t) ≈ 1/tα+1 (1 < α < 2) resulting in the significant divergences δΩα −3 at fundamental frequencies of the two bath states and algebraic cross-peak relaxation proportional to equation M309.234 The CTRW model further allows one to treat the aging of random walks with diverging average waiting time (0 < α < 1). Aging implies that the average mobility vanishes with time, because more and more paths are trapped in the long tails, resulting, e.g., in the anomalous relation between squared displacement and time left angle bracketx2 right angle bracket [proportional, variant]tα of the free diffusion. In addition, the individual paths σ do not uniformly track the space and the random walks are nonergodic.235 Clear signatures of aging can be seen in the third-order response function (with respect to the time interval from the start of the random walk to the first pulse).129 In Figure 25, the 2D absorptive line shape (eq 219) is shown to depend on the initial time t0 elapsed between the start of the random walk and the first interaction with laser pulse. The parameters correspond to fast bath limit at t0 = 0, with a single motionally narrowed peak. Diagonal peaks at the fundamental frequencies that grow with t0 correspond to immobilized particles and provide clear signatures of aging. Most notable is the coexistence of static diagonal peaks at fundamental bath frequencies and the central motional narrowing peak. This reflects the strongly inhomogenous (nonergodic) form of aging as described by the CTRW model. Measurements of the aging response would require the preparation of fresh samples before the application of each pulse sequence, making sure that the system has fully relaxed so that the age of the sample can be properly defined. 2D spectroscopy can thus provide a probe for glassy systems characterized by broad distributions of tunneling times.
Figure 24
Figure 24
Variation of the 2D absorptive signal (eq 219) for a single chromophore two-state-jump model with CTRW jump dynamics An external file that holds a picture, illustration, etc.
Object name is nihms212366ig32.jpg(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig33.jpg(t) = (θ(t))/(2π) ∫ dω eiωt[(1)/(1 + 1ω/[1 + (κα (more ...)
Figure 25
Figure 25
Aging effects in the 2D absorptive signal (eq 219). The SA signal for the nonstationary random walk model An external file that holds a picture, illustration, etc.
Object name is nihms212366ig32.jpg(t) = An external file that holds a picture, illustration, etc.
Object name is nihms212366ig33.jpg(t) = (θ(t))/(2π) ∫ dω eiωt[(1)/(1 + (κiω)α)] at t2 = 0, κ (more ...)
Chiral molecules and aggregates can be probed by pulse sequences that yield chirality-induced higher-resolution signals. The term “chirality” has been coined more than 100 years ago by Lord Kelvin. According to his definition, “any geometrical figure, or group of points is chiral if its image in a plane mirror, ideally realized, cannot be brought to coincide with itself”. Chiral systems, therefore, come in two equally probable mirror-image configurations.236 Mirror reflection with respect to the xy plane can be described by two successive symmetry operations: (i) parity P (space inversion P (x, y, z) = (−x, −y, −z)) and (ii) coordinate system rotation by π around the z-axis. Isotropic ensembles of randomly oriented molecules are invariant to an overall rotation. A racemic mixture of chiral molecules with opposite sense of chirality is invariant to parity, whereas an unequal mixture forms an isotropic chiral ensemble. The parity operation thus transforms between isotropic ensembles with an opposite sense of chirality. Aggregate chirality can arise from both the chemical structure of individual chromophores and from a chiral arrangement (even when the chromophores are nonchiral).
Since the successive application of two parity operations restores the system back into its original state, optical signals (like all system properties and physical observables) can be classified as either parity-odd (PFo =Fo) or even (PFe = Fe).237,238 The former change sign when the chirality is reversed and, therefore, must vanish for nonchiral systems and racemates (equal mixtures of chiral molecules with opposite chirality). Parity-even signals, in contrast, are not sensitive to the sense of chirality. We shall distinguish between two types of signals: chirality-induced (CI) and nonchiral (NC).
Circularly-polarized light is chiral since it has either left (L) or right (R) screw symmetry with respect to the propagation direction. The simplest chirality-induced linear optical technique, circular dichroism (CD), measures the difference between absorption of L and R circularly polarized light.93,239,240 Thanks to its high structural sensitivity, the technique has been extensively applied to structure determination using electronic transitions in the visible241,242 and the UV243245 or vibrational transitions in the IR.246248 Raman optical activity (ROA) is a closely related technique that measures the difference between resonant scattering of L and R circularly polarized light.237,249,250
CI techniques have been extended to the nonlinear regime. The chirality of liquid and solid interfaces where the surface layer of molecules has some degree of order is commonly studied using second-order techniques251 such as sum-frequency generation (SFG) and second-harmonic generation (SHG).252,253 CI tensor elements of the second-order susceptibility have been determined for specific molecular geometries.254256 Second-order techniques may also be used to study chirality in bulk samples; however, such signals are very weak due to phase-mismatch.89,238,257,258
In the bulk, third-order signal sensitivity to system chirality requires going beyond the dipole approximation. Chirality induced (CI) two-dimensional pump–probe, photon-echo, double-quatum correlation signals, have been simulated for vibrations in chiral polypeptides and electronic spectra of chromophore aggregates.79,81,82,9092,94,97,156,259263
Geometry enters optical signals through the transition dipoles and higher multipoles. All our applications so far were based on the dipole approximation for the coupling with the field. The molecule was assumed to be small compared to the wavelength of the light and can be treated as a point particle, neglecting the variation of the phase of the field across the molecule. This may not be generally justified for large aggregates. To take these phase shifts into account, we must go beyond the electric dipole approximation and include the system’s interaction with both electric and magnetic optical fields. This requires the introduction of a large number of electric, magnetic, and mixed-response functions. A more compact but equivalent description can be developed by describing the system–field interaction using the minimal coupling (pA) Hamiltonian,
equation M310
(233)
where A is the vector potential of the electromagnetic field.
equation M311
(234)
is the electric current operator, and
equation M312
(235)
is the charge density (rα and pα are the coordinate and momentum operators of electron α). The pA Hamiltonian has two coupling terms, which are linear and quadratic, respectively, in the field. For a time-domain impulsive experiment where each of the incoming laser fields interacts only once with the system and the fields do not overlap temporally, the AA term can be neglected. We thus adopt the simplified Hamiltonian of eq 1 with the system–field interaction:
equation M313
(236)
Similar to eq 16, we shall expand the vector potential in modes
equation M314
(237)
where equation M315 is the vector potential envelope corresponding to pulse j, centered at τj, with wavevector kj and carrier frequency ωj. Both electric and magnetic optical fields are represented by the vector potential since B = [nabla] × A and the Maxwell equations give E = − Å. We then have B = iΣj kj × Aj and E =iΣj ωjAj, where Aj is the vector potential of the jth pulse.
We now generalize eq 14 by introducing the nonlocal response function for the induced current,
equation M316
(238)
where tj are the delays between successive interactions with the fields and the third-order response function R(3) now relates the current density to the vector potentials. From the Maxwell equations and the charge continuity equation, it follows that the current is related to the polarization P and magnetization M vectors:
equation M317
(239)
We consider an isotropic ensemble of aggregates. We define the center of charge of the mth molecule rm and the displacement ρm [equivalent] rrm relative to the rm. The total current operator is given by a sum of operators localized around the chromophore charge centers, J(r) = Σm δ(rρmrm)Jm(ρm). Each |ρm| < l0, where l0 is the molecular size because Jm(ρm) = 0 for |ρm| ≥ l0. Performing the Fourier transform Jm(k) = ∫ dr eik(ρ), we can recast eq 238 for current amplitudes in a form similar to eq 18:
equation M318
(240)
Here, ks = u1k1 + u2k2 + u3k3 and ωs = u1ω1 + u2ω2 + u3ω3. We shall make the RWA for all interactions with the k1, k2, and k3 fields. Equation 240 is the spatially nonlocal generalization of eq 18.
The current-density operator for the m’th molecule will be written in the form
equation M319
(241)
Fourier transform with respect to r, Jm(k) = ∫ dr exp(−ikr)Jm(r), provides a direct connection with the optical signals. When the aggregate size is small (but not negligible) compared to the optical wavelength, we can expand the transition current densities to first order in the wavevectors k, thereby connecting them to the transition dipole, μ, the transition quadrupole q, and the magnetic transition dipole m for chromophore m:
equation M320
(242)
All molecular properties have been calculated with respect to the origin Rm. The kRm contribution is required for molecular aggregates, where each chromophore has its own reference point. Because of exciton delocalization, this enters as an additional contribution of the single-exciton eigenstate quadrupole (a collection of dipoles has higher multipoles: quadrupole, etc.). The total quadrupole contribution can then be represented as a sum of two terms: the transition intramolecular quadrupole amplitude of each chromophore equation M321 and a nonlocal interchromophore contribution originating from spatial distribution of chromophores of the aggregate Σn ψneRnμn, where Rn is the location (origin) of dipole μn. Thus, the total contribution of the e exciton to the quadrupole operator is
equation M322
(243)
We note that, while in the local molecular basis, we have equation M323; this is no longer the case for the delocalized excitons basis, where equation M324, 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
equation M325
(244)
The QP third-order response function in the coherent regime for the kI =−k1 + k2 + k3 (eq 89) technique similarly reads
equation M326
(245)
Here,
equation M327
(246)
and
equation M328
(247)
is the orientationally averaged product: Jj(k) is the transition current density of transition j, kj is the jth-field wavevector, and νj is its polarization. The summation over molecules describes an isotropic orientationally averaged ensemble.
The induced currents in isotropic ensembles (eq 247) can be expressed in terms of orientationally averaged products of the various multipoles:
equation M329
(248)
13.1. The Dipole Approximation
The first term in eq 248 gives the orientationally-averaged response functions in the dipole approximation (k = 0)
equation M330
(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: −i4 ωsu3ω3u2ω2u1ω1R(3) = R(3).
Orientational averaging of the transition dipoles performed using eq 424 in Appendix L requires the calculation of the vectors equation M331 in eq 429. The elements of this vector define the three linearly independent configurations of transition dipoles, equation M332, or equation M333, or equation M334. Assuming that all laser beams propagate along z, the independent polarization configurations, which survive the orientational averaging, are listed in Table 1; the labels in (ν4ν3ν2ν1) are ordered chronologically: ν1 is the polarization vector of the first pulse, ν2 is the second, etc.
Table 1
Table 1
Linearly Independent NC Third-Order Experiments for kI, kII, and kIII Techniques (There Is No Restriction on the Wavevector Configuration in This Case)
We next go beyond the dipole approximation, where the response functions become wavevector-dependent.
13.2. Beyond the Dipole Approximation; Collinear Chiral Techniques
We assume that all laser beams propagate along z and have the same carrier frequency. This configuration has been used experimentally.76 The various phase-matching signals can be separated using phase-cycling (combining various experiments defined with different phases of the fields) as is commonly done in NMR.58,264 Note that (u1, u2, u3) = (−1, 1, 1) for kI, (1, −1, 1) for kII, and (1, 1, −1) for kIII. From eq 248, we get
equation M335
(250)
where k = |k| is the wavevector amplitude and εγβα is the antisymmetric Levi–Civita tensor (εxyz = εyzx = εzxy = 1, any other permutation of indices changes sign, and εγβα = 0 when at least two indices coincide). The field polarizations are now restricted to νj = x, y, and we assumed m = im″ to be pure imaginary, which is the case for a real basis set.
Upon examination of the vector F(5) in eq 430, as was done for the dipole approximation, we find the three independent techniques listed in Table 2. These involve products of three transition dipoles and either one quadrupole or magnetic transition dipole. Since these products change sign upon parity transformation, the signals must be induced by system chirality and we denote them chirality-induced (CI).
Table 2
Table 2
Linearly Independent CI Collinear Third-Order Experiments for kI, kII, and kIII Techniques (We Assume Collinear Configuration Where All Wavevectors Are Along z)
13.3. Noncollinear Chiral Techniques
The resonant response functions in noncollinear configurations depend on both pulse polarizations and wavevectors, R (k4,k3,k2,k1; ν4,ν3,ν2,ν1). Using eq 240, we decompose all vectors into elementary components and obtain the tensor expression of the response function for an arbitrary wavevector configuration:
equation M336
(251)
Here, ν4ν3ν2ν1(k4k3k2k1) defines an elementary field polarization and wavevector configuration. Using the formalism of ref 265, we find that only few of these configurations are linearly independent and should be calculated. Within the dipole approximation, we only had three linearly independent configurations. To first order in the wavevector, the column vector equation M337 in eq 430 defines six linearly independent techniques. These are obtained by specifying equation M338, (0, 0, 0, 0, 1, 0)T, etc. The six noncollinear wavevector and polarization configurations that lead to this equation M339156 are listed in Table 3.
Table 3
Table 3
Linearly Independent CI Noncollinear Third-Order Experiments for kI, kII, and kIII Techniques with Pulse Configurations Denoted as ν4ν3ν2ν1(k4k3k2k1) and [alpha] [equivalent] −α
13.4. Chirality-Induced Signals in the FMO Complex
The simulations presented in section 10 used the quasi-particle expressions of the response function within the dipole approximation. Orientational averaging was described in section 13.1. FMO is a small complex with only 7 single-exciton states and 21 double-exciton states; the higher-level line-shape expressions of section 11 may be used in this case. Orientational averaging can be performed beyond the dipole approximation as described in sections 13.2 and 13.3.
In Figure 26, we compare NC (xxxx) and CI (xxxy) signals calculated using the doorway–window expressions of section 11.2 for different delay times t2. The xxxx (dipole approximation) is displayed in the top row. At short time (t2 ≤ 500 fs), it shows strong blue (negative) features along the diagonal (GSB and ESE contributions). Positive (green) peaks above the diagonal show ESA. Slow redistribution of peaks observed at longer picosecond timescales reflects population transport to lowest-energy state. The second row depicts the corresponding xxxy collinear CI signal. At short delay times, the fine structure of the spectrum is much more detailed than in the NC signal. The peaks along the diagonal now resemble the CD spectrum, while off-diagonal features provide information about couplings. The variation with t2 can now be followed more closely since separate well-resolved peaks appear in addition to peak shoulders seen in (xxxx). These simulations use the realistic line width.
Figure 26
Figure 26
kI signal of the FMO photosynthetic complex calculated using the CGF approach of section 11.2, with population transport for various t2 (broad band pulses with uniform An external file that holds a picture, illustration, etc.
Object name is nihms212366ig43.jpg(ω) spectral shape were assumed; the signal coincides with the response function). (more ...)
Signals specifically designed to probe coherent and dissipative dynamics81 are shown in Figure 27. The combination equation M340 selectively probes exciton-intraband coherences and their dynamics, since the contributions from populations exactly cancel out. On the other hand, the signal equation M341 highlights population dynamics: since exciton-coherence contributions from the ESE diagrams are canceled by quantum interference. Both are CI signals that capture structural information with high sensitivity. Since population contributions have been canceled, the off-diagonal peak signatures are now much better resolved and show the coherent exciton evolution. This is connected to coherent quantum effects in the system such as exciton delocalization, excitonic nature, and relaxation mechanisms.
Figure 27
Figure 27
Optimized kI signals B1 and C1 (for details, see the main text) of the FMO complex calculated using the same approach as in Figure 26. B1 shows only contributions from excitonic coherences, and C (more ...)
The C1 signal shows a much slower dynamics. At zero delay, the signal is simply the inverse of B1, but the differences rapidly develop at longer delays. The C1 peak structure is richer than B1, and specific peaks can be attributed to coherences in ESE and ESA processes by comparing both B1 and C1: the peaks that coincide in B1 and C1 are from ESA, and the ones missing in C1 are from ESE. Peaks in C1, which are absent in B1, come from populations. Note that coherences decay within 150 fs. Thus, the longer time dynamics is now solely due to populations. The C1 signal vanishes in the absence of population transfer in the ESE diagram. Thus, the peaks in C1 are particularly sensitive to the population-transport pathways (note that ground-state populations in GSB do not evolve in our model and only give static contributions).
In complex biological systems such as photosynthetic aggregates, coherent and dissipative dynamics occur on similar ultrafast timescales and line broadening leads to congested, weak cross-peaks. Global data analysis may be used to extract dynamical information out of time-resolved signals.50 The technique allows one to unravel spectra of excited species covered under the background of other signals. This method has been used for data analysis of ultrafast pump–probe (transient absorption) spectroscopy to characterize energy-transfer dynamics in light-harvesting antenna,47,266 to identify a photoprotective energy dissipation mechanism in higher plants,51 and for other applications. The interpretation is, however, model-dependent.
The active control of ultrafast processes can provide detailed information on the underlying vibrational and electronic dynamics.267270 Coherent-control pulse-shaping techniques have been used to drive quantum systems into a desired state,271275 to manipulate excitons in multidimensional spectroscopy,96,97,262,263,276282 and to control the function of biological systems.283286
A number of NMR control techniques involving simultaneous and shaped pulses, composite pulses, refocusing schemes, and effective Hamiltonians have been proposed.287 These developments are expected to find future applications in quantum information and computation. The same ideas may be used to design and simplify 2D signals of excitons.
2D spectroscopy of photosynthetic excitons provides valuable insights into coherent and dissipative excited-state dynamics.40,41,80,173,288,289 Specific energy-transfer pathways may be directly observed through the temporal evolution of cross-peaks, and their oscillations provide information about the long-lived electronic quantum coherences. These provide signatures of quantum effects in primary biological events.290 Analysis of these 2D spectra is complicated by spectral overlaps and weak cross-peaks. Specific pulse polarization sequences designed to highlight the off-diagonal spectral features were described in section 13.81,288 In spite of this progress, resolving weak congested cross-peaks and providing general methods for simplifying 2D optical spectra as in 2D NMR remain an open challenge.
Below, we demonstrate how pulse-shaping and coherent-control algorithms may be used to simplify the multidimensional spectra of photosynthetic excitons, revealing information about dynamical correlations between exciton states. In the weak-field (third-order) regime considered in this review, the shaped laser pulses do not modify the response functions of the system but allow one to selectively amplify and resolve the targeted spectral features by manipulating the constructive and destructive interferences among Liouville space pathways.
We shall represent the optical electric field in the form
equation M342
(252)
Compared to eq 16, here we have added the temporal phase function ϕjv (tτj)68,291 and introduced the pulse polarizations ν explicitly.
The electric field can be alternatively recast in the frequency domain
equation M343
(253)
where equation M344 is the Fourier transform of the complex envelope function equation M345 and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig24.jpg = [An external file that holds a picture, illustration, etc.
Object name is nihms212366ig25.jpg]*. Note that An external file that holds a picture, illustration, etc.
Object name is nihms212366ig26.jpg (ω) and [phi with tilde](ω) are not the Fourier transforms of An external file that holds a picture, illustration, etc.
Object name is nihms212366ig43.jpg(t) and ϕ (t).
The signal (eq 40) will now be recast in the form
equation M346
(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 An external file that holds a picture, illustration, etc.
Object name is nihms212366ig44.jpg(t) and the superposition of various tensor components. Later in this section we demonstrate how the pulse phases and envelopes may be manipulated in order to control the interferences among pathways and simplify the 2D signals.
14.1. Phase, Amplitude, and Polarization Pulse Shaping
The field of coherent control with laser pulse shaping has grown in many directions since the original proposals in 1980s292295 and the first experimental realizations in the early 1990s.296299 Shaped femtosecond laser pulses provide numerous “control knobs” for manipulating the outcome of light-matter interactions. Experimentally this is achieved with pulse shapers, which can independently modify the phase, amplitude, and polarization of each frequency component of a laser pulse. Various pulse-shaper designs have been introduced to create arbitrary pulse profiles over broad spectral ranges.300304
It is possible to perform single-parameter coherent-control protocols, where a pulse width, chirp, subpulse separation, or other single control knobs are varied. One can also perform free optimization control by modulating a large number of parameters simultaneously. Pulse-shaping may be performed in the time or frequency domain, by controlling either the spectral (eq 253) or temporal (eqs 252) phase and amplitude profiles. Shaping the polarization state of the electric field (degree of ellipticity and orientation angle) provides a different class of control parameters.275,305308 This technique, pioneered by Brixner, Gerber, and co-workers, opens up a new avenue of control based on manipulating vector properties of light propagation with promising applications to time-dependent chirality.309 In this method, the spectral phases of two orthogonal polarization components are modulated, and their interference leads to complex polarization-shaped field profiles. Various pulse-shaper designs provide full control of the polarization of ultrashort laser pulses by means of phase and amplitude shaping.310312
Early work had focused on the a priori design of optical pulse shapes and was limited to simple model systems and a few control parameters. In 1992, Judson and Rabitz proposed an iterative-control scheme, where the optimal pulse is found in a self-learning, adaptive loop that does not require any prior knowledge of the system.313 This paved the way for practical applications to complex systems with multiple control parameters. The approach involves the following steps: (i) generating an input trial laser pulse, (ii) applying the pulse to the sample and observing the signal, and (iii) using a genetic learning algorithm with a feedback loop to generate new pulse shapes based on the outcome of the previous trials.314 These steps are repeated iteratively until the desired target is met. Genetic algorithms perform a parallel search on an entire “population” of pulses, using payoff (cost function) information, rather than derivatives or auxiliary knowledge, and employ probabilistic, rather than deterministic, rules.
Adaptive control was demonstrated in 1997.315317 Various pulse-shaping schemes with genetic optimization were studied by Motzkus and co-workers.318 They investigated the stability and effectiveness of evolutionary algorithms, emphasizing the influence of steering parameters, number of configurations in search space, and noise. One of the most attractive applications of these closed-loop optimizations is the adaptive control of chemical reactions in gas and liquid phases, a subject that has been reviewed recently.269,319 Polarization pulse shaping has been applied to control the ionization of iodine320 and potassium dimers321 and was found to be more efficient than the shaping of linearly polarized pulses. Brixner and co-workers have also used polarization pulse shaping for the adaptive control of nano-optical fields.204 By controlling the interferences between near fields in nanostructures, they achieved subwavelength localization of electromagnetic energy, which can be used for nanoscopic ultrafast space-time-resolved spectroscopy.322 Silberberg and co-workers had used polarization pulse shaping to eliminate the nonresonant background and selectively excite closely lying Raman modes in coherent anti-Stokes Raman spectroscopy (CARS).307 The same technique was then applied to control angular-momentum states of atoms, demonstrating that it can access the three-dimensional vectorial properties of the system and excite states not accessible by linearly polarized control.323
Polarization pulse shaping has been applied to simulate the control of exciton localization in an ensemble of energetically disordered and randomly oriented FMO complexes278 and to control the exciton dynamics in FMO276 and in larger photosynthetic complexes.277 The excitation energy was localized on one of the chromophores by varying the spectral and temporal profiles of linearly polarized laser pulses. Linearly polarized pulses and simple forms of the spectral phase function have been applied toward the control of the energy flow in the light-harvesting complex LH2.286 Prokhorenko et al.283 varied the isomerization quantum yield of retinal in Bacteriorhodopsin by ± 20% using a genetic algorithm with a feedback loop and low-intensity shaped laser pulses.
Femtosecond pulse shaping has been applied to perform phase-coherent multidimensional spectroscopy.324327 Nelson et al. constructed a 2D femtosecond pulse-shaping apparatus325327 that provides control over the delay times, phases, and temporal shapes of all four interacting pulses and can be used for the coherent control of multidimensional spectroscopy. Zanni’s group demonstrated the use of a pulse shaper to generate a phase-stable collinear pulse pair that can be used to obtain 2D electronic spectra in a partially collinear configuration.324 These developments make it possible to combine coherent control with pulse shaping with 2D optical spectroscopy.
Here, we demonstrate how adaptive control can be used to simplify 2D optical spectra of complex systems and enhance desired spectral features by shaping the spectral, temporal, and polarization pulse profiles. Various pulse-shaping algorithms may be applied to improve the resolution of nonlinear optical spectra.96,97,279 Simulations demonstrate how the pump–probe spectrum of a model helical pentamer may be better resolved by using pure-phase polarization pulse shaping with genetic and iterative Fourier transform algorithms.279 The state of light was manipulated by varying the phases of two perpendicular polarization components of the pump, holding its total spectral and temporal intensity profiles fixed. Genetic and iterative Fourier transform algorithms were used to search for pulse phase profiles that optimize the ratio of the signal at two frequencies and allow for the observation of new features.
We adopt the description of polarization pulse shaping developed by Brixner et al.305,328 The time-evolution of the electric field vector within a single optical cycle around time t is represented by an ellipse (Figure 28B). A complete characterization of polarization-shaped laser pulses is provided by the quasi-3D representation obtained by specifying the temporal intensity ij(t), total phase ϕj(t), orientation of the ellipse θj(t), and ellipticity εj(t).
Figure 28
Figure 28
Phase, amplitude, and polarization pulse shaping. (A) Laser pulse sequence and time variables used in eqs 254 and 261 to describe a heterodyne-detected two-pulse 2D photon-echo experiment. (B) Polarization ellipse representation of the electric field (more ...)
We consider a pulse propagating along z (eq 252) with two polarization components ν = x and y. The field will be described using the following variables. An auxiliary angle χj(t) ε [0, π/2] represents the ratio of the amplitudes An external file that holds a picture, illustration, etc.
Object name is nihms212366ig27.jpg(t):
equation M347
(255)
A second angle δj(t) ε [-π, π] is given by the difference between the temporal phase modulations:
equation M348
(256)
[theta w/ tilde]j(t) ε [-π/4, π/4] is a third angle
equation M349
(257)
The orientation angle of the ellipse θ j(t) ε [-π/2, π/2] is then given by
equation M350
The ellipticity angle εj(t) ε [-π/4, π/4] is
equation M351
(258)
Finally, the total phase ϕj(t) of the laser-field oscillation with respect to the perihelion of the momentary light ellipse is defined as
equation M352
(259)
The instantaneous pulse frequency ωj(t) is given by the time derivative of ϕj(t)
equation M353
(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).
Figure 29
Figure 29
Pulse shaping control of excitons in the Soret band of a chiral porphyrin dimer. (A) Structure of the complex. (B) The exciton level scheme: the ground state, four single-exciton states, and six two-exciton states. (C) and (D) are the linear absorption (more ...)
We shall compare four laser pulse-shaping protocols. The first (ST||) involves spectral and temporal shaping of both x- and y-components, creating parallel linearly polarized pulses. The real and imaginary components of the electric field are independently varied. The second form (ST[perpendicular]) involves spectral and temporal shaping with perpendicular linearly polarized pulses. The third (P) involves pure-phase polarization pulse shaping constructed using the iterative Fourier transform (IFT) algorithm.279,332 The temporal envelopes An external file that holds a picture, illustration, etc.
Object name is nihms212366ig28.jpg(t) and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig29.jpg (t) are varied holding the total temporal intensity profile equation M354 as well as the spectral amplitudes of both electric field components An external file that holds a picture, illustration, etc.
Object name is nihms212366ig30.jpg(ω) and equation M355 fixed to the initial values. The temporal (ϕjv(t)) and spectral ([phi with tilde]jv(ω);) phases were calculated using the IFT algorithm. To obtain An external file that holds a picture, illustration, etc.
Object name is nihms212366ig31.jpg(t), we varied the auxiliary angle χj(t). The fourth (STP) is the most general and uses the entire parameter space involving spectral, temporal, and polarization pulse shaping.
We apply these pulse-shaping algorithms to control the amplification of weak spectral features in heterodyne-detected 2D two-pulse photon-echo (PE) signals with a Gaussian linearly polarized first pulse with wavevector k1, local oscillator ks pulse, and a second pulse k2 = k3 (Figure 28A) tuned with adaptive pulse shaping. The signal is generated in the direction ks = −k1 + 2k2 by varying t1 and t3 and holding t2 = 0. The signals were calculated by performing multiple time integrations and taking Fourier transforms with respect to the time delays between the first and second laser pulses, τ, and between the second and the local oscillator pulse, t,52
equation M356
(261)
where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig31.jpg(t) are the complex laser pulse envelopes.
The simulated spectra are displayed in Figure 29E. The reference 2D PE spectrum obtained with two linearly polarized Gaussian pulses (top row) has two major diagonal peaks, D1 and D2, and weak off-diagonal peaks, C1 and C2. The major contributions to D1 and D2 are the one-exciton states at equation M357 and equation M358, respectively (Figure 29C). The chosen control target was to amplify the off-diagonal (C1 and C2) relative to the diagonal (D1 and D2) peaks. To specify the target, we define the integrated intensity of the jth peak.
equation M359
(262)
The ratios of the integrated diagonal and cross-peaks T1 = (ID1 + ID2)/IC1 and T2 = (ID1 + ID2)/IC2 were used as cost functions and minimized using a genetic algorithm (GA). The optimal 2D PE spectra P and STP are displayed in rows 2 and 3 of Figure 29E. Peaks D1, D2, C1, and C2 are marked by circles. Arrows mark the amplified target peaks.
Complete resolution and amplification of the target peaks was only achieved using the STP (target T1) and STP′ (target T2). The resulting polarization-shaped pulses are displayed using the quasi-3D electric-field representation based on specifying the temporal intensity Ij(t), total phase ϕj(t), orientation of the ellipse θj(t), and ellipticity εj(t).305,328 In Figure 29F, we show these parameters as a set of instantaneous light ellipses at different times along the electric-field propagation axis. The ellipse sizes represent the intensity, and their shapes provide instantaneous snapshots of the polarization state. The projections represent the intensities of individual electric-field components. The variation of the total phase (chirp) is shown by color. Weak cross-peaks were enhanced compared to the intense diagonal peaks. The optimized electric fields with different polarizations select different Liouville space pathways of the tensor components of the response functions, highlighting several tensor elements with distinct features. It is clear by comparison of the four pulse-shaping strategies that modulation of the polarization profiles was necessary in order to achieve the complete resolution and amplification of the cross-peaks.97
14.2. Adaptive Optimization of Pulse-Polarization Configurations
The multiparameter coherent-control schemes described above involve elaborate pulse-shaping and characterization techniques. We now describe a different approach for simplifying 2D spectra by taking linear combinations of various tensor components of the response function.
It was shown in section 13 that the third-order 2D signals depend on various tensor components of the response function Rks4ν3ν2ν1.79 The nonchiral (NC) tensor components (dipole approximation) include ν4ν3ν2ν1 [equivalent] xxyy, xyxy, and xyyx (note that xxxx = xxyy + xyxy + xyyx). We label them equation M360. The nine linearly independent chirality-induced (CI) kI = −k1 + k2 + k3 response function tensor components are labeled equation M361. These were defined in Tables 13 in section 13.
We have constructed the following superposition of linearly independent response function tensor components:
equation M362
(263)
where n = 3 for the NC and n = 9 the CI techniques. The complex coefficients equation M363 were varied using a genetic algorithm in order to optimize the ratios of the amplitudes of selected peaks in the porphyrin dimer described above (Figure 29).263 The interference of CI tensor components was used to suppress the strong diagonal peaks and amplify weak cross-peaks (Figure 30). The initial nonchiral signals at delay time t2 = 0 were simulated with a homogeneous broadening γ = 500 cm−1. Only two diagonal peaks, P1 and P5, and two cross-peaks, P2 and P4, were resolved. In the chirality-induced T1T9 spectra, the remaining four cross-peaks P3, P6, P7, and P8 were only partially resolved. Using a genetic algorithm, we found optimal linear combinations of CI tensor components that amplify these peaks and reveal information about couplings of the x- and y-polarized transitions. This information was not available from the absorption spectra or the initial 2D NC spectra. The Soret-band in porphyrin aggregates is the strongest and is connected with interesting photophysical processes such as S2-exciton emission and dynamic intensity borrowing.333 The optimized cross-peaks provide additional information that is not revealed by superpositions of the nonchiral tensor components.
Figure 30
Figure 30
Coherent control of 2D signals of a porphyrin dimer by a linear combination of 2D CI tensor components. Shown are absolute magnitude at delay time t2 = 0: four nonchiral tensor components xxxx, xxyy = xyxy, xyyx; nine chirality-induced tensor components (more ...)
We next demonstrate how coherent control may be applied to manipulate exciton-transport pathways262 in the FMO complex as revealed by 2D experiments. The simulated xxxx spectrum shown in Figure 31 (left column) fits the experimental data.40,80 Thin lines on both axes in the 2D spectra mark the exciton 1–7 positions in order of increasing energy. The cross-peaks (2,4) and (2,5) corresponding to exciton states 4 and 5 on the Ω1-axis, and state 2 on the Ω3-axis, overlap with the diagonal peaks of the same sign, forming an L-shape signal.
Figure 31
Figure 31
Coherent control of energy transfer in FMO by adaptive optimization of a linear combination of 2D CI tensor components. Absolute magnitude of 2D signals: nonchiral xxxx and optimal 2D CI (fast) and (slow) linear combinations corresponding to the optimizations (more ...)
Using coherent control combined with a genetic learning algorithm, different energy-transfer channels were separated by optimizing a linear superposition of CI tensor components. The control targets were aimed at manipulating various peaks in these 2D spectra by selecting the ratios of the integrated cross-peaks in the absolute values of the 2D signals at long delay time t2 = 5 ps. Two cross-peaks (1,5) and (1,7) contributing to predicted fast and slow pathways40,80 are marked in Figure 31 by red and green circles, respectively. These cross-peaks represent energy transfer from higher to lower energy states, and grow as a function of t2. The (1,7) cross-peak was not resolved in the xxxx spectrum.
In the optimization of the fast energy-transfer pathway, the ratio of (1,5) to (1,7) is maximized. The red circle and arrow mark the cross-peaks (1,5) and (1,2), respectively. The optimal coefficients obtained from the optimizations at t2 = 5 ps were then used to obtain the 2D spectra for other t2 delay times. The exciton states 1, 2, and 5 are involved in the fast pathway, and the corresponding cross-peaks are enhanced. For the slow pathway, the ratio of (1,7) to (1,5) is maximized. Green circle and arrow indicate cross-peaks (1,7) and (1,3), respectively. As a result, both pathways are separately optimized and well-resolved. In addition, the fast pathway shows a node (a region where the signal vanishes) between the arrow and the red circle, corresponding to the elimination of the (1,3) cross-peak. Similarly, in the optimization of the slow pathway, we see nodes corresponding to the elimination of the (1,2), (1,4), and (1,5) cross-peaks of the fast pathway. Optimal laser-pulse polarization configurations obtained using coherent control and genetic algorithms enhance the chirality-induced spectral features, revealing a specific energy-transfer pathway that was not resolved in the xxxx spectra. These methods may be applied to other 2D techniques. They have also been applied, for example, in vibrational 2D spectroscopy to resolve the structure of amyloid fibrils.334
In conclusion, pulse shaping combined with coherent control can increase the sensitivity of multidimensional spectroscopy of molecular aggregates by unraveling weak cross-peaks from otherwise congested spectra.
Femtosecond nonlinear optical spectroscopy had gradually evolved from simple pump–probe (transient absorption) applications in the early 1980s to the elaborate coherent multidimensional techniques surveyed here. The field has been driven by both ultrafast pulse-laser source technology (down to attoseconds)268 and progress in the theoretical modeling of complex molecular assemblies with dissipation.
The response-function formalism provides a unified picture of nonlinear optical techniques in systems driven by moderately strong optical fields. The signals are given as convolutions of response functions to various orders in the incoming fields, with the field envelopes. The response functions carry all necessary molecular information for calculating the optical signals.
We have surveyed the key aspects of the theory, focusing on coherent third-order heterodyne-detected signals. This holographic phase-sensitive detection gives the signal field itself (both amplitude and phase). Early third-order optical measurements used homodyne detection, which only gives the signal intensity. For completeness, we briefly describe these signals using the response-function formalism.
The simplest third-order technique, pump–probe, is performed with two optical pulses with variable delay time τ. The first pulse, the pump, perturbs the system, and the second, the probe, detects the time evolution of this perturbation. The technique can be performed in various configurations: the pump and probe can have same/different carrier frequencies, have same/different bandwidths, and be in completely different optical regions ranging from THz to the UV. Instead of using the response function, this signal can be alternatively described as an induced absorption/emission events from a nonstationary state. Using the third-order response function when the pump and the probe are temporally well-separated, the pump is associated with the first two interactions, while the probe induces the third interaction and the self-heterodyned detection. Since the first two interactions occur within the pulse envelope and are indistinguishable, the signal is the sum of the kI and kII techniques, which only differ by the order of the first two interactions (Figures 4, ,5).5). The delay time t1 is controlled by the pump-field, while t3 is controlled by the probe. The delay between pump and probe pulses then coincides with the time variable t2. The pump–probe signal is thus given by
equation M364
(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 k1 and k2 are applied simultaneously; the probe comes after a delay τ and the homodyne signal is measured in the k1k2 + kpr direction. Usually ωpu and ωpr are tuned to the same optical transition and the 1D time-domain signal is recorded versus τ. The pulse bandwidth selects the relevant spectral region. This homodyne signal is given by
equation M365
(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 kI signal as a function of two delay times t1 and t2.
equation M366
(266)
This signal is usually displayed as a family of 1D plots vs t1, for various values of t2. For an ensemble of two-level systems, the signal along t1 grows for short t1, reaches a maximum at t1 = τM, and eventually decays. The value of τM and its dependence on t2 reflects the bath correlation function (eq 343). Thus, the dependence τM(t2) (the 3PEPS signal) reveals direct information on the system–bath couplings.4244,48,335,336 This interpretation is limited to a single two-level chromophore.
We have introduced two types of response functions: one (R) connecting the induced polarization with the incoming electric fields (eq 14), and the other (R) relates the induced current with the vector potential of the field (eq 238). The relation between them can be established using the Maxwell equations. In the dipole approximation where the signals calculated using both approaches must coin-cide, we find that R → (ω1ω2ω3ω4)−1R, where ωj is the carrier frequency of the jth laser field. We derived R in the dipole approximation, ignoring the magnetic dipole and the electric quadrupole contributions. This is adequate for molecular complexes of nonchiral molecules. The response function R may be extended to include the magnetic corrections. However, it is much more convenient to calculate R instead.
We have used the Frenkel-exciton model in the boson representation. Two types of parameters are responsible for the nonlinear optical signals. The quartic potential Umnkl controls the exciton scattering when the field is turned off, while equation M367 represents anharmonic system–field interaction. Hard-core bosons are a special case of this model, as shown in Appendix C. They have Ummmm = ∞. This model developed here for electronic excitations also applies to vibrational excitations. We then model the excitations as intrinsically boson-like with a finite Ummmm and equation M368 (soft-core bosons). This approach has been applied to coupled vibrations in polypeptides.67,100,337
The supermolecule and quasiparticle scattering are two conceptually different approaches for calculating the response functions of aggregates. The sources of nonlinearities in these descriptions look very different. In the supermolecular representation, the double-exciton states are included explicitly, and interband transitions do not follow a harmonic oscillator model. In the quasiparticle representation, the final expressions involve only single-exciton states. Double-exciton resonances are recovered from the frequency-dependent exciton-scattering matrix. The relation and equivalence of the two representations can be established when supermolecule double-exciton states are expanded in the single-exciton product basis for an isolated system. Differences remain since different approximations are made for the system–bath interaction.
In this review, the quasiparticle representation has been used in conjunction with simple Lorentzian line shapes in the Markovian limit. More elaborate line-shape theories could be developed (equivalent to the CGF or SLE); however, this is usually not necessary in large aggregates where the lineshapes are often controlled by inhomogeneities and spectral congestions. Each peak is then composed of many overlapping transitions. The detailed microscopic description of bath dynamics is crucial for small systems (~10 or less coupled oscillators), where the line shapes are better resolved.
The marriage of coherent-control pulse-shaping techniques with 2D spectroscopy offers numerous future possibilities. The rich information obtained through snapshots of electronic dynamics can then be further refined by manipulating the interference of various Liouville space pathways. These techniques may also be used to shift and isolate specific resonances, facilitating the assignment of electronic transitions.
Chirality-induced signals may lead to useful applications of multidimensional coherent spectroscopy to structure determination. X-ray crystallography and NMR determine structures of molecular complexes with atomic resolution. Optical spectroscopy can be used as a complimentary technique for further structure refinement.
Electronic nonlinear spectroscopy is most valuable for tracking dynamical events in molecules and their complexes. Unlike traditional 1D transient-absorption and transient-grating techniques, in multidimensional coherent techniques, a large amount of information is gathered in a single shot. Numerical data processing and analysis techniques need to be developed to help extract system information (e.g., Hamiltonian, bath, structural, and relaxation parameters) from the measured signals.
The various methods and techniques reviewed here are based on a generic exciton and supermolecule approaches within the RWA (sections 2, 5, 6) and may be applied to other frequency regimes: from THz to X-ray. 2D techniques, originally developed in NMR,53 were extended during the past decade to IR laser spectroscopy.60,338343 Simulation protocols based on molecular dynamics (MD) strategies have been reviewed recently in ref 67. These applications use the same formalism; however, the vibrations are bosons. Electrostatic parametrization of vibrational Hamiltonian parameters is currently a useful technique for including fluctuations, simulated by MD trajectories. Developing electrostatic maps for electronic Hamiltonians is necessary for affordable modeling of solvation effects.
The theoretical apparatus presented in this review has been implemented in the computer software package “SPECTRON”, which has been used for simulating coherent nonlinear optical responses in IR,79,91,92,94,95,155,344,345 visible,81,97 and UV82,156 regions. The code includes the super-molecule and quasiparticle approaches as well as the spectral lineshape-broadening models described in sections 2, 4, 9, and 11.
Figure 33
Figure 33
Feynman diagrams of the quasiparticle dynamics in kII technique. The left-most diagram is for the μ(0) = 0 case, while others include μ(0) ≠ 0. The symbols are the same as in Figure 32.
Acknowledgments
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.
Biographies
An external file that holds a picture, illustration, etc.
Object name is nihms212366b1.gif Object name is nihms212366b1.gif
Darius Abramavicius was born in 1974 in Alytus, Lithuania. He received his M.S. degree from Vilnius University in 1998 and his Ph.D. degree under the supervision of Prof. Leonas Valkunas from the Institute of Physics at Vytautas Magnus University, Lithuania, in 2002. He then joined the group of Professor Shaul Mukamel at University of Rochester (New York, U.S.A.). He moved to the University of California Irvine in 2003. During 2005–2006, he was a temporary lecturer in Theoretical Physics Department of Vilnius University in Lithuania. Currently he is a visiting scientist at the University of California, Irvine, in the group of Prof. S. Mukamel. His research interests are in the field of optical and dynamical properties of molecular complexes.
An external file that holds a picture, illustration, etc.
Object name is nihms212366b2.gif Object name is nihms212366b2.gif
Benoit Palmieri was born in 1978 in Montréal, Canada. He received his Ph.D. degree under the direction of Professor David Ronis from the chemistry department of McGill University in 2007. The main subject of his Ph.D. thesis was the theoretical study of diffusion of gas or liquid molecules inside channeled structures (zeolites). He joined the group of Professor Shaul Mukamel in 2007, where he worked on the nonlinear optical response of photosynthetic systems and on the importance of interband excitonic couplings. He is now working at the Canadian Space Agency under the direction of Dr. Marcus Dejmek, developing theoretical models of crystal growth.
An external file that holds a picture, illustration, etc.
Object name is nihms212366b3.gif Object name is nihms212366b3.gif
Dmitri V. Voronine was born in 1979 in Moscow, Russia, and is currently a Postdoctoral Researcher in the group of Prof. Tobias Brixner at the Institut für Physikalische Chemie, Universität Würzburg, Germany. He received his Ph.D. in 2004 from the Center for Photochemical Sciences, Bowling Green State University, under the supervision of Prof. Michael A. J. Rodgers, and did postdoctoral work with Prof. Shaul Mukamel at the Department of Chemistry, University of California, Irvine. He has been the recipient of the McMaster Fellowship (BGSU) and the RCCM Excellent Young Investigator award (U. Würzburg). His research interests in theoretical and experimental physical chemistry and nanobiophysics include applying the methods of coherent control with laser pulse shaping to femtosecond 2D spectroscopy of complex systems such as biomolecules, photosynthetic complexes, chiral and fractal molecular aggregates, porphyrins and carotenoids; and ultrafast nano-optics of metal nanostructures for space-time-resolved spectroscopy with nanometer spatial and femtosecond temporal resolution.
An external file that holds a picture, illustration, etc.
Object name is nihms212366b4.gif Object name is nihms212366b4.gif
František Šanda was born in 1976 in Domažlice, Czech Republic. He received his Ph.D. in Theoretical Physics (supervisor Vladislav Čápek) in 2003 from Charles University, Prague. From 2004 to 2005, he had a postdoctoral appointment at the University of California, Irvine, in the group of Prof. Shaul Mukamel. He is an Assistant Professor at Institute of Physics, Faculty of Mathematics and Physics, Charles University in Prague since 2006. His research interest focuses on stochastic modeling of dynamic processes in Liouville space, with emphasis on applications in spectroscopy.
An external file that holds a picture, illustration, etc.
Object name is nihms212366b5.gif Object name is nihms212366b5.gif
Shaul Mukamel, currently the Chancellor Professor of Chemistry at the University of California, Irvine, received his Ph.D. in 1976 from Tel Aviv University. Following postdoctoral appointments at MIT and the University of California, Berkeley, he has held faculty positions at Rice University, the Weizmann Institute, and the University of Rochester. He is the recipient of the Sloan, Dreyfus, Guggenheim, and Alexander von Humboldt Senior Scientist award and the Lippincort Award and is a fellow of the American Physical Society and the Optical Society of America. His interests focus on developing computational techniques for the design of novel ultrafast laser pulse sequences for probing electronic and vibrational dynamics in molecules. Biophysical applications include folding and dynamical fluctuations in proteins and DNA, hydrogen bonding, long-range electron and energy transfer in photosynthetic complexes, and signatures of chirality. Other areas are attosecond x-ray spectroscopy, excitons in semiconductor nanostructures, many-body effects in quantum optics, and photon and electron statistics in single-molecule spectroscopy and open molecular junctions. He is the author of over 600 publications in scientific journals and the textbook, Principles of Nonlinear Optical Spectroscopy (Oxford University Press), 1995.
17. Appendix A: Survey of Linear Optical Techniques
Linear optical signals are generated by a single interaction with the incident optical field. Three types of linear signals are most widely used: absorption, where the transmission of the incident beam is measured; linear dichroism, the difference between two polarized absorptions (longitudinal, L, and transverse, T) in oriented samples; and circular dichroism, the difference between the absorptions of two circularly polarized beams (left, L, and right, R).
Using the response function formalism, the linear absorption may be interpreted as follows. The incident beam interacts once with the molecule and induces the linear polarization. This polarization then creates the outgoing optical field, which interferes with the incoming beam. This is a self-heterodyne detection. The detector measures the net optical field, which is attenuated by the interference. The mathematical formulation is very nicely described in the book by Loudon.117
The absorption of the field is defined by following the field intensity in the medium, which follows
equation M387
(267)
The absorption coefficient κa is given by
equation M388
(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
equation M389
(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:
equation M390
(270)
Here, we introduced the transformation tensor from the molecular (μα) to the laboratory frame (μi), μi = ΣαTiαμα We assume that the molecules are oriented along z and take the molecular z-axis to coincide with the laboratory z-axis. Orientational averaging in the xy plane then amounts to rotating the molecular xy plane around the z-axis. This can be easily calculated since, in this case, only five elements of the transformation tensor are nonzero: Txx = cos [var phi], Txy = −sin [var phi], Tyx = -sin [var phi], Tyy = cos [var phi], and Tzz = 1, and the orientational averaging amounts to an integration over the rotation angle [var phi]:
equation M391
(271)
Thus, when μ2 = μ1 =μ, we obtain the following possible field configurations. For ν2 = ν1= z:
equation M392
(272)
and for ν2 = ν1 = x (the same holds for y):
equation M393
(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
equation M394
(274)
where κL is the longitudinal absorption given by eq 269 with the transition amplitude eq 272:
equation M395
(275)
whereas κT is the transverse absorption (eq 269) together with eq 273:
equation M396
(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:
equation M397
(277)
where “− (+)” sign corresponds to R(L) circular polarization.
Using the response function, we get the induced polarization generated by the R/L beam:
equation M398
(278)
The circular dichroism signal is commonly defined in terms of the ellipticity angle,
equation M399
(279)
where ER/L are the transmitted field amplitudes of the R and L beams given by eq 282. For very small angles, we use tan(θ) ≈ θ and get
equation M400
(280)
Instead of ellipticity angle, the molar ellipticity angle and the rotational strength are often used.93
18. Appendix B: Homodyne-Detected versus Heterodyne-Detected Coherent Signals
The induced polarization serves as a source for the signal field, which is eventually detected. Several detection schemes are possible. Neglecting the magnetic response, the electric field and polarization satisfy the Maxwell equation:
equation M401
(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:
equation M402
(282)
where n(ωs) is the dielectric constant at frequency ωs, l is the interaction length, and Δk = ksu1k1u2k2u3k3 is the wavevector mismatch between the polarization and the optical field; sinc(x) = sin(x)/x. Thus, the strongest amplitude occurs at the phase-matching condition Δk =0. The other important conclusion is that the signal field is directly proportional to the polarization amplitude, but is π/2 out of phase because of the i factor.
Two detection schemes for the signal are most commonly employed. These are known as homodyne and heterodyne detection. In homodyne detection (HOD), the detector is placed at the phase-matching direction and measures the integrated intensity of the signal field, or counts photons. The field intensity is proportional to the absolute value of the signal field. Thus, the homodyne signal is given by
equation M403
(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
equation M404
(284)
In heterodyne detection (HET), an additional laser pulse (local oscillator, EH) is applied in the direction of the signal field, ES. It interferes with the signal, and the resulting field is then detected by homodyne detector. We thus get
equation M405
(285)
The local oscillator (second term) can be readily subtracted. The |Es(t)|2 term is typically weak since Es [double less-than sign] EH, and thus, the HET approach detects the signal amplitude
equation M406
(286)
This measurement can be repeated with varying phase of EH to reveal the induced polarization Ps(t) using eq 282.
19. Appendix C: Bosonization Schemes for Frenkel Excitons
Transformations between fermion and boson representations widely used in many-body problems may be employed for describing excitations in semiconductors346353 and molecular aggregates.92,136,137,354,355
Consider molecular aggregates made of multilevel molecules. Molecule m has a set of transitions a from its ground state. Each transition acts as a two-level system, where the molecule can be either in the excited state a or the ground state. The molecule is now viewed as a collection of two-level systems. Within the space of a single excitation and double excitations, an aggregate made of such molecules is described by the Frenkel-exciton Hamiltonian
equation M407
(287)
where m, n label different molecules and equation M408 is a bosonic ( equation M409) exciton creation operator on molecule m of transition a (this can be brought into the form of eq 57 by relabeling pairs ma, nb, kc and ld as single indices). Vma,nb is the exciton resonant coupling term, and equation M410 is the two-boson interaction potential. When W [double less-than sign] V (weak scattering regime), the excitations are soft-core bosons, where two particles have a small penalty (anharmonicity), when occupying the same state in the molecule. The requirement that two excitations cannot reside on the same molecule is satisfied by shifting the unphysical state energies out of the physical window of interest. In the Hamiltonian equation (eq 287), this is done by setting equation M411. The electronic excitations now become hard-core bosons. Note that, by varying equation M412, we can continuously turn soft-core into hard-core bosons.
As a simple example, consider two-level molecules. The bosons do not scatter when they are on different molecules. To describe this model, we set equation M413 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,
equation M414
(288)
where all other elements vanish. Note that it involves the inversion of N × N matrix equation M415 (all other elements are neglected). These are given by
equation M416
(289)
The response is given by eqs 127129 using the transformation to the single-exciton basis with eq 120. In some cases (for instance, in the MFA), the response functions diverge. The scattering matrix of Appendix H can then be used.
An alternative bosonization scheme is possible in some cases.137 Consider a single anharmonic oscillator described by the Hamiltonian:
equation M417
(290)
where the eigenstate energy ladder consists of only d energy levels so that left angle bracket [B with circumflex] [B with circumflex]right angle bracket = 0, 1, •••, d – 1. In this case, [B with circumflex] are not boson operators; their commutation relations are as follows:
equation M418
(291)
The system thus behaves as harmonic only for low excitation levels: equation M419 when j < d; the higher-energy states do not exist.
This system can de described in terms of boson operators [b, b] = 1 using the following transformation,
equation M420
(292)
and
equation M421
(293)
where Z = exp(− b b).
This transformation can be extended for a system of coupled truncated oscillators described by the Hamiltonian (the zero-point energy (1/2 in eq 290) is ignored),
equation M422
(294)
where the operators B [B with circumflex] m refer to the mth molecule and satisfy the commutation rules:
equation M423
(295)
The bosonization transformation is then similar to eqs 292 and 293, but each operator has an additional index m. The resulting new Hamiltonian will be given in terms of various products of equation M424 and bm
As an example, consider a molecular aggregate of two-level molecules. The commutation relations (d = 2 in eq 295, hard-core bosons) are136,353,354
equation M425
(296)
In third-order spectroscopy, it is sufficient to consider up to three-particle operators; all terms involving four or more creation and annihilation operators contribute only to higher orders in the field and can be neglected. It is, therefore, sufficient to consider the Frenkel-exciton Hamiltonian of the form
equation M426
(297)
For the dipole operator, we have
equation M427
(298)
The transformation eqs 292 and 293 within the single and double excitons reduce to
equation M428
(299)
equation M429
(300)
The first term in the transformation is the single-exciton term, while the second term describes the Pauli exclusion. By substituting eqs 299 and 300 into eq 297, we get the transformed Hamiltonian,
equation M430
(301)
where umn,kl =hml(δmnδnk + δnkδkl) + Umn,kl. For the dipole operator, we similarly get
equation M431
(302)
We thus find equation M432 in this case (see eq 59 in equation M433).
20. Appendix D: Contributions of Dipole-Operator Nonlinearity to the Third-Order Response Function
In this section, we generalize eqs 115117 to include μ(2) contributions to the response function. In the NEE, this corresponds to several new sources to the third-order polarization, whose contributions are visualized in the diagrams shown in Figures 3234. The left-most diagram in each figure represents the contributions coming solely from μ as given in the main text (eqs 115117). Various configurations of interactions with μ(2) are shown in the other diagrams.
Figure 32
Figure 32
Feynman diagrams for quasiparticle dynamics in the kI technique. The left-most diagram is for μ(0) = 0 case, while others include μ(0) ≠0. Each interaction with the field is displayed by a solid dot. Time propagation is from the (more ...)
Figure 34
Figure 34
Feynman diagrams of the quasiparticle dynamics in kIII technique. The left-most diagram is for the μ(0) = 0 case, while others include μ(0) ≠ 0. The symbols are the same as in Figure 32.
The rules for these diagrams are as follows: vertical lines represent exciton propagators. As in the main text, a single line represents G, a double line stands for An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpg(both lines pointing up) or An external file that holds a picture, illustration, etc.
Object name is nihms212366ig8.jpg (lines pointing up and down), whereas a triple line represents An external file that holds a picture, illustration, etc.
Object name is nihms212366ig52.jpg. Solid dots represent μinteractions with the field, while red-circled dots represent μ(2) interactions. Since μ(2) is a local chromophore property in real space, all Green’s functions must converge into one point where the interaction takes place. The horizontal dashed line represents V. It also introduces an integral over the t3 interval.
Using these rules, we obtain the following additional contributions to eq 115 for the response function for kI:
equation M434
(303)
equation M435
(304)
equation M436
(305)
equation M437
(306)
Note that, for two-level molecules when equation M438, eq 303 cancels some contributions from eq 115 and eq 304 cancels eq 305.
For kII additional to eq 116, we similarly have
equation M439
(307)
equation M440
(308)
equation M441
(309)
equation M442
(310)
For kIII, we have the following additional contributions to eq 117:
equation M443
(311)
equation M444
(312)
equation M445
(313)
equation M446
(314)
equation M447
(315)
21. Appendix E: Coherent Single- And Two-Exciton Green’s Functions: The Bethe Salpeter Equation
In this section, we describe the coherent exciton dynamics in the absence of the bath and the field. The single-exciton dynamics is described by the Green’s function G(t) by definition B(t) = G(t)B(0) and is a solution of eq 73. We introduce the frequency-domain Green’s functions:
equation M448
(316)
with the inverse transform
equation M449
(317)
This gives, for single-exciton dynamics,
equation M450
(318)
where η→ +0 guarantees causality. The double-exciton dynamics is described by Green’s function An external file that holds a picture, illustration, etc.
Object name is nihms212366ig4.jpggiven by eq 84. In the frequency domain, we then get
equation M451
(319)
We next describe the two-exciton propagation in terms of single-exciton scattering. For this, we use h(Y) = h(0) + V, which gives
equation M452
(320)
and we define a Green’s function with respect to h(0)
equation M453
(321)
The two-exciton Green’s function can be written as
equation M454
(322)
Using the operator identity Â−1 = [B with circumflex] −1 + [B with circumflex] −1([B with circumflex] –  )  −1 and taking  = ωh(Y) + and [B with circumflex] = ωh(0) + , we obtain [B with circumflex] –  = VandarriveattheDysonequation:68,79,104,108,326
equation M455
(323)
We next define the boson-scattering matrix Γ(ω) by the equation
equation M456
(324)
Combining eqs 323 and 324 gives
equation M457
(325)
Γ, like V, is short-range, making it computationally tractable. Substituting eq 324 into eq 323 and transforming back to the time domain, we obtain the Bethe Salpeter equation
equation M458
(326)
where the first term, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig54.jpg (t), is the free-exciton propagator, with matrix elements equation M459. The second term represents the contribution of exciton scattering. We note that the scattering matrix has different units than the Green’s function: in the frequency domain, G has dimensions of time, while Γ has dimensions of angular frequency; in the time domain, G(t) is dimensionless, while Γ(t) has units of frequency squared.
Equation 326 recasts the two-exciton evolution in terms of the free quasi-particle dynamics and their scattering in real space. A doubly excited state in this representation is described by two quasi-particles residing on the same site. Thus, the quasiparticle picture applies even for a single chromophore (N= 1): two particles are then restricted to be on the same site, and their scattering appears as the diagonal anharmonicity, Δm.
For two-level chromophores, two excitons cannot reside on the same site. These are, therefore, known as hard-core bosons. By assuming that U couplings vanish, the hard-core boson scattering matrix is obtained from the boson-scattering matrix in the limit Umm,mm → ∞. The scattering matrix is then given by equation M460 with all other elements vanishing. More complicated expression for the scattering matrix (eq 369) can be obtained for hard-core bosons without performing bosonization as in Appendix H (we then also have equation M461 coming from Pauli blocking).
Note that, in the MFA limit for the response functions in eqs 8183, the response can be recast in the form of eqs 8991 using the scattering matrix, equation M462 or in the frequency domain,
equation M463
(327)
Comparing this with eq 323 shows that the MFA scattering matrix corresponds to the first-order correction of the two-exciton Green’s function (eq 323) in the scattering potential V. In the frequency domain, ΓMFA(ω) = V, and we can write
equation M464
(328)
22. Appendix F: Real Space Description of Exciton Dephasing and Transport
In this section, we derive the dephasing matrices of the NEE given by eqs 99102. For this purpose, we consider the case where the optical field is switched off.
The system and the bath are described by eqs 99102. The exciton dynamics is described by the Heisenberg equation of motion as in section 6. However, because of the bath, the equations for the exciton variables will be coupled with the bath variables. The bath contributions are calculated perturbatively to second order in the system–bath coupling, where averages of bath and exciton variables can be factorized as equation M465.107 For a harmonic bath at thermal equilibrium, we have
equation M466
(329)
where
equation M467
(330)
is the mean boson occupation number at temperature T and β = (kBT)−1. Using this factorization, we arrive at the following NEE:
equation M468
(331)
equation M469
(332)
equation M470
(333)
equation M471
(334)
We start with the single-exciton dynamics B(1) described by eq 331 and neglecting the Z variable. The bath-related term in the equation of single-exciton variables is as follows:
equation M472
(335)
The Heisenberg equation of motion for bath-assisted variables equation M473 and left angle bracket[B with circumflex](1)âαright angle bracket reads
equation M474
(336)
equation M475
(337)
A formal solution of these equations gives the following:
equation M476
(338)
equation M477
(339)
and where G(c)(t) = θ(t) exp(−iht), where h is the matrix and hmn is the zero-order Green’s function corresponding to eq 331 in the absence of the bath. Substituting these expressions into eq 335, we get
equation M478
(340)
k(B) is the relaxation memory kernel of the B variable:
equation M479
(341)
and
equation M480
(342)
is the spectral density of the coupling with the bath.
This definition is intimately related to the correlation function of bath coordinates Q defined by
equation M481
(343)
Using the symmetry
equation M482
(344)
we split it into real (which is even) and imaginary (odd) parts
equation M483
(345)
and get
equation M484
(346)
and
equation M485
(347)
Here, {A, B} = AB + BA and [A, B] = AB – BA.
The Fourier transform of the correlation function
equation M486
(348)
is real due to eq 344. We further split equation M487 along the lines of eq 345:
equation M488
(349)
equation M489
(350)
equation M490 is even, and equation M491 is odd. These two quantities are related by the fluctuation-dissipation theorem
equation M492
(351)
For the harmonic bath model (eq 95), these correlation functions can be calculated analytically and we get the following:
equation M493
(352)
This function is independent of temperature. We refer to it as the spectral density of bath coordinates and denote it equation M494 to distinguish from spectral densities of the Hamiltonian parameters used in eq 342.
Comparing eqs 342 and 352, we get
equation M495
(353)
Equation 340 may be simplified further by assuming a short-time memory kernel (the Markovian approximation) and replacing B(1)(tt′) in eq 340 by G(0)†(t′)B(1)(t). Note that even though the exact time propagation of B(1)(t) is described by the full Green’s function, G(t′), we can still use G(0)(t′) since we have already neglected higher than second order corrections in the system-bath coupling. This gives the following:
equation M496
(354)
where the one-exciton relaxation rate matrix is given by
equation M497
(355)
Dephasing of the second-order Y(2) variables can be calculated similarly from eq 332, since eq 332 is analogous to eq 331 for B. Following the procedure described above, we obtain the following:
equation M498
(356)
where k(Y) is the tetradic relaxation memory kernel of the Y variables:
equation M499
(357)
and An external file that holds a picture, illustration, etc.
Object name is nihms212366ig36.jpg(t) describes the free time evolution of Y in the absence of the bath.
In the Markovian approximation, we obtain
equation M500
(358)
where the two-exciton relaxation matrix is as follows:
equation M501
(359)
Exciton dephasing affects the B and Y variables whose Green’s functions involve coherences and oscillate at high, optical, frequencies. The Green’s function of the N variables is different: it describes the correlated dynamics of the conjugate variables [B with circumflex] and [B with circumflex] and contains low-frequency intraband oscillations related to intermolecular couplings.
Following the same procedure used in the derivation of eqs 340 and 356, we obtain
equation M502
(360)
The relaxation kernel, k(N), calculated to second order in the coupling with the bath is given by
equation M503
(361)
where equation M504 is the coherent N-variable Green’s function calculated in the absence of the bath. Invoking the Markovian approximation, we define the N-variable relaxation rates equation M505 where
equation M506
(362)
Kmn,kl [equivalent] Knm,lk is the Redfield relaxation superoperator.
The relaxation operator for Z variable is approximately constructed from γ and γ(Y) matrices by assuming that Z dephasing is equivalent to that of the product B*Y. The resulting NEE eqs 99102 contain no bath variables and include only the dephasing and relaxation matrices. The continuous spectral density bath models are widely used as described in Appendix G.
23. Appendix G: Exciton Dephasing and the Redfield Relaxation Kernel
In this section, we derive expressions for dephasing and relaxation rates of single-exciton eigenstates using the results of Appendix F when each chromophore is coupled to its own statistically independent bath with a continuous spectral density. Fluctuations of different chromophores are uncorrelated, and they all have the same spectral density so that
equation M507
(363)
where equation M508 is now independent of the system indices. The spectral density has the symmetry C″(ω)=C″(−ω) . Specific models for the spectral density, e.g., Ohmic, white noise, Lorentzian, etc., are often used for calculating the nonlinear responses of molecular systems.40,97,224 We assume that the bath has several collective modes α and recast the spectral density in the following form:
equation M509
(364)
where Cα(ω) is a dimensionless spectral density of mode α and λα is the coupling strength between. We define the auxiliary functions as follows:
equation M510
(365)
By performing a one-sided Fourier transform,
equation M511
(366)
we get
equation M512
(367)
equation M513
(368)
where PP stands for the principal part, i.e., integration where the singular point is removed (we have assumed that (1/ω)Cα(ω) is not singular at ω → 0); and we used πδ(x) = limη→0 η(x2 + η2)−1. We note the symmetry equation M514. The real part of [M with macron] is responsible for dephasing and transport rates, and the imaginary part represents spectral shifts.
Using this function, we can calculate the dephasing rate given by eq 355. We first transform this expression into the single-exciton basis, and then we neglect off-diagonal dephasing matrix elements (exciton coherence transfer). We get the dephasing rate for exciton e in terms of the [M with macron] function:
equation M515
(369)
where ωe′e = (εe′εe) and Φe4e3e2e1 = Σn ψne4 ψne3ψne2ψne1 (we assume that the exciton wave functions ψ are real). For the Y variables, we do not use their Green’s function explicitly. Instead, we use the exciton-scattering matrix, which depends only on single-exciton dephasings.
The Redfield relaxation rate, K (eq 362), for the single-exciton eigenstates is
equation M516
(370)
The expression for the population-transport rates is a special case (Ke′e′,ee elements)
equation M517
(371)
while the dephasing of single-exciton coherences equation M518. The pure dephasing (used in SOS expressions) can be simply obtained by equation M519. These compact expressions can be easily computed for an arbitrary spectral density.
We next adopt the overdamped Brownian oscillator spectral density,
equation M520
(372)
where λ corresponds to the system coupling strength with the collective bath mode and Λ is its timescale. Using eqs 364, 367, and 368, we get144
equation M521
(373)
with the auxiliary functions
equation M522
(374)
equation M523
(375)
νn = 2πnkBT/[variant Planck's over 2pi] are the Matsubara frequencies. The relaxation rates now are given in units of λ. This completes our description of the standard Redfield equations using second-order perturbation theory.
Assuming Gaussian statistics, an exact expression for diagonal bath fluctuations can be obtained for calculating the population (Pauli) relaxation block of the Redfield equations using the second order cumulant expansion. This expansion is correct to infinite order in the fluctuations. System-bath coupling of off-diagonal Hamiltonian elements is included perturbatively to second order. Exponentiation of the second-order correlation functions of diagonal bath fluctuations, which is exact for the harmonic bath, is then performed. Diagonal bath fluctuations, which only modulate the transport rate through the energy gap, can thus be included exactly, while off-diagonal fluctuations are treated by second-order perturbation theory. This gives for the population relaxation rate185
equation M524
(376)
equation M525
(378)
where equation M526 is the line-shape function. Here, a dot denotes the time derivative and
equation M527
(379)
equation M528
(380)
where λe4e3e2e1 = −limτ→∞[Im ġ e4e3e2e1 (τ)]. The modified Red-field theory eq 376 generalizes the population relaxation rate eq 371 to include an arbitrary fluctuation time scale. However, when the fluctuations are too slow, the transport can no longer be defined by simple Markovian rate equations. The calculated dephasing and relaxation rates can be used in calculations of Green’s functions in the signal expressions eqs 127129.
24. Appendix H: Exciton Scattering in Nonbosonic Systems
In section 6, we derived equations of motion for exciton variables using boson commutation relations. The bosonization procedure is described in Appendix C. In some applications, it may be preferable to use variables with the nonboson commutation relations. These should be incorporated in the equations of motion.
A generic nonbosonic system in the space of single and double excitons is defined by the commutation relation:
equation M529
(381)
We assume a linear transition dipole so that μ(2) = 0. The Heisenberg equation of motion (eq 63) truncated at the third order gives the following NEE:
equation M530
(382)
equation M531
(383)
equation M532
(384)
equation M533
(385)
where now Vmn,kl = 2Umn,kl – 2Σq Pmn,qkhql – 2Σqp Pmn,pqUpqkl.
The response functions and the signals are obtained in the form of eqs 127129. The scattering matrix for the coherent dynamics limit assumes the following form:
equation M534
(386)
In the MFA, this reduces to
equation M535
(387)
The corresponding scattering matrices of bosons are now obtained by simply taking P = 0, which gives eqs 325 and 327.
These commutation relations and the corresponding scattering matrix describe hard-core bosons, Pmn,kl = δmnδnkδkl, as well as Wannier excitons (electron–hole pairs) in semiconductors, where the commutation relations are more complex.108,114,354
25. Appendix I: Exciton-Scattering in Infinite Periodic Structures
Molecular crystals are periodic three-dimensional structures. Some molecular aggregates can be represented as infinite one-dimensional systems. These include J-aggregates and cylindrical dye aggregates.5,99,355,356 The quasiparticle scattering approach can be applied to the infinite systems by transforming from real space to the momentum (q) space using the Bloch states.
Consider a periodic D-dimensional system (D = 1, 2, 3) of linear extension L (the number of unit cells) made of cells with An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg; sites per unit cell. Thus, the number of chromophores is N= An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg × LD. We use periodic boundary conditions to represent the translational invariance. We label the sites by indices α, and each cell is identified by its position vector R. For a translationary invariant system, the intermode coupling Jmm [equivalent] JRα,R′α ′ = Jα,α ′ (R′– R) is a function of the distance between the cells R′– R. Our starting Hamiltonian is given by eq 57; however, for an infinite system, each oscillator m in real space is represented by a pair of indices Rα.
The one-exciton states of this system are the Bloch states. Each eigenstate ν has a pair of quantum numbers qλ, where λ denotes different Davydov’s sub-bands in the one–exciton band and the momentum q can have values 2π/(La)(–L/2, –L/2 + 1, ••• , L/2 – 1) (including 0) in each dimension where the lattice constant is a. The one-exciton wave functions are
equation M536
(388)
where An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpg = (La)D is the volume of the system and equation M537 are the eignstates of the cell,
equation M538
(389)
and
equation M539
(390)
Here, Jα,α′ (q) denotes the An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg × An external file that holds a picture, illustration, etc.
Object name is nihms212366ig9.jpg matrix with indices α and α′ where each matrix element depends parametrically on the vector q.
The one-exciton Green’s function in the one-exciton basis now reads
equation M540
(391)
where
equation M541
(392)
equation M542
(393)
We next assume a simplified form for the U matrix, Umn,mn = (Δm,n/4)(δmm δnn′+ δmn′ δnm′) (see eq 57), which reflects the boson scattering, i.e., the energy of the system in the real space when two molecules m and n are excited is εm + εn + Δmn. Similarly as for quadratic coupling, we use two indices for the infinite system Δmm′ = ΔRα,Rα = Δα,α (RR) and assume that it is nonzero only for short distances smaller than some cutoff parameter |RR| < lc.
The two-exciton Green’s function in real space is
equation M543
(394)
To calculate the Green’s function when Rn = Rm + r1 and equation M544 when r1 and r2 are within interaction radius lc, we define the reduced Green’s function equation M545 and expand it in the one-exciton basis,
equation M546
(395)
where
equation M547
(396)
is the unit cell’s Green’s function and
equation M548
(397)
Taking into account the translational invariance with respect to Rm and equation M549, we can transform the two-exciton Green’s function as
equation M550
(398)
By introducing the An external file that holds a picture, illustration, etc.
Object name is nihms212366ig55.jpgR × An external file that holds a picture, illustration, etc.
Object name is nihms212366ig55.jpgR matrix D,
equation M551
(399)
the exciton-scattering matrix may be recast as
equation M552
(400)
Since the optical wavelength is typically much longer than the system size, only Bloch states with zero lattice vector are optically active. The scattering matrix of these states then becomes with and
equation M553
(401)
with
equation M554
(402)
and
equation M555
(403)
Transforming this into zero-momentum Bloch state, we get
equation M556
(404)
This scattering matrix can be substituted into eqs 127129, where the summations over eigenstates are replaced by summations over zero-momentum Bloch states. For infinite systems, the index q is a continuous variable and all sums (1/An external file that holds a picture, illustration, etc.
Object name is nihms212366ig3.jpgq turn to be integrals (1/(2π)D) ∫ dq.
26. Appendix J: Doorway and Window Wavepackets
The doorway and window functions describe ground-state absorption and induced absorption/emission properties in the third-order response. These functions are defined by eq 170 in section 11.2. Expanded in the system eigenstates (see eq 176), we have185
equation M557
(405)
equation M558
(406)
with
equation M559
(407)
Here, equation M560 is the creation operator of either a single-exciton eigenstate a = e or a double-exciton eigenstate a = f.
equation M561
(408)
equation M562
(409)
The bleaching contribution is given by
equation M563
(410)
equation M564
(411)
The doorway and window functions may use the cumulant expansion. This gives
equation M565
(412)
equation M566
(413)
where the bath reorganization energy is λνν′, defined by eq 178.
27. Appendix K: Population Transfer with Slow Bath Fluctuations
In this appendix, we present the correlated contribution of the population-transfer pathways to the response functions described in section 11.3. The density matrix dynamics represented by the population-transfer diagrams now involves the characteristic single diagram shown in Figure 35. This diagram describes correlated propagation during the three time intervals and is given by the correlation function
Figure 35
Figure 35
Generalized diagram with population relaxation showing three density matrix propagation intervals. By a proper choice of the b and c states, this diagram can represent the density matrix evolution of all diagrams in Figure 7.
equation M567
(414)
Here, we assumed that the bath dynamics for different exciton states is correlated; a simpler case is given in ref 24. An external file that holds a picture, illustration, etc.
Object name is nihms212366ig37.jpg is evaluated by the cumulant expansion, and we get
equation M568
(415)
where
equation M569
(416)
Here, Im denotes the imaginary part. This expression contains line-shape functions that represent correlations during the time intervals t1 and t3.
The contribution of population-conserving terms to the response function now coincides with the coherent description. The contribution of the population-transfer diagram can be calculated using the cumulant expansion of
equation M570
(417)
where p denotes the summation over ia, iiia, iva, and via pathways: for ia, we have (–1)p [equivalent] 1, νp [equivalent] g, c [equivalent] e′, b [equivalent] g, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig38.jpg [equivalent] An external file that holds a picture, illustration, etc.
Object name is nihms212366ig39.jpg, left angle bracket... right angle bracket p left angle bracket ... right angle bracket *; for iiia, we have (–1)p [equivalent] –1, νp [equivalent] f, c [equivalent] f, b = e′,An external file that holds a picture, illustration, etc.
Object name is nihms212366ig38.jpg [equivalent] An external file that holds a picture, illustration, etc.
Object name is nihms212366ig40.jpg, left angle bracket... right angle bracketp [equivalent] left angle bracket... right angle bracket*; for iva, we have (–1)p [equivalent] 1, νp [equivalent] g, c [equivalent] g, b [equivalent] e′, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig38.jpg [equivalent] An external file that holds a picture, illustration, etc.
Object name is nihms212366ig39.jpg, left angle bracket ... right angle bracketp [equivalent] left angle bracket ... right angle bracket; for via, we have (–1)p [equivalent] –1, νp [equivalent] f, c [equivalent] e′, b [equivalent] f, An external file that holds a picture, illustration, etc.
Object name is nihms212366ig38.jpg[equivalent] An external file that holds a picture, illustration, etc.
Object name is nihms212366ig41.jpg, left angle bracket... right angle bracketp [equivalent] left angle bracket... right angle bracket. We next use μab = μba and the correlation function of eq 414. These yield in the following expression
equation M571
(418)
The dynamics of An external file that holds a picture, illustration, etc.
Object name is nihms212366ig37.jpg solely comes from the doorway and window functions. The phase functions in eq 179 are obtained from eq 418 by sorting out all Liouville space pathways.
28. Appendix L: Orientationally Averaged Molecular Properties
Here, we present expressions for third-order signals in isotropic ensembles of molecules. As shown in eq 251, the response function has several linearly independent configurations that differ by the polarization ν and wavevector k orientations in the laboratory coordinate frame. The complete expressions for this pulse-polarization configuration for the linear response are
equation M572
(419)
For third-order techniques, we set kI [equivalent] (u, v, w) = (–1, 1, 1), for kII [equivalent] (u, v, w) = (1, –1, 1), and for kIII [equivalent] (u, v, w) = (1, 1, –1). We then get
equation M573
(420)
where εγαβ is the antisymmetric Levi–Civita tensor. These can be obtained from orientational averaging of tensor products of arbitrary sets of vectors di and second-rank tensors ti
equation M574
(421)
equation M575
(422)
equation M576
(423)
equation M577
(424)
equation M578
(425)
equation M579
(426)
equation M580
(427)
equation M581
(428)
equation M582 are column vectors that depend on the laboratory-frame unit vectors ν:
equation M583
(429)
equation M584
(430)
The matrices M are system-independent:
equation M585
(431)
equation M586
(432)
Finally, V are column vectors that depend on the system properties:
equation M587
(433)
equation M588
(434)
equation M589
(435)
equation M590
(436)
equation M591
(437)
where d[middle dot in circle]t [equivalent] Σα3α 2α 1εα 3α 2α 1dα 3tα 2α 1, (td)α [equivalent] Σα′tαα′dα′, (dt)α [equivalent] Σα′ dα′tα ′α, and t [equivalent] Σαtαα; consequently, equation M592. By comparing with eqs 419 and 420 and replacing the vectors d with μ and m and the tensors t with quadrupoles q, we can calculate the orientationally averaged currents.
For the linear response in the exciton basis, we have equation M593 and the signal must propagate in the same direction of the incoming field. We can thus set κ = z. Orientational averaging for exciton e gives
equation M594
(438)
The orientational factors given by eqs 424437 relate the signals to molecular geometry. The CI terms can be identified by examining the products of transition dipoles. The electric transition dipole μξ and coordinate rm are parity-odd, while the magnetic transition dipole mξ and tensors qξ are even. By performing the parity operation on the response functions and inspecting the parity symmetry, we can connect the susceptibility with the chirality. All zero-order in wavevector contributions to odd (linear and cubic) response must be NC, since they involve either two or four electric transition dipoles. The corresponding first-order contributions in wavevectors are CI.
1. Steed JW, Atwood JL. Supramolecular Chemistry. Wiley; New York: 2000.
2. Kobayashi T. J-aggregates. World Scientific; Singapore: 1996.
3. Mobius D, Miller R, editors. Organized Monolayers and Assemblies: Structure, Processes and Function. Elsevier; New York: 2002.
4. van Amerongen H, Valkunas L, van Grondelle R. Photosynthetic Excitons. World Scientific; Singapore: 2000.
5. Agranovich VM, La Rocca GC. International School of Physics Enrico Fermi. Vol. 149 IOS Press; Amsterdam: 2002. Organic Nanostructures: Science and Applications.
6. Möbius D, Kuhn H. J Appl Phys. 1988;64:5138.
7. Zhang Q, Atay T, Tischler JR, Bradley MS, Bulovic V, Nurmikko AV. Nature Nanotechnol. 2007;2:555. [PubMed]
8. Forrest SR, Thompson ME. Chem Rev. 2007;107:923.
9. McDermott G, Prince SM, Freer AA, Hawthornthwaite-Lawless AM, Papiz MZ, Cogdell RJ, Isaacs NJ. Nature. 1995;374:517.
10. Kristein S, von Berlepsch H, Böttcher C, Burger C, Ouart A, Reck G, Dähne S. ChemPhysChem. 2000;3:146. [PubMed]
11. Didraga C. Dissertation. 2004. Optical Properties of Frenkel Excitons in Cylindrical Molecular Aggregates.
12. Pugzlys A, Hania PR, Augulis R, Duppen K, van Loosdrecht PHM. Int J Photoenergy. 2006 ID 29623.
13. Tomalia DA, Naylor AM, Goddard WA., III Angew Chem, Int Ed. 1990;29:138.
14. Love JC, Estroff LA, Kriebel JK, Nuzzo RG, Whitesides GM. Chem Rev. 2005;105(4):1103. [PubMed]
15. Olson J. Photosynth Res. 2004;80:181. [PubMed]
16. Jordan P, Fromme P, Witt H, Klukas O, Saenger W, Krauss N. Nature. 2001;411:909. [PubMed]
17. Davydov AA. Theory of Molecular Excitons. Mc.Graw-Hill; New York: 1962.
18. Rashba EI, Sturge MD, editors. Excitons. Elsevier; New York: 1987.
19. Mukamel S, Chemla DS, editors. Chem Phys. Vol. 210. Elsevier; New York: 1996. pp. 1–388.
20. Markovitsi D, Small G, editors. Chem Phys. Vol. 275. Elsevier; New York: 2002. pp. 1–398.
21. Pisliakov AV, Mančal T, Fleming GR. J Chem Phys. 2006;124:234505. [PubMed]
22. Mančal T, Valkunas L, Fleming GR. Chem Phys Lett. 2006;432:301.
23. Kjellberg P, Brüggemann B, Pullerits T. Phys Rev B: Condens Matter. 2006;74:024303.
24. Abramavicius D, Valkunas L, Mukamel S. Europhys Lett. 2007;80:17005.
25. Jelley EE. Nature. 1937;139:631.
26. Scheibe G. Angew Chem. 1937;50:212.
27. Ohta M, amd Yang K, Fleming GR. J Chem Phys. 2001;115:7609.
28. Fidder H. Chem Phys. 2007;341:158.
29. Heijs DJ, Dijkstra AG, Knoester J. Chem Phys. 2007;341:230.
30. Knapp EW. Chem Phys. 1984;85:73.
31. Gross M, Haroche S. Phys Rep. 1982;93:301.
32. Weitekamp DP, Duppen K, Wiersma DA. Phys Rev A. 1983;27:3089.
33. Spano FC, Kuklinski JR, Mukamel S. Phys Rev Lett. 1990;65:211. [PubMed]
34. Meier T, Zhao Y, Chernyak V, Mukamel S. J Chem Phys. 1997;107:3876.
35. Zhao Y, Meier T, Zhang W, Chernyak V, Mukamel S. J Phys Chem. 1999;103:3954.
36. Chernyak V, Meyer T, Tsiper E, Mukamel S. J Phys Chem A. 1999;103:10294.
37. Gobets B, van Grondelle R. Biochim Biophys Acta. 2001;1507:80. [PubMed]
38. Muller M, Niklas J, Lubitz W, Holzwarth A. Biophys J. 2003;85:3899. [PubMed]
39. Novoderezhkin VI, Rutkauskas D, van Grondelle R. Biophys J. 2006;90:2890. [PubMed]
40. Brixner T, Stenger J, Vaswani HM, Cho M, Blankenship RE, Fleming GR. Nature. 2005;434:625. [PubMed]
41. Engel GS, Calhoun TR, Read EL, Ahn TK, Mančal T, Cheng YC, Blankenship RE, Fleming GR. Nature. 2007;446:782. [PubMed]
42. Mukamel S, Loring RF. J Opt Soc Am B. 1986;3:595.
43. Fried LE, Mukamel S. Adv Chem Phys. 1993;84:435.
44. de Boeij WP, Pshenichnikov MS, Wiersma DA. Annu Rev Phys Chem. 1998;49:99. [PubMed]
45. Cho M, Scherer NF, Fleming GR, Mukamel S. J Chem Phys. 1992;96:5618.
46. Kovalevskij V, Gulbinas V, Piskarskas A, Hines MA, Scholes GD. Phys Status Solidi B. 2004;241(8):1986.
47. Doust AB, van Stokkum IHM, Larsen DS, Wilk KE, Curmi PMG, van Grondelle R, Scholes GD. J Phys Chem B. 2005;109:14219. [PubMed]
48. Cho M, Fleming GR. J Chem Phys. 2005;123:114506. [PubMed]
49. Berera R, Herrero C, van Stokkum IHM, Vengris M, Kodis G, Palacios RE, van Amerogen H, van Grondelle R, Gust D, Moore TA, Moore AL, Kennis JTM. Proc Natl Acad Sci USA. 2006;103:5343. [PubMed]
50. van Stokkum IHM, Larsen DS, van Grondelle R. Biochim Biophys Acta. 2004;82:1657.
51. Ruban AV, Berera R, Ilioaia C, van Stokkum IHM, Kennis JTM, Pascal AA, van Amerongen H, Robert B, Horton P, van Grondelle R. Nature. 2007;450:575. [PubMed]
52. Mukamel S. Principles of Nonlinear Optical Spectroscopy. Oxford University Press; New York: 1995.
53. Ernst RR, Bodenhausen G, Wokaun A. Principles of Nuclear Magnetic Resonance in One and Two Dimensions. Clarendon; Oxford, U.K: 1987.
54. Tanimura Y, Mukamel S. J Chem Phys. 1993;99:9496.
55. Mukamel S, Piryatinski A, Chernyak V. Acc Chem Res. 1999;32:145.
56. Milne CJ, Li Y, Miller RJD. Two Dimensional 5th Order Raman Spectroscopy. In: Torre R, editor. Time-Resolved Spectroscopy in Complex Liquids. Chapter 2. Springer-Verlag; New York: 2007. pp. 205–259.
57. Scheurer C, Mukamel S. J Chem Phys. 2001;115:4989.
58. Scheurer C, Mukamel S. J Chem Phys. 2002;116:6803.
59. Scheurer C, Mukamel S. Bull Chem Soc Jpn. 2002;75:989.
60. Chung H, Khalil M, Smith AW, Ganim Z, Tokmakoff A. Proc Natl Acad Sci USA. 2005;102:612. [PubMed]
61. Kim YS, Liu L, Axelsen PH, Hochstrasser RM. Proc Natl Acad Sci USA. 2008;105:7720. [PubMed]
62. Kwak KW, Park S, Fayer MD. Proc Natl Acad Sci USA. 2007;104:14221. [PubMed]
63. Eaves JD, Loparo JJ, Fecko CJ, Roberts TJ, Tokmakoff A, Geissler PL. Proc Natl Acad Sci USA. 2005;102:13019. [PubMed]
64. Kraemer D, Cowan ML, Paarmann A, Huse N, Nibbering ETJ, Elsasser T, Miller RJD. Proc Natl Acad Sci USA. 2008;105:437. [PubMed]
65. Cho M. Chem Rev. 2008;108:1331. [PubMed]
66. Hamm P, Helbing J, Bredenbeck J. Annu Rev Phys Chem. 2008;59:291. [PubMed]
67. Zhuang W, Hayashi T, Mukamel S. Angew Chem. 2009 in press.
68. Chernyak V, Zhang WM, Mukamel S. J Chem Phys. 1998;109:9587.
69. Zhang W, Chernyak V, Mukamel S. J Chem Phys. 1999;110:5011.
70. Mukamel S. Annu Rev Phys Chem. 2000;51:691. [PubMed]
71. Jonas DM. Annu Rev Phys Chem. 2003;54:425. [PubMed]
72. Nemeth A, Milota F, Sperling J, Abramavicius D, Mukamel S, Kauffman HF. Chem Phys Lett. 2009;469:130.
73. Zhang TH, Kuznetsova I, Meier T, Li XQ, Mirin RP, Thomas P, Cundiff ST. Proc Natl Acad Sci USA. 2007;104:14227. [PubMed]
74. Li X, Zhang T, Borca CN, Cundiff ST. Phys Rev Lett. 2006;96:057406. [PubMed]
75. Yang L, Mukamel S. Phys Rev B. 2008;77:075335.
76. Tian P, Keusters D, Suzaki Y, Warren WS. Science. 2003;300:1553. [PubMed]
77. Brixner T, Stiopkin IV, Fleming GR. Opt Lett. 2004;29:884. [PubMed]
78. Chernyak V, Mukamel S. Phys Rev Lett. 1995;74:4895. [PubMed]
79. Abramavicius D, Mukamel S. J Chem Phys. 2006;124:034113. [PubMed]
80. Cho M, Vaswani HM, Brixner T, Stenger J, Fleming GR. J Phys Chem B. 2005;109:10542. [PubMed]
81. Abramavicius D, Voronine DV, Mukamel S. Biophys J. 2008;94:3613. [PubMed]
82. Li Z, Abramavicius D, Zhuang W, Mukamel S. Chem Phys. 2007;341:29. [PMC free article] [PubMed]
83. Li Z, Abramavicius D, Mukamel S. J Am Chem Soc. 2008;130:3509. [PMC free article] [PubMed]
84. Oskouei AA, Bräm O, Cannizzo A, van Mourik F, Tortschanoff A, Chergui M. Chem Phys. 2008;350:104.
85. Bressler C, Chergui M. Chem Rev. 2004;104:1781. [PubMed]
86. Tanaka S, Chernyak V, Mukamel S. Phys Rev A. 2001;63:063405.
87. Schweigert IV, Mukamel S. Phys Rev Lett. 2007;99:163001. [PubMed]
88. Cho M, Fleming GR, Mukamel S. J Chem Phys. 1993;98:5314.
89. Cho M. J Chem Phys. 2002;116:1562.
90. Cho M. J Chem Phys. 2003;119:7003.
91. Abramavicius D, Mukamel S. Chem Phys. 2005;318:50.
92. Abramavicius D, Mukamel S. J Chem Phys. 2005;122:134305. [PubMed]
93. Berova N, Nakanishi K, Woody RW, editors. Circular Dichroism Principles and Applications. 2. John Wiley & Sons, Inc; New York: 2000.
94. Zhuang W, Abramavicius D, Mukamel S. Proc Natl Acad Sci USA. 2006;103:18934. [PubMed]
95. Zhuang W, Abramavicius D, Mukamel S. Proc Natl Acad Sci USA. 2005;102:7443. [PubMed]
96. Abramavicius D, Mukamel S. J Chem Phys. 2004;120:8373. [PMC free article] [PubMed]
97. Voronine DV, Abramavicius D, Mukamel S. J Chem Phys. 2007;126:044508. [PubMed]
98. Wright J, editor. Molecular Crystals. 2. Cambridge University Press; New York: 1995.
99. Knoester J, Agranovich VM. Thin Films Nanostruct. 2003:31.
100. Mukamel S, Abramavicius D. Chem Rev. 2004;104:2073. [PubMed]
101. Spano FC, Mukamel S. Phys Rev A. 1989;40:5783. [PubMed]
102. Spano FC, Mukamel S. Phys Rev Lett. 1991;66:1197. [PubMed]
103. Knoester J, Mukamel S. Phys Rep. 1991;205:1.
104. Leegwater JA, Mukamel S. Phys Rev A. 1992;46:452. [PubMed]
105. Mukamel S. Many-Body Effects in Nonlinear Susceptibilities: Beyond the Local-Field Approximation. In: Zyss J, Kelley P, Liao P, editors. Molecular Nonlinear Optics. Elsevier; New York: 1993.
106. Chernyak V, Wang N, Mukamel S. Phys Rep. 1995;263:213.
107. Mukamel S, Oszwaldowski R, Abramavicius D. Phys Rev B. 2007;75:245305.
108. Breuer HP, Petruccione F. The theory of open quantum systems. Oxford University Press; Oxford, U.K: 2002.
109. Tanimura Y. J Phys Soc Jpn. 2006;75:082001.
110. Weiss U. Quantum Dissipative Systems. World Scientific; River Edge, NJ: 2006.
111. van Kampen NG. Stochastic Processes in Physics and Chemistry. 3. North-Holland; The Netherlands: 2007.
112. Mukamel S. Phys Rev A. 1983;28:3480.
113. Jansen TlC, Hayashi T, Zhuang W, Mukamel S. J Chem Phys. 2005;123:114504. [PubMed]
114. Axt VM, Mukamel S. Rev Mod Phys. 1998;70:145.
115. Haug H, Koch SW. Quantum Theory of the Optical and Electronic Properties of Semiconductors. 4. World Scientific Publishing Company; River Edge, NJ: 2004.
116. Mukamel S, Oszwaldowski R, Yang L. J Chem Phys. 2007;127:221105. [PMC free article] [PubMed]
117. Loudon R, editor. Quantum theory of light. 3. Oxford University Press; New York: 2000.
118. Kubo R, Toda M, Hashitsume N. Statistical Physics II: Nonequilibrium Statistical Mechanics. 2. Springer; New York: 1991.
119. Zwanzig R. Nonequilibrium Statistical Mechanics. Oxford University Press; New York: 2001.
120. Breuer H-P, Petruccione F. The theory of open quantum systems. Oxford University Press; New York: 2002.
121. Lindblad G. Commun Math Phys. 1976;48:119.Gorini V, Kossakowski A, Sudarshan ECG. J Math Phys. 1976;17:821.
122. Stenholm S. Quant Semiclass Opt. 1996;8:297.
123. Hahn EL. Phys Rev. 1950;80:580.
124. Aartsma TJ, Wiersma DA. Phys Rev Lett. 1976;36:1360.
125. Samartsev VV. J Appl Spectrosc. 1979;30:401.
126. Prokhorenko VI, Holzwarth AR, Nowak FR, Aartsma TJ. J Phys Chem B. 2002;106(38):9923.
127. Schweigert IV, Mukamel S. Phys Rev A. 2008;77:033802.
128. Spano FC, Agranovich V, Mukamel S. J Chem Phys. 1991;95:1400.
129. Šanda F, Mukamel S. J Chem Phys. 2007;127:154107. [PubMed]
130. Meier T, Chernyak V, Mukamel S. J Phys Chem B. 1997;101:7332.
131. Kuhn O, Chernyak V, Mukamel S. J Chem Phys. 1996;105:8586.
132. Tretiak S, Middleton C, Chernyak V, Mukamel S. J Phys Chem B. 2000;104:4519.
133. Tretiak S, Middleton C, Chernyak V, Mukamel S. J Phys Chem B. 2000;104:9540.
134. Knox RS, Small G, Mukamel S. Chem Phys. 2002;281:1.
135. Mukamel S, Berman O. J Chem Phys. 2003;119:12194.
136. Agranovich VM, Toshich BS. Zh Eksp Teor Fiz. 1967;53:149.Sov Phys JETP. 1968;26:104.
137. Ilinskaia AV, Ilinski KN. J Phys A: Math Gen. 1996;29:L23.
138. Kopietz P. Bosonization of Interacting Fermions in Arbitrary Dimensions. Springer-Verlag; New York: 1997. [PubMed]
139. Gogolin AO, Nersesyan AA, Tsvelik AM. Bosonization and Strongly Correlated Systems. Cambridge University Press; New York: 2004.
140. Cheng Y, Silbey R. J Phys Chem B. 2005;109:21399. [PubMed]
141. Pollard W, Friesner R. J Chem Phys. 1994;100:5054.
142. Kohen D, Tannor D. J Chem Phys. 1997;107:5141.
143. Yan YJ, Xu RX. Annu Rev Phys Chem. 2005;126:114102.
144. Chernyak V, Minami T, Mukamel S. J Phys Chem. 2000;112:7953.
145. Kirkwood JC, Scheurer C, Chernyak V, Mukamel S. J Chem Phys. 2000;114:2419.
146. Davies E. Commun Math Phys. 1974;39:91.
147. Daffer S, Wódkiewicz K, Mciver J. J Mod Opt. 2004;15:1843.
148. Suárez A, Silbey R, Oppenheim I. J Chem Phys. 1992;97:5101.
149. Gaspard P, Nagaoka M. J Chem Phys. 1999;111:5668.
150. Renger T, May V, Kühn O. Phys Rep. 2001;343:137.
151. Michalet X, Kapanidis AN, Laurence T, Pinaud F, Doose S, Pflughoefft M, Weiss S. Annu Rev Biophys Biomol Struct. 2003;32:161. [PubMed]
152. Kong X, Nir E, Hamadani K, Weiss S. J Am Chem Soc. 2007;129:4643. [PMC free article] [PubMed]
153. Dexter DL. J Chem Phys. 1953;21:836.
154. Damjanovic A, Thorsten R, Schulten K. Phys Rev E. 1999;59:3293.
155. Zhuang W, Abramavicius D, Hayashi T, Mukamel S. J Phys Chem B. 2006;110:3362. [PMC free article] [PubMed]
156. Abramavicius D, Zhuang W, Mukamel S. J Phys B: At Mol Opt Phys. 2006;36:5051.
157. Frank HA, Young AJ, Britton G, Cogdell RJ. The Photochemistry of Carotenoids. Kluwer Academic Publishers; Norwell, MA: 1999.
158. Law CJ, Cogdell RJ. The Primary Processes of Photosynthesis. Royal Society Of Chemistry; London: 2008. The light-harvesting system of purple anoxygenic photosynthetic bacteria; pp. 205–259.
159. van Grondelle R, Dekker J, Gillbro T, Sundstrom V. Biochim Biophys Acta. 1994;1187:1.
160. Sugisaki M, Fujii R, Cogdell R, Hashimoto H. Photosynth Res. 2008;95:309. [PubMed]
161. Sugisaki M, Fujiwara M, Yanagi K, Cogdell R, Hashimoto H. Photosynth Res. 2008;95:299. [PubMed]
162. Holt N, Kennis J, Dall’Osto L, Bassi R, Fleming G. Chem Phys Lett. 2003;379:305.
163. Papagiannakis E, Kumar Das S, Gall A, van Stokhum IHM, Robert B, van Grondelle R, Frank H, Kennis JTM. J Phys Chem B. 2003;107:5642.
164. Holt N, Kennis J, Fleming G. J Phys Chem B. 2004;108:19029.
165. Dreuw A, Fleming G, Head-Gordon M. Biochem Soc Trans. 2005;33:858. [PubMed]
166. Holt NE, Zigmantas D, Valkunas L, Li XP, Niyogi KK, Fleming GR. Science. 2005;307:433. [PubMed]
167. Cogdell RJ, Mullineaux C. Photosynth Res. 2008;95:117. [PubMed]
168. Fenna R, Matthews B. Nature. 1975;258:573.
169. Li YF, Zhou W, Blankenship R, Allen J. J Mol Biol. 1997;271:456. [PubMed]
170. Louwe RJW, Vrieze J, Hoff AJ, Aartsma TJ. J Phys Chem B. 1997;101(51):11280.
171. Vulto S, de Baat M, Louwe R, Permentier H, Neef T, Miller M, van Amerongen H, Aartsma T. J Phys Chem B. 1998;102:9577.
172. Lee H, Cheng YC, Fleming G. Science. 2007;316:1462. [PubMed]
173. Zigmantas D, Read EL, Mančal T, Brixner T, Gardiner AT, Cogdell RJ, Fleming GR. Proc Natl Acad Sci USA. 2006;103:12672. [PubMed]
174. Sundström V, Pullerits T, van Grondelle R. J Phys Chem B. 1999;103:2327.
175. Pullerits T, Sundström V. Acc Chem Res. 1996;29:381.
176. Koepke J, Hu X, Muenke C, Schulten K, Michel H. Structure. 1996;4:581. [PubMed]
177. Hu X, Ritz T, Damjanovic A, Autenrieth F, Schulten K. Q Rev Biophys. 2002;35:1. [PubMed]
178. Schröder M, Kleinekathöfer U, Schreiber M. J Chem Phys. 2006;124:084903. [PubMed]
179. Schröder M, Schreiber M, Kleinekathöfer U. J Lumin. 2007;125:126.
180. Jang S, Silbey R. J Chem Phys. 2003;118:9312.
181. Jang S, Silbey R. J Chem Phys. 2003;118:9324.
182. Dalhbom M, Pullerits T, Mukamel S, Sundström V. J Phys Chem B. 2001;105:5515.
183. Chernyak V, Mukamel S. J Chem Phys. 1996;105:4565.
184. Meier T, Chernyak V, Mukamel S. J Chem Phys. 1997;107:8759.
185. Zhang WM, Meier T, Chernyak V, Mukamel S. J Chem Phys. 1998;108:7763.
186. Pullerits T, Chachisvilis M, Sundström V. J Phys Chem. 1996;100:10787.
187. Trinkunas G, Herek J, Polívka T, Sundström V, Pullerits T. Phys Rev Lett. 2001;86:4167. [PubMed]
188. Yang M, Damjanović A, Vaswani H, Fleming G. Biophys J. 2003;85:140. [PubMed]
189. Liu Z, Yan H, Wang K, Kuang T, Zhang J, Gui L, An X, Chang W. Nature. 2004;428:287. [PubMed]
190. Biesiadka J, Loll B, Kern J, Irrgang KD, Zouni A. Phys Chem Chem Phys. 2004;6:4733.
191. Vaitekonis S, Trinkunas G, Valkunas L. Photosynth Res. 2005;86:185. [PubMed]
192. Novoderezhkin V, Salverda J, van Amerongen H, van Grondelle R. J Phys Chem B. 2003;107:1893.
193. Novoderezhkin V, MAP, van Amerongen H, van Grondelle R. J Phys Chem B. 2005;109:10493. [PubMed]
194. Raszweski G, Saenger W, Renger T. Biophys J. 2005;88:986. [PubMed]
195. Novoderezhkin V, Dekker J, van Amerongen H, van Grondelle R. Biophys J. 2007;93:1293. [PubMed]
196. Cho M, Brixner T, Stiopkin I, Vaswani H, Flemig G. J Chin Chem Soc. 2006;53:15.
197. Yang M, Fleming GR. Chem Phys. 2002;282:163.
198. Pomyalov A, Tannor D. J Chem Phys. 2005;123:204111. [PubMed]
199. Kühn O, Sundström V, Pullerits T. Chem Phys. 2002;275:15.
200. Polivka T, Pullerits T, Herek J, Sundström V. J Phys Chem B. 2000;104:1088.
201. Campillo AJ, Shapiro SL, Kollman VH, Winn KR, Hyer RC. Biophys J. 1976;16:93. [PubMed]
202. Brüggemann B, Herek JL, Sundstrom V, Pullerits T, May V. J Phys Chem B. 2001;105:11391.
203. Brüggemann B, Kjellberg T, Pullerits P. Chem Phys Lett. 2007;444:192.
204. Aeschlimann M, Bauer M, Bayer D, Brixner T, Garcia de Abajo FJ, Pfeiffer W, Rohmer M, Spindler C, Steeb F. Nature. 2007;446:301. [PubMed]
205. Mukamel S, Tortschanoff A. Chem Phys Lett. 2002;357:327.
206. Abramavicius D, Voronine DV, Mukamel S. Proc Natl Acad Sci USA. 2008;105:8525. [PubMed]
207. Thouless DJ. J Phys C: Solid State Phys. 1972:77.
208. Piryatinski A, Asher SA, Mukamel S. J Phys Chem A. 2002;106:3524.
209. Dennison D. Rev Mod Phys. 1940;12:175.
210. Abramavicius D, Palmieri B, Mukamel S. Chem Phys. 2009;357:79. [PMC free article] [PubMed]
211. Yan Y, Mukamel S. Phys Rev E. 1990;41:6485. [PubMed]
212. Zhang WM, Meier T, Chernyak V, Mukamel S. Philos Trans R Soc London, Ser A. 1998;356:405.
213. Kubo R. J Math Phys. 1963;4:174.
214. Gamliel D, Levanon H. Stochastic Processes in Magnetic Resonance. World Scientific; New York: 1995.
215. Kubo R. A stochastic theory of line-shape and relaxation. In: ter Haar D, editor. Fluctuation, relaxation and resonance in magnetic systems. Oliver & Boyd; Edingburgh, U.K: 1962. p. 23.
216. Risken H, Frank T. The Fokker-Planck Equation: Methods of Solutions and Applications. 2. Springer; New York: 1996.
217. Anderson PW. J Phys Soc Jpn. 1954;9:935.
218. Anderson PW. J Phys Soc Jpn. 1954;9:316.
219. Zhao Y, Chernyak V, Mukamel S. J Phys Chem A. 1998;102:6614.
220. Šanda F, Mukamel S. J Chem Phys. 2006;125:014507. [PubMed]
221. Zheng J, Kwak K, Asbury J, Chen X, Piletic IR, Fayer MD. Science. 2005;309:1338. [PubMed]
222. SŠanda F, Mukamel S. J Phys Chem B. 2008;112:14212. [PMC free article] [PubMed]
223. Jansen TLC, Zhuang W, Mukamel S. J Chem Phys. 2004;121:10577. [PubMed]
224. Caldeira AO, Leggett AJ. Physica A. 1983;121:587.
225. Zurek WH. Rev Mod Phys. 2003;75:715.
226. Kubo R, Toyozawa Y. Prog Theor Phys (Kyoto) 1955;13:162.
227. Balian R, Berezin E. Il Nuovo Cimento, B. 1969;64:37.
228. Osad’ko IS, editor. Springer Series in Chemical Physics. Vol. 69 Springer; Berlin: 2002.
229. Egorov SA, Rabani E, Berne BJ. J Phys Chem B. 1999;103:10978.
230. Zusman D. Chem Phys. 1980;49:295.
231. Garg A, Onuchic J, Ambegaokar V. J Chem Phys. 1985;83:4491.
232. Klafter J, Shlesinger M, Zumofen G. Phys Today. 1996;49:33.
233. SŠanda F, Mukamel S. Phys Rev E. 2006;73:011103. [PubMed]
234. SŠanda F, Mukamel S. Phys Rev Lett. 2007;98:080603. [PubMed]
235. Bel G, Barkai E. Phys Rev Lett. 2005;94:240602.
236. Lough WJ, Wainer IW, editors. Chirality in Natural and Applied Science. CRC Press; Boca Raton, FL: 2002.
237. Barron LD. Molecular Light Scattering and Optical Activity. 2. Cambridge University Press; New York: 2004.
238. Buckingham AD, Fischer P. Optical Response of a Chiral Liquid. In: Hicks JM, editor. Chirality Physical Chemistry. Oxford University Press; New York: 2002. pp. 119–129.
239. Nafie LA. Annu Rev Phys Chem. 1997;48:357. [PubMed]
240. Wagersreiter T, Mukamel S. J Chem Phys. 1996;105:7997.
241. Somsen OJG, van Grondelle R, van Amerogen H. Biophys J. 1996;71:1934. [PubMed]
242. Georgakopoulou S, van Grondelle R, van der Zwan G. J Phys Chem B. 2006;110:3344. [PubMed]
243. Chin DH, Woody RW, Rohl CA, Baldwin RL. Proc Natl Acad Sci USA. 2002;99:15416. [PubMed]
244. Sreerama N, Venyaminov SY, Woody RW. Anal Biochem. 2000;287:243. [PubMed]
245. Sreerama N, Woody RW. Anal Biochem. 2000;287:252. [PubMed]
246. Eker F, Cao X, Nafie LA, Huang Q, Schweitzer-Stenner R. J Phys Chem B. 2003;107:358.
247. Kubelka J, Gangani RA, Silva D, Keiderling TA. J Am Chem Soc. 2002;124:5325. [PubMed]
248. Hilario J, Kubelka J, Keiderling TA. J Am Chem Soc. 2003;125:7562. [PubMed]
249. Hecht L, Barron LD. Raman Optical Activity. John Wiley & Sons Ltd; New York: 1996.
250. Barron LD, Blanch EW, McColl IH, Syme CD, Hecht L, Nielsen K. Spectroscopy. 2003;17:101.
251. Andrews DL, Allcock P. Optical Harmonics in Molecular Systems. Wiley-VCH; New York: 2002.
252. Byers JD, Yee HI, Hicks JM. J Chem Phys. 1994;101:6233.
253. Hicks JM, Petralli-Mallow T, Byers JD. Faraday Discuss. 1994;99:341. [PubMed]
254. Simpson GJ. J Chem Phys. 2002;117:3398.
255. Burke BJ, Moad AJ, Polizzi MA, Simpson GJ. J Am Chem Soc. 2003;125:9111. [PubMed]
256. Moad AJ, Simpson GJ. J Phys Chem B. 2004;108:3548.
257. Fischer P, Buckingham AD, Albrecht AC. Phys Rev A. 2001;64:053816.
258. Fischer P, Beckwitt K, Wise FW, Albrecht AC. Chem Phys Lett. 2002;352:463.
259. Choi JH, Cho M. J Phys Chem A. 2007;111:5176. [PubMed]
260. Choi JH, Cho M. Chem Phys. 2007;341:57.
261. Lee K-K, Oh K-I, Lee H, Joo C, Han H, MC ChemPhy-sChem. 2007;8:2218. [PubMed]
262. Voronine DV, Abramavicius D, Mukamel S. Ultrafast Phenomena XVI 2009. Springer; New York: in press.
263. Voronine DV, Abramavicius D, Mukamel S. J Chem Phys. 2006;125:224504. [PubMed]
264. Keusters D, Tan H, Warren W. J Phys Chem A. 1999;103:10369.
265. Andrews DL, Thirunamachandran T. J Chem Phys. 1977;67:5026.
266. Yeremenko S, van Stokkum IHM, Moffat K, Hellingwerf KJ. Biophys J. 2006;90:4224. [PubMed]
267. Kobayashi T, Okada T, Kobayashi T, Nelson K, De Silvestri S, editors. Ultrafast Phenomena XIV. Springer; New York: 2005.
268. Miller RJD, Weiner AM, Corcum P, Jonas DM, editors. Ultrafast Phenomena XV. Springer; New York: 2007.
269. Nuernberger P, Vogt G, Brixner T, Gerber G. Phys Chem Chem Phys. 2007;9:1.
270. Dantus M, Lozovoy VV. Chem Rev. 2004;104:1813. [PubMed]
271. Rice S, Zhao M. Optical Control of Molecular Dynamics. Wiley; New York: 2000.
272. Shapiro M, Brumer P. Principles of the Quantum Control of Molecular Processes. John Wiley and Sons; Hoboken, NJ: 2003.
273. Warren W, Rabitz H, Dahleh M. Science. 1993;259:1581. [PubMed]
274. Rabitz H. Science. 2003;299:525. [PubMed]
275. Brixner T, Kiefer B, Gerber G. Chem Phys. 2001;267:241.
276. Brüggemann B, May V. J Phys Chem B. 2004;108(29):10529.
277. Brüggemann B, May V. Chem Phys Lett. 2004;400:573.
278. Brüggemann B, Pullerits T, May V. J Photochem Photobiol A. 2006;180:322.
279. Voronine D, Abramavicius D, Mukamel S. J Chem Phys. 2006;124:034104. [PubMed]
280. Pshenichnikov MS, de Boeij WP, Wiersma DA. Phys Rev Lett. 1996;76:4701. [PubMed]
281. Lozovoy VV, Dantus M. ChemPhysChem. 2005;6:1970. [PubMed]
282. Dantus M. Annu Rev Phys Chem. 2001;52:639. [PubMed]
283. Prokhorenko VI, Nagy AM, Waschuk SA, Brown LS, Birge RR, Miller RJD. Science. 2006;313:1257. [PubMed]
284. Vogt G, Krampert G, Niklaus P, Nuernberger P, Gerber G. Phys Rev Lett. 2005;94:068305. [PubMed]
285. Herek JL, Wohlleben W, Cogdell RJ, Zeidler D, Motzkus M. Nature. 2002;417:533. [PubMed]
286. Wohlleben W, Buckup T, Herek JL, Motzkus M. ChemPhy-sChem. 2005;6:850.
287. Vandersypen LMK, Chuang IL. Rev Mod Phys. 2004;76:1037.
288. Read EL, Engel GS, Calhoun TR, Mančal T, Ahn TK, Blankenship RE, Fleming GR. Proc Natl Acad Sci USA. 2007;104:14203. [PubMed]
289. Cheng YC, Engel GS, Fleming GR. Chem Phys. 2007;341:285.
290. Nagy A, Prokhorenko V, Miller DRJ. Curr Opin Struct Biol. 2006;16:654. [PubMed]
291. Albrecht A, Hybl J, Faeder S, Jonas D. J Chem Phys. 1999;111:10934.
292. Brumer P, Shapiro M. Chem Phys Lett. 1986;126:541.
293. Shapiro M, Hepburn JW, Brumer P. Chem Phys Lett. 1988;149:451.
294. Tannor DJ, Rice SA. J Chem Phys. 1985;83:5013.
295. Tannor DJ, Kosloff R, Rice SA. J Chem Phys. 1986;85:5805.
296. Chen C, Yin YY, Elliott DS. Phys Rev Lett. 1990;64:507. [PubMed]
297. Baumert T, Grosser M, Thalweiser R, Gerber G. Phys Rev Lett. 1991;67:3753. [PubMed]
298. Baumert T, Buehler B, Grosser M, Thalweiser R, Weiss V, Wiedenmann E, Gerber G. J Phys Chem. 1991;95:8103.
299. Potter ED, Herek JL, Petersen S, Liu Q, Zewail AH. Nature. 1992;355:66.
300. Weiner AM, Leaird DE, Wiederrecht GP, Nelson KA. Science. 1990;247:1317. [PubMed]
301. Weiner AM, Leaird DE, Patel JS, Wullert JR. IEEE J Quantum Electron. 1992;28:908.
302. Weiner AM. Rev Sci Instrum. 2000;71:1929.
303. Goswami D. Phys Rep. 2003;374:385.
304. Feurer T, Vaughan JC, Nelson KA. Science. 2003;299:374. [PubMed]
305. Brixner T, Krampert G, Niklaus P, Gerber G. Appl Phys B. 2002;74:S133.
306. Brixner T, Damrauer N, Krampert G, Niklaus P, Gerber G. J Opt Soc Am B. 2003;20:878.
307. Oron D, Dudovich N, Silberberg Y. Phys Rev Lett. 2003;90:213902. [PubMed]
308. Wu R, Sola IR, Rabitz H. Chem Phys Lett. 2004;400:469.
309. Silberberg Y. Nature. 2004;430:624. [PubMed]
310. Polachek L, Oron D, Silberberg Y. Opt Lett. 2007;31:631. [PubMed]
311. Plewicki M, Weber SM, Weise F, Lindinger A. Appl Phys B. 2007;86:259.
312. Plewicki M, Weise F, Weber SM, Lindinger A. Appl Opt. 2006;45:8354. [PubMed]
313. Judson RS, Rabitz H. Phys Rev Lett. 1992;68:1500. [PubMed]
314. Goldberg DE. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley; Reading, MA: 1989.
315. Bardeen CJ, Yakovlev VV, Wilson KR, Carpenter SD, Weber PM, Warren WS. Chem Phys Lett. 1997;280:151.
316. Yelin D, Meshulach D, Silberberg Y. Opt Lett. 1997;22:1793. [PubMed]
317. Baumert T, Brixner T, Seyfried V, Strehle M, Gerber G. Appl Phys B. 1997;65:779.
318. Zeidler D, Frey S, Kompa KL, Motzkus M. Phys Rev A. 2001;64:023420.
319. Brixner T, Pfeifer T, Gerber G, Wollenhaupt M, Baumert T. In: Femtosecond Laser Spectroscopy. Hannaford P, editor. Chapter 9. Springer; 2007. p. 225.
320. Suzuki T, Minemoto S, Kanai T, Sakai H. Phys Rev Lett. 2004;92:133005. [PubMed]
321. Brixner T, Krampert G, Pfeifer T, Selle R, Gerber G, Wollenhaupt M, Graefe O, Horn C, Liese D, Baumert T. Phys Rev Lett. 2004;92:208301. [PubMed]
322. Brixner T, Garcia de Abajo FJ, Schneider J, Pfeiffer W. Phys Rev Lett. 2005;95:093901. [PubMed]
323. Dudovich N, Oron D, Silberberg Y. Phys Rev Lett. 2004;92:103303. [PubMed]
324. Grumstrup EM, Shim SH, Montgomery MA, Damrauer NH, Zanni MT. Opt Express. 2007;15:16681. [PubMed]
325. Gundogdu K, Stone KW, Turner DB, Nelson KA. Chem Phys. 2007;341:89.
326. Vaughan JC, Hornung T, Stone KW, Nelson KA. J Phys Chem A. 2007;111:4873. [PubMed]
327. Hornung T, Vaughan JC, Feurer T, Nelson KA. Opt Lett. 2004;29:2052. [PubMed]
328. Brixner T. Appl Phys B. 2003;76:531.
329. Kim D, Osuka A. Acc Chem Res. 2004;37:735. [PubMed]
330. Hajjaj F, Yoon ZS, Yoon MC, Park J, Satake A, Kim D, Kobuke Y. J Am Chem Soc. 2006;128:4612. [PubMed]
331. Won Y, Friesner R, Johnson M, Sessler J. Photosynth Res. 1989;22:201. [PubMed]
332. Hacker M, Stobrawa G, Feurer T. Opt Express. 2001;9:191. [PubMed]
333. Kano H, Kobayashi T. J Chem Phys. 2002;116:184.
334. Zhuang W, Abramavicius D, Voronine DV, Mukamel S. Proc Natl Acad Sci USA. 2007;104:14233. [PubMed]
335. Joo T, Jia Y, Yu JY, Lang MJ, Fleming GR. J Chem Phys. 1996;104:6089.
336. Everitt KF, Geva E, Skinner JL. J Chem Phys. 2001;114:1326.
337. Piryatinski A, Chernyak V, Mukamel S. Chem Phys. 2001;266:311.
338. Ge NH, Zanni MT, Hochstrasser RM. Local Structure and Dynamics of Liquid Acetone by Heterodyned 2D IR Spectroscopy. In: Murnane MM, Scherer NF, Miller RJD, Weiner AM, editors. Ultrafast Phenomena XIII. Spring-Verlag; New York: 2002.
339. Ge NH, Zanni MT, Hochstrasser RM. J Phys Chem B. 2002;106:962.
340. Demirdoven N, Khalil M, Tokmakoff A. Phys Rev Lett. 2002;89:237401. [PubMed]
341. Demidöven N, Cheatum CM, Chung HS, Khalil M, Knoester J, Tokmakoff A. J Am Chem Soc. 2004;126:7981. [PubMed]
342. Fang C, Wang J, Kim YS, Charnley AK, Barber-Armstrong W, Smith AB, III, Decatur SM, Hochstrasser RM. J Phys Chem B. 2004;108:10415.
343. Tang J, Mei E, Green C, Kaplan J, DeGrado W, Smith A, III, Hochstrasser R. J Phys Chem B. 2004;108:15910.
344. Abramavicius D, Zhuang W, Mukamel S. J Phys Chem B. 2004;108:18034.
345. Hayashi T, la Cour Jansen T, Zhuang W, Mukamel S. J Phys Chem A. 2005;109:64. [PubMed]
346. Usui T. Prog Theor Phys. 1960;27:787.
347. Hanamura E. J Phys Soc Jpn. 1970;29:50.
348. Sheboul MI, Ekardt W. Phys Status Solidi B. 1976;73:165.
349. Hiroshima T. Phys Rev B. 1989;40:3862. [PubMed]
350. Chernyak V, Mukamel S. J Opt Soc Am B. 1996;13:1302.
351. Cam HN. J Russ Laser Res. 2004;25:412.
352. Combescot M, Betbeder-Matibet O, Combescot R. Phys Rev B. 2007;75:114305.
353. Combescot M, Pogosov W. Composite boson many-body theory for Frenkel excitons. Phys Rev B. submitted for publication.
354. Agranovich VM, Galanin MD. Electronic excitation energy transfer in condensed matter. North-Holland; Amsterdam, The Netherlands: 1982.
355. Bassani GF, Agranovich VM, editors. Electronic Excitations in Organic Based Nanostructures. Elsevier; New York: 2003.
356. Ueta M, Kanzaki H, Kobayashi K, Hanamura H. Excitonic processes in solids. Springer; Berlin: 1986.