Two-Dimensional Impulsively Stimulated Resonant Raman Spectroscopy of Molecular Excited States

Giuseppe Fumero , Christoph Schnedermann , Giovanni Batignani , Torsten Wende, Matz Liebel , Giovanni Bassolino, Carino Ferrante, Shaul Mukamel , Philipp Kukura , and Tullio Scopigno 1,6,† Dipartimento di Fisica, Sapienza Universitá di Roma, Piazzale Aldo Moro 5, Roma, I-00185, Italy Dipartimento di Scienze di Base e Applicate per l’Ingegneria, Sapienza Universitá di Roma, Via Antonio Scarpa 14/16, Roma, I-00161, Italy Physical and Theoretical Chemistry Laboratory, South Parks Road, Oxford, OX1 3QZ, United Kingdom Cavendish Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom ICFO -Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, Barcelona, 08860 Castelldefels, Spain Istituto Italiano di Tecnologia, Center for Life Nano Science @Sapienza, Roma, I-00161, Italy Department of Chemistry and Physics and Astronomy, University of California, Irvine, California 92697, USA


I. INTRODUCTION
The investigation of light-induced processes is essential to the understanding of a variety of complex phenomena at the interface between physics, chemistry, and biology, in which excited-state dynamics cause the transient reconfiguration of atomic positions and electronic phases. Following the absorption of a photon, the behavior of the system is dictated by the interplay between vibrational and electronic degrees of freedom of the system and the environment [1][2][3][4]. The evolution typically occurs along a multidimensional landscape, which is effectively represented by vibrationally structured potential energy surfaces (PESs) and determines the competition between all the available radiative and nonradiative relaxation channels.
Ultrafast spectroscopy exploits tailored sequences of laser pulses to photoexcite and subsequently probe these channels, with the aim of unveiling the dynamics and the underlying vibronic structure. This goal requires the correct identification of the excited-state PESs involved in the photoinduced process as well as the mapping of their relative orientations and displacements. A number of different pulse schemes and strategies have been developed to meet these tasks, but the unambiguous identification of vibronic properties, such as quasiparticle couplings, mode-mixing and nonadiabatic effects, remains challenging [5,6]-in particular, on electronically excited states after the system has left the Franck-Condon (FC) regiondue to overlapping signal contributions arising from different physical processes. Linear vibrational techniques, such as infrared and spontaneous Raman spectroscopy, cannot monitor vibrational coherences on the excited states, whereas UV-visible absorption spectroscopy usually lacks the desired structural sensitivity. Considerable efforts have been aimed at the development of multidimensional spectroscopic techniques in order to separate the photoinduced response over additional spectral dimensions [7][8][9][10][11][12][13][14][15][16][17][18][19], which arise from the parameters tunable during the experiment, such as the time delays between multiple excitation pulses. The correlations of features on different dimensions offer a different perspective on the vibrational landscape, providing a connection between the structure of the system and its dynamics.
The development of a 2D-ISRS scheme aimed to probe electronically excited-state vibrations would disclose the conformation of excited energy landscapes, providing direct information on specific molecular properties, such as the geometrical configuration and orientations of different PESs. This information is encoded in the harmonic fifth-order Raman response, which vanishes in the off-resonant regime. Critically, resonant 2D-ISRS demands a substantial increase in complexity, regarding both the experimental layout [48] and the data interpretation, in order to establish a protocol able to disentangle the genuine excited-state contributions from ground-state features and to assign the measured 2D peaks to the corresponding molecular origin.
Here, we address these challenges by presenting resonant excited-state 2D-ISRS. In striking contrast with the time-resolved ISRS approach-which employs a photopump to create an electronically excited-state population and monitors the temporal evolution of vibrational frequencies during a photoreaction [30,32,49]-by exploiting three femtosecond pulses for stimulating Raman coherences, we induce and probe vibronic correlations on the electronically excited manifold. Building on the resonance Raman enhancement, in our realization we tune the optical wavelengths of the pulses used in 2D-ISRS in resonance with the static and transient electronic absorption transitions to isolate contributions pertaining to a targeted electronic state. Importantly, our pulse configuration additionally suppresses undesired lower-order cascade effects, which typically mask the 2D Raman features [50,51].
We decipher the complex multidimensional response by developing an analytical method to read out the properties of vibronic coupling, which allows us to map out the multidimensional PES from the intensities and locations of diagonal and off-diagonal peaks in the 2D-ISRS spectra. We subsequently apply our approach to study the wild-type green fluorescent protein (GFP) [52,53], a prototypical system used in the fluorescence bioimaging community, during its initial photoinduced relaxation process. By collecting 2D-ISRS spectra over the full vibrational fingerprint region, we project the initially induced coherences on a separate temporal dimension and then probe their correlations as they evolve out of the FC region. We rationalize our experimental results by means of a harmonic treatment of the molecular Hamiltonian, demonstrating that this simple model is able to trace the origin of the different couplings. We further illustrate how a careful design of the experimental conditions enables the observation of dark or weak vibrational modes as coupling peaks, PES displacements along normal coordinates, and signatures of harmonic mode mixing in the excited states beyond the approximation of linearly displaced potentials.

A. Three-color 2D-ISRS technique
The pulse sequence of three-color 2D-ISRS, along with an illustrative sketch of its working principle, is presented in Figs. 1(a) and 1(b), respectively. A femtosecond, frequency-tunable actinic pulse E a , resonant with the S 0 -S 1 absorption maximum, promotes the system into an excited electronic state S 1 , and it impulsively generates vibrational coherences due to its broad spectral bandwidth, provided that the pulse duration is shorter than the vibrational period [27]. Following a variable time delay T 1 , a femtosecond impulsive Raman pulse E p , resonant with the S 1 -S n excited-state absorption, induces additional vibrational coherences on S 1 . After a second variable time delay T 2 , a broadband white-light continuum (WLC) probe pulse records the temporal evolution of the vibrational coherences through spectrally-resolved transient absorption. The nonlinear polarization induced by the actinic and Raman pulses oscillates during T 1 and T 2 at the frequencies of the stimulated Raman-active modes, thereby modulating the transmitted broadband probe pulse. The spectral content of these coherent oscillations is retrieved via a 2D Fourier transformation along T 1 and T 2 , which reveals correlations between Raman modes appearing as cross and combination peaks. By tuning the actinic and Raman pulses into resonance with different excited-state transitions and selecting the probe wavelength region corresponding to the stimulated emission, the proposed scheme allows us to select S 1 as the PES on which the vibrational coherences are created and probed. We refer the reader to Fig. S2 of the Supplemental Material (SM) [54], where we report on the modification of the 2D-ISRS map upon changing the resonance condition of the probe pulse.
Specifically, we extract the oscillatory part of the detected signal at each probe wavelength along T 2 for every recorded point in time along T 1 , which yields the signal SðT 1 ; T 2 Þ [Figs. 1(c) and 1(d)]. Note that SðT 1 ; T 2 Þ (reported in Fig. S1 of the SM [54]) is subsequently 2D Fourier transformed and averaged over the spectral region of interest to provide the final 2D-ISRS map SðΩ 1 ; Ω 2 Þ [ Fig. 1(e); see the Appendix C for further details]. Considering that the 2D Fourier transform of a real time-domain signal is centrosymmetric, we present only the first and fourth spectral quadrants, which are associated with signal components of the same and the opposite sign along the two frequency axes, respectively [55].
Importantly, since the Raman pulse must interact with a molecule that was previously excited by the actinic pulse in order to contribute to the excited-state 2D-ISRS signal, the resonant scheme renders our technique free from lower-order background involving cascades between two different molecules that generally affect fifth-order spectroscopies [50,51,56,57].

B. Experimental measurements of wild-type GFP
To explore the capabilities of our technique, we measured resonant 2D-ISRS spectra of wild-type GFP during the early stages of its photodynamics. Under physiological conditions, the GFP chromophore exists predominantly in a neutral, protonated A 0 form [ Fig. 2(a), gray structure]. Photoexcitation at about 397 nm [ Fig. 2(b), gray absorption spectrum] promotes the chromophore to the first excited state A* [ Fig. 2(a), blue structure], which subsequently undergoes an excited-state proton transfer (ESPT) of the phenolic proton from Tyr 66 to Glu 222, across a hydrogen bonding network to yield the highly fluorescent, deprotonated I* form [ Fig. 2(a), green structure] [58].
We focus on the rapid initial (∼1 ps) relaxation out of the FC region, which precedes ESPT, occurring on a longer timescale (2-10 ps) [59,60]. As the photoinduced evolution  (c) Transient absorption signals recorded along T 2 for a given T 1 over all probe wavelengths λ s , and the vibrational coherence contribution, which is isolated by fitting and subtracting the electronic background. (d) At each λ s , we collect the coherence for each T 1 to build SðT 1 ; T 2 Þ, which is (e) subsequently 2D Fourier transformed to yield SðΩ 1 ; Ω 2 Þ, after averaging the WLC probe wavelengths λ s along the desired resonance region.
is accompanied by a large Stokes shift, the transient absorption spectrum after 1 ps from photoexcitation at 397 nm exhibits the characteristic stimulated emission band at 509 nm known for GFP, along with a broad photoinduced absorption band at higher wavelengths [720 nm, Fig. 2(b), green spectrum]. Beyond transient absorption, we perform 1D-and 2D-ISRS measurements on GFP using a 15-fs actinic pulse in combination with a 9-fs Raman pump pulse, whose spectral envelopes are shown in Fig. 2(b) (blue and orange spectra), by following the resonant strategy described previously and sampling delays of 1 ps along the T 1 and T 2 dimensions. In Fig. 2(c), we compare the 1D spectra obtained by two-pulse ISRS in the absence of the Raman pulse, probed immediately after the photoexcitation [ Fig. 2(c), blue spectrum] and after relaxation from the FC region [ Fig. 2(c), orange spectrum, averaged over 300 fs along T 1 ], which correspond to the A* excited-state Raman spectrum. The comparison shows that the "fingerprint" Raman band positions up to 1200 cm −1 are not affected by the FC relaxation on A*. This observation is in line with previous observations suggesting that the proton motion leading to the formation of the I* form is not dominant in the GFP subpicosecond transient dynamics [61].
As a consequence of the negligible changes in the chromophore configuration during the initial relaxation process, the corresponding 2D-ISRS spectrum, shown in Fig. 3(a), presents pronounced peaks along the principal diagonal at Fig. 3(a), black dashed diagonal line] that match frequencies obtained in the corresponding excited-state A* Raman spectrum [compare to Fig. 2(c), orange spectrum]. In addition, we observe several prominent off-diagonal peaks that occur only along vertical lines [i.e., parallel to Ω 1 , black dotted vertical lines in Fig. 3(a)], while horizontal correlations (parallel to Ω 2 ) appear to be missing. Off-diagonal peaks are gathered mainly in three regions: a vertical stripe for Ω 2 ¼ 1147 cm −1 and two subdiagonals for The spectral separation between the ground-state absorption, stimulated emission, and excited-state absorption guarantees the absence of ground-state contamination to the excited-state 2D maps. This result is further confirmed by the lack of contributions of ground-state Raman frequencies in the measured 2D spectra.

C. Origin of the 2D-ISRS couplings within a harmonic model
The interpretation of 2D spectra is generally hampered by similar spectral features that may arise from different and concurring physical processes [62]. In order to extract the structural information contained in these measurements, we derive the 2D-ISRS signal using a perturbative framework based on the density matrix expansion [63]. We then use the obtained analytical expression (see Eq. (D1) in the Appendix) to fit the experimental data, as shown in Fig. 3  Photoexcitation of the neutral A 0 chromophore at 397 nm prepares A*, which decays with biphasic dynamics (2-3 and 8-10 ps) to the fluorescently active deprotonated I* chromophore [59,60]. (b) Normalized absorption spectrum (gray) and transient absorption spectrum at 1 ps (green) and employed actinic and Raman pulse spectra (blue and orange). (c) Comparison of FC spectrum (blue) and A* Fourier amplitude spectrum binned in T 1 to 300 fs to average out any effect of oscillatory modulations along T 1 (orange). The absolute amplitude of the FC spectrum was scaled by 0.2 for clarity. Raman spectra were averaged over a probe wavelength region from 575 to 615 nm (stimulated emission), and coherent oscillations were Fourier transformed over a time delay of 1 ps (see Appendix C for further details).
presented in the Appendix D, and in Secs. c and d of the SM [54]. Briefly, the physical observables can be related to the nth order nonlinear optical polarization P ðnÞ , which consists of the convolution between matter correlation functions and the electromagnetic fields [63]. The radiation-matter interaction is treated perturbatively, and the density matrix is expanded in power of the fields, applying many-body Green function techniques in Liouville space. Diagrammatic representations are exploited to isolate all the relevant terms in the expansion and calculate nonequilibrium expectation values of the correlation functions [64,65]. In particular, 2D-ISRS signals originating from the fifth-order nonlinear polarization P ð5Þ are represented by the Feynman diagrams shown in Fig. 4(a).
The 2D-ISRS signal depends parametrically on the frequencies ω i and lifetimes γ i of the modes and the dipole matrix elements μ ij , in a way determined by the choice of the model molecular Hamiltonian, H 0 . As H 0 is selected, these parameters can be quantitatively computed by fitting the experimental observations to the chosen model. To describe the experimental data, we assume the harmonic oscillator (HO) approximation for H 0 , in which the vibrational manifold associated with each electronic state can be depicted as an n-dimensional parabola. The dipole matrix elements are given by the FC overlap integrals between the initial and final wave functions, whose PESs are displaced along the normal mode coordinates with respect to each other. It is worth stressing that deviations from the harmonic regime manifest themselves mainly in offresonance conditions, when the harmonic contribution vanishes, allowing us to pinpoint the otherwise smaller contribution of single-mode and multimode anharmonicities and polarizability nonlinearities [42]. Conversely, under resonant conditions, the harmonic response is dominant. To connect the structural information with the 2D-ISRS maps, we first consider a linearly displaced harmonic oscillator (LDHO) model with two vibrational modes ω l and ω h and three electronic manifolds. Within this minimal scenario, two displacements per vibrational mode are required to determine the 2D-ISRS signal. The effect of each displacement d i along the mode ω i is to enhance or suppress a given component of the FC progression relative to that mode, gradually shifting the maximum from the 0 → 0 transition to higher overtones. Thus, the signal depends on the four displacements between the electronic surfaces along each vibrational coordinate d h;l 1 and d h;l 2 . The displacement d 1 determines the probability of a transition to a specific vibrational state after the interaction with the actinic or the probe pulse since the two laser fields are resonant with the same electronic transition S 0 -S 1 . We remark that the fifth-order signal is emitted as a result of a six-wave mixing process upon an additional light-matter interaction (usually called free induction decay), which is due to the relaxation of the nonequilibrium polarization , which lead to a diagonal peak at ω l and two combination bands at Ω 1 ¼ ω h AE ω l and Ω 2 ¼ ω l . Red (blue) horizontal and vertical lines highlight the ω l (ω h ) frequency. created by the external pulses. Despite the different nature of the two interactions involving the probe-pulse modes, shown as rainbow-colored arrows in Fig. 4(a), the dipole moment associated with the free induction decay is still governed by d 1 , which gives the probability for a relaxation to a specific vibrational level of the electronically ground state. Similarly, d 2 controls the probability of a transition to a determined state mediated by the Raman pulse, involving the different electronic resonance S n -S 1 . From the expressions of the FC integrals, we can determine the role of each parameter: d i 1 ≠ 0 implies FC activity of mode i, which can be excited by the actinic pulse and probed by the WLC as a diagonal peak in the map. In order to observe combination peaks, both d 1 and d 2 of the involved modes must be finite.
In real systems, excited-state dynamics often modify this simple scenario, even in the case of the LDHO picture. For example, in the presence of a dynamic shift of the vibrational frequencies due to coupling to a thermal bath, peaks are asymmetrically broadened with respect to the principal diagonal [44,66]. Here, we consider a different situation in which the dynamics are exhausted within the finite duration of the femtosecond pulses. In this case, the displacements recorded by the actinic pulse are different from those recorded by the probe pulse, as shown by Fig. 4(b), and three displacements per mode are required to appropriately describe the 2D-ISRS spectra. For these reasons, 2D-ISRS is able to also access dynamics that are too fast to be resolved as a shift in the Raman peaks or that do not modify the frequencies of the modes but only the transition dipole moments. A common situation is that upon photoexcitation by the actinic pulse, the system evolves to a different region of the S 1 PES, and the harmonic potential describing the PES near the new minimum is shifted from the initial value corresponding to S FC 1 to a new value corresponding to a potential S 0 1 , which we refer to as the dynamic LDHO. We thus described the subpicosecond dynamics of GFP through the dynamic shift of the S 1 potential induced by the actinic pulse, in line with the observed Stokes shift [see Fig. 2(b)], and we fitted the experimental data with the dynamic LDHO model [ Fig. 3(b)]. Notably, at odds with the simple picture for the static two-mode LDHO described above, off-diagonal peaks are possible even if d 1 or d 3 is vanishing for one of the modes. This case is shown in Fig. 4(c), in which we present the simulated 2D-ISRS signal of the dynamic LDHO depending on three displacements for the three typical contributions that may appear in the 2D maps: diagonal peaks (ω i , ω i ), cross peaks (ω i , ω j ), and combination peaks (ω i AE ω j , ω i ). The corresponding dominant Liouville pathways are reported in Sec. e of the SM [54]. Specifically, setting two of the six available parameters equal to zero, it is possible to understand the effect of each displacement. A vanishing d 2 along both the vibrational cordinates results in diagonal peaks at the two vibrational frequencies ω l and ω h [Fig. 4(c), left panel]. Cross peaks between the two modes arise from switching off d 1 for one mode and d 3 for the other. For example, the cross peaks at Ω 1 ¼ AEω l , Ω 2 ¼ ω h in the central panel of Fig. 4(c) can be obtained with d h 1 ¼ 0 and d l 3 ¼ 0. Finally, combination bands at Ω 1 ¼ ω h AE ω l ¼ ω AE and Ω 2 ¼ ω l as shown in the right panel of Fig. 4(c) are isolated by setting d h 3 and d l 2 equal to zero. Building on such considerations based on the dynamic LDHO model, we can summarize the conditions that generate the different contributions in the 2D-ISRS maps: (i) A diagonal peak at (ω i , ω i ) results from d i 1 , d i 3 ≠ 0. (ii) A cross peak at (AEω i , ω j ) indicates that d 2 is not vanishing, and it is comparable for the two modes. (iii) A combination peak at (ω i AE ω j , ω i ) implies a nonvanishing and comparable d 1 for the two modes. Similarly, a combination band at (ω i , ω i AE ω j ) implies a nonvanishing and comparable d 3 . Both additionally require d j 2 ≠ 0. In general, all these processes contribute to the total signal, with weights depending on the relative magnitude of the displacements. Since the signal emission requires the state of the density matrix after the free induction decay to be diagonal (see Sec. c of the SM [54] and Ref. [63]), an offdiagonal 2D-ISRS peak is obtained if the interaction with at least one of the pulses changes the total vibrational quantum number n by two so that Δn ¼ P i jΔn i j ¼ 2. For the cross and combination peaks analyzed above, this process is mediated by the Raman and the actinic pulses, respectively. Even if, in general, many vibrational modes are present, correlations between more than two modes are negligible in the limit of small displacements because they involve multiple higher-order processes with Δn > 1. Thus, the scheme we built for the two-mode LDHO can be applied to multimode systems by considering all possible pairs.
Finally, we note that, by exploring more than two different electronic surfaces, resonant 2D-ISRS is sensitive to the sign of displacements, a key advantage over lower-order 1D techniques only capable of resolving their magnitude. Indeed, if only one electronic transition is probed (for instance, S 0 → S 1 ), the signal is totally symmetric with respect to a change of the sign of d 1 . This result happens because, even if transitions with an odd difference of quantum numbers in the initial and final states scale linearly with the displacement, the displacement appears squared in the signal since two transitions between the same states are required to go back to a diagonal state of the density matrix. Probing an additional excited state breaks this symmetry. Even if no change occurs when all the excited-state potentials are displaced in the same direction, a different sign between pairs of d i and d j will modify the spectrum. Linear and also thirdorder techniques usually probe resonances between two electronic states, while resonant 2D-ISRS is able to specifically probe three PESs at the same time.

III. DISCUSSION
Building on the results of the LDHO model, we can now interpret the 2D-ISRS signal of GFP. The peaks shown in the experimental map and retrieved by the theoretical TWO-DIMENSIONAL IMPULSIVELY STIMULATED RESONANT … PHYS. REV. X 10, 011051 (2020) 011051-7 simulation (Fig. 3)  ≠ 0. Here, we discuss the implications of the predominant features in the 2D map, while we refer to Sec. f of the SM [54] for a table summarizing the origin of all the diagonal and off-diagonal bands obtained by the fit. The presence of off-diagonal peaks along vertical and diagonal rather than horizontal lines suggests that combination peaks at (ω i AE ω j , ω i ) are primarily observed and, consequently, that most of the modes are displaced along d 1 since combination bands require a nonvanishing d 1 for both of the modes involved in the signal generation. The 1010 − cm −1 mode shows high FC activity on the S n state, testified by the subdiagonal at Ω 1 ≈ Ω 2 − 1010 cm −1 (Fig. 3, red  arrow), which corresponds to a series of combination bands with other peaks due to its large d 2 . In contrast, the absence of any features at 1010 cm −1 on the principal diagonal shows that d 3 ¼ 0 for this mode. These observations illustrate the capability of 2D-ISRS of uncovering dark or weak bands by boosting small FC displacements via another transition.
We further observe coupling between the 1147-cm −1 phenolic C-H bend and a low-frequency mode at about 112 cm −1 (Fig. 3, green dashed arrow), which has recently been under discussion fueled by results from femtosecond stimulated Raman measurements on GFP in both the frequency and time domain [21,61]. While the functional importance of this coupling for the ESPT is under debate, the observation of oscillatory modulation of the excitedstate Raman spectrum is generally attributed to anharmonic vibrational coupling. Conversely, related studies on other molecules suggest that anharmonic couplings are challenging to isolate and interpret in fifth-order experiments due to competitive cascade processes [51]. Our study instead highlights that the coupling between the 1147and 112-cm −1 modes can also be explained within a harmonic model by the displacement of the excited-state PES along these vibrational coordinates, as supported by the agreement between the combination band in our 2D measurements and fit. We further note that the intensity of this combination band is stronger above the principal diagonal than below, suggesting additional contributions due to a cross peak between the 1147-cm −1 and the 1248-cm −1 mode. Within our framework, single coupling peaks in the 2D map do not necessarily imply a functional importance for the ESPT reaction coordinate, but the determination of the displacements and PES orientations retrieved from the 2D-ISRS features allows for identifying the coordinates more involved in the relaxation during the first picosecond after the photoexcitation. In the case of GFP, this timescale involves the relaxation of the molecular structure towards the optimal geometry to support the ESPT [61]. In particular, the values of the displacements for the mode at 1010 cm −1 point to a relaxation of the system upon photoexcitation along this normal coordinate, which is initially in an out-ofequilibrium configuration (d 1 ≠ 0) and fully relaxed at the end (d 3 ¼ 0).
Beyond the 112-cm −1 mode, we observe additional couplings of the C-H phenol mode at 1147 cm −1 (Fig. 3, vertical dotted line at 1147 cm −1 ). According to the model discussed so far, the series of peaks along the vertical at Ω 2 ¼ 1147 cm −1 indicates a large value of d 1 and d 3 for this mode, while the absence of a subdiagonal at Ω 1 ≈ Ω 2 − 1147 cm −1 points to a small value of d 2 . Critically, the subdiagonal Ω 1 ¼ Ω 2 − 2294 cm −1 (highlighted by the blue solid arrow in Fig. 3) is not captured by this model. Within the LDHO model, these features would indicate a strong d 2 displacement for the 1147 − cm −1 mode, such that both the fundamental and the overtone become FC active on the S n state. This hypothesis is, however, at odds with the absence of any peaks for Ω 1 ¼ Ω 2 − 1147 cm −1 in the measured map [ Fig. 3(a), blue dashed arrow].
A possible explanation of this discrepancy is the presence of a mode mixing between high-and low-frequency modes in the excited state, arising in the presence of two electronic states that have different equilibrium geometries with nonparallel corresponding normal modes, commonly denoted as Duschinsky rotation [67][68][69]. Under such circumstances, the normal coordinates on the excited state (Q 0 ) can be expressed as a function of the ground state (Q) as Q 0 ¼ JQ þ D, where J is the orthogonal Duschinsky matrix which depends on the rotation angles Θ D and D is the displacement vector. Even if the excited state is not displaced, the transition 0 → n with n ≠ 0 can be strong, if the mode is coupled by the Duschinsky rotation (J ≠ I) to a displaced one. In a similar way, the FC activity of a displaced mode can decrease due to the coupling to other modes.
In order to test this hypothesis, we evaluate the 2D-ISRS response of GFP incorporating in the model a mixing between the 1147-cm −1 and 80-cm −1 modes, with a fixed Θ D ¼ 90°Duschinsky angle. In particular, by fitting the experimental data considering the modes at 80, 822, 888, and 1147 cm −1 , which are those involved in the generation of the off-diagonal peaks at Ω 1 ¼ Ω 2 − 2294 cm −1 , we are able to reproduce the experimental features missing in the LDHO model, as shown in Figs. 5(a) and 5(b). It is now worth dissecting the implications of the Duschinsky mechanism on the relative intensities between the fundamental and overtone contributions to combination bands. To this aim, we isolate the effect of Duschinsky rotation in Fig. 5 by considering the minimal scenario of a single mode (1147 cm −1 ), coupled to a nondisplaced low-frequency mode (80 cm −1 ), detected on a "spectator" mode (822 cm −1 ). In the absence of any mixing, i.e., for a vanishing Duschinsky angle Θ D ¼ 0, 65% for all the values of Θ D in the range 65°-115°(see Fig. S4 of the SM [54]).
We note that the correct identification of the origin of all the off-diagonal peaks-and, in particular, the assignment of peaks involving the overtone of the 1147 cm −1 -is required to prevent wrong interpretations of the spectra; more importantly, it demonstrates the sensitivity of this technique to the Duschinsky mixing. Since the Duschinsky rotation affects the relative weights between vibronic transitions, mode mixing has a large impact on the kinetics and efficiency of charge transfer processes [70][71][72], making 2D-ISRS particularly suitable for studying this class of samples.

IV. CONCLUSIONS
We introduced a resonant 2D-ISRS technique for multidimensional Raman spectroscopy. By implementing a three-color experimental scheme, we demonstrated how to decipher the correlation between vibronic modes in electronically excited molecules with electronic-state selectivity, using a perturbative sum-over-states approach. Specifically, we showed that resonant 2D-ISRS is sensitive to such correlations also in the absence of anharmonicities and efficiently suppresses lower-order cascade signals, which affect other fifth-order Raman techniques. Moreover, using a WLC probe enabled us to further isolate the electronic states from which 2D-ISRS features originate, allowing us to study the vibrational manifolds on electronically excited states. Upon identifying the vibronic origin of each peak in 2D-ISRS maps, we elucidated how different mechanisms such as linear displacements along specific normal coordinates and mode mixing contribute to the experimental signal. Notably, while lower-order techniques can only record one-dimensional projections of potential energy surfaces, 2D-ISRS is able to efficiently map complex PESs, determining FC overlaps directly over multiple dimensions.
As a proof of concept for the experimental scheme and theoretical model, we applied the technique to study the subpicosecond FC relaxation in GFP. We revealed that high-frequency mode correlations can be sustained even by fully harmonic interactions and showed that the model is capable of reproducing the previously observed coupling of the 112-and 1147-cm −1 modes. In particular, we found large FC activity on the excited state of a mode that is dark in the stimulated emission transitions and an enhancement of the mode overtones, which we linked to the presence of Duschinsky mixing. In the presence of ultrafast dynamics on the electronically excited state, the selectivity of 2D-ISRS can be exploited to directly access the structural conformation on the state in which the dynamics originate, disclosing the initial stages of the reaction. We anticipate that comparing mode displacements between electronic states involved for the reactant and product may be relevant to identify the key reaction coordinates implied at different stages of the photodynamics.

APPENDIX A: SAMPLE PREPARATION
Plasmids containing wild-type GFP were transformed in an E. coli BL21(DE3) cell line [73]. Briefly, 0.5 μl of plasmid were added in 25 μl of competent BL21(DE3) cells and incubated on ice for 25 min, heat shocked for 45 s at 42°C, and incubated further for 2 min on ice, followed by the addition of 150 μl of prewarmed at 37°C super optimal broth with catabolite repression media and incubation at 37°C with shaking for one hour. In a Lysogeny broth (LB)media Petri dish supplemented with 100 μg=mL kanamycin, 50 μl of the cell culture were plated and incubated overnight at 37°C. A few colonies were selected and allowed to inoculate 100 μl of LB media in the presence of 100 μg=mL of the appropriate antibiotic and allowed to grow overnight.
The overnight cell culture was used to inoculate 8 L (1 × 8 L in 2-L flasks) of LB media with 100 μg=ml antibiotic. Cells were grown at 37°C to an OD 600 of 0.6-0.8 upon which 0.5 mM IPTG was added. GFP was expressed overnight at 20°C. Cells were harvested through centrifugation at 10 000 xg for 10 min, and the pellets were resuspended in cell lysis buffer (300 mM NaCl, 50 mM Tris, pH 7.4) supplemented with a protease inhibitor tablet (Roche). Cells were lysed by microfluidizer, and the cell debris were removed by centrifugation at 20 000 xg for 20 min. In the supernatant, imidazole was added to a final concentration of 20 mM, and the first purification step was achieved by immobilized metal ion affinity chromatography (IMAC) utilizing a 20-mL IMAC column equilibrated in buffer A (300 mM NaCl, 20 mM Tris, and 20 mM imidazole, pH 7.4). After four column volume washes with buffer A, the sample was eluded with buffer B (300 mM NaCl, 20 mM Tris, and 500 mM imidazole, pH 7.4). The sample was concentrated to a final volume of 10 mL and loaded in a size exclusion chromatography (SEC) column 26=60 Superdex 75 previously equilibrated in SEC buffer (150 mM NaCl and 20 mM Tris, pH 7.4). The selected fractions were pooled and concentrated to the final OD 397 ¼ 9.
During the spectroscopy experiments, the sample (∼10 ml) was continuously flowed through a 200-μm path-length flowcell at sufficient speed to ensure sample replenishment during consecutive shots. To minimize thermal degradation, the sample reservoir was ice cooled.

APPENDIX B: EXPERIMENTAL SETUP
All pulses were derived from a Yb:KGW amplifier laser system (Pharos, Lightconversion) providing 5-W, 180-fs pulses centered at 1030 nm at a 10-kHz repetition rate. The Raman pulse [800 nm, 9 fs, 160 nJ, sample beam diameter of 55-μm full width at half maximum (FWHM)] was generated by a noncollinear optical parametric amplifier (NOPA) pumped by the second harmonic of the laser (515 nm) [74]. The pulse energy was adjusted to keep twophoton absorption of ground-state molecules below 1%. The actinic pulse (400 nm, ∼15 fs, 150 nJ, sample beam diameter of 80-μm FWHM) was generated by frequency doubling the output of a second Raman NOPA (800 nm, 10 fs, 6 μJ) in a 25-μm beta-barium-borate crystal (θ ¼ 29°) placed near the sample. The fundamental was removed by two reflective harmonic separators (Eksma Optics). Probe pulses (sample beam diameter: 35-μm FWHM) were generated via WLC generation in a 3-mm sapphire crystal and detected in a home-built single-shot prism spectrograph [27]. The actinic and Raman pulses were modulated at 2.5 and 5.0 kHz by mechanical choppers, respectively, to remove lower-order contributions to the signal. T 1 and T 2 were sampled in 10.1-and 5.36-fs steps (Thorlabs-LNR50S/M and PhysikInstrumente-M-230.10, respectively). The exact step size for each translation stage was determined by reference ISRS measurements on toluene with the same pulse parameters. All pulses were vertically polarized.

APPENDIX C: DATA ANALYSIS
The data analysis is analogous to previously reported time-domain Raman studies involving an actinic and Raman pulse [30]. Briefly, for each time point along T 2 , we record three wavelength-resolved transient absorption maps as a function of T 2 , corresponding to the Ramanpulse-induced differential absorbance in the presence and absence of the actinic pulse as well as the actinic-pulseinduced differential absorbance in the absence of the Raman pulse. For all three transient absorption maps, we discard time delays prior to T 2 ¼ 100 fs, due to cross-phase modulation [75,76], and describe the electronic background by a global fit consisting of a sum of three exponential decay functions to isolate the underlying vibrational coherence. For each probe wavelength, we subsequently construct the two-dimensional signal SðT 1 ; T 2 Þ and apply a 2D mono-exponentially modified Gaussian window function [23-fs rise time (sigma), 1526-fs decay constant in T 2 , and 6-fs rise time (sigma), 422-fs decay constant in T 1 ] prior to zero padding and 2D fast Fourier transformation. To selectively probe only S 1 vibrational coherences, we only included data greater than 536 fs in T 2 , thereby suppressing additional signals due to higher-lying excited electronic states. Importantly, both time delays covered a total of 1 ps to avoid possible distortions upon Fourier transformation. The final 2D-ISRS map is generated by averaging all detected probe wavelengths over the red side of the stimulated emission.
The FC spectrum [ Fig. 2(c), blue] was obtained by Fourier transformation of the vibrational coherence in the absence of the Raman pulse. The A* spectrum [ Fig. 2(c), orange] was generated by first subtracting the Raman-probe vibrational coherences in the absence of the actinic pulse from the vibrational coherences in its presence prior to Fourier transformation along T 2 to yield the signal SðT 1 ; Ω 2 Þ and averaging over the first 300 fs along T 1 . The 2D-ISRS Fourier map in Fig. 3(a) was generated by applying the same subtraction procedure for consistency. We remark, however, that the ground-state contributions do not evolve along T 1 in the absence of the actinic pulse, resulting in similar 2D-ISRS maps with and without a ground-state subtraction procedure.

APPENDIX D: SIGNAL DERIVATION
We consider, as a reference model, a three-electroniclevel system with a ground state and two excited electronic states, S 0 , S 1 , and S n , as depicted in Fig. 1(c), with the associated vibrational manifolds. The total field E acting on the sample consists in the three delayed pulses described above.
The 2D-ISRS signal is given by a convolution between the field amplitudes and the matter correlation function F , derived from the four Feynman pathways shown in Fig. 4(a): Note that F can be expanded as a sum over states (SoS) [64] corresponding to the eigenfunctions of free molecule Hamiltonian H 0 , obtaining S ð5Þ ðω;T 1 ;T 2 Þ ¼ X g;g 0 ;f e;e 0 ;e 00 KðgÞIm μ g 0 e 00 μ Ã g 0 e 0 E Ã s ðωÞE s ðω − ω e 00 e 0 Þ 2ðω −ω e 00 g 0 Þ − μ g 0 e 00 μ Ã e 0 g 0 E s ðωÞE Ã s ðω þ ω e 00 e 0 Þ 2ðω þω g 0 e 0 Þ × e −iω ee 0 T 1 e −iω e 00 e 0 T 2 μ e 00 f μ Ã fe W p ðω e 00 e 0 − ω ee 0 ;ω e 00 e 0 −ω fe 0 Þμ e 0 g μ Ã eg ½W a ðω ee 0 ;ω ee 0 −ω eg Þ − W a ðω ee 0 ;ω ge 0 Þ ; ðD1Þ whereω ij ¼ ω i − ω j − iΓ ij , μ ij are the matrix elements of the dipole operator V, and the weighted pulse spectral densities W k for k ¼ a, p are defined as and rule the resonant bandwidth accessible by the finite widths of the actinic and Raman pulses. The subscripts e, e 0 , e 00 run over the vibrational levels of S 1 , while f and g, g 0 run over those of S n and S 0 , respectively. The initial thermal population is ruled by the Boltzman factor KðgÞ that depends on the temperature. At room temperature, vibrational modes above a few hundred wave numbers are initially populated only in the vibrational ground state; i.e., g in Eq. (D1) is fixed. The displacements d n , n ¼ 1, 2, 3, and the Duschinsky angle Θ D enter in the signal expression via the dipole matrix elements μ ij as detailed in Sec. d of the SM [54]. Specifically, we note that the functional dependence of the signal on each displacement is the same since the signal is proportional to a product of six dipoles and each d n appears in two of them. Finally, a Fourier transformation over the two delays and the integration on ω lead to the 2D frequency correlation map: dT 1 dT 2 e iΩ 1 T 1 þiΩ 2 T 2 S ð5Þ ðω; T 1 ; T 2 Þ: Using a WLC probe, the effective domain of integration over ω is restricted by the bandwidth of the pulse EðωÞ. In our experimental implementation, the measurement is spectrally resolved in ω, and a selective average over a specific spectral region can be exploited to isolate the resonant contributions from the stimulated emission or the excited-state absorption of the system. In particular, here the signal has been averaged over the tail of the stimulated emission region of the GFP, from 575 to 615 nm, while the corresponding maps averaged on the excited-state absorption are presented in Fig. S2 of the SM [54].
In the SoS picture, the signal is determined by the contributions originating from different pathways in the Liouville space, which corresponds to the different permutation of the SoS indexes. In this way, the contribution of selected modes to the signal can be easily isolated, overcoming interpretative issues due to interference between different nonlinear optical effects, in addition to giving the advantage of speeding up the calculations by removing the contributions to an unobserved region of the spectrum. In the simulations, we did not include the contributions to the signal lying on one of the two axes, originating from pathways in which populations instead of coherences evolve during T 1 or T 2 , since they provide the same information of third-order techniques and are suppressed by the experimental analysis routine. Further details on the signal derivation are provided in Secs. c and d of the SM [54].