Frequency comb generation by Bloch gain induced giant Kerr nonlinearity

Optical nonlinearities are known to provide a coherent coupling between the amplitude and phase of the light, which can result in the formation of periodic waveforms. Lasers that emit such waveforms are referred to as optical frequency combs. Here we show that Bloch gain - a nonclassical phenomenon that was first predicted in the 1930s - plays an essential role in comb formation in quantum cascade lasers (QCLs). We develop a self-consistent theoretical model that considers all aspects of comb formation: bandstructure, electron transport, and cavity dynamics. It reveals that Bloch gain gives rise to a giant Kerr nonlinearity and serves as the physical origin of the linewidth enhancement factor in QCLs. Using a master equation approach, we explain how frequency modulated combs can be produced in Fabry-P\'{e}rot QCLs over the entire bias range. In ring resonators, Bloch gain triggers phase turbulence and the formation of soliton-like patterns.

Bloch and Zener predicted charge oscillations in a periodic potential under an applied constant electric field in the 1930s 1,2 , a phenomenon which is commonly referred to as Bloch oscillations. It attracted researchers ever since due to the property of oscillating charges to couple with electromagnetic waves, potentially offering new sources of radiation 3 . In condensed-matter theory, the motion of electrons in a periodic crystal lattice is governed by the energy-momentum relation within a Brillouin zone 4 . A constant electric field accelerates the electrons towards the edge of the Brillouin zone, where they experience Bragg reflection, resulting in an oscillatory motion. The width of the Brillouin zone in bulk crystals is large and thus electrons cannot reach the edge before they scatter. However, in semiconductor superlattices (Fig. 1a), made of alternate semiconductor layers 5 , it is significantly narrower and electrons can complete multiple oscillation cycles within their lifetime 6,7 . Ktitorov et al. 8 predicted tunable optical Bloch gain arising from these oscillations, which was subsequently verified in a GaAs/AlGaAs superlattice 9 . The gain is present even without population inversion, a necessary ingredient in the analysis of a classical harmonic oscillator (Fig. 1b). Moreover, the Bloch gain possesses an S-shaped profile (Fig. 1a), referred to as the dispersive gain 9 . This unique spectral response, sharply contrasted with the well-known symmetric Lorentzian gain of a harmonic oscillator (Fig. 1b), serves as the fingerprint feature of the Bloch gain. Analogous observations of this ubiquitous phenomenon are reported in other physical systems * nikola.opacak@tuwien.ac.at, † benedikt.schwarz@tuwien.ac.at as well, e.g. Josephson junctions 10 , Bose-Einstein condensates 11 , complex potentials with PT symmetry 12 and in optical 13 and acoustic 14 waves.
More recently, pure quantum-mechanical treatments of the Bloch gain were developed in the density matrix formalism 15 and the Green's function formalism 16 . They generalized the concept of the Bloch gain and showed that it is not exclusive to superlattices, but appears also between any two states (subbands) in semiconductor heterostructures, such as quantum cascade lasers (QCLs). QCLs are unipolar laser sources 17 , which emit in the midinfrared 18 and terahertz 19 spectral regions by nanoscale engineering of the conduction-band profile (Fig. 1c). The gain bandwidth of QCLs is broadened by elastic scattering processes beyond its natural limit defined with the carrier lifetimes 20,21 . An accompanying effect of these processes, neglected by most researchers so far, is the occurrence of scattering-assisted optical transitions between subbands 15,22 . They connect electronic states with nonidentical wavevectors and give rise to dispersive Bloch gain. The total gain is comprised of the Bloch contribution and the usual Lorentzian gain generated by the harmonic oscillator (Fig. 1c).
In this work, we conduct a rigorous theoretical and numerical study of the Bloch gain and its influence on the laser dynamics. A meticulous simulation tool is developed which models and self-consistently couples every aspect of QCL operation -from electronic band structure and charge transport to the light spatio-temporal evolution within the laser cavity. We show that a dominant Bloch gain contribution is present in any operating QCL and causes a giant Kerr nonlinearity at the laser wavelength. The induced nonlinearity plays an essential role in the laser cavity dynamics as it is a requirement for self-starting optical frequency combs 23 . Bloch gain is not only the reason why frequency modulated (FM) comb formation is predominantly found in dispersion compensated cavities 24 , but it also allows tuning the laser into the phase turbulence regime. This can trigger the generation of soliton-like structures 25,26 , establishing a bridge between semiconductor QCL lasers and Kerr microresonators 27 .
The spectral response of the laser active region is fully captured by its complex susceptibility χ = χ R + iχ I . The optical gain is defined as g = ωχ I /n r c with ω being the frequency, n r the refractive index and c the speed of light. The susceptibility that arises from any two subbands u and their complex susceptibilities χ = χ R + χ I . Optical gain, which is proportional to χ I , has a symmetric Lorentzian shape and χ R has a dispersive shape in the case of the harmonic oscillator. For a Bloch oscillator, the shape profiles of χ R and χ I are exchanged compared to the harmonic oscillator, due to a π/2 phase shift in the spectral response. This results in the dispersive shape of the Bloch gain. (c) Schematic of the QCL bandstructure with the laser levels and the optical transition. The complex susceptibility can be represented as a sum of both Bloch and harmonic contributions.
l in a semiconductor heterostructure is calculated as 15 : The total susceptibility in Eq. (1) comprises two components. The usual harmonic contribution is given by the first term in the square brackets in Eq. (1). It depends on the population inversion f u (k) − f l (k) and yields a Lorentzian gainshape. The dipole matrix element is µ ul , ε 0 is the vacuum permittivity, f (k) and γ(k) are the electron distribution and broadening at wavevector k and ∆W (k) = W u (k) − W l (k) is the resonant transition energy, where ∆W 0 = ∆W (k = 0) = ω 0 . The highlighted second term in Eq. (1) is more intriguing. It generates the Bloch gain by allowing optical transitions between states with different wavevectors. Introduced notations are level broadenings γ u,l , where γ = γ u + γ l 21 , and in-plane momenta of the final states, defined as . A thorough analysis is given in the Supplementary section 1.
While Eq. (1) provides an exact treatment of the Bloch gain, the origin of the dispersive spectral shape is not well understood. It is not rare in physics to opt for a simpler model that provides an intuitive understanding over the exact one. Bearing this in mind, Eq. (1) significantly reduces its complexity by assuming subband electron distributions in the frame of Boltzmann statistics. The electron concentration is usually low enough so that the Fermi-Dirac distribution reduces to the Boltzmann distribution and the carrier-carrier interaction is sufficiently large to enforce carrier thermalization 28 . Following the derivation presented in the Supplementary section 1.1, we analyti-cally obtain a simplified definition of χ: where L p is the QCL period length, k B is the Boltzmann constant, T is the temperature and n u,l are the electron sheet densities of subbands u, l. Eq. (2) provides an understanding of the origin of Bloch gain, which is proportional to the highlighted terms. Contrary to the harmonic susceptibility, it is not dependent on the population inversion (n u − n l ) but rather on the population sum (n u + n l ). The dispersive gainshape appears due to the imaginary value of the highlighted terms. They induce a π/2 phase shift and exchange the shapes of χ R and χ I (Fig. 1a & b). The factor b in Eq. (2) captures the impact of the Bloch gain. It deviates the total gainshape from a Lorentzian curve and causes spectral asymmetry. Subband nonparabolicity, which is known to induce similar behavior, has a weaker effect. Most importantly, Eq.
(2) allows straightforward implementation of Bloch gain in any carrier transport model, unlike Eq. (1), which requires k-space resolved approaches.
With the aim of quantitatively assessing the influence of the Bloch gain on the laser dynamics, our model is employed to a reference QCL device 29 . For details about our band structure and charge transport model, see the Methods section. The calculated conduction-band profile with probability densities of the states and the electron density are shown in the Supplementary section 4.1. With the knowledge of the electron population, the calculation of the optical gain for the lasing transition follows from Eq. (1). Its unsaturated value is shown in Fig. 2a. The Bloch gain induces an asymmetric total gainshape and a redshift of the peak. However, the unsaturated gain asymmetry conveys only a fraction of what happens above the laser threshold. In the usual harmonic description, the emission of light depletes the population inversion until the gain saturates to the threshold value, while χ R remains zero at the gain peak (Fig. 2b). On the other hand, the Bloch gain is independent of the population inversion and thus remains mostly unaffected. As the harmonic gain fades away with stronger light intensity, the Bloch contribution prevails and results in an increasingly asymmetrical gain accompanied with a red-shift and non-zero χ R at the gain peak (Fig. 2c). Intriguingly, a negative global population inversion is required to completely diminish the total gain (Supplementary section 4.2). Fig. 2d shows the saturated gain for three different current densities J. The gain peak is blueshifted due to the quantum-confined Stark effect and, more importantly, the asymmetry increases towards the dispersive shape.
Gain asymmetry has historically been analyzed in the context of the laser linewidth broadening. It was treated with the empirical linewidth enhancement factor (LEF) 30 , defined as LEF = −( ∂χ R /∂N)/( ∂χ I /∂N), where N is the carrier population. In interband lasers, the gain asymmetry and LEF dominantly originate from the opposite curvature of the valence and conduction band. Since both laser levels in a QCL have similar curvatures, the LEF was expected to vanish. Interestingly, non-zero experimental values were obtained mostly between -0.5 and 1.5 31,32 . We explain this with the gain asymmetry in QCLs that is dominantly caused by the Bloch gain. Furthermore, the LEF is frequency dependent, which yields a large range of values shown in Fig. 2e. Elimination of the Bloch gain in QCLs yields a symmetric gain profile and vanishing LEF (Supplementary secion 4.2). Fig. 2f shows the sim-ulated light-current-voltage (LIV) characteristic with the lasing threshold at around J = 1.6 kA/cm 2 and rollover at J = 5.5 kA/cm 2 29 . The calculated values of the LEF and factor b at the gain peak for the entire range of the current density J from the LIV are shown in Fig. 2g. Although the population inversion is clamped, the population sum increases with the current in Eq. (2). This leads to a linear dependence of the LEF and factor b on J, which matches observations found in literature 31,32,33 . The impact of the gain saturation is underlined yet again, as the saturated values notably break off from the unsaturated ones.
The gain asymmetry causes changes of both the gain and the refractive index of the active region. Their alterations are induced by variable electron population, as was described by Agrawal 34 . This is closely related to a dependence of the gain and refractive index on the intensity (Fig. 2c), which gives rise to a Kerr nonlinearity. Although the bulk nonlinearity of a semiconductor crystal is small, the resonant contribution from the asymmetric nature of the gain yields a giant Kerr nonlinearity due to ultrafast dynamics in QCLs 23 . Based on the saturation analysis of χ in Supplementary section 3, we calculate the resonant Kerr contribution due to Bloch gain to be in the range of 10 −15 m 2 /W, which is two orders of magnitude larger than the highest bulk values 35 .
Optical nonlinearities couple the amplitude and the phase of the intracavity laser field and give rise to coherent processes such as frequency comb formation. Frequency combs are lasers whose spectra consist of equidis- tant modes with a fixed phase relation 36,37 . Although historically their formation relied on the emission of short pulses 38 , recently a new type of frequency modulated (FM) combs is blossoming. They are self-starting and appear in numerous Fabry-Pérot laser types such as QCLs 39 , interband cascade lasers 40 , quantum dot lasers 41 and laser diodes 42 . The fascinating property of FM combs, which distinguishes them from other frequency combs, is an almost constant intensity accompanied with a linear frequency chirp 43 . This unique behavior was explained in 23 as a result of the group velocity dispersion (GVD) or more importantly a Kerr nonlinearity, thus bringing the role of the Bloch gain in FM comb formation to the foreground.
In order to quantitatively study the frequency comb dynamics, we conduct spatio-temporal simulations of the intracavity field based on a master equation approach 23 . We describe the gain shape asymmetry through the parameter b from Eq. (2). Its population dependence can be accurately modeled as a function of the current density J and the laser field intensity I (Supplementary equation (42)). This allows a self-consistent implementation into the master equation to include the Bloch gain: where E± are the right and left propagating field envelopes, T 1 , T 2 and T g are the recovery times of the gain, polarization and the population grating, α w is the waveguide loss, g is the saturated gain, I sat is the saturation intensity and I = |E + | 2 +|E − | 2 the normalized intensity. The Bloch gain enters the equation through terms b, ξ(b) and T 2 (b). A detailed derivation is presented in the Supplementary section 2.3, along with the analysis for interband lasers with slower dynamics. The numerical results for a Fabry-Pérot QCL are shown in Fig. 3. Using Eq. (3), we simulate 30 000 roundtrips of the electric field evolution to ensure that a steady state has been reached. Temporal evolution of the light intensity for one bias point is shown in Fig. 3a. The inclusion of the Bloch gain leads to a periodic waveform after 6000 roundtrips and the formation of an FM comb, which is fully characterized in the Supplementary section 4.2. Conversely, the intensity evolves chaotically in the absence of a locking mechanism provided by the Bloch gain induced Kerr nonlinearity. By extracting the scattering rates from the transport model, we are able to accurately simulate the intracavity dynamics from the laser threshold to rollover (Fig. 3b). The laser is in the single mode regime near the threshold and significantly broadens its spectrum with the current increase. The key role of the Bloch gain is clearly visible, as it leads to an FM comb operation over the entire bias range. This is indicated by the autocorrelation value equal to one in Fig. 3b. In sharp contrast, the pure harmonic gain results in unlocked states with the autocorrelation smaller than unity and chaotic spectra. This validates the Bloch gain induced giant Kerr nonlinearity as an efficient locking mechanism and explains why FM combs in QCLs have mostly been found in GVD compensated cavities 24,43,44 . The interplay with a non-zero GVD yields an unlocked state for most of the bias range (Supplementary section 4.4), in accordance with literature 44 .
Linking the physics of FM combs to Bloch gain induced giant Kerr nonlinearity suggests a connection to the Kerr combs in microresonators 27 . They represent passive me- Bloch gain leads to a multi-mode instability through phase turbulence 26 . A sech 2 envelope is fitted to the spectrum. Elimination of Bloch gain yields single-mode operation. (c) Temporal evolution of the intensity shows an initial turbulent regime that forms a frequency comb after 510 000 roundtrips. dia, where pumping is achieved through an external injection of a monochromatic laser and the gain stems from the Kerr nonlinearity of the bulk crystal. Through a cascaded parametric process, the injected wave induces the appearance of side-modes giving rise to phase-locked frequency combs in the form of temporal solitons 45 . QCL combs in ring cavities (Fig. 4a) have recently been shown to possess several similarities with Kerr microresonators 33,25 . Within the framework of the Ginzburg-Landau formalism 46 , it was demonstrated that a single mode operation is destabilized by the phase turbulence. The latter is controlled with the laser nonlinearity to induce multi-mode emission with a sech-type spectrum (Fig. 4b). This can trigger the formation of localized structures in the waveform (Fig. 4c), which are related to dissipative Kerr solitons. In the absence of the Bloch gain, the laser operates in a single-mode regime with constant intensity. The study in 47 demonstrated that the Kerr microresonators and ring QCLs can both be analyzed within the same theoretical framework and predicted the emission of temporal solitons from a ring QCL with a suitable nonlinearity. As we now know that the Kerr nonlinearity dominantly stems from the Bloch gain, by using our model to carefully optimize the gainshape for soliton emission, new methods of QCL comb formation with wider spectral bandwidth could be realized.
In conclusion, we have shown that a substantial Bloch gain contribution is present in QCLs due to scatteringassisted optical transitions and that it provides a coherent locking mechanism for frequency comb formation. A selfconsistent simulation model was built to study QCL operation altogether, including bands tructure, carrier transport and cavity field dynamics. It showed that the saturated gain considerably deviates from a Lorentzian curve towards a dispersive asymmetric shape due to Bloch gain.
We connect the gain asymmetry to the LEF and explain its experimentally obtained values. The nonclassical nature of the Bloch gain is captured with a single populationdependent parameter b. It is implemented in the master equation to study the spatio-temporal evolution of the laser field. In accordance with this, we discover that FM comb formation in Fabry-Pérot QCLs is triggered by a Bloch gain induced giant Kerr nonlinearity. The Bloch gain therefore acts as a efficient locking mechanism for the entire range of the bias current, which explains why FM combs were experimentally observed mostly in GVD compensated cavities. In a ring resonator, the impact of Bloch gain is particularly strong, due to the low cavity losses and the stronger saturation. The induced Kerr nonlinearity destabilizes the single mode operation through phase turbulence and can result in comb formation and the emission of localized structures. This paves the way towards broadband active Kerr combs in the mid-infrared range. By careful design of the laser active region and cavity, the Bloch contribution to the total gain can be controlled in order to tailor the induced LEF and Kerr nonlinearity. This would allow us to further optimize QCL frequency combs for broadband emission and discover new states of light.

Methods
Band structure and electron transport: The band structure of the device is calculated using an envelope function formalism in the two-band k · p model 48 . Electron transport between the states is computed by employing a one-dimensional density matrix approach 21 , which includes the energy-preserving tunneling 49,50 . The electron density is calculated self-consistently and accounts for the most relevant scattering processes via longitudinal optical phonons, acoustic phonons, interface roughness and alloy scattering 51,52 .
Laser cavity model: In order to study frequency comb dynamics, we use the master equation (3). We simulate 30 000 roundtrips of spatio-temporal laser field evolution in Fabry-Pérot cavities and 720 000 roundtrips in ring cavities to ensure that the laser is in a frequency comb regime. Such long simulation times are allowed by a highly efficient numeric implementation using the CUDA library 53 . The code is parallelized and runs on 2000 threads on the graphics processing unit (GPU) within a PC, which resulted in a speed-up factor of 500 compared to the implementation on a central processing unit. We used an NVIDIA GeForce GTX 1070 Ti GPU. As an example, the spatio-temporal simulation shown in Fig. 4c, which consists of 720 million time steps, took 32 minutes to run. The values of the parameters that are used in the cavity model are listed in the Supplementary section 5.