Center Line Slope Analysis in Two-Dimensional Electronic Spectroscopy

Center line slope (CLS) analysis in 2D infrared spectroscopy has been extensively used to extract frequency–frequency correlation functions of vibrational transitions. We apply this concept to 2D electronic spectroscopy, where CLS is a measure of electronic gap fluctuations. The two domains, infrared and electronic, possess differences: In the infrared, the frequency fluctuations are classical, often slow and Gaussian. In contrast, electronic spectra are subject to fast spectral diffusion and affected by underdamped vibrational wavepackets in addition to Stokes shift. All these effects result in non-Gaussian peak profiles. Here, we extend CLS-analysis beyond Gaussian line shapes and test the developed methodology on a solvated molecule, zinc phthalocyanine. We find that CLS facilitates the interpretation of 2D electronic spectra by reducing their complexity to one dimension. In this way, CLS provides a highly sensitive measure of model parameters describing electronic–vibrational and electronic–solvent interaction.


INTRODUCTION
Ultrafast laser technology 1 in the last two decades progressed to the point where heterodyne-detected four wave mixing (FWM) experiments can fully characterize the third-order optical response. 2 In such experiments, three excitation pulses, separated by two time delays t 1 and t 2 , induce a coherent signal in a molecular sample, which is emitted during signal time t 3 . Such FWM signals are best displayed as ω 1 vs ω 3 correlation plots between absorption events during the coherence time t 1 and signal emissions during time t 3 . Such plots are called two-dimensional (2D) optical spectra 3−5 in close analogy to 2D methods used throughout the NMR community. 6 2D optical spectroscopy, first developed in the infrared (IR) domain, 7−10 has been brought to the visible 11 and other domains 12−16 over the past decade.
Analysis of 2D spectrograms is focused either on peak intensity and position or on the peak shape, where peak intensity and position are usually much simpler to interpret. For instance, at short waiting times t 2 ≈ 0 cross peak intensities bear information on the relative angle between the involved transition dipole moments. 2D IR was thus used to elucidate molecular 17 or protein 18,19 structure, as well as chemical exchange. 20,21 Similarly, the evolution of cross peak magnitudes are instrumental for tracking excitation and relaxation pathways in multichromophoric systems such as natural 22,23 and artificial light harvesting complexes. 24 Line shapes of individual peaks in 2D spectra are related to fluctuations of transition frequencies as caused by spectral diffusion processes. 25 In some instances, peak shapes reflect molecular structure or energy deactivation networks, 26 but more typically, 2D line shapes are heavily influenced by the environment of the probed molecular system. Although peak positions and intensities are explained quantitatively, detailed calculations of 2D line shapes are relatively costly and comparison between simulation and experiment is often left to visual inspection. This rationalizes the need for a simple quantitative characterization of 2D line shapes in the electronic domain. Ideally, such measures would allow for a simple parametrization of standard microscopic models of spectral diffusion. Moreover, these measures should also be well-defined for more complicated line shapes, i.e., for atypical parameter regimes or microscopic models.
There are two limiting cases for 2D peak shapes: (i) Fast and homogeneous spectral fluctuations induce Lorentzian star-like profiles, as typically found in 2D NMR. And (ii) slow Gaussian spectral diffusion produces characteristic 2D IR Gaussian peaks, the contours of which are approximately elliptic. 27 Elongation of such elliptic peaks along the diagonal is a measure for inhomogeneous disorder and its waiting time (t 2 ) evolution shows the extent of relaxation during spectral diffusion processes described by the frequency−frequency correlation function (FFCF). Several measures of 2D peaks have been introduced to characterize spectral diffusion, 28 such as contour eccentricity, nodal lines of the dispersive part of the spectrum or center lines slopes (CLS) of the absorptive parts. 29 The ratio of a peak's diagonal Δ D and antidiagonal Δ AD width is connected to the FFCF for Gaussian line shapes. 30 This measure is useful when the diagonal elongation is caused predominantly by static disorder. However, the concept has no reasonable extension beyond the very limit of slow Gaussian spectral diffusion, which is not always obeyed for electronic spectra and it is also difficult to identify if the 2D spectrum has deviated from a Gaussian line shape.
The nodal line in the dispersive spectrum is a rather phenomenological concept. 31 Nodal line dynamics reflect the relaxation of the line shape and can thus be used to infer relaxation time scales. However, we have no analytical expression relating the nodal line slope to other model parameters.
In contrast, the center line (CL) as obtained by maximization of a 2D signal along the ω 1 axis at fixed waiting times t 2 is a more robust concept. For slow, approximately Gaussian, 32 spectral diffusion, typical for vibrational spectra measured in the IR domain, the CL is linear and its slope manifests the frequency−frequency correlation function (FFCF). 29 In such a scenario, the dynamics of a Gaussian coordinate is entirely defined by the FFCF, and thus spectral diffusion is completely characterized by the dependence of CLS on the waiting time t 2 .
Dynamics of electronic transitions are often different. Intermediate time scales of spectral diffusion, significant bath reorganization effects such as Stokes shift and vibrational progression of electronic spectra induce complex, non-Gaussian line shapes. Consequently, the CL becomes nonlinear. However, we will demonstrate that it still carries information about dynamics and helps with parametrization of a spectral diffusion model.
In the present paper we will analyze the CLS concept in the context of 2D electronic spectroscopy (2D ES). After explaining all factors influencing CLS in 2D ES, we apply our findings on experimental data of zinc−phthalocyanine (ZnPc), 33 where the main diagonal peak shows different sections with linear or curved CL, a behavior rather dissimilar from that for the paradigmatic IR or NMR 2D spectra.
The paper is organized as follows. In section 2 we interpret 2D IR spectra as dynamical maps of a diffusive coordinate and review the concept of CLS. In section 3 we introduce spin− boson dynamics as a convenient model of Gaussian fluctuations of electronic transitions. It will be related to the stochastic picture of section 2 for a certain parametric limit. We will show how more general parametrization allows us to analyze CLS in the presence of typical phenomena for the electronic domain such as Stokes shift. We will also discuss limitations imposed by the spin−boson model and investigate the effects of non-Gaussian spectral diffusion in section 4. In section 5 we apply our models to experimental electronic 2D spectra of ZnPc. In section 6, we conclude by evaluating CLS analysis as a tool for parametrizing dynamics of electronic transitions.

CENTER LINE SLOPE OF GAUSSIAN 2D IR LINE
SHAPES We start this section with a short review of 2D spectroscopy. In FWM experiments, three short laser pulses impinge upon the sample at time delays t 1 and t 2 and the resultant third-order signal field generated at a delay t 3 is mixed with a local oscillator field. Consider a two-level chromophore, whose transition frequency ω(t) undergoes stochastic spectral fluctuations. The signal is a function of the three time intervals and reads 25 e +i corresponds to the so-called rephasing signal (R) and e −i to the nonrephasing signal (N). Angular brackets ⟨ ⟩ represent averaging over realizations of the stochastic process ω(t). Equation 1 shows clearly why ω 3 , t 2 , ω 1 is the standard domain for 2D spectroscopy. We will focus on the real part of the spectrum, e.g., The total (absorptive) signal R combines rephasing and nonrephasing contributions 35 In the infrared, spectral diffusion is typically a slow process leading to stable transition frequencies ω 1 and ω 3 during the t 1 a n d t 3 i n t e r v a l s . W e c a n t h e n a p p r o x i m a t e The total signal (eq 3) is then interpreted as two-time joint densities 36 of ω(t) Spectral fluctuations are most often of Gaussian type and entirely characterized by the two-point correlation function 32,37 Substituting the standard prescription for two-point joint density of a Gaussian process with vanishing mean ⟨ω(t)⟩ = 0 into eq 4 yields The Gaussian line shape R g delivers elliptical contours as shown in Figure 1 for several waiting times t 2 . Line shapes develop from highly eccentric contours at short waiting times t 2 ≈ 0 (left panel) to the circular limiting case for t 2 → ∞ (right panel) and depend on the extent of relaxation given by t ( ) 2 . The obvious task is to extract the correlation function t ( ) from the experimental 2D line shapes.
The black curve in each panel of Figure 1 depicts the center line 29 ω 1 BCL (ω 3 ). It was obtained by finding position ω 1 BCL of maximal signal R(ω 3 ,t 2 ,ω 1 ) along axis ω 1 at given ω 3 As far as we consider the 2D line shape to be a joint probabilistic density of a classical stochastic variable as in eq 4, ω 1 BCL (ω 3 ) is the most probable value of ω(t) at t = 0 given that the frequency at a later time t = t 2 is ω 3 . The center line thus represents backward-time evolution and we abbreviate this center line as BCL (backward-time center line). The BCL line of a Gaussian peak (eq 5) is linear because BCL with a slope given by the inverse of the normalized correlation function Alternatively, one can search for the maximal signal along the ω 3 axis. 34 The associated red line ω 3 FCL (ω 1 ) ( Figure 1) represents a typical spectral diffusion trajectory ω(t), which starts from ω(0) = ω 1 and has the most probable position ω(t 2 ) = ω 3 FCL (ω 1 ) after waiting time t 2 . Hence, ω 3 FCL (ω 1 ) represents the ordinary time-forward picture of a stochastic process and will be abbreviated FCL (forward-time center line). For a Gaussian peak (eq 5), FCL is again linear with a slope given directly by the normalized correlation function C(t) The symmetry between BCL and FCL evident from comparing eqs 6 and 7 has a fundamental background. The joint density of equilibrium stochastic processes is subjected to the microscopic time reversibility The line shapes are symmetric along the diagonal, BCL can be deduced from FCL and vice versa because they are axial images; FCL and BCL carry similar information in this case, as illustrated in Figure 1. Experimentally, it is inequivalent to determine the BCL and FCL, as ω 3 is typically measured directly in frequency domain via spectrally dispersed detection, whereas ω 1 is obtained in postprocessing from time-domain data. Routinely produced 2D spectra has often better resolved BCL. Equation 8 then provides a meaningful relation between experimentally accessible BCL and the ordinary forward-time picture of FCL. Beyond the assumptions underlying eq 4, however, the BCL and FCL may become dissimilar, so we continue to discuss both of them throughout this paper.
Here we note that the Gaussian joint density (eq 5) forms the backbone of most experimental line shape measures of the FFCF. For instance, the ratio of a peak's diagonal Δ D and antidiagonal width Δ AD within the limit of eq 5 is connected to C(t) as For Gaussian line shapes C(t) can thus be extracted similarly well from BCLS, FCLS or the eccentricity Δ AD /Δ D . However, eccentricity is limited to elliptical contours and cannot be easily generalized beyond Gaussian line shapes. In contrast, FCL and BCL allow the detection of non-Gausianities by deviations from linearity. In the following, we develop a theoretical framework for BCLS and FCLS, going beyond the limit of Gaussian line shapes.

SPECTRAL FEATURES OF ELECTRONIC TRANSITIONS
We now turn to 2D spectroscopy in the electronic domain. Although modulation of vibrational transition frequencies originates in the solvent and with few exceptions 38 can be considered classical, electronic transitions in molecular systems incorporate a more complex modulation of vibrations and show clear signatures of quantum behavior (e.g., Stokes shift).
We start with reinterpreting the transition frequency in eq 1 from the classical diffusive variable to a quantum coordinate ω(t) → ω eg + Q̂(t) and introduce the standard dynamical model of an electronic transition. Electronic transitions in molecules are always modulated by vibrations. By approximating them by a set of harmonic oscillators which are eventually displaced in the electronic excited state , we adopted the spin−boson Hamiltonian, which is well documented in literature. 39 Thus, each electronic level is accompanied by a handful of vibrational states forming a band, and there is a number of transitions between them which are overlapping and merge into a single peak (or few peaks) at room temperature. We thus understand the electronic 2D line shape as a map of a transition between a single electronic ground and excited level with the transition (gap) frequency dynamically modulated by vibrations. Within the standard Condon approximation 40,41 the transition frequency represents the difference between ground and excited state surfaces, i.e., ℏQ̂= Ĥe − Ĥg − ℏω eg = ∑ j m j Ω j 2 d j xĵ. Its time profile is found by switching into the Dirac picture Q̂(t) = e iĤgt/ℏ Q̂e −iĤgt/ℏ . It represents quantum Gaussian fluctuations at arbitrary time scales around the mean of the transition frequency ω eg .
In this quantum case the excited and the ground state dynamics are different and the total 2D signal of a single transition is defined by four Liouville space pathways depicted in Figure 2.
obeying the fluctuation−dissipation relation between its real t ( ) and imaginary ′ t ( ) components (see eq 47 for Fourier transform conventions) where β is the Boltzmann factor, i.e., the inverse temperature β = 1/kT. At high-temperatures ℏβω → 0 the real part dominates, corresponding to the classical case of eq 5. Spin− boson dynamics are solvable and the third-order response function can be calculated exactly using the second cumulant (eqs 42−45) that is reviewed in Appendix A. We next address several phenomena typical for electronic spectra that are rare for vibrational transitions in the infrared. First, potential surfaces of electronic ground and excited states are significantly different, and thus frequencies of emitted photons are lower than those absorbed. This Stokes shift commonly represents a simple displacement between ground and excited state harmonic potential surfaces. Differences between these surfaces, however, can be more dramatic. For instance, differences between the curvature of electronic ground and excited states induces nonlinear electronic−solvent coupling that results in non-Gaussian spectral diffusion, 42 which will be treated along with anharmonicity of potential surfaces in the section 4.
In many cases the electronic transition is coupled to underdamped vibrations that appear as side bands in the absorption spectrum in a vibronic progression. Moreover, this coupling modulates both the amplitude and peak shape of the primary peak as a function of time. This is another phenomenon which is rarely seen in the infrared.
Last, when the transition frequency is altered rapidly during t 1 and t 3 , Gaussian frequency fluctuations do not translate into Gaussian line shapes. Rather the line shapes are motionally narrowed and the slow fluctuation approximation introduced above (eq 4) fails. Although motional narrowing is not unknown in the infrared, 29 electronic transitions quite typically exhibit some of these signatures and thereby limit the use of the straightforward analysis of section 2.
Analysis within the present section will neglect the interference between levels of multilevel and multichromophoric systems, namely pathways of excited state absorption and cascading processes whose complex effects escape simple classification. Instead, we focus on a single transition between electronic ground and excited state and the effects of Stokes shift, finite fluctuation time-scales and vibrational structure on center lines will be addressed in the coming sections 3.1−3.3.
In the following we will outline quantum-classical correspondence and formally reproduce the classical case of Figure 1. In principle, the FFCF ℏ t ( ) can be chosen arbitrarily, only bounded by eq 13. The most common spectral diffusion model 2,39,43 (see Appendix A for details) describes overdamped quantum motion (diffusion) in a harmonic potential at moderately high temperatures (tanh(ℏβω/2) → ℏβω/2 in eq 13), with relaxation rate Λ, and coupling λ.
Model 14 thus refers to Figure 1 with magnitude of fluctuations The overdamped FFCF of eq 14 is capable of describing both the Stokes shift and effects of diffusion time scale. The Gaussian shape of eq 5 will be recovered in the slow fluctuation limit, which is specified by comparing the relaxation rate to the magnitude of fluctuations, i.e., Λ ≪ Δ. Following the procedure of ref 29 (see Appendix A), we approximate eqs 42−45 by where g is the auxiliary line broadening function The classical, real valued FFCF implies Im g(t) = 0. In this case, we recover the Gaussian line shape of eq 5 by insertion of eqs 16−19 into eqs 2 and 11. The correspondence between dynamics generated by the spin−boson Hamiltonian and classical Gaussian fluctuations is thus established. Our considerations will be demonstrated in the slow limit introduced by eqs 16−19. Line shapes of GSB pathways are Gaussians (eq 5) around transition frequency ω eg . The SE contributions are different due to the factor e −4i Img(t 2 )t 3 in eqs 17 and 19. The Gaussian line shapes defined by eq 5 are retained, but the peak is shifted along the ω 3 axis below the diagonal by −4 Im g(t 2 ): The total 2D signal (eq 11) is thus a combination of the two shifted Gaussian peaks. The waiting time evolution (Figure 3) shows two combined effects: (i) correlation loss similar to that in Figure 1 (i.e., change of contours from elongated ellipse to circle) and (ii) the development of Stokes shift for SE contributions. At short times t 2 ≈ 0, the Stokes shift Im g(t 2 ) is small with no apparent influence on shape. With increasing t 2 , the Stokes shift becomes significant (asymptotically Im g(∞) → −λ) and deforms the elliptic peak shapes. We now examine the extent to which the FFCF can be measured by CLS in the presence of the Stokes shift. Three distinct regions with different BCLS were uncovered in Figure  3, right panel. For ω 3 ≫ ω eg the line shape is dominated by GSB contribution and the BCLS follows For ω 3 ≪ ω eg − 4 Im g(t 2 ), the BCLS follows the SE contribution and thus BCL eg 1 2 This indicates that FFCF can be measured on the peak's periphery. Unfortunately, as will be demonstrated below, the peripheral BCLS is often distorted by motional narrowing or anharmonic effects. In the central region the relation between FFCF and BCLS is complex, because the Taylor expansion around the peak center suggests The BCLS is thus inapplicable as a direct measure of the FFCF in the presence of Stokes shift. We next turn to the FCLS, maximizing the signal along ω 3 . In a typical situation shown in the left and central panels of Figure 3 GSB and SE contributions are strongly overlapped. Their sum exhibits only a single maximum, located halfway between the maxima of GSB and SE signals The FCL is linear, shifted by −2 Im g(t 2 ) from the diagonal, but its slope remains unchanged following both the GSB and SE. We thus conclude that, from a theoretical point of view, the FCLS is a better measure for the FFCF. For large Stokes shift λ ≈ Δ the SE and GSB peaks start to separate, and the line shape along ω 3 is flat. The global maximum becomes ill defined, as demonstrated in the right panel of Figure 3 around ω 1 = 1.8Δ. For very large Stokes shift λ ≫ Δ SE and GSB peaks become separated, and center lines are better considered as local measures for each observed peak.

Effects of Rapid Spectral Diffusion.
Exact response functions (eqs 42−45) of the spin−boson model apply to both slow and fast spectral diffusion. Line shapes in the regime of fast spectral diffusion Λ ≫ Δ will be analyzed in the Λ → ∞ limit of eq 14, i.e., approximating the FFCF by The peaks in 2D electronic spectra are Lorentzian in this limit, with purely vertical (horizontal) BCL (FCL) and no waiting time dynamics, defined by

The Journal of Physical Chemistry
where Φ ≡ (2λ)/(ℏβΛ) and represents the rate of pure dephasing.
Finite fluctuation time scales shall be treated numerically, Δ/ Λ is varied in Figure 4. Contours in the intermediate regime Λ ∼ Δ change significantly from the edges to the center of the peak. The edges of a peak at ω − ω eg ≫ Λ reflect deviations far from equilibrium, where relaxation is fast and the contours assume star-like structures (eq 23) of fast (homogeneous) processes. The BCL is then aligned almost vertically to the ω 3 axis, even for slow spectral diffusion shown in the right panel of Figure 4. Near the center of the peak at ω − ω eg ≪ min(λ,Δ), the relaxation is slow and the Gaussian structures of eq 5 can be observed, with elliptic contours and with a BCL and FCL closer to the diagonal representing FFCF along with eqs 6 and 7 in Figure 1. The BCLS is sensitive to waiting time t 2 primarily in the central region where the BCLS can be deduced. Note that time scale effects do not modify time-reversal symmetry; line shapes with insignificant Stokes shift Im g(t) ≈ 0 obey a symmetry relation similar to eq 8 beyond the slow fluctuation limit. 44 The FCL thus still approximates a mirror image of the BCL.
The form of eq 23, with dephasing rate Φ as the only line shape-determining parameter, is reminiscent of effects of spontaneous emission on line shapes. Indeed, radiative dephasing Γ can be accounted for by adding a Gaussian coordinate with δ = Γ Γ t t ( ) ( ) to the correlation function. Numerical modeling of radiative dephasing is thus easily accomplished. However, the relation between CLS and FFCF (eq 6 or 7) becomes less straightforward. 29 Because the radiative rate Γ is often known, one can speculate that 2D spectra could be more easily interpreted, when the effect of Γ on line shapes is removed by post processing of the 2D data in the time domain, i.e., by multiplying the response function by e Γ(t 1 +t 3 ) . Solving practical difficulties of such a procedure is, however, left for future work.
3.3. Effects of Underdamped Vibrations. Electronic transitions are often modulated by underdamped vibrations. 23,33,45−47 Vibrations appear in the absorption spectrum as additional displaced peaks, referred to as vibrational (vibronic) progression. In 2D, besides a rich dynamical peak structure, vibrations also induce waiting time oscillatory dynamics of the principal peak, which will be the focus of this section.
We separate the environmental effects on an electronic transition into a solvent-related response R S and a vibrational response R V , e.g., and similarly for the other pathways of Figure 2. For the solvent response, 48,49 we assume the Gaussian line shapes of the slow where Ω is the vibrational frequency, λ V is vibrational coupling, and γ ≪ Ω is the damping rate. Equations 14, 24, and 25 fully define the model and were used to simulate the results shown in Figures 5 and 6. For a better understanding, we will analyze some important limits of this model. Simulated 2D electronic spectra depend on parametrization. We adopt a regime related to experiments on ZnPc, 33 where the principal peak at (ω eg , ω eg ) is spectrally resolved (Ω > Δ, λ V ) from peaks in the vibronic progression at (ω eg ± nΩ, ω eg ± mΩ). We can then expand the response function in powers of λ V /Ω and neglect components oscillating in t 1 and t 3 intervals ∝ sin(Ωt 1,3 ), cos(Ωt 1,3 ) . The only harmonic variation relevant for the principal peak occurs during the waiting time interval ∼cos Ωt 2 , sin Ωt 2 .
To first order in λ V /Ω, using eqs 42−45 and 49, we thus approximate We next use eq 24 and combine the vibrational contribution with the solvent response R S in the slow Λ ≪ Δ limit as expressed in eqs 16−19. Transformation into the frequency domain eq 2 is made with the use of the convolution theorem. For example, the nonrephasing contribution is a convolution of Gaussian line shape (eq 5) with the Fourier image of the characteristic function of the t 1 , where p.v. stands for the (Cauchy) principal value. After summing R and N signals and taking the real part (eqs 2 and 3), we obtain for GSB responses where ω̃i = ω i − (ω eg − 2λ V ). The SE response functions read

The Journal of Physical Chemistry
with ω̃1 = ω 1 − (ω eg − 2λ V ) and ω̊3 = ω 3 − (ω eg + 4 Im g(t 2 ) − 2λ V ). The standard Gaussian solvent line shape of absorptive nature represented by the first line in eq 26 is modulated by the second, harmonic term ∝ cos(Ωt 2 ) of dispersive nature. The ∝ sin(Ωt 2 ) terms cancel in the total signal, when the solvent Stokes shift vanishes, Im g(t 2 ) → 0, rendering it negligible in most realistic cases. The numerical simulations based on full cumulant expressions (eqs 42−45) are shown in Figure 5. Line shapes of the principal peak are similar to those in Figure 4, but the BCLS, FCLS is oscillating with the waiting time t 2 . In Figure 5 the BCLS of total, rephasing and nonrephasing signals were extracted from the central (linear) area of the peak. We noted three significant dynamical time scales: (i) the overall trend of BCLS represents solvent relaxation on a Λ −1 time scale ( Figure  5, top panel), (ii) vibrational oscillations of the BCLS on a Ω −1 time scale mainly derive from the nonrephasing part of signal which are (iii) damped away with a rate of γ −1 . In the present simulation we separated time scales in a realistic parametric regime of Ω ≫ γ > Λ (i.e., underdamped vibration, damped still before the solvent is relaxed). We shall next analyze the BCLS quantitatively and discuss the possibility of extracting Ω, γ, Λ, λ V , λ, etc. from experiments.
The effects of the oscillatory ∝ cos (Ωt 2 ) term on BCL are complicated, but they can be circumvented by averaging the signal over a period of 2π/Ω. The solvent FFCF can thus be measured by averaging the BCLS over the period from which the solvent parameters λ, Λ can be deduced according to section 2. The relation between the oscillations of the BCLS and the vibrational FFCF t ( ) V is complicated. We can straightforwardly measure the vibrational frequency Ω. In certain parametric regimes, an exponential decay of oscillatory amplitude may be observed and attributed to vibrational damping γ. Our simulations show (bottom panel of Figure 5), however, that such an approach for determining γ works surprisingly poorly at early t 2 times where the peaks are narrow and the BCL (FCL) is almost static. The oscillations first increase before being damped.
We have addressed this behavior in Appendix C, eq 52, where line shapes were analyzed around the peak center. Neglecting solvent Stokes shift, the BCLS has been approximated by  (27) and FCLS as its axial image. At short times, C(t 2 ) ∼ 1 − Λt 2 , the magnitude of oscillations is modulated by a singular prefactor ∝ t 2 3/2 , which thwarts direct extraction of γ using a simple fit of an exponential form e −γt 2 . Instead, using the form of eq 27, one can substantially improve the experimental determination of γ.
In Figure 6 we compare methods for deducing γ out of the 2D spectrum. One approach is based on CLS analysis. In particular, we will study oscillations of BCLS and its correction obtained by multiplying BCLS by t −3/2 to eliminate the singular factor ∝ t 3/2 (eq 27), as just disscused. The other approach is based on measuring peak volumes in electronic 2D spectra. In particular, we measured the oscillations of the volume of the principal peak (DP1) and for the first diagonal peak of the vibrational progression (DP2). 50 Oscillatory amplitudes shown at Figure 6 were defined as the difference between consecutive local maximum and minimum of the BCLS oscillatory curve of Figure 5. The same definition applies for the corrected BCLS and the peak volumes. A typical result is shown in Figure 6 in the semilog scale, where the exponential decays are linear lines. The exponential damping e −γt is plotted for comparison. We conclude that the peak volume of the higher lying diagonal peak DP2 and the corrected BCLS reproduce the correct decay rate, whereas the untreated CLS and the volume of DP1 fail to do so. We note that differences between the DP1 and DP2 peaks is due to the constructive (destructive) interference of R and N signals for DP2 (DP1). 50 The CLS correction factor retrieved from eq 27 has been shown to be essential for extracting the correct decay rate. We note that the need of the correction factor can be lesser beyond the slow fluctuation limit.
We also investigated the role of the parametric regime (λ V ≪ Ω, Λ ≪ Δ, γ ≪ Ω, Im ġ≈ 0) used in our analysis. We compared eqs 26 and 27 with the full simulations used in Figures 5 and 6. We found that the quadratic term (λ V /Ω) 2 should not be completely neglected in real situations (whenever the vibronic modulation is apparent), so one should not rely on eqs 26 and 27 quantitatively. However, all the features discussed above, such as the amplitude of CLS oscillations, their damping, and the short time t 3/2 modulation, are kept beyond the linear regime and are in fact even more clearly pronounced. Our analysis is thus qualitatively correct even far beyond formal validity.

ANHARMONIC SPECTRAL DIFFUSION AND NONLINEAR ELECTRONIC−VIBRATIONAL COUPLING
The elementary understanding of 2D line shapes outlined by eq 4 suggests additional explanations for the emergence of atypical (nonlinear) CL such as observed in ZnPc. A linear CL is a direct consequence of harmonic potential surfaces within the spin−boson model. Non-Gaussian spectral diffusion will result in a curved CL. 51 A simple example would be diffusion on an anharmonic surface. Another example emerges when the transition frequency is a nonlinear function of a harmonic The Journal of Physical Chemistry A Article coordinate. Both cases generalize beyond the spin−boson Hamiltonian, where a full quantum-mechanical treatment is difficult, and researchers mostly approximate by calculating few higher cumulants. 42 Truncating cumulant expansions, however, tends to distort line shapes heavily. We thus prefer to approximate the 2D line shape by a joint distribution of electronic transition frequencies (eq 28) following a classical diffusive coordinate Q. 4.1. Classical Diffusive Coordinate on Potential Surfaces. A two-level system modulated by a classical diffusive anharmonic coordinate should be analyzed by means of stochastic quantum dynamics (SQD), which has been elaborated in detail elsewhere. 52−54 To bring SQD to finite temperatures, the ground state ℏV g (GSB diagram) and the excited state ℏV e (SE diagram) state potential surfaces for the spectral diffusion are allowed to be different. This modification of classical SQD was introduced by ref 55 and, for harmonic surfaces, it is equivalent to the Kubo−Tanimura hierarchy. 56,57 Within the slow diffusion limit (eq 4) the GSB can be approximated as a joint probability of transition frequency following equilibrium spectral diffusion on the ground state potential surface The SE contribution is represented by a joint distribution of transition frequencies following nonequilibrium spectral diffusion on the excited state potential The diffusion of a classical coordinate Q on the potential surface is described by a Smoluchowski equation 58 (30) where β is reciprocal temperature, and D the diffusion constant. The potential surface ℏV i is different for the electronic ground V g and excited states V e . Diffusion can be equivalently modeled by an asymmetric random walk 59 with step length Δ Dt 2 , sampling time t Δ and probabilities to step forward which is the ground state Boltzmann distribution for both GSB and SE diagrams. For the same reason, the Franck−Condon principle, the transition frequency ω ≡ ω(Q) reflects gap between electronic ground and excited state The joint distribution of the Q coordinates ′ Q Q t ( , , ) i Q should be transformed into frequency joint distributions ω ω′ t ( , , ) (33) 4.2. Relation to Overdamped Quantum Brownian Motion. In the simplest case the Q coordinate is linear to the fluctuations of the transition frequency and can be rescaled to represent them directly

Equation 33 then reads
eg . The overdamped Brownian oscillator correlation function (eq 14) used corresponds to diffusion on harmonic potential surfaces where α is the force constant related to coupling in eq 14 by λ = 1/(2α). The joint distribution (eq 30) for the GSB is a Gaussian The maximum of the Gaussian packet (eq 36) Q g (t) started from Q g (t=0) = Q(0) and approaches the center of potential at Q = 0 as Q g (t) = e −2ℏβαDt Q(0). This shows directly that the FCL is linear with the slope e −2ℏβαDt . The time-reversal symmetry (eq 8) guarantees for the 2D line shape axial symmetry along the diagonal R GSB (ω 1 ,t 2 ,ω 3 ) = R GSB (ω 3 ,t 2 ,ω 1 ). The BCLS thus assumes the slope e 2ℏβαDt , in perfect agreement with the results of section 2. Similar analysis of diffusion on excited state surface recovers line shapes of SE pathway.
4.3. Spectral Difussion in Anharmonic Potential. We next discuss CL in the presence of the anharmonic potential surfaces. As the additional effect of Stokes shift would not affect our conclusions on anharmonic spectral diffusion, we adopt the high temperature limit where the Stokes shift becomes much smaller than the peak width λ Δ = → The above potential is in part harmonic, but the linear back force has different force constants α ± for positive and negative Q. Far from the center, the line shapes are approximately Gaussian and related to Figure 1 by mapping Λ ± = 2α ± ℏβD and Δ ± = βα ℏ ± 1/ 2 for low (−) and (+) high frequencies, respectively. FCL follows a typical trajectory Q(t) = e −2α + ℏβDt Q(0) for Q ≫ 0, and Q(t) = e −2α − ℏβDt Q(0) for Q ≪ The Journal of Physical Chemistry A Article 0, where D is the diffusion constant. Figure 7 shows simulation results of the diffusion model outlined by eqs 30, 31, and 37 with the predicted asymptotic properties, i.e., a CL with two linear sections and a transition region around ω ∼ ω eg .
Similar behavior can be observed in Figure 8 for a more realistic and smooth form of anharmonicity, which is achieved after adding a cubic term to the harmonic ground state potential The curvature of BCL is proportional to α 3 and evolves with waiting time. We thus conclude that spectral diffusion in an anharmonic potential is represented by a curved CL in electronic 2D spectra.
The transition frequency in the Condon approximation (eq 32) becomes a nonlinear function of the harmonic coordinate Q The Q-coordinate during diffusion on both excited and ground state is Gaussian, however, the transition frequency is not. A nonlinear transformation, eq 40, between Gaussian Q and non-Gaussian ω distributions must be simulated with full use of the general rule, eq 33. The resulting 2D spectra are shown in Figure 9, where we observe a moderate curvature of the CL that is t 2 dependent and ∝ α e − α g proportional. We note that in this model the CL is visually only moderately curved unless differences between α e and α g are rather large (>10%). These results can be understood by an elementary analysis of the model. The diffusive motion of the Gaussian coordinate Q still follows Q g (t) = e −2α g ℏβDt Q g (0) in the ground state (eq 36). The initial position Q(0) shall be obtained from inverting eq 40 where we assumed 4(ω − ω eg )(α e − α g ) < 1. In the present subsection ≈ stands for Taylor expansion to second order in ω(0) − ω eg . The maxima of Q and ω distributions may be slightly different due to the Jacobian in eq 33, but the difference is usually small and we shall neglect it for the present analysis. The frequency evolution in ground state is then and the FCL is represented by graphing ω(0) vs ω g (t). The BCL can be obtained by axial symmetry (eq 8).
The joint probability for SE pathways (eq 30) is The center of this Gaussian distribution approaches the excited state potential minimum as  GSB and SE line shapes are usually similar, being only slightly shifted along ω 3 . We can thus conclude that CL of eq 39 are curved ∝ α e − α g , which was confirmed by the full simulation shown in Figure 9.
5. CENTER LINE OF 2D SPECTRA OF ZINC−PHTHALOCYANINE As a simple test case for the methodology developed in the previous sections, we analyze the 2D spectra of zinc− phthalocyanine (ZnPc). ZnPc is a rigid, planar, and square symmetric (D 4h ) molecule (see Figure 10 for its molecular structure). The absorption spectrum (bottom panel) shows a narrow electronic transition at 14 850 cm −1 . 60 Peaks of lower intensity shifted by Ω a ≈ 700 cm −1 and Ω b ≈ 1600 cm −1 are readily identified as vibronic side bands. 61 The experimental details of the 2D ES measurements of ZnPc ( Figure 10) have been published previously. 33 Briefly, excitation pulses tunable throughout the visible spectral range are provided by a home-built noncollinear optical parametric amplifier, 62 pumped by a regenerative titanium-sapphire amplifier system (RegA 9050, Coherent Inc.) at 200 kHz repetition rate. Pulse spectra were chosen to overlap with ZnPc's absorption spectrum (Figure 10, lower panel) and compressed to a width of sub-8 fs that was determined using intensity autocorrelation. The pulses were attenuated by a neutral density filter to yield 8.5 nJ per excitation pulse at the sample. This corresponds to a fluence of less than 3.0 × 10 14 The Journal of Physical Chemistry A Article photons/cm 2 per pulse. The setup employed for 2D ES relies on a passively phase stabilized setup with a transmission grating 11,63 and has a temporal resolution of 0.67 fs for t 1 , and 5.3 fs for t 2 . A detailed description was given in ref 64. The emerging third-order signal was spectrally resolved in ω 3 by a grating-based spectrograph and recorded with a CCD camera. At given t 1 and t 2 delays, spectra were recorded by integration over approximately 10 5 shots per spectrum.
Sample circulation was accomplished by a wire-guided drop jet 65 with a flow rate of 20 mL/min and a film thickness of approximately 180 μm. ZnPc was dissolved in benzonitrile, and the concentration was set to obtain an optical density of 0.13 OD measured directly in the flowing sample jet. All measurements were performed under ambient temperatures (295 K).
An example 2D spectrum is shown at delay time t 2 = 96 fs in Figure 10. The signal has positive and negative features, the latter indicating excited state absorption (ESA). 57 We therefore consider an energy level structure involving three states, i.e., electronic ground, one first excited and one doubly excited state. Vibrational modes identified in the absorption spectra are considered as two underdamped harmonic oscillators as described in section 3.3.
We next analyze the shape and the center line of the principal peak DP1 of the ZnPc spectrum ( Figure 10) using the models presented in the previous sections. Focusing on the BCL (black line), the main experimental observations are as follows: (i) The center line is curved and the slopes within the high and low frequency regions are somewhat different. (ii) The larger slope in low frequency periphery can be attributed to motional narrowing described in sub section 3.2. In the high frequency periphery no signatures of motional narrowing were apparent due to interference with other peaks such as the ESA. (iii) The slope shows fast periodic modulation on time scales similar to the period of 700 cm −1 vibrational mode, in line with the discussion of sub section 3.3. (iv) The effect of the Stokes shift Figure 11. Comparison of experimental (top row) and simulated (bottom row) 2D ES of ZnPc at waiting times t 2 = 24−96 fs. Simulations assume a three-level chromophore modulated by an overdamped solvent mode (eq 14) and two underdamped vibrational modes (eq 25). The color scale is the same as in Figure 10. Parameters for simulation were retrieved as follows: Ω a = 700 cm −1 , Ω b = 1600 cm −1 , ω eg = 15000 cm −1 . Solvent parameters: λ = 80 cm −1 , Λ = 5 cm −1 . Vibrational parameters: λ V,a = 40 cm −1 , λ V,b = 160 cm −1 , γ a = 10 cm −1 , γ b = 150 cm −1 . Parameters for ESA pathway (eq 46) ω ef = 15300 cm −1 , g ff = g ee . Homogenous dephasing Γ = 30 cm −1 . Temperature T = 300 K. The Journal of Physical Chemistry A Article on the CLS of DP1, as discussed in section 3.1, is negligible given its relatively small value of 175 cm −1 and its large associated time scale of up to 2.5 ps 66 compared to 96 fs of our experiment. And finally (v) Quantum chemical studies suggest that anharmonicity of the involved potential energy curves is negligible for ZnPc 33,60 and thus anharmonic corrections discussed in section 4 are not relevant to this discussion. Thus, we reach a qualitative understanding of ZnPc using the BCL.
Given the multitude of influencing factors stated above, no simple expression for BCL of FCL can be given for ZnPc. Instead, both quantities must be retrieved from modeled spectra. In an effort to quantitatively describe the experimental ZnPc spectra, we therefore implemented the spin−boson model with two underdamped oscillators with V prescribed as eq 25; ESA contributions were included using eq 46 and homogeneous dephasing 67 as the fast component of relaxation δ = Γ Γ t ( ) . The conventional approach to evaluate the quality of a fit between model and experiment is to visually compare simulated and measured 2D spectra. Such a qualitative assessment can be made in Figure 11, where experimental and simulated 2D spectra of ZnPc are plotted. The experimental data were collected with a waiting time step of 12 fs from 0 to 96 fs covering two periods of a slower (700 cm −1 ) vibrational mode with period of 48 fs. Visually, the experimental and simulated DP1 peak seems quite similar and the CL also compares well. Slight differences between the CL's can be observed in the low frequency periphery of DP1 and are attributed to the lower stability of CL in peripheral regions as well as interference from a low energy cross peak in the simulation that is not present or, at best, not well resolved in the experimental data. As described, visual inspection of 2D plots can often be misleading; we thus focus on center lines and slopes as these lower-dimensional objects are easier to compare. Figure 12 compares the slopes of BCL for simulated and experimental data. The BCLS were measured by linear regression using only the central part of BCL for both experimental and simulated 2D spectra. It was found that the central part is more linear than the periphery and thus better for regression analysis without practical difficulties. For ZnPc, practical values for the linear region of the BCL slope are for maxima determined for peak values greater than 0.7−0.8 of the peak intensity (with respect to DP1 maximum).
The CLS dynamics for ZnPc within the experimental time window show a dominant oscillatory modulation with the expected period of 48 fs associated with the underdamped 700 cm −1 mode. The simulation parameters, ω eg , Ω a , and Ω b were estimated from the absorption spectrum. Global properties of 2D spectra (such as antidiagonal width of peaks) determined homogeneous rate Γ. Couplings and vibrational relaxation time scales (λ, λ V , γ) were obtained by fitting the BCLS dynamics in Figure 12. Solvent relaxation was determined to be negligible (Λ < 5 cm −1 ) though due to the temporal resolution of the experiment, this estimate is uncertain. Note that the retrieved coupling constant of the dominant mode λ V,a = 40 cm −1 , fitted form BCLS corresponds to the Huang−Rhys factors , i.e, λ V,a = 39.2 cm −1 reported from analysis of the absorption spectrum. 61 Generally, the experimental and simulated BCLS show similar trends; the discrepancy at early waiting times arguably result from pulse overlap artifacts that were not included in the model. We note that moderate changes of the model parameters had a pronounced effect on BCLS. This is in contrast to visual control of 2D peaks, where only rather extreme changes are visually recognizable. Comparison via CLSs (Figure 12) is thus a better indicator of successful simulation than comparison of 2D plots.

CONCLUSIONS
Here we summarize the merits of CLS analysis of electronic 2D spectra. Compared to 2D IR spectra, where the CLS is often a straightforward measure of the FFCF, our interpretation of typical electronic transition met several obstacles. We analyzed physical processes and regimes typical for electronic transitions such as Stokes shift, finite fluctuation time scales, and anharmonicity of electronic surfaces. These processes result in non-Gaussian line shapes, where the linear character of the CL breaks down, and thereby result in a less direct relation between CLS and FFCF. We still, however, identified specific regions of the CL, where the FFCF can be extracted in a simple manner and in other cases connected the geometric characteristics of the CL with microscopic parametrizations. We also highlighted differences between BCL and FCL and their dynamical implications.
To assess the applicability of CL, we analyzed experimental data of ZnPc, where 2D ES showed atypical curved (or multilinear) BCLS. We explained this feature using the spin− boson Hamiltonian in a parametric regime where (i) motional narrowing affects the low frequency periphery of the peak and (ii) the high frequency region is affected by interference with ESA pathways. CLS analysis proved to be a sensitive test for line shape oscillations originating from underdamped vibrations. These oscillations were superimposed on a slow exponential decay of CLS caused by solvent spectral diffusion. This provided a correct measure of vibrational damping and other parameters that would otherwise be difficult to access. Although CL analysis of electronic 2D spectra is certainly less straightforward and less powerful than in the IR, the method presented here provides a benchmark to estimate the quality of simulated data. CL analysis in 2D ES thereby helps to characterize all system parameters, even for electronic transitions, and serves as a one-dimensional reduction of complex two-dimensional electronic spectra and thereby facilitates the quality assessment of simulated data.

■ A. CUMULANT SOLUTION TO SPIN−BOSON HAMILTONIAN
In this Appendix we review the response for the spin−boson Hamiltonian, a standard microscopic representation of quantum Gaussian fluctuation. Its optical response functions of all orders can be calculated exactly using second cumulant. 2 Second cumulants are also a common approximation even beyond the spin−boson model linking their FFCF with nonlinear response. They were used extensively for describing a whole class of older spectroscopic techniques for measurement ot the FFCF such as photon echo peak shift. 48 Equations 42−45 are used in simulations. Expanding g to second order in t 1,3 , e.g., g(t 1 + t 2 ) ≈ g(t 2 ) + g(t 2 ) t 1 + g(t 2 )t 1 2 /2 and neglecting the imaginary part of the quadratic terms =* = g t g t t ( ) ( ) ( ) we recover the slow limit and obtain eqs 16−19.
The second cumulant is also used to calculate excited state absorption required to simulate the ZnPc spectra in section 4 After substituting into eq 51, we finally arrive at The Journal of Physical Chemistry A Article ω ω