Dynamics of a broad-band Quantum Cascade Laser. From chaos to coherent dynamics and mode-locking

The dynamics of a multimode Quantum Cascade Laser, is studied in a model based on effective semiconductor Maxwell-Bloch equations, encompassing key features for the radiationmedium interaction such as an asymmetric, frequency dependent, gain and refractive index as well as the phase-amplitude coupling provided by the Henry factor. By considering the role of the free spectral range and Henry factor, we develop criteria suitable to identify the conditions which allow to destabilize, close to threshold, the traveling wave emitted by the laser and lead to chaotic or regular multimode dynamics. In the latter case our simulations show that the field oscillations are associated to self-confined structures which travel along the laser cavity, bridging mode-locking and solitary wave propagation. In addition, we show how a RF modulation of the bias current leads to active mode-locking yielding high-contrast, picosecond pulses. Our results compare well with recent experiments on broad-band THz-QCLs and may help understanding the conditions for the generation of ultrashort pulses and comb operation in Mid-IR and THz spectral regions


Introduction
The multimode dynamics of Mid-IR and THz Quantum Cascade Lasers (QCLs) became a focus of interest for the realization of pulsed regimes and frequency combs for a number of applications including time-resolved measurements, frequency mixing, high-precision spectroscopy etc. [1]- [12]. Such phenomena can be retraced to the capability of realizing coherent locking of modes in a multimode emission regime. In literature, different models have been proposed to theoretically interpret the experimental evidences [13]- [21]. A first class of models are essentially based on Maxwell-Bloch equations (MBEs) (in2 or 3-level medium approximation) [13,14,16,17,[19][20][21], where in some cases additional effects (e.g. cubic nonlinearity, saturable absorption, cavity dispersion) have been phenomenologically introduced to provide a proper interpretation of the experimental findings. A common understanding emerging thereof was the role of Spatial Hole Burning (SHB) (introduced to describe the effect of counter-propagating fields in a Fabry-Perot (FP) scheme of the laser cavity) in reducing the instability threshold leading to the mode-locking (see e.g. [20]). In other approaches adopting a closer semiconductor (s.c.) optical response [15,18], a standard adiabatic elimination of the macroscopic semiconductor polarization has been introduced, leading to the commonly termed rate equation model. While this approximation proved successful in describing the dynamics of (especially single mode) optical nonlinear systems, it is known that it entails unphysical dynamical effects in multimode lasers because it corresponds to an infinitely broad gain/dispersion [22]. In order to circumvent this limitation, spectral filtering terms were added in the description of broad-band QCLs [18]. On a more fundamental basis, other works report on nonstandard techniques for the adiabatic elimination of the polarization which formally amount to introducing diffusive terms in the equations for the multimode laser field [23]. In order to provide a fundamentally self-consistent modelization, it is desirable to take into account a frequency dependent, asymmetric dispersion/gain line and the Linewidth Enhancement Factor (LEF) that are known to play a critical role in coherent multimode regimes of s.c. lasers [24,25]. Thus, in this work we adapt to a QCL, a s.c. laser model whose basics have already been proved quite effective in the description of multimode instabilities in bipolar lasers [24][25][26]. Also, in order to reduce the complexity of the model, we retain the transition dynamics in the 2-level framework [1] 1. Specifically, we introduce an asymmetric gain/dispersion line in the linear-gain approximation fitting the gain peak position and width to the experimental data and account for the LEF without resorting to first principle derivation of the full polarization evolution. In the model, we consider a unidirectional ring cavity configuration; this hypothesis allows us to show that phase-amplitude coupling (associated to the LEF) and mode competition can provide, even in absence of SHB typical of Fabry-Perot (FP) configurations considered e.g. in [13,20], a radically lower threshold for coherent multimode dynamics than what predicted by the Risken-Nummedal-Graham-Haken (RNGH) instability [28][29][30], which is typical of two-levels laser. In this regard, we observe that unidirectional lasing in ring cavities in QCLs has been demonstrated in both a monolithically integrated cavity and an external one [31,32]. Recently, in the latter configuration, working in both unidirectional and bi-directional emission regimes, the first experimental evidence of active mode-locking in QCLs has been reported [10]. On a farther perspective we might envisage that compact QCL ring resonators would share a similar future with the Near-IR ring resonators that nowadays represent key elements in integrated photonics circuits [33]. The extension of the present model to a FP resonator, where other physical mechanisms such as SHB may play a role in the system dynamics and possibly hinder the formation of stable pulses via mode-locking [16,32] is left to a future work.
In Section 2 we (i) introduce the s.c. laser model, (ii) derive an analytical expression for the singlemode laser solutions in the free-running regime in the form of traveling waves (TWs) and (iii) study their respective thresholds when the LEF and free spectral range (FSR) are varied. As expecEFd, the TW closest to the gain peak has the lowest threshold, but at higher pump currents, we predict competition among several laser modes. In this case a "rule of thumb" can be obtained to gauge the laser operating regimes towards irregular (spatio-temporal chaos) or coherently locked emission, the latter leading to regular oscillations of the output field. In Section 3 our simulations show that increasing the bias current leads to the destabilization of the lowest threshold TW and drives the laser through a sequence of alternating irregular regimes and regular, mode-locked, regimes in agreement with recent experiments results on the multimode dynamics of an ultra-broad-gain THz-QCL [8,34]. Also, we find that a larger LEF (more typical of a Mid-IR QCL [1]) favors the TW destabilization towards multimode regimes. In the regular, mode-locked regimes we interpret the pulsed dynamics with the formation inside the cavity of self-confined structures, with a variable number of intensity peaks, which can also explain the experimentally observed disappearance of the fundamental beat note (BN) in the spectrum. Such structures have a clear solitonic character, not unlike those recently observed in long cavity (i.e. "quasi class-A") s.c.lasers [26]. They appear as low-contrast pulsed emission on a CW background and therefore represent an interesting example of a spontaneous mode-locking in a QCL laser. While such a small peak contrast is a limitation for applications, in Section 4 we show the generation of picosecond pulses, introducing a RF modulation in a portion of the cavity which might have a relevant impact on applications such as time resolved measurements, high-precision spectroscopy, multispectral sensing and imaging. In Section 5 we summarize our main results and trace the pathways to future investigations.

Effective Semiconductor Bloch Equations
We consider the ring resonator sketched in Fig. 1.a where a semiconductor active medium of length l ∼ 1 mm is placed in a unidirectional ring cavity of total length L ≥ l. We suppose that two mirrors have transmissivity T 0, while the others are perfect reflectors (T = 0). This configuration choice excludes the onset of the spatial hole burning multimode instability [16]. In the hypothesis of linearly polarized field, we writẽ where k 0 = ω 0 /v = ω 0 n/c, n = √ b is the background refractive index and ω 0 is the angular frequency of an empty cavity mode and it will be taken in the following as the reference frequency. Analogously, the medium macroscopic polarization can be written as: boundary condition reads explicitly where we have taken into account that ω0 = 2πc/λ0 and Λ is the optical length of the cavity which includes the presence of the background refractive index nB (but not the contribution to the refractive index from the two-level atomic system) In addition, the transmitted field is given by where L is the distance between mirror 4 and mirror 1 (see Fig. 8.4) and the reflected field corresponds to By passing to the normalized envelopes In the slowly varying envelope and rotating wave approximation, the radiation-matter interaction is then described by the following nonlinear partial differential equations: where N is the carrier density in the upper laser level (in the hypothesis of instantaneous relaxation of the lower state (τ 2 = 0) [1,14,16]), τ e = τ 32 is the carrier density nonradiative decay time, I and V are the pump current and the sample volume respectively, g = i N p ω 0 Γ c 2 0 nc , Γ c is the overlap factor between the optical mode and the active region and N p is the number of cascading gain stages. The Fourier components of P are proportional to the intracavity field via a complex susceptibility:P which can be referred to the resonant medium in a general way by assuming; so that, anti-transforming, we get the following dynamical equation for P: By considering complex quantities A = Re(A)+iIm(A) and B = Re(B)+iIm(B), the susceptibility takes the form: Note that it must be Re(B) > 0 in order for χ to represent an analytical function of ω in the upper half plane, and in particular we will take where τ d is the dipole dephasing time and Γ and δ are reasonably assumed independent from N [24]. As illustrated by Fig.2.a, the (asymmetric) gain and the refractive index curves associated to Re( χ) and Im( χ) now lack even or odd parity and moreover, the zero of the latter does not coincide with the maximum of the former. The the canonically symmetrical dispersion curve and Lorentzian gain of the two-level model can be readily recovered by assuming Re(A) = 0. i.e.: Figure.2.b evidences this difference. From Eq. (7) it follows that the maximum of Im( χ), and hence of the gain curve is given by: Therefore, in the limit Im(A)/Re(A) >> 1, (or in the two-level case, Re(A) = 0), one has ω M δ/τ d and the gain linewidth (FWHM) is approximately Γ/τ d , so that we recover the intuitive picture for the susceptibility behavior ( Fig.2

.b).
A further simplification can be introduced for moderate carrier densities (as it will be valid in the following sections), by modeling a linear dependence of the s.c. susceptibiity on N at the reference frequency (see e.g. Eq. (2.115) of [35] as follows: where α is the LEF. Hence using Eq.(11) we get: Setting δ = −αΓ (so that Re(A) = −2 f 0 NαΓ/τ d ) implies the reasonable assumption that the gain maximum coincides with the reference cavity mode (from Eq. (10): ω M = 0 ). In this way we recover the two-level description by setting α = 0. The real constants f 0 , α, Γ can be obtained by fitting the experimentally measured gain spectra reported for example in [8,36] around their maxima [24]. By using relations (8) and (11) in Eqs. (3)-(4) we get the Effective Semiconductor Maxwell-Bloch Equations (ESMBEs) [37] which correctly incorporate the peculiarities of the s.c. laser nonlinearities: We also remark that (1) the adiabatic elimination of the macroscopic polarization P obtained by setting ∂P ∂t = 0 in Eqs. (13)-(15) corresponds to assuming a flat gain curve and it leads to the rate equation model analogous to that adopted in [15] (in the three-level case) and (2) the adiabatic elimination of the macroscopic polarization P supplemented by the inclusion of a spectral filtering term in Eq. (13) to model the finite gain linewidth, leads to a modified rate equation model as adopted in [18] (in the three-level case).
A more compact form can be achieved by introducing the following variable changes and scaled parameters: The final formal steps are (i) the introduction of the boundary conditions modelling the ring cavity, which allows us to isolate effects of s.c. medium and mode spectra from SHB due to standing wave gratings, (ii) the application of the low transmission limit, which makes the concept of cavity mode well defined, and (iii) the assumption that the medium completely fills the cavity (L = l in Fig.1a) which is coherent with the focus on millimetric monolithic devices, although our model is apt to describe extended cavity configurations as those in [10,16]. Since the procedure is rather standard [30], we refer the reader to Appendix A where Eqs. (16)- (18) are worked in the form: whereẼ = E,P = AP,D = AD andμ = A µ with A =gl/T. Time is normalized to the fastest timescale τ d so that σ = τ d /τ p , b = τ d /τ e with τ p = nl/cT, and the longitudinal coordinate to the cavity length l.
Although our results might also apply to Mid-IR-QCLs, we adopt in the following typical values for a THz-QCL by setting τ d = 0.1ps, τ p = 10ps, τ e = 1ps (that give σ = 0.01 and b = 0.1). The considered value of τp is corrected to account for the total losses inside the laser cavity [19].

Free running traveling wave solutions
Due to the peculiarities of the s.c. laser we considered, the stationary solutions differ in intensity and frequency from the known ones pertaining to 2-level lasers. The solution is sought among TWs of longitudinal wavevector k and frequency Ω in the form: Setting to zero the RHS of Eqs. (19)-(21), using expressions (22) and dropping the tilde we have: and hence the stationary solutions for P s and D s are: where The stationary field solution from Eq. (23) is then: Two real equations follow, which constitute the generalization of the stationary lasing intensity and the (implicit) dispersion relation Ω(k) of the TW angular frequency for the considered QCL:

Traveling wave selection at threshold
A general stability analysis for the TW solutions derived above is left to a future work, while in this instance we will focus on the destabilization of the TW emitted at threshold and to the multimode regimes appearing beyond it. In particular we are interested in the role played by the mode separation and the LEF. The instability of the TW solutions in QCL models based on MBEs for 2-level systems in case of FP and ring configurations has been linked to the RNGH instability, associated with the parametric gain of the cavity modes in resonance with the Rabi frequency [13,20]. In our case, though, the phase-amplitude coupling points towards a different character, namely the Benjamin-Feir instability [38], that appears to be peculiar of the multimode s.c. laser dynamics as assessed in [24,39] and is triggered by the growth of modes having wavevectors larger than a critical value that depends on boundary conditions. In particular in [39] the authors, using to a codimension 2 bifurcation analysis and the hypothesis of instantaneous medium response ("class A" laser), show how a set of effective semiconductor Bloch equations analogous to those used here can be reduced close to threshold to the prototypical cubic Complex Ginzburg Landau Equation (CGLE) that describes the behavior of a large class of spatially extended nonlinear systems [40]. The "phase instability" (or Benjamin-Feir instability) of the TW solutions of the CGLE leads to a rich variety of spatio-temporal complexity encompassing phaseand defect-mediated turbulence (where respectively small or large amplitude modulations of the unstable TW occur) and coherent phenomena as the formation of localized structures. Because of the small carriers relaxation time that makes QCLs "quasi class A" lasers, we may thus expect to observe a similar dynamical scenario. We note that with respect to the RNGH instability the Benjamin-Feir instability has a much lower threshold, in agreement with the experimental evidences. Now, in order to determine the conditions where an easier destabilization of the single mode TW 0 emitted just above the lasing threshold, or a stronger modal competition emerging therefrom, can occur we proceed to derive the analytical pump threshold curve µ th of the TW lasing solutions (setting |E s | = 0 in Eq.(29)) and study it for different values of the LEF α and cavity length l (i.e. of the FSR). As it turns out, the curve shape and the mode spectrum can provide us with a first insight about the stability of the TW 0 . In Fig. 3.a and Fig. 3.b we report µ th as a function of k, considered as a continuous variable, for different values of α and l. We set Γ = 1.1 that corresponds to a gain linewidth of 1.75T Hz which is close to that of the broad-band QCL studied in [8,34]. Clearly, an increment of α or l lowers the lasing threshold for an increasing number of TWs. An easier destabilization of TW 0 and a more complex multimode dynamics can thus be expected. In Fig. 3.c we plot the pump threshold µ th against Ω for α = 1.5 [41,42] and l = 1mm. As indicated by the arrow, in this case the TW 0 angular frequency is Ω 0 = 0.015 that corresponds in physical units to a frequency separation of 24GHz from the gain peak. This will be the case study for the following sections.

Numerical results. Coherent multimode dynamics and chaotic regimes.
In this section we study the laser behavior emerging from the destabilization of the TW0 by numerical integration of the Eqs. (19)-(21) using a split-step method based on a second order Runge-Kutta and a FFTW algorithm. In the case of Fig. 3.c, we simulated the laser dynamics for increasing values of the pump parameter µ and we plotted the temporal evolution of the field intensity, the optical spectrum (OS) and the beat note spectrum (BNS) for specific values of the normalized pump parameter. In Fig. 4.a, just above threshold (µ = 1.002µ th with µ th = µ th (k = 0) = 1), we find as expected the emission of TW 0 , affected by a self-oscillation, whose signature is seen in the BNS peak at 0.062 (quite close to the normalized FSR = (2πc/nl)τ d = 0.0628 that corresponds to a frequency separation of 100GHz in physical units) and the presence of higher angular frequencies separated by FSR. Consistently, in the OS of Fig. 4b, we observe the main peak at the emission angular frequency Ω 0 = 0.015 of the solution TW 0 , according to Eq. 30, weaker peaks at Ω 0 + FSR and FSR − Ω 0 that are the signature of nonlinear parametric processes. The corresponding dynamical regime consists of regular small amplitude oscillations around an almost constant intensity.
Beyond the immediate proximity of the threshold (µ > 1.002µ th ), a multimode instability takes place and brings the system in a turbulent spatio-temporal regime. An example of the temporal plot of the field intensity taken at z = 1/2 is shown in Fig. 5 for µ = 1.2µ th . The BNS of Fig. 4.a shows the broadening of the first BN and the appearance of higher beatnotes separated by FSR. The corresponding OS of Fig. 4.b shows the same, but the broadening of the QCL modes is more pronounced and the wealth of modes setting off is more evident. A similar irregular behavior is found for µ = 1.4µ th . The next value of µ = 1.6µ th brings the system back to a regular regime whose field intensity dynamics is shown in Fig. 6, along with the OS. At steady state, both the plot of the intracavity field intensity at z = 1/2 versus time and the longitudinal intensity profile along the resonator at a fixed time reveal the presence of a self-confined, two-peaked structure that indefinitely travels in the ring resonator with constant speed and shape, as shown in Fig. 7. It represents a phenomenon of spontaneous mode-locking, although it involves a limited number of modes. The corresponding OS is in fact characterized by only 7 active modes within two decades (see Fig. 6.b). The presence of two traveling peaks is of course associated to an intensity spectrum showing a periodicity at half the round trip time, i.e. at twice the BN angular frequency. This feature can thus explain the disappearance of the fundamental BN observed in Fig. 4.a and in certain dynamical regimes reported in [34] (see for example Fig.4 of [34] in the region between 400 and 420 mA). We could a posteriori verify that this traveling structure is also present in the field profile just above threshold (µ = 1.002µ th ) where, of course, the low absolute emission  Color scale variations across these regions are due to graphical interpolation. The red rectangle highlights the region of the BNS where analogous experimental data are available [34]. Parameters as in Fig. 3.c. leaves a very low-contrast. We also observe that these results are stable upon inclusion of an additive, stochastic noise process in the simulations. At µ = 1.8µ th and at µ = 2.0µ th the dynamics is still regular in both cases (the number of locked modes in the OS in the two regimes is 9 and 10 within two decades, respectively). A proof in support of the solitonic character of the travelling structures described above is provided at µ = 2.0µ th . Figure. 8 shows how two different initial conditions can evolve at steady state into a single-peaked self-localized structure or to a two-peaked one propagating at a constant speed, a phenomenon of multistability, typical of dissipative solitons in extended systems [43]. Up to N = 3 stable structures could be obtained at regime from different initial conditions.
We may thus conclude that a good qualitative agreement exists with the results reported e.g. in Fig. 4 of [34], as our model properly reproduces (though for different values of the the bias-to-threshold current ratio) a sequence of narrow lines at multiples of the FSR, and broad bands.  Fig. 8. Intensity profile at a given instant of time obtained starting from two different initial conditions at µ = 2.0µ th : a single peak structure (down) and a two-peaked structure (up). Other parameters as in Fig. 4.
The absence of the multibeatnote regimes reported in [34] (see Fig. 5 in [34] at pump 440 mA) can be explained by the homogeneous gain line assumed by our model, as opposed to the inhomogeneous line of the THz-QCL in [34] which integrates three different active regions in the same waveguide. The different gain peaks probably introduce different BN with the common resonator. Moreover, our numerical simulations show that the extension of the regions of regular and irregular system dynamics, as well as the number of competing longitudinal modes can be controlled by acting on the QCL intrinsic parameters such the LEF α, the carrier decay time τ e , the cavity length l and the FWHM of the gain curve (∼ Γ). As an example, in Fig. 9 we report the BNS for a value of Γ almost three times smaller than that used to produce Fig. 4 and hence comparable with the gain linewidth of the single stack active region considered in [8,34].
Although, as expected, the number of competing modes becomes smaller, we observe a dynamical scenario characterized by a stable TW emission up to µ = 1.04µ th , and alternating windows of regular oscillations and chaotic dynamics for 1.04 < µ/µ th < 1.14.

Effect of RF modulation
An interesting result reported in [34] was the line narrowing in the BNS caused by the application of a RF modulation to the current pump. The narrowing was interpreted as the onset of comb coherence cased by the external "forcing". Our model proved capable to reproduce this effect, though only qualitatively, when an external RF modulation is added to the pump current to just a small portion of the cavity in the form µ(t) = µ M (t) = µ 0 + µ A cos(Ω M t), with µ A = 0.1µ 0 and where Ω M is close to the first beatnote observed in absence of modulation. In the remaining part of the cavity a DC bias a µ = µ 0 was maintained. As shown in Fig. 10, the first evidence is that the addition of the RF modulation causes the fundamental BN to reappear for pump values where it was absent (in agreement with experimental evidences, as in Fig. 10 of [34]). The second effect, better illustrated in Fig.  11, is the important reduction of the BN linewidth with a transition to coherent dynamics. The phase-noise reduction is in fact a clear manifestation on an increased coherence in the multimode emission. Starting from a chaotic state at µ 0 = 1.4µ th (whose spectrum is the red line in Fig.11) we observed a transition to a regularly oscillating (multimode) state when the RF is turned on (black line in Fig.11). The qualitative transition is reproduced, tough the line shrinking is much less pronounced than that reported in experiments [34]. In this respect we note that, at difference from the experiment, a RF applied to the whole cavity does not yield the same effect and it leaves the system in the chaotic state with little effect on the BNS. A possible explanation of this apparent discrepancy may lie in microwaves attenuation, estimated of few dB/mm [44], that removes spatial uniformity in the experiments. No RF modulation RF modulation Beat note spectrum ω Fig. 11. Beat note spectrum of the electric field emitted by the QCL in presence (black line) and absence (red line) of RF pump modulation. In the simulations we used µ 0 = 1.4µ th and a RF modulation, applied only to a just 2/5 of the whole laser cavity, with an amplitude equal to 35% of the mean value and Ω M = 0.062. Other parameters as in Fig. 4.

Active mode-locking
Finally, we checked that our model could robustly reproduce active mode-locking with the generation of picosecond pulses, as experimentally verified [5,6,9,12]. To this purpose we applied a strong current modulation at the round trip frequency ω M to a short cavity section. In particular we set in 1/7 of the cavity length, while in the remaining part a standard DC bias (µ = µ b < µ th ) is maintained. In Figure 12 we plot the field intensity at z = 1/2( Fig. 12.a) and the intracavity field intensity ( Fig. 12.b) versus time at steady state for µ b = 0.6µ th , µ 0 = 1.5µ th and µ A = 4µ th . Fig. 12 shows the formation of pulses with FWHM of 1.5 ps, a repetition rate Ω M and a contrast S = (Max(I) − Min(I))/Max(I) = 1. The corresponding OS and BNS are shown in Fig. 13 where the first is characterized by 13 modes in the first two decades. Robustness was assessed by reaching the same regime with varied parameters: e.g. for a DC bias slightly above the lasing threshold (µ b = 1.1µ th ) and we also verified that a variation of few percents in the pulse width can be obtained by varying the modulation amplitude µ A . We checked that active mode-locking is also achieved in case of faster carriers (down to subpicosecond time scale), and broader gain linewidth (up to 10T Hz range). As expected, slower carriers are associated to longer pulses and above the picosecond time scale, the effect disappears because of the medium inertia. An experimental validation of our numerical predictions would pave the way towards a number of fascinating applications based of actively mode-locked THz-QCL.

Conclusion
In conclusion we have established a unidirectional broad-area QCL model describing multimode regimes which focuses on features typical of active s.c. media, such as an asymmetric dispersion/gain line and a LEF. We have analytically derived an expression for the TW solutions and their lasing thresholds. While retaining only a minimum of mechanisms typical of the complex QCL device, our model proves capable of predicting low threshold multimode regimes, coherent multimode dynamics and it is capable to reproduce a number of experimentally confirmed features typical of broad-band THz-QCLs, such as the alternance between broad and narrow BN, but also of coherent regimes where we found that a narrow BN corresponds to stable self-confined structures that travel in the cavity at the group velocity of light and whose number depends on initial conditions. From a fundamental point of view, in spite of their low-contrast due to the small number of active modes, these coherent structures represent an interesting and unique phenomenon of spontaneous mode-locking in a QCL. Also, we could validate how the effect of a radio frequency modulation of the pump facilitates a coherent dynamics and leads to active mode-locking where the laser emits picoseconds pulses. These evidences might have a serious impact on applications such as high-precision spectroscopy and time resolved measurements in the Mid-IR and THz windows of the electromagnetic spectrum.