Two-Dimensional Terahertz Spectroscopy of Condensed-Phase Molecular Systems

Nonlinear terahertz (THz) spectroscopy relies on the interaction of matter with few-cycle THz pulses of electric field amplitudes up to megavolts/centimeter (MV/cm). In condensed-phase molecular systems, both resonant interactions with elementary excitations at low frequency such as intra- and intermolecular vibrations and nonresonant field-driven processes are relevant. Two-dimensional THz (2D-THz) spectroscopy is a key method for following nonequilibrium processes and dynamics of excitations to decipher the underlying interactions and molecular couplings. This article addresses the state of the art in 2D-THz spectroscopy by discussing the main concepts and illustrating them with recent results. The latter include the response of vibrational excitations in molecular crystals up to the nonperturbative regime of light-matter interaction and field-driven ionization processes and electron transport in liquid water.


I. INTRODUCTION
Molecular systems in the condensed phase display a variety of low-energy excitations in the frequency range from some 10 gigahertz (GHz) to 30 terahertz (THz). Collective molecular motions, inter-and intramolecular vibrations and coupled nuclear-electronic degrees of freedom such as soft modes give rise to dipole-allowed absorption and/or Raman bands, most of which exhibit complex and partly overlapping line shapes. While linear dielectric, far-infrared, and Raman spectroscopies have been developed over decades and reached a high level of sophistication, 1,2 nonlinear THz methods have opened a new field of research, in which the nonlinear low-frequency response, the intrinsic dynamics of excitations, their couplings, and the interaction with a thermal bath are addressed most directly. [3][4][5][6][7] This rapidly developing field has benefitted from major progress in generating THz few-cycle pulses with electric field amplitudes reaching the MV/cm range [8][9][10][11] and from the introduction of new nonlinear spectroscopies in which THz pulses induce a nonlinear response of a solid or a molecular ensemble.
The concept of two-dimensional correlation spectroscopy goes back to nuclear magnetic resonance 12 and has been transferred to the optical domain to study vibrational excitations in the mid-infrared range and/or electronic and vibronic excitations in the visible and ultraviolet. Both pump-probe and three-pulse photon-echo methods have been applied to follow the nonlinear response of a molecular system in third order in the optical field, mostly under conditions of resonant excitation. Detailed accounts of 2D third-order spectroscopy have been given in Refs. [13][14][15][16]. Two-dimensional THz (2D-THz) spectroscopy, introduced some 10 years ago, [17][18][19] goes well beyond the third-order limit. It allows for mapping nonlinearities up to the regime of nonperturbative light-matter interaction, where the coupling to the external THz electric field is comparable to or even stronger than interactions within the a) Electronic mail: reimann@mbi-berlin.de b) Electronic mail: woerner@mbi-berlin.de c) Electronic mail: elsasser@mbi-berlin.de molecular system. 7 Both resonant and nonresonant interactions as well as field-driven processes of charge transport have been studied by 2D-THz spectroscopy, which is inherently phase-resolved and gives access to absolute rather than relative optical phases.
Early work in 2D-THz spectroscopy has focused on method development and applications to electronic excitations in solids. Here, the large transition dipoles, e.g., of intersubband transitions in semiconductor quantum wells, 17,18,20 facilitated the generation of a quasi-resonant nonlinear response with moderate peak electric fields in the range of 1 to 50 kV/cm. More recent applications of 2D-THz spectroscopy in spectroscopic and transport studies of solids, metamaterials, and semiconductor devices have been reported in Refs. 21-34. The rapid development of generation schemes for THz pulses with MV/cm amplitudes has strongly widened the application range of 2D-THz methods, now including vibrational and/or phononic excitations and, in particular, a broad range of nonresonant nonlinear interactions. This article gives an account of the current state of this field with a focus on molecular excitations and processes in the condensed phase. It combines a presentation of the method and experimental aspects with a discussion of recent prototypical applications.
The content of this article is organized as follows. Section II summarizes methods for THz pulse generation in laser-driven sources (Sec. II A) and field-resolved detection (Sec. II B), including a short description of linear THz spectroscopy. Section III introduces the basic concepts of 2D-THz spectroscopy and the experimental implementation (Sec. III A), complemented by a discussion of data handling and analysis (Sec. III B). A comparison to related 2D methods such as nonlinear 2D Raman-THz and multicolor 2D spectroscopy is included as well (Sec. III C). A prototypical study of softmode excitations in molecular crystals of aspirin is presented in Sec. IV, followed by a discussion of ionization and electron transport in water driven by strong THz fields (Sec. V). Conclusions and an outlook on future developments are given in Sec. VI. Nonlinear THz spectroscopy relies on ultrashort THz pulses with an electric field amplitude sufficient to induce a nonlinear response in the sample under study. Such pulses have been generated with electrically biased photoconductive switches and antenna structures or by nonlinear optical frequency conversion of optical pulses from the visible and near-infrared to the THz range. Irradiation of a photoconductive switch made from GaAs, InP or other semiconductors by a femto-to picosecond optical pulse induces short-lived transient currents, which radiate electric fields at THz frequencies. The maxi-mum amplitude of the THz field is set by the DC bias field of the switch and reaches values on the order of 10 kV/cm for large-area emitters. While photoconductive switches driven by femtosecond laser oscillators are widely applied in commercial devices for linear THz spectroscopy, their applicability in nonlinear experiments is limited because of the comparably small peak electric fields.
Plasmas generated with intense femtosecond laser pulses in dilute gases emit strong few-cyle THz pulses. 8,37,38 Twocolor generation of a plasma in air or nitrogen, e.g. by the fundamental and second-harmonic frequency of an amplified femtosecond laser pulse, generates a unidirectional transient plasma current of electrons, which radiates a very broad frequency spectrum including pronounced THz components (up to 70 THz). 39 Peak electric fields of up to 400 kV/cm have been generated at a kilohertz repetition rate with 25-fs pulses from an amplified Ti:sapphire laser system. 8 This type of source has been applied in a number of nonlinear THz experiments, although plasma instabilities result in fluctuations of the electric field amplitudes and the spatial beam parameters.
Optical difference frequency mixing or rectification in nonlinear crystals and/or at surfaces with a high second-order nonlinearity χ (2) is the predominant method for generating few-cycle THz pulses with high electric field amplitudes. Here, the difference frequency of two synchronized incoming pulses or between different components in the broad spectrum of a single femtosecond pulse is generated, preferentially in a phase-matched conversion process, in which the phase velocity of the THz pulse needs to match the group velocity of the pump pulses. Beyond the THz frequency range discussed here, this method has been widely applied to generate pulses in the mid-infrared up to frequencies on the order of 100 THz.
Inorganic semiconductors with a zincblende crystal structure such as GaAs, GaP, InP, CdTe, and ZnTe and the ferroelectric LiNbO 3 are standard materials to cover a frequency range up to some 6 THz with electric field amplitudes 7 up to several hundred kV/cm. The introduction of noncollinear rectification in LiNbO 3 using pump pulses with a tilted wavefront 10,40 allows for reaching peak electric fields of several MV/cm at frequencies between 0.2 and 1.5 THz. The time-dependent electric field of a THz pulse generated with this method in shown in Fig. 1(a), the corresponding intensity spectrum in Fig. 1(b). Drawbacks of this method are the spatial and temporal break-up of the optical pump due to angular dispersion, which limits the conversion efficiency, and the limited quality of the generated THz beam. Very recent work has addressed such issues both by numerical simulations of different interaction geometries and experiments. 41 It should be noted that few-cycle pulses with electric field amplitudes up to 100 MV/cm have been generated 9 in the frequency range from 10 to 70 THz.
Organic nonlinear crystals consisting of molecular chromophores have received substantial interest for THz generation because of their high second-order nonlinearities. 42,43 Figures 1(c) and (d) display the time-resolved electric field and the intensity spectrum of pulses generated by driving a 400-µm thick DSTMS crystal (DSTMS: 4-N,N-dimethylamino-4'-N'-methylstilbazolium 2,4,6-trimethylbenzenesul-fonate) with 25-fs pulses from an amplified Ti:sapphire laser. 44 The THz spectrum covers the range from 0.1 to 30 THz. The narrow features in the spectrum originate from radiation on vibrational resonances and correspond to the long-lived oscillations of the electric field in the time-domain. The frequency conversion process is only partially phasematched. However, due to the very high optical nonlinearity, an appreciable conversion efficiency is possible even without phase-matching. The comparably low optical damage thresholds of organic nonlinear crystals at higher pulse repetition rates eventually limit the attainable electric field amplitudes.

B. Field-resolved detection and measurement
Phase-resolved detection of THz pulses, i.e., the measurement of the electric field as a function of time, is a key ingredient of 2D-THz spectroscopy. The standard method is freespace electrooptic (EO) sampling, which exploits the change of refractive index induced by the THz field in an electrooptic crystal. 45 This change is mapped by the polarization change of a femtosecond probe pulse, the duration of which needs to be much smaller than the THz period. Measuring the polarization change for different delay times between the THz transient and the probe pulse allows for reconstructing the time-dependent THz field and, in calibrated setups, the absolute electric field amplitudes. Widely used electrooptic materials include ZnTe, GaP, GaSe, LiTaO 3 , and LiNbO 3 . A detection bandwidth on the order of 100 THz has been demonstrated with 10 to 50 µm thick ZnTe and GaSe crystals. 46,47 The smallest electric fields measured by EO sampling are on the order of 1 V/cm, as shown in studies of vacuum fluctuations of the electric-field. 48,49 Further details and literature can be found in Ref. 7. A setup for EO sampling is presented in Sec. III A (Fig. 4).
Time-domain THz methods have found broad application in linear mid-and far-infrared spectroscopy. A field-resolved detection scheme allows for deriving both the real and imaginary part of the linear dielectric function from the THz field transmitted through a sample. Moreover, the fact that electric fields are detected rather than intensities as in conventional far-infrared spectroscopy makes the study of optically thick samples possible. The electric field amplitude decreases only proportional to the square root of the intensity, e.g., a reduction of intensity to 1/100 corresponds to a reduction of 1/10 in electric field amplitude only. This aspect is relevant for aqueous solutions, e.g., of biomolecules where the background absorption of water represents a major contribution to the overall absorption spectrum. To illustrate the potential of fieldresolved spectroscopy, Fig. 1 value of approximately 2.8, corresponding to a very small intensity transmission of 1.5 × 10 −3 , which is typically beyond the detection range of conventional Fourier spectrometers.

A. Concept and experimental implementation
Two-dimensional THz spectroscopy is based on the interaction of a pulse sequence consisting of two or three THz pulses with a sample and a measurement of the total transmitted THz electric field in amplitude and phase. If the response of the sample is perfectly linear, the total transmitted field E lin ABC of, e.g., three pulses (A, B, C) is equal to the sum of the electric fields of the individual pulses: Accordingly, 2D-THz spectroscopy only yields additional information in case of a nonlinear response of the sample. In 2D-THz spectroscopy, one exploits the nonlinearities resulting from interaction with all pulses applied. For two pulses [ Fig. 2(a)], the nonlinear response determines the electric field E NL (τ,t) given by 17 In this equation, τ is the time delay between the two pulses A and B, t the real time, i.e., the time axis along which the electric field is measured by electrooptic sampling, E A (τ,t) the transmitted electric field when only pulse A is incident on the sample, E B (t) the transmitted field when only pulse B is incident on the sample, and E AB (τ,t) the transmitted field with both pulses present. The single-pulse nonlinearities are contained in E A (τ,t) and E B (t), so that they are removed by Eq. (2). For three pulses [ Fig. 2(b)], one has to subtract both oneand two-pulse nonlinearities to recover E NL (τ, T,t) according to 50 The delay between pulses B and C is T , often called waiting time. The meaning of E ABC (τ, T,t) etc. is analogous to the two-pulse case. For two pulses, the nonlinear susceptibilities contributing to E NL are χ (2) , χ (3) , and higher, for three pulses they are χ (3) , χ (4) , and higher. In general, an n-pulse scheme of this type will show nonlinearities of order n and higher. Nonlinearities of even and odd orders can be distinguished by looking at the spectrum of the nonlinear field. For input fields of a frequency ω even-order nonlinearities will result in nonlinear fields with frequencies 0, 2ω, 4ω, . . ., whereas odd-order nonlinearities generate fields with ω, 3ω, 5ω, . . .. Thus, a nonlinear signal with the frequency ω in a two-pulse experiments can not result from a second-order nonlinearity but originates from a third-or higher odd-order nonlinearity. Most third-order nonlinearities can be measured with only two pulses.
The different contributions to the total nonlinear signal can be analyzed and separated with density matrix theory describing the nonlinear light-matter interaction. 51 The signal components correspond to different pathways in Liouville space. Such pathways have been visualized by diagrammatic techniques such as double-sided Feynman diagrams, 51 multi-level schemes, 52 or sequences of arrows in 2D frequency space. 19 Reproducible femtosecond delays between the THz pulses and synchronization with the read-out pulse for electrooptic sampling are essential for 2D-THz spectroscopy. Typically, all pulses are derived from a single mode-locked laser oscillator. There are various schemes to generate the two or three THz pulses for 2D-THz spectroscopy as schematically illustrated in Fig. 3 for two pulses. In scheme I, a single THz pulse is generated, split in two replicas, and recombined after variable time delays. In scheme II, the beam from the pump laser is split and then sent to two separate THz generators G. In this way it is possible to have different frequencies of the two beams for two-color 2D spectroscopy. Scheme III requires only a single THz generator G, but has the disadvantage that the generation process leads to nonlinear electric fields according to the definitions (2) and (3), in particular during the time overlap of the pump pulses.
Key components of schemes I and II are the THz beam splitters S and D in I and D in II. An ideal beam splitter should generate a reflected pulse with E r (t) = a E(t − δ r ) and a trans- Here, δ t and δ r are time delays upon reflection and transmission. Components coming close to this ideal are thin pellicle beam splitters, plate-type beam splitters, and wire grid polarizers. The thin plastic substrate of a pellicle introduces weak dispersion-induced changes of the transmitted transients, but may display considerable absorption at certain THz frequencies. Disadvantages of pellicles are their fragility and their high sensitivity to acoustic noise and air currents. Thicker plate-type beam splitters need to possess low THz dispersion and absorption. Suitable materials are undoped diamond, silicon, and germanium. Unwanted reflections from the back surface of beamsplitter plates can be suppressed by antireflection coatings. [53][54][55][56] Another solution is to use the beam splitter under Brewster's angle for p polarization and add a metallic coating on one side for the required reflectivity. 57 Wire grid polarizers consist of parallel thin metallic wires. For wavelengths large compared to the distance between neighboring wires, they transmit electric fields perpendicular to the direction of the wires and reflect electric fields parallel to the wires. After a wire grid beamsplitter D, the two beams will be orthogonally polarized. To obtain the same polarization, one can introduce an additional polarizer between D and the sample.
In both schemes I and II, the beam splitters D can be replaced by having the two beams parallel and close to each other (lower part of Fig. 3). Because of the generally high divergence of THz beams, they will overlap after some propagation length. For instance, an initial beam waist of 5 mm (typical aperture of nonlinear crystals, e.g., GaSe) at a frequency of 1 THz (wavelength λ = 300 µm) results in a minimum divergence of 38 mrad. After a propagation length of 30 cm the two beams overlap nearly perfectly. This method can easily be extended to three beams arranged, e.g., in a triangular pattern (Fig. 4). In such geometries, essentially no losses occur, i.e., in scheme I for two (three) beams each beam has, upon hitting the sample, one half (one third) the energy of the initially generated THz pulse outside the time overlap of the two pulses. In contrast, the pulse energies on the sample is only one quarter (one ninth) of the initial energy for ideal beam splitters D.
A setup for three pulses using scheme II is shown in   50 The pump pulses are generated with a Ti:sapphire oscillator-amplifier system. This particular setup includes two multipass amplifiers, both seeded from the same oscillator to ensure synchronization. In this scheme, the shape of the two amplified pulses can be set independently with acousto-optic pulse shapers, permitting both phase and amplitude changes. The settings are chosen for optimum generation of the required THz pulses, 58 a feature particularly important for twocolor measurements. The setup of Fig. 4 uses pulses from the oscillator as readout pulses in electrooptic sampling. They are shorter than the amplified pulses (12 fs versus 20-25 fs) and, thus, allow for measuring THz pulses at higher frequency (40 THz vs. 25 THz). 7,59 Moreover, the smaller pulse-to-pulse fluctuations improve the signal-to-noise ratio. While the temporal jitter between oscillator and amplifier output is typically on the order of a few femtoseconds, temporal drifts of the read-out relative to the THz pulse represent a challenge in this scheme. Due to the long optical path length through the amplifiers (on the order of 10 m), even a very small temperature difference between the path of the oscillator pulse and the path of the amplified pulse can lead to significant drifts. For example, a temperature difference of 0.1 K between parts of the optical table will lead to a change of the distance of 10 µm, equivalent to a time difference of 33 fs. It is, however, possible to correct for such drifts after the measurement (see Sec. III B).
The read-out and the THz pulse overlap in a thin (thickness ≈ 10 µm) electrooptic ZnTe crystal of (110) orientation.
The thin crystal is mounted onto a 0.5 mm-thick (100) ZnTe crystal, which shows no electrooptic effect. 60 In this way, the thin electroptic crystal is mechanically more stable, and, even more important, the time range over which reflections from the back surface of the crystal are absent, is extended. The THz electric field changes the dielectric tensor of the (110) ZnTe crystal, which in turns leads to a change of the polarization state of the initially linearly polarized read-out pulse. This change is converted by the quarter-wave plate and the polarizer into a difference between the signals of two balanced photodiodes. The difference signal is proportional to the THz electric field at this instant of time. It is digitized for every laser shot and stored in a PC, which also detects the settings of all choppers (see below) and moves the various delays according to a predetermined sequence.
The THz generation stages, the sample, and the EO sampling setup are placed in a closed chamber, which is evacuated with a scroll and a turbomolecular pump for removing gases such as CO 2 and water vapor. In this way, absorption and temporal distortions of the THz pulses are suppressed. At a pressure of 10 −3 mbar, the gas concentrations are reduced from their values at atmospheric pressure by six orders of magnitude. At pressures of the order of 10 −6 mbar it is possible to perform measurements at low temperatures by mounting the sample on the cold finger of a helium-flow cryostat.
Mechanical choppers synchronized to the 1-kHz repetition rate of the Ti:sapphire amplifiers are introduced in all beams to distinguish the THz pulses in the pulse sequence. The chopper in beam A works with half the repetition rate (500 Hz), B with one quarter (250 Hz), and C with one eighth (125 Hz). In this way chopper A transmits one pulse and blocks one, B transmits two and blocks two, and C transmits four and blocks four (lower part of Fig. 4). In this way, all possible combinations of the three pulses A, B, and C are measured. Varying the delays T , τ, and t one obtains the single-pulse transients E A (T, τ,t), E B (T,t), and E C (t), the two-pulse transients E AB (T, τ,t), E BC (T,t), and E AC (τ,t), the three-pulse transient E ABC (T, τ,t), and the nonlinear signal [Eq. (3)].

B. Data analysis
Two-dimensional THz spectra are derived from the electric field E NL (τ, T,t) at a fixed waiting time T by a 2D Fourier transform along τ and t, giving the signal E NL (ν τ , ν t ) as a function of the excitation frequency ν τ and the detection frequency ν t . For illustrating this concept and discussing key aspects of data analysis, we use a 2D-THz data set recorded in the three-pulse setup of Fig. 4 with the narrow-gap semiconductor InSb. 50,61 InSb is a direct-gap III-V semiconductor with a band gap of E g = 0.17 eV, corresponding to a frequency of ν = E g /h = 41 THz. 62 There is a TO two-phonon resonance at 10 THz, as observed in the higher-order Raman spectrum. The THz pulses with a center frequency of 21 THz are neither resonant to the band gap nor to the two-phonon resonance. The THz peak field strength E ≈ 50 kV/cm and the extraordinarily large interband transition dipole d cv = e · (4 nm) at the Γ . In particular, the nonperturbative regime allows for multiple interactions of a single THz pulse with the InSb crystal. In Refs. 50 and 61, we reported the ultrafast dynamics of two-phonon coherences and signals related to multiple two-photon interband excitation of electron-hole pairs. In the following, we concentrate on the dynamics of two-TO-phonon coherences in InSb.
In Fig. 5(a), the sum of electric field transients E A (τ, T,t) + E B (T,t)+E C (t) transmitted through the InSb sample is shown in a contour plot as a function of the coherence time τ and real time t for the waiting time T = 35 fs. The corresponding nonlinear signal E NL (τ, T,t) derived from Eq. (3) is plotted in panel (b). In accordance with causality, a nonvanishing E NL (τ, T,t) sets in only with or after the last pulse in the timing sequence. For better orientation, the center of pulse A is indicated by the orange dashed line, which intersects the horizontal τ = 0 line (black) at t = −T . In Fig. 5(c), the contour plot of the amplitude |E NL (ν τ , ν t )| is displayed for T = 35 fs.
The circles of different size and color indicate the positions of relevant signals in the 2D frequency space. As a function of detection frequency ν t , strong nonlinear signals occur in the spectral range of the driving pulses 15 THz < ν t < 25 THz (large circles) and at the two-phonon resonance ν t = 10 THz (small circles). The former signals have been discussed in detail in Ref. 61. In the following we focus on the small circles at the two-phonon resonance ν t = 10 THz. As a function of excitation frequency ν τ we observe significant nonlinear signals for ν τ = 0 (blue circles), ν τ = −ν t (black circles), and ν τ = +ν t (red circles). The signal at (ν τ , ν t ) = (0, 10) THz is a nonlinear signal of at least 11th order in the THz field, as has been discussed in detail in Ref. 50. For τ > 0, i.e., for the pulse sequence (A(τ, T,t), B(T,t), C(t)) pulse A induces a complete two-photon absorption event, i.e., it creates an electron-hole pair by four interactions with its electric field. After the coherence time τ, pulse B creates a two-phonon coherence via an impulsive third-order Raman excitation, which evolves during the waiting time T . Pulse C generates a second electron-hole pair, again requiring four interactions and thereby transferring the two-phonon coherence to the final electronic state. After the third pulse, the two-phonon coherence evolves along the real time t and emits radiation via its optical transition dipole.
Selection of the peak at (ν τ , ν t ) = (0, 10) THz by a 2D-Gaussian filter and subsequent Fourier back-transform into the time-domain gives the signal field E 2-phonon (τ,t), which is plotted in Fig. 5(d) as a function of τ and t. Panel (f) shows a cut of this signal field for τ = 500 fs. The Fourier filtering by the 2D Gaussian results in an extraordinarily low noise level below < 0.1 kV/cm. For comparison, a cut of the total signal field E NL (τ,t) for τ = 500 fs and T = 35 fs is shown in Fig. 5(e). Here, a noise level of ≈ 2 kV/cm is estimated from the signal at t < 0, which should vanish for causality reasons. This analysis shows that 2D-Fourier filtering can considerably increase the signal-to-noise ratio of narrow signal peaks in the 2D-THz spectrum. To further illustrate the power of 2D-Fourier filtering, we added artificially white noise to the measured E NL (τ,t) shown in panel (b), which results in the grey dots in panels (e) and (f). Here, the total nonlinear signal E NL (τ = 500 fs,t) is barely detectable. In contrast, the 2D-filtered nonlinear two-phonon coherence E 2-phonon (τ = 500 fs,t) has a noise level of ≈ 0.3 kV/cm [grey dots for t < 0 in panel (f)] allowing for a clear experimental characterization of this nonlinear contribution.
We now address the correction of long-term delay drifts between the three THz pulses A, B, and C derived from the output of the Ti:sapphire amplifiers in Fig. 4 and the read-out pulse from EO sampling which is derived from the modelocked oscillator. In Fig. 5(g), the nonlinear signal field E NL (τ set ,t set ) is plotted as a function of the uncorrected times τ set and t set , as given by the computer controlling the delay stages. The time drifts in this raw signal lead to slight distortions of the shape and the period of the signal wavefronts. Such distortions can easily be corrected as the timing of the individual three pulses is measured simultaneously during a complete 2D scan, applying the pulse sequences shown in the lower part of Fig. 4. Using a 2D linear interpolation, one can reconstruct the nonlinear signal on a regular equidistant (τ,t) 2D grid resulting in the 2D signal E NL (τ,t) shown in Fig. 5(b).
The separation of nonlinear signals with respect to their order in the THz field and the assigment to different pathways in Liouville space (cf. Sec. III A) requires distinct nonoverlapping signal peaks in the frequency-domain 2D spectrum. In the nonperturbative regime of both resonant and nonresonant light-matter coupling, this condition may not be fulfilled. Resonantly excited Rabi flops, e.g., of intersubband transitions, 17 include signal components up to a very high nonlinear order and display overlapping 2D signal peaks at high driving fields. Electrons coherently driven in the conduction band of a semiconductor by a strong nonresonant THz field display a nonperturbative response and emit a multitude of harmonics of the fundamental driving frequency which are difficult to separate. 64,65 In general, an interplay of resonant and nonresonant contributions to the overall nonlinear response results in highly complex 2D-THz spectra. Here, variation of the phase and polarization of the THz pulses and, in anisotropic samples, measurements with different sample orientations can help to decipher the 2D spectra.

C. Related 2D Methods
There is a number of related 2D methods which involve at least one THz pulse in a sequence of ultrashort pulses and/or rely on the detection of coherent THz emission from an excited sample. In the present context, 2D Raman-THz, 66,67 2D THz-Raman, [68][69][70] and THz-infrared-visible spectroscopy 71 are relevant, all so far being applied in the third-order limit of light-matter interaction.
Two-dimensional Raman-THz spectroscopy of liquids gives insight in the dynamics and couplings of intermolecular degrees of freedom. The method involves a femtosecond pulse in the visible or near-infrared for the Raman interaction, which is second order in the electric field, and a THz pulse which induces a coherent polarization by a single interaction with the transition dipole of a low-frequency excitation. The THz emission of the sample represents the third-order nonlinear signal, which is measured as a function of real time t 2 in amplitude and phase by EO sampling. The second relevant time axis is the temporal separation t 1 of the Raman and the THz pulse.
Two interaction sequences have mainly been applied in experiments on liquid water and aqueous salt solutions. 66,67,[72][73][74][75] For positive t 1 , The Raman pulse interacts with the sample first and generates a Raman coherence, which is translated to a dipole-induced coherence by an interaction with the delayed THz pulse. The latter coherence gives rise to THz dipole emission, the nonlinear signal. For negative t 1 , the THz pulse comes first and generates a dipole coherence, which is switched to a Raman coherence by two interactions with the delayed Raman pulse. The THz signal field is emitted at a time t 2 after the THz pulse.
Two-dimensional Raman-THz spectroscopy of water and aqueous salt solutions has revealed echo-like signals in the range of water-water hydrogen-bond modes. A detailed account on such results has been given in Ref. 67. The ultrafast structural fluctuations of aqueous systems set a time window, in which a particular local structure exists and rephasing of the vibrational excitations of different environments is possible. As a result, 2D Raman-THz spectroscopy represents a direct probe of the complex intermolecular dynamics in liquids. A correct theoretical description of the nonlinear Raman-THz response requires polarizable water models, which can be benchmarked by comparison with 2D spectra. [76][77][78][79] Two-dimensional THz-Raman spectroscopy implies two interactions with THz pulses and one interaction with the near-infrared Raman pulse. 69,80 In a three-level system, interaction with the first THz pulse generates a coherence on the 1-2 transition, which the second THz pulse at a delay t 1 transfers to the 2-3 transition. This coherent excitation is read out after a time interval t 2 by a near-infrared Raman pulse on the 3-1 transition. The coherences during the periods t 1 and t 2 allow for generating 2D spectra as a function of t 1 and t 2 or-after a 2D Fourier transform-as a function of the corresponding frequencies ν 1 and ν 2 . Anharmonic couplings of low-frequency modes in liquid CHBr 3 , CCl 4 , and CBr 2 Cl 2 have been studied by 2D THz-Raman spectroscopy. 68,69,80 In part, these results have been re-analyzed in Ref. 81.
Third-order THz-infrared-visible spectroscopy 71 aims at elucidating couplings between low-and high-frequency vibrational modes and involves a THz pulse, a femtosecond midinfrared pulse, and a sub-picosecond near-infrared or visible pulse. The THz and the mid-infrared pulse are both in resonance with vibrational transitions, so that each generates a vibrational coherence. The third interaction with the visible pulse transfers the coherent polarization to the visible spectral range and, thus, generates a coherent emission at the frequencies (ν VIS + ν IR ± ν THz ). There is a delay t between the THz pulse and the mid-infrared and visible pulses, the latter with a fixed temporal delay to each other. The emitted electric field is heterodyned with a local oscillator for sensitive phase-resolved detection. From this signal, 2D spectra are calculated as a function of ν 1 , which is in the THz range and derived from a Fourier transform along t, and as a function of the mid-infrared frequency ν 2 , which is obtained by subtracting ν VIS from the frequency of the nonlinear signal. Details of data processing and an analysis have been given in Ref. 82. THz-infrared-visible spectroscopy has been applied to phonons in solids and to liquid water.
We note that there are a number of 2D experiments in which a coherent THz emission has been recorded after interaction of the sample with pulses at higher frequencies. A recent example is a study of TO phonon coherences in bulk GaAs. 83

IV. SOFT-MODE EXCITATIONS IN POLAR MOLECULAR CRYSTALS
Soft modes are lattice vibrations with a particularly strong coupling to the electronic system of a polar and/or ionic crystal. Excitation of a soft mode is connected with a pronounced relocation of electronic charge, thus directly affecting the electric polarization of the material. Vice versa, the coupling to electronic degrees of freedom enhances the vibrational oscillator strength significantly. Prototypical systems displaying soft modes are displacive ferroelectrics, e.g., with a perovskite lattice structure, and other ionic crystals. The interplay of soft-mode excitations and changes of electronic charge density has been revealed in femtosecond x-ray diffraction experiments following the transient charge densities in space and time. 85,86 Much less is known on the intrinsic nonlinear response of soft modes, 87 in particular in molecular materials. In the following, we discuss results from 2D-THz spectroscopy which give detailed insight in soft mode nonlinearities in aspirin crystals. 84 The crystal structure of aspirin consists of layers, in which pairs of aspirin molecules are arranged as cyclic hydrogenbonded dimers [ Fig. 6(a)]. The layers are connected via hydrogen bonds between the acetyl groups of neighboring molecules. Theoretical calculations of the electronic and vibrational structure of aspirin have shown that a number of vibrational modes in the THz range display a significant coupling to the electronic system. 88 Among them is the rotation of the methyl groups at a frequency of 1.1 THz, much lower than that of a free CH 3 rotator. Figure 6(b) shows the linear THz absorption of a polycrystalline aspirin sample at a temperature of 80 K as measured with weak THz pulses. The methyl rotation gives rise to the weak absorption shoulder around 1 THz. A 2D-THz experiment was performed with two pulses A and B, the spectra of which are shown in Fig. 6(b). The maximum electric field amplitudes of pulses A and B are 25 and 50 kV/cm. The 2D-THz spectrum shown in Fig. 6(c) exhibits different pump-probe and photon-echo signals, among which the A pump -B probe signal at ν τ = 0 is the strongest one (oval boundary). Applying a 2D Fourier filter (see Sec. III B) and transforming the signal back to the time domain gives the contour plot of Fig. 6(d), which shows an oscillation along t with a frequency below 2 THz, the maximum of the pulse spectra in Fig. 6(b). Moreover, there is a phase shift compared to the field of pulse B (not shown). To assess this transient in more detail, it was Fourier-transformed along the real time t to generate the spectrally resolved pump-probe signal as a function of ν t and the delay time τ between the two pulses [ Fig. 6(e) and (f)]. The symbols in panel (e) give the pump-probe signal at a delay time τ = 1.6 ps, which shows a bleaching (increased transmission) for frequencies ν t < 1.4 THz and an enhanced absorption (decreased transmission) at higher frequencies. This behavior persists up to delay times of several picoseconds [ Fig. 6(f)].
Pulse A in the 2D experiment induces both a population change of the v=0 and v=1 states of the vibrations within the pump spectrum and a nonlinear polarization. The population changes result in a third-order (χ (3) ) nonlinear response with a transmission increase on the v=0 to 1 transition and a transmission decrease on the anharmonically redshifted v = 1 to 2 transition. The spectral envelope of the measured pump-probe signal in Figs. 6(e) and (f) with a blueshifted transmission decrease demonstrates that such population-induced signals play a minor role here. Instead, the signal is dominated by a pumpinduced blueshift of the v = 0 to 1 transition of the soft mode and represents a response beyond the χ (3) regime.
The observed behavior reflects the transient local-field response of the aspirin crystal. Due to the strong enhancement of the soft-mode transition dipole by coupling to the electronic system, there is a strong dipole-dipole coupling between different aspirin sites, a concomitant electric polarization and a strong Lorentz electric field. In thermal equilibrium, the electrostatic energy of the crystal is minimized, resulting in a pronounced redshift of the soft-mode frequency compared to its value without a Lorentz field. The transition frequency of 1.1 THz observed in linear THz absorption spectrum is a result of this redshift. Excitation by pulse A in the 2D-THz experiment reduces the contribution of the excited molecules to the overall polarization. This nonlinear saturation of electric polarization reduces the local electric field and shifts the soft mode back to a higher frequency, i.e., a blueshift arises. This nonperturbative mechanism governs the line shapes of the pump-probe signals in Figs. 6(e) and (f) and strongly dominates over the population-induced signals. This interpretation is supported by an analysis of the photon-echo signals in the 2D-THz spectrum [ Fig. 6(c)], which has been presented in Ref. 84. The photon-echo signals display pronounced components at negative delay times τ, which originate from non-instantaneous contributions of Lorentz fields with the time structure of a free-induction decay.
The 2D-THz results reveal a nonlinear optical response of the aspirin soft mode in the nonperturbative regime, even for moderate THz driving fields of some 50 kV/cm. This fact reflects the sensitivity of collective electric properties, here the local field, to a comparably weak external stimulus. A similar behavior is expected for a large range of vibrational and/or phononic excitations with a pronounced coupling to the electronic system. The interplay of vibrational excitations and electronic properties is highly relevant for optically induced phase transitions of crystal structure and for transiently steering macroscopic electric properties of functional materials.

V. FIELD-INDUCED IONIZATION AND ELECTRON TRANSPORT IN LIQUID WATER
Intermolecular and/or collective molecular motions together with librations, i.e., hindered rotations, give rise to the complex absorption spectrum of liquid water in the frequency range from a few GHz to some 30 THz [ Fig. 1(e)]. The broad unstructured envelopes and the substantial spectral overlap of different absorption bands make an analysis of linear spectra in terms of the underlying degrees of freedoms highly demanding. This fact has generated a quest for nonlinear timeresolved spectroscopies, by which the different excitations can be discerned in a better way and couplings between modes be identified.
A first proposal for nonlinear 2D-THz spectroscopy of liquid water has been published in the 1990s. 89 Recent work in 2D Raman-THz and THz-infrared-visible spectroscopy has started to address the properties of intermolecular vibrations (see Sec. III C). There has been related work studying the nonlinear Kerr effect in a similar frequency range. [90][91][92][93][94][95] Until very recently, however, there have been no genuine 2D-THz results on liquid water, partly due to the insufficient experimental sensitivity and comparably small transition dipole moment of intermolecular and collective modes. In the following, we discuss very recent results of a 2D-THz study of liquid water, which applies THz pulses with peak electric field amplitudes in a range from 100 kV/cm to 2 MV/cm to study the nonperturbative response of the liquid. 96 The experimental geometry is sketched in Fig. 7(a). A phase-locked pair of pulses A and B separated by a delay τ interacts with a 50-µm thick water jet. As the water jet cannot be operated in vacuum, it was placed in an atmosphere of dry nitrogen, together with the 2D-THz setup. The transmit- ted pulses are recorded by EO sampling (EOS). A transmitted pulse pair detected this way is shown in Fig. 8(a). The pulses cover a spectral range from 0.2 to 2 THz. The nonlinear signal E NL (τ,t) is given by Eq. (2). In Fig. 8(b), the field E AB (τ,t) is plotted as a function of τ and t, while panel (d) shows the nonlinear signal field E NL (τ,t). In the 2D frequency-domain [ Fig. 8(c)], the strong pump-probe signal generated with the stronger pump pulse A and the weaker probe pulse B (A-B signal) is complemented by the much weaker B-A pumpprobe signal. This observation points to a highly nonlinear dependence of E NL on the pump field, which has been studied in a series of measurements with pump pulses A of different peak field strengths E A . Data for different values of E A and a pump-probe delay of τ = 7 ps are summarized in Fig. 8(e).
Here, the orange lines show the transmitted probe pulse B measured without excitation and the blue lines the nonlinear signal field.
The results of Fig. 8(e) demonstrate a threshold behavior of the nonlinear pump-probe signal. For field amplitudes of pulse A below the threshold E th = 225 kV/cm, the nonlinear signal (blue lines) vanishes for all delay times τ. Above the threshold, E NL (τ = 7ps,t) shows a dispersive shape. It reaches considerable amplitudes up to 100 kV/cm. As seen from the total transmitted probe field (dashed lines), the nonlinear response leads to an increase in amplitude and a phase shift along the t axis, i.e., there is a nonlinear change of both the real and imaginary part of the refractive index of the excited water jet. A quantitative analysis discussed in detail in Ref. 96 gives a decrease of the real and imaginary part by some 10 %, with a flat spectral envelope distinctly different from the absorption spectrum shown in Fig. 1(e).
The strong changes of the real and imaginary parts of the refractive index together with their spectra demonstrate that excitations of intermolecular vibrations giving rise to the comparably weak absorption bands of Fig. 1(e) (molar extinction coefficient ε(1 THz) = 1.67 M −1 cm −1 ) play a minor role for the observed nonlinear response. The same holds for the nonlinear Kerr effect in the THz range which is connected with changes of the refractive index several orders of magnitude smaller than found here. [96][97][98] Instead, the nonlinear response is governed by the generation of solvated electrons in the presence of the THz field, strongly modifiying both the real and imaginary part of the refractive index in a wider frequency range.
Water molecules in the liquid phase have an electric dipole moment d ≈ 2.9 Debye, which gives rise to a very strong electric field in the hydrogen bond network of H 2 O. Due to thermally activated molecular motions at ambient temperature, such fields display fluctuations in a time range from approximately 50 fs to several picoseconds. An electric field trajectory from molecular dynamics simulations is shown in Fig. 7(b), where the electric field projected on the 3a 1 electron orbital of a water molecule is plotted as a function of time. Individual fluctuation peaks reach amplitudes on the order of 100 to 200 MV/cm.
At such high external fields, the binding potential of electrons in the water molecule is strongly distorted and tunneling ionization becomes possible, as schematically shown in Fig. 7(c). As both amplitude and direction of the fluctuating electric field change on an ultrafast time scale, recombination processes of a released electron with its H 2 O + parent ion at the ionization site strongly dominate over successful ionization events, resulting in a very low concentration of released electrons on the order of some 10 −7 M. This scenario changes dramatically if one superimposes a directed THz field on the fluctuating electric field. Now, the THz field can drive the released electron away from the ionization site and, thus, induce a persistent charge separation. After charge separation and transport, electrons are equilibrating in their own solvation shell on a time scale of a few picoseconds [ Fig. 7(d)]. 99 The presence of such solvated electrons changes the dielectric function or refractive index in the THz range and, thus, causes the observed strong nonlinear response. This response is in the regime of nonperturbative light-matter interaction under nonresonant conditions. Independent evidence for the generation of solvated electrons comes from the observation of the characteristic absorption band of solvated electrons around 700 nm after interaction of the water sample with the strong THz field. 96 The strength of this absorption with a known molar extinction coefficient 100 allows for estimating the electron concentration c e . One derives c e = 5 × 10 −6 M for a THz peak field of 500 kV/cm and c e = 2 × 10 −5 M for 1.9 MV/cm.
The results in Fig. 8(e) demonstrate a threshold electric field E th = 225 kV/cm of pulse A for the generation of solvated electrons. This threshold originates from the requirement of a persistent spatial separation of the electron from its parent ion by the transport process along the THz field. The electron needs to take up a ponderomotive (or kinetic) energy exceeding the ionization energy of a water molecule of approximately 11 eV. This energy is provided by the local THz field in the liquid which is, due to the polarity of the water molecules, roughly two times higher than the external field. Thus, an external field of some 250 kV/cm corresponds to a local field of 500 kV/cm at which the ponderomotive energy reaches a value of 11 eV, the water ionization energy. Below this threshold, recombination of the electron with its parent ion prevails, while for even higher fields the electron either travels larger distances in the liquid or generates secondary electrons by impact ionization of water molecules.
The present 2D-THz study of water in the nonperturbative regime reveals a novel aspect of the strong fluctuating electric fields in polar liquids, namely their potential for inducing tunneling ionization processes. In combination with strong THz fields, mobile electrons can be generated and their transport and localization processes be steered by tailoring the THz fields in amplitude and time. Such an approach may allow for studying a new regime of charge transport in liquids.

VI. CONCLUSIONS
Two-dimensional THz spectroscopy has undergone a rapid development and represents an important addition to the portfolio of multidimensional spectroscopies. It allows for studying resonant excitations and field-driven processes up to arbitrary nonlinear orders, i.e., clearly exceeds the third-order limit. In particular, the response of a molecular ensemble under conditions of a nonperturbative light-matter interaction becomes accessible.
There is a number of directions for future developments. First, a more systematic application of 2D-THz spectroscopy and related 2D methods to intermolecular modes in liquids and low-frequency phonons in crystalline and disordered solids appears promising. Secondly, the investigation of charge transport processes in molecular systems driven by nonresonant THz fields holds potential for elucidating frictional forces and the interplay between delocalized and localized states of electrons, protons, and ions. Third, fieldinduced changes of vibrational transition frequencies, i.e., a vibrational Stark effect induced by strong THz fields, will allow for mapping electric fields in liquids and biomolecules at their intrinsic THz fluctuation frequencies and, thus, go beyond steady-state electrostatics. Work along those lines will benefit from the ongoing rapid development of THz sources and detection systems, and the implementation of more complex pulse sequences in multidimensional THz spectroscopy.

ACKNOWLEDGMENTS
We thank our colleagues who were involved in the research reviewed here, in particular Klaus Biermann, Benjamin P. Fingerhut, Giulia Folpini, Christos Flytzanis, Ahmed Ghalgaoui, and Carmine Somma. We acknowledge funding by the Deutsche Forschungsgemeinschaft (WO 558/14-1) and by the European Research Council (ERC) under the European Unions Horizon 2020 Research and Innovation Program (grant agreement no. 833365).

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