Interference among Multiple Vibronic Modes in Two-Dimensional Electronic Spectroscopy

Vibronic coupling between electronic and vibrational states in molecules plays a critical role in most photo-induced phenomena. Many key details about a molecule’s vibronic coupling are hidden in linear spectroscopic measurements, and therefore nonlinear optical spectroscopy methods such as two-dimensional electronic spectroscopy (2D ES) have become more broadly adopted. A single vibrational mode of a molecule leads to a Franck–Condon progression of peaks in a 2D spectrum. Each peak oscillates as a function of the waiting time, and Fourier transformation can produce a spectral slice known as a ‘beating map’ at the oscillation frequency. The single vibrational mode produces a characteristic peak structure in the beating map. Studies of single modes have limited utility, however, because most molecules have numerous vibrational modes that couple to the electronic transition. Interactions or interference among the modes may lead to complicated peak patterns in each beating map. Here, we use lineshape-function theory to simulate 2D ES arising from a system having multiple vibrational modes. The simulations reveal that the peaks in each beating map are affected by all of the vibrational modes and therefore do not isolate a single mode, which was anticipated.

Two-dimensional electronic spectroscopy (2D ES) is the ultimate four-wave mixing spectroscopy method because it can provide maximum resolution of the underlying nonlinear optical signal pathways [20][21][22][23]. In the most prevalent implementation of 2D ES measurements, three femtosecond laser pulses focus on a sample, from which they produce a four-wave mixing signal, which is a laser beam emitted by the sample that is subsequently detected by a spectrometer [24]. Heterodyning with a weak local-oscillator reference pulse allows researchers to extract both amplitude and phase information from the intensity-level measurement via the spectral interferometry method [25]. 2D ES is analogous to many more common methods from 2D nuclear magnetic resonance [26], and indeed the two techniques have similar interpretations regarding diagonal and cross peaks. Furthermore, 2D ES is equivalent to spectrally resolved transient absorption spectroscopy with the additional benefit of resolution of excitation frequencies [27]. In contrast to spectrally resolved transient absorption spectroscopy, however, 2D ES does not suffer from a time-frequency bandwidth constraint because the excitation frequency resolution arises from varying the time interval between the two pump pulses and a subsequent Fourier transformation.
The study of vibronic coupling using 2D ES dates back to some of the earliest 2D ES method-development research efforts of the late 1990s and early 2000s [28][29][30][31]. Since that time, numerous research groups have continued to pursue theoretical and experimental studies [32][33][34][35][36][37][38][39][40]. The result of these research efforts is that 2D ES measurements can be conducted with ultrabroad bandwidths adequate to resolve the complete manifold of vibronic peaks, allowing 2D ES to begin to be used as a mechanistic tool for photochemical reactions [41].
2D ES signals can be described theoretically using a perturbative treatment of the nonlinear optical signal [20]. In this method, the electronic system has a linear coupling to a manifold of harmonic oscillators to represent the surrounding solvent. The solvent causes dynamic fluctuations of the system's electronic energy gap, and these fluctuations are captured mathematically by a lineshape function, g(t), which is inserted into each third-order response function, R (3) (τ 1 , τ 2 , τ 3 ), to model the nonlinear optical signal that is measured in the laboratory. The lineshape function modulates the electronic transition during the three time intervals-τ 1 , τ 2 , and τ 3 -between the laser pulse interactions. This general approach has been used to model a variety of four-wave mixing signals including transient absorption and two-pulse photon echo, among others [42]. Indeed many 2D ES studies have used this theoretical approach [32,[43][44][45][46]. Despite this large body of literature, the interaction among multiple vibronic modes in 2D ES measurements has been observed [7,12,[47][48][49], but not studied in detail. Therefore, in this contribution, we examine explicitly how multiple vibrational modes influence any given beating map constructed from a sequence of 2D spectra, in a context without the many possible nuclear modes and extensive coupling interactions that often obscure this relationship in laboratory measurements. The findings reveal that beating maps-although created by isolating the oscillations at a specific frequency during the waiting time interval-are indeed influenced by all of the vibrational modes of the system.

Theory: 2D ES Simulations
We aim to model the optical response of a molecular system due to excitation by light, which usually takes the form of laser pulses in the case of 2D ES. We therefore use a two-level electronic system as depicted in Figure 1 to describe the sample molecule, with the initial system in the ground state, |g . Franck-Condon excitation enables electronic transitions between |g and the displaced excited state, |e , which are separated by a (thermally-averaged) energy gap with a frequency of ω eg . Solving for the first-order optical response function produces the linear spectrum resulting from excitation of the system, which can be conveniently described by the lineshape of the spectrum [20]. The absorption and fluorescence spectra, respectively, are given by, where i = √ −1. Here, g(t) is the function for the line-broadening component of the lineshape, and λ is a coupling strength parameter, sometimes referred to as the reorganization energy [50]. In Equation (2), the factor 2λ produces a redshift of the fluorescence peak with respect to the absorption peak, also known as the Stokes shift [20].
To simulate a 2D electronic spectrum, we calculate the third-order nonlinear optical response as a function of the three time intervals, τ 1 , τ 2 , and τ 3 . From the second-order cumulant expansion, we obtain the conventional third-order response functions, (3) Figure 1. Two-level system of harmonic oscillators used to model a sample with electronic states |g and |e (green) in the presence of a solvent environment (blue). Excitation induces an electronic transition (red) between |g and |e , which are separated by a frequency gap of ω eg . Nuclear modes, or vibrational transitions (light green) with characteristic frequencies indicated by ω 0 , also occur. Vibronic coupling takes place as coupling between the vibrational transitions and the electronic transition.
Both the linear and nonlinear response functions include the lineshape function, g(t), which takes the general expression [20], M(t) represents the electronic frequency correlation function, an analytic expression for incorporation into the nonlinear response functions. The correlation functions represent the model system and incorporate the necessary components for vibronic coupling. The electronic transitions induced by optical excitation also involve transitions between the many vibrational (sub-)levels within each electronic level-vibronic coupling refers to this coupled relationship between electronic transitions and nuclear vibrational motion. We apply the Brownian oscillator model for the correlation function to incorporate the vibrational modes of both the molecular sample and its surrounding solvent. In this context, a nuclear vibration within a molecule acts as a coherent harmonic mode, which can be modeled as an underdamped Brownian oscillator [20,51]. We use the following correlation function for an underdamped oscillator with frequency ω 0 : Here, γ is a vibrational relaxation rate parameter. Relaxation yields the reduced vibrational frequency, ω 0 , of the original nuclear mode with frequency ω 0 : In practice, the solvent surrounding the sample molecule also contributes nuclear vibrations that influence the total spectroscopic response. Transfer of vibrational energy from the primary nuclear mode to the vibrations of the surrounding solvent promotes relaxation and dephasing of the nuclear oscillation. To incorporate these effects on the optical response, we model the solvent as a group of low-frequency, overdamped Brownian oscillators [20,50]: In this correlation function, the parameter γ s represents the vibrational relaxation due to coupling to the solvent. Therefore, evaluating Equation (4) with the correlation functions for the nuclear modes (5) and solvent oscillations (7) yields their corresponding lineshape functions [50]: where we have set β = 1/ (k B T) consistent with convention. These lineshape functions arise from a displaced harmonic oscillator model in which the frequencies of the ground and excited electronic states are identical and the Condon approximation is used. Furthermore, the model itself is widely used but has some physical flaws discussed by others [31,52] that make it a tractable theoretical model but prevent it from being used to make direct comparison to measured spectra. Lastly, the total lineshape function is the sum of the lineshapes for the discrete nuclear modes, indexed by n, and the solvent: The complete lineshape function g(t) is used to evaluate the desired spectrum for the system. For example, the total g(t) can be inserted into the third-order response functions from Equations (3) to compute the 2D spectrum. Evaluation over the τ 1 and τ 3 ranges followed by Fourier transformation produces the 2D spectrum across the ω 1 (excitation frequency) and ω 3 (detection frequency) axes, at a specific τ 2 time that has elapsed since excitation of the sample with the pump laser pulses. A subsequent Fourier transform over a sequence of 2D spectra for varying τ 2 values produces a 3D spectrum [53], and slicing through the 3D spectrum at a specific ω 2 value produces a so-called 'beating map' of the oscillation amplitude and phase at this frequency [7,34,54]. Peaks arising in the beating maps indicate electronic transitions that induce oscillations at the ω 2 frequency-therefore, selecting ω 2 = ω 0 means that, if peaks appear in the beating map, some will be spaced apart with a distance equivalent to ω 2 corresponding to the frequency gap of vibrational transitions occurring in addition to the electronic transition. Prior studies of beating maps have not specifically examined the peak patterns that result when n > 1 for the underdamped vibrational modes, likely because the anticipated result is that each beating map would isolate the contribution of the vibrations at that specified frequency. However, as seen in Figures 4 and 5 in Ref. [55] and in Figures 3 and S5 of Ref. [56], peaks do appear in beating maps at spacings other than the selected ω 2 value. In this work, we examine this effect in detail by comparing beating maps computed using single-mode and multiple-mode systems that share vibrational frequencies.

Linear Response: Spectral Density and Absorption and Fluorescence Spectra
We selected vibrational frequencies of 9 THz and 40 THz as nuclear modes, and examined the spectra produced by the individual modes and the combination of both modes using the theoretical methods described above. Importantly, we selected the mode frequencies to be neither factors nor multiples of each other to distinguish their spectral peaks: This avoids misidentification of any potential overtone peaks from one mode as peaks due to interference between two modes. We also implemented an energy gap between the ground and excited electronic states with frequency ω eg = 500 THz, used a temperature high enough to satisfy the high-temperature limit [29], and applied consistent solvent parameters λ s = 60 THz and γ s = 2 THz to maintain the extent of homogeneous broadening and relaxation across all simulations [57].
First, we computed the linear response functions and produced 1D absorption and fluorescence spectra from the lineshapes in Equations (1) and (2). The spectral density, C (ω), is the 1D total correlation function in the frequency domain [20]. As shown in Figure 2, the amplitude of the spectral density is centered at the frequency of the nuclear mode(s), ω 0 , selected in each simulation. The Ohmic solvent contribution is visible at frequencies less than 5 THz. Figure 2 also shows the absorption and fluorescence spectra obtained from the 1D lineshapes for each simulation, with the maximum absorption frequency and the maximum fluorescence frequency separated by a Stokes shift of 2λ ≡ 20 THz in all simulations, which we kept consistent through our choice of parameter values.

Nonlinear Response: 2D Electronic Spectra
Next, we present the 2D electronic spectra for the same single-mode and multi-mode systems described above, which are shown in Figure 3. As expected, at early τ 2 time points ( 1 ps) after simulated pulsed excitation, the 2D spectra show a single positive peak that is elongated along the diagonal. The elongated, ellipse-shaped peak results from the line-broadening included in g(t), and is a typical feature observed in experimental spectra at early times due to spectral diffusion [58]. By 1 ps, the diagonal line-broadening has diminished and the main peak has dispersed. As relaxation occurs, another peak composed of the stimulated emission signal splits from the main ground-state bleach peak, while moving from the maximum absorption frequency toward the maximum fluorescence frequency, which is redshifted from the maximum absorption by the Stokes shift. In all simulations, relaxation completes by 5 ps, as evidenced by the separation of the two peaks by the full Stokes shift with the ground-state bleach signal remaining at the maximum absorption frequency and the stimulated emission signal shifted to the maximum fluorescence frequency along the axis of detection (ω 3 ).

Vibronic Spectra and Beating Maps from 3D Electronic Spectra
Finally, we examine the spectroscopic signatures of vibronic coupling by assessing 1D vibronic spectra, and beating maps measured at the frequencies of the simulated nuclear modes. A vibronic spectrum is obtained by taking the sum of the 3D spectral data in the frequency domain along the ω 1 and ω 3 axes, which yields a plot of spectral amplitude versus ω 2 frequency. As shown in Figure 4, a vibronic spectrum reveals the frequencies at which the vibrational wavepacket oscillates in the electronic states during time period τ 2 . Peaks indicate the frequencies at which vibrational transitions occur as a result of the optical excitation, with the peak intensities being determined by the coupling strength, λ. As expected, in simulations of one nuclear mode at 9 THz or 40 THz, the resulting vibronic spectra are dominated by a single main peak at the mode frequency ω 0 (Figure 4a). Notably, both spectra simulated with one mode show a clear additional peak that is much smaller in amplitude and located at a frequency of 2ω 0 . This peak corresponds to the overtone of the nuclear mode. When the system is simulated with two nuclear modes, the vibronic spectrum contains both primary peaks at their mode frequencies, the individual overtones, and combination-band peaks that were not present in the spectra of the individual modes (Figure 4b). The combination-band peaks have amplitudes on the same scale as the overtones and appear at 49 THz, corresponding to the sum of the two mode frequencies, and 31 THz, corresponding to their difference. The appearance of the combination-band peaks indicates that the vibrational wavepacket dynamics occur along both vibrational coordinates simultaneously, and that signatures of this complicated oscillatory motion should also appear in the 2D beating maps. To further examine the origin of the combination-band peaks, we performed a series of simulations in which we held the parameters of the solvent and 40 THz mode constant but varied the coupling strength, λ, of the 9 THz mode. We quantified the relative intensity by extracting the intensity of the sum and difference frequency peaks in the 1D vibronic spectra and normalizing them to the intensity of the 40 THz peak. Based on a series of simulations in which λ for the 9 THz mode varied from 0.4× to 1.8× its original value, we observed that a greater coupling strength enhances the intensity of the combination-band peaks. As follows, we surmise that increasing any component of the total coupling strength parameter will increase the intensity of any spectral signatures of vibronic coupling in the beating maps as well.
Next, we investigate the vibronic coupling in more detail by examining two-dimensional beating maps, which are slices of the 3D spectrum at selected ω 2 frequencies where peaks arise from electronic transitions that involve vibrational transitions occurring at this ω 2 frequency. Simulations of one nuclear mode display spectroscopic signal in the beating map magnitude (Figure 5a) measured at the ω 2 frequency that matches the nuclear mode frequency, ω 0 . All the maps presented here show clear signatures of oscillations of the vibrational modes as a result of the optical excitation, indicating vibronic coupling between the electronic transition and the nuclear modes. Furthermore, the beating maps show unique spectroscopic signatures for the different individual modes. As expected, the peaks appear with spacings that correspond to their nuclear mode frequency: For the 40 THz mode, all the peaks are spaced 40 THz apart in the ω 3 direction, and along the ω 1 axis they are also 40 THz apart, while the two rows of peaks on either side of the frequency of maximum absorption are 40 THz away from this frequency, where a node is located. Similarly in the beating map of the 9 THz mode, a node is also found at approximately the frequency of maximum absorption along ω 1 , with the main peak located 9 THz above it. The beating maps from single-mode simulations reveal multiple broad peaks, which, comparing to the respective 2D spectra in Figure 3, indicate that the oscillations occur across the entirety of the peaks. Furthermore, the horizontal node at the excitation frequency corresponding to ω max abs indicates that peaks above and below this value oscillate out of phase.
In simulations of a system with two nuclear modes, the resulting beating maps (Figure 5b) also show clear evidence of vibronic coupling, yet are strikingly different from those of the same modes studied individually. Notably, two-mode beating maps seem to share some spectral features with the beating map of one individual nuclear mode, and some with the beating map of the other nuclear mode. The two-mode beating maps contain a combination of the peaks of their two individual modes, and this can be observed at ω 2 frequencies equal to either of the two nuclear mode frequencies: The beating map of two modes at 9 THz and 40 THz (Figure 5b) shows significant oscillatory signal when examined at both ω 2 = 9 THz and ω 2 = 40 THz. In addition, both of these two-mode beating maps are noticeably different from the corresponding one-mode beating maps (Figure 5a), 9 THz at ω 2 = 9 THz, and 40 THz at ω 2 = 40 THz. Furthermore, we find peak spacings that correspond to both nuclear mode frequencies in the 9 THz + 40 THz beating maps. In the beating map measured at ω 2 = 9 THz, the two peaks above the frequency of maximum absorption are 40 THz apart along the ω 1 axis, and the displacement of the three peaks along the ω 3 axis is 49 THz, corresponding to the combination of both modes. In the beating map at ω 2 = 40 THz, we observe peak spacings of 31 THz, corresponding to the difference between the two modes, and 18 THz along the ω 1 axis, and both 31 THz and 49 THz along the ω 3 axis. Clearly, as observed in the vibronic spectra as combination-band peaks, there is also evidence of the oscillations of both modes present in the beating maps, measured at either frequency. Overall, these results indicate that the spectral features of multiple nuclear modes contribute additional signal arising from interference to the 2D electronic spectra.
An additional observation is that the beating maps from the same two-mode simulation vary distinctly when observed at the different ω 2 frequencies of the two modes. Again demonstrated in Figure 5, the two-mode beating map with 9 THz and 40 THz modes, for example, is more similar to the one-mode beating map of the 9 THz mode when observed at ω 2 = 9 THz, but retains more of the spectral signatures of the individual 40 THz mode beating map when observed at ω 2 = 40 THz. This is logical based on the assumption that beating maps should isolate only the oscillations contributed by the specified ω 2 frequency. However, the variance of the beating maps with multiple modes from the beating map of a single mode at the same ω 2 frequency demonstrates that the oscillations corresponding to this frequency are not completely isolated-this indicates that interactions between different modes are also evident in beating maps. To understand the origin of these interactions, we can re-examine the expressions for g(t) and the response functions, from Equations (8) and (3). The response functions include many terms, but the terms most relevant to the interesting peaks in the beating maps are the oscillatory terms including cosine and sine functions [51], contributed from the expression for g vib,n (t) as cos(ω 0,n t) + sin(ω 0,n t). Because the response functions always involve the exponent of g(t), when there are multiple modes present, this yields terms of the form e cos(ω 0,a t) cos(ω 0,b t) for the cosine term, and the parallel expression for the sine term. Using standard trigonometric identities, these terms can be recast in terms of the sum and difference frequencies. For example, 2 cos(ω 0,a t) cos The presence of both the sum and difference of the two modes in the response functions leads directly to the generation of the combination bands in the vibronic spectra and to the interference and mixing between oscillatory signals observed in the beating maps of two nuclear modes. The occurrence of overtones and combination bands as a result of vibrational coherences between multiple fundamental modes is well-documented in vibrational spectroscopy [59,60]. Interactions between coupled vibrations have been demonstrated to produce combinations of the frequencies and new cross peaks in other nonlinear spectroscopies as well [61][62][63][64]. Overall, the various couplings and coherences can be revealed as additional signals and peak splittings in multidimensional spectroscopy techniques, such as the beating map measurements demonstrated here. The transformation of 2D spectra into beating maps provides more detailed information in multidimensional electronic spectroscopy, and reveals mixed oscillatory signals arising from all nuclear modes that are vibronically-coupled to electronic transitions.
Finally, it is worth envisioning extending these results using a more physically accurate model. In the lineshape-function model used here, the vibrational frequencies are the same in the ground and excited electronic states. In real molecules, however, the normal vibrational modes of the excited and ground states can be slightly different because the bonds have different force constants. This is known to affect the vibronic structure that is observable in linear absorption and fluorescence spectra [65,66]. Due to completeness, one can write a transformation matrix between the two normal-mode bases. Such a matrix is known as a Duschinksy rotation matrix [67]. This means, for example, that one could write a normal vibrational mode of the excited electronic state as a linear combination of some or all of the normal vibrational modes of the ground electronic state. This effect has been studied in the context of time-resolved four-wave mixing measurements previously [68], and it is intriguing to consider how this effect could be incorporated into the lineshape-function model and what effect it would have. The main expected difference is that each fundamental peak in the vibronic spectrum would become a doublet, one frequency for each electronic state. The relative intensities would be related to the coupling strengths. This increase in total number of peaks would then lead to additional interference terms. We anticipate that the Duschinksy effect could be incorporated into this model using an extended method developed by Tanimura and Mukamel [69] and recently used by our group to model Herzberg-Teller coupling in 2D ES beating maps [45].

Conclusions
In this article, we have simulated 2D electronic spectra and the beating maps resulting from a system with multiple nuclear modes. The interactions between nuclear modes in 2D ES have not been investigated in detail in previous work, although their spectral effects have been observed in laboratory measurements. We demonstrate that beating maps, showing oscillatory signals occurring at the frequency (ω 2 ) at which the map is measured, do not isolate purely the oscillations corresponding to the ω 2 frequency as might be expected. Instead, if there are multiple modes, interference between their oscillatory signals arises in the beating maps as combinations of the spectral signatures from the individual modes. This ability to resolve signals of vibrational coherence and vibronic coupling through 2D ES beating maps is relevant to most molecular systems, which usually have numerous vibrational modes, and especially to many-atom materials like crystalline solids, which will continue to be of interest for applications in semiconductors, quantum computers, and photovoltaics. Further extensions of this model to include additional electronic states or doubly excited electronic states would be relevant for studies of electronic energy-transfer or many-body interactions, respectively. In such cases, beating maps can be especially useful for deciphering the coupling interactions that are possible between multiple electronic and vibrational transitions. The spectral interferences observed in this work will be important for interpreting measurements of these systems, especially as multidimensional spectroscopy is used to study a wider range of samples and more complicated excited-state processes in the near future.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: 2D ES-Two-dimensional electronic spectroscopy