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 t1 and t2, induce a coherent signal in a molecular sample, which is emitted during signal time t3. Such FWM signals are best displayed as ω1 vs ω3 correlation plots between absorption events during the coherence time t1 and signal emissions during time t3. 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 t2 ≈ 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 (t2) 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 t2 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 t2.

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 we interpret 2D IR spectra as dynamical maps of a diffusive coordinate and review the concept of CLS. In section 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 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 . In section we apply our models to experimental electronic 2D spectra of ZnPc. In section , 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 t1 and t2 and the resultant third-order signal field generated at a delay t3 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 shows clearly why ω3, t2, ω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 t1 and t3 intervals. We can then approximate and . The total signal (eq ) 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], where Δ̃ω(t) ≡ ω(t) – ⟨ω(t)⟩. Substituting the standard prescription for two-point joint density of a Gaussian process with vanishing mean ⟨ω(t)⟩ = 0 into eq yieldsThe Gaussian line shape Rg delivers elliptical contours as shown in Figure for several waiting times t2. Line shapes develop from highly eccentric contours at short waiting times t2 ≈ 0 (left panel) to the circular limiting case for t2 → ∞ (right panel) and depend on the extent of relaxation given by . The obvious task is to extract the correlation function from the experimental 2D line shapes.

Figure 1

Center lines BCL (black) and FCL (red) for Gaussian 2D line shapes (eq ) with , and Δ/Λ = 4000. Waiting time increases from left to right Λt2 = 0.01, 1, 10. Throughout the text, 2D plots are normalized to maximal value and for Figures , 3, 4, 7–9 contour lines are plotted at 5% steps, and the same color scale applies.

Figure 3

Effects of the time dependent (Λt2 = 0.1, 0.5, 1 from left to right) Stokes shift on BCL (black) and FCL (red). Parameters are λ/Δ = 1.5, ℏβΛ = 0.1.

Figure 4

Interpolating from slow bath (Δ/Λ = 10, left panel) to fast bath limit (Δ/Λ = 0.5 in the right panel) at waiting time Λt2 = 0.5, with ℏβΛ = 0.001.

Figure 7

Waiting time evolution of a 2D-signal for a potential that is harmonic in parts, eq , shown at times Λ–t2 = 0.1, 0.4, 0.8 for parameters α+/α– = Λ+/Λ– = 0.2. We used .

Figure 9

Effect of nonlinear coupling (eq ) when the ground and excited electronic potentials have different curvatures αg/αe = Λg/Λe = 1.03. Line shapes are shown at Λet2 = 0.2, 0.5, 1.0. The ground state width is defined as . The shift between ground and excited state is characterized by 1/(2α) = 1.46Δ.

The black curve in each panel of Figure depicts the center line[29] ω1BCL(ω3). It was obtained by finding position ω1BCL of maximal signal R(ω3,t2,ω1) along axis ω1 at given ω3As far as we consider the 2D line shape to be a joint probabilistic density of a classical stochastic variable as in eq , ω1BCL(ω3) is the most probable value of ω(t) at t = 0 given that the frequency at a later time t = t2 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 ) is linear becausewith 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 ω3FCL(ω1) (Figure 1) represents a typical spectral diffusion trajectory ω(t), which starts from ω(0) = ω1 and has the most probable position ω(t2) = ω3FCL(ω1) after waiting time t2. Hence, ω3FCL(ω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 ), 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 and 7 has a fundamental background. The joint density of equilibrium stochastic processes is subjected to the microscopic time reversibilityThe 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 . 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 then provides a meaningful relation between experimentally accessible BCL and the ordinary forward-time picture of FCL. Beyond the assumptions underlying eq , 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 ) 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 is connected to C(t) asFor 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 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 oscillatorswhich 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 = ∑mΩ2dx̂. Its time profile is found by switching into the Dirac picture Q̂(t) = eiQ̂e–i. 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 .where ground state bleaching diagrams (GSB) propagate within the electronic ground state during waiting time t2, and stimulated emission diagrams (SE) through the excited state.

Figure 2

Liouville space pathways for the third-order responses of a two-level chromophore.

As in the classical case, all properties of quantum Gaussian noise are described by the FFCF, which is now a complex-valued quantityobeying the fluctuation–dissipation relation between its real and imaginary components (see eq 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 . Spin–boson dynamics are solvable and the third-order response function can be calculated exactly using the second cumulant (eqs –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 .

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 t1 and t3, 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 ) 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 .

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

In the following we will outline quantum-classical correspondence and formally reproduce the classical case of Figure . In principle, the FFCF can be chosen arbitrarily, only bounded by eq . 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 ), with relaxation rate Λ, and coupling λ.Model 14 thus refers to Figure with magnitude of fluctuationsThe overdamped FFCF of eq is capable of describing both the Stokes shift and effects of diffusion time scale.

The Gaussian shape of eq 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 –45 bywhere g is the auxiliary line broadening functionThe classical, real valued FFCF implies Im ġ(t) = 0. In this case, we recover the Gaussian line shape of eq by insertion of eqs –19 into eqs and 11. The correspondence between dynamics generated by the spin–boson Hamiltonian and classical Gaussian fluctuations is thus established.

Effects of Stokes Shift

The Stokes shift is manifested in SE pathways. Here, the photon emission follows the evolution during waiting time t2 on the excited state potential surface. According to the Franck–Condon principle[40,41] the waiting time evolution of Q̂ starts from the ground state equilibrium. If the excited state equilibrium is different, microscopic reversibility (eq ) is violated and the energy is dissipated, resulting in the Stokes shift. In contrast, the GSB pathways exhibit similar frequencies of absorbed and emitted photons. The ground state evolution obeys eq and remains much closer to the stochastic line shapes of section .

Our considerations will be demonstrated in the slow limit introduced by eqs –19. Line shapes of GSB pathwaysare Gaussians (eq ) around transition frequency ωeg.

The SE contributions are different due to the factor e–4i Im in eqs and 19. The Gaussian line shapes defined by eq are retained, but the peak is shifted along the ω3 axis below the diagonal by −4 Im ġ(t2):The total 2D signal (eq ) is thus a combination of the two shifted Gaussian peaks. The waiting time evolution (Figure ) shows two combined effects: (i) correlation loss similar to that in Figure (i.e., change of contours from elongated ellipse to circle) and (ii) the development of Stokes shift for SE contributions. At short times t2 ≈ 0, the Stokes shift Im ġ(t2) is small with no apparent influence on shape. With increasing t2, the Stokes shift becomes significant (asymptotically Im ġ(∞) → −λ) 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 , right panel. For ω3 ≫ ωeg the line shape is dominated by GSB contribution and the BCLS followsFor ω3 ≪ ωeg – 4 Im ġ(t2), the BCLS follows the SE contribution and thusThis 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 suggestsThe 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 GSB and SE contributions are strongly overlapped. Their sum exhibits only a single maximum, located halfway between the maxima of GSB and SE signalsThe FCL is linear, shifted by −2 Im ġ(t2) 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 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 –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 , 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 bywhere Φ ≡ (2λ)/(ℏβΛ) and represents the rate of pure dephasing.

Finite fluctuation time scales shall be treated numerically, Δ/Λ is varied in Figure . 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 ) 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 . Near the center of the peak at ω – ωeg ≪ min(λ,Δ), the relaxation is slow and the Gaussian structures of eq can be observed, with elliptic contours and with a BCL and FCL closer to the diagonal representing FFCF along with eqs and 7 in Figure . The BCLS is sensitive to waiting time t2 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 ġ(t) ≈ 0 obey a symmetry relation similar to eq beyond the slow fluctuation limit.[44] The FCL thus still approximates a mirror image of the BCL.

The form of eq , 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 to the correlation function. Numerical modeling of radiative dephasing is thus easily accomplished. However, the relation between CLS and FFCF (eq 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Γ(. Solving practical difficulties of such a procedure is, however, left for future work.

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 RS and a vibrational response RV, e.g.,and similarly for the other pathways of Figure . For the solvent response,[48,49] we assume the Gaussian line shapes of the slow diffusion limit eq –19, i.e., RS ≡ R determined by the FFCF of eq . The response of an underdamped vibration RV can still be calculated in the framework of the spin–boson Hamiltonian; i.e., we can use the cumulant expressions eqs –45, but with the FFCF of a weakly damped harmonic oscillator (see Appendix B) approximated bywhere Ω is the vibrational frequency, λV is vibrational coupling, and γ ≪ Ω is the damping rate. Equations , 24, and 25 fully define the model and were used to simulate the results shown in Figures and 6. For a better understanding, we will analyze some important limits of this model.

Figure 5

Waiting time evolution of center lines for a vibrationally modulated transition as described in eq in the parametric regime consistent with eq . Top panel: evolution of BCLS for a vibrationally modulated two-level system at the solvent relaxation time scale Λ–1. The rephasing, nonrephasing, and total absorptive signal are shown in blue, red, and black, respectively. Bottom panel: the same at the vibrational relaxation time scale Ω–1. Other parameters are Δ/Λ = 100, γ/Λ = 3.75, Δ/Ω = 0.1, λV/Ω = 0.0025, ℏβΩ = 0.01. Simulations were carried out using the full cumulant solution eqs , 25, 50, and 42–45.

Figure 6

Damping rate γ retrieved by different measures. Exponential decays appear as straight lines in the linear-log scale. Black line: vibrational decay (to be retrieved) e–. Black dots: BCLS. Blue dots: BCLS with correction factor t23/2 as explained in the text. Green dashed line: peak volume of principal (DP1) peak. Red dashed line: peak volume of DP2 peak (diagonal peak representing first vibration). Parameters are Δ/Λ = 100, γ/Λ = 2.5, Δ/Ω = 0.1, λV/Ω = 0.0025, ℏβΩ = 0.01. Oscillatory amplitudes were retrieved from the difference between consecutive local maxima and minima along t2, followed by normalization at γt2 = 0.52 to e–γ.

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 t1 and t3 intervals ∝ sin(Ωt1,3), cos(Ωt1,3) . The only harmonic variation relevant for the principal peak occurs during the waiting time interval ∼cos Ωt2, sin Ωt2.

To first order in λV/Ω, using eqs –45 and 49, we thus approximate

We next use eq and combine the vibrational contribution with the solvent response RS in the slow Λ ≪ Δ limit as expressed in eqs –19. Transformation into the frequency domain eq is made with the use of the convolution theorem. For example, the nonrephasing contribution is a convolution of Gaussian line shape (eq ) with the Fourier image of the characteristic function of the t1, t3 > 0 quadrant where p.v. stands for the (Cauchy) principal value. After summing R and N signals and taking the real part (eqs and 3), we obtain for GSB responseswhere ω̃ = ω – (ωeg – 2λV). The SE response functions readwith ω̃1 = ω1 – (ωeg – 2λV) and ω̊3 = ω3 – (ωeg + 4 Im ġ(t2) – 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(Ωt2) of dispersive nature. The ∝ sin(Ωt2) terms cancel in the total signal, when the solvent Stokes shift vanishes, Im ġ(t2) → 0, rendering it negligible in most realistic cases.

The numerical simulations based on full cumulant expressions (eqs –45) are shown in Figure . Line shapes of the principal peak are similar to those in Figure , but the BCLS, FCLS is oscillating with the waiting time t2. In Figure 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 , 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 (Ωt2) 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 .

The relation between the oscillations of the BCLS and the vibrational FFCF 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 ), however, that such an approach for determining γ works surprisingly poorly at early t2 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 , where line shapes were analyzed around the peak center. Neglecting solvent Stokes shift, the BCLS has been approximated byand FCLS as its axial image. At short times, C(t2) ∼ 1 – Λt2, the magnitude of oscillations is modulated by a singular prefactor ∝ t23/2, which thwarts direct extraction of γ using a simple fit of an exponential form e–. Instead, using the form of eq , one can substantially improve the experimental determination of γ.

In Figure 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 ∝ t3/2 (eq ), 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 were defined as the difference between consecutive local maximum and minimum of the BCLS oscillatory curve of Figure . The same definition applies for the corrected BCLS and the peak volumes. A typical result is shown in Figure in the semilog scale, where the exponential decays are linear lines. The exponential damping e– 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 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 and 27 with the full simulations used in Figures 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 and 27 quantitatively. However, all the features discussed above, such as the amplitude of CLS oscillations, their damping, and the short time t3/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 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 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 ) following a classical diffusive coordinate Q.

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 ℏVg (GSB diagram) and the excited state ℏVe (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 ) the GSB can be approximated as a joint probability of transition frequency following equilibrium spectral diffusion on the ground state potential surfaceThe 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]where β is reciprocal temperature, and D the diffusion constant. The potential surface ℏV is different for the electronic ground Vg and excited states Ve. Diffusion can be equivalently modeled by an asymmetric random walk[59] with step length , sampling time tΔ and probabilities to step forwardand backward

Coordinate Q is not changed during excitation (Franck–Condon principle). The Smoluchowski equation (eq ) thus should be solved with the initial conditionwhich 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 should be transformed into frequency joint distributions using

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 directlyEquation then reads . The overdamped Brownian oscillator correlation function (eq ) used corresponds to diffusion on harmonic potential surfaceswhere α is the force constant related to coupling in eq by λ = 1/(2α). The joint distribution (eq ) for the GSB is a Gaussianwith FFCF and obeys eq . Correspondence to Figure is established with Λ ≡ 2ℏβαD and Δ2 ≡ 1/(2ℏβα).

The maximum of the Gaussian packet (eq ) Qg(t) started from Qg(t=0) = Q(0) and approaches the center of potential at Q = 0 as Qg(t) = e–2ℏQ(0). This shows directly that the FCL is linear with the slope e–2ℏ. The time-reversal symmetry (eq ) guarantees for the 2D line shape axial symmetry along the diagonal RGSB(ω1,t2,ω3) = RGSB(ω3,t2,ω1). The BCLS thus assumes the slope e2ℏ, in perfect agreement with the results of section . Similar analysis of diffusion on excited state surface recovers line shapes of SE pathway.

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 and thus RSE ≈ RGSB in the following simulations. We start with a simple form of an anharmonic potential

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 by mapping Λ± = 2α±ℏβD and Δ± = for low (−) and (+) high frequencies, respectively. FCL follows a typical trajectory Q(t) = e–2αQ(0) for Q ≫ 0, and Q(t) = e–2αQ(0) for Q ≪ 0, where D is the diffusion constant. Figure shows simulation results of the diffusion model outlined by eqs , 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 for a more realistic and smooth form of anharmonicity, which is achieved after adding a cubic term to the harmonic ground state potentialThe 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.

Figure 8

Waiting time evolution of a 2D-signal for a cubic anharmonic potential, eq , at Λt2 = 0.1, 0.4, and 0.8. The extent of anharmonicity is characterized by the ratio of harmonic and cubic terms at the equilibrium half-width = 3.75. Parameters Λ and Δ were defined by eq .

Nonlinear Electronic–Vibrational Coupling

We shall now analyze the CL when the ground and excited state surfaces are harmonic but have a different curvature αe ≠ αgThe transition frequency in the Condon approximation (eq ) becomes a nonlinear function of the harmonic coordinate QThe Q-coordinate during diffusion on both excited and ground state is Gaussian, however, the transition frequency is not. A nonlinear transformation, eq , between Gaussian Q and non-Gaussian ω distributions must be simulated with full use of the general rule, eq . The resulting 2D spectra are shown in Figure , where we observe a moderate curvature of the CL that is t2 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(t) = e–2αQg(0) in the ground state (eq ). The initial position Q(0) shall be obtained from inverting eq 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 , but the difference is usually small and we shall neglect it for the present analysis. The frequency evolution in ground state is thenand the FCL is represented by graphing ω(0) vs ωg(t). The BCL can be obtained by axial symmetry (eq ).

The joint probability for SE pathways (eq ) isThe center of this Gaussian distribution approaches the excited state potential minimum asThe evolution of transition frequency readsGSB and SE line shapes are usually similar, being only slightly shifted along ω3. We can thus conclude that CL of eq are curved ∝ αe – αg, which was confirmed by the full simulation shown in Figure .

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 (D4) molecule (see Figure 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 Ω ≈ 700 cm–1 and Ω ≈ 1600 cm–1 are readily identified as vibronic side bands.[61]

Figure 10

Bottom right: molecular structure of ZnPc. Top right: level structure of ZnPc. Three electronic levels S0, S1, S relevant for 2D ES. Bottom left: absorption spectrum of ZnPc dissolved in benzonitrile (solid blue line). The red line shows the laser pulse spectrum. Top left: 2D electronic spectrum of ZnPc dissolved in benzonitrile measured at room temperature and t2 = 96 fs. The DP1 peak corresponds to the S0 → S1 electronic transition. DP2 is the first vibrational peak. ESA denotes excited state absorption S1 → S. BCL and FCL are indicated as black and red lines, respectively.

The experimental details of the 2D ES measurements of ZnPc (Figure ) 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 , 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 × 1014 photons/cm2 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 t1, and 5.3 fs for t2. 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 t1 and t2 delays, spectra were recorded by integration over approximately 105 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 t2 = 96 fs in Figure . 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 .

We next analyze the shape and the center line of the principal peak DP1 of the ZnPc spectrum (Figure ) 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 . 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 . (iv) The effect of the Stokes shift on the CLS of DP1, as discussed in section , 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 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 prescribed as eq ; ESA contributions were included using eq and homogeneous dephasing[67] as the fast component of relaxation . 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 , 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 11

Comparison of experimental (top row) and simulated (bottom row) 2D ES of ZnPc at waiting times t2 = 24–96 fs. Simulations assume a three-level chromophore modulated by an overdamped solvent mode (eq ) and two underdamped vibrational modes (eq ). The color scale is the same as in Figure . Parameters for simulation were retrieved as follows: Ω = 700 cm–1, Ω = 1600 cm–1, ωeg = 15000 cm–1. Solvent parameters: λ = 80 cm–1, Λ = 5 cm–1. Vibrational parameters: λV, = 40 cm–1, λV, = 160 cm–1, γ = 10 cm–1, γ = 150 cm–1. Parameters for ESA pathway (eq ) ωef = 15300 cm–1, gff = gee. Homogenous dephasing Γ = 30 cm–1. Temperature T = 300 K.

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

Figure 12

Inverse of BCLS extracted from experimental (black line) and simulated (red line) 2D ES of ZnPc. The slope was measured at the central part of the DP1 peak. Parameters are the same as in Figure .

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, Ω, and Ω 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 . 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, = 40 cm–1, fitted form BCLS corresponds to the Huang–Rhys factors , i.e, λV, = 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 ) 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.