Analytical expression for Risken-Nummedal-Graham-Haken instability threshold in quantum cascade lasers

We have obtained a closed-form expression for the threshold of Risken-Nummedal-Graham-Haken (RNGH) multimode instability in a Fabry-P\'erot (FP) cavity quantum cascade laser (QCL). This simple analytical expression is a versatile tool that can easily be applied in practical situations which require analysis of QCL dynamic behavior and estimation of its second threshold. Our model for a FP cavity laser accounts for the carrier coherence grating and carrier population grating as well as their relaxation due to carrier diffusion. In the model, the RNGH instability threshold is analyzed using a second-order bi-orthogonal perturbation theory and we confirm our analytical solution by a comparison with the numerical simulations. In particular, the model predicts a low second threshold in QCLs. This agrees very well with experimental data available in the literature.


Introduction
Nowadays quantum cascade lasers (QCLs) emitting in the mid-infrared (MIR) spectral range find widespread applications in health care, gas sensing, telecommunications and broadband spectroscopy [1]. While many applications require QCLs operating strictly in a single mode regime, QCLs producing frequency combs are proven to be very efficient in spectroscopic applications [2]. However, the majority of QCL based frequency combs do not exhibit phase locking relationship between individual modes. So far, a very moderate progress was made towards achieving ultra-short pulse production regimes in QCLs operating under room temperature conditions [3]. At the same time, ultra-short optical pulse generation by QCLs would enable applications in other important domains such as time resolved pump-probe MIR spectroscopy of rotovibrational modes of molecules or the environmental remote sensing of molecules in the mid-infrared spectral range.
One very interesting proposal on how to produce short MIR pulses in a QCL [4] has emerged when a multimode Risken-Nummedal-Graham-Haken (RNGH) like instability [5,6] (or so-called "second threshold" of a laser) was experimentally observed in these lasers. Although the interpretation of measured autocorrelation traces in particular QCLs was not straightforward, one can expect that the lasing dynamics in a QCL might be tailored so as to produce regular RNGH self-pulsations.
Interestingly, low-threshold RNGH instability was observed both in ridge waveguide QCLs and in buried heterostructure QCLs [4,7,8]. These experimental observations strongly disagree with the second threshold predicted by the RNGH theory [5,6]. Although the role of spatial hole burning (SHB) in lowering the second threshold was recognized, a conclusion was drawn in [7] that RNGH instability does not occur just as a result of the induced grating of carrier population. In [7], an additional assumption was made about a built-in saturable absorber in the cavity of QCL, which reduces its second threshold from the 9-fold excess as predicted in [5,6] to about 1.1 times above the lasing threshold as it was experimentally measured. As a matter of fact, the saturable absorber has been shown to lower the instability threshold in a laser [9]. However, the nature of saturable absorption in QCLs has never been fully clarified [1]. In particular, in [7], the Kerr lensing effect has been proposed as a possible mechanism responsible for the saturable absorption. In case of narrow ridge waveguide lasers, due to the overlap of the waveguide mode tails and the metal contact deposited on the waveguide, the Kerr lensing may lead to a saturable absorption. However it cannot produce the saturable absorption effect of the same strength in buried heterostructure QCLs. At the same time low-threshold RNGH instability was observed in both ridge waveguide and buried heterostructure QCLs.
In our recent study [10] we have analyzed conditions for multimode RNGH instability in a Fabry-Pérot cavity laser and found that a set of inaccuracies was admitted in the previous treatment [7]. As a result we propose an alternative mechanism responsible for low second threshold in QCLs. It also favors the second threshold lowering by the SHB effect but does not require a saturable absorber. More specifically we show that a combined effect of the carrier coherence grating, carrier population grating and relaxation processes due to carrier diffusion leads to the onset of multimode instabilities at Rabi flopping frequency. We consider such spectral behavior as an evidence of RNGH instability. Our approach to calculate the second threshold has shown a convincing agreement with the available literature data and with the results of numerical simulations. From a practical point of view, this can be converted into a useful engineering tool enabling to examine particular QCL design and predict its dynamic behavior or to tailor QCL design for a specific application. In this paper we obtain a simple closed-form analytical expression for the second threshold in a QCL with monolithic FP cavity. The main result of this paper is summarized by Eq. (14). Technically, it is obtained by using the second-order bi-orthogonal perturbation theory applied to the Lyapunov stability analysis of the single-mode lasing regime in a QCL.

Theory of the 2 nd threshold
The RNGH instability in a continuous wave (CW) single-mode laser arises due to Rabi splitting of the lasing transition induced by the lasing mode [11]. As a result of such spectral broadening and reshaping of the gain curve, the laser can provide sufficient optical gain for other longitudinal modes. If below the RNGH instability threshold, the optical mode is fully controlled by the laser cavity, at above the second threshold, the dynamic behavior is very different. The multimode RNGH instability is related to excitation of rapid coupled oscillations of the medium polarization P and population inversion N [10]. The buildup of Rabi oscillations at above the second threshold indicates that now a macroscopic polarization (coherence) has the major impact on the optical field in the laser cavity [6]. The medium polarization thus exhibits the opposite behavior to that which is known in the literature as the adiabatic following approximation [12].
RNGH instability is thus a multi-longitudinal-mode amplitude instability of a laser [13] which occurs on a very fast time scale. It cardinally differs from other possible instabilities in semiconductor lasers such as transversal mode instability in broad area lasers [14] or relatively slow phase instabilities in semiconductor laser [15,16], which admit adiabatic approximation for the medium polarization.
In order to find the gain curve for instabilities in initially non-lasing cavity modes, one applies the linear stability analysis to the Maxwell-Bloch equations (MBE) describing a single-mode CW laser [5,6]. In Appendix A we briefly present the original model equations, while further elaboration is given in [10]. Here we provide a brief account of these. Starting from MBE and accounting for the standing wave pattern of the optical field in the cavity, we first get a truncated system of coupled-mode equations for the field, carriers and medium polarization as well as their spatial harmonics. We then find a stationary solution for the single-mode CW lasing regime and perform its linear stability analysis with respect to small perturbations to the lasing mode. Then using the ansatz from [5], we recast the perturbations in the form of propagating wavelets exp( / ) g in z c t    , where Λ can be interpreted as Lyapunov exponent, Ωn g /c is the detuning of the propagation constant from the main lasing mode and n g is the group velocity index (for simplicity, we call the parameter Ω as a frequency offset [10]). The resulting linearized system of differential equations is defined by a 9×9 block-diagonal matrix M 9×9 [10]. The dispersion of the instability increment Re(Λ) as a function of the detuning Ω then follows from the eigenvalue problem of M 9×9 : 9 9 det( ( ) ) 0 M I      ( 1 ) where I is the identity matrix. It appears that up to 40-fold excess of the pump rate above the lasing threshold of QCL, only one eigenvalue in Eq. (1) may have a positive real part Re(Λ) [10].
The purpose of this paper is thus to find an analytic expression for the second threshold condition when Re(Λ)>0 and leads to a multimode RNGH instability. Note that the matrix M 9×9 is composed of two main diagonal blocks of 4×4 and 5×5 sizes and the eigenvalue with a positive real part is related to the smaller size (4×4) block, which we show below in Eq. (2). (The remaining 5×5 block exhibits only stable solutions in the range of pump rates that are used in practice [10].) For the reason we explain later, we split this 4×4 matrix into two matrices M (0) and M (1) : where the coefficient √2 appears due to the standing wave pattern of the optical field in the cavity. Here p is the pump rate normalized to the lasing threshold in the absence of SHB, while parameter ν 0 (p) accounts for the SHB effect, see also [10]. It increases the (effective) lasing threshold and reduces the slope efficiency:  T  TT  T  TT  T  TT  p  p  p  p  T  TT  T  TT  T Note that ν 0 is independent of the photon lifetime in the cavity. Figure 1(a) displays the effective threshold behavior in QCLs emitting at 10 and 8 µm (QCL 1 and QCL 2) as compared to a typical quantum well laser diode (QW LD) operating at 820 nm wavelength (green curve). For all considered lasers we assume a simple single-section Fabry-Perot (FP) cavity design and cavity length of 4 mm. In Table 1 we summarize the model parameters [7,17] for the three types of lasers considered here.
In the case of LD, the influence of the SHB effect on the effective lasing threshold Eq. (3) is negligible in the considered range of pump rates (p<3). This can be attributed to the gain recovery time T 1 in LDs being three orders of magnitude higher than T 2 . Our model curves for the effective threshold in QCLs plotted in Fig. 1(a) take into account large dispersion in the relaxation time constants due to possible design variations in the active region of QCLs. In both considered QCL examples, the effect of SHB on the effective lasing threshold is significant. As expected, at the lasing threshold (p=1), when the optical field is weak, the SHB does not affects the output characteristics of QCL and ν 0 =1. However at higher pump rates (p>1), the SHB effect sets in, yielding an increase of the effective threshold ν 0 . At very high pump rates (p>>1, not shown in the Fig. 1(a)), ν 0 asymptotically approaches the value of . Note that despite the fact that longitudinal and transversal relaxation times in the two QCL examples differ by a factor of two, the difference in ν 0 behavior is marginal. Figure 1(b) illustrates the effect of SHB on the Rabi oscillation frequencies (Eq. (4)) linked to the amplitudes of the optical fields in the cavities of considered LD and two QCL examples. The Rabi oscillation frequencies are plotted with (solid curves) and without (dashed curves) accounting for the SHB effect. In both QCLs, the SHB reduces the Rabi flopping frequency due to the effect on the lasing threshold ν 0 . Thus at the normalized pump rate p=3, the reduction in Rabi frequencies of QCLs is about 10% (compare the solid and dashed curves in Fig. 1(b)). However as a consequence of ν 0 behavior in LD depicted in Fig  1(a), the reduction of the Rabi frequency in a LD is small (red solid and dashed curves overlap in the figure).
In Figure the instability this example, is 100 µm. T separation of 2(a), the inst round-trip gai non-lasing cav above the las the dispersion modes. There function of th equency Ω is n rresponds to t s and the maxi maximum of t abi (further deta behavior as a s established fo e for some of pump rate is at tion grating an g. 2(a), the ga p=1.1. Howev ccur at such pu the spectral be r stability matr me that the cavi use of the large ed example. In p time, yieldin urves for instab om p=1.1 to p= 2) does not ac condition of th es by plotting normalized to the lue curves) and for coherence grating cies are calculated ngth is 4 mm in all normalized to t the initial lasin imum of the g the gain for m ails age given an evidence o for a CW unidi f the non-lasin least 9 times a nd the carrier c ain for instabili ver from the n ump rate.
havior of rix (2). In ity length e spectral n the Fig.  ng Fig.  2(a). We use the travelling wave (TW) rate equation model that is similar to the one in Ref. [18] but does not assume a saturable absorber in the cavity and accounts for the diffusion of carriers and coherences [10]. Time domain simulations are performed for the QCL being initially in non -lasing state. The TW model incorporates Langevin force terms that seed the spontaneous polarization noise into the system. We vary polarization noise power over six orders of magnitude, from eight to two orders below the polarization noise level in homogeneously broadened ensemble of two-level quantum oscillators (see Eq. (5.5.22) in [19]). In all cases, we have not seen the impact of the noise power on the lasing regime. Only the delay time to the onset of emission was changing. For each pump rate p, we perform a series of numerical simulations [10] and count the number of occurrences of RNGH selfpulsations as opposed to the CW lasing regime (for details see Figs. 4. and 5. in [10]). The criterion to determine whether the particular simulation evolves into RNGH self-pulsations is the observation of oscillations in the output power (Fig. 2(c)) and RF power spectrum (Fig. 2  (d)) as well as the peculiar shape of the system attractor in the P-N plane (polarization-carrier density variables) reaching high coherence values P on the system trajectory, that is with the maximum of P approaching a half maximum of N. The P and N variables used in Fig. 2 (e) are introduced following the approach of [18]. In particular, the variable P measures the order parameter in the system (the carrier density in the coherent state) and has the same units as the carrier density N. In Fig. 2(e), both variables P and N are normalized on the transparency carrier density. More details can be found in Ref. [10].
The self-pulsations in a short-cavity QCL are observed at frequencies very close to the cavity round-trip frequency (or its harmonic) and when it (or its harmonic) approaches Ω Rabi [10]. In Fig 2(b) we show the average probabilities to obtain multimode RNGH selfpulsations. Each point is the result of 80 realizations so as the uncertainty is of ±0.056 (shown as error bars). Note that out of 80 realizations at p=1.5, no occurrence is detected with RNGH self-pulsations. The RNGH instability starts to develop at p between 2 and 2.5, with p=2.5 being the first point for which such probability exceeds the uncertainty range. We clearly see that this is a higher pump rate than one would expect from the approach of [5,6] based on the positive instability increment for single mode CW lasing [the gain curve at p=1.1 in Fig.  2(a)].
The positive increment for instability of single mode CW lasing regime in Eq. (1) indicates the instability of the lasing mode but it does not imply that RNGH self-pulsations will occur. The laser may just switch to another mode, which is also unstable. The increment for instability of the mode is small [e.g. see the cyan curve in Fig.2 (a) obtained at p=1.3] so the switching process will be slow and may be driven by technical noise. A classic example is a ring laser with reciprocal cavity for counter-propagating modes and absence of mode coupling via backscattering in the cavity: it shows sporadic switching of the lasing direction (spontaneous mode-jumps) at times much longer than the characteristic times T 1 , T 2 and τ [20]. A FP cavity laser may also exhibit such spontaneous mode-jumps [16] and the linear stability analysis in Eq. (1) does not exclude such possibility. Other possibility is related to multimode phase instability in semiconductor laser which admit adiabatic approximation for the medium polarization [15,16]. (Note that since RNGH instability is the multimode amplitude instability [13], our numeric model utilizes slowly varying envelope approximation (SVEA) for the field amplitude and therefore it cannot reproduce spontaneous mode-jumps or slow phase instabilities).
What distinguishes RNGH instability from all other instabilities is that, when the regime of RNGH self-pulsations sets in, the medium polarization does not follow adiabatically the optical field in the cavity. However, at the initial stage when the CW single-mode lasing regime becomes unstable, there is only a small-amplitude perturbation of the optical field circulating in the cavity. In the most general case it has the form of an optical pulse. The interaction of the active medium with such a small-amplitude perturbation is still governed by adiabatic approximation [12] and does not necessarily give rise to a non-adiabatic behavior at a later stage. In order to define when the non-adiabatic behavior develops, we apply the pulse area theorem [19,21]. According to it, a perturbation in the form of an optical pulse is unstable and the pulse area grows if Ω Rabi ·τ p >π, where τ p is the characteristic pulse width. Otherwise, the medium polarization just tracks the optical field behavior. The most general form of a small-amplitude perturbation that initially circulates in the cavity is such that its characteristic time is roughly a half of the cavity round-trip time Ln g /c. Using also the fact that the maximum gain for instability is shifted from the initial (main) lasing mode by the Rabi frequency [ Fig. 2(a)], we recast the pulse area theorem into the form of instability condition at the maximum gain frequency max / 2 / 2 g c Ln Fig.   2(a), the first case where the gain maximum is on the first non-lasing cavity mode occurs at p~2.5, when the numerical simulations do reveal the RNGH self-pulsations [see Fig. 2(b)]. The dashed vertical line in Fig. 2(b) shows the exact location of the second threshold obtained from the unstable pulse area condition. Because of the large frequency separation between the cavity modes, the second threshold raises considerably in short-cavity devices. At the same time, the threshold condition originating from the pulse area theorem has a minor impact on QCLs with long cavities, where mode separation frequency c/2Ln g is small. Thus for QCLs with cavity lengths of 3-4 mm, our linear stability analysis predicts that the second threshold is at a few percent above the lasing threshold [10], in agreement with available experimental data [4,7,8]. In addition, our numerical simulations based on TW rate equations model confirm the RNGH multimode operation at small excess above the lasing threshold [10]. At first glance, it is quite straightforward to solve the eigenproblem of the matrix from Eq. (2) and to find an analytic expression for the pump rate p th2 at which one fulfills the second threshold condition: for the spectrum of the positive Lyapunov exponent. However the characteristic equation for the Lyapunov exponent is of the fourth order. Therefore we use a perturbation-theory approach. Since the matrix in Eq. (2) is not Hermitian, its eigenvectors are not orthogonal. For that reason we use the second-order bi-orthogonal perturbation theory [22,23] and split the linear stability matrix into the zero-order approximation M (0) and perturbation M (1) matrices, as indicated in Eq. (2). An example of the eigenvalue spectra of the initial matrix M (solid black curves) and its zero-order approximation M (0) (red dashed curves) for a QCL pumped at p=1.5 times above the lasing threshold is depicted in Fig.3. The eigenvalue with the largest real part Λ max is associated with the upper 2×2 block in M (0) , which is independent of Rabi  . In the zero order approximation, the increment ) Re( ) 0 ( max Λ just barely increases with the pump rate p due to the SHB effect accounted for in the parameter ν 0 : where τ=n g /cl 0 is the photon lifetime in the cavity. Note that 0 at the lasing threshold (p=1) and zero detuning from the main lasing mode ( 0 Ω  ). The second eigenvalue originating from this 2×2 block has the largest negative real part, of ~0 Fig. 3). In contrast to this behavior [See Eq.(7)], the two eigenvalues originating from the lower 2×2 block in M (0) are independent of the spectral detuning Ω but exhibit splitting that varies with Rabi  . At the lasing threshold (p=1) the splitting is of   Table 1. It can be seen from the example in Fig. 3 plotted for p= 1.5, that the eigenvalues of the matrices M and M (0) are close to each other as compare to the splitting. The perturbation matrix M (1) has no diagonal elements. Therefore the first-order correction terms to the eigenvalues of M (0) vanish. Hence we apply the second order of bi-orthogonal perturbative expansion in Eq. (2)   which provides a realistic approximation for the spectrum of RNGH instability gain ) Re( max Λ in the vicinity of its maxima at Ω~± Rabi  (see Fig. 4, which we will discuss below). In Eq.
(8), the diagonal term due to the first-order correction vanishes while the second-order correction terms are of the order of g  Rabi  eff  Rabi   T  T   ,  2  2  ,  2  2 ,   , which is much smaller than the initial spectral splitting indicated in Fig. 3. This justifies our decomposition onto the zeroorder approximation M (0) and perturbation M (1) matrices in Eq.(2). In (8) where coefficients C i (p) and A(p) are independent of the offset frequency Ω (they are shown in the Appendix B). The first two terms in Eq. (9) originate from 0 ( 1 0 ) Using the expressions from Appendix B for parameters A and 0 1 / C C  , we finally obtain that the highest gain for multimode RNGH instability is at the following offset frequencies: This expression is valid when , that is true in all cases of practical interest considered here. The frequency max  in Eq. (11) monotonically increases with the pump rate (e.g. see Fig. 4). Using Eq. (4), we find a reciprocal relationship to Eq. (11): ( 1 2 ) In Appendix C we obtain a general solution of Eq. (12) with respect to p. It appears that the instability increment (9) is negative and multimode instability is impossible at very small pump rates, when is the pump excess at which 0 2 max   in Eq. (12). For the pump rates above this value, the increment (9) is positive at the frequency (10) [see Fig. 2], which may lead to multimode instability. Finally to obtain the second threshold p th2 , we substitute in Eq. (12) our condition (6) for the max  and we get The second threshold reduces towards min p [Eq. (13)] with increasing cavity length L. The meaning of min p corresponds to the notion of the second threshold in the original RNGH theory. Nevertheless by no means Eq. (13) can be directly compared with the second threshold for a unidirectional ring laser obtained in [5,6] (see also discussion in Section 3.3). However the shape of the multimode instability gain curve (9) can be assessed against the one in a unidirectional ring laser. For this purpose we provide an approximate expression for the maximum of increment (9) at the offset frequency (11), which we will use later:

Validation of the analytic solution for the second threshold
In Fig. 4 we compare the instability gain curves obtained from the analytic expression (9) (red curves) and by solving numerically the eigenvalue problem of the linear stability matrix (2) (black curves). The results are plotted as the round-trip gain 2Re(Λ)n g L/c vs the offset frequency Ω/2π. In this example we use a QCL with the cavity length of 4 mm and other parameters from the first column in Table 1. For the pump rate p from 1 to 2.5 times above the lasing threshold Fig. 4. The spectral behavior of the round-trip gain coefficient for multimode instabilities in QCL with the cavity length L = 4 mm (The cavity mode separation c/2Ln g =11 GHz) at different pump rates p in the range from 1 to 2.5 times above the lasing threshold. The red curves show the gain spectra obtained from the analytical expression (9) and the black curves are obtained by solving numerically the eigenproblem of the matrix (2). The parameters used in simulations are listed in the column "QCL 1" in Table 1. and in the large spectral range, our analytical expression (9) is in reasonable agreement with the RNGH round-tip gain obtained by solving numerically for the eigenvalues of the matrix M in Eq. (2). It is interesting to compare the instability gain curves for long -cavity QCL in Fig. 4 (the cavity mode separation is 11 GHz) with the ones obtained in Fig. 2(a) for a QCL with the cavity length of 100 µm (the intermodal frequency is 450 GHz). For concreteness we take the gain curves at p=2.5. Both gain curves are peaked at the same offset frequency of Ω max /2π~450 GHz. Our Eq. (11) also shows that max  is independent of the photon lifetime in the cavity. Figure 5 provides a comparison between the frequency Ω max /2π obtained from the analytical expression (11) (red curve) and from the numerical solution of the eigenproblem of the matrix (2) (black curve). Our analytical expression for Ω max is in reasonable agreement with the one obtained numerically. The Ω max frequency approaches the Rabi oscillation frequency Ω Rabi at high p (blue curve). In order to illustrate the second threshold condition (6) we show a horizontal line at Ω/2π=c/2Ln g for the cavity length of L=100 μm. Its crossing point with Ω max /2π curve defines the second threshold, which is indicated by a vertical line. Our analytic approach provides a realistic approximation for the second threshold (the crossing point of the horizontal line with the red curve). exact, p=2.5 perturb., p=2.5 exact, p=2 perturb., p=2 exact, p=1.5 perturb., p=1.5 exact, p=1.1 perturb., p=1.1 Fig. 5. Peak gain frequency of RNGH instability calculated numerically from our linear stability matrix (black curve) and from our analytic expression Eq. (10) (red curve) is plotted vs. pump normalized to lasing threshold p. It is shown in comparison with the behavior of Rabi oscillations Ω Rabi /2π (blue curve). According to Eq. (6), the intersection between horizontal dotted line at c/2Ln g and the peak gain frequency defines the value of second threshold (vertical dotted line) in case of the 100 µm long cavity. Fig. 6. Second threshold (represented as a relative pump excess above the lasing threshold p th2 -1) vs. cavity length. We compare p th2 calculated from our Eq. (14) (curves) and calculated by numerical solving of the eigenproblem of the matrix (2) (squares and triangles). The parameters for QCL 1 and QCL 2 are listed in Table 1. Figure 6 shows variation of the second threshold with the length of QCL cavity. (We plot the relative excess above the lasing threshold p th2 -1 in logarithmic scale.) We consider two sets of QCL parameters from Table 1. In both cases, our analytical expression (14) (curves) shows good agreement with the numerical solutions for the eigenvalues of the linear stability matrix (2) (points). Note that the second threshold behavior predicted in Fig. 6 is in agreement with the results of numerical simulations based on the TW model [10] (see also Fig. 2) and with the numerous experimental reports indicating low second threshold in QCLs with the cavity length of 2-4 mm [4,7,8].
Our small signal analysis is based on a truncated system of coupled-mode equations [10] and our analytical expressions are valid in the limited range of the pump rates, when  Table 1) this condition implies that p<<10 4 . However since we account only for the first nine coupledmode equations in (1) and do not analyze the truncation error at high pump rates, our analysis . Instead, it shows a linear growth which is depicted in Fig. 1(a) by the dashed blue curve. As a result, both the Rabi frequency and the maximum increment frequency max  (11) decrease. We obtain the following approximate expressions for these two frequencies when 0