2D THz spectroscopic investigation of ballistic conduction-band electron dynamics in InSb

Using reflective cross-polarized 2D THz time-domain spectroscopy in the range of 1-12 THz, we follow the trajectory of the out-of-equilibrium electron population in the low-bandgap semiconductor InSb. The 2D THz spectra show a set of distinct features at combinations of the plasma-edge and vibration frequencies. Using finite difference time domain simulations combined with a tight binding model of the band structure, we assign these features to electronic nonlinearities and show that the nonlinear response in the first picoseconds is dominated by coherent ballistic motion of the electrons. We demonstrate that this technique can be used to investigate the landscape of the band curvature near the Gamma-point as illustrated by the observation of anisotropy in the (100)-plane.


Introduction
The development of ultrafast laser systems and light conversion techniques has enabled the generation of few-cycle terahertz (THz) pulses with high electric fields and has enabled a new range of nonlinear optical studies. For carrier dynamics in semiconductors in particular, a variety of high-field effects have been demonstrated using ultrafast THz spectroscopy. These include interband tunneling in GaAs [1], impact ionization and intervalley scattering [2][3][4][5], Bloch oscillations in GaSe [6] and ballistic transport of electrons in GaAs and InGaAs [7,8]. For the latter, carriers are freely accelerated by electric fields along the conduction band, leading to a coherent optical response which contains information about the band structure of the material [7,8].
THz multidimensional spectroscopy provides additional spectral information compared to conventional ultrafast THz spectroscopy techniques. By probing the nonlinear response of the system along two independent frequency components, it enables studying the coupling between different degrees of freedom (electrons, lattice, etc.) and also purely electronic phenomena such as the ballistic transport of electrons. In the study of condensed matter systems, 2D THz spectroscopy has so far been used in the time domain to investigate different types of excitations and their couplings, such as the correlations of electronic and lattice excitations in GaAs/AlGaAs quantum wells [9], two-phonon coherences in bulk InSb [10] and magnon coherence and correlations in YFeO3 [11]. For most of these studies relatively narrowband THz radiation is tuned to be close to resonance with intrinsic material frequencies (intersubband resonances or phonon/magnon resonances for example). 2D THz spectroscopy has also been implemented to investigate nonlinear dynamics of electrons in various low-bandgap systems [12][13][14][15]. In many materials there are a wide range of frequencies in the THz range that are relevant for studies of electron dynamics and coupling, making broadband THz pulses combined with 2D THz spectroscopy a highly attractive technical goal.
In this letter, we present cross-polarized 2D THz spectroscopy measurements over a broad 1-12 THz range to investigate the time-resolved electronic band nonlinearities of the low-bandgap semiconductor InSb during the first few picoseconds after excitation by a strong, nearly single-cycle THz field. With this method and precise controls over the field polarization, we can perform a complete parity analysis and thus are able to distinguish nonlinear contributions from physically distinct mechanisms. Although electron scattering effects are known to be sub-100 fs, we show that the coherent ballistic motion of field-driven electrons dominates the nonlinear response in the first few picoseconds, as previously demonstrated in GaAs films [8]. Our experimental observations are supported by finite-difference timedomain (FDTD) simulations [16] of the ballistic response of InSb to the THz electric fields. These enable us to identify the conduction band curvature characteristics that give rise to the observed nonlinearities.

Experimental methods
In order to create the THz pulses used in our measurements, we start with an amplified femtosecond laser system. The output of a Ti:Sapphire amplifier (λ= 800 nm, 100-fs pulse, repetition rate: 1 kHz) seeds two three-stage opticalparametric-amplifiers (OPAs) resulting in two output beams tuned at 1.3 m and 1.5 m, with 1 mJ/pulse and 1.5 mJ/pulse respectively. These are used to generate broadband THz pulses from two separate sources: a two-color plasma source [17] providing up to 100 kV/cm peak field (amplitude) at the sample position with a 1-12 THz bandwidth and a source based on optical rectification (OR) in an organic crystal (DSTMS) [18] providing up to 250 kV/cm peak field with 1.5-4.5 THz spectral content. Two independent choppers on the plasma source and OR source paths are set to allow transmission of the light pulses at a quarter and half of the repetition rate of the laser, respectively, in order to excite the sample with either the "probe" field Eplasma, the "pump" EOR or the combined fields EOR+p. The inset of figure   1a shows a schematic of the temporal pulse sequence.
A key feature of this experiment is the careful control of the two THz beam polarizations to obtain cross-polarized electric fields to excite the sample. The THz beam from plasma generation is horizontally polarized, which is maintained until it interacts with the sample. The THz beam generated from OR is also initially horizontally polarized, but is rotated to vertical polarization using a sequence of 2 wire-grid polarizers (WGs). In this scheme the first polarizer is oriented in order to pass electric fields polarized at an angle θ with respect to the horizontal, and the second is oriented to allow only vertically polarized electric fields to pass. By rotating θ between values of -45° and +45°, we can precisely control both the amplitude and sign of the electric field pulses transmitted by the second WG. The two pulses are combined on the sample at normal incidence using another WG polarizer that transmits the horizontallypolarized plasma-generated pulse and reflects the vertically-polarized OR-generated pulse, as shown in Fig. 1a. The horizontally-polarized plasma-generated pulse reflected by the sample gets transmitted through the combining WG and then is partly reflected on a Si beamsplitter towards the detection cell. There, the THz field is detected using broadband Air-Biased Coherent Detection (ABCD) [19] with a split-off sampling beam at 1.3 m, mostly sensitive to horizontally polarized beams. As the cross-polarizations are well-defined, the part of the OR-generated pulse which is reflected by the sample is mostly reflected back to the generation crystal at the combining WG so that despite the collinear incidence geometry, the detection is nearly insensitive to the OR-generated field.
Examples of the time-and spectral-dependence of the field Eplasma generated by the plasma source (red) and that of the field EOR generated by the OR source (green) are shown in figures 1b and 1c. These traces were measured using the broadband ABCD method at the detection position. To obtain accurate electric field measurements at the sample position, we temporarily replaced the sample by electro-optic sampling detection, using a 100 m-thick GaP crystal.
Due to its limited spectral sensitivity towards higher frequencies, this measure gives a lower estimate of the actual field strength for Eplasma. In the following, the indicated peak field amplitudes refer to these EO-sampling measurements at the sample position. The delays of the sampling infrared beam for ABCD and the pumping infrared beam for the OR source are independently controlled, and in the following will be referred to as "detection delay t" and "excitation delay ", respectively. For the 2D measurements, both delays are scanned and the nonlinear signal ENL, defined as can thus be extracted, for every delay t and . The detection specificity to horizontallypolarized fields is further enhanced by placing an additional WG polarizer before the detection cell.
For 1D reflectivity measurements, the setup is used as a broadband time-domain spectrometer (TDS) with the plasma field only. In order to get a reference signal at the exact position of the sample, an intense 1500-nm-laser pulse illuminates the semiconducting sample, inducing a transient state with a plasma frequency far beyond the observed frequency range. Using this photo-screening method, we briefly create a mirror-like reference at the exact position of the sample. The sample equilibrium reflectivity can be determined by comparing the reference with the semiconducting state response measured with the following laser pulse 1 ms later [20]. Fig. 1. a) Measurement geometry of reflective 2D THz spectroscopy on InSb at normal incidence. The horizontally-polarized plasma field, Eplasma, is transmitted back and forth through the wire-grid polarizer (WG), while the vertically-polarized field from OR in the organic crystal, EOR, is mostly reflected. The nonlinear response is retrieved using the chopping scheme depicted in the inset. b) Eplasma (red) and EOR (green) normalized electric fields, measured with broadband ABCD method. c) Corresponding Eplasma (red) and EOR (green) normalized spectra.

Figure 2a
shows the reflectivity amplitude (black circles) and phase (green circles) measured on the (100)-cut bulk nominally-undoped InSb sample by TDS using a strongly attenuated broadband plasma source and the photo-screening method to reference the data. The reflectivity data can be fitted using the combined Drude-Lorentz model for the dielectric function, given as: where TO and LO are transverse (TO) and longitudinal optical (LO) phonon frequencies respectively, ph is the corresponding phonon damping, the plasma frequency of Drude carriers with damping p, and ∞ represents the high frequency dielectric function. Fitting the complex reflectivity, we extract all parameters: the plasma-edge frequency , TO and LO-phonon frequencies TO = 5.3 THz and LO = 5.7 THz, ∞ = 15.8 and the damping p = 0.3 THz and ph = 0.6 THz. All are in good agreement with previously reported data [21], despite common variations in growth-process that affects doping and thus the Drude coefficients. Figure 2b shows the 2D temporal trace of the nonlinear signal acquired on InSb varying detection delay t and excitation delay  with peak electric fields of Eplasma  25 kV/cm and EOR  65 kV/cm. The field Eplasma has a constant phase at a given detection delay t, while EOR evolves following the diagonal with a delay (dashed line). For delays , the nonlinear signal shows a vertical structure almost completely independent of the excitation delay, which can be understood as a persistent change of the reflectivity of Eplasma induced by EOR. For InSb, these changes can be attributed to carrier multiplication induced by impact ionization [3]. In the first picosecond, where both THz pulses temporally overlap, the nonlinear signal shows a clear dependence on both detection delay t and excitation delay . This nonlinear contribution amounts to about 15% of the incoming peak field Eplasma.
The different signal modulations observed in the time domain transform into distinct peaks in frequency space, as shown by the 2D spectrum presented in figure 2c. We can assign these features to combinations of plasma-edge frequency p and LO phonon frequency LO that were measured with TDS ( fig. 2a). We note here that the peaks do not appear exactly at the plasma-edge frequency previously measured (2.2 THz), but rather at the maximum of the phase shift of the reflectivity (1.9 THz). Expressed as (t, ), the main observed peaks here are (p, 0), (p, 2p), (LO, 0) and In order to distinguish the various contributions to the nonlinear response in the 2D frequency domain and to separate coherent and incoherent contributions, we can invert the polarity of EOR (phase shift of 180°) and then extract the odd (E -) and the even (E + ) components of the nonlinear signal. A similar approach was followed for 1D measurements by Ho et al. [22]. For 2D experiments, odd and even components are a combination of 5 different fields, defined as: where OR indicates the inverted polarity of the OR field.
The even component E + contains nonlinear terms with even orders of the OR field, including incoherent contributions from carrier multiplication proportional to the field intensity. The odd signal, conversely, contains nonlinear terms with odd orders of the OR field and is then expected to show only coherent responses that have a fixed phase relation to EOR. Because our detection scheme does not completely eliminate contributions from EOR, we might expect some contributions in the "odd" signal arising from an influence of Eplasma on EOR. These contributions are, however, very close to the noise floor of our measurement.

Simulations
To better understand the features of the 2D spectra observed experimentally, we model the response of the InSb bulk sample in the temporal overlap region beyond the perturbation limit. Here we model the electronic motion in the first few picoseconds as a coherent wave packet motion of conduction band electrons [8,23]. Indeed, these systems have a long electronic coherence time [8] and the probability of tunneling from the valence to the conduction band is reduced by the low joint density of states at the -point. In the following, we will hence neglect electronic decoherence and approximate the dynamics as those of carriers moving in a two-dimensional conduction band. Our model therefore describes the ballistic motion of electrons, driven by THz electric fields in the conduction band but does not include interband tunneling or impact ionization or any scattering channels apart from those in the Drude-Lorentz model.
To implement the physics of electronic motion and the interaction of polar phonons with few-cycle THz pulses, we solve the time dependent Maxwell's equations in the presence of an interacting medium, using the FDTD method [16,23]. This algorithm, together with a description of ballistic motion of conduction band carriers for a polar semiconductor in response to electric fields, enables modeling of the response in the temporal overlap region beyond the perturbation limit. To allow for arbitrary crystal orientations and two dimensional carrier motion, we include the full band structure, obtained from tight binding calculations [24], in our FDTD simulations. In the absence of a magnetic field the spin-split conduction bands are equally populated and we therefore approximate the conduction band as an average of the spin-split bands. Following the work by Yu et al. [23], the local electric displacement in the semiconductor is given by: where is the driving THz electric field, ℎ and correspond to the phonon and electron polarization densities respectively and the dielectric background is given by the scalar ∞ . The phonon polarization density ℎ is implemented as a Lorentz oscillator with given strength, damping and frequencies, previously introduced in equation (1). The temporal evolution of the electronic polarization density is then given by: Although a change in scattering rate is known to occur for higher energy electrons, we keep the damping p at its equilibrium value and do the same for the carrier density ne. The equilibrium values, listed above, were found by fitting the measured linear reflectivity shown in figure 2a, so that the simulation has no free parameters.
Using the FDTD method with the previously described assumptions, we model the response of bulk (100)

Discussion
The nonlinear peaks observed in the 2D spectrum in fig. 2c  due to small changes in the plasma frequency that can be attributed to two distinct effects taking place in the InSb system under high THz field excitation. The first effect is impact ionization, the electron-hole pair creation induced by energetic carriers scattering, which is known to occur in InSb under high THz electric field excitation [3]. Impact ionization induces an increase in the carrier density ne, which leads to a change in plasma frequency. The second effect stems from non-parabolic features in the conduction band. When the conduction band electrons are driven by the THz electric fields under ballistic transport conditions, they dynamically explore the conduction band. If the band is nonparabolic, electrons will experience different effective masses for different k positions and therefore changes in plasma frequency. We thus expect that both impact ionization and ballistic transport contribute to the nonlinear peak (p, 0).
Large-amplitude driving of electrons in a non-parabolic conduction band can also lead to nonlinear effects at harmonics of the plasma frequency, which result in the (p, 2p) and (LO, 2p) peaks in the 2D spectrum. We remark here that second order nonlinearities, which would appear at (0,p) and (2p, p) (corresponding to difference and sum frequency generation respectively, observed in [11]), are absent from the spectrum. This is expected given that the only second order nonlinear susceptibility tensor components of the zincblende structure of InSb that are non-zero are those with i ≠ j ≠ k. The second order polarization therefore vanishes at normal incidence on a (100)-surface and the observed peaks can be assigned to third or higher order nonlinearities.
We now focus on the nonlinearities observed at the LO phonon frequency LO, appearing as the (LO, 0) and (LO, 2p) peaks in both the experimental and simulated 2D spectrum. A modification of the dielectric function due to a change in the plasma frequency leads to a redistribution in spectral weight which affects the nearby modes such as the LO-phonon mode, even in the absence of direct coupling. Then a THz field driven modification of the plasma frequency, through either carrier density or effective mass changes, leads to a change in phonon frequency [25] and then to the observed (LO, 0) and (LO, 2p) peaks.
Impact ionization is an incoherent process due to the photo-excitation of electrons and evolves independently of the EOR phase. Ballistic transport, on the contrary, is a coherent process since the electrons are directly driven by both electric fields. These two effects can be partially differentiated by separating odd and even parity signals. In the even signal, the strong feature at (p, 0) contains incoherent changes due to impact ionization as well as a coherent contribution from the pump-probe effect of EOR on Eplasma. In figure 3, the peaks at (LO, 0) and (LO, 2p), previously observed in figure 2c, are weak and so they are not discussed further. In contrast, the odd parity signal of fully crosspolarized beams contains only coherent process contributions. The peak (p, p) observed in the odd signal of fig. 3 for  = 22.5° and  = 67.5° emerges from a coherent process and is then a nonlinear signature of the non-parabolic (2) ijk  conduction band experienced by the ballistically-driven electrons. Its intensity dependence with different  angles (it disappears for  = 45°) leads us to conclude that this feature originates from the local band curvature experienced by the driven electrons for different directions within the conduction band. The weaker feature at frequency (p, 0) can be attributed to leakage from intense feature at (p, 0) in the even signal, due to non-perfect balancing of the polarity inversion.
The simulated spectra, from the last sets of fig. 3a and 3b, are in good agreement with the experimental results both in terms of intensity and angular dependence. In the simulated even signal, the (p, 0) amplitude is close to 50% of the maximum of Eplasma spectrum, which is noticeably higher than the 30% found experimentally. The three secondary peaks observed experimentally, at (p, 2p), (LO, 0) and (LO, 2p) (see fig. 2c) are well reproduced in the simulated even signal but again with higher relative intensities than in the experiment. These discrepancies in intensity could be possibly due to impact ionization, which occurs in the experiments but is not included in the simulations. In addition, the simulation is based on a plane-wave approximation with homogeneous field-strength, while in the experiment, the sample is placed at the focus of the THz beam, where the peak field strength is only reached at the very center. In the simulated odd spectrum, the relative amplitude of the (p, p) peak is near 2%, comparable with the experiment. The simulated response for an angle =45° (spectrum not shown here) shows no feature in the odd signal, also in good agreement with experimental observations.
The agreement between experimental data and simulations that include only ballistic transport indicate that ballistic transport dominates the THz nonlinear response in the first few picoseconds after excitation. This agreement is particularly good for the odd spectra where we were able to successfully isolate coherent ballistic transport signatures.
We can then conclude that the odd spectral feature (p, p) is a manifestation of the band curvature symmetry along kOR, and orthogonal band curvatures i.e. along kplasma. These differences in band curvature lead to different nonlinear responses. Figure 4b shows the evolution of the (p, p) peak intensity from the odd signal with  angles between 0° and 90°. Results from FDTD simulations (blue line) are in good agreement with experimental values (red markers) extracted from the spectra in figures 3a and 3b. In particular, they show a strong dependence on the  angle, with maxima at 22.5° and 67.5° and minima at 0°, 45° and 90°. At this stage, it seems evident that the nonlinear signals in the simulations arise from a deviation of the conduction band dispersion from the parabolic free-energy model: with . This means that our technique can be used to investigate anharmonicity and anisotropy characteristics of the band curvature. In order to do so, we simulate the response for two different anharmonic conduction band models given by: First, we will consider the fully parabolic band model given by eq. (2). Here our simulations show no nonlinearities.
This is expected since the effective mass of a parabolic conduction band is constant for any k-vector (we recall that ). For the model Eband2 with anharmonicity, we observe nonlinearities but only in the even channel. These are similar to features in the even spectrum previously simulated with the actual InSb band structure Although these analytical models are too simple to accurately characterize the conduction band of InSb, they provide the right band curvature landscape within the k-range of the ballistic dynamics and explain qualitatively the experimentally measured nonlinearities. Specifically, we showed that the nonlinear signal of even parity is sensitive to the deviation from band harmonicity and the odd signal captures the anisotropy of the band structure. A more numerically efficient implementation of the FDTD calculation including higher order analytical models would enable approaching the conduction band curvature properties by fitting to the experimental data directly.

Conclusion
In conclusion, we studied the nonlinear dynamics of electrons in InSb by cross-polarized 2D THz spectroscopy  spectra and their symmetries contain information regarding the local band curvature so that we can identify the symmetry of the valley in which the carriers reside. Our technique could in the future be applied to study a broad range of systems including the band structure of strongly nonparabolic electronic systems such as Dirac and Weyl semimetals in bulk, providing a possible complement to more direct measures of band structure that are typically, however, limited to near-surface states.