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⁴ 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₃₂E_z (t) / ℏ. Here, $-e{z}_{32}=-e〈3|\widehat{z}|2〉$ 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 Δ₃₂ = ℏω₀.

 

Fig. 1 Schematic diagram of a simple three-level LO-phonon depopulation THz QCL, where the upper laser level is populated via resonant tunneling.

Download Full Size | PDF

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}=〈i|\widehat{\rho }|j〉$ denote the corresponding density matrix elements and we have set the zero energy at (E_1′ + E₃)/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

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,

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₀ ≈ 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, |LLL1〉 for the higher energy level from a pair of lower laser levels, |LLL2〉 for the lower energy level from the same pair, |DEP1〉 for the higher energy level from a doublet of depopulation levels, and finally |DEP2〉 for the lowest energy level. Furthermore, since the structure is periodic, the depopulation levels from the previous period will be referred to as |IN J1〉 and |IN J2〉, respectively.

 

Fig. 2 (a) The moduli squared of the wave functions of the THz QCL in [2] within the tight-binding approximation. (b), The anticrossing coupling strengths between the pair of injector levels, |IN J1〉, |IN J2〉, and the upper laser level |ULL〉, computed via the method outlined in [31]. (c), Dipole elements calculated for the same laser, obtained with the extended basis Hamiltonian, Ĥ_ext.

Download Full Size | PDF

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_{i j}|, 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 |IN J2〉 and |ULL〉, whereas |IN J1〉 and |ULL〉 are separated by approximately Δ_INJ1,ULL ≈ 4.3 meV (not shown in the figure). Even though the anticrossing energies ℏΩ_INJ1,ULL ≈ 1.18 meV and ℏΩ_INJ2,ULL ≈ 1.38 meV are of comparable strength, the strong resonance condition between the pair |IN J2〉, |ULL〉 enhances the tunneling probability, i.e. reduces the tunneling time, between those levels [33], as compared to the tunneling transition between |IN J1〉 and |ULL〉. This means that the majority of the tunneling electrons will prefer the |IN J2〉 ↔ |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〉, |LLL1〉, and thus we can assign the state |2〉 in our equations to be subband |LLL1〉 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 J2〉. Finally, we can map |3〉 to subband |ULL〉. The remaining states |IN J1〉 and |LLL2〉 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_j 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.

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 − ω_ct)] }, where f is the (slowly varying) envelope function and k_c 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 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′ (ω) + in″ (ω), 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 (ω − ω_c) exp(ik_c L) = F_in (ω − ω_c) exp iω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.

 

Fig. 3 Simulated spectral gain profile (blue curve, left y-axis) together with the normalized group velocity v_gn₀/c (red curve, right y-axis). (Left inset) The higher order phase acquired by the test pulse. (Right inset) The second derivative of the wave number with respect to the angular frequency.

Download Full Size | PDF

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⁻¹ 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(ω)/∂ω]⁻¹ with k(ω) = n′ (ω)ω/c, and normalized to the central frequency’s phase velocity c/n₀. 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_gn₀/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, Ψ_fit, 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. 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 ω₁ and ω₂ < ω₁. 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 ω₁ and ω₂ were varied. In Fig. 4(a), we chose ω₁ and ω₂ to be separated by the free spectral range, Δω = 2πc/ (n₀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; ω₁, −ω₂, ω₁) and χ(−ω_s; ω₂, −ω₁, ω₂) [37], where the pump modes ω₁ and ω₂ mix to produce a signal at the anti-Stokes and Stokes frequencies ω_a and ω_s, respectively. This process is schematically illustrated in Fig. 4(d). In Figs. 4(a) and 4(b), ω_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.

 

Fig. 4 (a) Optical spectrum obtained from THz-TDs simulations when the seed frequencies, ω₁ and ω₂, are separated by the free spectral range Δω. (b) Same as (a), however this time ω₁ − ω₂ = 2Δω. (c) Optical spectrum from simulations where the seed frequencies are distributed under both peaks of the spectral gain. (d) Schematic representation of degenerate FWM [37], where two pump modes ω₁,ω₂ combine to produce sidebands at frequencies ω_a = 2ω₁ − ω₂ and ω_s = 2ω₂ − ω₁.

Download Full Size | PDF

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

Plugging in typical values of L = 5 mm, n₀ = 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.

 

Fig. 5 Optical power spectra (left column) and beatnotes (right column) from experiment and simulations of the device in [2]. (a), (b) Experimental data for the dispersion compensated laser, at driving current of 0.9 A. Simulation results without dispersion compensation (c)–(d) and with dispersion compensation (e)–(f). The experimentally detected beatnote in (b) has a FWHM of approx. 0.553MHz, whereas the simulated beatnote in (f) has a resolution limited FWHM of 1.66MHz. The strongest peak in (d) has a FWHM of 3.04 MHz.

Download Full Size | PDF

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⁴ 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_g, 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.

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.

 

Fig. 6 Time domain separation of the optical field into high and low frequency lobe components. (a) Measured intensity over time from [32] for the same THz comb device as in [2], driven with 0.9 A current. (b) Simulated intensity over time for a dispersion compensated QCL evolved for ~15000 round trips. (c) Time dependence of the electron populations in the injector and upper laser level (i.e. ρ_1′1′, ρ₃₃) as well as the dressed states (i.e. ρ₊₊, ρ₋₋).

Download Full Size | PDF

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 ρ₃₃, ρ₊₊, ρ₋₋ 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].

 

Fig. 7 Simulated high (blue) and low (red) frequency pulses for a numerically dispersion compensated (left column) and an uncompensated (right column) THz QCL. Data from round trip 4615-4618 (top) and 4710-4713 (bottom) are shown.

Download Full Size | PDF

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

(10)$$E_z\left(x,t\right)=\frac{1}{2}\left\{f_{+}\left(x,t\right)\exp \left[\mathrm{i}\left(k_{c}x-{\omega }_{c}t\right)\right]+f_{-}\left(x,t\right)\exp \left[-\mathrm{i}\left(k_{c}x+{\omega }_{c}t\right)\right]+c.c.\right\},$$

where ”c.c.” denotes the complex conjugate of the preceding expression, the + and − signs specify the forward/backward propagating waves’ envelopes, respectively, and ω_c and k_c are the field’s carrier frequency and wave number, related by k_c = n₀ω_c/c. 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 x, a phenomenon also known as spatial hole burning. Thus for the diagonal elements of the density matrix we make the following ansatz

(11)$${\rho }_{ii}\left(x,t\right)={\rho }_{ii}^0\left(x,t\right)+{\rho }_{ii}^{++}\left(x,t\right)\exp \left(2\mathrm{i}{k}_{c}x\right)+{\rho }_{ii}^{-}\left(x,t\right)\exp \left(-2\mathrm{i}{k}_{c}x\right),$$

where ${\rho }_{ii}^{+}={\left({\rho }_{ii}^{-}\right)}^{*}$ are the inversion grating’s amplitudes [14]. Lastly, we decompose the coherences of the density matrix as

(12a)$${\rho }_{32}\left(x,t\right)={\eta }_{32}^{+}\left(x,t\right)\exp \left[\mathrm{i}\left(k_{c}x-{\omega }_{c}t\right)\right]+{\eta }_{32}^{-}\left(x,t\right)\exp \left[-\mathrm{i}\left(k_{c}x+{\omega }_{c}t\right)\right],$$

(12b)$${\rho }_{1^{\prime }2}\left(x,t\right)={\eta }_{1^{\prime }2}^{+}\left(x,t\right)\exp \left[\mathrm{i}\left(k_{c}x-{\omega }_{c}t\right)\right]+{\eta }_{1^{\prime }2}^{-}\left(x,t\right)\exp \left[-\mathrm{i}\left(k_{c}x+{\omega }_{c}t\right)\right],$$

(12c)$${\rho }_{1^{\prime }3}\left(x,t\right)={\rho }_{1^{\prime }3}^0\left(x,t\right)+{\rho }_{1^{\prime }3}^{+}\left(x,t\right)\exp \left(2\mathrm{i}{k}_{c}x\right)+{\rho }_{1^{\prime }3}^{-}\left(x,t\right)\exp \left(-2\mathrm{i}{k}_{c}x\right).$$

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:

1. Propagation equations

   (13)$$\frac{n_0}{C}{\partial }_t{f}_{\pm }\pm {\partial }_x{f}_{\pm }=-\mathrm{i}\frac{N\mathrm{\Gamma }e{z}_{32}{k}_c}{{\mathit{ϵ}}_0{n}_0^2}{\eta }_{32}^{\pm }-l_0{f}_{\pm },$$
2. Population densities

   (14a)$$\frac{d{\rho }_{1^{\prime }{1}^{\prime }}^0}{dt}=\mathrm{i}{\mathrm{\Omega }}_{1^{\prime }3}\left({\rho }_{1^{\prime }3}^0-{\rho }_{31^{\prime }}^0\right)+\left(\frac{1}{{\tau }_{31^{\prime }}}+\frac{1}{{\tau }_{31}}\right){\rho }_{33}^0+\left(\frac{1}{{\tau }_{21^{\prime }}}+\frac{1}{{\tau }_{21}}\right){\rho }_{22}^0-\frac{{\rho }_{1^{\prime }{1}^{\prime }}^0}{{\tau }_{1^{\prime }}},$$

   (14b)$$\frac{d{\rho }_{1^{\prime }{1}^{\prime }}^{+}}{dt}=\mathrm{i}{\mathrm{\Omega }}_{1^{\prime }3}\left({\rho }_{1^{\prime }3}^{+}-{\rho }_{31^{\prime }}^{+}\right)+\left(\frac{1}{{\tau }_{31^{\prime }}}+\frac{1}{{\tau }_{31}}\right){\rho }_{33}^{+}+\left(\frac{1}{{\tau }_{21^{\prime }}}+\frac{1}{{\tau }_{21}}\right){\rho }_{22}^{+}-\left(\frac{1}{{\tau }_{1^{\prime }}}+4k_c^{2}D\right){\rho }_{1^{\prime }{1}^{\prime }}^{+},$$

   (14c)$$\frac{d{\rho }_{33}^0}{dt}=\mathrm{i}{\mathrm{\Omega }}_{1^{\prime }3}\left({\rho }_{31^{\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},$$

   (14d)$$\frac{d{\rho }_{33}^{+}}{dt}=\mathrm{i}{\mathrm{\Omega }}_{1^{\prime }3}\left({\rho }_{31^{\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}+4k_c^{2}D\right){\rho }_{33}^{+},$$

   (14e)$$\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}},$$

   (14f)$$\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}+4k_c^{2}D\right){\rho }_{22}^{+},$$
3. Coherence terms

   (15a)$$\frac{d{\rho }_{1^{\prime }3}^0}{dt}=-\mathrm{i}\mathit{ϵ}{\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,$$

   (15b)$$\frac{d{\rho }_{1^{\prime }3}^{\pm }}{dt}=-\mathrm{i}\mathit{ϵ}{\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}+4k_c^{2}D\right){\rho }_{1^{\prime }3}^{\pm },$$

   (15c)$$\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 },$$

   (15d)$$\frac{d{\eta }_{1^{\prime }2}^{\pm }}{dt}=\mathrm{i}\left({\omega }_c-{\omega }_0-\mathit{ϵ}\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₀, and in Eqs. (14b), (14d), (14f) and (15b) a diffusion term $4k_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²/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₀. 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

(16)$$\frac{\partial f_{\pm }}{\partial t}=\mp c\frac{\partial f_{\pm }}{\partial x}+{\eta }_{\pm }\left(x,t\right)+k{f}_{\pm }.$$

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

(17)$$\begin{array}{l}{f}_{\pm }\left(m,n+1\right)=f_{\pm }\left(m\mp 1,n\right)+\mathrm{\Delta }t\left[{\eta }_{\pm }\left(m,n\right)+k{f}_{\pm }\left(m,n\right)\right]\\ +\frac{\mathrm{\Delta }{t}^2}{2}\left\{{\left[\frac{\partial {\eta }_{\pm }}{\partial t}\right]}_m^n\mp c{\left[\frac{\partial {\eta }_{\pm }}{\partial x}\right]}_m^n\mp 2kc{\left[\frac{\partial f_{\pm }}{\partial x}\right]}_m^n+k{\eta }_{\pm }\left(m,n\right)+k^2{f}_{\pm }\left(m,n\right)\right\},\end{array}$$

where the time step is chosen as Δt = Δ xn₀/c, with n₀/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⁻¹.

Table 1. Simulation parameters for a THz QCL, modelled after the device in [2]. The modal overlap factor and the facet reflectivities are selected based on [43] for a metal-metal waveguide with thickness 10 μm and width of 20 μm.

View Table | View all tables in this article

Table 2. Total scattering rates between each pair of relevant subbands of the device in [2] for a bias of 10.8 kV/cm. The rates are presented in ps⁻¹.

View Table | View all tables in this article

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]