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

J. Riepl, J. Raab, P. Abajyan, H. Nong, J. R. Freeman, L. H. Li, E. H. Linfield, A. G. Davies, A. Wacker, T. Albes, et al.

► To cite this version:

J. Riepl, J. Raab, P. Abajyan, H. Nong, J. R. Freeman, et al. Field-resolved high-order sub-cycle nonlinearities in a terahertz semiconductor laser. Light: Science and Applications, 2021, 10 (1), pp.246. 10.1038/s41377-021-00685-5. hal-03510075

HAL Id: hal-03510075

https://hal.sorbonne-universite.fr/hal-03510075

Submitted on 4 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Field-resolved high-order sub-cycle nonlinearities in a terahertz semiconductor laser

J. Riepl¹, J. Raab¹, P. Abajyan², H. Nong², J. R. Freeman¹,³, L. H. Li¹,³, E. H. Linfield¹,³, A. G. Davies¹,³, A. Wacker⁴, T. Albes⁵, C. Jirauschek¹,⁵, C. Lange¹,⁶, S. S. Dhillon², and R. Huber¹

## Abstract

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.

## Introduction

The terahertz (THz) window of the electromagnetic spectrum is at the heart of an ongoing revolution in optical sciences and technology. Ultrashort THz pulses have enabled resonant access to a broad variety of low-energy elementary excitations in condensed matter, such as plasmons, phonons, magnons, and Higgs modes in superconductors¹⁻⁷. Strong THz fields have paved the way to extreme nonlinearities²,⁸⁻¹¹ such as multi-wave mixing¹²⁻¹⁵, high-harmonic generation⁹,¹⁶, lightwave electronics⁸, and ultrafast nanoscopy¹⁰,¹¹,¹⁷,¹⁸. Moreover, a wide range of novel photonics applications has employed THz radiation including high-resolution spectroscopy¹,¹⁹, near-field microscopy²⁰, and telecommunication concepts²¹. Thus, powerful, efficient, and compact sources of THz pulses and

broadband frequency combs are in high demand, for which THz quantum cascade lasers (QCLs)²² have presented themselves as a particularly promising solution. These electrically pumped semiconductor lasers based on inter-subband transitions and resonant tunneling in quantum wells have undergone considerable development over the last two decades. Important milestones include high-temperature²³,²⁴, high-power²⁵, tunable²⁶, and frequency comb operation²⁷,²⁸, as well as short pulse generation²⁹⁻³¹.

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 pump-probe experiments³²⁻³⁶, 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, four-wave mixing has been observed in time-integrated experiments³⁷. In THz pump-probe measurements using weak pulses, incoherent gain dynamics of THz QCLs³⁸⁻⁴⁰

<figure><img src="image_1.png" /></figure>

<footnote>Correspondence: S. S. Dhillon (Sukhdeep.Dhillon@phys.ens.fr)
¹Department of Physics, University of Regensburg, Regensburg, Germany
²Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, Paris, France
Full list of author information is available at the end of the article</footnote>
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 free-running THz QCL. Phase-stable, single-cycle THz transients with $kV \text{ 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$ hetero-structure (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 wet-chemical 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

<figure><img src="image_2.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_5.png" /></figure>

**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 electro-optically. 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 $\tilde{E}(\nu)$ of the waveforms in **b** and lasing spectrum $\tilde{A}(\nu)$ 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 $\text{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⁻¹ 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 hetero-structure 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, $\tilde{E}_b(v)$, features an increase over the unbiased situation, $\tilde{E}_u(v)$. 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³⁴.

For a more quantitative analysis of the gain, Fig. 1d shows the ratio of transmission, $\tilde{E}_b(v)/\tilde{E}_u(v)$ (red curve), measured with an incident THz peak field of 1.4 kV cm⁻¹ (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 $\alpha_{E_0} \approx 2.2 \text{ 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 ($\alpha_{E_0/2} \approx 3.4 \text{ 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³⁹. 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%⁴⁵ 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¹²,¹⁴,⁴⁶,⁴⁷ (see “Materials and methods” section). To this end, two identical THz waveforms $E_A(t)$ and $E_B(t, \tau)$, which are mutually delayed by a variable time $\tau$, 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, \tau)$ 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, \tau)$. Synchronous mechanical chopping of the two incident THz beams $E_A(t)$ and $E_B(t, \tau)$ isolates their individual contributions, and by subtracting them from the transmitted signal $E_{AB}(t, \tau)$, we are left with the nonlinear signal $E_{NL}(t, \tau) = E_{AB}(t, \tau) - E_A(t) - E_B(t, \tau)$.

Figure 2b shows the resulting 2D field map $E_{NL}(t, \tau)$ 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, $\tau$. The strongest nonlinearities occur near $\tau = 0$ ps, where both pulses overlap temporarily. By varying $\tau$, the total incident field is strongly modulated leading to a pronounced variation of $E_{NL}$. Even for delay times $|\tau| > 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 $\tau$. The nonlinear response manifests as oscillations at the transition frequency of the unbiased QCL at $\nu_{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,
<figure><img src="image_6.png" /><figcaption>**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, \tau)$, of the unbiased QCL as a function of the delay times $t$ and $\tau$. The electric field is multiplied by a factor of 10 and color-coded as in **c**. **c** $E_{NL}(t, \tau)$, as in **b**, but for biased QCL with a current of $I_B = 840$ mA. **d** 2D-amplitude spectrum $\tilde{E}_{NL}(\nu_t, \nu_\tau)$ of the nonlinear response ($I_B = 0$) in **b**. The spectrum is multiplied by a factor of 7 for $\nu_t < -5$ THz and $\nu_\tau > 2.8$ THz and color-coded as in **e**. **e** 2D-amplitude spectrum $\tilde{E}_{NL}(\nu_t, \nu_\tau)$ 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 $\tilde{E}_{NL}(\nu_t, \nu_\tau)$</figcaption></figure>

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 $\nu_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 $|\tau|$, 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 $|\tau|$ 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 $\Delta\varphi = \tau\nu_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 non-linear processes, we perform a 2D Fourier transform of the time-domain data of Fig. 2c. By comparing the

resulting amplitude spectrum $\tilde{E}_{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 $\tilde{E}_{NL}$ are located at $(\nu_D, \nu_\tau) = (\nu_L, n \times \nu_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 $(\nu_L, 0)$, and $(\nu_L, -\nu_L)$, respectively. Starting from the origin (0,0), every maximum of $\tilde{E}_{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 $(\nu_L, 0)$ and $(\nu_L, -\nu_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 $(\nu_L, \nu_L)$, and $(\nu_L, -2\nu_L)$ are caused by four-wave mixing (4WM$_1$ and 4WM$_2$) processes originating from
<figure><img src="image_7.png" /><figcaption>Fig. 3 Quantitative analysis of the coherent and incoherent dynamics in the QCL. a Back-transformed electric field corresponding to the pump-probe signal $E_{PP1}$, b Slice through $E_{PP1}(t, \tau)$ at a delay time $t = 1.0$ ps (see black dashed line in a). Blue dashed line: numerical fit of function $A_1 \exp(-\tau/T_{gr}) + C$. c Gain recovery time $T_{gr}$ of the QCL as a function of the bias current, extracted from the PP1 signal (red dots) and the full nonlinear time domain signal (gray). The error bars display the uncertainty of the numerical fit. Light blue dashed line: $T_{gr}$ as a function of the pump rate $G$, extracted from the rate equation model. d Back-transformed electric field corresponding to the four-wave mixing signal $E_{4WM1}$. e Slice through $E_{4WM1}(t, \tau)$ at an example delay time $t = 1.0$ ps (see black dashed line in d). Blue dashed line: numerical fit of function $A_2 \sin(2\pi\nu_L\tau + \phi) \exp(-\tau/T_{gr})$. f Decay time of the coherent population $T_2$ as a function of the bias current. The error displays the uncertainty of $T_2$ for different delay times between 0.45 and 1.75 ps</figcaption></figure>

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} \text{ pm}^2 \text{ V}^{-2}$ (see Supplementary Information for details). Remarkably, the QCL’s strongly nonlinear response also includes six- and even eight-wave-mixing signals at ($v_L, -3v_L$) and ($v_L, -4v_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 PP1 signal into the time domain (Fig. 3a). The exponential decay of the corresponding nonlinear field component along the $\tau$ axis takes the form

$$E_{PP1}(t = 1.0 \text{ ps}, \tau) = A_1 \exp\left(-\frac{\tau}{T_{gr}}\right) + C \quad (1)$$

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 $\tau$ 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
current, $I_B$ (Fig. 3c), by rapidly measuring one-dimensional 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} \le 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

$$E_{4WM1}(t = 1.0 \text{ ps}, \tau) = A_2 \sin(2\pi\nu_L\tau + \phi) \exp\left(-\frac{\tau}{T_2}\right) \quad (2)$$

we extract the polarization decay time, $T_2$. Here, $A_2$ is the signal amplitude and $\phi$ 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 $\rho$ coupled to the self-consistently calculated THz cavity electric field, $E(t)$, leading to a Hamiltonian of the form:

$$H = \begin{pmatrix} 0 & \mu_{12}E(t) \\ \mu_{12}E(t) & h\nu_L \end{pmatrix} \quad (3)$$

Here, $\mu_{12}$ and $\nu_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:

$$\dot{\rho}_{11} = \frac{i}{\hbar}\mu_{12}E(t)(\rho_{12} - \rho_{12}^*) - \frac{1}{T_{gr}}(\rho_{11} - \rho_{11}^0) \quad (4a)$$

$$\dot{\rho}_{22} = \frac{i}{\hbar}\mu_{12}E(t)(\rho_{12}^* - \rho_{12}) - \frac{1}{T_{gr}}(\rho_{22} - \rho_{22}^0) \quad (4b)$$

$$\dot{\rho}_{12} = \frac{i}{\hbar}\mu_{12}E(t)(\rho_{11} - \rho_{22}) - i2\pi\nu_L\rho_{12} - \frac{\rho_{12}}{T_2} = \dot{\rho}_{21}^* \quad (4c)$$

In the stationary state without THz excitation, we assume the QCL to exhibit a slight population inversion $w_0 = \rho_{22}^0 - \rho_{11}^0$. 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, $\rho_{12}$:

$$E(t) = E_{THz} - \Gamma\mu_0 \frac{c}{2n}\dot{\rho}_{12} \quad (5)$$

where $\mu_0$ denotes the vacuum permeability, $c$ is the speed of light, and $n$ is the refractive index of the material. $\Gamma$ 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 = \rho_{22} - \rho_{11}$ and the polarization $\rho_{12} = u + iv$ as a function of $t$ and $\tau$.
<figure><img src="image_8.png" /><figcaption>Fig. 4 Microscopic picture of the QCL dynamics. a Temporal evolution of the QCL's density matrix, shown in a reduced Bloch sphere with a radius of 0.2 for a delay time of $\tau = 0$ ps and a bias current of $I_B = 780$ mA. b Close-up view of the trajectory of the Bloch vector shown in a. c As b, but for $I_B = 1080$ mA. d Simulated band diagram of the QCL under an electric field of $F = 1.1$ kV cm$^{-1}$. The blue and orange lines represent the squared modulus of the wave function of the upper and lower laser level ($|\Psi_2|^2$ and $|\Psi_1|^2$, respectively). Additionally, the highest miniband state $|\Psi_{mb1}|^2$ as well as the other miniband states $|\Psi_{mb2}|^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$.</figcaption></figure>

Figure 4a shows the evolution of the QCL's state in a reduced Bloch sphere representation for $\tau = 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
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⁻¹ (Fig. 4d), the squared wave functions of the miniband states (e.g., |Ψ<sub>mb1</sub>|², black curve) are delocalized over almost the full injector period whereas they spatially separate with increasing bias fields (1.7 kV cm⁻¹, 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 free-
running THz QCL revealed strong nonlinearities includ-
ing 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 dimen-
sion of electron dynamics in QCLs. While our experi-
ments 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 trans-
port for dispersion management and single-cycle pulse
formation. Moreover, highly efficient and tunable multi-
wave mixing opens up exciting perspectives for intracavity
frequency conversion, multiplexing, potentially even
quantum squeezing, and strong-field light-matter cou-
pling in a single compact electrically pumped device.

**Materials and methods**

Experimental setup

Phase-locked THz pulses are generated by optical rec-
tification 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⁻¹ 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 beams-
plitter 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 col-
linearly. The pulse pair is coupled into the facet of the
QCL within the cryostat and the transmitted THz radia-
tion 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 non-
linear processes inside the QCL by a Liouville path ana-
lysis which relates each nonlinear signature in the two-
dimensional spectra to a specific wave vector combination
of the incident fields⁴⁷. For example, χ⁽³⁾ 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₁ at frequencies of (νₜ, ντ) =
(νₗ, 0), a third-order polarization

$$
P_{PP_1}^{(3)}(\nu_L) = \epsilon_0 \chi^{(3)}(\nu_L, \nu_L, -\nu_L) \tilde{\mathcal{E}}_A(\nu_L) \tilde{\mathcal{E}}_B(\nu_L) \tilde{\mathcal{E}}_B^*(-\nu_L) \quad (6)
$$

is created by one photon from the incident pulse A and
two photons from pulse B, represented by the complex-
valued 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<sup>(3)</sup><sub>PP₁</sub>.
Therefore, PP₁ is incoherent with respect to the phase of
pulse B, whereas PP₂, 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₁ signal
where (ν<sub>t</sub>, ν<sub>τ</sub>) = (ν<sub>L</sub>, ν<sub>L</sub>), two photons from pulse A and
one photon from pulse B of inverted phase contribute to
the resulting third-order polarization:

$$
P_{4\text{WM}_1}^{(3)}(\nu_L) = \epsilon_0 \chi^{(3)}(\nu_L, \nu_L, -\nu_L) \tilde{\mathcal{E}}_A(\nu_L) \tilde{\mathcal{E}}_A(\nu_L) \tilde{\mathcal{E}}_B^*(-\nu_L) \quad (7)
$$

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 sig-
nature can be reached by a path consisting of an integer
linear combination of the wave vectors **k**<sub>A</sub> = (ν<sub>L</sub>, 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 $k_{4WM1} = (v_L, v_L)$ and $k_{4WM2} = (v_L, -2v_L)$ are decomposed to $k_{4WM1} = k_A + k_A - k_B$ and $k_{4WM2} = k_B + k_B - k_A$.

### 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:

$$ \frac{dN_1}{dt} = -\frac{N_1}{\tau_1} + \frac{N_2}{\tau_2} + B_{21}(N_2 - N_1) \cdot \rho_0 N_{ph} \quad (8) $$

$$ \frac{dN_2}{dt} = (Ne - N_2) \cdot G - \frac{N_2}{\tau_2} - B_{21}(N_2 - N_1) \cdot \rho_0 N_{ph} \quad (9) $$

$$ \frac{dN_{ph}}{dt} = -\frac{N_{ph}}{\tau_{cav}} + B_{21}(N_2 - N_1) \cdot \rho_0 N_{ph} + \beta \cdot \delta(t - t_{pert}) \quad (10) $$

Here, $\tau_1$ and $\tau_2$ describe the lifetimes of the lower and upper laser level. The losses of the cavity are considered by a corresponding cavity lifetime, $\tau_{cav}$. To fit the experimental data and in general agreement with literature values$^{36,48,49}$, we choose $\tau_1 = 4.3$ ps, $\tau_2 = 5.5$ ps and $\tau_{cav} = 1.2$ ps. Stimulated emission is taken into account by the Einstein coefficient in the dipole approximation, $B_{21} = \frac{4\pi^2}{3\hbar^2} \left(\frac{1}{4\pi\epsilon_0}\right) |\langle\Psi_2|e\vec{r}|\Psi_1\rangle|^2$, and the spectral energy density, $\rho_0 \cdot N_{ph} = \frac{h\nu_L}{V\Delta\nu_L} \cdot N_{ph}$. Here, $\Delta\nu_L = \frac{1}{2\pi}(\frac{1}{\tau_1} + \frac{1}{\tau_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 $\mu_{21} = \langle\Psi_2|e\vec{r}|\Psi_1\rangle = 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 \cdot 10^9$.

At delay time $t = t_{perturb} = 0$ ps, we increase the photon population using a short, delta-shaped term $\beta \cdot \delta(t - t_{perturb})$. Here, $\beta = 5.7 \cdot 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 inter- and 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_i)/2 + \gamma_p$. Here, $T_u$ and $T_i$ denote the upper and lower laser level lifetimes obtained from the corresponding intersubband scattering rates, and $\gamma_p$ is the pure dephasing rate due to intrasubband processes, which is also self-consistently calculated in our simulation approach$^{50}$.

### Acknowledgements

We acknowledge funding from the European Union under the Horizon 2020 research and innovation programs FET-Open grant ULTRAQCL 665158, the French National Research Agency (ANR-18-CE24-0013-02 - "TERASEL"), the Swedish Research Council (project 2017-04287) and the EPSRC (UK) program grant "HyperTerahertz" (EP/P021859/1).

### Author details

¹Department of Physics, University of Regensburg, Regensburg, Germany.

²Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, Paris, France. ³School of Electronic and Electrical Engineering, University of Leeds, Woodhouse Lane, Leeds, UK. ⁴Mathematical Physics and NanoLund, Lund University, Lund, Sweden. ⁵Department of Electrical and Computer Engineering, Technical University of Munich, Munich, Germany. ⁶Department of Physics, TU Dortmund University, Dortmund, Germany

### Author contributions

C.L., S.S.D., and R.H. conceived the study. J.R.F., L.H.L., E.H.L., and A.G.D. designed, epitaxially grew, and fabricated the quantum cascade laser. H.N. and S.S.D. characterized the QCL. J.Ra., J.Ri., C.L., and R.H. carried out the 2D spectroscopy. J.Ri., J.Ra., C.L., S.S.D., and R.H. analyzed the experimental results. A.W. developed the model for the gain dynamics. J.Ri., C.L., and R.H. implemented the rate equation model. P.A. simulated the band structure. J.Ri., J.Ra., and C.L. set up and carried out the density matrix formalism. T.A. and C.J. simulated the band structure and the microscopic carrier transport in the QCL. J.Ri., S.S.D., and R.H. wrote the manuscript with contributions from all authors.

### Funding

Open Access funding enabled and organized by Projekt DEAL.
**Competing interests**

The authors declare no competing interests.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41377-021-00685-5.

Received: 22 June 2021 Revised: 18 November 2021 Accepted: 23 November 2021

Published online: 20 December 2021

**References**

1. Dhillon, S. S. et al. The 2017 terahertz science and technology roadmap. *J. Phys. D Appl. Phys.* **50**, 043001 (2017).

2. Kampfrath, T., Tanaka, K. & Nelson, K. A. Resonant and nonresonant control over matter and light by intense terahertz transients. *Nat. Photonics* **7**, 680–690 (2013).

3. Makihara, T. et al. Ultrastrong magnon-magnon coupling dominated by antiresonant interactions. *Nat. Commun.* **12**, 3115 (2021).

4. Vaswani, C. et al. Light quantum control of persisting Higgs modes in iron-based superconductors. *Nat. Commun.* **12**, 258 (2021).

5. Tonouchi, M. Cutting-edge terahertz technology. *Nat. Photonics* **1**, 97–105 (2007).

6. Matsunaga, R. et al. Light-induced collective pseudospin precession resonating with Higgs mode in a superconductor. *Science* **345**, 1145–1149 (2014).

7. Jepsen, P. U., Cooke, D. G. & Koch, M. Terahertz spectroscopy and imaging—modern techniques and applications. *Laser Photonics Rev.* **5**, 124–166 (2011).

8. Reimann, J. et al. Subcycle observation of lightwave-driven Dirac currents in a topological surface band. *Nature* **562**, 396–400 (2018).

9. Schmid, C. P. et al. Tunable non-integer high-harmonic generation in a topological insulator. *Nature* **593**, 385–390 (2021).

10. Peller, D. et al. Sub-cycle atomic-scale forces coherently control a single-molecule switch. *Nature* **585**, 58–62 (2020).

11. Plankl, M. et al. Subcycle contact-free nanoscopy of ultrafast interlayer transport in atomically thin heterostructures. *Nat. Photonics* **10**, 594–600 (2021).

12. Kuehn, W. et al. Two-dimensional terahertz correlation spectra of electronic excitations in semiconductor quantum wells. *J. Phys. Chem. B* **115**, 5448–5455 (2011).

13. Zaks, B., Liu, R. B. & Sherwin, M. S. Experimental observation of electron-hole recollisions. *Nature* **483**, 580–583 (2012).

14. Maag, T. et al. Coherent cyclotron motion beyond Kohn's theorem. *Nat. Phys.* **12**, 119–123 (2016).

15. Lu, J. et al. Coherent two-dimensional terahertz magnetic resonance spectroscopy of collective spin waves. *Phys. Rev. Lett.* **118**, 207204 (2017).

16. Hafez, H. A. et al. Extremely efficient terahertz high-harmonic generation in graphene by hot Dirac fermions. *Nature* **561**, 507–511 (2018).

17. Yoshioka, K. et al. Real-space coherent manipulation of electrons in a single tunnel junction by single-cycle terahertz electric fields. *Nat. Photonics* **10**, 762–765 (2016).

18. Klarskov, P. et al. Nanoscale laser terahertz emission microscopy. *ACS Photonics* **4**, 2676–2680 (2017).

19. Hübers, H. W., Richter, H. & Wienold, M. High-resolution terahertz spectroscopy with quantum-cascade lasers. *J. Appl. Phys.* **125**, 151401 (2019).

20. Pogna, E. A. A. et al. Unveiling the detection dynamics of semiconductor nanowire photodetectors by terahertz near-field nanoscopy. *Light Sci. Appl.* **9**, 189 (2020).

21. Freeman, W. Decoherence in high barrier quantum cascade structures. *J. Appl. Phys.* **122**, 045701 (2017).

22. Köhler, R. et al. Terahertz semiconductor-heterostructure laser. *Nature* **417**, 156–159 (2002).

23. Kainz, M. A. et al. Thermoelectric-cooled terahertz quantum cascade lasers. *Opt. Express* **27**, 20688–20693 (2019).

24. Khalatpour, A. et al. High-power portable terahertz laser systems. *Nat. Photonics* **15**, 16–20 (2021).

25. Li, L. H. et al. Terahertz quantum cascade lasers with >1 W output powers. *Electron. Lett.* **50**, 309–311 (2014).

26. Curwen, C. A., Reno, J. L. & Williams, B. S. Broadband continuous single-mode tuning of a short-cavity quantum-cascade VECSEL. *Nat. Photonics* **13**, 855–859 (2019).

27. Burghoff, D. et al. Terahertz laser frequency combs. *Nat. Photonics* **8**, 462–467 (2014).

28. Rösch, M. et al. Octave-spanning semiconductor laser. *Nat. Photonics* **9**, 42–47 (2015).

29. Barbieri, S. et al. Coherent sampling of active mode-locked terahertz quantum cascade lasers and frequency synthesis. *Nat. Photonics* **5**, 306–313 (2011).

30. Wang, F. H. et al. Short terahertz pulse generation from a dispersion compensated modelocked semiconductor laser. *Laser Photonics Rev.* **11**, 1700013 (2017).

31. Hillbrand, J. et al. In-phase and anti-phase synchronization in a laser frequency comb. *Phys. Rev. Lett.* **124**, 023901 (2020).

32. Eickemeyer, F. et al. Ultrafast coherent electron transport in semiconductor quantum cascade structures. *Phys. Rev. Lett.* **89**, 047402 (2002).

33. Kröll, J. et al. Phase-resolved measurements of stimulated emission in a laser. *Nature* **449**, 698–701 (2007).

34. Choi, H. et al. Femtosecond dynamics of resonant tunneling and superlattice relaxation in quantum cascade lasers. *Appl. Phys. Lett.* **92**, 122114 (2008).

35. Choi, H. et al. Gain recovery dynamics and photon-driven transport in quantum cascade lasers. *Phys. Rev. Lett.* **100**, 167401 (2008).

36. Kuehn, W. et al. Ultrafast phase-resolved pump-probe measurements on a quantum cascade laser. *Appl. Phys. Lett.* **93**, 151106 (2008).

37. Friedli, P. et al. Four-wave mixing in a quantum cascade laser amplifier. *Appl. Phys. Lett.* **102**, 222104 (2013).

38. Green, R. P. et al. Gain recovery dynamics of a terahertz quantum cascade laser. *Phys. Rev. B* **80**, 075303 (2009).

39. Bacon, D. R. et al. Gain recovery time in a terahertz quantum cascade laser. *Appl. Phys. Lett.* **108**, 81104 (2016).

40. Deml, C. G. et al. Gain dynamics in a heterogeneous terahertz quantum cascade laser. *Appl. Phys. Lett.* **113**, 181102 (2018).

41. Markmann, S. et al. Two-dimensional coherent spectroscopy of a THz quantum cascade laser: observation of multiple harmonics. *Opt. Express* **25**, 21753–21761 (2017).

42. Markmann, S. et al. Two-dimensional spectroscopy on a THz quantum cascade structure. *Nanophotonics* **10**, 171–180 (2020).

43. Luo, C. W. et al. Phase-resolved nonlinear response of a two-dimensional electron gas under femtosecond intersubband excitation. *Phys. Rev. Lett.* **92**, 047402 (2004).

44. Worrall, C. et al. Continuous wave operation of a superlattic quantum cascade laser emitting at 2 THz. *Opt. Express* **14**, 171–181 (2006).

45. Li, L. H. et al. The MBE growth and optimization of high performance terahertz frequency quantum cascade lasers. *Opt. Express* **23**, 2720–2729 (2015).

46. Raab, J. et al. Ultrafast terahertz saturable absorbers using tailored intersubband polaritons. *Nat. Commun.* **11**, 4290 (2020).

47. Mornhinweg, J. et al. Tailored subcycle nonlinearities of ultrastrong light-matter coupling. *Phys. Rev. Lett.* **126**, 177404 (2021).

48. Tzenov, P. et al. Passive and hybrid mode locking in multi-section terahertz quantum cascade lasers. *N. J. Phys.* **20**, 053055 (2018).

49. Wang, F. H. et al. Ultrafast response of harmonic modelocked THz lasers. *Light Sci. Appl.* **9**, 51 (2020).

50. Jirauschek, C. Density matrix Monte Carlo modeling of quantum cascade lasers. *J. Appl. Phys.* **122**, 133105 (2017).