Terahertz molecular motions are often probed by high-frequency molecular oscillators in different types of non-linear vibrational spectroscopy. Recently developed two-dimensional terahertz-infrared-visible spectroscopy allows direct measuring of this coupling and, thus, obtaining site-specific terahertz vibrational spectrum. However, these data are affected by the intensity and phase of the employed laser pulses. In this work, we develop a method of extracting sample response—representing solely physical properties of a material—from experimental spectra. Using dimethyl sulfoxide as a model molecule to verify this method, we measure the coupling between C–H stretch vibration of its methyl groups and terahertz intramolecular twist and wagging modes.

## INTRODUCTION

Molecular low-frequency vibrational modes (LFMs) in the terahertz (THz) frequency range can play a role in (bio-) chemical processes in condensed phases at room temperature.^{1–4} Non-linear mid-infrared (MIR) spectroscopy techniques (transient absorption and 2D IR) measure high-frequency vibrational modes (HFMs) to obtain insight into these motions at locations of specific chemical groups.^{5–10} Such measurements probe terahertz dynamics indirectly via LFM-HFM coupling, by exciting and probing only a high-frequency oscillator. Sophisticated relations between the HFM frequency/transition dipole moment and the LFM coordinate complicate determining the relevant LFM modes from such measurements. Recently, we have developed two-dimensional terahertz-infrared-visible (2D TIRV) spectroscopy to measure LFM-HFM coupling directly.^{11,12} This technique can identify the relevant LFMs using mid-infrared spectroscopy selectivity, thereby disentangling the typically broad featureless low-frequency spectral response.^{13–16} Previously, we have used this technique to study the coupling between intra- and inter-molecular motions in liquid water and the coupling between phonon vibrations and N–H stretch vibrations in hybrid perovskites.^{11,12} However, taking full advantage of the 2D TIRV spectroscopy requires extracting the system’s complex-valued third-order non-linear response function $S3$ from the measured data.

A measured 2D TIRV spectrum is determined by a product of $S3$ and spectra of laser pulses. The latter can be considered as an instrument response function (IRF), which affects both the intensity and phase of the 2D spectrum. Deriving $S3$, which contains the physical properties of a sample, thus requires eliminating IRF from the experimental data. This requirement is analogous to data processing in other non-linear spectroscopy techniques where the signal depends on the intensities and phases of the employed laser pulses. For example, the progress in 2D IR and 2D electronic spectroscopy was enabled by a method of determining the phase of a spectrum.^{17–21} The method utilizes the equality of 2D spectrum projection and transient absorption spectrum. Another technique, phase-resolved sum-frequency generation (PR SFG) spectroscopy, uses the response of reference material to obtain the second-order response function of a sample.^{22,23} Here, we develop a similar method for 2D TIRV spectroscopy.

We have recently reported a procedure that allows separating different excitation pathways in 2D TIRV spectroscopy.^{24} This was a critical step toward deriving material response function from experimental data because different excitation pathways produce signals via different mixing of the terahertz and infrared fields: sum- and difference-frequency mixing. Here, we report a method to eliminate the effect of the laser fields and obtain the spectrum of sample response function in 2D TIRV spectroscopy. To this end, we utilize a thin (1 *µ*m) silicon nitride (SiN_{x}) membrane as the reference material. Its response gives a spectrum that directly reflects the distribution of intensities and phases of all the laser pulses used in 2D TIRV spectroscopy to generate and detect signal field. We demonstrate the approach on a model sample, liquid dimethyl sulfoxide (DMSO), which has intense low- and high-frequency modes that are sufficiently strongly coupled to generate a substantial 2D TIRV signal.

## PRINCIPLES OF 2D TIRV SPECTROSCOPY

The theoretical formalism and experimental implementation of 2D TIRV spectroscopy have been discussed in detail elsewhere.^{11,24,25} In brief, in a 2D TIRV spectroscopy measurement, four-wave mixing of a narrowband visible (VIS), and broadband terahertz (THz) and mid-infrared (IR) pulses induce the emission of a signal field by a sample (Fig. 1). The resonant interactions of THz and IR pulses with the sample excite its LFMs and HFMs, respectively, with time delay *τ* between the pulses. VIS pulse adds third, non-resonant interaction between electromagnetic field and material to enable signal generation from isotropic media. Also, using VIS pulses leads to the signal being emitted in the visible frequency range, typically 620–720 nm, for which sensitive detectors are commercially available. Doubly resonant generation of signal—via excitation of both LFMs and HFMs—requires a two-quantum transition at one of the excitation steps and, thus, coupling between the modes. The signal field is dispersed in a spectrometer and measured using its interference with a local oscillator (LO) pulse at different time delays *τ* between the THz pulse and IR/VIS pulse pair. Fourier transformation of the time-domain data over *τ* provides the THz frequency axis (*ω*_{1}), representing LFMs. The second frequency axis, *ω*_{2}, representing HFMs, is obtained by subtracting the frequency of the VIS pulse from the spectrometer frequency. A constant time delay *δ* ≈ 4.1 ps between the signal and LO pulses is used to distinguish positive and negative *ω*_{2} frequencies and separate signals generated by rephasing (VIS + IR−THz difference-frequency mixing) and non-rephasing (VIS + IR + THz sum-frequency mixing) excitation pathways (quadrant separation).^{24} Rephasing and non-rephasing signals appear in the II/IV and I/III quadrants of the 2D spectrum, respectively.

^{24}

Here, $S3$ is the third-order response function of a sample, which reflects its physical properties; E_{VIS} is the complex-valued spectral amplitude of the VIS electric field with frequency Ω_{VIS}; and E_{THz}, E_{IR}, and E_{LO} are complex-valued electric field spectra of THz, IR, and LO pulses, respectively. The first and second terms in Eq. (1) represent different interaction sequences between the THz and IR fields and the sample. In the first (second) term, the interaction of the sample with THz (IR) precedes interaction with IR (THz). Equation (1) shows that the intensity and phase of a measured 2D TIRV spectrum are affected by the spectral shapes and phases of the employed laser pulses. We note that the THz, IR, and LO pulse phases can have complex dependence on frequency (due to, e.g., dispersion). Hence, a measured spectrum can be considered a product of the sample response function and instrument response function (IRF), which is given by the product of electric fields.

^{22}In this technique, a measured spectrum Γ

_{sample}of a sample is divided by a non-resonant spectrum Γ

_{ref}of a reference material (typically quartz). For non-resonant interactions, a material response function of any order is real-valued and constant (frequency-independent).

^{26}Thus, the non-resonant spectrum is linearly proportional to the IRF of an experimental setup and $\Gamma sample/\Gamma ref\u221dSsample3$. This approach can be employed using a resonant reference material as well, given that its response function $Sref3$ is known. In this case,

We demonstrate that in 2D TIRV spectroscopy, a thin silicon nitride (SiN_{x}) membrane is a suitable reference material to derive the complex-valued spectrum of the sample response function.

## RESULTS AND DISCUSSION

We use liquid dimethyl sulfoxide (DMSO) as a model sample. The absorption spectrum of DMSO in the mid- and far-infrared frequency range is shown in Fig. 2. In the terahertz frequency range, this spectrum contains two sharp vibrational modes at 333 and 383 cm^{−1} labeled 2 and 3. These vibrations are assigned to intramolecular CH_{3}–CH_{3} twist (*v*_{23}) and wagging (*v*_{11}) modes, respectively.^{27–29} There is also a weak shoulder at about 308 cm^{−1} (labeled 1), which is also produced by intramolecular motion and is assigned to C-S-C bending (*v*_{12}).^{29} Besides intramolecular vibrations, the spectrum contains intermolecular libration motion at frequencies below 150 cm^{−1} (labeled 0).^{27}

Because LFMs 2 (333 cm^{−1}) and 3 (383 cm^{−1}) of DMSO displace CH_{3} groups of the molecule, one can expect substantial coupling between LFMs and the C–H stretch vibrations. The C–H stretch symmetric and asymmetric vibrational modes have frequencies of 2913 and 2996 cm^{−1}, respectively (labeled with s and a in Fig. 2). Figure 3(a) shows the time-domain 2D TIRV data of DMSO measured in this IR frequency range using a femtosecond IR pulse with a bandwidth of 250 cm^{−1} and a central frequency of 2780 cm^{−1} (Fig. 2). Fourier transformation over time axis and quadrant separation yields 2D spectrum Γ_{DMSO}. Figure 3(b) shows the first and second quadrants of the absolute-value 2D TIRV spectrum ($\Gamma DMSO$). Sum-frequency mixing of the THz and IR fields produces a signal in the first quadrant at *ω*_{2}-frequencies of CH_{3} stretch vibrations.^{24} The offset of the IR pulse frequency from CH_{3} modes to red favors this process. The similar signal in the second quadrant is weak because it is produced by difference-frequency mixing of the IR and THz fields, which is hampered by the IR frequency offset. Instead, the second quadrant contains a signal parallel to the diagonal and intersecting *ω*_{2}-frequency axis at the CH_{3} stretch frequencies.

_{DMSO}, respectively. Their shape is affected by the spectral intensities and phases of the THz, IR, and LO pulses. Inter alia, an offset of time zero in time-domain data [Fig. 3(a)] and time delay

*δ*between signal and LO results in oscillations along

*ω*

_{1}and

*ω*

_{2}axes, respectively. To eliminate this effect of IRF, we measure the 2D TIRV spectrum Γ

_{SiN}of a reference material, a ∼1

*µ*m thick SiN

_{x}membrane [Figs. 3(e)–3(h)]. Because SiN

_{x}has no resonances at mid-infrared and visible frequencies, interaction with IR and VIS pulses is non-resonant for this material. In contrast, in the far-infrared frequency range, SiN

_{x}has intense vibrational resonances (Fig. 4). Hence, the signal can be generated by its non-resonant and resonant interaction with the THz pulse. In the former case, SiN

_{x}response function $SSiN3$ is real-valued and constant, whereas in the latter case, it is complex-valued and depends on the

*ω*

_{1}-frequency,

_{x}. Figure 4 shows the $SSiN1$ spectrum calculated using its reported permittivity

*ɛ*

_{SiN},

^{30}

_{0}is the number density of absorbers (we assume ρ

_{0}= 1). In the 0-700 cm

^{−1}THz frequency range of our 2D TIRV experiment, the $SSiN1$ amplitude varies by about 67% and its phase by up to 0.8 rad. This dispersion is primarily caused by the transverse optical phonon resonance around 825 cm

^{−1}.

^{31}

_{x}as a reference material, we test the two possible extreme cases of interaction with the THz pulse. Figures 5(a)–5(c) and 5(d)–5(f) show $SDMSO3$ in the first quadrant obtained assuming fully non-resonant (case I) and partly resonant (for the THz pulse interaction, case II) SiN

_{x}signal generation, respectively,

*ω*

_{1}≲ 400 cm

^{−1}. For higher

*ω*

_{1}-frequencies, the difference of the line shape becomes more significant because of the ∼0.8 rad accumulated phase of $SSiN1$. A 2D spectrum with the correct phase should fulfill the Kramers–Kronig relations (KKR) calculated parallel to one of the frequency axis,

^{26}for instance,

Figure 6 shows real and imaginary parts of the spectra in Fig. 5 taken at the symmetric stretch frequency (*ω*_{2} = 2913 cm^{−1}, mode s in Fig. 2). For both cases I and II, there is a significant discrepancy in the frequency range below $\u223c100cm\u2212\u20091$ between the imaginary part inferred from the real part using the KKR and experimental data [Figs. 6(b) and 6(d)]. This can be explained by two factors: low signal intensity in this spectral range and the limited spectral range (0–600 cm^{−1}) used in the calculations. The 2D TIRV signal from SiN_{x} reduces by a factor of ten and more at *ω*_{1} ≲ 20 cm^{−1}, thus making unreliable intensities of the derived $SDMSO3$ at these *ω*_{1}-frequencies. Because the KKR is not local, this can affect the inferred imaginary part in Fig. 6 in the frequency range significantly beyond 20 cm^{−1}. In addition, the correct KKR also requires integration over negative frequencies (quadrant II). In the calculations, we use data only from quadrant I (0–600 cm^{−1}) because of the low signal-to-noise ratio in quadrant II. This approximation can influence the inferred imaginary part, particularly at low *ω*_{1}-frequencies. Above 100 cm^{−1}, the inferred imaginary part agrees reasonably with the experimental data for case II, given that the calculated data in 500–600 cm^{−1} range can also be affected by the limited spectral range used in the KKR. For case I, the agreement between the calculated and experimental spectra is worse [Fig. 6(b)]. Thus, the 2D TIRV signal of SiN_{x} is generated predominantly via the resonant interaction with the THz pulse, and the correct spectrum of $SDMSO3$ is given by case II.

For each of the C–H stretch vibrational modes s and a, the $SDMSO3$ spectrum contains three peaks at *ω*_{1} ≈ 300, 326 and 376 cm^{−1} (s1, s2, s3 and a1, a2, a3 in Fig. 5) reflecting their coupling with the intramolecular LFMs 1, 2, and 3. We attribute the additional prominent peaks s4 and a4 to an LFM with a weak transition dipole moment and strong coupling with C–H stretch vibrations. This weakly absorbing mode is obscured in the absorption spectrum in Fig. 2. To determine which of the THz-IR and IR-THz interaction sequences produces these signals, we perform a Fourier transformation of the spectrum over the *ω*_{1}-frequency axis back into the time domain [Fig. 7(a)]. Because the data processing eliminates phases of the laser pulses, $SDMSO3$ has a maximum at *t* = 0 (note that here the *t*-axis is uniquely determined by the phase of the frequency-domain spectrum). We use a square window function to select data at negative (*t* < 0) and positive (*t* > 0) time and perform their Fourier transformation separately to derive the signal generated by the IR-THz and THz-IR interaction sequences, respectively. The spectrum of the signal at *t* < 0 contains peaks that are broad over the *ω*_{1} axis and narrow over the *ω*_{2} axis [Figs. 7(b)–7(d)]. Thus, the IR-THz interaction sequence is resonant for the IR interaction and presumably non-resonant for the THz interaction. In contrast, the spectrum of *t* > 0 contains a signal generated by the coupled intramolecular LFMs and C–H stretch vibrations [Figs. 7(e)–7(g)]. Hence, the doubly resonant signal is produced by DMSO molecules interacting first with the THz pulse, followed by interaction with the IR pulse.

Cross-peaks between CH_{3} stretch vibrations and intra-molecular LFMs 1, 2, and 3 are at frequencies of *ω*_{1} ≈ 300, 326, and 376 cm^{−1}. These frequencies are systematically lower by about 5 cm^{−1} than LFMs locations in the absorption spectrum (accuracy of the *ω*_{1} and *ω*_{2} frequencies is about 1 cm^{−1} each; see the supplementary material for details). This deviation cannot be explained by different frequencies of LFMs in the excited state of the C–H oscillator because the THz pulse excites LFMs in the ground vibrational state of the C–H stretch. Because the first excited state of the DMSO low-frequency intramolecular vibrations is populated at room temperature, the THz pulse excites molecules also from the first to the second excited state. One can hypothesize a significantly stronger coupling between the second excited state of THz modes and C–H stretch vibrations, which can dominate the 2D TIRV spectrum. Thus, the $\u22485cm\u2212\u20091$ offset of resonances in the 2D TIRV spectrum along the *ω*_{1}-axis can stem from the own mechanical anharmonicity of DMSO intramolecular vibrations. To verify this hypothesis, we measure the absorption spectrum of the DMSO in a 1 mm thick sample cell to identify possible overtone transitions of the molecule in the THz frequency range (Fig. S4 in the supplementary material). A weak resonance at about 765 cm^{−1} can be tentatively assigned to the overtone transition of the 383 cm^{−1} vibrational mode, which implies $\u22481cm\u2212\u20091$ own mechanical anharmonicity of this LFM. An alternative explanation of the red-shifted response can be a slight inhomogeneous broadening of the LFMs. For instance, LFM frequencies can have a weak dependence on the instantaneous conformation of DMSO molecules, such as the rotational orientation of its methyl groups around S–C bonds. Assuming LFM-HFM interaction strength also depends on instantaneous DMSO conformation, this can lead to stronger coupling between red-shifted LFM oscillators and C–H stretch vibrations. Hence, an offset of resonance frequencies in the 2D TIRV spectrum can occur as compared to linear absorption. However, a comprehensive explanation of this phenomenon requires accurate theoretical calculations.

At even lower *ω*_{1} frequency, the $SDMSO3$ spectra display coupling between the C–H stretch vibration and intermolecular librational motions [s0 and a0 in Fig. 5(b)]. The librational motion changes the spatial orientation of the C–H stretch transition dipole moment, and because the intermolecular forces on the CH_{3} groups are typically rather weak, this coupling is likely produced by predominantly electrical anharmonicity. The maximum intensity in the 2D spectrum is substantially below a THz frequency of 100 cm^{−1}—the maximum of linear absorption (Fig. 2). We hypothesize that this can be due to the larger amplitude (angle) of libration at lower frequencies and, thus, stronger LFM-HFM coupling.

In summary, we present a method of deriving a complex-valued sample response function in experimental 2D TIRV spectroscopy. The sample response free from the instrument response function reflects its physical properties and can be directly compared with theoretical calculations. For a model liquid DMSO sample, we directly measure the coupling between high-frequency C–H stretch modes and THz molecular motions, both intramolecular vibrations and intermolecular libration. These results pave the way for the utilization of 2D TIRV spectroscopy for the quantitative study of vibrational dynamics in liquids and soft materials important in different fields of chemistry.

## METHODS

The experimental setup is described in detail elsewhere.^{24} Here, we provide a brief overview. The output of a femtosecond Ti: sapphire regenerative amplifier (Astrella, Coherent) centered at 800 nm with a repetition rate of 1 kHz, and FWHM of ∼60 nm is split into three beams. The first beam (≈1 mJ/pulse) pumps a commercial TOPAS (traveling wave optical parametric amplification of superfluorescence) coupled with an NDFG (non-collinear difference-frequency generation) unit to generate tunable mid-infrared (IR) pulses with FWHM of ∼250 cm^{−1} (spectrum of IR pulse centered at 2780 cm^{−1} is shown in Fig. 2). We use a 4 f pulse shaper to produce a narrowband visible pulse (VIS, centered at about 800 nm, FWHM ≈ 10 cm^{−1}) from the second beam (≈0.8 mJ/pulse). IR and VIS beams are made collinear using a dichroic beam combiner. The IR/VIS pulse pair is then routed to a displaced Sagnac interferometer to generate a local oscillator (LO) by sum-frequency generation in a 10 *μ*m thick type 1 β-barium borate (BBO) crystal. The third beam (≈1.1 mJ/pulse) is used to generate a broadband terahertz (THz) pulse by two-color mixing of laser field in air plasma (spectrum of THz pulse is shown in Fig. 2). IR, VIS, and LO pulses are focused by a parabolic mirror and overlapped with horizontally polarized terahertz (THz) pulse at the sample to generate the third-order signal. After the sample, the signal and LO are collimated with a CaF_{2} lens (20 cm focal length) and passed through a polarizer. The signal and LO are aligned to a spectrometer (SpectraPro HRS-300, Princeton Instruments) and detected with an electron-multiplying charge-coupled device (EMCCD) detector (Newton 971, Andor). The energies of IR and VIS pulses at the sample are 1.5 and 4.5 *μ*J/pulse, respectively. We use parallel polarization of THz, IR, VIS, and signal fields.

We measure time-domain data by varying the time delay between the THz and IR/VIS pulse pair in steps of 10 fs using a motorized delay stage [V-551.2B, Physik Instrumente (PI)]. The starting position of the delay stage is the same for DMSO and SiN_{x} membrane. We scan time delay across a total time interval of 7.33 and 1.33 ps for DMSO and SiN_{x}, respectively. Data acquisition for one spectrum takes about 25 min for DMSO and 5 min for the SiN_{x} membrane.

The sample cell for DMSO consists of a front and a back window separated by a 1 mm thick Viton o-ring. The front window is made of a 200 *μ*m thick silicon frame with a low-stress SiN_{x} membrane of 50 nm thickness and 0.5 × 0.5 mm^{2} aperture (NX5050A, Norcada). The back window is a 2 mm thick CaF_{2} window. The 2D TIRV signal from the front window is about 30 times weaker than the DMSO signal and, thus, is neglected in data interpretation and analysis (see Fig. S5 in the supplementary material).

To measure the absorption spectrum in far- and mid-infrared, we use a commercial Fourier transform infrared (FTIR) spectrometer (Bruker Vertex 70), which has a spectral range of 50–6000 cm^{−1}. Liquid DMSO is placed between two diamond windows (400 *µ*m thick) with a 25 *µ*m thick PTFE spacer.

## SUPPLEMENTARY MATERIAL

See the supplementary material for calibration of ω1 and ω2 frequencies, infrared absorption spectrum of 1 mm thick DMSO sample, and comparison of 2D TIRV signal intensity of DMSO and 50 nm thick SiNx membrane.

## ACKNOWLEDGMENTS

We acknowledge the Max Planck Society for financial support.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Pankaj Seliya**: Data curation (lead); Formal analysis (lead); Investigation (lead); Visualization (lead); Writing – original draft (equal). **Mischa Bonn**: Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Methodology (supporting); Project administration (supporting); Supervision (supporting); Visualization (supporting); Writing – original draft (equal). **Maksim Grechko**: Conceptualization (lead); Formal analysis (supporting); Investigation (supporting); Methodology (lead); Project administration (lead); Supervision (lead); Visualization (supporting); Writing – original draft (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

_{2}O

_{2}O

*X*

^{(3)}for doubly vibrationally enhanced four wave mixing spectroscopy

*Springer Series in Optical Sciences*

*d*

_{6}

_{x}:H Prepared by plasma-enhanced chemical-vapor deposition