1. 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–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 YFeO₃ [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–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 time-domain (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.

2. 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 optical-parametric-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 E_plasma, the “pump” E_OR or the combined fields E_{OR + p}. The inset of Fig. 1(a) shows a schematic of the temporal pulse sequence.

 

Fig. 1 a) Measurement geometry of reflective 2D THz spectroscopy on InSb at normal incidence. The horizontally-polarized plasma field, E_plasma, is transmitted back and forth through the wire-grid polarizer (WG), while the vertically-polarized field from OR in the organic crystal, E_OR, is mostly reflected. The nonlinear response is retrieved using the chopping scheme depicted in the inset. b) E_plasma (red) and E_OR (green) normalized electric fields, measured with broadband ABCD method. c) Corresponding E_plasma (red) and E_OR (green) normalized spectra.

Download Full Size | PDF

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 horizontally-polarized plasma-generated pulse and reflects the vertically-polarized OR-generated pulse, as shown in Fig. 1(a). 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 E_plasma generated by the plasma source (red) and that of the field E_OR generated by the OR source (green) are shown in Figs. 1(b) and 1(c). 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 E_plasma. 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 E_NL is defined for every delay t and τ as $E_{NL}=E{\text{'}}_{OR+p}-E{\text{'}}_{plasma}-E{\text{'}}_{OR}$, where E’_plasma, E’_OR and E’_{OR + p} are respectively the plasma-generated, the OR-generated and the combined fields detected after interaction with the sample. The detection specificity to horizontally-polarized fields is enhanced by placing an additional WG polarizer before the detection cell, so that OR-generated field leakage is kept at the minimum in E’_{OR + p}.

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].

3. Results

Figure 2(a) 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:

 

Fig. 2 a) Reflectivity amplitude (black circles) and phase (green circles) of InSb measured by TDS. Solid lines show fits to the complex reflectivity data, based on Eq. (1). Dashed red lines indicate the maxima in the phase-shift of plasma-edge and LO phonon frequencies at 1.9 THz and 5.7 THz, respectively. b) Nonlinear response E_NL(t,τ) of InSb measured with polarizations E_plasma and E_OR along (010) and (001), respectively, as a function of detection delay t and excitation delay τ, for E_plasma ∼25 kV/cm and E_OR ∼65 kV/cm. The color scale is normalized to the maximum of E’_plasma 2D measurement. c) 2D spectrum obtained from the Fourier transform of E_NL(t,τ), plotted as a function of detection frequency ν_t and excitation frequency ν_τ. Values are given with respect to the maximum of E’_plasma 2D spectrum.

Download Full Size | PDF

Figure 2(b) shows the 2D temporal trace of the nonlinear signal acquired on InSb varying detection delay t and excitation delay τ with peak electric fields of E_plasma ∼25 kV/cm and E_OR ∼65 kV/cm. The field E_plasma has a constant phase at a given detection delay t, while E_OR evolves following the diagonal with a delay t + τ = 0 (dashed line). For delays t + τ > 1.5 ps, 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 E_plasma induced by E_OR. For InSb, these changes can be attributed to carrier generation induced by impact ionization [3] and interband tunneling [1]. 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 linear response E’_plasma.

The different signal modulations observed in the time domain transform into distinct peaks in frequency space, as shown by the 2D spectrum presented in Fig. 2(c). We can assign these features to combinations of plasma-edge frequency ν_p and LO phonon frequency ν_LO that were measured with TDS (Fig. 2(a)). 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 (ν_LO, 2ν_p).

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 E_OR (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: $E^{\pm }=\left(E{\text{'}}_{OR+p}-E{\text{'}}_{plasma}-E{\text{'}}_{OR}\right)\pm \left(E{\text{'}}_{\overline{OR}+p}-E{\text{'}}_{plasma}-E{\text{'}}_{\overline{OR}}\right)$ where$\overline{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 generation 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 E_OR. Because our detection scheme does not completely eliminate contributions from E_OR, we might expect some contributions in the “odd” signal arising from an influence of E_plasma on E_OR. These contributions are, however, very close to the noise floor of our measurement.

In order to investigate any possible anisotropy in the nonlinear response of the sample, we studied the dependence of the 2D spectrum on the angle between the THz electric field and a reference axis of the InSb sample. The polarizations of E_plasma and E_OR are kept orthogonal, as described in Fig. 1(a), but the sample is rotated so that the THz electric fields excite different directions. We define the angle φ as the angle between the E_plasma polarization and the in-plane InSb axis (010) which in reciprocal space corresponds to the angle between E_plasma and the Γ→X direction in the Brillouin zone.

Figures 3(a) and 3(b) show odd and even 2D spectra for two sets of measurements. Two different φ angles were acquired per set, φ = 45° and φ = 22.5° for Fig. 3(a) and φ = 45° and φ = 67.5° for Fig. 3(b). In both cases, we see that odd and even spectra present different features. The even signal E⁺ is similar for all angles and presents the main characteristics previously pointed out in Fig. 2(c), with the strongest peak at (ν_p, 0) (relative intensity of 30%). The odd signal E^- shows a weak but well-defined feature on the diagonal at (ν_p, ν_p) (relative intensity near 2%). This feature appears for φ angles of 22.5° and 67.5° (on Figs. 3(a) and 3(b) respectively) but is at the noise level for φ = 45° (in both figures) and thus indicates an anisotropic response.

 

Fig. 3 a) Odd and even parity spectra of InSb for E_plasma ∼35 kV/cm and E_OR ∼65 kV/cm. The first two sets show the experimental results from the same measurement at two different angles, φ = 45° and φ = 22.5°. The last set shows simulated spectra (FDTD) for an angle φ = 22.5°. b) Odd and even parity spectra for similar electric fields. First two sets: experimental results for the angles φ = 45° and φ = 67.5°; last set: simulated spectra at φ = 67.5°. Intensities are given with respect to the maximum of the E’_plasma spectrum. Simulated spectra at φ = 45° (not shown here) lead to a vanishing odd signal.

Download Full Size | PDF

4. 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:

Since the electrons explore only a few percent of the Brillouin zone, modifying only slightly the plasma frequency, the scattering rate is assumed to be energy independent and the damping γ_p is kept at its equilibrium value. As we do not account for carrier generation in this model, the carrier density n_e is also kept constant. The equilibrium values, listed above, were found by fitting the measured linear reflectivity shown in Fig. 2(a), so that the simulation has no free parameters.

Using the FDTD method with the previously described assumptions, we model the response of bulk (100)-InSb excited by two THz electric fields. The THz electric field waveforms used in the simulations come from a fit to the experimental plasma and OR temporal traces and we set E_plasma = 35 kV/cm and E_OR = 65 kV/cm. The simulated linear response to the plasma field excitation alone (not shown) is in very good agreement with the measured reflectivity spectrum presented in Fig. 2(a). Under the excitation of both E_OR and E_plasma, the nonlinear signal obtained in the temporal overlap is about 15% of the linear response E’_plasma, comparable to the experimental nonlinear signal presented in Fig. 2(b).

By Fourier transforming the temporal responses, we obtain simulated odd and even spectra for different configurations. The simulated spectra for φ = 22.5° and for φ = 67.5° are included in Fig. 3 for comparison with the experimental data. As expected from the fourfold symmetry of the potential in this plane, both angles show very similar results in the simulations. The even spectrum shows the features already pointed out in Fig. 2(c), at (ν_p, 0), (ν_p, 2ν_p), (ν_LO, 0) and (ν_LO, 2ν_p). The main peak (ν_p, 0) intensity is close to 50% of the E’_plasma spectrum maximum. The odd spectrum shows mainly the peak (ν_p, ν_p), with a relative intensity near 2%. For φ = 45° (not shown), the simulated spectra lead to a similar even signal and to a vanishing odd signal.

5. Discussion

The nonlinear peaks observed in the 2D spectrum in Fig. 2(c) are combinations of the plasma frequency and the phonon frequency appearing at (ν_p, 0), (ν_p, 2ν_p), (ν_LO, 0) and (ν_LO, 2ν_p). Nonlinearities appearing at the plasma frequency ${\Omega }_p=\sqrt{e^2{n}_e/{\epsilon }_0{\epsilon }_{\infty }{m}^{*}}$ (with n_e the electron density and m* the effective mass), especially the main peak (ν_p, 0), are 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 high field induced carrier generation, due to impact ionization and interband tunneling. The induced increase in the carrier density n_e 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 non-parabolic, electrons will experience different effective masses for different k positions and therefore changes in plasma frequency. We thus expect that both carrier generation processes 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${\chi }_{ijk}^{(2)}$ 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.

Carrier generation is an incoherent process due to the photo-excitation of electrons and evolves independently of the E_OR 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 carrier generation processes as well as a coherent contribution from the pump-probe effect of E_OR on E_plasma. In Fig. 3, the peaks at (ν_LO, 0) and (ν_LO, 2ν_p), previously observed in Fig. 2(c), are weak and so they are not discussed further. In contrast, the odd parity signal of fully cross-polarized 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 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 Figs. 3(a) and 3(b), 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 the E’_plasma 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. 2(c)) 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 carrier generation processes, which occur in the experiments but are 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 experienced by the electrons, which are driven in different directions in the conduction band by the orthogonal THz electric fields.

We now further examine the local properties of the conduction band (band curvature) explored by the driven electrons. With FDTD calculations, one can follow step by step in the temporal domain the ballistic motion of electrons driven along the conduction band and the resulting nonlinearities. In Fig. 4(a), we show the reconstructed electron trajectory along the conduction band from FDTD simulations evaluated at the sample surface, for E_plasma = 35 kV/cm, E_OR = 65 kV/cm and φ = 0°. This plot shows the fraction of the conduction band explored by the excited electrons. For these electric field amplitudes, the electrons span up to 5% of the Brillouin zone in k_x and k_y, mostly along the direction of E_OR. Depending on φ and on the polarity of E_OR, the electrons may encounter different lateral band curvatures, i.e. along k_OR, and orthogonal band curvatures i.e. along k_plasma. These differences in band curvature lead to different nonlinear responses. Figure 4(b) 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 Figs. 3(a) and 3(b). 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°.

 

Fig. 4 a) Reconstructed electron trajectory along the conduction band for E_plasma = 35kV/cm, E_OR = 65kV/cm and φ = 0°. Inset: φ is defined as the angle between E_plasma and Γ→X, E_plasma and E_OR are kept orthogonal. b) Maximum of the (ν_p, ν_p) peak intensity as a function of φ angle, FDTD simulations using the InSb band structure calculated with a tight binding model are plotted as a blue line, experimental results for φ = 22.5°, φ = 45° and φ = 67.5° are indicated as red markers. c) Simulated odd and even spectra using analytical model E_band3 for the conduction band (cf. text).

Download Full Size | PDF

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:

6. Conclusion

In conclusion, we studied the nonlinear dynamics of electrons in InSb by cross-polarized 2D THz spectroscopy highlighting the contribution of coherent ballistic transport in the first few picoseconds. The broadband nature of the plasma electric field enabled the observation of nonlinearities at combinations of plasma and phonon frequencies. With the support of FDTD-simulations, we described the temporal overlap regime and demonstrated that the nonlinear observations result from ballistic transport of electrons along an anharmonic and anisotropic conduction band. The 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.