Field-resolved high-order sub-cycle nonlinearities in a terahertz semiconductor laser

The exploitation of ultrafast electron dynamics in quantum cascade lasers (QCLs) holds enormous potential for intense, compact mode-locked terahertz (THz) sources, squeezed THz light, frequency mixers, and comb-based metrology systems. Yet the important sub-cycle dynamics have been notoriously difficult to access in operational THz QCLs. Here, we employ high-field THz pulses to perform the first ultrafast two-dimensional spectroscopy of a free-running THz QCL. Strong incoherent and coherent nonlinearities up to eight-wave mixing are detected below and above the laser threshold. These data not only reveal extremely short gain recovery times of 2 ps at the laser threshold, they also reflect the nonlinear polarization dynamics of the QCL laser transition for the first time, where we quantify the corresponding dephasing times between 0.9 and 1.5 ps with increasing bias currents. A density-matrix approach reproducing the emergence of all nonlinearities and their ultrafast evolution, simultaneously, allows us to map the coherently induced trajectory of the Bloch vector. The observed high-order multi-wave mixing nonlinearities benefit from resonant enhancement in the absence of absorption losses and bear potential for a number of future applications, ranging from efficient intracavity frequency conversion, mode proliferation to passive mode locking. Two-dimensional terahertz strong-field spectroscopy reveals wave-mixing processes up to eighth order in a free-running quantum cascade laser, unraveling its sub-cycle gain dynamics and nonlinearities in a regime of negative absorption. Two-dimensional terahertz strong-field spectroscopy reveals wave-mixing processes up to eighth order in a free-running quantum cascade laser, unraveling its sub-cycle gain dynamics and nonlinearities in a regime of negative absorption. Two-dimensional terahertz strong-field spectroscopy reveals wave-mixing processes up to eighth order in a free-running quantum cascade laser, unraveling its sub-cycle gain dynamics and nonlinearities in a regime of negative absorption.

In order to meet the growing demands of nonlinear optics and ultrashort pulse generation and to explore novel application fields, a detailed understanding and control of coherent and incoherent electron dynamics in operational THz QCLs is indispensable. In seminal pumpprobe experiments [32][33][34][35][36] , gain relaxation oscillations and recovery times ranging from 0.2 ps to a few 10s ps have been reported for mid-infrared (MIR) QCLs. Also, fourwave mixing has been observed in time-integrated experiments 37 . In THz pump-probe measurements using weak pulses, incoherent gain dynamics of THz QCLs [38][39][40] have been studied, and gain recovery times varying from 5 to 50 ps have been found. A recent two-dimensional (2D) spectroscopy approach 41 based on injection seeding and RF switching (i.e., not free running) has yielded a gain recovery time of~10 ps as well as the emergence of different nonlinear optical processes. Finally, 2D spectroscopy using strong THz pulses 42 has allowed the investigation of incoherent absorption dynamics in an unbiased quantum cascade structure without a resonator.
Here we combine strong THz fields with a QCL, and show the first field-resolved 2D spectroscopy of a freerunning THz QCL. Phase-stable, single-cycle THz transients with kV cm −1 peak fields strongly saturate the QCL gain by stimulated emission on ultrashort time scales, allowing for the direct observation of the gain recovery dynamics over the entire current operating range with a sub-cycle resolution by a subsequent THz pulse. Remarkably, we record coherent nonlinearities up to eight-wave mixing, which qualifies the QCL as an efficient nonlinear optical medium and enables a direct extraction of the dephasing time as a function of the bias current. This situation represents a rare scenario where resonantly enhanced high-order nonlinearities occur in a system with population inversion, bypassing the otherwise inevitable absorption losses. A density-matrix theory reproduces both incoherent and coherent multi-wave mixing signals and maps the internal QCL dynamics to the trajectory of the corresponding Bloch vector 43 . This approach may be particularly helpful in tailoring lasers for efficient intracavity frequency conversion, mode proliferation, and self-starting mode locking.

Results
Our QCL is based on an Al 0.1 Ga 0.9 As/GaAs heterostructure (described in detail elsewhere 44 ) that shows laser action at 2.2 THz above a threshold bias current of I th = 910 mA. The laser structure was processed, by wetchemical etching, into a laser ridge, 250 µm wide and with a thickness of 200 µm, which includes both the active region and the substrate (Fig. 1a). Contact electrodes were  Fig. 1 High-field spectroscopy of a THz QCL. a Experimental arrangement; gray: GaAs substrate, red: active medium, gold: waveguide and electrical contacts. The THz waveform E in (blue) is focused onto the QCL facet and the transmitted waveform (E out , orange) is detected electrooptically. The rectangular modulation of the bias voltage allows the current-induced change of the nonlinear response to be extracted. b Measured THz waveform as a function of delay time t, after the transmission through the QCL with (red) and without (blue) bias current. c Amplitude spectrum E ν ð Þ of the waveforms in b and lasing spectrumÃ ν ð Þ of the QCL (gray). d Modulation of the transmission through the QCL as a function of frequency, defined as the ratio of the transmission spectra of the biased and the unbiased QCL (see c) for two peak electric fields E p = E 0 = 1.4 kV cm −1 (red) and E p = E 0 /2 (green, broken curve) deposited on the top and the sides of the laser ridge for current injection. The sample was cleaved into a laser resonator of length 2.9 mm, soldered onto a gold-coated copper base, and finally mounted on the cold finger of a liquid-helium cryostat. The bias voltage provided by an external power supply was modulated at a frequency of 500 Hz and synchronized with the THz probe pulses, allowing for a direct extraction of the bias-induced change in transmission. Single-cycle THz pulses with a tunable peak field strength of up to 3 kV cm −1 and a repetition rate of 1 MHz generated by tilted-pulse-front optical rectification of intense near-infrared (NIR) laser pulses in lithium niobate (see "Materials and methods" section) are polarized parallel to the growth direction of the heterostructure and coupled into the QCL waveguide through the front facet (Fig. 1a). The THz transient propagates through the active medium and is partially transmitted through the end facet of the QCL. We fully resolve the amplitude and the absolute phase of the transmitted waveform by electro-optic sampling in a GaP crystal, as a function of the delay time t. This field-resolved detection allows us to distinguish unequivocally the transmitted THz fields, which have a fixed phase relation with the incident pulses, from the steady-state THz radiation emitted by the free-running QCL. This way, we can study nonlinearities both below and above the laser threshold with excellent contrast despite strong QCL emission. Figure 1b displays the transmitted THz waveform after a single pass through the structure for the biased (red) and the unbiased QCL (blue). The field transmitted through the biased laser is~30% higher and shows a clear oscillation for delay times t > 1 ps owing to the amplification of the incident THz pulse by stimulated emission in the inverted active medium. The Fourier transform of the two waveforms ( Fig. 1c) provides the spectral characteristics of the QCL gain. At a frequency of 2.2 THz, the spectral amplitude of the biased device,Ẽ b (ν), features an increase over the unbiased situation,Ẽ u (ν). This frequency agrees very well with the designed lasing spectrum of the QCL (Fig. 1c). The overall redistribution of charge carriers inside the heterostructure upon biasing leads to a broadband background of increased transmission 34 .
For a more quantitative analysis of the gain, Fig. 1d shows the ratio of transmission,Ẽ b (ν)/Ẽ u (ν) (red curve), measured with an incident THz peak field of 1.4 kV cm −1 (see Supplementary Information for the corresponding phase spectrum). The dominant gain maximum in the QCL active medium at a frequency of 2.2 THz reaches values of up to 37% above the baseline. The latter exceeds unity because of the overall enhancement of transmission by electrical pumping. Assuming a spatially independent linear gain over the cavity length, we obtain a gain coefficient of α E 0 % 2:2 cm À1 . Interestingly, the gain spectrum depends sensitively on the incident THz field strength.
A transmission spectrum measured at half the peak incident electric field of the THz pulse shows a similar base level (Fig. 1d, green curve), yet the resonant gain (α E 0 =2 % 3:4 cm À1 ) reaches values of up to 64% above the baseline-almost twice the value of the high-field case. This drastic change attests to a strong saturation of the QCL gain 39 . Moreover, the linewidth increases from 0.35 to 0.39 THz as the incident THz field doubles, which we attribute to power broadening and a reduced lifetime of the upper laser state under intense excitation. Inhomogeneous broadening effects play only a minor role, since the relative variation of the well thickness of order 1% 45 corresponds to a broadening of the laser transition of less than~40 GHz. In this strongly nonlinear regime, the overall response may result from a variety of simultaneous competing or cooperating nonlinear optical elementary processes. Yet, their individual contributions, which could provide key insights into the microscopic electron dynamics in principle, remain hidden in single-pulse transmission studies.
Field-resolved two-dimensional THz spectroscopy, in contrast, represents a powerful strategy to disentangle all relevant coherent and incoherent nonlinear processes 12,14,46,47 (see "Materials and methods" section). To this end, two identical THz waveforms E A (t) and E B (t, τ), which are mutually delayed by a variable time τ, are focused onto the QCL (Fig. 2a). The first pulse induces a polarization or population change inside the active medium, which is subsequently interrogated by the second pulse. The total transmitted field E AB (t, τ) is recorded with absolute phase and amplitude using electro-optic sampling and consists of a linear superposition of the incident fields and the nonlinear contributions, E NL (t, τ). Synchronous mechanical chopping of the two incident THz beams E A (t) and E B (t, τ) isolates their individual contributions, and by subtracting them from the transmitted signal E AB (t, τ), we are left with the nonlinear signal Figure 2b shows the resulting 2D field map E NL (t, τ) for the unbiased QCL as a function of the delay time t of the electro-optic sampling process and the delay time between the two pulses, τ. The strongest nonlinearities occur near τ = 0 ps, where both pulses overlap temporally. By varying τ, the total incident field is strongly modulated leading to a pronounced variation of E NL . Even for delay times |τ| > 2 ps and t > 0 ps, where the two THz pulses no longer overlap in time, nonlinear optical signals are observed. On the electro-optic time axis, t, they set in with the second pulse, which is pulse A or pulse B (Fig. 2b, broken lines) depending on τ. The nonlinear response manifests as oscillations at the transition frequency of the unbiased QCL at ν I¼0 = 1.7 THz ( Fig. 2d and Supplementary Information). In sharp contrast, passing a bias current of 840 mA, slightly below the laser threshold, through the QCL leads to a qualitatively different response. First, E NL increases by more than one order of magnitude, up to almost 10% of the amplitude of the transmitted THz pulses, E AB (Fig. 2c). Second, the coherent modulation following the second pulse is more long-lived and blue shifted to the laser resonance of ν L = 2.2 THz. Third, the modulation of E NL along t, observed only near time zero for the unbiased structure, now persists for much larger delays |τ|, signifying the presence of coherent nonlinear processes.
To first order, the observed nonlinear signal may be understood considering the saturation of the active QCL gain discussed above. The first THz pulse to arrive at the structure depletes the population of the upper laser level by stimulated emission, decreasing the transmission of the second pulse through the active medium within the gain recovery time. As |τ| increases, the gain is restored by the pump current. Moreover, the first THz pulse imprints its phase on the oscillating intersubband polarization. The relative phase of Δφ = τν L between the induced coherent polarization and the second pulse determines whether the latter can further reduce or restore the population inversion by stimulated emission or absorption.
To expand E NL systematically into its constituent nonlinear processes, we perform a 2D Fourier transform of the time-domain data of Fig. 2c. By comparing the resulting amplitude spectrumẼ NL (Fig. 2e) with the spectrum of the unbiased case (Fig. 2d), we observe that the nonlinearities do not only drastically increase, but also higher-order multi-wave mixing signals appear. The distinct maxima ofẼ NL are located at (ν t , ν τ ) = (ν L , n × ν L ), where n is an integer number. Given the time order of the 2D scans, the linear polarization of the incident waves A and B resonant with the laser transition would be located at (ν L , 0), and (ν L , −ν L ), respectively. Starting from the origin (0,0), every maximum ofẼ NL can be navigated by adding or subtracting a suitable number of these two vectors representing photons from pulses A and B (see example path in Fig. 2e). This Liouville path analysis allows us to assign uniquely the number and the phase of the interacting photons contributing to the nonlinear optical process that gives rise to a specific maximum 12 .
The dominant structures at (ν L , 0) and (ν L , −ν L ) are incoherent pump-probe (PP 1 and PP 2 ) contributions 12,14,46,47 , for which the phase of either pulse A (PP 2 ) or pulse B (PP 1 ) is irrelevant (see "Materials and methods" section for details). Therefore, they are directly linked to population dynamics emerging from the modulation of the QCL gain by saturation and pump-induced restoration of the population inversion. In contrast, the signatures at (ν L , ν L ), and (ν L , −2ν L ) are caused by four-wave mixing (4WM 1 and 4WM 2 ) processes originating from Fig. 2 Field-resolved 2D THz spectroscopy of a free-running QCL. a Experimental arrangement. The THz-waveforms E A and E B (blue) are focused onto the QCL facet and the transmitted waveform (E AB , orange) is detected electro-optically. b Electric field emitted by the nonlinear polarization, E NL (t, τ), of the unbiased QCL as a function of the delay times t and τ. The electric field is multiplied by a factor of 10 and color-coded as in c. c E NL (t, τ), as in b, but for biased QCL with a current of The spectrum is multiplied by a factor of 7 for ν τ < −5 THz and ν τ > 2.8 THz and color-coded as in e. e 2D-amplitude spectrumẼ NL ν t ; ν τ ð Þof the nonlinear response (I B = 840 mA) in c. For better visibility of the 6WM and 8WM signals, the corresponding part of the spectrum is multiplied by a factor of 7. f Corresponding simulated 2D-amplitude spectrumẼ NL ν t ; ν τ ð Þ purely coherent nonlinear polarization dynamics which depend on the phase of all participating photons (see "Materials and methods" section for details). The relative strength of the four-wave mixing signals amounts to as much as~1% of the amplitude of the transmitted THz pulses, E AB , indicating an unusually large third-order susceptibility of at least 10 10 pm 2 V −2 (see Supplementary  Information for details). Remarkably, the QCL's strongly nonlinear response also includes six-and even eightwave-mixing signals at (ν L , −3ν L ) and (ν L , −4ν L ), respectively, corresponding to the interaction of as many as eight photons from both THz pulses. To the best of our knowledge, these observations mark the highest-order nonlinearity reported for QCLs so far. This wealth of information enables us to extract valuable parameters of the microscopic electron dynamics in the free-running QCL. Since the nonlinearities studied here occur within a single pass through the QCL cavity, propagation effects do not play a dominant role and spatially averaged dynamics provide a meaningful description (Supplementary Information). As a first test, we quantitatively retrieve the gain recovery time, T gr , for a series of values of I B from below to above the laser threshold, which for THz QCLs has only been available for weak saturation to date 41 . For strong gain depletion, we proceed by selectively back-transforming the PP 1 signal into the time domain (Fig. 3a). The exponential decay of the corresponding nonlinear field component along the τ axis takes the form when evaluated at a fixed delay time, such as t = 1.0 ps (dotted black line in Fig. 3a). Here, A 1 is the signal amplitude and C is an offset accounting for thermal effects. Alternatively, T gr can be extracted from the measured time-domain signal E NL directly since E PP1 is the only contribution that does not oscillate in the τ direction (Fig. 3b, see Supplementary Information for details). The excellent agreement of both approaches allows us to determine T gr reliably as a function of the bias   (Fig. 3c), by rapidly measuring onedimensional sets of E NL for constant t. Starting at a value of T gr = 1.3 ps, below the laser threshold (I B = 680 mA), T gr grows with increasing I B until it reaches a maximum of (2.0 ± 0.2) ps at the laser threshold current of I th = 910 mA, followed by a drop upon further increase of I B .
The overall time constants are remarkably short given that the laser emission occurs far below the optical phonon energy and the lattice temperature is kept below 10 K. Microscopically, the characteristic current dependence of T gr can be understood within a rate equation model that goes beyond the picture of a pure quantum cascade amplifier by considering the electrically pumped laser transition in a two-level picture and the feedback by the laser cavity (see "Materials and methods" section). Below the threshold, increasing pump currents boost the population inversion w and hence the losses by amplified spontaneous emission, such that the recovery of w takes longer. When I B reaches the laser threshold, gain clamping limits a further increase of w and T gr is maximal. With yet higher pump currents, T gr finally drops because of faster repopulation of the upper laser level, as observed in Fig. 3c. Although longer relaxation times may still exist owing to the complex transport in a QCL band structure, the predominant gain relaxation time in our experiment is extremely fast (T gr ≤ 2 ps). Going beyond population dynamics, we can also analyze the nonlinear polarization of the free-running THz QCL for the first time. By back-transforming the four-wave mixing signal 4WM 1 (Fig. 3d) and fitting the experimental data at a fixed delay time (e.g., t = 1.0 ps) with the exponentially decaying oscillatory function we extract the polarization decay time, T 2 . Here, A 2 is the signal amplitude and ϕ is a constant phase. Given the excellent agreement of the above fit function with the experimental data (Fig. 3e), faithful values of the decay time of the nonlinear polarization can be extracted for pump currents below and above the threshold. Interestingly, T 2 monotonically increases from T 2 = 0.9 ps at I B = 780 mA by a factor of almost 1.6 to T 2 = 1.5 ps at I B = 1080 mA (Fig. 3f).

Discussion
Although analyzing selected nonlinear optical processes provides access to key information associated with the population and polarization dynamics, our access to a large number of different high-order nonlinearities (Fig. 2e) facilitates a much more comprehensive insight into the microscopic sub-cycle dynamics. As a first demonstration, we utilize the optical Bloch equations to track the full THz-induced coherent trajectory of the Bloch vector of the resonant laser transition, and explain the emergence of the various nonlinear optical processes all at once. The gain medium is described by a density matrix ρ coupled to the self-consistently calculated THz cavity electric field, E(t), leading to a Hamiltonian of the form: Here, μ 12 and ν L represent the dipole moment and the frequency of the laser transition, respectively. The dynamics are calculated beyond the rotating-wave approximation by solving the von Neumann equation. The equations of motion for the density matrix elements are given as: In the stationary state without THz excitation, we assume the QCL to exhibit a slight population inversion w 0 ¼ ρ 0 22 À ρ 0 11 . The electric field E(t) is calculated by considering the incident THz field, E THz , as well as the electric field reradiated by the electronic polarization, ρ 12 : where μ 0 denotes the vacuum permeability, c is the speed of light, and n is the refractive index of the material. Γ is a constant accounting for the coupling of the electron ensemble to the cavity field (Supplementary Information). All parameters are taken from the literature or auxiliary measurements as far as available. Only T 1 and T 2 are used as fitting parameters. Analogously to the experiment, we derive the nonlinear field E NL from calculations of the individual configurations including both or only one of the THz pulses. The resulting 2D-amplitude spectrum (Fig. 2f) reproduces the emergence of all nonlinear signals including pump-probe, four-, six-, and eight-wave-mixing signals. Furthermore, the simulation reveals the underlying temporal evolution of the microscopic QCL dynamics including the population inversion w = ρ 22 -ρ 11 and the polarization ρ 12 = u + iv as a function of t and τ. Figure 4a shows the evolution of the QCL's state in a reduced Bloch sphere representation for τ = 0 ps. Unlike in an equilibrium two-level system, the initial state of the Bloch vector is not given by the south pole of the Bloch sphere, but an incoherent state (u = v = 0) located slightly above w = 0 corresponding to a moderate population inversion of the gain medium, prepared by the pump current. The intense THz field drives the Bloch vector onto a coherent trajectory, which, unlike in a conventional equilibrium two-level system, points downwards as the pump pulse coherently depletes the population inversion by stimulated emission within a single cycle of the THz carrier wave. Subsequently, electrical pumping restores the population inversion while the free induction decay of the induced polarization causes a spiraling trajectory back to the stationary incoherent starting point. Because the THz field is resonant to the laser transition, these dynamics and all multi-wave mixing nonlinearities associated with them are resonantly enhanced. Yet, resonant enhancement in a running laser is associated with gain rather than absorption, making QCLs a highly attractive nonlinear optical medium. Importantly, the coherent Bloch trajectory sensitively depends on the bias current. The close-up in Fig. 4b highlights the dynamics for a low current of I B = 780 mA, where the coherent polarization given by the excursion in the u and v directions of the Bloch trajectory rapidly collapses. In contrast, for a bias current of I B = 1080 mA, which is above threshold (Fig. 4c), the polarization maintains a much larger amplitude until the population inversion is fully restored. To establish a solid microscopic understanding of these bias-dependent coherent dynamics, we simulated the underlying electronic transport with a density matrix ensemble Monte Carlo approach coupled to a Schrödinger-Poisson solver. This allows us to extract the relevant scattering, tunneling, and dephasing rates self-consistently. We accounted for eight states per QCL period. Figure 4d, e displays the resulting band diagram of the QCL under two typical bias fields of 1.1 kV cm −1 and 1.7 kV cm −1 together with the computed dephasing time, T 2 , of the lasing transition (Fig. 4e, inset). The theory clearly confirms a monotonic rise of T 2 with I B as observed experimentally (Fig. 3f). In our simulations, this behavior The blue and orange lines represent the squared modulus of the wave function of the upper and lower laser level (|Ψ 2 | 2 and |Ψ 1 | 2 , respectively). Additionally, the highest miniband state |Ψ mb1 | 2 as well as the other miniband states |Ψ mb | 2 are shown in black and gray. e As in d but for a bias voltage of F = 1.7 kV cm −1 . Inset: Calculated dependence of T 2 on the bias current I B is mainly caused by the decreasing overlap of the lower laser level with other levels with increasing bias. As seen in Fig. 4d, e, a stronger bias induces a more pronounced energy spread of the minibands, which reduces the rate of electron-electron scattering between them and limits scattering out of the lower laser level. This behavior is supported by the bias-dependent spatial overlap of the miniband states. For a low bias field of 1.1 kV cm −1 (Fig. 4d), the squared wave functions of the miniband states (e.g., |Ψ mb1 | 2 , black curve) are delocalized over almost the full injector period whereas they spatially separate with increasing bias fields (1.7 kV cm −1 , Fig. 4e), suppressing phase destroying scattering events in the miniband states. Since our model does not explicitly consider back action by the laser cavity, we restricted our simulations to currents below the threshold. The new 2D spectroscopy data hence provide an important benchmark for future quantum theories treating both electronic transport processes and intracavity light-matter interaction on equal sub-cycle footing.
In conclusion, the first 2D spectroscopy of a freerunning THz QCL revealed strong nonlinearities including up to eight-wave mixing, on a sub-cycle time scale. Besides extremely fast gain recovery times down to 1.2 ps, our approach reveals the coherent trajectory of the Bloch vector itself, unlocking a previously inaccessible dimension of electron dynamics in QCLs. While our experiments were carried out with a prototypical design of a THz QCL, we expect this new comprehensive access to coherent and incoherent response functions to be broadly applicable to all kinds of innovative device architectures. In particular, mode-locked THz QCLs will depend on an in-depth understanding of the carrier-wave dynamics of the gain, the coherent polarization, and electronic transport for dispersion management and single-cycle pulse formation. Moreover, highly efficient and tunable multiwave mixing opens up exciting perspectives for intracavity frequency conversion, multiplexing, potentially even quantum squeezing, and strong-field light-matter coupling in a single compact electrically pumped device.

Experimental setup
Phase-locked THz pulses are generated by optical rectification of near-infrared laser pulses of a duration of 220 fs from an Yb:KGW amplifier (repetition rate, 1 MHz) in a cryogenically cooled (T = 77 K) lithium niobate crystal in a tilted pulse-front configuration. The THz waveforms ( Supplementary Fig. S1) feature peak fields of up to 3 kV cm −1 and a duration of 0.84 ps of the intensity envelope (FWHM), corresponding to 1.3 optical cycles of the carrier frequency of 1.6 THz ( Supplementary Fig. S1, inset). A Michelson interferometer with a silicon beamsplitter creates two identical copies of the THz pulses which are subsequently delayed with respect to each other by a variable delay time τ, before being recombined collinearly. The pulse pair is coupled into the facet of the QCL within the cryostat and the transmitted THz radiation is focused onto a GaP crystal of a thickness of 1 mm for electro-optic detection, using a small portion of the laser pulse as a gate. In each of the two interferometer arms, a mechanical chopper is placed in an intermediate focus, enabling independent modulation of both THz pulse trains.

2D THz spectroscopy
We employ 2D THz spectroscopy to disentangle nonlinear processes inside the QCL by a Liouville path analysis which relates each nonlinear signature in the twodimensional spectra to a specific wave vector combination of the incident fields 47 . For example, χ (3) nonlinearities, such as pump-probe and four-wave mixing processes, combine three photons of the two incident fields. In case of the pump-probe signal PP 1 at frequencies of (v t , v τ ) = (v L , 0), a third-order polarization is created by one photon from the incident pulse A and two photons from pulse B, represented by the complexvalued spectral amplitudesẼ A andẼ B , respectively. Since the field componentsẼ B ðν L Þ andẼ Ã B ðÀν L Þ are each other's complex conjugate with mutually inverted phases (see inverted sign of their frequencies), the phase of field B cancels. Only the phase of the probe photon from pulse A remains imprinted in the nonlinear polarization P 3 ð Þ PP 1 . Therefore, PP 1 is incoherent with respect to the phase of pulse B, whereas PP 2 , where A and B switch roles, is incoherent with respect to the phase of pulse A. In contrast, four-wave mixing nonlinearities preserve the phases of both pulses. Considering the 4WM 1 signal where (v t , v τ ) = (v L , v L ), two photons from pulse A and one photon from pulse B of inverted phase contribute to the resulting third-order polarization: In conclusion, nonlinear pump-probe signals relate to the incoherent population dynamics of the system, while four-wave mixing or higher-order multi-wave mixing signals represent the coherent nonlinear polarization.
A conventional graphical representation of this Liouville path decomposition plots the contributing wave vectors in 2D frequency space, as shown in Fig. 2e. Here, each signature can be reached by a path consisting of an integer linear combination of the wave vectors k A = (v L , 0) and k B = (v L , −v L ), each representing one photon from pulse A or B, respectively. Correspondingly, the signatures at k PP1 = (v L , 0) and k PP2 = (v L , −v L ) belong to the Liouville paths k PP1 = k A + k B − k B and k PP2 = k A − k A + k B , respectively, while four-wave mixing signals occurring at

Rate equation model for the QCL gain dynamics
To model the gain recovery time of the QCL we assume a two-level system with an electron reservoir and a drain for the upper and lower level, respectively, located inside the QCL cavity. The electron number N 1 and N 2 in the lower and upper laser level as well as the photon number N ph in the cavity are described by the following rate equations: Here, τ 1 and τ 2 describe the lifetimes of the lower and upper laser level. The losses of the cavity are considered by a corresponding cavity lifetime, τ cav . To fit the experimental data and in general agreement with literature values 36,48,49 , we choose τ 1 = 4.3 ps, τ 2 = 5.5 ps and τ cav = 1.2 ps. Stimulated emission is taken into account by the Einstein coefficient in the dipole approximation, B 21 ¼ 4π 2 3 h 2 1 4πε 0 hΨ 2 er j jΨ 1 i j j 2 , and the spectral energy density, ρ 0 Á N ph ¼ hv L V Δv L Á N ph . Here, Δv L ¼ 1 2π 1 τ 1 þ 1 τ 2 is the bandwidth of the laser transition at v L = 2.2 THz and V = A Laser •d describes the mode volume. The transition dipole moment μ 21 ¼ hΨ 2 er j jΨ 1 i ¼ 8:5 nm Á e with the elementary charge, e, is extracted from a band structure simulation of the QCL. The injection of electrons into level 2 is determined by the pump rate G and limited by the pre-factor (N e − N 2 ) implementing Pauli blocking by including a maximal electron number in the upper level of N e = 0.9•10 9 .
At delay time t = t perturb = 0 ps, we increase the photon population using a short, delta-shaped term β⋅δ(t − t perturb ). Here, β ¼ 5:7 Á 10 7 ps À1 describes the strength of the perturbation and is chosen to fit the experimental data and to increase N ph by approximately an order of magnitude. The subsequent dynamics of the cavity photon number back to quasi-equilibrium are fitted with a decaying exponential function in order to extract the gain recovery time.

Microscopic carrier transport model of the QCL
The dephasing time T 2 of the lasing transition in the inset of Fig. 4e has been computed based on self-consistent stationary carrier transport simulations, using a density matrix ensemble Monte Carlo approach 50 . Here, the energy quantization due to electron confinement in growth direction as well as the in-plane electron wave vector are considered, allowing us to include both interand intrasubband processes. Space charge effects are considered by coupling the carrier transport simulations to a Schrödinger-Poisson solver, providing the wave functions and eigenenergies of the quantized electron states (Fig. 4d, e). Our carrier transport model accounts for tunneling through the thick injection barrier as well as all relevant scattering mechanisms, such as electron-electron interactions beyond mean-field theory, scattering with acoustic and longitudinal-optical phonons (also accounting for nonequilibrium phonons), and elastic scattering due to impurities and interface roughness. The dephasing time is modeled as 1/T 2 = (1/T u + 1/T l )/2 + γ p . Here, T u and T l denote the upper and lower laser level lifetimes obtained from the corresponding intersubband scattering rates, and γ p is the pure dephasing rate due to intrasubband processes, which is also self-consistently calculated in our simulation approach 50 .