## 1.

## Introduction

Stability of pulse generation from midinfrared (MIR) quantum-cascade lasers (QCLs)^{1}^{,}^{2} has been an active field of research owing to the fact that stable MIR pulses are desired for free-space communication, medical diagnosis, and environmental sensing applications,^{3} and QCLs are important sources for such pulses. However, a full picture of stability analysis of QCLs is inevitably complicated when compared to conventional lasers because of the unique combination of ultrafast carrier scatterings and gain recovery,^{4}5.^{–}^{6} significant nonlinearities,^{7} and dispersion effect^{8} in a QCL medium.

There have been several discussions on the stability of typical MIR QCLs both theoretically and experimentally;^{9}10.11.^{–}^{12} however, none of these discussions specifically addresses the effect from group-velocity dispersion (GVD). In Ref. 9, observations of the coherent Risken–Nummedal–Graham–Haken (RNGH)-like^{13}^{,}^{14} multimode instability in QCLs have been reported. It is found that the instability threshold is significantly lower than that of the original RNGH instability, which is attributed to the saturable absorber effect in the lasing medium. Thereafter, more comprehensive theoretical and experimental investigations on the multimode operation regimes in QCLs have been documented in Ref. 10. It is shown that the fast gain recovery of QCLs exhibits two kinds of instabilities: one is associated with the spatial hole burning (SHB) effect and the other is RNGH-like instability. The QCLs’ stability analysis is based on the linearization of the Maxwell–Bloch formalism. The RNGH instability of QCLs exhibits two differences from that of the conventional RNGH instability, namely, lower instability threshold and missing of the central spike in the optical spectrum. According to the analysis reported in Ref. 10, the first difference is owing to the presence of a saturable absorber effect in the QCL cavity and the second difference is owing to the SHB effect.

The saturable absorber effect originates from the Kerr nonlinearity^{15} in QCLs, and it is also called the transverse Kerr effect. Accompanying the saturable absorber effect, there is the longitudinal Kerr effect, which induces self-phase-modulation (SPM) of lasers. A discussion of the effect of SPM on the QCLs’ instability is reported in Ref. 11, where it is found that the SPM could significantly broaden the unstable domain of the RNGH instability but it only has a trivial effect on the instability threshold. Moreover, the instability of the QCLs is viewed at a different angle,^{12} where the instability mechanisms are decoupled into the amplitude and phase domains based on the symmetry and antisymmetry of the propagating fields in the lasing cavity. The amplitude instability is of the RNGH-type multimode instability and the phase instability is of a single-mode nature. Both saturable absorber and SPM effects are accounted for in the analysis. The saturable absorber effect could lower the threshold of both amplitude and phase instabilities, while the SPM effect could broaden the frequency domain of both instabilities.

Similar to the Kerr nonlinearity, the GVD also originates from the intersubband transitions of the QCLs and it is embedded in the same structure as that of the Kerr nonlinearity. However, so far, its interaction with either saturable absorber or SPM has never been accounted for in the stability analysis of a typical QCL structure. A measurement of GVD in the MIR QCLs has been reported in Ref. 8. An ultrafast upconversion technique based on the sum-frequency generation was applied to investigate the total GVD by measuring the wavelength-dependent pulse propagation delay.^{8} The measurement on a $\sim 5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{m}$ QCL with an InGaAs/InAlAs active region shows an anomalous GVD with the coefficient of ${\beta}_{2}\sim -4.6\times {10}^{-6}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{ps}}^{2}/\mu \text{m}$.^{8} It is stated in Ref. 16 that the low GVD of the QCLs tends to lock the longitudinal modes of the QCL designed for the MIR frequency comb. The interaction between the GVD and the saturable absorber in the self-induced transparency (SIT) modelocking of QCLs is discussed in Ref. 17, where analysis based on the pulse shape evolution demonstrates that the stability limits for saturable gain and saturable loss with the presence of GVD depends on the ratio of the SIT-induced gain and absorption to the linear loss.

In our current work, we focus our discussion on the effect of the GVD on the stability of a typical MIR QCL without the mode-locking mechanism. The novelty of this paper is in the inclusion of the GVD effect in the analysis of instability mechanisms of MIR QCLs. We particularly investigate the interaction between the GVD and the saturable absorber on the QCL stability. At this point, we have not involved the SPM effect, as we are planning to take it in our next work. These works are toward our goal of investigating the possibility of forming an optical soliton in the QCLs for MIR pulse generation.

The rest of the paper is organized as follows. In Sec. 2, the linearization of nonlinear Maxwell–Bloch equations with the GVD effect and the decoupling of amplitude and phase stability mechanisms is described. Section 3 presents the solution to quadratic eigenvalue problems that result from the GVD effect in the amplitude and phase stability analysis. In Sec. 4, results are presented and discussed. Conclusions are drawn in Sec. 5.

## 2.

## Linear Stability Analysis with Second-Order GVD

The dynamic behavior of QCLs is modeled based on the two-level Maxwell–Bloch formalism. The QCL gain medium is assumed to be built with the Fabry–Perot cavity, thus, the SHB effect, brought by the standing waves, are accounted in the model. Both the second-order GVD and saturable absorber effects are included as well. We have the Maxwell–Bloch equations as

## (1)

$$\frac{\partial {E}_{+}}{\partial t}=-\frac{c}{n}\frac{\partial {E}_{+}}{\partial z}-\frac{c}{n}{P}_{+}-\frac{c}{2n}({l}_{o}-\tilde{\gamma}{|E|}^{2}){E}_{+}-i\frac{c}{2n}{\beta}_{2}\frac{{\partial}^{2}{E}_{+}}{\partial {t}^{2}},$$## (2)

$$\frac{\partial {E}_{-}}{\partial t}=+\frac{c}{n}\frac{\partial {E}_{-}}{\partial z}-\frac{c}{n}{P}_{-}-\frac{c}{2n}({l}_{o}-\tilde{\gamma}{|E|}^{2}){E}_{-}-i\frac{c}{2n}{\beta}_{2}\frac{{\partial}^{2}{E}_{-}}{\partial {t}^{2}},$$## (3)

$$\frac{\partial {P}_{+}}{\partial t}=-\frac{1}{2}({D}_{o}{E}_{+}+{D}_{2}{E}_{-})-\frac{{P}_{+}}{{T}_{2}},$$## (4)

$$\frac{\partial {P}_{-}}{\partial t}=-\frac{1}{2}({D}_{o}{E}_{-}+{D}_{2}^{*}{E}_{+})-\frac{{P}_{-}}{{T}_{2}},$$## (5)

$$\frac{\partial {D}_{o}}{\partial t}=\frac{{D}_{p}-{D}_{o}}{{T}_{1}}+({E}_{+}^{*}{P}_{+}+{E}_{-}^{*}{P}_{-}+c.c.),$$## (6)

$$\frac{\partial {D}_{2}}{\partial t}=({E}_{+}^{*}{P}_{-}+{E}_{-}{P}_{+}^{*}+c.c.)-\frac{{D}_{2}}{{T}_{g}},$$^{8}${D}_{p}$ is the steady-state normalized population inversion proportional to the pumping factor $p$ ($p\ge 1$), ${T}_{1}$ and ${T}_{2}$ are the longitudinal and transverse relaxation times, respectively, ${T}_{g}$ is the parameter that describes the strength of the SHB and it ranges from 0 (no SHB) to ${T}_{1}$ (the strongest SHB), ${l}_{0}$ is the linear cavity loss; $\tilde{\gamma}={\hslash}^{2}\gamma /{\mu}^{2}$, where $\gamma $ is the coefficient that describes the strength of the saturable absorber and $\mu $ is the transition dipole moment, and the parameter ${\beta}_{2}$ is the GVD coefficient.

The linearization of the dynamic Eqs. (1)–(6) is achieved through the addition of a small perturbation to the steady-state solution. Then, a set of equations with respect to the perturbations are obtained as

## (7)

$$\frac{\partial \delta {E}_{+}}{\partial t}=\frac{c}{n}[-\frac{\partial \delta {E}_{+}}{\partial z}-\delta {P}_{+}-\frac{1}{2}({l}_{0}-3\tilde{\gamma}{|\overline{E}|}^{2})\delta {E}_{+}-i\frac{{\beta}_{2}}{2}\frac{{\partial}^{2}\delta {E}_{+}}{\partial {t}^{2}}],$$## (8)

$$\frac{\partial \delta {E}_{-}}{\partial t}=\frac{c}{n}[+\frac{\partial \delta {E}_{-}}{\partial z}-\delta {P}_{-}-\frac{1}{2}({l}_{0}-3\tilde{\gamma}{|\overline{E}|}^{2})\delta {E}_{-}-i\frac{{\beta}_{2}}{2}\frac{{\partial}^{2}\delta {E}_{-}}{\partial {t}^{2}}],$$## (9)

$$\frac{\partial \delta {P}_{+}}{\partial t}=-\frac{1}{2}[\overline{E}(\delta {D}_{0}+\delta {D}_{2})+{\overline{D}}_{0}\delta {E}_{+}+{\overline{D}}_{2}\delta {E}_{-}]-\frac{\delta {P}_{+}}{{T}_{2}},$$## (10)

$$\frac{\partial \delta {P}_{-}}{\partial t}=-\frac{1}{2}[\overline{E}(\delta {D}_{0}+\delta {D}_{2}^{*})+{\overline{D}}_{0}\delta {E}_{-}+{\overline{D}}_{2}\delta {E}_{+}]-\frac{\delta {P}_{-}}{{T}_{2}},$$## (11)

$$\frac{\partial \delta {D}_{0}}{\partial t}=-\frac{\delta {D}_{0}}{{T}_{1}}+\overline{P}(\delta {E}_{+}+\delta {E}_{-}+c.c.)+\overline{E}(\delta {P}_{+}+\delta {P}_{-}+c.c.),$$## (12)

$$\frac{\partial \delta {D}_{2}}{\partial t}=-\frac{\delta {D}_{2}}{{T}_{g}}+\overline{P}(\delta {E}_{+}^{*}+\delta {E}_{-})+\overline{E}(\delta {P}_{+}^{*}+\delta {P}_{-}).$$Without loss of generality, we assumed that all the variables are complex numbers except the population inversion ${D}_{o}$. The steady-state solution for Eqs. (1–6) is obtained as

## (13)

$$\overline{E}={\overline{\rho}}_{E}=\sqrt{\frac{{l}_{0}{T}_{2}(2{T}_{1}+{T}_{g})-\tilde{\gamma}-\sqrt{{[\tilde{\gamma}+{l}_{0}{T}_{2}(2{T}_{1}+{T}_{g})]}^{2}-4{l}_{0}p\tilde{\gamma}{T}_{2}(2{T}_{1}+{T}_{g})}}{2{T}_{2}(2{T}_{1}+{T}_{g})\tilde{\gamma}}},$$## (16)

$${\overline{D}}_{0}={\overline{\rho}}_{0}=\frac{{l}_{0}^{\prime}}{{T}_{2}}-{\overline{\rho}}_{2},$$In order to decouple the behaviors of these variables into amplitude and phase domains, we adopt the polar expression for each variable, i.e., $X=\rho {\mathrm{e}}^{j\theta}$, where $X$ is any of the variables in Eqs. (1–6). Through some mathematical manipulations, the small perturbation is expressed in the polar form as

## (17)

$$\delta X(z,t)=[\delta \rho (z,t)+j\overline{\rho}(z)\delta \theta (z,t)]{\mathrm{e}}^{j\overline{\theta}(z)},$$Bringing Eqs. (17) and (18–21) into the differential equations in Eqs. (7–12) of the perturbations, the perturbations in amplitude and phase are decoupled into two equation sets as

## (22)

$$\frac{\partial \delta {\rho}_{+}^{F}}{\partial t}=\frac{c}{n}[-\frac{\partial \delta {\rho}_{+}^{F}}{\partial z}-\delta {\rho}_{+}^{P}-\frac{1}{2}({l}_{0}-3\tilde{\gamma}{|\overline{F}|}^{2})\delta {\rho}_{+}^{F}-i\frac{{\beta}_{2}}{2}\frac{{\partial}^{2}\delta {\rho}_{+}^{F}}{\partial {t}^{2}}],$$## (23)

$$\frac{\partial \delta {\rho}_{-}^{F}}{\partial t}=\frac{c}{n}[\frac{\partial \delta {\rho}_{-}^{F}}{\partial z}-\delta {\rho}_{-}^{P}-\frac{1}{2}({l}_{0}-3\tilde{\gamma}{|\overline{F}|}^{2})\delta {\rho}_{-}^{F}-i\frac{{\beta}_{2}}{2}\frac{{\partial}^{2}\delta {\rho}_{=}^{F}}{\partial {t}^{2}}],$$## (24)

$$\frac{\partial \delta {\rho}_{+}^{P}}{\partial t}=-\frac{1}{2}\overline{F}(\delta {\rho}_{0}+\delta {\rho}_{2})-\frac{1}{2}{\overline{\rho}}_{0}\delta {\rho}_{+}^{F}-\frac{1}{2}{\overline{\rho}}_{2}\delta {\rho}_{-}^{F}-\frac{\delta {\rho}_{+}^{P}}{{T}_{2}},$$## (25)

$$\frac{\partial \delta {\rho}_{-}^{P}}{\partial t}=-\frac{1}{2}\overline{F}(\delta {\rho}_{0}+\delta {\rho}_{2})-\frac{1}{2}{\overline{\rho}}_{0}\delta {\rho}_{-}^{F}-\frac{1}{2}{\overline{\rho}}_{2}\delta {\rho}_{+}^{F}-\frac{\delta {\rho}_{-}^{P}}{{T}_{2}},$$## (26)

$$\frac{\partial \delta {\rho}_{0}}{\partial t}=2\overline{F}(\delta {\rho}_{+}^{P}+\delta {\rho}_{-}^{P})+2\overline{P}(\delta {\rho}_{+}^{F}+\delta {\rho}_{-}^{F})-\frac{\delta {\rho}_{0}}{{T}_{1}},$$## (27)

$$\frac{\partial \delta {\rho}_{2}}{\partial t}=\overline{F}(\delta {\rho}_{+}^{P}+\delta {\rho}_{-}^{P})+\overline{P}(\delta {\rho}_{+}^{F}+\delta {\rho}_{-}^{F})-\frac{\delta {\rho}_{2}}{{T}_{g}}$$## (28)

$$\frac{\partial \delta {\theta}_{+}^{F}}{\partial t}=\frac{c}{n}[-\frac{\partial \delta {\theta}_{+}^{F}}{\partial z}-\frac{\overline{P}}{\overline{F}}\delta {\theta}_{+}^{P}-\frac{1}{2}({l}_{0}-3\tilde{\gamma}{|\overline{F}|}^{2})\delta {\theta}_{+}^{F}],$$## (29)

$$\frac{\partial \delta {\theta}_{-}^{F}}{\partial t}=\frac{c}{n}[\frac{\partial \delta {\theta}_{-}^{F}}{\partial z}-\frac{\overline{P}}{\overline{F}}\delta {\theta}_{-}^{P}-\frac{1}{2}({l}_{0}-3\tilde{\gamma}{|\overline{F}|}^{2})\delta {\theta}_{-}^{F}],$$## (30)

$$\frac{\partial \delta {\theta}_{+}^{P}}{\partial t}=-\frac{1}{2}\frac{\overline{F}}{\overline{P}}({\overline{\rho}}_{0}\delta {\theta}_{+}^{F}+{\overline{\rho}}_{2}\delta {\theta}_{-}^{F}+{\overline{\rho}}_{2}\delta {\theta}_{2})-\frac{\delta {\theta}_{+}^{P}}{{T}_{2}},$$## (31)

$$\frac{\partial \delta {\theta}_{-}^{P}}{\partial t}=-\frac{1}{2}\frac{\overline{F}}{\overline{P}}({\overline{\rho}}_{0}\delta {\theta}_{-}^{F}+{\overline{\rho}}_{2}\delta {\theta}_{+}^{F}-{\overline{\rho}}_{2}\delta {\theta}_{2})-\frac{\delta {\theta}_{-}^{P}}{{T}_{2}},$$## (32)

$$\frac{\partial \delta {\theta}_{2}}{\partial t}=\frac{\overline{P}\overline{F}}{{\overline{\rho}}_{2}}(\delta {\theta}_{-}^{F}+\delta {\theta}_{-}^{P}-\delta {\theta}_{+}^{F}-\delta {\theta}_{+}^{P})-\frac{\delta {\theta}_{2}}{{T}_{g}}.$$Equations (22–27) describe the amplitude perturbations, while Eqs. (28–32) are phase perturbations. Owing to the symmetry and antisymmetry relationships between each pair of forward and backward propagating variables, such as $({\rho}_{+}^{E},{\rho}_{-}^{E})$, $({\theta}_{+}^{E},{\theta}_{-}^{E})$, $({\rho}_{+}^{P},{\rho}_{-}^{P})$, and $({\theta}_{+}^{P},{\theta}_{-}^{P})$, the number of equations involved in the stability analysis can be reduced. In the Fabry–Perot cavity $(0,L)$, the symmetry and antisymmetry relationships about the cavity center $z=L/2$ are expressed as

From Eqs. (22–27) and (28–32), it can be found that the pairs $({\rho}_{+}^{E},{\rho}_{-}^{E})$ and $({\rho}_{+}^{P},{\rho}_{-}^{P})$ possess symmetry about the cavity center and the pairs $({\theta}_{+}^{E},{\theta}_{-}^{E})$ and $({\theta}_{+}^{P},{\theta}_{-}^{P})$ possess antisymmetry about the cavity center. Subsequently, those two sets of equations are reduced to Eqs. (35–38) and (39–41), respectively,

## (35)

$$\frac{\partial \delta {\rho}_{+}^{E}}{\partial t}=\frac{c}{n}[-\frac{\partial \delta {\rho}_{+}^{E}}{\partial z}-\delta {\rho}_{+}^{P}-\frac{1}{2}({l}_{0}-3\tilde{\gamma}|\overline{E}{|}^{2})\delta {\rho}_{+}^{E}-i\frac{{\beta}_{2}}{2}\frac{{\partial}^{2}\delta {\rho}_{+}^{E}}{\partial {t}^{2}}],$$## (36)

$$\frac{\partial \delta {\rho}_{+}^{P}}{\partial t}=-\frac{1}{2}\overline{E}(\delta {\rho}_{0}+\delta {\rho}_{2})-\frac{1}{2}({\overline{\rho}}_{o}+{\overline{\rho}}_{2})\delta {\rho}_{+}^{E}-\frac{\delta {\rho}_{+}^{P}}{{T}_{2}},$$## (37)

$$\frac{\partial \delta {\rho}_{0}}{\partial t}=4\overline{P}(\delta {\rho}_{+}^{E})+4\overline{E}(\delta {\rho}_{+}^{P})-\frac{\delta {\rho}_{0}}{{T}_{1}},$$## (38)

$$\frac{\partial \delta {\rho}_{2}}{\partial t}=2\overline{P}(\delta {\rho}_{+}^{E})+2\overline{E}(\delta {\rho}_{+}^{P})-\frac{\delta {\rho}_{2}}{{T}_{g}}$$## (39)

$$\frac{\partial \delta {\theta}_{+}^{E}}{\partial t}=\frac{c}{n}[-\frac{\partial \delta {\theta}_{+}^{E}}{\partial z}-\frac{\overline{P}}{\overline{E}}\delta {\theta}_{+}^{P}-\frac{1}{2}({l}_{0}-3\tilde{\gamma}{|\overline{E}|}^{2})\delta {\theta}_{+}^{E}],$$## (40)

$$\frac{\partial \delta {\theta}_{+}^{P}}{\partial t}=-\frac{1}{2}\frac{\overline{E}}{\overline{P}}[({\overline{\rho}}_{0}-{\overline{\rho}}_{2})\delta {\theta}_{+}^{E}+{\overline{\rho}}_{2}\delta {\theta}_{2}]-\frac{\delta {\theta}_{+}^{P}}{{T}_{2}},$$## (41)

$$\frac{\partial \delta {\theta}_{2}}{\partial t}=\frac{\overline{P}\overline{E}}{{\overline{\rho}}_{2}}(-2\delta {\theta}_{+}^{E}-2\delta {\theta}_{+}^{P})-\frac{\delta {\theta}_{2}}{{T}_{g}}.$$The amplitude stability analysis is based on Eqs. (35–38), while phase stability analysis is based on Eqs. (39–41). Both equation sets involve a second-order differential term that is brought by the GVD effect, which results in the quadratic eigenvalue problem in the stability analysis.

## 3.

## Stability Analysis Through the Quadratic Eigenvalue Problem

Dynamic behaviors of amplitude and phase are described by equation sets (35–38) and (39–41), respectively. Examination of stability is conducted by solving the eigenvalues for each equation set. In this case, we need to deal with the quadratic eigenvalue problem that results from the GVD effect. We take equation set (35–38) for amplitude stability as the example to illustrate our solution procedure for the quadratic eigenvalue problem.

We consider transforming the second-order differential equation to a first-order differential equation; thus, the quadratic eigenvalue problem will be degenerated to a linear eigenvalue problem that we are familiar with. For equation set (35–38), by defining an additional variable, $[\partial (\delta {\rho}_{+}^{E})]/\partial t$, the second-order differential equation set (35–38) is linearized (or reduced to the first-order differential equation) in the matrix form as

## (42)

$$\frac{\partial}{\partial t}\left[\begin{array}{c}\begin{array}{l}\frac{\partial (\delta {\rho}_{+}^{E})}{\partial t}\end{array}\\ \begin{array}{l}\delta {\rho}_{+}^{E}\end{array}\\ \begin{array}{l}\delta {\rho}_{+}^{P}\\ \end{array}\\ \begin{array}{l}\delta {\rho}_{o}\\ \end{array}\\ \delta {\rho}_{2}\end{array}\right]=\left[\begin{array}{ccccc}\begin{array}{l}\frac{2n}{{\beta}_{2}c}i\\ \end{array}& \begin{array}{l}\frac{2i}{{\beta}_{2}}[\frac{\partial}{\partial z}+\frac{1}{2}({l}_{0}-3\tilde{\gamma}{|\overline{E}|}^{2})]\end{array}& \begin{array}{l}\frac{2i}{{\beta}_{2}}\\ \end{array}& \begin{array}{l}0\\ \end{array}& \begin{array}{l}0\\ \end{array}\\ \begin{array}{l}1\end{array}& \begin{array}{l}0\\ \end{array}& \begin{array}{l}0\end{array}& \begin{array}{l}0\\ \end{array}& \begin{array}{l}0\\ \end{array}\\ 0& -\frac{{\overline{\rho}}_{o}+{\overline{\rho}}_{2}}{2}& -\frac{1}{{T}_{2}}& -\frac{\overline{E}}{2}& -\frac{\overline{E}}{2}\\ 0& 4\overline{P}& 4\overline{E}& -\frac{1}{{T}_{1}}& 0\\ 0& 2\overline{P}& 2\overline{F}& 0& -\frac{1}{{T}_{g}}\end{array}\right]\left[\begin{array}{c}\begin{array}{l}\frac{\partial (\delta {\rho}_{+}^{E})}{\partial t}\\ \end{array}\\ \begin{array}{l}\delta {\rho}_{+}^{E}\\ \end{array}\\ \begin{array}{l}\delta {\rho}_{+}^{P}\\ \end{array}\\ \begin{array}{l}\delta {\rho}_{o}\\ \end{array}\\ \delta {\rho}_{2}\end{array}\right].$$The stability of the dynamic equation set (35–38) is studied by finding the eigenvalues of the matrix in Eq. (42). Assume that $\lambda $ is any of the eigenvalues and with some mathematical manipulations based on Eq. (42), we arrive at

where## (44)

$${f}_{a}(\lambda )=-\frac{{\beta}_{2}}{2}i{\lambda}^{2}-\frac{n}{c}\lambda -\frac{1}{2}({l}_{o}-3\overline{\gamma}{|\overline{E}|}^{2})+\frac{({\overline{\rho}}_{o}+{\overline{\rho}}_{2})(\lambda +{T}_{1}^{=1})(\lambda +{T}_{g}^{-1})+2\overline{E}\overline{P}(\lambda +{T}_{1}^{=1})+4\overline{E}\overline{P}(\lambda +{T}_{g}^{-1})}{2[(\lambda +{T}_{1}^{=1})(\lambda +{T}_{2}^{-1})(\lambda +{T}_{g}^{-1})+{\overline{E}}^{2}(\lambda +{T}_{1}^{=1})+2{\overline{E}}^{2}(\lambda +{T}_{g}^{-1})]}.$$The value for ${f}_{a}(\lambda )$ can be figured out from the boundary conditions of the Fabry–Perot cavity. As the mirror reflection at the two facets has been accounted for in the linear loss ${l}_{o}$ in Eqs. (1) and (2), the boundary conditions become formally identical to those ideal Fabry–Perot resonators with perfectly reflecting mirrors, i.e., ${\rho}_{+}^{E}(0)={\rho}_{-}^{E}(0)$ and ${\rho}_{+}^{E}(L)=\phantom{\rule{0ex}{0ex}}{\rho}_{-}^{E}(L)$. Combining these boundary conditions with the symmetry relationship between the forward and backward propagations stated in Eq. (33), we have ${\rho}_{+}^{E}(0)={\rho}_{+}^{E}(L)$, and then,

After substituting Eq. (46) into Eq. (45), we get the transcendental equation for the entire spectrum of eigenvalues ${\lambda}_{k}$ as

where $k$ is any integer number, i.e., $k=0,\pm 1,\pm 2,\pm 3,\dots $.The stability of dynamic behaviors in the amplitude described in Eq. (42) is characterized by the eigenvalues ${\lambda}_{k}$. Each integer number $k$ represents a transverse wave mode with an offset frequency of $\mathrm{\Delta}f$ from the central frequency ${f}_{o}$, and $\mathrm{\Delta}f=(c/2nL)k$, where it implies that the separation between the two adjacent wave modes is the round trip frequency along the cavity. For each mode, if the maximum real part of all the eigenvalues (the parametric gain) is negative, then the mode is stable; otherwise, the system is unstable when the parametric gain is positive or the system is marginally stable with zero parametric gain.

We follow a similar procedure to analyze phase stability. To avoid redundancy, here, we only present the major steps in the derivation procedure. The linearized dynamic equation for the phase domain based on Eqs. (39–41) is

## (48)

$$\frac{\partial}{\partial t}\left[\begin{array}{c}\begin{array}{l}\frac{\partial \delta {\theta}_{+}^{E}}{\partial t}\end{array}\\ \begin{array}{l}\delta {\theta}_{+}^{E}\end{array}\\ \begin{array}{l}\delta {\theta}_{+}^{E}\end{array}\\ \delta {\theta}_{2}\end{array}\right]=\left[\begin{array}{cccc}\frac{2n}{{\beta}_{2}c}i& \frac{2i}{{\beta}_{2}}[\frac{\partial}{\partial z}+\frac{1}{2}({l}_{0}-3\tilde{\gamma}|\overline{E}{|}^{2})]& \frac{2\overline{P}}{{\beta}_{2}\overline{E}}& 0\\ 1& 0& 0& 0\\ \begin{array}{l}0\end{array}& \begin{array}{l}\frac{({\overline{\rho}}_{2}-{\overline{\rho}}_{o})\overline{E}}{2\overline{P}}\end{array}& \begin{array}{l}-\frac{1}{{T}_{2}}\end{array}& \begin{array}{l}-\frac{{\overline{\rho}}_{2}E}{2\overline{P}}\end{array}\\ 0& -\frac{2\overline{E}\overline{P}}{{\overline{\rho}}_{2}}& -\frac{2\overline{E}\overline{P}}{{\overline{\rho}}_{2}}& -\frac{1}{{T}_{g}}\end{array}\right]\left[\begin{array}{c}\begin{array}{l}\frac{\partial \delta {\theta}_{+}^{E}}{\partial t}\\ \end{array}\\ \begin{array}{l}\delta {\theta}_{+}^{E}\\ \end{array}\\ \begin{array}{l}\delta {\theta}_{+}^{E}\\ \end{array}\\ \delta {\theta}_{2}\end{array}\right].$$With any eigenvalue $\lambda $ of the above matrix, we have the following relationship about the phase angle of the electric field at the two boundaries of the cavity,

where## (50)

$${f}_{p}(\lambda )=-\frac{{\beta}_{2}}{2}i{\lambda}^{2}-\frac{n}{c}\lambda -\frac{1}{2}({l}_{o}-3\gamma {|\overline{E}|}^{2})-\frac{({\overline{\rho}}_{2}-{\overline{\rho}}_{o})(\lambda +{T}_{g}^{-1})+2\overline{E}\overline{P}}{2[(\lambda +{T}_{2}^{-1})(\lambda +{T}_{g}^{-1})-{\overline{E}}^{2}]}.$$The above transcendental equation is solved by incorporating the boundary conditions of the phase angle. The following relationship holds because of the reflection of the phase angle at each of the boundaries and the antisymmetry between the forward and backward propagations of phase angles

Thus, we haveThe full spectrum of eigenvalues ${\lambda}_{k}$ for phase stability is solved based on Eq. (52) for each mode identified by the integer $k$. The stability criteria are determined by the parametric gain as explained previously.

## 4.

## Results and Discussion

We aim to investigate the amplitude and phase stabilities under the GVD and saturable absorber effects. Our ultimate goal is to tailor these effects to design a QCL for MIR pulse generation, so we vary the GVD and saturable absorber strength within a certain range by referring to the parameters provided in the literature^{8}^{,}^{10} to examine the interaction between the GVD and the saturable absorber. The saturable absorber strength $\gamma $ of the QCL structure (wafer No. 3251) emitting at 8.47 *μ*m is between $1\sim 2\times {10}^{-9}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{cm}/{\mathrm{V}}^{2}$ when the active region width is in the range of 3–7.5 *μ*m.^{10} The saturable absorber strength is estimated based on ${n}_{2}=2\times {10}^{-8}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{2}/\mathrm{W}$.^{10} We take this structure as our simulation prototype and set the variation range of the saturable absorber strength $\gamma $ to be $1\sim 4\times {10}^{-9}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{cm}/{\mathrm{V}}^{2}$. The waveguide structures of the QCL in Refs. 8 and 10 (the waveguide structure of Ref. 10 is detailed in Ref. 18) are similar. They are both grown by metal organic vapor-phase epitaxy (MOVPE) and made by an InGaAs/InAlAs core and highly Fe-doped InP cladding layers. Thus, for our simulation, we can adopt the GVD coefficient of the structure in Ref. 8 as the baseline of our GVD variation. As reported in Ref. 8, the measurement of a $\sim 5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{m}$ QCL with an InGaAs/InAlAs active region shows an anomalous GVD coefficient ${\beta}_{2}\sim -4.6\times {10}^{-6}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{ps}}^{2}/\mu \text{m}$.^{8} We set the variation range of the GVD coefficient ${\beta}_{2}$ to be $(-10\sim -1)\times {10}^{-6}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{ps}}^{2}/\mu \text{m}$. For other parameters, we follow the parameters provided for the Fabry–Perot cavity in Ref. 10.

We studied the amplitude and phase instabilities with both the GVD and saturable absorber existing in the cavity under various pumping strengths. Results are presented in Fig. 1 for amplitude instability and in Fig. 2 for phase instability. We want to note that the line type and color assignment shown in the insert of Fig. 1(a) also apply to all other figures. In each case, the stability is demonstrated by the parametric gain versus detuning frequency. For reference, in each figure, we also present the cases when the GVD or saturable absorber is absent, i.e., ${\beta}_{2}=0$ or $\gamma =0$. In our simulation procedure, we found that the amplitude (or phase instability) responds almost indifferently to the strength of the GVD when ${\beta}_{2}$ has the range of $(-10\sim -1)\times \phantom{\rule{0ex}{0ex}}{10}^{-6}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{ps}}^{2}/\mu \text{m}$, which is a reasonable range for the QCLs (Refs. 8 and 17). Thus, we only show the results with ${\beta}_{2}={\beta}_{20}=-5\times {10}^{-6}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{ps}}^{2}/\mu \text{m}$ to represent the cases with the GVD effect.

The amplitude instability without the GVD effect exhibits as that of an RNGH-type instability. From Fig. 1, it can be found that, under each pumping level, the presence of the GVD alters the influence of the saturable absorber on the amplitude instability. That is, an increase in the saturable absorber strength can drive the amplitude to be more stable when the GVD is present. Inversely, a stronger saturable absorber boosts the instability by lowering the instability threshold when the GVD is absent.^{10}^{,}^{12} Thus, the lasing cavity with a stronger saturable absorber effect is less affected by the presence of the GVD. The interaction between the GVD and the saturable absorber on the amplitude instability is also affected by the pumping strength $p$. Figure 1(a) tells that, under low pumping strength, the GVD tends to destabilize the amplitude and flattens the parametric gain curve throughout all frequency modes. When the pumping strength gets stronger, as shown in Figs. 1(b), 1(c), and 1(d), the effect of the GVD can be suppressed by the saturable absorber effect. Hence, the instability without the GVD resembles that of the instability with the GVD, especially for the frequency portion that is closely around the central frequency. Apparently, to achieve such suppression, the weaker the saturable absorber is, the higher is the pumping strength needed. This is because of the fact that the suppression is related to the nonlinear portion of the loss $\tilde{\gamma}{|\overline{E}|}^{2}$. We notice that, in Figs. 1(c) and 1(d), for the case with the highest saturable absorber strength under study, the symmetry of the gain spectrum about the central frequency is broken, which is mathematically caused by the emergence of an imaginary part in the steady-state solution.

The phase instability is of a single-mode nature. The GVD has a lesser effect on the phase instability than it does on amplitude instability, especially, the GVD effect is almost suppressed at a higher pumping level. At a lower pumping level, the destabilization of the phase by the GVD only happens at frequencies that are away from the central frequency.

## 5.

## Conclusions

We present an analysis procedure on the amplitude and phase instabilities of MIR QCLs with the presence of GVD and saturable absorber effects in the lasing cavity. The significance of this work is in accounting for the GVD effect, whose presence in MIR QCLs has been substantiated, in the discussion of QCL instability mechanisms.

Stability analysis was performed through the linearization of dynamic equations based on the Maxwell–Bloch formalism. Because of the second-order differential term brought by the GVD, the quadratic eigenvalue problem results from the linearization and it is treated by a degeneration technique. The instability mechanisms are decoupled into amplitude and phase domains through the symmetry or antisymmetry propagation of fields along the cavity.

Simulation results show that the effect of the GVD on the stability of QCLs is largely dependent on the saturable absorber and pumping strength. The GVD could significantly destabilize the QCL pulse amplitude under a lower pumping level and weaker saturable absorber strength. With the increase of pumping and saturable absorber strength, the GVD effect could be suppressed. The GVD has only a slight effect on the phase instability, and the effect decreases as the frequency approaches the central frequency. The results obtained from this research will contribute to our further investigation on the topic of soliton dynamics of QCLs for MIR pulse generation.

## Acknowledgments

This work was supported in part by the National Science Foundation (NSF) by Grant No. ECCS-1232273. The authors would like to thank Dr. Qijie Wang in Nanyang Technological University, Singapore, for helpful discussions.

## References

## Biography

**Jing Bai** is an associate professor in the Department of Electrical Engineering at the University of Minnesota Duluth (UMD). She received her PhD degree and MSc degree in electrical and computer engineering at Georgia Institute of Technology in 2007 and 2003, respectively. She also earned the MEng degree in Nanyang Technological University, Singapore, in 1999 and the BEng degree in Tsinghua University, China, in 1996, both in mechanical engineering. Her studies focus on theoretical investigations on intracavity nonlinearities and dynamic behaviors of nanoscale optoelectronic structures. She also extends her research interests to the biomedical applications of optoelectronic technologies.

**Debao Zhou** joined the Faculty of Mechanical and Industrial Engineering at UMD in January 2009. Prior to that, he was a McKnight professor in the Department of Electrical Engineering at UMD. From January 2005 to January 2008, he was with Georgia Tech Research Institute as a postdoctoral fellow. He received his PhD degree in mechanical and production engineering from Nanyang Technological University, Singapore, in 2004. He earned his MEng degree and BEng degree from Tsinghua University, China, both in mechanical engineering. His research interests include system dynamics and control, distributed force sensor design for biomedical applications, and image processing techniques for intelligent transportation systems.