## Abstract

The generation of frequency combs in the mid-infrared and terahertz regimes from compact and potentially cheap sources could have a strong impact on spectroscopy, as many molecules have their rotovibrational bands in this spectral range. Thus, quantum cascade lasers (QCLs) are the perfect candidates for comb generation in these portions of the electromagnetic spectrum. Here we present a theoretical model based on a full numerical solution of Maxwell-Bloch equations suitable for the simulation of such devices. We show that our approach captures the intricate interplay between four wave mixing, spatial hole burning, coherent tunneling and chromatic dispersion which are present in free running QCLs. We investigate the premises for the generation of QCL based terahertz combs. The simulated comb spectrum is in good agreement with experiment, and also the observed temporal pulse switching between high and low frequency components is reproduced. Furthermore, non-comb operation resulting in a complex multimode dynamics is investigated.

© 2016 Optical Society of America

## 1. Introduction

Quantum cascade lasers (QCLs) promise to be efficient, cheap and compact generators of frequency combs in the mid- and far-infrared portions of the electromagnetic spectrum. QCL based combs in both spectral regions have been experimentally demonstrated [1–5], but their spectral coverage has been limited to a fraction of their central frequency. From an application point of view, frequency combs with bandwidths spanning an octave are highly desirable, since then the carrier offset frequency can be readily identified via a self-referencing scheme based on heterodyne detection [6]. However, achieving comb stability over such a broadband frequency range has proven to be difficult [7], namely due to the distorting effect of chromatic dispersion. Nevertheless, also more narrowband QCL combs are useful for practical applications, as demonstrated using the so-called dual-comb spectroscopy technique [8,9].

The experimentally demonstrated QCL frequency combs are based on free running lasers [1–3,7]. Here, high order nonlinear optical processes, in particular four-wave mixing (FWM), have been identified as the main mode proliferation mechanisms that contribute to comb formation [10, 11]. In contrast, it has been argued that group velocity dispersion (GVD) leads to unstable multimode operation and thus limits the full exploitation of the gain bandwidth of the material [12]. In the terahertz (THz) regime, the two widest comb generating devices demonstrated so far have shown a strong variation of the beatnote’s linewidth with changing injection current, indicating that comb operation comprises only a fraction of the whole dynamic range of these lasers [2,7].

This paper addresses fully time dependent simulations of QCL comb operation based on the semi-classical Maxwell-Bloch (MB) laser equations [13]. Based on numerical solutions of the MB equations, mode-locking in QCLs and the emergence of coherent optical instabilities have been analyzed [14–19]. In those works, **k**-space dependence has not been included as it significantly increases the computational load involved and limits the simulation time to only short intervals [20]. A direct numerical solution of the MB equations is even much more intricate for comb operation, since the optical field must be propagated over some 10^{4} cavity round trips to obtain converging results, which in turn poses strong requirements on the accuracy and numerical efficiency of the simulation approach. This problem had been circumvented by perturbatively solving the MB equations under the assumption of a discrete mode spectrum [11,12], confirming FWM as the main comb forming mechanism and GVD as a major detrimental effect. In this paper, we derive extended MB equations which also include coherent resonant tunneling between the injector and upper laser level, as required for a realistic modeling of QCL-based THz combs. These equations are coupled to self-consistent ensemble Monte Carlo (EMC) carrier transport simulations which provide the non-radiative transition rates between the energy levels, eliminating the need to use empirical electron lifetimes. Our simulation approach considers the full transient dynamics of the system, directly solving the extended MB equations without invoking any further assumptions. In this way, also regimes of non-comb or imperfect comb operation can be analyzed, including cases where spectrally separated sub-combs exist or only part of the spectral lines are phase-locked and participate in comb operation.

The paper is organized as follows: In Sec. 2 we introduce our theoretical model. This approach is used in Sec. 3 to investigate the active region of a longitudinal optical (LO) phonon depopulation THz QCL design used for frequency comb generation [2]. Here we focus on effects relevant for frequency comb operation, in particular the gain characteristics, GVD and FWM. Lastly, in Sec. 4 we present simulation results for comb operation and unstable behavior, respectively, and show that the theoretical results agree well with experimental data.

## 2. Theoretical model

We base our investigations on the Maxwell-Bloch (MB) laser equations, which present a suitable model for the description of light-matter interaction in microscopic systems. Specifically, here we use extended Bloch equations to describe the optical transition between the upper and lower laser level, additionally accounting for coherent resonant tunneling between the injector and upper laser level to obtain a realistic description of the investigated THz QCL design. We employ the standard rotating wave and slowly varying envelope approximations to reduce the numerical load of the solution [13–16]. The effect of nonradiative scattering mechanisms onto the carrier dynamics is captured phenomenologically via a rate equations approach [18, 21, 22]. The level eigenenergies, dipole moments and anticrossing strengths are obtained from a Schrödinger-Poisson solver [23]. The MB approach is coupled to a well-established ensemble Monte Carlo (EMC) carrier transport simulation code for QCLs [24–27], which provides the scattering rates. Thus, our approach is self-consistent and does not require empirical parameters, with the exception of pure dephasing rates.

First, in order to introduce some notation, let us take a simple toy model for a resonant phonon THz QCL, depicted in Fig. 1. In this configuration, we consider four relevant levels |1′〉, |3〉, |2〉 and |1〉 which are the injector level, the upper and lower laser levels and the depopulation level, which also serves as the injector level for the next period. The state |1〉 couples to |3〉 via the anticrossing energy ℏΩ_{1′3}, whereas the upper and lower laser levels interact via the instantaneous Rabi frequency Ω* _{L}*(

*t*) = −

*ez*

_{32}

*E*(

_{z}*t*) / ℏ. Here, $-e{z}_{32}=-e\u30083|\widehat{z}|2\u3009$ is the dipole matrix element,

*E*(

_{z}*t*) is the electric field component along the growth direction

*z*, and

*e*is the elementary charge. Lastly, we assume that the energy separation between |1′〉 and |3〉 is Δ

_{1′3}= ℏ

*ϵ*and that between the upper and lower laser levels Δ

_{32}= ℏ

*ω*

_{0}.

In THz QCLs, resonant tunneling plays an important role especially for thick barriers where the subbands involved in the electron transport have a small energy separation, which particularly applies to the injection barriers of various QCL designs [28–30]. This effect is usually treated in the tight-binding (TB) approximation [31], where the wave functions are calculated for a single isolated period, and the electron subbands in adjacent modules are simply obtained by adequate translation in energy and space by the module length *L _{p}* [28]. The energetic coupling between levels spanning the intramodule barrier is modelled by including the corresponding anticrossing energies in the Hamiltonian matrix [29,30]. In our case these are the injector and upper laser level, with the coupling energy ℏΩ

_{1′3}. In this tight-binding basis, the time evolution is governed by the von Neumann equation which, with phenomenologically included scattering rates, reads

In Eq. (1) the
${\rho}_{ij}=\u3008i|\widehat{\rho}|j\u3009$ denote the corresponding density matrix elements and we have set the zero energy at (*E*_{1′} + *E*_{3})/2. Furthermore
${\tau}_{ij}^{-1}$ denotes the outscattering rate from level *i* to level *j*, *τ _{i}* is the lifetime of level

*i*with ${\tau}_{i}^{-1}={\displaystyle {\sum}_{j\ne i}{\tau}_{ij}^{-1}}$, and

*i*and

*j*, containing lifetime broadening and a ”pure” dephasing time ${\tau}_{ij}^{pure}$ [28].

Notice that in Eq. (1) we have omitted the time evolution related to state |1〉, because it serves as the injector level for the next period, which allows us to eliminate it from the model by employing ”periodic” boundary conditions in the scattering rates matrix. It is vital for the simulation that these periodic boundary conditions are implemented correctly, because our model is formulated in a manner which does not phenomenologically include the injection current density *J* into the equations. Instead, we assume a periodic system where all carriers that reach the depopulation level |1〉 are immediately re-injected into the system through level |1′〉. In such a configuration the overall carrier density has to be conserved, which means that the relation ∑* _{j} dρ_{j j}/dt* = 0 has to be satisfied at all times. Then, the current injected into the system is essentially equal to the current entering into the depopulation level, |1〉. From rate equation analysis, the current density for a volume-averaged carrier concentration

*N*is given by [29]

The coupling of the microscopic density matrix equations to the macroscopic Maxwell’s equations is usually done via the incorporation of the polarization term in Maxwell’s equations as the expectation value of the quantum mechanical dipole moment operator,

*x*denotes the propagation direction,

*c*is the speed of light in vacuum,

*n*

_{0}denotes the refractive index of the material, and

*ϵ*

_{0}is the permittivity of free space. The full system of equations which we developed and solve in practice, is presented in Appendix A, and the numerical scheme is discussed in Appendix B.

## 3. Theoretical analysis of an LO phonon depopulation THz QCL with a strong injector anticrossing

In this section we extend our three level model, outlined in Sec. 2, to simulate the QCL design in [2]. This laser is based on a resonant LO phonon depopulation active region lasing at the central frequency *f*_{0} ≈ 3.5 THz, and has been shown to produce a stable comb spectrum containing more than 70 equidistant longitudinal modes in a free running regime of operation. At optimum bias, field autocorrelation measurements indicate multimode lasing, and radio frequency (RF) measurements produce a very strong and narrow beatnote with a minimum reported full width at half-maximum (FWHM) linewidth of approximately 1.53 kHz, indicating that at least part of the modes are phase-locked [1–3, 7]. Furthermore, a novel comb coherence detection technique, shifted-wave interference Fourier-transform (SWIFT) spectroscopy, has been developed to prove the stability of the THz comb over a large number of round trips [32]. To correctly describe the device under test, we first determine operating regimes where the ”single resonant tunneling and single optical transition” assumption, made in Sec. 2, is justified.

Figure 2(a) illustrates the calculated wave functions at a bias of 11 kV/cm obtained with the tight-binding approximation. From graphical inspection we see that there are in total five relevant levels per period, which we will refer to, according to their assumed role, as |*ULL*〉 for the upper laser level, |*LLL*1〉 for the higher energy level from a pair of lower laser levels, |*LLL*2〉 for the lower energy level from the same pair, |*DEP*1〉 for the higher energy level from a doublet of depopulation levels, and finally |*DEP*2〉 for the lowest energy level. Furthermore, since the structure is periodic, the depopulation levels from the previous period will be referred to as |*IN J*1〉 and |*IN J*2〉, respectively.

Figures 2(b) and 2(c) depict the calculated coupling strengths ℏΩ* _{i j}* between the injector states and the upper laser level for different biases, as well as the magnitudes of the dipole matrix elements, |

*z*|, between the upper laser level and the doublet of lower laser levels. The anticrossing (AC) energies were calculated via the method described in [31], and the numerical values were verified by diagonalization of the tight-binding Hamiltonian. Focusing on bias values at around 11 kV/cm, our calculations show that there is almost perfect energetic alignment between |

_{i j}*IN J*2〉 and |

*ULL*〉, whereas |

*IN J*1〉 and |

*ULL*〉 are separated by approximately Δ

_{INJ}_{1,}

*≈ 4.3 meV (not shown in the figure). Even though the anticrossing energies ℏΩ*

_{ULL}

_{INJ}_{1,}

*≈ 1.18 meV and ℏΩ*

_{ULL}

_{INJ}_{2,}

*≈ 1.38 meV are of comparable strength, the strong resonance condition between the pair |*

_{ULL}*IN J*2〉, |

*ULL*〉 enhances the tunneling probability, i.e. reduces the tunneling time, between those levels [33], as compared to the tunneling transition between |

*IN J*1〉 and |

*ULL*〉. This means that the majority of the tunneling electrons will prefer the |

*IN J*2〉 ↔ |

*ULL*〉 transport channel and hence our model, which includes only a single tunneling transition, ought to correctly capture the microscopic dynamics of the real device. From dipole moment calculations in Fig. 2(c), we see that the optical transition at 11 kV/cm is most likely to occur between the pair |

*ULL*〉, |

*LLL*1〉, and thus we can assign the state |2〉 in our equations to be subband |

*LLL*1〉 from the system under investigation. Similarly, due to the discussion above, we can set the injector level in our reduced model, i.e. |1′〉, to state |

*IN J*2〉. Finally, we can map |3〉 to subband |

*ULL*〉. The remaining states |

*IN J*1〉 and |

*LLL*2〉 are then treated within a rate equations approach and included into the scattering rates matrix by adequately extending Eq. (1).

#### 3.1. Gain and dispersion characterization

In order to extract the spectral gain profile as well as the strength of the chromatic dispersion induced by the active region design, we have performed simulations emulating the THz time-domain spectroscopy (THz-TDS) technique, often used for gain characterization of THz QCLs [34–36].

We apply our model to a ring cavity configuration with length *L* = 2.5 mm and neglect spatial hole burning effects which are not expected to play a role for this simulation. A weak unchirped Gaussian pulse is propagated inside the cavity for one round trip of length *L*. At each time step *t _{n}* of our simulation, as well as at different points

*x*along the cavity length, we record the electric field envelope ${f}_{j}^{n}$ for further data processing. From this data the real and imaginary parts of the refractive index can be calculated in a straightforward manner, as briefly discussed below.

_{j}In a ring cavity, we have only forward propagating waves (no standing waves) and thus the electric field can be written as *E*(*x, t*) = ℜ { *f* (*x, t*) exp [i(*k _{c} x* −

*ω*)] }, where

_{c}t*f*is the (slowly varying) envelope function and

*k*and

_{c}*ω*are as defined in Appendix A. We will denote the electric field of the injected seed pulse at the input facet of our cavity as

_{c}*E*(

_{in}*t*) and the field of the detected pulse as

*E*(

_{out}*t*). Their Fourier transforms are

*E*(

_{in}*ω*) and

*E*(

_{out}*ω*), with the angular frequency

*ω*. The corresponding envelope functions are

*f*(

_{in}*t*) and

*f*(

_{out}*t*), with the Fourier transforms

*F*(

_{in}*ω*) and

*F*(

_{out}*ω*), respectively. For a weak seed pulse, the light-matter interaction is linear, and the active region can be described by its complex refractive index

*n*(

*ω*) =

*n*′ (

*ω*) + i

*n*″ (

*ω*), with the amplitude gain coefficient given by

*g*(

*ω*) = −

*n*″ (

*ω*)

*ω/c*. The dependence between the seed and output fields is then given by

*E*(

_{out}*ω*) =

*E*(

_{in}*ω*) exp i

*ωnL/c*, which corresponds to

*F*(

_{out}*ω*−

*ω*) exp(i

_{c}*k*) =

_{c}L*F*(

_{in}*ω*−

*ω*) exp i

_{c}*ωnL/c*. From this, we obtain

Figure 3 shows the results from numerical THz-TDS simulations at a bias of 10.8 kV/cm for the active region in [2] for a single round trip of the seed THz pulse. All relevant simulation parameters, as well as the scattering times are given in Appendix C.

Since our aim was to probe the unsaturated gain profile, for these numerical experiments we used a weak Gaussian seed pulse with an instantaneous Rabi frequency amplitude of 0.5 ns^{−1} and FWHM duration of 0.707 ps, corresponding to a FWHM bandwidth of 623 GHz for a transform limited pulse. In Fig. 3 the blue curve illustrates the simulated spectral amplitude gain, obtained from Eq. (7). We clearly can observe a pronounced splitting of the gain spectra into two frequency lobes, one centred around 3.55 THz and another one around 4.21 THz. The red curve in Fig. 3 depicts the frequency resolved group velocity, calculated from *v _{g}* = [

*∂k*(

*ω*)

*/∂ω*]

^{−1}with

*k*(

*ω*) =

*n*′ (

*ω*)

*ω/c*, and normalized to the central frequency’s phase velocity

*c/n*

_{0}. Due to the strong resonances at 3.55 THz and 4.21 THz, the low and high frequency components are delayed with respect to each other as

*v*

_{g}n_{0}

*/c*approaches 0.9511 and 0.9326 for the low and high frequency gain peaks, respectively, where this ratio depends on the strength of the corresponding transition. The upper left inset of Fig. 3 depicts the calculated higher order phase Ψ =

*k*(

*ω*)

*L*(with the linear part removed) together with a third order polynomial fit to it, Ψ

*, and shows that the dispersion relation within the spectral range of interest (i.e. from 3.5 THz to 4.2 THz) is approximately cubic. Lastly, the top right inset illustrates the second derivative of the wave number with respect to frequency, i.e.*

_{fit}*k*″ (

*ω*), which is a measure of GVD.

From the above analysis we can conclude that even in the absence of bulk or waveguide dispersion, the doubly peaked resonant nature of the transition will cause significant dispersion in the cavity which is expected to deteriorate the comb performance. Here a question naturally arises: ”How does the presence of strong chromatic dispersion impact the mode proliferation process?” If the multimode behaviour of free running QCLs is due to FWM, as suggested in [10, 11], then one would intuitively expect that such a high GVD will violate the phase matching condition and thus yield this nonlinear process ineffective. However, from experiment [2, 7], multimode operation of both mid-infrared and THz QCLs could be observed, even in the presence of strong dispersion. To investigate further this question, we analyze the nature of this mode generation mechanism and estimate the amount of phase mismatch induced by the resonant transition.

#### 3.2. Four wave mixing

We have conducted numerical experiments, similar to the THz-TDs technique, where we pumped the 2.5 mm ring cavity laser, outlined above, with two frequencies *ω*_{1} and *ω*_{2} < *ω*_{1}. For 50 round trips we collected the signal and calculated the resulting power spectrum.

Figures 4(a)–4(c) illustrate the results from these simulations when the seed frequencies *ω*_{1} and *ω*_{2} were varied. In Fig. 4(a), we chose *ω*_{1} and *ω*_{2} to be separated by the free spectral range, Δ*ω* = 2*πc/* (*n*_{0}*L*), and to reside within the high frequency lobe of the spectrum. We can clearly observe the formation of side-modes separated by Δ*ω*, populating this whole frequency lobe. Similarly in Fig. 4(b), our seed frequencies were chosen to be separated by 2Δ*ω* and again side-modes can be seen in the spectrum, however this time with mode spacing of 2Δ*ω*. This kind of behaviour is a trademark of degenerate FWM described by susceptibilities *χ*(−*ω _{a}*;

*ω*

_{1}, −

*ω*

_{2},

*ω*

_{1}) and

*χ*(−

*ω*;

_{s}*ω*

_{2}, −

*ω*

_{1},

*ω*

_{2}) [37], where the pump modes

*ω*

_{1}and

*ω*

_{2}mix to produce a signal at the anti-Stokes and Stokes frequencies

*ω*and

_{a}*ω*, respectively. This process is schematically illustrated in Fig. 4(d). In Figs. 4(a) and 4(b),

_{s}*ω*

_{1,2}are distributed around 4.2 THz, and there is no activation of modes under the low-frequency lobe part of the spectrum. Such dynamics differs from simulations at the same bias, but within a Fabry-Perot cavity, where we observe lasing from both lobes of the gain, as will be discussed in the next section. This means that in order for the FWM process to start, some kind of seeding mechanism is necessary. Since in these experiments we considered a ring-cavity laser, multimode instabilities such as spatial hole burning [15] were not included, and lasing only started in the lobe excited by the seed. Lastly, Fig. 4(c) shows simulation results when both pump modes are chosen to be in resonance with the corresponding gain peaks, and again we observe the familiar formation of side modes due to FWM.

The phase mismatch for the anti-Stokes component of the degenerate FWM process, depicted in Fig. 4(a), is [37]

*L*, we have Δ

*ω*=

*πc*/(

*n*

_{0}

*L*), and Taylor-expanding around

*ω*

_{1}up to third order in Δ

*ω*yields, with

*ω*

_{2}=

*ω*

_{1}– Δ

*ω*and

*ω*= 2

_{a}*ω*1 −

*ω*2 =

*ω*1 + Δ

*ω*,

Plugging in typical values of *L* = 5 mm, *n*_{0} = 3.6 and
${\partial}_{\omega}^{2}k\approx 2{\text{ps}}^{2}/\text{mm}$ as an upper estimate for GVD, we obtain a phase mismatch of *L* |Δ*k*| = 0.0274 rad, which is negligible. This means that despite significant dispersion present in the cavity, under the favourable conditions of broadband gain, strong third order nonlinearity and some kind of multimode instability mechanism, such lasers can potentially emit a multitude of longitudinal modes even in a free running regime of operation, as reported in numerous experiments [1–3,7]. This, however, does not mean that such GVD is not strong enough to hamper the comb formation over the full spectral bandwidth. We will elaborate further in Sec. 4.1 and Sec. 4.2 on the detrimental effect of GVD onto the comb coherence, where we simulate devices with and without taking care of GVD compensation.

## 4. Simulation of comb operation and comparison to experiment

In this section, we present simulation results for comb operation of the device in [2], now considering counter-propagating waves in the Fabry-Perot cavity with a length *L* = 5 mm, giving rise to spatial hole burning. Our simulations only use the device specifications and well known material parameters as an input, with the exception of empirical dephasing times, which were adapted from [30]. Again, the electron wavefunctions and eigenenergies are computed with a Schrödinger-Poisson solver, and the carrier lifetimes are extracted from Monte Carlo carrier transport simulations. The used model parameters are listed in Appendix C. We propagate over ~15000 round trips to obtain results close to steady state.

#### 4.1. Frequency domain

In Fig. 5, we compare the calculated spectra and beatnotes from simulations without dispersion compensation, Figs. 5(c) and 5(d), and with dispersion compensation, Figs. 5(e) and 5(f), to experimental data, Figs. 5(a) and 5(b). Here, the simulation bias is set slightly below resonance, i.e., at 10.8 kV/cm. In both the dispersion compensated and uncompensated case, we observe reasonable agreement with the experimental optical power spectra, Fig. 5(a), however substantial differences arise when we consider the corresponding beatnotes.

Let us first discuss the simulation results in Figs. 5(c) and 5(d). While experimental data show a strong and narrow beatnote with an FWHM linewidth of approximately 0.55 MHz, Fig. 5(b), in Fig. 5(d) we observe a multi-beatnote signal distributed around 8.14 GHz. This deviates from the experimental value of around 6.8 GHz because the experimental implementation of GVD compensation in the cavity results in an increased effective cavity length. In the RF spectrum of Fig. 5(d), we can distinguish two peaks, one at 8.149 GHz and another at 8.151 GHz, both with a FWHM of approximately 3 MHz (Fourier transformation of 10^{4} round trips yields a frequency resolution of 0.83 MHz, where at least two frequency intervals are required to resolve the peak). From the group velocity delay plot in Fig. 3 as well as from the spectra in Fig. 5(c) we can deduce that since the high frequency lobe’s spectral components carry more power and have lower group velocity *v _{g}*, the stronger peak in the beatnote belongs to those frequency components, whereas the weaker one is associated with the lower frequency modes. This means that there ought to be two co-propagating pulses, albeit with a different

*v*, which interfere in the beatnote measurement to produce the chaotic RF spectrum in Fig. 5(d). The amount of noise in this non-comb regime of operation can be therefore traced mainly to the timing jitter induced by the difference in group velocities of those pulses. This will be elaborated further below, when we consider our results in the time domain.

_{g}The situation drastically changes when we employ dispersion compensation in our simulations. To cancel the accumulated higher order phase components, after each round trip we Fourier-transform the electric field envelope, subtract from the resulting phasor the phase Ψ, extracted analogously to the simulation in Fig. 3, and inverse Fourier-transform the result. The obtained spectral power density and RF spectra are plotted in Figs. 5(e) and 5(f). We see that such a procedure equilibrates the difference in the group velocities of different lasing modes, which results in a single, strong and narrow beatnote with a linewidth corresponding to our numerical frequency resolution, as expected from the experimental results for comb operation in Fig. 5(b).

#### 4.2. Time domain

In [32] it was shown that the strong injector anticrossing leads to a splitting of the emitted comb spectra in frequency domain and to an effect which was coined as ”temporal hole burning” in time domain. To compare our simulation results to experiment, we applied a low-pass and high-pass finite impulse response filter to the simulated electric field, in order to separate the low and high frequency lobe components of the signal, respectively. Since the experimental design is dispersion compensated, we first consider this case in our simulation, corresponding to the results in Figs. 5(e) and 5(f). The smoothed experimental and simulated instantaneous intensity are depicted in Figs. 6(a) and 6(b), respectively, with a smoothing length of 10 ps in experiment and 3 ps in simulation.

We note that when the population dynamics is extremely fast, memory effects become relevant, which can be taken into account by using a non-Markovian approach [38,39], however at the cost of considerably increased numerical complexity. Since the modeling of frequency comb operation requires simulations over many 1000 round trips as discussed above, we rely on our model Eqs. (13)–(15) (see Appendix A). This approach already yields decent agreement between experiment and simulation as both data demonstrate a pulse switching behaviour between the high and low frequency lobe components of the spectra. The high frequency lobe pulse (blue curve) has a longer time duration than the low lobe signal (red curve). We believe this results from the difference in relative strengths of the corresponding spectral lobes, as the high lobe components experience more gain, as shown in Fig. 3. Furthermore, in the experimental trace an additional oscillatory substructure emerges, the reason for which remains to be clarified. However, a possible explanation could be the presence of windowing effects of the Fourier transform spectrometer.

The dependence of the relative strengths of both lobes on the bias has already been discussed in [30], however its implications on the time domain behaviour of the laser were not considered in greater detail. To summarize the conclusions in [30], the inclusion of a strong injector anticrossing leads to a splitting of the injector |1′〉 and upper laser level |3〉 into a doublet of so-called ”dressed states”. The higher energy level is denoted by |+〉 (also known as the ”anti-symmetric” state), and the other level by |−〉 (i.e. the ”symmetric” state). These states are separated by approximately the anticrossing energy 2ℏΩ_{1′3} (depending on the detuning), which can be readily confirmed from experiment [2]. The relative radiative coupling strength of states |+〉 and |−〉 with the lower laser level depends on the detuning from resonance. Assuming low intracavity intensity, it can be shown that below the resonant bias, i.e., for ℏ*ϵ* < 0, the high frequency lobe of the gain dominates the transition, whereas above resonance the low lobe does [30]. As our simulations and experimental data show, this results in a longer high frequency pulse.

Figure 6(c) displays the calculated density matrix elements *ρ*_{33}, *ρ*_{++}, *ρ*_{−−} and *ρ*_{1′1′} as a function of time. For all four terms, we observe dampened Rabi-oscillations occurring upon a pulse switching event, which manifest the intramode beating of the high and low lobe. The period of these oscillations is approximately 2.84 ps, giving us a frequency of 350 GHz, which is namely the beat frequency between the high and low lobes of the dispersion compensated laser’s spectrum (see Fig. 5). We can also clearly observe that whenever the high frequency lobe lases, state |+〉 gets more strongly depleted, whereas whenever the low lobe is switched on, state |−〉 saturates.

To illustrate the implications of the existence of co-propagating high and low frequency lobe pulses in the presence of GVD, the simulated instantaneous intensities with (left column) and without (right column) dispersion compensation are shown in Fig. 7 for different temporal intervals. Temporal hole burning, or pulse switching, is clearly present in both cases. With dispersion compensation, the periodicity is preserved over many round trips, as can be seen by comparing the intensity traces for the two displayed temporal intervals. By contrast, in the presence of GVD the high and low frequency lobes constantly compete with each other, which leads to a complex non-periodic temporal shape of the signal. Thus, while the signal from a low GVD QCL has a minimal timing jitter, which results in a narrow and strong beatnote, the one from a dispersion uncompensated device operates in a multi-pulse regime with a complicated temporal and spectral profile, and thus never quite reaches a steady state. From this it follows that coherent comb cannot be formed over the full spectrum, but either no comb or sub-combs can be observed [3,7].

One can naturally extrapolate these conclusions to a situation where the effects of material and waveguide dispersion add to a much more complicated dispersion profile, as compared to the one in Fig. 3. In such a case, the competition between electric field signals propagating with various different group velocities will lead to a chaotic and unstable beatnote, characteristic for multimode non-comb behaviour of QCLs [3,7].

Lastly, we would like to add that in our simulations, the temporal hole burning effect could be reproduced at different values of the applied bias ranging between 10.7 and 11.4 kV/cm (not shown here). For lower or alternatively higher biases, lasing occurred only in the dominant frequency lobe, due to the large difference between the corresponding components of the gain, and thus no pulse switching could be observed. Also we have found that the temporal overlap increases as the spectral gap between the high and low lobes is decreased, i.e. with decreasing anticrossing energy ℏΩ_{1′3}, as well as with broadening of the radiative transition’s linewidth.

## 5. Conclusion

We have presented a theoretical model suitable for the simulation of THz QCLs in comb and non-comb regimes of operation, based on the full numerical solution of the Maxwell-Bloch equations in rotating wave approximation. We have shown that our approach correctly captures the complicated dynamics between four wave mixing, coherent tunneling and spectral gain splitting, group velocity dispersion, as well as spatial hole burning, as it delivers simulation results in good agreement with experimental data. We have used this model to characterize the spectral gain profile and the corresponding group velocity delay, and to investigate the FWM as the dominant comb proliferation mechanism. Simulations of comb operation yielded a spectral comb shape in good agreement with experiment, and the temporal switching behaviour between the high and low frequency components reported in experiment could also be reproduced. Furthermore, we have demonstrated that in contrast to earlier theoretical approaches our model can also capture non-ideal comb operation and transitions between comb and non-comb operating regimes.

## A. Maxwell-Bloch equations for a Fabry-Perot resonator

Let us consider a Fabry-Perot resonator of a certain length *L*, and employ the rotating wave (RWA) and the slowly varying envelope approximations (SVEA) [13, 15]. The optical field inside the cavity can be written as a superposition of forward and backward propagating waves

*ω*and

_{c}*k*are the field’s carrier frequency and wave number, related by

_{c}*k*=

_{c}*n*

_{0}

*ω*. Since the superposition of two counter-propagating waves forms a standing wave, this will lead to the formation of an inversion grating along the propagation direction

_{c}/c*x*, a phenomenon also known as spatial hole burning. Thus for the diagonal elements of the density matrix we make the following ansatz

Notice that Eqs. (12a) and (12b) follow the electric field ansatz since the corresponding transition energies are in close resonance with the optical field. In contrast, the ansatz Eq. (12c) is chosen in analogy to Eq. (11) for the diagonal elements since the term *ρ*_{1′3} is expected to oscillate only slowly due to the small energetic spacing between subbands 1′ and 3.

Finally, plugging Eqs. (10)–(12) into Eqs. (1), (4) and (5) and invoking the rotating wave and slowly varying amplitude approximations, we obtain our model in its final form:

- Population densities$$\frac{d{\rho}_{{1}^{\prime}{1}^{\prime}}^{0}}{dt}=\mathrm{i}{\mathrm{\Omega}}_{{1}^{\prime}3}\left({\rho}_{{1}^{\prime}3}^{0}-{\rho}_{3{1}^{\prime}}^{0}\right)+\left(\frac{1}{{\tau}_{3{1}^{\prime}}}+\frac{1}{{\tau}_{31}}\right){\rho}_{33}^{0}+\left(\frac{1}{{\tau}_{2{1}^{\prime}}}+\frac{1}{{\tau}_{21}}\right){\rho}_{22}^{0}-\frac{{\rho}_{{1}^{\prime}{1}^{\prime}}^{0}}{{\tau}_{{1}^{\prime}}},$$$$\frac{d{\rho}_{{1}^{\prime}{1}^{\prime}}^{+}}{dt}=\mathrm{i}{\mathrm{\Omega}}_{{1}^{\prime}3}\left({\rho}_{{1}^{\prime}3}^{+}-{\rho}_{3{1}^{\prime}}^{+}\right)+\left(\frac{1}{{\tau}_{3{1}^{\prime}}}+\frac{1}{{\tau}_{31}}\right){\rho}_{33}^{+}+\left(\frac{1}{{\tau}_{2{1}^{\prime}}}+\frac{1}{{\tau}_{21}}\right){\rho}_{22}^{+}-\left(\frac{1}{{\tau}_{{1}^{\prime}}}+4{k}_{c}^{2}D\right){\rho}_{{1}^{\prime}{1}^{\prime}}^{+},$$$$\frac{d{\rho}_{33}^{0}}{dt}=\mathrm{i}{\mathrm{\Omega}}_{{1}^{\prime}3}\left({\rho}_{3{1}^{\prime}}^{0}-{\rho}_{{1}^{\prime}3}^{0}\right)+\mathrm{i}\frac{e{z}_{32}}{2\hslash}\left({f}_{-}^{*}{\eta}_{32}^{-}+{f}_{+}^{*}{\eta}_{32}^{+}-c.c\right)+\frac{1}{{\tau}_{{1}^{\prime}3}}{\rho}_{{1}^{\prime}{1}^{\prime}}^{0}+\frac{1}{{\tau}_{23}}{\rho}_{22}^{0}-\frac{{\rho}_{33}^{0}}{{\tau}_{3}},$$$$\frac{d{\rho}_{33}^{+}}{dt}=\mathrm{i}{\mathrm{\Omega}}_{{1}^{\prime}3}\left({\rho}_{3{1}^{\prime}}^{+}-{\rho}_{{1}^{\prime}3}^{+}\right)+\mathrm{i}\frac{e{z}_{32}}{2\hslash}\left[{f}_{-}^{*}{\eta}_{32}^{+}-{f}_{+}{\left({\eta}_{32}^{-}\right)}^{*}\right]+\frac{{\rho}_{{1}^{\prime}{1}^{\prime}}^{+}}{{\tau}_{{1}^{\prime}3}}+\frac{{\rho}_{22}^{+}}{{\tau}_{23}}-\left(\frac{1}{{\tau}_{3}}+4{k}_{c}^{2}D\right){\rho}_{33}^{+},$$$$\frac{d{\rho}_{22}^{0}}{dt}=-\mathrm{i}\frac{e{z}_{32}}{2\hslash}\left({f}_{-}^{*}{\eta}_{32}^{-}+{f}_{+}^{*}{\eta}_{32}^{+}-c.c\right)+\frac{1}{{\tau}_{{1}^{\prime}2}}{\rho}_{{1}^{\prime}{1}^{\prime}}^{0}+\frac{1}{{\tau}_{32}}{\rho}_{33}^{0}-\frac{{\rho}_{22}^{0}}{{\tau}_{21}},$$$$\frac{d{\rho}_{22}^{+}}{dt}=-\mathrm{i}\frac{e{z}_{32}}{2\hslash}\left[{f}_{-}^{*}{\eta}_{32}^{+}-{f}_{+}{\left({\eta}_{32}^{-}\right)}^{*}\right]+\frac{1}{{\tau}_{{1}^{\prime}2}}{\rho}_{{1}^{\prime}{1}^{\prime}}^{+}+\frac{1}{{\tau}_{32}}{\rho}_{33}^{+}-\left(\frac{1}{{\tau}_{2}}+4{k}_{c}^{2}D\right){\rho}_{22}^{+},$$
- Coherence terms$$\frac{d{\rho}_{{1}^{\prime}3}^{0}}{dt}=-\mathrm{i}\mathit{\u03f5}{\rho}_{{1}^{\prime}3}^{0}+\mathrm{i}{\mathrm{\Omega}}_{{1}^{\prime}3}\left({\rho}_{{1}^{\prime}{1}^{\prime}}^{0}-{\rho}_{33}^{0}\right)+\mathrm{i}\frac{e{z}_{32}}{2\hslash}\left({f}_{+}^{*}{\eta}_{{1}^{\prime}2}^{+}+{f}_{-}^{*}{\eta}_{{1}^{\prime}2}^{-}\right)-{\tau}_{||{1}^{\prime}3}^{-1}{\rho}_{{1}^{\prime}3}^{0},$$$$\frac{d{\rho}_{{1}^{\prime}3}^{\pm}}{dt}=-\mathrm{i}\mathit{\u03f5}{\rho}_{{1}^{\prime}3}^{\pm}+\mathrm{i}{\mathrm{\Omega}}_{{1}^{\prime}3}\left({\rho}_{{1}^{\prime}{1}^{\prime}}^{\pm}-{\rho}_{33}^{\pm}\right)+\mathrm{i}\frac{e{z}_{32}}{2\hslash}{f}_{\mp}^{*}{\eta}_{{1}^{\prime}2}^{\pm}-\left({\tau}_{||{1}^{\prime}3}^{-1}+4{k}_{c}^{2}D\right){\rho}_{{1}^{\prime}3}^{\pm},$$$$\frac{d{\eta}_{32}^{\pm}}{dt}=\mathrm{i}\left({\omega}_{c}-{\omega}_{0}\right){\eta}_{32}^{\pm}+\mathrm{i}\frac{e{z}_{32}}{2\hslash}\left[{f}_{\pm}\left({\rho}_{33}^{0}-{\rho}_{22}^{0}\right)+{f}_{\mp}\left({\rho}_{33}^{\pm}-{\rho}_{22}^{\pm}\right)\right]-\mathrm{i}{\mathrm{\Omega}}_{{1}^{\prime}3}{\eta}_{{1}^{\prime}2}^{\pm}-{\tau}_{||32}^{-1}{\eta}_{32}^{\pm},$$$$\frac{d{\eta}_{{1}^{\prime}2}^{\pm}}{dt}=\mathrm{i}\left({\omega}_{c}-{\omega}_{0}-\mathit{\u03f5}\right){\eta}_{{1}^{\prime}2}^{\pm}+\mathrm{i}\frac{e{z}_{32}}{2\hslash}\left({f}_{\pm}{\rho}_{{1}^{\prime}3}^{0}+{f}_{\mp}{\rho}_{{1}^{\prime}3}^{\pm}\right)-\mathrm{i}{\mathrm{\Omega}}_{{1}^{\prime}3}{\eta}_{32}^{\pm}-{\tau}_{||{1}^{\prime}2}^{-1}{\eta}_{{1}^{\prime}2}^{\pm}.$$

Notice that in Eq. (13) we have phenomenologically added a linear amplitude loss coefficient *l*_{0}, and in Eqs. (14b), (14d), (14f) and (15b) a diffusion term
$4{k}_{c}^{2}D$ describing the diffusion rate at which carriers diffuse away from peaks of the population grating. Here *D* denotes the diffusion constant, which for GaAs/AlGaAs systems is 46 cm^{2}/s [19, 40]. The fields *f*_{±} in Eq. (13) satisfy the boundary conditions *f*_{+} (0) = *r f*_{−} (0) and *r f*_{+} (*L*) = *f*_{−} (*L*) [14], where we set the reflection coefficient *r* = 1 and instead include mirror outcoupling into the loss coefficient *l*_{0}. In Appendix B, the numerical methods used to solve Eqs. (13)–(15) are briefly discussed.

## B. Numerics

Eqs. (13)–(15) comprise a system of two partial and thirteen ordinary differential equations, which we solve numerically. The field propagation equations, Eq. (13), are a pair of inhomogeneous hyperbolic equations the accurate numerical solution of which is far from trivial. From the area of computational fluid dynamics [41], it is known that a simple central differences discretization scheme for Eq. (13) will be highly unstable due to the introduction of strong *numerical* dispersion near sharp edges or discontinuities of the solution. A finite difference scheme that does not generate such spurious oscillations is called monotonicity preserving [41] and its usage is essential for the correct interpretation of simulation results, especially when one tries to quantify the amount of *physical* dispersion present. Without getting too much into detail, we present a second order linear finite difference discretization scheme, possessing the monotonicity preserving property (in the case when *η*_{±} ∝ *f*_{±}), for the model equation

We take an equidistant spatio-temporal grid with grid size Δ*x* and time step Δ*t* and set the values of the grid variables at spatial point *x _{m}* =

*m*Δ

*x*and time

*t*=

_{n}*n*Δ

*t*as

*f*

_{±}(

*m, n*), where the ± index denotes the forward/backward propagating envelopes, respectively, and

*η*

_{±}denotes the corresponding source terms. The time stepping scheme we use is based on a second order upwind discretization, first introduced by Risken and Nummedal [42], and is given by

*t*= Δ

*xn*

_{0}/

*c*, with

*n*

_{0}/

*c*being the velocity of light in the medium. The evaluation of the time derivative of

*η*

_{±}(

*x*,

*t*) can be computed analytically from the density matrix equations, Eq. (15). The terms in rectangular brackets, ${\left[\partial {\eta}_{\pm}/\partial x\right]}_{m}^{n}$ and ${\left[\partial {f}_{\pm}/\partial x\right]}_{m}^{n}$, can be computed via forward/backward finite differences, depending on the propagation direction of the field.

The density matrix equations, Eqs. (14) and (15), form a system of ordinary differential equations and the usage of any out of the box numerical solver ought to suffice. In our experience a suitable method which is high-order accurate and also preserves the normalization property of the density matrix, i.e. Tr(*ρ*) = 1, is given by the fifth order linear multi-step Adams-Bashforth method.

## C. Simulation parameters

In Table 1 we summarize the parameter set used for our simulations and in Table 2, the scattering rates between each pair of the active region subbands calculated with our ensemble Monte Carlo simulation code [21] is shown. The code includes all relevant scattering mechanisms reported to play a role in quantum cascade lasers, including longitudinal optical phonons, acoustic phonons, interface roughness and impurity scattering as well as electron-electron scattering. The calculated rates include all these mechanisms and are presented in units of ps^{−1}.

## Funding

This work was supported by the German Research Foundation (DFG) within the Heisenberg program (JI 115/4-1) and under DFG Grant No. JI 115/9-1.

## References and links

**1. **A. Hugi, G. Villares, S. Blaser, H. Liu, and J. Faist, “Mid-infrared frequency comb based on a quantum cascade laser,” Nature **492**, 229–233 (2012). [CrossRef]

**2. **D. Burghoff, T.-Y. Kao, N. Han, C. W. I. Chan, X. Cai, Y. Yang, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, “Terahertz laser frequency combs,” Nature Photon. **8**, 462–467 (2014). [CrossRef]

**3. **M. Wienold, B. Röben, L. Schrottke, and H. Grahn, “Evidence for frequency comb emission from a Fabry-Pérot terahertz quantum-cascade laser,” Opt. Express **22**, 30410–30424 (2014). [CrossRef]

**4. **Q. Lu, M. Razeghi, S. Slivken, N. Bandyopadhyay, Y. Bai, W. Zhou, M. Chen, D. Heydari, A. Haddadi, R. McClintock, M. Amanti, and C. Sirtori, “High power frequency comb based on mid-infrared quantum cascade laser at λ ~9 μ m,” Appl. Phys. Lett. **106**, 051105 (2015). [CrossRef]

**5. **F. Cappelli, G. Villares, S. Riedi, and J. Faist, “Intrinsic linewidth of quantum cascade laser frequency combs,” Optica **2**, 836–840 (2015). [CrossRef]

**6. **J. Ye, *Femtosecond Optical Frequency Comb: Principle, Operation and Applications* (Springer Science & Business Media, 2005). [CrossRef]

**7. **M. Rösch, G. Scalari, M. Beck, and J. Faist, “Octave-spanning semiconductor laser,” Nature Photon. **9**, 42–47 (2015). [CrossRef]

**8. **G. Villares, A. Hugi, S. Blaser, and J. Faist, “Dual-comb spectroscopy based on quantum-cascade-laser frequency combs,” Nat. Commun. **5**, 5192 (2014). [CrossRef]

**9. **Y. Yang, D. Burghoff, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, “Terahertz multiheterodyne spectroscopy using laser frequency combs,” Optica **3**, 499–502 (2016). [CrossRef]

**10. **P. Friedli, H. Sigg, B. Hinkov, A. Hugi, S. Riedi, M. Beck, and J. Faist, “Four-wave mixing in a quantum cascade laser amplifier,” Appl. Phys. Lett. **102**, 222104 (2013). [CrossRef]

**11. **J. Khurgin, Y. Dikmelik, A. Hugi, and J. Faist, “Coherent frequency combs produced by self frequency modulation in quantum cascade lasers,” Appl. Phys. Lett. **104**, 081118 (2014). [CrossRef]

**12. **G. Villares and J. Faist, “Quantum cascade laser combs: effects of modulation and dispersion,” Opt. Express **23**, 1651–1669 (2015). [CrossRef]

**13. **R. W. Boyd, *Nonlinear Optics* (Academic, 2003).

**14. **C. Y. Wang, L. Diehl, A. Gordon, C. Jirauschek, F. X. Kärtner, A. Belyanin, D. Bour, S. Corzine, G. Höfler, M. Troccoli, J. Faist, and F. Capasso, “Coherent instabilities in a semiconductor laser with fast gain recovery,” Phys. Rev. A **75**, 031802 (2007). [CrossRef]

**15. **A. Gordon, C. Y. Wang, L. Diehl, F. X. Kärtner, A. Belyanin, D. Bour, S. Corzine, G. Höfler, H. C. Liu, H. Schneider, T. Maier, M. Troccoli, J. Faist, and F. Capasso, “Multimode regimes in quantum cascade lasers: From coherent instabilities to spatial hole burning,” Phys. Rev. A **77**, 053804 (2008). [CrossRef]

**16. **V.-M. Gkortsas, C. Wang, L. Kuznetsova, L. Diehl, A. Gordon, C. Jirauschek, M. Belkin, A. Belyanin, F. Capasso, and F. Kärtner, “Dynamics of actively mode-locked quantum cascade lasers,” Opt. Express **18**, 13616–13630 (2010). [CrossRef]

**17. **M. A. Talukder and C. R. Menyuk, “Self-induced transparency modelocking of quantum cascade lasers in the presence of saturable nonlinearity and group velocity dispersion,” Opt. Express **18**, 5639–5653 (2010). [CrossRef]

**18. **Y. Wang and A. Belyanin, “Active mode-locking of mid-infrared quantum cascade lasers with short gain recovery time,” Opt. Express **23**, 4173–4185 (2015). [CrossRef]

**19. **N. Vukovic, J. Radovanovic, V. Milanovic, and D. Boiko, “Multimode RNGH instabilities of Fabry-Pérot cavity QCLs: impact of diffusion,” Opt. Quant. Electron. **48**, 1–10 (2016). [CrossRef]

**20. **C. Weber, F. Banit, S. Butscher, A. Knorr, and A. Wacker, “Theory of the ultrafast nonlinear response of terahertz quantum cascade laser structures,” Appl. Phys. Lett. **89**, 091112 (2006). [CrossRef]

**21. **C. Jirauschek and T. Kubis, “Modeling techniques for quantum cascade lasers,” Appl. Phys. Rev. **1**, 011307 (2014). [CrossRef]

**22. **R. C. Iotti and F. Rossi, “Microscopic theory of semiconductor-based optoelectronic devices,” Rep. Prog. Phys. **68**, 2533 (2005). [CrossRef]

**23. **C. Jirauschek, “Accuracy of transfer matrix approaches for solving the effective mass Schrödinger equation,” IEEE J. Quantum Electron. **45**, 1059–1067 (2009). [CrossRef]

**24. **C. Jirauschek, G. Scarpa, P. Lugli, M. S. Vitiello, and G. Scamarcio, “Comparative analysis of resonant phonon THz quantum cascade lasers,” J. Appl. Phys. **101**, 086109 (2007). [CrossRef]

**25. **C. Jirauschek and P. Lugli, “Monte-Carlo-based spectral gain analysis for terahertz quantum cascade lasers,” J. Appl. Phys. **105**, 123102 (2009). [CrossRef]

**26. **C. Jirauschek, “Monte Carlo study of intrinsic linewidths in terahertz quantum cascade lasers,” Opt. Express **18**, 25922–25927 (2010). [CrossRef]

**27. **C. Jirauschek, “Monte Carlo study of carrier-light coupling in terahertz quantum cascade lasers,” Appl. Phys. Lett. **96**, 011103 (2010). [CrossRef]

**28. **H. Callebaut and Q. Hu, “Importance of coherence for electron transport in terahertz quantum cascade lasers,” J. Appl. Phys. **98**, 104505 (2005). [CrossRef]

**29. **S. Kumar and Q. Hu, “Coherence of resonant-tunneling transport in terahertz quantum-cascade lasers,” Phys. Rev. B **80**, 245316 (2009). [CrossRef]

**30. **E. Dupont, S. Fathololoumi, and H. Liu, “Simplified density-matrix model applied to three-well terahertz quantum cascade lasers,” Phys. Rev. B **81**, 205311 (2010). [CrossRef]

**31. **G. Bastard, *Wave Mechanics Applied to Semiconductor Heterostructures* (John Wiley & Sons, 1990).

**32. **D. Burghoff, Y. Yang, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, “Evaluating the coherence and time-domain profile of quantum cascade laser frequency combs,” Opt. Express **23**, 1190–1202 (2015). [CrossRef]

**33. **B. S. Williams, “Terahertz quantum-cascade lasers,” Ph.D. thesis, Massachusetts Institute of Technology (2003).

**34. **D. P. Burghoff, “Broadband terahertz photonics,” Ph.D. thesis, Massachusetts Institute of Technology (2014).

**35. **N. Jukam, S. Dhillon, Z. Y. Zhao, G. Duerr, J. Armijo, N. Sirmons, S. Hameau, S. Barbieri, P. Filloux, C. Sirtori, X. Marcadet, and J. Tignon, “Gain measurements of THz quantum cascade lasers using THz time-domain spectroscopy,” IEEE J. Sel. Top. Quantum Electron. **14**, 436–442 (2008). [CrossRef]

**36. **M. Martl, J. Darmo, C. Deutsch, M. Brandstetter, A. M. Andrews, P. Klang, G. Strasser, and K. Unterrainer, “Gain and losses in THz quantum cascade laser with metal-metal waveguide,” Opt. Express **19**, 733–738 (2011). [CrossRef]

**37. **P. N. Butcher and D. Cotter, *The Elements of Nonlinear Optics* (Cambridge University, 1991).

**38. **S. Butscher, J. Förstner, I. Waldmüller, and A. Knorr, “Ultrafast electron-phonon interaction of intersubband transitions: Quantum kinetics from adiabatic following to Rabi-oscillations,” Phys. Rev. B **72**, 045314 (2005). [CrossRef]

**39. **I. Knezevic and B. Novakovic, “Time-dependent transport in open systems based on quantum master equations,” J. Comput. Electron. **12**, 363–374 (2013). [CrossRef]

**40. **C. Y. Wang, L. Kuznetsova, V. M. Gkortsas, L. Diehl, F. X. Kärtner, M. A. Belkin, A. Belyanin, X. Li, D. Ham, H. Schneider, P. Grant, C. Y. Song, S. Haffouz, Z. R. Wasilewski, H. C. Liu, and F. Capasso, “Mode-locked pulses from mid-infrared quantum cascade lasers,” Opt. Express **17**, 12929–12943 (2009). [CrossRef]

**41. **P. Wesseling, *Principles of Computational Fluid Dynamics* (Springer Science & Business Media, 2009).

**42. **H. Risken and K. Nummedal, “Self-pulsing in lasers,” J. Appl. Phys. **39**, 4662–4672 (1968). [CrossRef]

**43. **S. Kohen, B. S. Williams, and Q. Hu, “Electromagnetic modeling of terahertz quantum cascade laser waveguides and resonators,” J. Appl. Phys. **97**, 053106 (2005). [CrossRef]