Four-dimensional coherent electronic Raman spectroscopy

The physical properties of molecular systems are strongly dependent on the intra- and inter-molecular interactions among the electronic, nuclear, and bath degrees of freedom. Revealing these interactions is challenging, in large part, because of the vastly different energy, and hence time scales, of the associated quantum-mechanical modes. Spectroscopic methods are often employed to tease out these interactions but interpretation of the signals is oftentimes hampered by congested spectral features. This, in turn, leads to uncertainty in the physical origin of signals free from potentially biased microscopic models. Spreading the spectral bands across multiple dimensions has proven to be an extremely successful approach for revealing the underlying physical structure of many complex systems. Perhaps, there is no better example of this than multi-dimensional nuclear magnetic resonance (nD-NMR) spectroscopy,¹ which systematically maps networks of coupled nuclear spins, thereby revealing molecular structure in solution. The optical analog, multi-dimensional coherent optical spectroscopy (nD-COS) aims to reveal the structure and dynamics of molecules, which may be thought of as complex networks of oscillators and dipoles. Specifically, nD-COS aims to provide quantitative information on intra- and inter-molecular coupling among electronic, vibrational, and bath degrees of freedom and the effects of these couplings on dynamics on an ultrafast time scale.^2,3 As with nD-NMR, nD-COS requires the use of multiple and short (relative to the dephasing time) electromagnetic pulses in a specific timed sequence that encodes the various interactions as multi-frequency correlations.

One of the most widespread nD-COS methods is 2D Fourier Transform (2DFT) spectroscopy, which may be performed in the UV/visible/near-infrared (referred to as 2D electronic spectroscopy) and mid-infrared (referred to as 2D-IR spectroscopy) regions of the spectrum. Two-dimensional electronic spectroscopy (2D-ES) reveals which electronic transitions share a common ground state. For instance, electronic coupling may be manifested as off-diagonal features in the 2D spectrum. Information gleaned from the lineshapes of the transitions reveals the interplay between homogeneous and inhomogeneous broadening, as well as interactions of the system with the fluctuating bath. As the 2D-ES spectrum is tracked as a function of the waiting time (analogous to a pump-probe delay), the appearance and disappearance of features indicate dynamical changes such as those arising from energy or charge transfer processes or those caused by rapid chemical transformations. In the infrared, the dynamics reflect processes that occur on the ground electronic state, thereby indicating both the static and dynamic couplings between different vibrational degrees of freedom. 2D-IR allows ultrafast chemical reactions to be tracked by measuring changes in the coupling between different infrared-active modes of the reactant and product. Irrespective of the spectral region probed, 2DFT methods are third-order in the applied electric field and therefore fall under the category of four-wave mixing (FWM) techniques.⁴ While multi-dimensional vibrational and electronic spectroscopies may be considered as complimentary, for systems that have strongly coupled quantum-mechanical degrees of freedom, there may no longer exist a simple separation of energy scales. While there have been attempts to measure electronic-vibrational interactions directly,⁵ the origin of the couplings are typically obscured by the lack of spectral resolution inherent in any third-order method.

Recently, our group has explored high-order⁶ and high-dimensionality⁷ approaches to directly interrogate couplings between electronic and vibrational states. Specifically, we have experimentally demonstrated a four-dimensional coherent spectroscopy, which we call GRadient-Assisted Multi-dimensional Electronic Raman Spectroscopy (GAMERS). The GAMERS pulse sequence is shown in Figure 1. A key and unique feature of GAMERS is the combined use of resonant and non-resonant excitation. In essence, the non-resonant excitation amounts to impulsive Raman scattering, which places the system in a non-equilibrium coherent state on the ground-state potential surface. This state, rather than the thermally populated ground state, acts as the initial state for subsequent resonant 2D ES. In total, the pulse sequence creates a pair of zero-quantum coherences (ZQCs) and a pair of single-quantum coherences (SQCs). ZQCs, as their name suggests, are coherences between states with the same number of electronic quanta, while SQCs are those with a difference of one electronic quantum. Correlations between all of these coherences simultaneously encode information on electronic-electronic (2D SQC), vibrational-vibrational (2D ZQC), and electronic-vibrational (2D SQC-ZQC) interactions. Formally, GAMERS is a six-wave mixing (6WM) process as five electric fields act to generate a time-dependent polarization in the system, which then acts as a source for the emitted sixth field. In this contribution, we develop the theory behind GAMERS and compare it to other recently (and not so recently) introduced methods such as 2D Resonance Raman (2D RR)⁸ and 2D non-resonant Raman (2D NRR)⁹ spectroscopies. We highlight the unique aspects of GAMERS to reveal couplings that have thus far remained hidden in lower-dimensionality and lower-order experiments. We end by briefly discussing the prospects of GAMERS for studying quantum effects in which the degrees of freedom strongly couple. Such effects are thought to be responsible for the remarkable efficiency of photosynthetic pigment-protein complexes, perovskite solar cells, and singlet fission materials to name a few.

In GAMERS, the first two field-matter interactions originate from a time-coincident pulse pair (i.e., simultaneous Raman pump and Stokes), which creates a ground-state coherence that subsequently evolves on the electronic ground state for some time $τ1$⁠. This impulsive Raman scattering process is followed by three electronically resonant pulses, with inter-pulse delays $τ2$ and $τ3$⁠, which cause transitions between the ground state, S₀, and the first electronically excited state, S₁. The signal is emitted a time, $τ4$⁠, after pulse 3. Other higher lying excited states may be accessible within the pulse bandwidth (e.g., $S2→Sn$⁠). For the sake of simplicity, we will ignore these states which manifest in the spectrum as excited-state absorption (ESA). The signal is formally fifth-order in the electric field because the interactions caused by the first pulse pair depend on the electronic polarizability, while the subsequent resonant pulses interact through the dipole moment operator. Therefore, the typical perturbation theory does not apply. Further, in the experiment, the non-resonant pulse is set at a significantly higher intensity than the resonant pulse. Theoretically, we treat the interactions by first applying the Raman pulses to create an initial, non-equilibrium density operator, and then subsequently treat the resonant interactions by third-order perturbation theory.

This section is organized as follows: (A) We treat the non-resonant interaction by introducing the polarizability operator and discuss differences between resonant and non-resonant excitation on the density matrix. (B) We derive the optical response functions (ORFs), which are then combined to reproduce the signal emitted in a specific phase-matched direction. We discuss the specific pathways diagrammatically and their manifestation in the 4D spectrum. (C) Using the Franck-Condon approximation, we calculate the overlap integrals as a function of the displacement and relative curvature of excited and ground state potentials. (D) We discuss how the polarizability and anharmonicity conspire to create ground-state coherences between anharmonic wave functions. Using perturbation theory, we derive some analytical results which are useful for understanding the results of the simulations. (E) We discuss in which cases excited-state coherences may be distinguished from ground-state coherences.

To treat the non-resonant Raman interactions, we only need to consider the ground-state potential and the interaction of the system with the electric field.¹⁰ While this is an approximation, it is accurate in the limit of large detuning, especially when the frequency of the external field, $ω0$⁠, is far below the resonance of the $S0→S1$ transition, $ωge$⁠. If the Raman pulses are sufficiently below resonance, we may ignore excited state populations and coherences. This is valid as long as the off-resonant detuning, $Δω≡ωeg−ω0$⁠, is sufficiently large so that the excited-state population times are short, $τelec≈Δω−1$⁠. In practice, the detuning is such that $τelec<5fs$⁠. Since coherences require populations, the coherences are also limited to this time scale. We then obtain the equation of motion for the ground-state terms only

where $Agg≡PgAPg$⁠, P_g is a projection operator, and the tilde corresponds to the rotating frame defined by the time-independent Hamiltonian, H_g, which may contain both harmonic and anharmonic terms (see  Appendix A).

We may then solve this equation to first-order in perturbation theory to obtain

where $Aμμ′≡⟨μ|Agg|μ′⟩$ and $Δρμμ′(eq)≡ρμμ(eq)−ρμ′μ′(eq)$ is the difference in the equilibrium population between ground vibrational states. $H′(t)$ is the time-dependent Hamiltonian describing the field-matter interaction, which we write in terms of the electronic polarizability. Note that we do not treat the Raman scattering process as two separate field-matter interactions; rather, we treat the simultaneous pump and Stokes fields as a single interaction, which we may be described as

where $α^(𝐪)$ is the electronic polarizability operator and $𝐪=(q1,q2,…,q3N−6)$ describes the vibrational normal modes. Substituting this interaction into the density matrix expression, we arrive at

where

The polarizability operator may then be expanded in terms of the vibrational coordinates

where $α^(n)≡∂(n)α^(𝐪)/∂q1…∂qn$⁠. The polarizability operator may be written in terms of annihilation and destruction operators (see  Appendix A).

If we examine the time-dependence in Equation (5), we see that E²(t) represents the intensity envelope of the Raman field. The integral does not vanish as long as $δtpulse<ωμμ′−1$⁠. That is, the pulse duration should be shorter than the period of the vibrational coherence that is being excited, while remaining below resonance with the $S0→S1$ transition within the bandwidth of the transform limited pulse.

It is prudent to compare the density matrix immediately following non-resonant and resonant Raman (near-impulsive) excitation. If the excitation is resonant, we may perform a similar analysis to that above. However, because the interaction in the dipole approximation is given by $H′(t)$ = −|d|E(t), the first order term in the expansion vanishes as the dipole operator does not allow for direct transitions between vibrational states within the same electronic state when the frequency of the applied field is electronically resonant, $ω0≈ωge≫ωμμ′$⁠. An infrared field resonant with vibrational transitions could be used instead of the Raman field, which opens up exciting opportunities to study the relationship between IR and Raman active modes, but this will not be considered further here. For the resonant case, the first non-vanishing term in the expansion is second-order in the field, which physically corresponds to the pump and Stokes field-matter interactions in resonance Raman scattering. Further, we will assume that the electronic state has negligible lifetime, which implies that population and vibrational coherence on the excited state may be ignored. Under most circumstances, this is not the case and these processes greatly complicate the interpretation of the spectra in fully resonant experiments; we will discuss these issues later on. For now, we are mainly interested in a qualitative comparison of ground-state coherences created by resonant and non-resonant excitation.

In the case of resonant excitation, invoking the Franck-Condon approximation, and with the conditions prescribed above, the resonant analog of Equation (4) to lowest non-vanishing order is given by

where

We use Greek and Roman letters to indicate ground vibrational and excited vibrational states, respectively. Rigorously, the matrix elements of the polarizability operator, $α^$⁠, are calculated by summing over all the excited states in second-order perturbation theory, so the result for the resonant and non-resonant cases is perfectly consistent. We could have started with the well-known Albrecht theory for resonant and non-resonant Raman scattering,¹¹ but the approach given here is more tractable for incorporating into high-order response theory by minimizing the computational cost. The equation above also ignores higher-order terms in the interaction, which arise when the dipole operator has a coordinate dependence.

It can be seen by Equation (7) that in the resonant case, a non-vanishing product of Franck-Condon factors is required to place the system into a ground-state coherence, while the summation occurs over all excited states, ${|a⟩}$⁠. This just represents the fact that ground-state coherences may be created through any electronic state where there is sufficient wave function overlap (and within the bandwidth of the excitation pulse). The term, g(t), reflects the effect of resonant excitation. It is non-vanishing when the field has frequency components resonant with transitions from the ground to excited state, i.e., $μ↔a$ or $μ′↔a$⁠. If we now compare Equations (4) and (7) , we see that both non-resonant and resonant excitations create ground-state coherences, but the former is dependent on the off-diagonal elements of the polarizability operator while the latter is dependent on the product of Franck-Condon overlap integrals (summed over all excited states). Physically speaking, the polarizability operator depends on the derivatives of the transition dipole moment with respect to normal modes or combinations of normal modes. It is interesting to note that the polarizability tensor is sensitive to non-trivial vibronic coupling, such as those that occur when the Born-Oppenheimer approximation breaks down. Sensitivity of GAMERS to such effects will be explored in future work.

We treat the response functions to third-order in perturbation theory with respect to the dipole interaction induced by the resonant pulses. The “initial” density matrix, $ρ(0)=limt→0ρ(t)$⁠, is given in Equation (4), which is the result of the non-resonant Raman pulse pair acting on the thermalized density matrix. The response is then dependent on the four time delays $τ→=(τ1,τ2,τ3,τ4)$ in the pulse sequence according to

where the superscript emphasizes that the response function is third-order with respect to the resonant field and second-order with respect to the non-resonant field. V and $V$ represent the dipole operator and $G$ is the Liouville space Green’s function. $⟨⟨⟩⟩$ indicates the trace over all states. In terms of commutators, we have

where the optical response functions (ORFs) R_i and $Ri*$ are formed by combinations of commutators acting on the right or the left of the density matrix.⁴ These terms are written out explicitly in  Appendix A. Note that we have ignored excited-state absorption (ESA) pathways, which act to create excited-state coherences but with an overall negative contribution to the signal. These pathways, while critical for a quantitative analysis of the spectra, do not qualitatively change the general features for a monomeric system. Their analysis is beyond the scope of this work. However, ESA pathways are critical to understanding the underlying physics in a dimer or multi-chromophoric system where doubly excited states encode information on electronic coupling. This will be the focus of future work.

In the rotating wave approximation, not all terms survive.¹² We focus primarily on the signal that arises in a particular phase matched direction given by

where $Δ𝐤0=𝐤P−𝐤S$ is the difference in wavevectors between the non-resonant Raman pump and Stokes pulses. The detected signal represents the photon-echo (PE) or rephasing pathway along the electronic frequency dimensions (⁠$ω4$⁠). We note that in the experiment, we observe the signal in the approximate photon echo direction $𝐤PE4(6)WM=−𝐤1+𝐤2+𝐤3$ because the term $±Δ𝐤0$ approximately cancels when the Stokes and pump pulses are collinear (note that for large bandwidths the signal direction may be significantly affected). Alternatively, these pathways may be separately detected if $𝐤P≠𝐤S$⁠. The advantage of a non-collinear geometry is that the signal is emitted in a background-free direction, enabling spatial filtering from the much stronger 4WM signals. In the phase-matched direction in Equation (11), the only terms that survive the rotating wave approximation (i.e., those terms that involve near-resonance of transitions within the bandwidth of the excitation light) are

The signal pathways may be represented diagrammatically as shown in Figure 2. For each pathway, the diagram depicts the fate of the term $|α⟩⟨α|$ in the equilibrium density operator. At $T=0K$ this is just the ground state of the system. If we examine the R₃ pathway, for instance, we see that the Raman interaction acts on the ket states (solid lines), creating a ZQC $|α⟩⟨γ|$ represented as a green circle. The system acquires phase during $τ1$ as it evolves on the ground-state potential. The phase is negative for this particular pathway because $Eα<Eγ$ at $0K$⁠, and this is represented as an open circle. The next pulse operates on the bra state (dashed lines), forming the SQC coherence $|a⟩⟨γ|$⁠, which evolves positive phase for $τ2$⁠. The next interaction operates on the ket and creates another ZQC, $|a⟩⟨b|$⁠, which is now a vibrational coherence on the excited state. The phase evolved during $τ3$ may be either positive or negative depending on the relative energies of a and b. The next interaction operates on the bra creating another SQC, $|β⟩⟨b|$⁠, that evolves for $τ4$ at which point the system emits back to a ground-state population, $|β⟩⟨β|$⁠. Notice that in pathways R₃ and R₄, the frequency of the signal along $ω1$ is typically negative (at higher temperatures, it could be positive), and that for all pathways, $ω2$ is positive, while $ω4$ is negative using the convention here. These sign differences, in some instances, offer clues to the origin of signals in the spectrum.

In the impulsive limit (pulse duration goes to zero), the Fourier transformation along all four time dimensions yields the four-dimensional spectrum

where $ω→≡(ω1,ω2,ω3,ω4)$⁠. Note that the forms of the response functions given in  Appendix A ignore population relaxation and dephasing, but these may be added phenomenologically by introducing vibrational and electronic lineshape functions. More sophisticated methods abound for introducing a variety of spectral density functions to more accurately simulate relaxation,¹³ but these are beyond the scope of this work. Finally, we should mention that, in principle, other signal pathways may be detected which give access to other ORFs. In addition, judicious combinations of signals corresponding to different pathways may help eliminate dispersive features seen in the spectra,¹⁴ thereby improving spectral resolution. We do not explore other pathways here but note that they may be of importance for specific applications.

The matrix elements in the ORFs may be written as

To evaluate the first term in the sum, $⟨a|α(0)⟩$⁠, we assume a displaced harmonic potential for the electronically excited state so that $|a⟩=|a(0)⟩$⁠. The harmonic parts of the ground-state and electronically excited-state are given by

where d_i is the displacement operator, $ωi,g$ is the energy of mode i in the ground state, $ωj,e$ is the energy of mode j in the excited state, and $Ni,e(g)≡ai,e(g)†ai,e(g)$ is the number operator. We further define

We may apply analytical expressions for Franck-Condon factors in terms of $βi$⁠, which is related to the ratio of curvatures of the ground- and excited-state potentials, and the Huang-Rhys factor $Si≡Di2/2$⁠. These expressions are found elsewhere.¹⁵ It should be noted that, in general, the normal coordinates of the ground and excited state may not be related by a simple displacement. Rather,

where $J^$ is a rotation matrix and $D→$ is the displacement vector. If $J=1$⁠, then

where the product is over the normal modes. However, if $J≠1$ then this separability no longer holds and mode mixing occurs (Duschinsky effect¹⁶). The Franck-Condon factors may still be calculated but become much more complicated and computationally expensive. Here, we ignore mode mixing, but note that for polyatomic systems, it may be necessary for quantitative analysis. The second term in Equation (14) is evaluated by solving for the anharmonic wave function (see  Appendix A).

To analyze the origin of different features in the 4D spectra, we need to carefully consider the effects of the polarizability and anharmonicity on the signal. In this section, we focus on developing intuition for understanding the effects of the Raman pulse pair on the creation of ground-state vibrational coherences. The Raman interaction on the ket states creates a coherence of the form $|γ⟩⟨α|$⁠, where the $|α⟩$ state is populated at thermal equilibrium. For the sake of discussion and in order to build physical intuition, we take $|α⟩≡|01(0),02(0),…,0m(0)⟩$⁠, where m is the number of vibrational modes and the superscript (0) indicates that the modes are harmonic. Of course, in the presence of anharmonicity, this state is not the ground state. However, the energy shift to second-order in perturbation theory for small diagonal anharmonicity is negligible relative to the resolution of the measurement. Therefore, in the limit of $T=0 K$⁠, we may take the ground state to be $|0→(0)⟩$⁠. In the simulation, we numerically diagonalize and use the partition function, so these approximations are not relevant. However, for the purposes here, these approximations allow some insight into the most salient features in the spectra. As discussed in Section I, the electronic polarizability operator creates a ground-state coherence. To second-order, the matrix element is given by

where we use the shorthand notation $|αi(0)βj(0)⟩≡|01(0),…,αi(0),…,βj(0),…,0m(0)⟩$⁠. We see immediately that if the states $γ$ are harmonic, then the Raman interaction creates coherences between the ground vibrational state and the states $|1i(0)⟩$ or the combination states $|1i(0)1j(0)⟩$⁠, with the former due to the first-order polarizability and the latter due to the second-order polarizability. Therefore, in the limit of no anharmonicity, peaks along $ω1$ at combination frequencies $ωi,g+ωj,g$ are indicative of second-order polarizability effects.

In the presence of anharmonicity, $|γ⟩$ is a super-position of the harmonic states, ${|α(0)⟩}$⁠. Assuming that the anharmonicity is small compared to the vibrational quanta, $g(3)≪ωi,g$⁠, we may apply perturbation theory to find the eigenfunctions of $Hυib$⁠,

where

To first order in perturbation theory

Therefore,

and

We may then evaluate the matrix elements of the anharmonic operator

Similarly,

Note that, if say j = k, then $|1j(0)1k(0)⟩≡|2j(0)⟩$⁠, etc.

Thus, the first-order polarizability creates coherences between the ground state and anharmonic states $|1i⟩$⁠, $|1j1k⟩$⁠, and $|2i1j1k⟩$⁠, while the second-order polarizability creates coherences between the ground state and the anharmonic states $|1k⟩$⁠, $|2i1k⟩$⁠, and $|2i2j1k⟩$⁠. The off-diagonal density matrix elements corresponding to these coherences are dependent on the anharmonic coupling, polarizability tensor elements, energies, and relative populations. For instance, the effect of diagonal anharmonicity in the presence of first-order polarizability is to create coherence between the ground-state and the two-quanta anharmonic state

If off-diagonal anharmonicity is present, then the anharmonic combination states may also be excited

The energies of these anharmonic states may be calculated by second-order perturbation theory. For instance, the anharmonic combination state above is shifted in energy from the harmonic combination state by

According to Equation (28), this is given by

Therefore, peaks along $ω1$ that differ from the harmonic combination bands indicate the presence of first-order polarizability and anharmonicity. The peak shift is then a direct measure of the off-diagonal anharmonicity between those modes. In general, many vibrational quanta may be excited by a combination of the polarizability and anharmonicity operators depending on the strength of the applied fields and anharmonic couplings.

It is evident from our discussion here and in Section II A that ground-state coherence created by the non-resonant Raman interactions is generally different from those created by resonant excitation. Due to the Franck-Condon factors, there may be modes in the ground-state that cannot be easily accessed through the excited state (e.g., by means of resonance Raman scattering). For instance, the excited-state displacement for some modes may be sufficiently large that, within the bandwidth of the pulse, the overlap integral vanishes even in the presence of anharmonicity. On the other hand, the non-resonant excitation is less sensitive to the displacement operator. This means that GAMERS may allow for unique opportunities to access regions of the excited state potential that would normally be off limits.

Since the first time period $τ1$ selectively encodes ground-state coherences (GSC), excited-state coherences (ESC) may only occur during $τ3$⁠. Specifically, these coherences occur in the R₃ and $R3*$ response functions, in which $|a⟩⟨b|$ coherences are generated. We reiterate that while ESA pathways also create such coherences, we will not consider them here. Here, we explore how these types of coherences may be manifested in the 2D Raman spectra and under what circumstances may they be spectrally distinguished.

We first note that at $T=0K$⁠, the pathways R₃ and R₄ are in separate quadrants of the $ω1/ω3$ spectrum than the $R3*$ and $R4*$ pathways. This is because in the former, the Raman interaction creates an $|γ⟩⟨α|$ coherence while the latter creates an $|α⟩⟨γ|$ coherence, a consequence of the interaction occurring on either the bra or ket states. We note that the experiment may also be performed in a fully phase-matched geometry that incorporates the Raman pulses to select either pathway separately. Within a given quadrant, the Franck-Condon prefactor for each pathway at a particular position in the 2D SQC spectrum is identical so that the amplitude of peaks is of no help in distinguishing ESC from GSC. However, these coherences may be distinguished if there is resolvable energy differences between $ωab$ and $ωαβ$⁠. These ZQCs may differ in energy under three possible scenarios: (1) the ground-state anharmonicity is sufficiently strong that the anharmonic energies shift away from the harmonic energies in the excited state (assuming no anharmonicity is present there), (2) the curvatures of the ground- and excited-state surfaces are sufficiently different (i.e., $βi≠1$⁠), or (3) the electronically excited states is not simply a displaced version of the ground-state (i.e., it involves a change in the normal mode coordinates so that $J≠1$⁠). We will only focus on the second possibility since we are in the limit of weak anharmonicity and no coordinate rotation as previously discussed. These will be interesting regimes, however, to explore in future work. The curvature differences in the excited-state are easier to analyze and measure because they lead to unique signature such as asymmetry between the absorption and emission spectra.¹⁷ 

The relationship between the ground- and excited-state energies is given by Equation (17). Therefore, we expect that the ZQC in period $τ3$ will differ by a factor of $Bi2$ if we assume the same effective mass for a given mode for each surface. In the limit that $βi= 1$⁠, diagonal features in the 2D Raman-like spectrum, $S(ω1,ω3)$⁠, will include both GSCs and ESCs. However, when $βi≠ 1$⁠, the ESCs should shift along $ω3$ by $βi2ωg,i$⁠. The shift may be observed if $βi$ is sufficiently large relative to the vibrational dephasing time of the mode i. It is important to realize that such spectral discrimination is not possible with all resonant excitation, in which the shifts would occur along both $ω1$ and $ω3$⁠. The Raman excitation restricts the coherence created in the first period exclusively to the ground state so that no excited state vibronic coherence can exist. Therefore, the spectrum along $ω1$ serves as a reference, which provides a direct means to distinguish the GSCs and ESCs.

There are other scenarios that allow GSCs and ESCs to be distinguished even if the vibrational energies for a particular mode on the ground and excited potential surfaces are identical. In some cases, for instance, R₃ and R₄ may be distinguished by taking advantage of the selectivity along $ω2$⁠. Consider the lower right quadrant of the 2D ZQC spectrum (⁠$ω1< 0,ω3> 0$⁠). In the presence of sufficient displacement, a cross peak may form at some location $(ω1,ω3)=(−ωg(e),i,ωg(e),j)$⁠. The use of the notation g(e) indicates that we do not know at this point if the states involved are on the ground or excited electronic state. In R₄, the fourth interaction is on the bra, thus forming a GSC during $τ3$⁠, while in R₃, the forth interaction is on the ket, thus forming a ESC during the same period. However, for some points along $ω2$⁠, the only way to form a cross-peak in the ZQC spectrum is to excite to a state below $0′0′$⁠, the lowest energy state on the S₁ surface (see red “X” marks in R₃ pathway in Figure 8). This is in contrast to R₄, which forms such cross peaks readily. In these cases, the signal for the R₃ pathway will mostly vanish. Therefore, resonance may also be used to distinguish some pathways under special circumstances. This example will be elaborated on in more detail in Sec. III.

In this section, we attempt to highlight the sensitivity of GAMERS to non-trivial effects concerning the nature and couplings within and between the ground and the electronically excited states. We first explore different 2D projections of the 4D spectrum for a simple reference model where all potentials are harmonic and vibronic effects are weak. These projections encode electronic-electronic, vibrational-vibrational, and electronic-vibrational correlations. Different coherence pathways are explored and used to explain salient features in the spectra. Specifically, we show how electronic and vibrational resolution is effectively enhanced by taking various projections or slices through the relatively sparse 4D spectrum. Next, we explore vibronic effects in which the displacement between the ground and excited potentials leads to subtle changes in the Franck-Condon factors that contribute to the signal. While such effects may be seen even in linear spectroscopies, we highlight how GAMERS isolates these interactions free from other dominant signal pathways, thereby taking direct advantage of its large dimensionality. The highly nonlinear nature of GAMERS amplifies subtle changes in the potentials which we directly exploit. Next, we consider anharmonic effects on the GAMER spectra. Again, we highlight signals that unambiguously identify even weak anharmonic effects, which are typically obscured by lower-dimensionality methods. We end by discussing the role of electronic coherence and how GAMERS may distinguish coherent dynamics occurring on different potential surfaces. All simulation parameters are summarized in Table I of  Appendix B.

We begin by considering the simplest model consisting of two harmonic modes. The curvature of the excited state is set to be the same as the ground state (i.e., $β1(2)=1$⁠). We set the bandwidth of each pulse (non-resonant and resonant) at $1200cm−1$⁠, and evaluate the response functions up to n_q = 10 quanta for each mode. While the response functions are calculated in the impulsive regime, we only evaluate terms which fall within the bandwidth of our pulses. This allows the response functions to be evaluated in the frequency domain, which greatly speeds up the calculations. Note that if the Huang-Rhys factor, S, is identically zero and $β= 1$⁠, there is no signal within the bandwidth of the excitation due to selection rules. This can be seen by examining the pathways in Figure 2: the Raman interaction places the system into a ground-state coherence, restricting transitions to be of the type $νμ↔ν′μ′$ where the prime refers to the excited state. This necessarily leaves no pathways by which the system is left in a population after signal emission. Therefore, we set the Huang-Rhys factor to a small, but non-vanishing value, S₁₍₂₎ = 0.01, corresponding to a displacement along each mode of $D1(2)≈0.14$⁠. We use Lorentzian lineshape functions to convolve the ORFs, with $Γυib=10cm−1$ and $Γele=100cm−1$⁠. The same dephasing rates for ZQCs in the ground and excited states are used unless otherwise specified. In order to generate initial GSC by means of impulsive Raman scattering, we set the first-order electronic polarizability to be nonzero, $α1(2)(1)=−100cm−1$ (note the actual numerical value is somewhat arbitrary as it only acts to scale the total signal). While we may calculate the 4D spectrum at any temperature, we limit ourselves to $T=0K$ to simplify the analysis and limit the computation to only one initial state. We only calculate spectra for a single molecule, thereby ignoring the effects of inhomogeneous broadening so that no photon echo is observed. The large dimensionality of the calculated spectrum (⁠$Nω1(3)=128$ and $Nω2(4)=64$ or nearly 10⁸ elements in total) necessitates the approximations given within the constrains of our current computational resources. This model serves as our reference 4D spectrum.

For the sake of comparison, we show the 4WM spectrum, $I3D(ω2,ω3,ω4)$⁠, calculated using the ORFs given in Equations (A15) and (A16), but now using the equilibrium density matrix, $ρ(0)=ρeq$ (thus replacing the $γ$ index with $α$⁠). In Figure 3, we compare the 4WM spectrum at $ω3=0$ and 6WM spectrum at $(ω1,ω3)=0$⁠. We refer to these projections as “DC” spectra. The remaining dimensions (⁠$ω2$ and $ω4$⁠) represent SQCs between the ground and electronically excited states. In the case of 6WM, we call this type of projection a 2D SQC spectrum. All parameters are identical for the simulations. Note that both the 4WM and 6WM 2D spectra are dispersive, rather than purely absorptive, because only the rephasing component (i.e., the signal in the direction $ks≈−k1+k2+k3$⁠) is calculated. The 4WM spectrum shows an elongated positive feature along the diagonal (⁠$ω2=−ω4$⁠) with a strong negative lobe above. These features are commonly observed in third-order 2D spectra and result from the product of matrix elements with different signs. The 6WM mixing spectrum shows a similar pattern with a positive central lobe flanked by negative features on either side with additional structure clearly resolved. Other peaks become apparent, especially along $ω2$⁠, where the first SQC is formed, that lie close to $ω2=ωge±ωg,i$⁠, indicated by the dashed lines in the figure. These features are a result of the first-order polarizability operator, $α(1)$⁠, which allows for excitation into additional modes compared to those allowed by the resonant transitions. These excitations evolve as ZQCs during $τ1$ but allow additional opportunities for transitions further downstream in the pulse sequence.

At each point, $(ω2,ω4)$ in the 2D SQC spectrum, we may extract $I4D(τ1,τ3)$⁠, which may be described as a 2D “beating” map. ZQCs evolve along both of these temporal dimensions as shown in the last panel of Figure 3. We show this representation to highlight the fact that in the experiment, the signal is acquired by parametrically scanning $τ1$ and $τ3$⁠, yielding many thousands of 2D beating maps, one for each spectral position at fixed $ω2$ (i.e., also known as the coherence axis) and $ω4$ (i.e., also know as the rephasing axis). Double Fourier transformation along these two axes yields the 2D ZQC spectra shown in Figure 4 at locations marked with various icons (circle, diamond, star, square) and also marked in the 2D SQC DC spectrum of Figure 3. For instance, at the center of the 2D SQC spectrum (labeled as a diamond), the 2D ZQC spectrum exhibits peaks at the two modes, $ωg,1(2)$ along the diagonal (red line) in the upper-right (+, +) and lower-left (−, −) quadrants. There are also strong peaks along the DC axis $ω3=0$ in both the positive and negative $ω1$ quadrants. Notice that there are no peaks along the other DC axis (⁠$ω1$⁠) nor are there strong cross peaks between the modes at locations, $(ω1,ω3)=(±ωg(e),1(2),±ωg(e),2(1))$⁠. The spectrum shows inversion symmetry but no mirror symmetry about $ω1(3)=0$⁠. The fact that no peaks are observed along $ω1=0$ is a consequence of the polarizability tensor, which only creates signal for coherences, rather than populations in the ground state (see $Δρeq$ term in Equation (4)). In the experiment, signal along this axis may arise from nominal excited-state population that is inevitably formed when there is slight overlap between the resonant and non-resonant pulses (see Figure 1) or when imperfectly subtracted 4WM signal contaminates the measurement. At different points in the 2D SQC spectrum, the 2D ZQC spectrum changes dramatically. For instance, at the position labeled by a square, prominent peaks are observed in the lower-left quadrant, while weaker signals appear in the upper-right quadrant. No peaks appear in the other two quadrants ((−, +), (+, −)) in any of the marked positions.

We may explain all of these features by considering the optical response functions. For sake of brevity, we only consider a relatively small number of pathways to develop intuition for how the signals arise. Consider once again the diamond position (⁠$ω2(4)=ωge$⁠). There are two dominant signal types, peaks along the diagonal at $(ω1,ω3)=(±ωg(e),i,±ωg(e),i)$ and those along the $ω1$ axes, $(ω1,ω3)=(±ωg,i,0)$⁠. If we restrict ourselves to the upper-right quadrant, we can rationalize these features by considering the pathway, $R4*$⁠. The diagonal feature arises because of the non-negligible Franck-Condon factors between states of the same vibrational quantum number such as $⟨00|0′0′⟩$⁠. The amplitudes of the peaks along the $ω1$ axes are even larger, however, because they have the fewest changes in vibrational quanta. For the star peak, the diagonal features in the lower-left quadrant are dominant. Again, we may rationalize this by considering the coherence pathways, this time arising from R₃. The dominant signals are those that involve transitions between states of the same quantum number. This is because the displacements for each mode in the reference model are small so that the largest matrix elements are those of the type $⟨ν|ν′⟩$⁠. Note that weak amplitude peaks are observed in the upper-right quadrant at the square peak location. Strictly speaking there should be no pathways that allow signals in this region. However, because the electronic transitions are significantly broadened, the signals have contributions from other regions such as circle and diamond. These small amplitude peaks, especially when studying systems at higher temperatures, make unambiguous assignments challenging.

Another useful means to visualize the 4D data is to plot the 2D SQC spectrum for a given point in the 2D ZQC spectrum. Examples of these types of projections are shown in Figure 5. As a reference, we plot the 2D DC SQC spectrum in the background to visualize how the 2D SQC spectra move relative to it. As we select different points in the 2D ZQC spectrum (see plot at the center of the figure), we retrieve markedly different 2D SQC spectra. This selectivity arises because the 2D ZQC spectra may only be generated by pathways that involve specific vibrational states. For instance, at the point marked in orange, $(ω1,ω3)=(−ωg(e),2,−ωg(e),2)$⁠, in the 2D ZQC spectrum, we observe a 2D SQC spectrum shifted to lower frequency along $ω2$⁠, with features at $−ω4=ωge$ and $−ω4=ωge+ωg(e),2$⁠. This is readily explained by considering the coherence pathways, in this case R₃, shown to the left of the figure. By fixing $ω1$ and $ω3$ for this pathway, the primary transitions in the case of small displacement are those which cause the frequencies along $ω2$ to be one vibrational quanta lower than $ωge$⁠, while there are two possible transitions along $ω4$⁠. One of these involves a $0′0′↔01$ transition on the bra followed by a $0′1′↔01$ transition on the ket. The other involves a $0′0′↔ 00$ transition on the bra followed by a $0′1′↔ 00$ transition on the ket. Note that the overall product of these transition moments is the same and therefore the two signals in the 2D SQC spectrum have roughly the same magnitude. We may think of such 2D SQC spectra as “super-resolution” 2D electronic spectra in the sense that they are necessarily more resolved than both the 4WM and 6WM DC spectra. Such spectra may be of great utility to accurately determine the origin of features in multi-dimensional electronic spectroscopy. For instance, mode-specific relaxation processes may be directly measured in this way.

Finally, we may analyze the 4D spectrum by considering a projection along a single ZQC and SQC axis. This effectively generates a type of 2D electronic-vibrational spectrum. We will highlight the importance and usefulness of such projections as we discuss vibronic and anharmonic effects in Sec. III B.

One of our goals is to analyze the sensitivity of GAMERS to vibronic effects. First, we will consider the effects of displacement between the ground and excited-state potentials on the 4D spectrum. To do this, we will compare the reference model to a system in which a single mode has a moderate Huang-Rhys factor, (S₁, S₂) = (0.2, 0.01). Comparing the 4WM and 6WM 2D SQC DC spectra in Figure 6, we observe that both are quite similar, showing the characteristic negative lobes on one or both sides of the positive central lobe oriented along the diagonal (⁠$ω2=−ω4$⁠) direction. Subtle differences may be seen, for instance, along $ω2$⁠, which again is a consequence of the differences in excitation induced by the dipole and polarizability operators. Further, these spectra are nearly indistinguishable from the corresponding reference spectra shown earlier. More striking differences are seen in the 2D ZQC spectra at the location marked with a star in the 6WM 2D SQC DC spectrum. A new peak at $(ω1,ω3)=(−ωg,1,−2ωe,1)$ is observed. We may attribute this peak to the R₃ pathway. The R₄ pathway (not shown) is very weak in this case because it would require the transition $0′↔ 3$ during the fourth field-matter interaction (via pulse 2). Therefore, we see that this signal is in fact indicative of an ESC. More generally, signals at harmonics of the fundamental vibrational frequencies are indicators of vibronic effects. They arise because the larger displacement operator permits transitions to vibrational quanta that would otherwise be forbidden. In this case, the $1→ 2′$ transition becomes favorable for mode 1, but remain largely forbidden for mode 2. This illustrates how GAMERS allows a direct means by which we can estimate the displacement per mode. Lower-order spectroscopies would project this information, either along $ω1$ in the case of non-resonant Raman, in which case the peak is completely hidden, or along $ω2$ in the case of resonant Raman, in which case there is no direct means to know if the peak is at a fundamental or harmonic frequency. The absence of the harmonic along the diagonal of the ZQC spectrum makes the assignment definitive.

It should also be noted that, as with the reference model, strong amplitude cross peaks do not appear anywhere in the 2D ZQC spectra. The reason for this is easy to justify by analyzing pathways that may give rise to cross peaks. Consider, for instance, pathway $R4*$⁠. In order to see a cross peak at $(ω1,ω3)=(−ωg,1,−ωg,2)$⁠, transitions of the type $1→0′$ must occur for mode 2. Since its displacement is small, the signals corresponding to these pathways are small. The intensity of transitions of the type $01→10$ caused by two field-matter interactions is therefore sensitive to the displacements along both modes since they involve changes of both vibrational quanta. The main differences between the reference spectrum and the single-mode displacement spectrum are in the relative intensities of the peaks, but without a unique spectral signature such as that arising from the harmonic pathways discussed above, it becomes exceedingly difficult to make accurate spectral assignments absent a microscopic model.

While the signature in the ZQC spectrum discussed above is indicative of mode-specific displacement, in practice it may prove difficult to isolate this peak in the presence of a strong diagonal peak or if the harmonic frequency is too close to another fundamental mode to readily resolve. It is, therefore, worth exploring other spectral projections that may better isolate the signals of interest. Out of the four spectrally resolved dimensions, $ω1→ω4$⁠, we may choose one of the six 2D projections. Besides the 2D ZQC and 2D SQC spectra discussed thus far, we may also examine the 2D ZQC-SQC spectra. In this case, we may select a 2D SQC spectrum at a fixed $(ω1,ω2)$⁠, $(ω1,ω4)$⁠, $(ω2,ω3)$⁠, or $(ω3,ω4)$⁠. Considering all of these possibilities is outside the scope of this work, but it is instructive to examine a few select projections where unique features may be isolated. As an example, we compare the 2D ZQC-SQC spectra at a fixed $(ω2,ω3)$ value. This leaves the ZQC-SQC correlation between $ω1$ and $ω4$⁠. As shown in the middle panel of Figure 7, slices along these directions easily reveal differences in the spectra that may be obscured by the 2D ZQC or 2D SQC spectra. In this particular example, the pathway was chosen so as to take advantage of the high nonlinearity of the measurement. For the R₃ pathway shown in the figure, there are three $0↔1′$ transitions for mode 1. When the Huang-Rhys factor is small, these transitions are strongly suppressed and the signal in the highlighted region is weak. As the HRS is increased to more moderate values, the peak intensity rapidly rises. For this pathway, the signal is proportional to the product of Franck-Condon factors, $|⟨00|1′0′⟩|2|⟨10|0′0′⟩|$⁠. As can be seen in the leftmost plot, these transitions are much stronger in the case of the displaced model than for the reference model. Therefore, the signal is strongly amplified in this region. This represents just one example of how different projections or slices of the 4D GAMERS spectrum may reveal differences otherwise hidden in lower-order or lower dimensionality spectroscopies.

Finally, we may also consider the case when both modes have moderate displacements, (S₁, S₂) = (0.2, 0.2). Again, we note that the 4WM and 6WM spectra are nearly identical to one another and, in general, look quite similar to the corresponding reference spectra. In other words, electronic spectra exhibit poor sensitivity to vibronic effects under these conditions. However, the ZQC spectra are markedly different. For instance, by analyzing the R₄ pathway, we see that cross-peaks in the 2D ZQC spectrum should now appear as both transitions $10↔0′0′$ and $0′0′↔01$ have appreciable oscillator strength. That is, there is a pathway by means of two Franck-Condon transitions that connects 10 and 01. If we select the diamond location in the 2D SQC spectrum in Figure 8, we observe cross peaks between $ωg(e),1$ and $ωg(e),2$ in both right quadrants of the 2D ZQC spectrum. The cross-peaks are necessarily weaker than the diagonal or zero-frequency peaks but are clearly visible. Unlike the reference case, strong signals now appear in the lower right quadrant because of the availability of these new transitions. The region highlighted in green may arise, in principle, from either the R₃ or R₄ pathway shown in the same figure. However, the R₃ pathway may be ruled out since it would require a transition to a state below $0′0′$⁠, the lowest state on the excited-state potential surface. Indeed we see that very little cross-peak amplitude exists at this location. The fact that some signal does appear is again a reflection of the poor separation of pathways from the broadened electronic spectra (i.e., the diamond position in the SQC DC spectrum has contributions from other points with ZQC cross-peak amplitude). We may therefore conclude that R₄ is the dominant pathway, which is born out by the calculations. This pathway involves transition to combination states, such as 11, which can only occur if displacements lie along both coordinates. It is important to remember that in the experiment signals from all pathways lie directly on top of one another. Understanding how the signals arise from different ORFs is therefore critical to interpreting the spectral features. Since R₄ is the dominant pathway, this implies that this signal does in fact correspond to a ZQC on the ground electronic state. We will discuss differences between ZQCs on different potentials in Secs. III C and III D, but here we want to highlight the power of GAMERS to distinguish these type of coherences even when the vibrational energies on the ground and excited states are identical (i.e., $βi=1$⁠). From this analysis, we conclude that the strength of the cross peak is highly sensitive to displacements along both coordinates, q₁ and q₂.

There are other features in the spectrum worth noting such as the peak inside the red square at $(ω1,ω3)=(ω1,g,ω2,g−ω1,g)$⁠. Considering the R₄ pathway, we immediately see its origin: the fourth interaction on the bra (pulse 2 in the pulse sequence) could cause a transition to the state 01 instead of the combination band 11. Pulse 3 would then need to act on the ket to cause a transition to $0′0′$⁠. Therefore, this spectral signature also is sensitive to the displacement along both coordinates as it requires a transition of the type $1′0′↔01$⁠.

In general, cross peaks between any two harmonic modes indicate that both modes have sufficient vibronic coupling, a consequence of using resonance to generate the non-linear signal. From our discussion in Section II D, we should expect that spectral signatures along $ω1$ indicate other non-trivial behavior, such as the presence of anharmonicity or higher-order polarizability. We first consider the case of diagonal anharmonicity, which is captured by the g_iii term in Equation (A6). Before analyzing the spectra, it is prudent to try and predict where some of the spectral features that encode information on the anharmonicity may reside. First, we may examine the effects of diagonal anharmonicity on the anharmonic wave functions, in particular, on the expansion coefficients in the harmonic basis. In Figure 9, we plot these expansion coefficients for the first 15 anharmonic states (in order of increasing energy). For instance, we see that diagonal anharmonicity mixes primarily the 00 and 10 states to form the anharmonic state $|0⟩anh$⁠, while the state $|1⟩anh$ is a mix primarily of the harmonic states, 00, 10, and 20. Because the anharmonicity is only diagonal, it cannot mix harmonic states from different modes. This mixing causes changes in the Franck-Condon factors as shown in second panel in the figure. However, these changes are very subtle for the lower-lying vibrational states. States with higher vibrational quantum number show far larger changes, but these are typically not accessible within the bandwidth of the pulse so we will not consider them here. These results emphasize that the changes to the Franck-Condon factors are unlikely to produce significant changes in the transition strengths involving the lower-quanta vibrational states. In essence, this is the challenge of using resonant spectroscopy as a sensitive probe of anharmonicity. On the other hand, GAMERS utilizes both resonant and non-resonant interactions, so it is worthwhile to examine differences in the initial coherence created in the anharmonic and reference models. From our discussion in Section II D, we see that in the presence of first-order polarizability, the only initial states that may be formed when anharmonicity is absent are between the ground state and the harmonic states 01 and 10. However, in the presence of diagonal anharmonicity, we should also observe coherences between the ground state and the harmonics of the modes. In the last column of the figure, we show the matrix elements of $|Δρα^(1)|$ on a logarithmic scale. We see that the initial coherences between the ground state, $|0⟩anh$⁠, and other anharmonic states do in fact change in the presence of even small amounts of anharmonicity. For instance, we observe a coherence between the ground state and $|2⟩anh$⁠, which is not present in the reference model as predicted based on first-order perturbation theory. This term is about 1/3 as large as the coherence between the ground-state and fundamental states (e.g., $|1⟩$⁠). However, it is important to note that while such coherences may be created initially, the nonlinear response, and hence the magnitude of the signal at a particular location in the 4D spectrum, still relies on Franck-Condon overlap integrals between these states.

To assess the overall influence of diagonal anharmonicity, we need to look at the overall optical response. We choose $g111=30cm−1$ and leave all other parameters identical to the reference model. In Figure 10, we show 2D ZQC spectra at two specific locations in the 2D SQC spectra marked with a circle and diamond. The 2D spectrum in the top middle panel of the figure shows the 2D ZQC spectrum at $(ω2,−ω4)=(ωge−2ωg,1,ωge+ωg,1)$⁠. This location was chosen based on the R₃ pathway (not shown) in which the initial ground-state ZQC is created between $|0⟩$ and $|2⟩$⁠. This is then followed by an interaction on the bra that drives a transition between $|0⟩$ and $0′0′$⁠, followed by an interaction on the ket which causes a transition to $2′0′$⁠. After pulse 3 and emission, the system is left in a 10 population. Overall, we expect this to be one of the strongest “unique” pathways in the sense that involves the product of relatively strong transitions that are effectively absent in the reference model. Indeed, we see that small amplitude peaks at the predicted location, $(ω1,ω3)=(−2ωg,1,−2ωg,1)$ are indeed observed. In fact, two peaks side-by-side are seen, corresponding to the small splitting induced by anharmonic interaction. Additional peaks arise when selecting other locations in the SQC spectrum, such as $(ω2,−ω4)=(ωge+ωg,2,ωge)$⁠, where we observe another small amplitude peak at $(ω1,ω3)=(ωg,1,2ωg,1)$⁠. However, in both these cases, the new peaks are quite weak and could easily be swamped by the tails of the much larger diagonal features. Therefore, it is fruitful to try and find other projections where the spectral signatures are relatively free from unwanted background signal. If instead we choose to examine the 2D SQC spectrum at the same location $(ω1,ω3)=(ωg,1,2ωg,1)$⁠, we observe more obvious difference in the signals between the reference and anharmonic models. Positive and negative lobe signals appear near $(ω2,−ω4)=(ωge+2ωg,1,ωge−ωg,1)$ which may be attributed primarily to the $R4*$ pathway. Unlike the previous example where we exploit the ground-state coherences formed in the presence of anharmonicity, in this case the spectral features arise primarily from subtle changes in the Franck-Condon factors, which promote pathways that are largely forbidden in the reference model. That is, the anharmonicity mixes the harmonic states in order to borrow oscillator strength from relatively strong transitions. This nonlinearity in the signal makes GAMERS exquisitely sensitive to relatively small couplings.

Next, we turn our attention to off-diagonal anharmonicity. In this case, we expect the creation of coherences involving anharmonic combination states that evolve during $τ1$⁠. Despite the creation of these new coherences, the 2D DC SQC spectra shown in Figure 11 is nearly the same as the reference and diagonal anharmonicity models. Again, this points to the insensitivity of purely electronic spectroscopy to identify weak vibrational couplings. On the other hand, the 2D ZQC spectrum shows the appearance of new features in the presence of off-diagonal anharmonic coupling (⁠$giij=30cm−1$⁠). These peaks are highlighted by the red and blue boxes in the figure, which are centered at $ω3=−(ωg,1+ωg,2)$ and at $ω1=−ωg,1$ or $ω1=−(ωg,1+ωg,2)$⁠. If we consider the latter (red box), then we may identify the dominant pathway as R₄. We see that in this case, a coherence is created between the ground state and the anharmonic combination state, $|11⟩$⁠, due to the combined effects of the first-order polarizability and the off-diagonal anharmonicity as predicted in Equation (30). We may also examine the 2D SQC spectrum at this location, where we see that the peaks have shifted relative to the reference case (black contour lines in the background). Even more obvious is the appearance of new peaks well-isolated from the background in the SQC-ZQC spectrum shown in the bottom right hand corner of the figure. Here we see the appearance of peaks along $ω1=−(ωg,1+ωg,2)$⁠, with one band at $−ω4=ωge$ and the other at $−ω4=ωge+ωg,1+ωg,2$⁠. Again, these features may be readily explained by examining the various ORF pathways. This highlights the multiple ways by which we can isolate the signals of interest that encode information on mode-specific vibrational coupling.