Analytical theory of frequency-modulated combs: generalized mean-field theory, complex cavities, and harmonic states

Levi Humbard; David Burghoff

doi:10.1364/OE.445570

1. Introduction

Frequency-modulated (FM) combs have long existed in lasers [1], and in particular tend to arise in semiconductor lasers, which can spontaneously enter self-FM regimes [2–10]. There are many applications of laser frequency combs, particularly dual-comb spectroscopy [11,12] and radiometry [13]. Frequently, these combs exhibit an interesting linearly-chirped behavior [5–7], which was not theoretically predicted. This behavior was first conclusively identified in quantum cascade lasers (QCLs) [5], which are notable for their fast gain recovery times (gain recovery times shorter than the cavity round trip time, typically a few picoseconds in mid-IR QCLs), but have also been identified in slow gain recovery time systems (system with gain recovery times longer than the cavity round trip time) as well [7–10].

Recently, we developed a mean-field theory (MFT) analogous to the Lugiato-Lefever equation that is essential for understanding the universality of this behavior [14]. The Lugiato-Lefever equation has been used to describe cubic nonlinear properties, such as Kerr combs, in various resonators and materials [15–18]. The basic idea is that if the changes to the field are small within a single round-trip, mean-field theories allow one to quickly predict the long term behavior of a system without taking small timesteps. The fundamental difficulty with lasers is that most lasers do not have small changes within a round trip: they often have large mirror losses that are balanced by gain. Thus, most simulations that had previously been developed relied either on explicit solutions to the Maxwell-Bloch equations [3,19–27] or simpler coupled master equations [28,29]. These methods require the simulation of many timesteps and in general are not analytically tractable. To address this limitation, we showed that the large changes within a laser could be eliminated by normalizing the electric field to a scaling function that essentially preconditions the master equations [14,30]. This scaling function is found using the steady-state power profile of the laser, i.e., the power produced by a laser without any phase-dependent terms, the continuous wave (CW) solution. Note that this normalization is essential for most Fabry-Perot semiconductor lasers (whose typical mirror reflectivities are around 30%), but it is not needed for ring lasers (whose coupling losses are much smaller). Thus, suitable mean-field theories can be derived for rings by setting this scale factor to unity [31].

Although active cavity mean-field theory can apply to several different systems [14,31,32], in this work we focus on the system shown in Fig. 1 (data for this figure and the others are available online, as we show in Code 2 (Ref. [45])). This setup describes a two-level laser gain medium inside a Fabry-Perot (FP) cavity. If one then makes some weak simplifying assumptions—a slowly varying envelope, an amplitude that closely mimics the CW solution, linear internal power growth, and infinite gain bandwidth (we later show how this method can be modified to account for weak gain curvature, i.e. limited gain bandwidth)—the analytical result follows in the form of a phase-dominated nonlinear Schrödinger equation (phase NLSE) [14],

(1)$$-i \frac{\partial F}{\partial T} = \frac{\beta}{2} \frac{\partial^{2} F}{\partial z^{2}} + \gamma \left|F\right|^{2} \left(\arg F-\langle \arg F \rangle\right) F +i r \left(\left|F\right|^{2}-P_0\right)F,$$

where $F$ is the electric field normalized to the unitless intracavity gain (proportional to the laser’s steady-state intracavity intensity), T is the slow time, $\beta$ is the normalized group velocity dispersion (GVD), $\gamma$ is the nonlinear cross-steepening, angle brackets represent an average over position, and $r$ represents amplitude relaxation to the steady state intensity at the left facet, $P_0$. This model is surprisingly robust, applying to both fast gain media lasers and to the CW-like solutions of slow gain media lasers. The position $z$ is taken within the artificially extended laser cavity with its origin ($z=0$) position fixed at the left facet—if $L_c$ is the physical cavity length, wave traveling in the positive direction are are mapped to $z \in \left [0,L_c\right )$, and waves traveling in the negative direction are mapped to $z \in \left [-L_c,0\right )$. Previously, it was shown that the fundamental solution to this equation can be written analytically as [14]

(2)$$F_1\left(z,T\right) =A_0 \exp\left[i \frac{\gamma \left|A_0\right|^{2}}{2 \beta} \left(z^{2}-\frac{1}{3}L_c^{2} \gamma \left|A_0\right|^{2} T\right) \right],$$

where $A_0=\sqrt {\frac {P_0}{1+\gamma /2r}}$. This solution was referred to as an extendon due to its spatially-extended nature. The higher-order solutions to the phase NLSE are harmonic states [27,33], and the $N$th one can be taken as [14]

(3)$$F_N\left(z,T\right) = F_1\left(z\bmod \frac{2L_c}{N},T\right),$$

where the modulo operator is taken over $z\in \left [-\frac {L_c}{N},\frac {L_c}{N}\right )$. Fabry-Perot lasers can sometimes prefer harmonic modes of operation, but even when the fundamental state is preferred the harmonic states detailed in Eq. (3) are useful for forming a basis set for further analysis.

Fig. 1. Diagram of a Fabry-Perot laser, showing the intensity of waves traveling in each direction as solid lines. The intensity closely mimics the steady-state power (i.e., the continuous wave solution), which are shown here as dashed lines. These quantities are plotted for the case of $R_1=0.2$ and $R_2=0.45$.

Download Full Size | PDF

In this work, we extend these results to derive several key results pertaining to FM combs. First, we generalize the previous mean-field theory from Ref. [14], showing how it can be applied to arbitrary models of the laser under consideration—two-level models, three-level models, systems with defects, etc. An algorithm is provided describing the steps one needs to follow to convert an arbitrary master equation into its corresponding mean-field theory.

Next, we generalize previous results to allow for arbitrary cavities. In order to derive the phase NLSE, it was previously assumed that one of the end mirrors of the Fabry-Perot cavity have unity reflectivity (R=100%). It was not shown that the phase NLSE applies to cavities with two non-unity reflectors, which includes most lasers. An analysis is presented here to expand this previous result. This is a useful pursuit because the mirror losses are integral for optimizing the the spectral bandwidth [34].

Then, we improve the accuracy of the extendon theory in the face of gain curvature (or equivalently, a situation with limited gain bandwidth). Though the analytical solution of the NLSE (Eq. (1)) matches the numerical solution while neglecting gain curvature (everywhere but the turnaround point in fast gain media, can have additional discrepancies in slow gain media), when more realistic gain curvature was applied the accuracy dips noticeably [14]. To address this, we use a perturbative analysis to show the relevant parameters can be modified when gain curvature is present. These results allow the numerical and analytical results to match closely in the case of moderate gain curvature, though at very high values the accuracy drops (these corrections are good for fast gain media, but will eventually breakdown for slow enough gain media). We use this result to analyze the effect of a linewidth enhancement factor, showing analytically and numerically that the chirp of the field is unchanged.

Lastly, we show that by calculating the diagonal elements of the matrix corresponding to Eq. (1), one can gain insight into the dynamics of the phase NLSE. These calculations provide the carrier envelope offset (CEO) of the comb—the offset of the frequency comb teeth. This method also provides limited insights into the stability of the various harmonic states. In particular, we show that the addition of reflective defects—which has been shown to be an effective way of modifying harmonic frequency combs [35]— adds to the imaginary component of these diagonal elements, which can improve the stability of certain eigenstates while weakening the others. For sufficiently large defect magnitudes, this results in the laser favoring harmonic states rather than fundamental states, or for moderate magnitudes can generate quasi-stable solutions which breathe periodically in time.

2. Analysis

2.1 Laser mean field theory for arbitrary sets of master equations

We wish to expand and generalize the previous work’s mean-field theory analysis. The prior theory was derived specifically for the set of laser master equations describing a two-level system, but other formulations are possible depending on the physics one wishes to consider. For example, in section 2.5 we will consider the case of adding reflective defects to the cavity. In this section, we will closely examine how one can convert arbitrary master equation terms into appropriate mean-field terms. To begin, consider a two-level system, whose master equations can mostly simply be expressed as

(4)$$\begin{aligned} \frac{n}{c} \frac{\partial E_\pm}{\partial t} \pm \frac{\partial E_\pm}{\partial z} &= \frac{g_0}{2}\left(1-\frac{1}{P_s}\left(\left|E_\pm\right|^{2}+2\left|E_\mp\right|^{2}\right) \right) E_\pm{-}\frac{\alpha_w}{2} E_\pm\\ & +i\frac{1}{2}k^{\prime\prime}\frac{\partial^{2} E_\pm}{\partial t^{2}} + \frac{g_0}{2}T_2^{2} \frac{\partial^{2} E_\pm}{\partial t^{2}} \\ &+\frac{g_0}{2 P_s}\left(\left(2T_1+3T_2\right) \frac{\partial E_\mp^{{\ast}}}{\partial t} E_\mp{+}\left(T_1+\tfrac{5}{2}T_2\right) E_\mp^{{\ast}} \frac{\partial E_\mp}{\partial t} \right)E_\pm, \end{aligned}$$

where $E_\pm$ is the envelope of the wave traveling in the $\pm z$ direction (see Fig. 1 for reference), $g_0$ is the small signal gain, $k''$ is the group velocity dispersion (GVD), $T_1$ and $T_2$ are the population and coherence lifetimes for a two-level system, and $n$ is the refractive index of the medium (see Appendix E for the full version that includes all nonlinear terms). The first line represents losses, gain, and gain saturation (an approximation valid near the lasing threshold), the second represents dispersion and gain curvature, and the final lines cross-steepening. This model is valid provided that $T_1 \partial _t<<1$, which automatically holds for fast gain media. It also holds for slow gain media as long as they are not producing pulses, for example through the addition of a saturable absorber.

As previously discussed, the most critical terms for FM comb generation are the cross-steepening terms, which arise due to temporal variation of the gain grating (the spatially-varying gain created by gain saturation, also known as spatial hole burning). Other effects can be included as well, but these serve only to alter the parameters shown here. (For example, Kerr nonlinearity and linewidth enhancement shift the effective dispersion [14].) Alternatively, the master equations of a three level system can be described by [36]

(5)$$\begin{aligned} \left(\frac{n}{c} \partial_t \pm \partial_z \right)E_\pm & = \frac{g}{2} \frac{1+ib}{1+i\xi} \left[1-\tilde{T_2} \partial_t+\tilde{T_2}^{2} \partial_t^{2} \right]E_\pm\\ & -\frac{g T_g}{T_1 I_\textrm{sat}} \frac{1+ib}{1+i\xi}\left[|E_\mp|^{2} E_\pm{-} \left(\tilde{T_2}+T_g\right)|E_\mp|^{2} \partial_t E_\pm \right.\\ & \left. -\left(\tilde{T_2}+T_1\right)E_\pm E_\mp \partial_t E_\mp^{*}-\tilde{T_2} E_\pm E_\mp^{*} \partial_t E_\mp\right]-\frac{\alpha_w}{2}E_\pm, \end{aligned}$$

where $T_1$, $T_2$, and $T_g$ are the recovery times of the gain, polarization, and population grating, $\alpha _w$ is the waveguide loss, $g = \frac {g_0}{1+\frac {I}{I_{sat}}}$ is the saturated gain, $I_\textrm {sat}$ is saturated intensity, and $I=|E_-|^{2}+|E_+|^{2}$ the normalized intensity. The parameter $b$ describes the gain shape asymmetry and $\xi \left (b\right )$ and $\tilde {T}_2\left (b\right )$ are functions of containing contributions of the Bloch gain, which gives rise to a linewidth enhancement factor.

In order to unify these descriptions, note that Eq. (4) can be written more generally as

(6)$$\frac{n}{c}\frac{\partial{E_\pm}}{\partial{t}} \pm \frac{\partial E_\pm}{\partial z}=\sum_i p_i(z)\widehat{Q}_i(E_\pm(z,t))\widehat{R}_i(E_\mp(z,t)),$$

where $p_i(z)$ is a time invariant function of position and $\widehat {Q}_i$ and $\widehat {R}_i$ are operators acting on $E_\pm$ and $E_\mp$, respectively. This is valid because every term in the right-hand side of (4) is separable. Though (5) is not separable, it can be written in the same way provided that one uses a series expansion to take $g={}g_0\left (1-\frac {I}{I_\textrm {sat}} +\cdots \right )$. This is the general form we will therefore consider.

Next, backward propagating waves are flipped and extended to $\left (-L_c,0\right )$, where $L_c$ is the cavity length. This artificial cavity extension reduces the master equations to a single forward propagating wave on $\left (-L_c, L_c\right )$, and periodic boundary conditions now apply. Mathematically, this corresponds to making the substitutions $E_\pm (z,t) \rightarrow E(z,t)$ and $E_\mp (z,t) \rightarrow E(-z,t)$. For convenience, we adopt the notation $E_-(z)\equiv E(-z)$. The only modification needed to the master equation (aside from changing the domain of the fields) is to replace explicit space dependence ($p_i(z$)) with an appropriately modified piecewise-defined function, $\tilde {p}_i(z)$. After this, the steady-state intensity $P\left (z\right )$ of the system with no phase dynamics is then calculated, taking into account gain saturation, waveguide loss, and mirror losses. This intensity is typically very similar to the master equation solutions (see Fig. 1), and so normalizing by it eliminates large changes due to the facets. Mathematically, we define the slow-intensity envelope $F(z,t)$ using $E\left (z,t\right )\equiv K^{1/2}\left (z\right ) F\left (z,t\right )$. Here $K\left (z\right )$ can be thought of as the dimensionless power gain that the wave experiences as it traverses the cavity. It is time-independent, and for a laser above threshold it is related to the steady-state intensity by $K\left (z\right )=P\left (z\right )/P_0$, where $P_0$ is the power at the left cavity facet. It can alternatively be expressed in terms of the effective gain as

(7)$$K(z) = \exp\left({\int_0^{z} g_\textrm{eff}(z')dz'}\right),$$

where $g_\textrm {eff}$ contains the effects of gain saturation, waveguiding loss, and mirror losses.

(8)$$\begin{aligned} \frac{n}{c}\frac{\partial{F}}{\partial{t}}+\frac{\partial F}{\partial z}&=\sum_{i} \frac{\tilde{p}_i(z)}{\sqrt{K(z)}}\,\widehat{Q}_i\left(\sqrt{K(z)}F(z,t)\right)\widehat{R}_i\left(\sqrt{K_-(z)}F_-(z,t)\right)\\ &- \sum_{i\in P} \frac{\tilde{p}_i(z)}{\sqrt{K(z)P_0}}\,\widehat{Q}_i\left(\sqrt{K(z) P_0}\right)\widehat{R}_i\left(\sqrt{K_-(z)P_0}\right)F(z,t) \end{aligned}$$

where the first line contains all terms from the master equation and the second contains just those terms contained within the power calculation. Most terms included in the power calculation will now vanish. For example, for (4) one obtains

(9)$$\begin{aligned} \frac{n}{c} \frac{\partial F}{\partial t} + \frac{\partial F}{\partial z} = &-\frac{g_0}{2P_s}\left(K \left(\left|F\right|^{2}-P_0\right)+2K_-\left(\left|F_-\right|^{2} - P_0\right) \right) F +i\frac{1}{2}k^{\prime\prime}\frac{\partial^{2} F}{\partial t^{2}} + \frac{g_0}{2}T_2^{2} \frac{\partial^{2} F}{\partial t^{2}}\\ &+\frac{g_0}{2 P_s}K_-\left((2T_1+3T_2) \frac{\partial F_-^{{\ast}}}{\partial t} F_-{+}(T_1+\tfrac{5}{2}T_2) F_-^{{\ast}} \frac{\partial F_-}{\partial t} \right)F. \end{aligned}$$

The mirror losses, waveguide losses, and small-signal gain have vanished; the gain saturation is now driving the system towards a power of $P_0$. Again, (8) can be compactly written as

(10)$$\frac{n}{c}\frac{\partial{F}}{\partial{t}}+\frac{\partial F}{\partial z}=\sum_i f_i(z)\widehat{G}_i(F(z,t))\widehat{H}_i(F_-(z,t)),$$

where $\widehat {G}_i$ and $\widehat {H}_i$ are operators acting on $F$ and $F_-$, respectively, and $f_i(z)$ is a time invariant function of position.

At this point, one must integrate Eq. (10) over a full round trip over the cavity to finish the mean field derivation. The left-hand side of Eq. (10) is the typical wave equation and is periodic, so it vanishes with averaging. Hence, the only relevant round trip change comes from the sum of the right-hand side of the equation over every point during a round trip of the cavity. Therefore, (10) becomes

(11)$$\frac{\partial F}{\partial T} = \frac{\Delta F}{T_r} =\frac{c}{2L_cn}\sum_i \int_{{-}L_c}^{L_c} f_i\left(z+z'\right)\widehat{G}_i\left(F\left(z\right)\right)\widehat{H}_i\left(F\left({-}z-2z'\right)\right)dz',$$

where $T_r$ is the round trip time (the total time it takes for the wave to travel the extended cavity), $T$ is a slow time, usually described in fractional units of $T_r$, and $z\in \left (-L_c,L_c\right )$. To understand why each function is offset differently in Eq. (11), one should consult Fig. 2. For this averaging, assume the origin is attached to a point on the forward traveling envelope ($F$). After a short time of $z'/v_g$ passes, stationary functions such as $f\left (z\right )$ shift to left $\left (f\left (z\right )\rightarrow f\left (z+z'\right )\right )$ and backward traveling envelopes change according to $F_-\left (z\right )=F\left (-z\right )\rightarrow F\left (-z-2z'\right )=F_-\left (z+2z'\right )$. The second case involving backward envelopes is not quite as intuitive as the first, but if one considers that the forward envelope (origin) moves to the right z’ and the backward envelope moves to the left $z'$, it is apparent the overall shift is $2z'$. Armed with this information, one just needs to form the convolution integral to deal with backward propagating waves. The convolution integral is derived from the integral in (11) by applying the substitution $u=z+z'$, noting the argument of f is flipped since we are considering backward propagating, and arguing that since the integration still covers the whole extended cavity the specific bounds are arbitrary provide they cover an interval of $2L_c$. With this tool it is then straightforward to apply the generalized mean-field theory. Defining this integral as

(12)$$\widetilde{f}\left[q\right](z)=\frac{1}{2L_c} \int_{0}^{2L_c}f\left({-}u\right) q\left(z-2u\right)du.$$

With Eq. (12), the final result can be written compactly as

(13)$$\begin{aligned} \frac{n}{c} \frac{\partial F}{\partial T}=\sum_i \widehat{G}_i\left(F\right)\widetilde{f}[\widehat{H}_i\left(F\right)]. \end{aligned}$$

For convenience, many common terms converted from master equations to mean-field theory are detailed in Table (1). If desired, the spatiotemporal substitution [17] $\partial _t \rightarrow -\frac {c}{n} \partial _z$ can be used to convert fast-time derivatives to spatial derivatives, or vice versa.

Fig. 2. How the forward traveling envelope ($F$) changes relative to the backward traveling envelope ($F_-$) and an arbitrary fixed time function ($f$). Note the origin is fixed to a point on the forward traveling wave; that is, the origin moves with the forward traveling wave (as indicated by the dashed line).

Download Full Size | PDF

Table 1. Conversion of common master equation terms to normalized extended cavity and mean-field theory. Terms in gray are associated with the steady-state power calculation and mostly vanish. $\widetilde{f}[q](z)$ is defined according to Eq. (12) and $\langle f(z) \rangle = \frac{1}{2L_c}\int_0^{2L_c} f(z) dz$ is an average over the cavity space.

View Table | View all tables in this article

For the rest of this work, we rely on the two-level system—for reference, this leads to an evolution equation of the form

(14)$$\begin{aligned} \frac{n}{c} \frac{\partial F}{\partial T} = &\left(\frac{c}{n}\right)^{2}\left(i\frac{1}{2}k^{\prime\prime} + \frac{g_0}{2}T_2^{2}\right) \frac{\partial^{2} F}{\partial z^{2}} -\frac{g_0}{2P_s}\left(\langle K \rangle \left(\left|F\right|^{2}-P_0\right)\right. \left.+2\left(\widetilde{K}\left[\left|F\right|^{2}\right] - \langle K \rangle P_0\right) \right) F\\ &-\left.\frac{c}{n}\right.\frac{g_0}{2 P_s}\left((2T_1+3T_2) \widetilde{K} \left[F \frac{\partial F}{\partial z}^{*}\right] + (T_1+\tfrac{5}{2}T_2) \widetilde{K} \left[F^{*} \frac{\partial F}{\partial z}\right] \right)F. \end{aligned}$$

In Appendix C, the effect of additional terms are investigated, including additional third-order nonlinearities and first-order linewidth terms. These terms include ones neglected in establishing (4) and one associated with the first time derivative related to the coherence lifetime, $T_2$—For instance, the second term inside the square brackets on the first line of (5). In general, these do not have major effects. Primarily, they either modify extendon formation by a few percent or slightly adjust its group velocity.

2.2 Phase potential NLSE of a variable reflectivity cavity

In previous work [14], it was shown that a Fabry-Perot cavity with one mirror of unity reflectivity obeyed a nonlinear Schrodinger equation with a phase potential, but this result did not generalize to arbitrary cavities. From numerical simulations it is apparently the case that the result actually should generalize, but it was not proven and no expression for the important chirp bandwidth was derived for the case where $R_1$ and $R_2$ is allowed to be variable. We perform this analysis here. While the details are more complicated, the final result is similar to the unity cavity result. Consider the cross-steepening terms of the two-level system mean-field theory, which ultimately give rise to a phase potential. They take the form

(15)$$\begin{aligned} \left(\frac{\partial F}{\partial T}\right)_\mathrm{XS} & ={-}\left(\frac{c}{n}\right)^{2}\frac{g_0}{2P_s}\left( \left(2T_1+3T_2\right) \tilde{K}\left[ F\tfrac{\partial F^{{\ast}}}{\partial z}\right] +\left(T_1+\tfrac{5}{2}T_2\right) \tilde{K}\left[F^{{\ast}} \tfrac{\partial F}{\partial z}\right] \right) F. \end{aligned}$$

In order to evaluate this for a general cavity, two assumptions are made. First, the averaged field envelope ($F$) is assumed to have a constant amplitude ($A$). This is typically a very good approximation, as fast gain saturation suppresses amplitude fluctuations. Doing so allows the amplitude to factor out of the convolution and leaves only the phase ($\phi$), since $\tilde {K}\left [F\tfrac {\partial F^{\ast }}{\partial z} \right ]=-i|A|^{2}\tilde {K}\left [\tfrac {\partial \phi }{\partial z}\right ]=-\tilde {K}\left [F^{\ast } \tfrac {\partial F}{\partial z}\right ]$. It is interesting to note that even when the field has significant amplitude fluctuations, this approximation still ends up being very good; see Appendix D for a detailed discussion of why this is. Secondly, the steady-state power profile is assumed to be piecewise linear (see Fig. 3). This also turns out to be a surprisingly good approximation when compared to using a more accurate numerically simulated steady-state power profile; the key to this feature is most of the relevant comb dynamics are governed by the discontinuities caused by the facets, not what happens in between. These discontinuities will ultimately each give rise to a delta function when (15) is evaluated. The system’s power profile and its derivative can be simply expressed as follows:

(16)$$P\left(z\right) = \begin{cases} P_0+\frac{z}{L_c} \Delta_1, & \text{for } 0\le z <L_c \\ \left(1-R_2\right)\left(P_0+ \Delta_1\right)+\frac{z-L_c}{L_c} \Delta_2, & \text{for } L_c\le z<2L_c \end{cases}$$

and

(17)$$\frac{\partial P}{\partial z} = \begin{cases} \frac{\Delta_1}{L_c} -\Delta P_2\delta\left(z\right), & \text{for } 0\le z <L_c \\ \frac{\Delta_2}{L_c} -\Delta P_1\delta\left(z-L_c\right), & \text{for } L_c\le z<2L_c \end{cases}$$

Within the constant amplitude assumption, the only terms that can change drastically with a general cavity are those containing a $\widetilde {K}\left [\partial \phi /\partial z\right ]$. This is because it is the main operator that relies on the steady-state power explicitly—the steady-state power also affects $\langle K \rangle$ but this quantity is very straightforward to calculate. To evaluate $\widetilde {K}\left [\partial \phi /\partial z\right ]$, one starts with (12) to obtain

(18)$$\widetilde{K}\left[\frac{\partial\phi}{\partial z}\right] =\frac{1}{2L_c} \int_{z}^{z+2L_c}K\left({-}u\right) \frac{\partial\phi}{\partial z}\left(z-2u\right)du$$

(19)$$=\frac{1}{2L_c} \left[K\left({-}u\right)\phi\left(z-2u\right)\right]_0^{2L_c}-\frac{1}{2L_c}\int_{z}^{z+2L_c}\frac{\partial K}{\partial z}\left({-}u\right)\phi\left(z-2u\right)du.$$

Now, recalling that $K(z)=P(z)/P_0$ and that the first term goes to zero (because $K\left (z\right )$ and $\phi \left (z\right )$ are both periodic with respect to $2L_c$), one obtains the following:

\begin{aligned} \widetilde{K}\left[\frac{\partial\phi}{\partial z}\right]={} & -\frac{1}{2L_c P_0} \int_{z}^{z+2L_c}\frac{\partial P}{\partial z}\left({-}u\right)\phi\left(z-2u\right)du\\ & ={-}\frac{1}{2L_c P_0} \int_{z+2L_c}^{z}\frac{\partial P}{\partial z}\left(u\right)\phi\left(z-2u\right)du\\ & ={-}\frac{1}{2L_c P_0} \left[\int^{0}_{L_c}\left(\frac{\Delta_1}{L_c} -\Delta P_2 \delta\left(u\right)\right)\phi\left(z-2u\right)du \right.\\ & \hspace{1.5cm}\left. +\int^{L_c}_{2L_c}\left(\frac{\Delta_2}{L_c} - \Delta P_1 \delta\left(u-L_c \right)\right)\phi\left(z-2u\right)du\right]\\ & ={-}\frac{1}{2L_c P_0} \left[\left(\Delta_1+\Delta_2\right)\langle \phi \rangle -\left(\Delta P_1+\Delta P_2\right)\phi\left(z\right)\right] \end{aligned}

Now, observe that as they are defined in Fig. 3, it is the case that $\Delta _1+\Delta _2=\Delta P_1+\Delta P_2\equiv \Delta P.$ Therefore, the final result is that

(20)$$\widetilde{K}\left[\frac{\partial\phi}{\partial z}\right]=\frac{1}{2L_c} \frac{\Delta P}{P_0} \left(\phi\left(z\right)-\langle \phi \rangle\right).$$

This is the same result obtained in Ref. [14] for a cavity with a unity reflector and leads to the same NLSE, the only difference being that there are two $\Delta P$ components. The mean-field result then takes the form

(21)$$-i \frac{\partial F}{\partial T} = \frac{1}{2}\left(\beta-i\frac{D_g}{2}\right) \frac{\partial^{2} F}{\partial z^{2}} + \gamma \left|F\right|^{2} \left(\arg F-\langle \arg F \rangle\right) F +i r \left(\left|F\right|^{2}-P_0\right)F,$$

where the normalized parameters are defined in Table 2 for reference.

Fig. 3. (a) Power profile sketch showing key features like initial power at the left facet, $P_0$ , and the corresponding power drops at the right facet ($\Delta P_1$) and the left facet ($\Delta P_2$). (b) The analytical power function, Eq. (22), versus the numerical power function. Here $R_1=0.2$ and $R_2=0.45$. Observe that the values for $\Delta P_1$ are from the analytical expression Eq. (23). Note that the data in this plot has not been normalized by $K(z)$. (c) Numerical and analytical phase function calculated using (29) with the same parameters as part (b), showing good agreement. Parameter values for (b) and (c) are included in Table 4. (d) Comparison of the mean-field theory, the analytical theory with several sets of mirror reflectivities, and mid-IR experimental data. The experimental and simulation data was taken in a cavity 6 mm long, with $J=1150$ $\textrm {A}/\textrm {cm}^{2}$, dispersion $k''=-4000 \textrm {fs}^{2}/\textrm {mm}$, and $R_1=R_2=0.3$ (the mirror reflectivities of the analytical solutions were changed according to the legend).

Download Full Size | PDF

Table 2. Properties of fundamental extendons in terms of physical laser parameters for a general cavity. Note $N_\textrm {d}$ is the atomic density, the dipole matrix element is for a two level laser system, and $w_0$ is the equilibrium population inversion. For numerical values refer to Table (3) in Appendix A.

View Table | View all tables in this article

To properly understand the result from (20), one must apply boundary conditions to obtain expressions for $\Delta P_1$ and $\Delta P_2$. For these derivations, one assumes the intracavity power profile is exponential and has gain balanced by mirror losses $\alpha _m\equiv \frac {1}{2L_c}\ln \left (\frac {1}{R_1R_2}\right )$ (a valid approximation near threshold):

(22)$$P\left(z\right) = \begin{cases} P_0 e^{\alpha_m z}, & \text{for } 0\le z <L_c \\ R_2 P_0 e^{\alpha_m z}, & \text{for } L_c\le z<2L_c, \end{cases}$$

With this assumption, it is straightforward to show the expressions for $\Delta P_1$ and $\Delta P_2$ are:

(23)$$\Delta P_1=\frac{P_0}{\sqrt{R_1R_2}}\left(1-R_2\right) \hspace{0.4cm}\textrm{and}\hspace{0.4cm} \Delta P_2=P_0\left(\frac{1}{R_1} -1\right).$$

Note that if $R_2=1$, then $\Delta P_1=0$ and $\Delta P=\Delta P_2$ as before.

In spite of this similarity, the initial power required at the left facet does change depending on the setup. To derive this result, note that in steady state the gain must approximately equal the loss,

(24)$$g_0 \left(1-\frac{1}{P_s} \left(P_+{+}2P_-\right)\right)=\alpha_w+\alpha_m.$$

Averaging both sides so that $\left (P_++2P_-\ \right )=3 \langle P \rangle$, the result for $\langle P \rangle /P_s$ is

(25)$$\frac{\langle P \rangle}{P_s}=\frac{1}{3}\left(1-\frac{1}{g_0} \left(\alpha_w+\alpha_m\right)\right).$$

Alternatively, $\langle P \rangle$ can be calculated as

(26)$$\langle P \rangle={} \frac{1}{2L_c}\int_{0}^{2L_c} P\left(z\right)dz =\frac{1}{2L_c} \left[\int_{0}^{L_c} P_0 e^{\alpha_m z} dz+\int_{L_c}^{2L_c}R_1 P_0 e^{\alpha_m z} dz\right] =\frac{\Delta P}{2L_c \alpha_m}$$

Combining Eq. (25) and Eq. (26), one gets

(27)$$\Delta P=\frac{2L_c\alpha_m}{3}P_s\left(1-\frac{1}{g_0}\left(\alpha_w+\alpha_m\right)\right).$$

Finally, using (23), one can rewrite (27) in terms of $P_0$ as

(28)$$P_0=\frac{1}{3}\frac{2L_c\alpha_mP_s}{\frac{1}{\sqrt{R_1R_2}}\left(1-R_2\right)+\left(\frac{1}{R_1}-1\right)} \left(1-\frac{1}{g_0}\left(\alpha_w+\alpha_m\right)\right).$$

Hence, an expression for the initial steady-state intensity at the left facet was derived in terms of physical parameters. Due to our formalism, Eq. (28) is not symmetric with respect to $R_1$ and $R_2$ — it matters which side one takes as a starting point.

Because the form of the results are the same, in the absence of gain curvature the chirp bandwidth of an extendon can be calculated from (1) as $BW=\frac {1}{2\pi }\frac {c}{n}A_0^{2} \frac {\gamma }{\beta }2L_c$. By assuming that $A_0^{2}\approx P_0$ and plugging in (27) and (28)), one can show that the chirp bandwidth formula is unchanged from previous work, resulting in

(29)$$BW \approx \frac{1}{24\pi} 2L_c \alpha_m \left(g_0-\alpha_w-\alpha_m\right) \frac{1}{k^{\prime\prime}} \left(T_1+\frac{1}{2}T_2\right).$$

As a way of verifying the results derived in this section, we use Fig. (3d) to compare the frequency chirp from the mean-field theory, this analytical theory (with different mirror reflectivities), and experimental data from a mid-infrared QCL. The model matches the mean-field theory and experimental data quite well. The comparison of the analytical data with $R_1=R_2=0.3$ and the analytical data with $R_1=0.09$ and $R_2=1$ reveals that the bandwidth result is unchanged (the chirp is the same) when the mirror reflectivity product is the same. This result is interesting because it implies that comb bandwidth is not impacted by the particular values of $R_1$ and $R_2$, provided their product, $R_1R_2$, remains constant. One can then derive a product that maximizes the comb bandwidth as

(30)$$R_1R_2=\exp{\left[{-}L_c\left(g_o-\alpha_w\right)\right]}.$$

Based on typical simulation parameters, (30) predicts $R_1R_2 \approx 0.2$ to maximize the comb bandwidth. This result supports findings from [34] which numerically predicted that high mirror losses could be beneficial for widening the bandwidth. Referring back to (28), it turns out one can minimize $P_0$ by making $R_1<<R_2$, thereby maximizing the extracted power. Combining this result and the mirror reflectivity product result in Eq. (30), produces a general prescription for modifying $R_1$ and $R_2$ to maximize chirp bandwidth and maximize the extracted power. One can implement this design strategy either with coatings or with double-chirped mirrors, which naturally make one side high-reflectivity [4,37].

For ease of reference, all of the main results from this section are tabulated in Table 2.

2.3 Harmonic states: the higher-order eigenstates of the phase NLSE

Another consequence of (29) is it demonstrates the relevance of higher-order linearly-chirped harmonic states (which have been demonstrated experimentally in [27]). Recall that for an FM comb state to form, its comb bandwidth cannot exceed the available gain bandwidth. Now note that for the Nth order harmonic state the bandwidth defined in (29) is divided by N (this will be shown later), and hence, higher order states have a much narrower maximum bandwidth. This implies these higher order harmonic states can form in situations with much lower group velocity dispersion ($k''$) than the fundamental state, a situation which can occur due to the dispersion’s dependence on bias [38].

As harmonic states are of experimental interest, we now perform a more thorough investigation of some of their properties and link them to the higher-order eigenstates of the phase NLSE. An interesting consequence of the mean-field theory is that it can reproduce both fundamental FM combs and harmonic combs. By initializing the electric field with a field that is periodic over 1/Nth a round trip, stable linearly-chirped harmonic combs are produced rather than fundamental combs. These turn out to be numerically unstable—eventually decaying to the fundamental due to machine precision-level asymmetry—but can be resymmetrized to remain stable indefinitely. The phase profiles of the first five harmonic states are plotted in Fig. 4. Note that they agree very well with the form predicted by Eq. (3), which constructs them in terms of the fundamental state.

Fig. 4. The first five harmonic states plotted in descending order. The analytical (red) data was obtained using Eq. (3) and the numerical (blue) data was obtained by numerically solving the mean-field theory. These phases have been arbitrarily offset by their own average for ease of comparison. The numerical harmonic states use gain curvature whereas the analytical harmonic states do not.

Download Full Size | PDF

In order to better understand the properties of these states, the phase NLSE can be reformulated as an eigenvalue problem expressed in a finite basis of the harmonic extendons. By computing matrix elements and applying the well-known LCAO method [39], the carrier-envelope offset frequencies and stability properties of the various extendons can be better-understood. First, the phase NLSE is rewritten as

(31)$$i\frac{\partial F}{\partial T}=\widehat{H}F,$$

where $\widehat {H}$ is the Hamiltonian of the phase NLSE, given by

(32)$$\widehat{H}={-}\frac{1}{2}\beta\frac{\partial^{2}}{\partial z^{2}}-\gamma \left|F\right|^{2}\left(\arg{F}-\langle\arg{F}\rangle\right)-ir\left(\left|F\right|^{2}-P_0\right).$$

At this point, the phase NLSE resembles a traditional time-evolution problem, albeit a nonlinear one. It is therefore instructive to consider the eigenvalues, defined by the $F_N$ that satisfy

(33)$$\widehat{H}F_N=\lambda_N F_N.$$

Functions that satisfy this equation will make up the stationary solutions or harmonic states of the system. The full solution can then be expressed as a sum of these stationary states multiplied by an appropriate complex exponential. Combining (31) and (33), one gets

(34)$$F\left(T,z\right)=\sum_{N}^{\infty} a_N F_N\left(z\right)e^{{-}i\lambda_N T},$$

where $F_N$ is the Nth harmonic state and $a_N$ is a weighting factor. Note that for completeness one should also include the trivial $N=0$ state, the continuous-wave (CW) solution with $F=\sqrt {P_0}$.

Because harmonic states have the same chirp as fundamental states, they are also solutions to the NLSE, and can also be taken as eigenstates. However, in truth neither the fundamental extendon nor the harmonic extendon calculated above are the true eigenstates, as they are not differentiable at the jump. Currently, we do not have a working form for the true eigenstates—they may not even be expressible in terms of elementary functions. Nevertheless, because the extendons previously proposed are very close to the true eigenstates in an integrated error sense, it is valuable to consider the expectation value of the Hamlitonian as the approximate eigenvalues.

To find an approximate expression for the $\lambda _N$, the diagonal elements of $\widetilde {H}$, the matrix corresponding to $\widehat {H}$, are calculated. We note that by following this procedure one obtains the same results compared to using the full matrices, but the interpretation is more straightforward. Finding these diagonal elements is a simple matter of calculating the inner product

(35)$$\widetilde{H}_{NN}=\langle F_N \left|H\right| F_N\rangle=\int_{{-}L_c}^{L_c}F_N^{{\ast}} \left(HF_N\right)dz,$$

and normalizing to the overlap

(36)$$\widetilde{S}_{NN}=\langle F_N | F_N \rangle=\int_{{-}L_c}^{L_c} F_N^{{\ast}} F_N dz.$$

Performing this calculation, one finds that

(37)$$\lambda_{N}=\frac{\widetilde{H}_{NN}}{\widetilde{S}_{NN}}=\frac{2 L_c^{2} r^{2} \gamma^{2} P_0^{2}}{3\beta\left(2r+\gamma\right)^{2}}\left(\frac{1}{N}+i\frac{3\beta}{\gamma A_0^{2} L_c^{2}}\right).$$

The imaginary component calculated above arises from the delta function that is generated by the discontinuiuty in the frequency at the boundary. One expects that the true eigenstates have imaginary parts that are exactly zero, but even if that is not the case they are numerically very small. Thus, for an undisturbed system with no gain curvature, the eigenvalues have the form

(38)$$\operatorname{Re}\lambda_N=\frac{\lambda}{N^{2}},$$

where $\lambda$ is used to denote the eigenvalue of the fundamental and $N$ is the index of the harmonic. This result holds true for all harmonic states of the form given in Eq. (3). Referring to (34), the evolution of each eigenstate has $e^{-i\lambda _N T}$ dependence, and so $\lambda _N$ has an additional interpretation as the carrier-envelope offset frequency of the extendon. Specifically,

(39)$$f_{N}^{\left(\textrm{ceo}\right)}=\frac{1}{2\pi}\frac{d}{dT}\arg{F_N}={-}\frac{1}{2\pi}\lambda_N.$$

Hence, the real part of these eigenvalues provides the CEO angular frequency its corresponding harmonic state.

Next, consider the imaginary component. This value is the constant for all of the harmonic states and manifests as a consequence of the discontinuity in the phase equation. This discontinuity ultimately arises from the boundary conditions present at the cavity’s mirrors. The impact of this term is essentially to multiply each $F_N$ with a term which grows exponentially with $T$. From this, one would expect the harmonic states to be equally stable—but one can verify numerically that the lower-energy states are actually more robust. We will explore the implications of this analysis extensively in the subsequent sections, and show how defects can modify this behavior.

Though transition between high order harmonic states and the fundamental state has been exhibited as a function of bias [40], we have found that under the same parameter space, both states can be supported. No tuning of the dispersion is actually needed, only a change in history. Although the fundamental is slightly preferred, a small defect in the cavity can favor a harmonic instead. Additionally, the dispersion of optical systems such as QCLs can be a strong function of bias [38]. Depending on the frequency of interest, one could find that the dispersion varies massively due to the complex interaction of several intersubband features, including but not limited to the lasing transition. The dispersion can even cross zero—such crossing points would naturally favor harmonic states.

Note also that while the results reported here use mean-field theory to show the existence of harmonic states, previous explanations of harmonic states relied on the concept of parametric gain induced by temporal population pulsations [40]. In fact, population pulsations ultimately explain their appearance in the mean-field theory as well. Mathematically, pulsations manifest through the five temporally-varying nonlinear terms (the last two lines of each model in Appendix E). The details behind their derivation were discussed in [30], but briefly the Maxwell-Bloch equations were solved for population, inserted into coherence, and linearized in time derivative. Essentially, population pulsations in the Maxwell-Bloch equations become non-instantaneous $\chi ^{(3)}$ nonlinearities in the master equations. While our emphasis here is on time-domain analysis, one could alternatively perform frequency-domain analysis of the kind shown in [40] and show the existence of coherent parametric gain that favors harmonic sidebands.

2.4 Effect of gain curvature, Kerr nonlinearity, and linewidth enhancement factor on stable extendons

Although the results in the previous sections suffice for guiding intuition, they fail to account for the effects of gain curvature, a necessary consequence of the finite linewidth of a real laser. While it is known that gain curvature plays a crucial role in destabilizing combs at low dispersions, for stable combs its role is typically small (though it still has impact in modifying parameters such as dispersion). With gain curvature, the phase NLSE begins to more closely resemble the Complex-Ginzburg Landau equation, albeit one with an additional phase potential. In order to examine these effects, one can adjust many of the previous results for the phase-dominated NLSE by perturbatively modifying Eq. (1) and its related results.

To accomplish this goal, one essentially takes $\beta$ to be complex, i.e., replacing it with $\beta -\frac {i}{2}D_g$. Then, the phase NLSE (Eq. (21)) is solved using a Gaussian ansatz

(40)$$F_{t}(z,T)\equiv \sqrt{R_0 e^{{-}az^{2}}} e^{i\left(q z^{2} -\lambda_t T\right)},$$

where a is the decay factor of the Gaussian, $R_0$ is its amplitude, q is the chirp parameter, and $\lambda$ is the system’s eigenvalue, related to the average phase. From here, one can perturbatively solve for the unknown variables $q$, $a$, $\lambda$, and $R_0$. After series expansions around $D_g=0$ and $\gamma =0$ and algebra, these can be approximated as:

(41)$$\begin{aligned} q&=\frac{\gamma A_0^{2} }{2 \beta }+D_g \frac{A_0^{4} \gamma ^{3} L_c^{2} }{24 \beta ^{3} r} \end{aligned}$$

(42)$$\begin{aligned}a&=D_g\frac{A_0^{2} \gamma ^{2}}{4 \beta ^{2} r} \end{aligned}$$

(43)$$\begin{aligned}R_{0}&=A_0^{2}-D_g\frac{A_0^{4} \gamma ^{3} L_c^{2}}{24 \beta ^{2} r^{2}} \end{aligned}$$

(44)$$\begin{aligned}\lambda_t&= \frac{A_0^{4} \gamma ^{2} L_c^{2}}{6 \beta }-D_g\frac{A_0^{2} \gamma}{4 \beta } \end{aligned}$$

The details are long and are contained in a supplemental notebook, as we show in Code 1 (Ref. [41]). The leftmost term in the above expressions represents the equivalent term unmodified by gain curvature. While only approximate, these terms are a reasonable approximation for moderate values of $D_g$, as demonstrated in Fig. 5(a). The most important effect of gain curvature is to introduce a non-zero amplitude decay, which causes frequency modulation to induce amplitude modulation, decreasing intensity away from the center of the FM. It also modifies the chirp, peak power, and eigenvalue by a few percent. We note that the results of this section will eventually breakdown in slow enough gain media, this is due to the expansion about $\gamma =0$ since $\gamma$ increases in slower gain media. At that point one should set $D_g=0$ and use the derived formulas as a decent approximation.

Fig. 5. (a) Gain curvature modifications (red) for the fundamental solution’s intensity and phase compared with the numerical solution (blue). (b) Effect of a linewidth enhancement factor $\alpha =0.5$ on numerical simulation, along with the analytical calculation. Although the frequency and peak gain are shifted, the chirp is unaffected.

Download Full Size | PDF

The extension to harmonic extendons is straightforward. The eigenvalues are modified by a factor of $N^{2}$, leaving

(45)$$\begin{aligned} \operatorname{Re} \lambda_{tN} \approx {} & \frac{1}{N^{2}} \left(\frac{A_0^{4} \gamma ^{2} L_c^{2}}{6 \beta }-D_g\frac{A_0^{2} \gamma}{4 \beta } \right), \end{aligned}$$

while the chirp bandwidth $BW=q\frac {2L_c c}{\pi n}$ is reduced by N, leading to

(46)$$BW_t = \left( \frac{\gamma A_0^{2} }{2 \beta }+D_g \frac{A_0^{4} \gamma ^{3} L_c^{2} }{24 \beta ^{3} r}\right)\frac{2L_cc}{\pi n N}.$$

For ease of reference, these results are summarized in Table 3.

Table 3. Summary of key extendon parameters, including intensity, phase, CEO frequencies, harmonic states, chirp bandwidth, Kerr nonlinearity, and effective dispersion.

View Table | View all tables in this article

The extendon theory can also be modified to capture the effects of Kerr nonlinearity and linewidth enhancement. Kerr nonlinearity induces a mean-field change of

(47)$$\frac{\partial F}{\partial T} \Bigg\rvert_{\textrm{Kerr}}={-}i \gamma_K \frac{c}{n}\left(\langle K \rangle|F|^{2}+2\tilde{K}[|F|^{2}]\right) F \approx{-}i3 \gamma_K \frac{c}{n}\langle K \rangle |F|^{2} F \equiv{-}i \gamma_K' |F|^{2} F$$

where $\gamma _K\equiv \frac {3 \chi ^{(3)} \omega _0}{2nc}$ and $\gamma _K'\equiv 3 \gamma _K \frac {c}{n}\langle K \rangle$. In the absence of gain curvature, this value is constant and only shifts the CEO frequency of the extendon. As was previously shown, gain curvature causes the power to sag parabolically according to $-az^{2}$, shifting the effective dispersion according to

(48)$$\begin{aligned} \beta_\textrm{eff} &= \beta - \frac{1}{2r}D_g \gamma_K'\\ k^{\prime\prime}_\textrm{eff} &= k^{\prime\prime} - 2 P_s T_2^{2} \gamma_K. \end{aligned}$$

Linewidth enhancement factors (LEFs) have also been proffered as an important aspect of comb generation. Though two-level systems do not naturally have LEFs, they can be added phenomenologically to the master equations by making the peak gain complex, effectively using $g_0 \rightarrow g_0\left (1+i\alpha \right )$. Note that this shifts the peak gain value and center frequency—LEFs result from gain asymmetries in models with more than two states—but we have found that it does not substantially affect the phase dynamics. To see why, consider that the effect of linewidth enhancement is to couple amplitude to phase and vice-versa. Most of the terms that appear with a $g_0$ in the mean-field theory are phase terms; the only amplitude-affecting terms are the gain curvature and power restoring terms. This affects the imaginary part of the mean-field equations as

(49)$$\begin{aligned}\operatorname{Im} \frac{\partial F}{\partial T} \Bigg\rvert_{\textrm{LEF}}&= \alpha \frac{g_0}{2}\frac{c}{n} F - \alpha \frac{g_0}{2P_s}\left( \frac{c}{n} \left(\langle K \rangle|F|^{2}+2\tilde{K}[|F|^{2}]-3\langle K \rangle P_0\right) +T_2^{2} \left(\frac{c}{n}\right)^{3} \frac{\partial^{2} }{\partial z^{2}}\right)F\\ &= \alpha \frac{g_0}{2}\frac{c}{n}\left(1+3\frac{P_0}{P_s}\langle P \rangle\right)F - \alpha \frac{g_0}{2P_s}\left( \frac{c}{n} \left(\langle K \rangle|F|^{2}+2\tilde{K}[|F|^{2}]\right) +T_2^{2} \left(\frac{c}{n}\right)^{3} \frac{\partial^{2} }{\partial z^{2}}\right)F. \end{aligned}$$

Close inspection reveals that the three terms represent additional CEO frequency, additional Kerr nonlinearity, and additional dispersion, respectively. The CEO shift is fairly minor, the Kerr nonlinearity $\gamma _K$ is shifted by $\frac {g_0 \alpha }{2P_s}$, and the dispersion $k''$ is shifted by $\alpha g_0 T_2^{2}$. Therefore, the total effective dispersion is

(50)$$\begin{aligned} k^{\prime\prime}_\textrm{eff} &= k^{\prime\prime} + \alpha g_0 T_2^{2} - 2 P_s T_2^{2} \left(\gamma_K+\frac{g_0 \alpha}{2P_s}\right)\\ &= k^{\prime\prime} - 2 P_s T_2^{2} \gamma_K, \end{aligned}$$

meaning that the LEF had no effect on the chirp. This is confirmed numerically in Fig. 5(b). While a LEF can have an effect on the output, it does so by changing the effective gain of the cavity (and its output power). When this gain is properly accounted for, there is no direct effect on phase and little effect on stability.

2.5 Analysis of reflective defects

The phase NLSE provides a natural framework for analyzing perturbations to simple cavities. Here, we explore the impact of adding reflective defects, which can create asymmetrical cavities and allow one to engineer the phase profile. This type of defect engineering was recently demonstrated using ring laser frequency combs [35], where it was shown that defects can induce a laser to act in a harmonic state. It is expected that similar behavior can be induced in Fabry-Perot cavities—in particular, we will demonstrate how adding a minor reflection defect in the center of the cavity causes the system to prefer the second harmonic state over the first.

To construct the effect of a minor defect inside our cavity, one can form a scattering matrix, $\widetilde {S}$, relating the incident, reflected, and transmitted fields to each other [42]. We use this information to formulate an appropriate field term for an arbitrary collection of reflection point-defects (see Appendix B for the full derivation). The resulting master equation term is

(51)$$\left(\frac{n}{c}\frac{\partial{E_\pm}}{\partial{t}}+\frac{\partial E_\pm}{\partial z}\right)_\textrm{defect}=\sum_j \frac{1}{2} \ln \left( \frac{1+r_j}{1-r_j}\right)\delta\left(z-z_j\right) E_\mp,$$

where $r_j$ is the field reflectivity of the defects and $z_j$ are their physical locations within the cavity. Converting to the normalized and extended cavity, one finds that

(52)$$\left(\frac{n}{c}\frac{\partial{F}}{\partial{t}}+\frac{\partial F}{\partial z}\right)_\textrm{defect}=\sum_j \frac{1}{2}\ln\left(\frac{1+r_j}{1-r_j}\right)\left[\delta\left(z+z_j \right)+\delta\left(z-z_j \right)\right]\frac{K_-^{1/2}\left(z\right)}{K^{1/2}\left(z\right)}F_-\left(z\right),$$

where the defect now appears twice to account for the fact that waves see it twice a round trip. Performing the mean-field averaging and employing the spatiotemporal substitution $\partial t\rightarrow \ -n/c\ \partial z$, one find that

(53)$$\left(\frac{\partial F}{\partial T}\right)_{\textrm{defect}}=\frac{c}{2L_cn}\sum_j\frac{1}{2} \ln\left(\frac{1+r_j}{1-r_j}\right)\left[\frac{K_-^{1/2}\left({-}z_j\right)}{K^{1/2}\left({-}z_j\right)}F\left(z+2z_j \right)+\frac{K_-^{1/2}\left(z_j \right)}{K^{1/2}\left(z_j\right)}F\left(z-2z_j\right)\right],$$

since the delta functions give rise to a shift. Note that the defect derivation is only valid when $|r_j|\ll 1$ if one wishes to make use of the above results for steady-state intensity, as a sufficiently large reflector will begin to violate the mean-field assumption of small internal changes. To make it valid for large values, one would need to directly modify the effective gain upon which $K\left (z\right )$ is based. This is necessary to account for the additional reflections. The impact of small $r_j$ on the effective gain is negligible. Thus, (53) provides a useful explanation of the dynamics added by a reflection defect.

Figure 6 shows the effects of a defect term on the ground state solution of a system with a defect located at $L_c/2$ and magnitude $r_1$. Applying Eq. (53), this leads to mean-field terms of the form $F\left (z-L_c\right )$ and $F\left (z+L_c\right )$, essentially coupling fields separated by a distance of $L_c$. As this matches the peak spacing for the second harmonic of the system (Fig. 4) but not the other harmonics, this leads to the second harmonic becoming the new ground state solution as the defect magnitude increases. It is worth noting that in between the regions where the first harmonic and the second harmonic are dominant, quasi-stable solutions still exist despite the chaotic appearance of the heat plot in that region (Fig. 6(a)). However, these solutions are not time invariant; they are breather states that fluctuate periodically with time. The middle plot of Fig. 6(b) shows the phase profile of one of these states at a fixed time.

Fig. 6. (a) Top: phase heatmap showing the transition between the fundamental harmonic and the second harmonic, as a function of defect size. These values were obtained after thousands of simulated round-trips; for moderate defect sizes, the system becomes unstable. Bottom: calculated diagonal elements corresponding to the first five harmonic states. (b) Phase plots at specific defect magnitudes ($r_1$ values)

Download Full Size | PDF

To provide a quantitative explanation of this behavior, the eigenvalue approach from Section 2.3 can be employed. Adding the following defect operator to the Hamiltonian established in (32), one obtains

(54)$$H_\textrm{defect} = i \frac{c}{2L_cn}\sum_j \frac{1}{2} \ln\left(\frac{1+r_j}{1-r_j}\right)\left[\frac{K_-^{1/2}\left({-}z_j \right)}{K^{1/2}\left({-}z_j \right)}T^{2z_j}+\frac{K_-^{1/2}\left(z_j\right)}{K^{1/2}\left(z_j\right)}T^{{-}2z_j}\right],$$

where $T^{t}=e^{t \partial /\partial z}$ is the shift operator $\left (T^{t}f\left (z\right )=f\left (z+t\right )\right )$. For the example detailed in Fig. 6, (54) reduces to

(55)$$H_\textrm{defect} = i \frac{c}{2L_cn}\frac{1}{2} \ln\left(\frac{1+r_1}{1-r_1}\right)\left[\frac{K_-^{1/2}\left(-\frac{L_c}{2}\right)}{K^{1/2}\left(-\frac{L_c}{2} \right)}T^{L_c}+\frac{K_-^{1/2}\left(\frac{L_c}{2}\right)}{K^{1/2}\left(\frac{L_c}{2}\right)}T^{{-}L_c}\right].$$

Denote the positive constants attached to $T^{L_c\ }$ and $T^{-L_c\ }$ as $C_+$ and $C_-$, respectively. The defect term’s contribution to $\widetilde {H}_{NN}$ is

(56)$$\begin{aligned} \left(\widetilde{H}_{NN}\right)_{\textrm{defect}}&=\langle F_N \left|i\left(C_+ T^{L_c}+C_- T^{{-}L_c}\right)\right| F_N \rangle\\ &=iC_+ \langle F_N\left(z\right)\left|F_N\left(z+L_c\right)\rangle+iC_- \langle F_N\left(z\right)\right|F_N\left(z-L_c\right)\rangle. \end{aligned}$$

Note that for even harmonic states (denoted by $F_e$), it must be the case that $F_e\left (z\right )=F_e\left (z\pm L_c\right )$. Therefore, it must be the case that

(57)$$\left(\widetilde{H}_{ee}\right)_{\textrm{defect}}=i2L_c\left(C_+{+} C_- \right)\left|F\right|^{2},$$

where $\left |F\right |^{2}$ is the amplitude squared of the harmonic state. By contrast, for odd diagonal elements, these matrix elements will be complex rather than purely imaginary. The amplitude of all the harmonic states is the same and the overlap is much worse for all other state configurations. Hence, it follows that the even index diagonal elements of $\widetilde {H}$ will receive the dominant contribution from the defects, as shown in Fig. 6(a). In other words,

(58)$$\operatorname{Im}\left(\widetilde{H}_{ee}\right)_{\textrm{defect}}>\operatorname{Im}\left(\widetilde{H}_{oo}\right)_{\textrm{defect}},$$

where $\widetilde {H}_{oo}$ represents the eigenvalue change of odd diagonal elements. For sufficiently small reflectivity, $\frac {1}{2} \ln \left (\frac {1+r_1}{1-r_1}\right ) \approx r_1$, so an approximate shift of the eigenvalue for an even state will be given by

(59)$$\lambda_{\textrm{defect}}^{\left(e\right)}=i \frac{c}{n} r_1 \left|F\right|^{2}\left(\frac{K_-^{1/2}\left(-\frac{L_c}{2}\right)}{K^{1/2}\left(-\frac{L_c}{2}\right)}+\frac{K_-^{1/2}\left(\frac{L_c}{2}\right)}{K^{1/2}\left(\frac{L_c}{2}\right)}\right).$$

This perturbation is purely imaginary, and so it will generate an evolution of $e^{-i \lambda _{\textrm {defect}}^{\left (e\right )} T}$. Because positive imaginary part corresponds to exponential growth, the presence of defects now favors even states. Once a defect is turned on, even states grow faster than odd states, causing them to be dominant. In steady-state a perturbation will always develop at the boundary to cancel this growth out; this perturbation is energetically easiest for the lowest-order state, and so the system favors the $N=2$ state.

The eigenvalue approach utilizing the harmonic basis set is quite informative, however, it does not fully explain why the system always prefers the first harmonic as its ground state solution in the absence of defects or why the ground state solution prefers to the second harmonic over other even states. The explanation is that the current formalism fails to adequately account for the boundary spike, which lifts the degeneracy of the eigenstates. Numerical investigation has revealed that a system initialized in a higher-order state always decays to the lowest possible harmonic (for an unmodified system, this would be the fundamental solution). Adding reflection defects changes this base system to prefer the lowest even state ($N=2$). Thus, to make more detailed predictions analytically, one would also need to account for boundary spike in the eigenstate. Of course, stability information can still be obtained numerically.

More generally, an important topic of interest is stability of FM combs. While an analytical theory of the stability of linearly-chirped FM combs is of major interest (as it would determine their theoretical bandwidth limits), currently no adequate stability criterion exists. Such a theory would require an exact solution for the phase NLSE, including at the boundary point. Even numerical analysis is fraught with danger—we have found that due to the fast change of the frequency at the jump point, when solving master equations the apparent stability range can vary fairly drastically with numerical choices. First order schemes introduce numerical diffusion (gain curvature), and second-order schemes introduce numerical dispersion [43]. The use of split-step methods mitigates this error, but still does not permit the determination of the all-important scaling laws.

Finding a full solution that incorporates the jump point is very challenging. Perturbative expansions have proven to have unacceptably large error, and it is not even known whether the phase NLSE admits exact solutions that can be expressed in terms of elementary functions. (Actually, it is unlikely—even the driven NLSE requires cnoidal functions to describe its solutions in a cavity [44].) However, by drawing an analogy to the complex Ginzburg-Landau equation (CGLE, whose stability criterion is well-known to be the Benjamin-Feir lines) one can make some observations. For example, transforming the phase NLSE with gain curvature into dimensionless units (replacing $F\rightarrow F/\sqrt {P_0}$, $z\rightarrow z/\sqrt {\frac {D_g}{4rP_0}}$, and $T\rightarrow T rP_0$) yields a phase NLSE of the form

\begin{aligned} \frac{\partial F}{\partial T} &= F + \left(i \frac{2\beta}{D_g} + 1\right) \frac{\partial^{2} F}{\partial z^{2}} -\left(1 - i \frac{\gamma}{r}\left(\arg F - \langle \arg F \rangle \right)\right)|F|^{2} F\\ &\equiv F + \left(i c_D + 1\right) \frac{\partial^{2} F}{\partial z^{2}} -\left(1 + i c_{NL}\right)|F|^{2} F . \end{aligned}

This generally has the same form as a standard CGLE, with a linear and nonlinear parameter, but there are two key differences. First, there is explicit dependence of the nonlinear coefficient $c_{NL}$ on the phase. Secondly, because the boundary effects determine stability, the dimensionless cavity length ($z_b\equiv L_c/\sqrt {\frac {D_g}{4rP_0}}$) plays a critical role. Thus, at least three dimensionless parameters are needed to determine stability: the ratio of dispersion and gain curvature, the ratio of cross-steepening to gain saturation, and the ratio of amplitude decay to steady-state power. Further work is needed to determine the existence and stability of exact solutions of the phase NLSE.

3. Conclusion

We have performed a comprehensive analytical treatment of FM comb behavior. We demonstrated how one can generalize active cavity mean-field theory to handle various laser systems, including systems with defects or other perturbations. We described how one can derive a phase-driven NLSE for a cavity with different facet reflectivities and derived equations for several key parameters in terms of known physical constants, and showed how it can be modified to account for nonidealities such as gain curvature, Kerr nonlinearity, and linewidth enhancement factor. Finally, we showed how a formulation of the phase-dominated NLSE as an eigenvalue problem, where harmonic states represent the higher-order eigenstates, allows one to analyze the stability these systems in the presence of defects.

Appendix A. Simulation parameters and values

For completeness, all simulation parameters used inside this work are included in Table 4. These values were always used unless otherwise specified.

Table 4. Parameters of a two-level system used during numerical analysis.

View Table | View all tables in this article

Appendix B. Derivation of reflection defect term

In this section, we show how a small reflective point-defect can be formulated in a spatially-extended way that is suitable for mean-field theory. For master equations the most straightforward approach is to use explicit boundary conditions, but for mean-field theory these are difficult to integrate. The essential idea is to convert the boundary conditions into a gain matrix that can act on both master equations at once. Consider the situation illustrated in Fig. 7, from this setup, we can form a scattering matrix, $\widetilde {S}$, relating the fields (the four shown in the figure) to each other as follows [46]:

(60)$$\begin{pmatrix} E_-^{\left(\textrm{left}\right)}\\ E_+^{\left(\textrm{right}\right)} \end{pmatrix} =\widetilde{S} \begin{pmatrix} E_-^{\left(\textrm{right}\right)}\\ E_+^{\left(\textrm{left}\right)} \end{pmatrix}$$

where

(61)$$\widetilde{S} = \begin{pmatrix} r & it\\ it & r \end{pmatrix}.$$

Here the defect reflection is given as $r$ and the transmission is given as $t$. Without loss of generality, one can shift the reference plane of the scattering matrix by multiplying left and right [42] by

(62)$$\widetilde{P} = \begin{pmatrix} -i & 0\\ 0 & 1 \end{pmatrix},$$

resulting in

(63)$$\widetilde{S}'=\widetilde{P}\widetilde{S}\widetilde{P}= \begin{pmatrix} -r & t\\ t & r \end{pmatrix}.$$

Next, this matrix is converted from a scattering matrix into a transmission matrix, in order to represent the change that occurs going from the left fields to the right fields. This is accomplished by solving for $E\pm ^{\left (\textrm {left}\right )}$ in terms of $E\pm ^{\left (\textrm {right}\right )}$, resulting in an expression that relates the fields on one side of the defect to the fields on the other side:

(64)$$\begin{pmatrix} E_+^{\left(\textrm{right}\right)} \\ E_-^{\left(\textrm{right}\right)} \end{pmatrix} = \begin{pmatrix} \frac{1}{t} & \frac{r}{t}\\ \frac{r}{t} & \frac{1}{t} \end{pmatrix} \begin{pmatrix} E_+^{\left(\textrm{left}\right)} \\ E_-^{\left(\textrm{left}\right)} \end{pmatrix} \equiv\widetilde{T}\begin{pmatrix} E_+^{\left(\textrm{left}\right)} \\ E_-^{\left(\textrm{left}\right)} \end{pmatrix}$$

We define the gain matrix $\widetilde {G}(z)$ as the matrix that satisfies

(65)$$\frac{\partial}{\partial z} \begin{pmatrix} E_+\\ E_- \end{pmatrix} =\frac{1}{2}\widetilde{G}(z) \begin{pmatrix} E_+\\ E_- \end{pmatrix}.$$

We defined this matrix because we want to insert the effects of the defect into the master equations—which requires expressing the defect in terms of its contribution to the system’s gain. Continuing the derivation, one finds that the formal solution is

(66)$$\begin{pmatrix} E_+(z_j+\epsilon)\\ E_-(z_j+\epsilon) \end{pmatrix} =\exp{\left(\frac{1}{2}\int_{z_j-\epsilon}^{z_j+\epsilon}\widetilde{G}(z')dz'\right)} \begin{pmatrix} E_+(z_j-\epsilon)\\ E_-(z_j-\epsilon) \end{pmatrix}.$$

Equating the exponential on the right-hand side of Eq. (66) with the transmission matrix $\widetilde {T}$, one can use the matrix logarithm to find

(67)$$\log \left(\widetilde{T}\right) = \log \begin{pmatrix} \frac{1}{t} & \frac{r}{t}\\ \frac{r}{t} & \frac{1}{t} \end{pmatrix}_{t=\sqrt{1-r^{2}}} = \begin{pmatrix} 0 & \frac{1}{2}\log \left(\frac{1+r}{1-r} \right)\\ \frac{1}{2}\log \left(\frac{1+r}{1-r} \right) & 0 \end{pmatrix} = \frac{1}{2}\int_{z_j-z}^{z_j-z}\widetilde{G}(z')dz'.$$

From this expression, it is clear that to satisfy the relationship between $\widetilde {G}$ and $\log \left (\widetilde {T}\right )$, one must take the defect’s contribution to the gain as

(68)$$\widetilde{G}_\mathrm{defect} = \begin{pmatrix} 0 & \log \left(\frac{1+r}{1-r} \right) \delta \left(z-z_j\right)\\ \log \left(\frac{1+r}{1-r} \right) \delta \left(z-z_j\right) & 0 \end{pmatrix}.$$

This causes counterpropagating waves to be coupled by a factor of $\log \left (\frac {1+r}{1-r} \right ) \delta \left (z-z_j\right )$, leading to the master equation

(69)$$\left(\frac{n}{c}\frac{\partial{E_\pm}}{\partial{t}}+\frac{\partial E_\pm}{\partial z}\right)_\textrm{defect}=\sum_j \frac{1}{2} \ln \left( \frac{1+r_j}{1-r_j}\right)\delta\left(z-z_j\right) E_\mp.$$

Fig. 7. Forward and backward traveling waves encountering a defect

Download Full Size | PDF

Finally, note that one must take some caution when constructing mean-field theory. Within the extended cavity formalism, defects actually appear twice—once for each pass through the cavity.

Appendix C. Effect of additional third-order nonlinearities: nonlinear group velocity shifts

Previous analysis [30] has discussed additional self-steepening terms as being a source of group delay, but this has not been studied as they are not key to understanding FM comb formation. However, they are worthy of investigation because they are a potential explanation for a commonly encountered discrepancy between the usual free-spectral range (FSR) of a cavity and the optimal frequency for injection locking [47]. To that end, we will derive an expression estimating this group delay in an idealized situation without gain curvature. The additional nonlinear and linear terms that can induce group velocity shifts are

(70)$$\begin{aligned} \left(\frac{\partial F}{\partial T}\right)_{\textrm{gd}}={} & \frac{g_0}{2} T_2 \left(\frac{c}{n}\right)^{2} \frac{\partial F}{\partial z}-\frac{g_0}{2P_s} \left[\left(T_1+\frac{5}{2} T_2 \right)\left( \widetilde{K}\left[\left|F\right|^{2}\right]+ \langle K \rangle \left|F\right|^{2} \right) \frac{\partial F}{\partial z}\right.\\ &\left.\hspace{3cm} +\left(T_1+\frac{3}{2} T_2 \right) \langle K \rangle F^{2} \frac{\partial F^{{\ast}}}{\partial z}\right], \end{aligned}$$

where the first term is the linear shift due to the system’s finite linewidth, and the latter three are the additional third-order nonlinearities. To find the group delay caused by the terms in (70), we need to get the prefactor attached to $\partial F/\partial z$ from each of these terms and average them over the whole round trip. To perform these averages, we will assume $F\left (z\right )$ is in the fundamental state. Then, all the averages from the first line of Eq. (70) are trivial to perform. Averaging the second line is more subtle though, so we will describe the steps in detail. First, note that

(71)$$F \frac{\partial F^{{\ast}}}{\partial z}=\frac{\partial\left|F\right|^{2}}{\partial z}-F^{{\ast}} \frac{\partial F}{\partial z}.$$

Now using Eq. (71) to rewrite the second line of Eq. (70) we obtain

(72)$$\left(T_1+\frac{3}{2} T_2\right) \langle K \rangle F\left(\frac{\partial\left|F\right|^{2}}{\partial z}-F^{{\ast}} \frac{\partial F}{\partial z}\right)=\left(T_1+\frac{3}{2} T_2\right) \langle K \rangle \left(F \frac{\partial\left|F\right|^{2}}{\partial z}-\left|F\right|^{2} \frac{\partial F}{\partial z}\right).$$

In Eq. (72) we have a term with an explicit dependence on $\partial F/\partial z$ which can be treated the same as the other terms in Eq. (70). The $\left (\partial \left |F\right |^{2}\right )/\partial z$ term on the other hand is still problematic. It can be thought of as an intensity (or amplitude) induced group delay source unconnected to the phase. For systems with negligible gain curvature the amplitude is practically constant (save the boundary spike) and hence, the adding the averaged results provides a reasonably accurate description of the group delay associated with these terms (see Fig. 8). However, when stronger gain curvature is present, the amplitude is curved in a quadratic fashion and this term can no longer be ignored. Due to its atypical form, it is difficult to analyze currently.

Fig. 8. Analytically calculated group delay from self-steepening terms and refractive index shift (red) compared to numerically simulated group delay using the same terms. The numerical version is not a function of dispersion, but the analytical version is slightly impacted. The shift in the system’s free spectral range (FSR) is also calculated.

Download Full Size | PDF

In the case of negligible gain curvature as discussed above, the group delay per round trip associated with the terms specified in Eq. (70), denoted by $T_d$, can be approximated as

(73)$$T_d \approx \frac{g_0}{2P_{\textrm{sat}}}\left[T_2 P_{\textrm{sat}}-\left(T_1+\frac{7}{2} T_2\right)\left|F\right|^{2} \frac{\langle P \rangle}{P_0}\right]\left(\frac{c}{n}\right).$$

Eq. (73) gives the group delay in seconds. We also calculated the percent FSR shift caused by this change in group delay, defined as

(74)$$\Delta \mbox{FSR}=\frac{1}{T_r}-\frac{1}{T_r+T_d}$$

The agreement between the theoretical and numerical results, shown in Fig. 8, is quite good. While the numerical value of the fractional FSR shift is somewhat smaller than what has been reported (0.3% here, versus the 2.7% measured in Ref. [47]—using an injected cavity FSR of $\Delta f_{\textrm {rep}} = 7.4$ MHz and a frequency shift of 200 kHz), it is conceivable that additional mechanisms (such as the temperature dependence of the refractive index) can contribute to these shifts.

Appendix D. Effect of amplitude fluctuations on analysis

At first glance, it may appear that the constant-amplitude assumption is fairly severe, particularly given that substantial amplitude fluctuations are typically observed even in FM combs. However, the analytical results that derive from it are actually much more accurate than one may initially expect. The reason for this is subtle, but ultimately comes from the fact that the nonlinear terms that appear due to population pulsations in the mean-field theory are convolutional in nature. For example, consider the key approximation of $\tilde {K}\left [F\tfrac {\partial F^{\ast }}{\partial z} \right ]\approx -i|A|^{2}\tilde {K}\left [\tfrac {\partial \phi }{\partial z}\right ]$, which is what enables one to derive the phase NLSE from the mean-field theory. Because the $\tilde {K}$ operator acts as a low-pass filter, amplitude fluctuations are heavily dampened. This is illustrated in Fig. 9, which compares $\textrm {Im} \,\tilde {K}\left [F\tfrac {\partial F^{\ast }}{\partial z} \right ]$ and $-\langle |A|^{2}\rangle \tilde {K}\left [\tfrac {\partial \phi }{\partial z}\right ]$ for various amplitude profiles. For typical mean-field simulations of FM combs (Fig. 9(a)), the minor perturbations due to the amplitude pulsation and due to gain curvature allow the constant-amplitude approximation to remain excellent. If the amplitude is now significantly modulated, as is shown in Fig. 9(b), the result is still fairly good. Though an offset has been introduced, the overall shape remains similar (approximately parabolic).

Fig. 9. Efficacy of constant amplitude approximation. a. For typical simulated amplitudes (top), the constant amplitude approximation is excellent (bottom), and the result is essentially parabolic. b. If significant amplitude modulation is now introduced (top), the approximation (bottom) remains suprisingly adequate. This explains why FM combs can be produced even in the face of rather severe amplitude modulation.

Download Full Size | PDF

Appendix E. Full model of a two-level system

For completeness, we include here the full model of a two-level Fabry-Perot system. The master equations are

(75)$$\begin{aligned}\frac{n}{c} &\frac{\partial E_\pm}{\partial t}\pm \frac{\partial E_\pm}{\partial z} ={-}\frac{\alpha_w}{2} E_\pm{+}i\frac{k^{\prime\prime}}{2}\frac{\partial^{2} E_\pm}{\partial t^{2}} -i \gamma_K (|E_\pm|^{2}+2|E_\mp|^{2}) E_\pm\\ &+ \frac{g_0}{2}\left(1+i \alpha\right)\bigg[ 1-\frac{1}{P_s}(\left|E_\pm\right|^{2}+2\left|E_\mp\right|^{2}) -T_2 \frac{\partial}{\partial t} +T_2^{2} \frac{\partial^{2} }{\partial t^{2}}\\ &\hspace{1cm}+\frac{1}{P_s}\bigg( (2T_1+3T_2) \frac{\partial E_\mp^{{\ast}}}{\partial t} E_\mp{+}(T_1+\tfrac{5}{2}T_2) E_\mp^{{\ast}} \frac{\partial E_\mp}{\partial t}\\ &\hspace{1.7cm}+(T_1+\tfrac{5}{2}T_2) E_\mp^{{\ast}} E_\mp \frac{\partial}{\partial t} + (T_1+\tfrac{5}{2}T_2) E_\pm^{{\ast}} \frac{\partial E_\pm}{\partial t} +(T_1+\tfrac{3}{2}T_2) \frac{\partial E_\pm^{{\ast}}}{\partial t} E_\pm \bigg) \bigg] E_\pm. \end{aligned}$$

where the parameters are defined in the main text, and the mean-field theory is described by

(76)$$\begin{aligned}\frac{\partial F}{\partial T} &= i\frac{k^{\prime\prime}}{2}\left(\frac{c}{n}\right)^{3}\frac{\partial^{2} F}{\partial z^{2}} -i \gamma_K \frac{c}{n}\left(\langle K \rangle|F|^{2}+2\tilde{K}[|F|^{2}]\right) F + i \alpha \frac{g_0}{2}\frac{c}{n}F\\ &+ \frac{g_0}{2}\left(1+i \alpha\right)\bigg[ -\frac{1}{P_s} \frac{c}{n} \left(\langle K \rangle|F|^{2}+2\tilde{K}[|F|^{2}]-3\langle K \rangle P_0\right) + T_2 \left(\frac{c}{n}\right)^{2} \frac{\partial}{\partial z} +T_2^{2} \left(\frac{c}{n}\right)^{3} \frac{\partial^{2} }{\partial z^{2}}\\ &\hspace{2.1cm}-\left(\frac{c}{n}\right)^{2}\frac{1}{P_s}\bigg( (2T_1+3T_2) \tilde{K}\left[\tfrac{\partial F^{{\ast}}}{\partial z} F\right] +(T_1+\tfrac{5}{2}T_2) \tilde{K}\left[F^{{\ast}} \tfrac{\partial F}{\partial z}\right]\\ &\hspace{2.8cm}+(T_1+\tfrac{5}{2}T_2) \left(\tilde{K}\left[|F|^{2}\right]+\langle K \rangle |F|^{2} \right) \tfrac{\partial}{\partial z} +(T_1+\tfrac{3}{2}T_2) \langle K \rangle \tfrac{\partial F^{{\ast}}}{\partial z} F \bigg) \bigg] F. \end{aligned}$$

The last two lines in each case contain the cross- and self-steepening nonlinearities, which arise as a consequence of the laser’s population pulsations.

Funding

National Defense Science and Engineering Graduate; Office of Naval Research (N00014-21-1-2735); National Science Foundation (ECCS-2046772); Air Force Office of Scientific Research (FA9550-20-1-0192).

Acknowledgments

The authors acknowledge support by the Department of Defense (DoD) through the National Defense Science and Engineering Graduate (NDSEG) Fellowship Program, the Air Force Office of Scientific Research Young Investigator Program (award number FA9550-20-1-0192), the National Science Foundation CAREER program (award number ECCS-2046772), and the Office of Naval Research Young Investigator Program (award number N00014-21-1-2735).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are available online, as we show in Code 2 (Ref. [45]).

References

1. A. Yariv, “Internal Modulation in Multimode Laser Oscillators,” J. Appl. Phys. 36(2), 388–391 (1965). [CrossRef]

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

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

4. 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,” Nat. Photonics 8(6), 462–467 (2014). [CrossRef]

5. M. Singleton, P. Jouy, M. Beck, and J. Faist, “Evidence of linear chirp in mid-infrared quantum cascade lasers,” Optica 5(8), 948 (2018). [CrossRef]

6. J. Hillbrand, D. Auth, M. Piccardo, N. Opačak, E. Gornik, G. Strasser, F. Capasso, S. Breuer, and B. Schwarz, “In-Phase and Anti-Phase Synchronization in a Laser Frequency Comb,” Phys. Rev. Lett. 124(2), 023901 (2020). [CrossRef]

7. L. A. Sterczewski, C. Frez, S. Forouhar, D. Burghoff, and M. Bagheri, “Frequency-modulated diode laser frequency combs at 2 μm wavelength,” APL Photonics 5(7), 076111 (2020). [CrossRef]

8. M. W. Day, M. Dong, B. C. Smith, R. C. Owen, G. C. Kerber, T. Ma, H. G. Winful, and S. T. Cundiff, “Simple single-section diode frequency combs,” APL Photonics 5(12), 121303 (2020). [CrossRef]

9. C. Kriso, A. Barua, O. Mohiuddin, C. Möller, A. Ruiz-Perez, W. Stolz, M. Koch, and A. Rahimi-Iman, “Signatures of a frequency-modulated comb in a VECSEL,” Optica 8(4), 458–463 (2021). [CrossRef]

10. M. Heck, A. Renault, E. Bente, Y.-S. Oei, M. Smit, K. Eikema, W. Ubachs, S. Anantathanasarn, and R. Nötzel, “Passively Mode-Locked 4.6 and 10.5 GHz Quantum Dot Laser Diodes Around 1.55 μm With Large Operating Regime,” IEEE J. Select. Topics Quantum Electron. 15(3), 634–643 (2009). [CrossRef]

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

12. 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(5), 499 (2016). [CrossRef]

13. D. J. Benirschke, N. Han, and D. Burghoff, “Frequency comb ptychoscopy,” Nat. Commun. 12(1), 4244 (2021). [CrossRef]

14. D. Burghoff, “Unraveling the origin of frequency modulated combs using active cavity mean-field theory,” Optica 7(12), 1781 (2020). [CrossRef]

15. S. Coen, H. G. Randle, T. Sylvestre, and M. Erkintalo, “Modeling of octave-spanning Kerr frequency combs using a generalized mean-field Lugiato–Lefever model,” Opt. Lett. 38(1), 37–39 (2013). [CrossRef]

16. L. A. Lugiato, F. Prati, M. L. Gorodetsky, and T. J. Kippenberg, “From the Lugiato–Lefever equation to microresonator-based soliton Kerr frequency combs,” Phil. Trans. R. Soc. A. 376(2135), 20180113 (2018). [CrossRef]

17. Y. K. Chembo and C. R. Menyuk, “Spatiotemporal Lugiato-Lefever formalism for Kerr-comb generation in whispering-gallery-mode resonators,” Phys. Rev. A 87(5), 053852 (2013). [CrossRef]

18. A. B. Matsko, A. A. Savchenkov, W. Liang, V. S. Ilchenko, D. Seidel, and L. Maleki, “Mode-locked Kerr frequency combs,” Opt. Lett., OL 36(15), 2845–2847 (2011). [CrossRef]

19. A. Gordon and D. Majer, “Coherent transport in semiconductor heterostructures: A phenomenological approach,” Phys. Rev. B 80(19), 195317 (2009). [CrossRef]

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

21. C. Silvestri, L. L. Columbo, M. Brambilla, and M. Gioannini, “Coherent multi-mode dynamics in a quantum cascade laser: Amplitude- and frequency-modulated optical frequency combs,” Opt. Express 28(16), 23846 (2020). [CrossRef]

22. Y. Wang and A. Belyanin, “Harmonic frequency combs in quantum cascade lasers: Time-domain and frequency-domain theory,” Phys. Rev. A 102(1), 013519 (2020). [CrossRef]

23. P. Tzenov, D. Burghoff, Q. Hu, and C. Jirauschek, “Analysis of Operating Regimes of Terahertz Quantum Cascade Laser Frequency Combs,” IEEE Trans. Terahertz Sci. Technol. 7(4), 351–359 (2017). [CrossRef]

24. M. Dong, S. T. Cundiff, and H. G. Winful, “Physics of frequency-modulated comb generation in quantum-well diode lasers,” Phys. Rev. A 97(5), 053822 (2018). [CrossRef]

25. N. Henry, D. Burghoff, Y. Yang, Q. Hu, and J. B. Khurgin, “Pseudorandom dynamics of frequency combs in free-running quantum cascade lasers,” Opt. Eng. 57(01), 1 (2017). [CrossRef]

26. N. Henry, D. Burghoff, Q. Hu, and J. B. Khurgin, “Temporal characteristics of quantum cascade laser frequency modulated combs in long wave infrared and THz regions,” Opt. Express 26(11), 14201 (2018). [CrossRef]

27. M. Piccardo, P. Chevalier, B. Schwarz, D. Kazakov, Y. Wang, A. Belyanin, and F. Capasso, “Frequency-Modulated Combs Obey a Variational Principle,” Phys. Rev. Lett. 122(25), 253901 (2019). [CrossRef]

28. N. Opačak and B. Schwarz, “Theory of Frequency-Modulated Combs in Lasers with Spatial Hole Burning, Dispersion, and Kerr Nonlinearity,” Phys. Rev. Lett. 123(24), 243902 (2019). [CrossRef]

29. M. Piccardo, B. Schwarz, D. Kazakov, M. Beiser, N. Opačak, Y. Wang, S. Jha, J. Hillbrand, M. Tamagnone, W. T. Chen, A. Y. Zhu, L. L. Columbo, A. Belyanin, and F. Capasso, “Frequency combs induced by phase turbulence,” Nature 582(7812), 360–364 (2020). [CrossRef]

30. D. Burghoff, Supplementary document for Frequency-modulated combs as phase solitons - 4926929.pdf, Optica (2020).

31. L. Columbo, M. Piccardo, F. Prati, L. A. Lugiato, M. Brambilla, A. Gatti, C. Silvestri, M. Gioannini, N. Opačak, B. Schwarz, and F. Capasso, “Unifying Frequency Combs in Active and Passive Cavities: Temporal Solitons in Externally Driven Ring Lasers,” Phys. Rev. Lett. 126(17), 173903 (2021). [CrossRef]

32. M. Jaidl, M. Jaidl, N. Opačak, M. A. Kainz, M. A. Kainz, S. Schönhuber, S. Schönhuber, D. Theiner, D. Theiner, B. Limbacher, B. Limbacher, M. Beiser, M. Beiser, M. Giparakis, M. Giparakis, A. M. Andrews, A. M. Andrews, G. Strasser, G. Strasser, B. Schwarz, B. Schwarz, J. Darmo, J. Darmo, K. Unterrainer, and K. Unterrainer, “Comb operation in terahertz quantum cascade ring lasers,” Optica 8(6), 780–787 (2021). [CrossRef]

33. A. Forrer, Y. Wang, M. Beck, A. Belyanin, J. Faist, and G. Scalari, “Self-starting harmonic comb emission in THz quantum cascade lasers,” Appl. Phys. Lett. 118(13), 131112 (2021). [CrossRef]

34. M. Beiser, N. Opačak, J. Hillbrand, G. Strasser, and B. Schwarz, “Engineering the spectral bandwidth of quantum cascade laser frequency combs,” Opt. Lett. 46(14), 3416–3419 (2021). [CrossRef]

35. D. Kazakov, N. Opačak, M. Beiser, A. Belyanin, B. Schwarz, M. Piccardo, M. Piccardo, F. Capasso, and F. Capasso, “Defect-engineered ring laser harmonic frequency combs,” Optica 8(10), 1277–1280 (2021). [CrossRef]

36. N. Opačak, S. D. Cin, J. Hillbrand, and B. Schwarz, “Frequency comb generation by Bloch gain induced giant Kerr nonlinearity,” Phys. Rev. Lett. 127(9), 093902 (2021). [CrossRef]

37. Y. Yang, D. Burghoff, J. Reno, and Q. Hu, “Achieving comb formation over the entire lasing range of quantum cascade lasers,” Opt. Lett. 42(19), 3888–3891 (2017). [CrossRef]

38. D. Burghoff, Y. Yang, J. L. Reno, and Q. Hu, “Dispersion dynamics of quantum cascade lasers,” Optica 3(12), 1362–1365 (2016). [CrossRef]

39. R. S. Mulliken, “Electronic Population Analysis on LCAO-MO Molecular Wave Functions. IV. Bonding and Antibonding in LCAO and Valence-Bond Theories,” The Journal of Chemical Physics 23(12), 2343–2346 (1955). [CrossRef]

40. T. S. Mansuripur, C. Vernet, P. Chevalier, G. Aoust, B. Schwarz, F. Xie, C. Caneau, K. Lascola, C.-e. Zah, D. P. Caffey, T. Day, L. J. Missaggia, M. K. Connors, C. A. Wang, A. Belyanin, and F. Capasso, “Single-mode instability in standing-wave lasers: The quantum cascade laser as a self-pumped parametric oscillator,” Phys. Rev. A 94(6), 063807 (2016). [CrossRef]

41. L. Humbard and D. Burghoff, Gain curvature perturbation analysis, figshare (2021) [retrieved 20 December 2021], https://doi.org/10.6084/m9.figshare.16783651.v1.

42. D. M. Pozar, Microwave Engineering (Wiley, 2004), 3rd ed.

43. P. B. Tzenov, Modeling and simulations of quantum cascade lasers for frequency comb generation, Ph.D. thesis, Technische UniversitÄt München (2018). 00000.

44. W. H. Renninger and P. T. Rakich, “Closed-form solutions and scaling laws for kerr frequency combs,” Sci. Rep. 6(1), 24742 (2016). [CrossRef]

45. L. Humbard and D. Burghoff, MATLAB scripts and data used for plot generation, figshare (2021) [retrieved 20 December 2021], https://doi.org/10.6084/m9.figshare.17203535.v4.

46. H. A. Haus, Waves and Fields in Optoelectronics (Prentice Hall, 1983).

47. J. Hillbrand, A. M. Andrews, H. Detz, G. Strasser, and B. Schwarz, “Coherent injection locking of quantum cascade laser frequency combs,” Nat. Photonics 13(2), 101–104 (2019). [CrossRef]

Name	Description
Code 1	A Mathematica notebook that details how to modify a solution to account for gain curvature.
Code 2	Data and Matlab scripts for generating the figures in this paper.

Master equation	Normalized extended cavity	Mean-field theory
$p \hat{Q} (E_{\pm}) \hat{R} (E_{\mp})$	$f \hat{G} (F) \hat{H} (F_{-})$	$\hat{G} (F) \tilde{f} [\hat{H} (F)]$
$\frac{1}{2} (- α_{w} + g_{0}) E_{\pm} - \frac{g_{0}}{2 P_{s}} (\| E_{\pm} \|^{2} + 2 \| E_{\mp} \|^{2})) E_{\pm}$	0 $- \frac{g_{0}}{2 P_{s}} (K ({\| F \|}^{2} - P_{0})$	0 $+ 2 K_{-} ({\| F_{-} \|}^{2} - P_{0})) F$
	$- \frac{g_{0}}{2 P_{s}} (⟨ K ⟩ ({\| F \|}^{2} - P_{0})$	$+ 2 (\tilde{K} [{\| F \|}^{2}] - ⟨ K ⟩ P_{0})) F$
$E_{\pm}$	$F$	$F$
$E_{\mp}$	$\frac{K_{-}^{1 / 2}}{K^{1 / 2}} F_{-} \equiv K_{f r a c} F_{-}$	${\tilde{K}}_{f r a c} [F]$
$\| E_{\mp} \|^{2} \frac{\partial E_{\pm}}{\partial t}$	$K_{-} \| F_{-} \|^{2} \frac{\partial F}{\partial t}$	$\tilde{K} [\| F \|^{2}] \frac{\partial F}{\partial t}$
$E_{\mp} {\frac{\partial E_{\mp}}{\partial t}}^{*} E_{\pm}$	$K_{-} F_{-} {\frac{\partial F_{-}}{\partial t}}^{*} F$	$\tilde{K} [F {\frac{\partial F}{\partial t}}^{*}] F$
$\| E_{\pm} \|^{2} \frac{\partial E_{\pm}}{\partial t}$	$K \| F \|^{2} \frac{\partial F}{\partial t}$	$⟨ K ⟩ \| F \|^{2} \frac{\partial F}{\partial t}$
$E_{\pm} {\frac{\partial E_{\pm}}{\partial t}}^{*} E_{\pm}$	$K F {\frac{\partial F}{\partial t}}^{*} F$	$⟨ K ⟩ F {\frac{\partial F}{\partial t}}^{*} F$

Name	Symbol	Expression
Dipole matrix element	$\| μ \|^{2}$	$\frac{z_{0}^{2} e J T_{1}}{L_{mod} w_{0}}$
Small signal gain	$g_{0}$	$\frac{N_{d} z_{0}^{2} ω_{0} J T_{1} T_{2}}{n ϵ_{0} c ℏ L_{mod}}$
Saturation intensity	$P_{s}$	${(\frac{4 T_{1} T_{2}}{ℏ^{2}} \| μ \|^{2})}^{- 1}$
Total intensity drop	$Δ P$	$\frac{2 L_{c} α_{m}}{3} P_{s} (1 - \frac{1}{g_{0}} (α_{w} + α_{m}))$
Average intensity	$⟨ P ⟩$	$\frac{Δ P}{2 L_{c} α_{m}}$
Left facet intensity	$P_{0}$	$\frac{1}{3} \frac{2 L_{c} α_{m} P_{s}}{\frac{1}{\sqrt{R_{1} R_{2}}} (1 - R_{2}) + (\frac{1}{R_{1}} - 1)} (1 - \frac{1}{g_{0}} (α_{w} + α_{m}))$
Normalized group velocity dispersion	$β$	${(\frac{c}{n})}^{3} k^{″}$
Nonlinear cross-steepening	$γ$	$\frac{g_{0}}{2 P_{s}} {(\frac{c}{n})}^{2} (T_{1} + \frac{1}{2} T_{2}) \frac{1}{4 L_{c}} \frac{Δ P}{P_{0}}$
Amplitude relaxation	$r$	$\frac{g_{0}}{2 P_{s}} \frac{c}{n} \frac{3 ⟨ P ⟩}{P_{0}}$
Gain curvature	$D_{g}$	$2 g_{0} T_{2}^{2} {(\frac{c}{n})}^{3}$
Chirp bandwidth	$B W$	$\frac{1}{24 π} 2 L_{c} α_{m} (g_{0} - α_{w} - α_{m}) \frac{1}{k^{''}} (T_{1} + \frac{1}{2} T_{2})$
Extendon amplitude	$A_{0}$	$\sqrt{\frac{P_{0}}{1 + \frac{γ}{2 r}}}$
Fundamental eigenvalue	$λ$	$\frac{γ {\| A_{0} \|}^{2}}{2 β} (\frac{1}{3} L_{c}^{2} γ {\| A_{0} \|}^{2})$

Key extendon parameters with gain curvature, Kerr nonlinearity, and linewidth enhancement factor
Intensity	$\| F_{t} \|^{2}$	$[A_{0}^{2} - D_{g} \frac{A_{0}^{4} γ^{3} L_{c}^{2}}{24 β^{2} r^{2}}] \exp (- D_{g} \frac{A_{0}^{2} γ^{2}}{4 β^{2} r} z^{2})$
Phase	$\arg F_{t}$	$(\frac{γ A_{0}^{2}}{2 β} + D_{g} \frac{A_{0}^{4} γ^{3} L_{c}^{2}}{24 β^{3} r}) z^{2}$
Carrier-envelope offset frequency	$- \frac{λ_{t N}}{2 π}$	$- \frac{1}{2 π N^{2}} (\frac{A_{0}^{4} γ^{2} L_{c}^{2}}{6 β} - D_{g} \frac{A_{0}^{2} γ}{4 β})$
Harmonic states	$F_{t N}$	$F_{t} (z mod \frac{2 L_{c}}{N}, \frac{T}{N^{2}})$
Chirp bandwidth	$B W_{t}$	$(\frac{γ A_{0}^{2}}{2 β} + D_{g} \frac{A_{0}^{4} γ^{3} L_{c}^{2}}{24 β^{3} r}) \frac{2 L_{c} c}{π n N}$
Kerr nonlinearity	$γ_{K}$	$\frac{3 χ^{(3)} ω_{0}}{2 n c}$
Dispersion with LEF	$k_{LEF}^{″}$	$k^{″} + α g_{0} T_{2}^{2}$
Kerr nonlinearity with LEF	$γ_{K, LEF}$	$γ_{K, LEF} + α \frac{g_{0}}{2 P_{s}}$
Total effective dispersion with Kerr nonlinearity and LEF	$k_{eff}^{″}$	$k^{″} - 2 P_{s} T_{2}^{2} γ_{K}$

Name	Symbol	Value
Population lifetime	$T_{1}$	$0.4 ps$
Coherence lifetime	$T_{2}$	$50 fs$
Refractive index	$n$	$3.3$
Waveguide losses	$α_{w}$	$4 {cm}^{- 1}$
Dipole moment	$z_{0}$	$2.3 nm$
Left mirror reflectivity	$R_{1}$	variable
Right mirror reflectivity	$R_{2}$	variable
Cavity length	$L_{c}$	$4 mm$
Module length	$L_{mod}$	$58 nm$
Wavelength	$λ_{0}$	$8 μ m$
Current density	$J$	$1100 \frac{A}{{cm}^{2}}$
Dispersion	$k^{″}$	$- 2000 \frac{{fs}^{2}}{mm}$
Kerr nonlinearity	$γ_{K}$	$0$
Reflection defect	$r_{1}$	variable
Reflection defect position	$z_{1}$	$2$ mm

Analytical theory of frequency-modulated combs: generalized mean-field theory, complex cavities, and harmonic states

Abstract

1. Introduction

2. Analysis

2.1 Laser mean field theory for arbitrary sets of master equations

2.2 Phase potential NLSE of a variable reflectivity cavity

2.3 Harmonic states: the higher-order eigenstates of the phase NLSE

2.4 Effect of gain curvature, Kerr nonlinearity, and linewidth enhancement factor on stable extendons

2.5 Analysis of reflective defects

3. Conclusion

Appendix A. Simulation parameters and values

Appendix B. Derivation of reflection defect term

Appendix C. Effect of additional third-order nonlinearities: nonlinear group velocity shifts

Appendix D. Effect of amplitude fluctuations on analysis

Appendix E. Full model of a two-level system

Funding

Acknowledgments

Disclosures

Data availability

References

Supplementary Material (2)

Data availability

Cited By

Figures (9)

Tables (4)

Equations (78)

Optics Express