Efficient computation of coherent multimode instabilities in lasers using a spectral approach

Coherent multimode instabilities are responsible for several phenomena of recent interest in semiconductor lasers, such as the generation of frequency combs and ultrashort pulses. These techonologies have proven disruptive in optical telecommunications and spectroscopy applications. While the standard Maxwell-Bloch equations (MBEs) encompass such complex lasing phenomena, their integration is computationally expensive and offers limited analytical insight. In this paper, we demonstrate an efficient spectral approach to the simulation of multimode instabilities via a quantitative analysis of the instability of single-frequency lasing in ring lasers, referred to as the Lorenz-Haken (LH) instability or the RNGH instability in distinct parameter regimes. Our approach, referred to as CFTD, uses generally non-Hermitian Constant Flux modes to obtain projected Time Domain equations. CFTD provides excellent agreement with finite-difference integration of the MBEs across a wide range of parameters in regimes of non-stationary inversion, including frequency comb formation and spatiotemporal chaos. We also develop a modal linear stability analysis using CFTD to efficiently predict multimode instabilities in lasers. The combination of numerical accuracy, speedup, and semi-analytic insight across a variety of dynamical regimes make the CFTD approach ideal to analyze multimode instabilities in lasers, especially in more complex geometries or coupled laser arrays.


I. INTRODUCTION
Lasers are complex dynamical systems the operation of which is enabled by the interaction between the field of an often-multimode cavity and a gain medium.For a considerable part of the history of lasers, experimental efforts have been focused on improving stability, monochromaticity, power performance and beam quality.This is largely due to specific applications where these attributes are desirable, such as fiber optical telecommunication, materials processing, and precision spectroscopy.Even in such high-performance regimes, nonlinear, multimode effects emerge via mode competition and spatial hole burning, which must be understood and suppressed [1,2] to maximize laser power and efficiency.More recently, however, rather than an undesired effect, coherent multimode lasing phenomena have been the focus of many studies [3][4][5][6][7][8].Such phenomena include ultrashort (subps) pulse formation [9,10] useful for precision machining [11] and probing of ultrafast processes [12,13].Active mode locking is another example traditionally linked with the generation of ultrashort pulses [14,15]; recent efforts have focused on achieving mode locking in lasers with fast gain recovery, such as quantum cascade lasers [16,17].Furthermore, soliton generation has, in recent years, proven to be a disruptive technology via passive dissipative Kerr combs [18,19] as well as active platforms [20].Optical frequency combs -whether or not characterized by ultrashort pulses -are also a prime example of a technology hinging on coherent instabilities that has been impactful in applications of spectroscopy [21], metrology [22,23], optical clocks [24] and optical telecommunication [25].
In contrast to standard lasing regimes marked by stationary population inversion in the gain medium, coherent mul-timode instabilities present in the examples above are mediated by nontrivial dynamics of the lasing medium itself, referred to as population pulsations in prior work [26].One of the earliest predicted examples of such dynamics are the instabilities of single-frequency lasing in ring lasers, often classed into Lorenz-Haken (LH, or 'single-mode' ) and Risken-Nummedal-Graham-Haken (RNGH, or 'multi-mode') instabilities [27,28].Either instability emerges past a second threshold where the gain overcomes loss for additional modes, leading to steady-state multi-frequency emission.This second threshold was predicted to be at high pump powers, rendering experimental verification elusive [29][30][31][32].However, in more recent years, low-threshold multimode instabilities have been investigated in ring quantum cascade lasers [33][34][35], as well as alternate geometries including Fabry-Pérot lasers [26,36], and the instability mechanism at play has been subject of much work.
The development of simulation tools to capture such instabilities is thus very timely.While the semi-classical Maxwell-Bloch equations (MBEs) have provided the standard description of spatio-temporal lasing dynamics across the aforementioned single and multi-mode lasing regimes, the powerful simulation methods that have been developed to simulate MBEs are computationally expensive [37,38], and provide only limited analytic insight into these lasing phenomena.As a result, alternative simulation methods to more efficiently analyze lasing phenomena in nonlinear dynamical regimes are highly desirable.Such schemes are also particularly relevant with increasing interest in lasers of more complex geometries or coupled laser arrays [39][40][41][42][43], for which simulation complexity will further increase.
In this paper, we develop and validate a spectral approach to simulating complex multimode laser instabilities to address this need.While spectral or modal approaches have been used previously to analyze laser instabilities [44], they use lossless, closed cavity modes as a basis, and often consider restricted parameter regimes to reduce the dimensionality of arXiv:2306.03290v1[physics.optics]5 Jun 2023 the MBEs.Our approach differs in these two keys aspects.First, our spectral basis of choice is spanned by cold modes defined by the lossy laser cavity: the so-called constant-flux (CF) modes of the laser [45], which satisfy a generally non-Hermitian boundary value problem with associated complex eigenvalues.For cavities where the internal loss dominates over loss due to openness and radiation, such as the original RNGH and LH models we compare to , the modal theory can be cast in terms of standard closed-cavity modes with complex eigenvalues denoting intrinsic loss.Secondly, and more importantly, we analyze the full MBEs, placing no constraints on the gain medium dynamics.This requires us to use the spectral basis to perform a spatial projection of the spatiotemporal inversion field, casting the modal theory in terms of a general time-dependent inversion matrix.This non-standard representation ultimately proves crucial to correctly capturing the dynamics of the inversion field and thus of instabilities mediated by population pulsations.
The resulting coupled-mode theory, referred to as CFTD in our earlier work [46], has previously been applied to specific situations where only brief violations of stationary inversion occur, such as synchronization.However, in this paper we test and demonstrate the utility of CFTD across a wide range of regimes in ring lasers characterized by nontrivial dynamics of the gain medium, encompassing both LH and RNGH instabilities, below, past, and far above the so-called second threshold.This spans dynamical phenomena ranging from decaying spatiotemporal oscillations around a single-mode lasing state, to stable broadband frequency combs, and even chaotic dynamics.Crucially, we provide a thorough benchmarking study of CFTD across these regimes using at least one of two standard methods: finite difference time domain (FDTD) integration of multimode lasing dynamics and a split-step Runge-Kutta method (SSRK) tailored to ring lasers.We find excellent qualitative agreement across the broad range of considered regimes, with very good quantitative agreement in specific regimes that we identify.Furthermore, the agreement is obtained via CFTD simulations several orders of magnitude faster than FDTD or SSRK simulations.
Finally, the CFTD approach goes beyond providing an efficient numerical tool: it also forms the foundation for a modal linear stability analysis (LSA) that we show can be used to predict multimode instabilities.The modal LSA explicitly describes the instability of discrete modes coupled via inversion matrix elements describing gain medium dynamics.This is in contrast to more standard perturbative approaches to ring lasers, which analyze the instability of a continuous perturbation of the spatiotemporal PDE [27].Not only does the modal LSA typically agree with the spatiotemporal (ST) LSA, it can at times provide additional information.In particular, the modal LSA derived from CFTD is natively aware of the laser cavity mode spectrum, unlike the ST LSA.As a result, we find it can correctly predict the absence of instabilities in parameter regimes when a newly generated frequency does not coincide with a cavity mode, in contrast to the ST LSA.The numerical accuracy and efficiency of CFTD simulations, together with predictive capabilities provided by the modal LSA, make it ideal to study multimode instabilities in more complex, coupled laser geometries.
This paper is organized as follows.In Sec.II we recount the standard description of lasing via MBEs, including their form within the slowly-varying envelope approximation typically employed in the analysis of ring laser instabilities.In Sec.III, we present the CFTD approach starting from the general MBEs, and obtain the set of time-dependent ordinary differential equations (ODEs) that constitute the CFTD equations.Sec.IV considers the single-mode lasing regime below the multifrequency instability threshold, applied in particular to spatially non-trivial initial field distributions involving multiple spatial modes.In Sec.V, we consider dynamics above the threshold of the ring laser instability leading to stable frequency comb formation, as a function of laser loss parameters.In Sec.VI we investigate the CFTD method in a parameter space that leads to chaotic behavior, and finally in Sec.VII we quantify the simulation time improvement that our method has over spatiotemporal schemes.

II. MAXWELL-BLOCH EQUATIONS FOR RING LASERS
Lasing dynamics of a variety of lasers have been very successfully described by MBEs for the electric field inside the laser cavity E(r, t) coupled to the polarization P(r, t) and inversion D(r, t) of the gain medium, where c = 1 √ µ0ϵ0 is the speed of light in vacuum.Mathematically, Eqs.(1a)-(1c) most generally describe a set of coupled partial differential equations (PDEs) for arbitrary lasers.They must be accompanied by boundary conditions, here set by the geometry of the cavity confining the electric field.In this paper, we will consider the specific case of a ring laser (radius R and length L = 2πR), depicted schematically in Fig. 1 (a), for which the boundary conditions must be periodic.This enables a simplification of the MBEs to a set of scalar PDEs, defined along the ring laser arc length coordinate, which we label x.Then, n is the effective refractive index of the cavity medium, which can be complex to incorporate cavity loss; in particular we write it as n = n R +in I .Phenomenological parameters γ ∥ and γ ⊥ represent the population and polarization decay rate, respectively, as shown in Fig. 1.Ω is the center frequency of the gain curve, g is the dipole moment, and D 0 describes incoherent pumping threshold necessary to achieve population inversion and lasing.For simplicity, we consider a spatially homogeneous pump, although the spectral approach we employ can also account for nontrivial spatial pump profiles.
We begin by considering the standard approach to Eqs. (1a)-(1c), namely employing the slowly-varying envelope approximation [47] to reduce the order of derivatives, reducing the wave equation for electric field evolution to an advection equation (or, in nonlinear optical media, the nonlinear Schrödinger's equation).

A. Slowly-varying envelope approximation
Within the standard slowly-varying envelope approximation, the form of the MBEs is simplified by explicitly extracting the spatiotemporal dependence at the frequency set by the atomic transition frequency.For ring lasers, this takes the form: where we have also extracted dimensionful factors of the physical fields for convenience: E(x, t), P (x, t) then describe the envelopes of the total electric and polarization fields respectively, which typically evolve at frequencies much slower than the large atomic transition frequency Ω that has been explicitly extracted.Substituting Eq. (2a)-(2c) into Eqs.(1a)-(1c) and dropping terms proportional to the second-order time-derivative of the envelope fields defines the slowly-varying envelope approximation.The final MBEs under the slowly-varying envelope approximation take the form of an advection equation for the electric field (instead of a second-order wave equation), coupled to ODEs for the polarization and inversion fields: where we have also introduced the dimensionless space, time, and frequency scales: and where D 0 = D 0 Dc , while κ = nI nR Ω describes the cavity loss rate proportional to the imaginary part of the refractive index.
Having introduced the various scaling transformations via Eqs.(3), (5), we will now drop the (•) notation in the remainder of this paper, for clarity of the presentation.All quantities from here on are therefore to be understood as dimensionless, unless otherwise noted.
Eqs. (4a)-(4c) and their associated periodic boundary conditions still describe a set of PDEs, and thus must be solved using a spatio-temporal integration scheme such as a FDTD method.While well-established, such schemes are computationally expensive and scale unfavourably with system size.In the following sections, we develop an efficient spectral approach to capturing dynamics of complex multimode lasers described by the MBEs, and then benchmark this approach against more standard numerical techniques for simulating Eqs.(4a)-(4c).

III. MULTIMODE CFTD APPROACH FOR LASER DYNAMICS
Any spectral approach to analyzing laser dynamics uses a spatial basis to project Eqs.(1a)-(1c), thereby yielding a set of ordinary differential equations.The distinguishing feature of the CFTD spectral approach, introduced in Ref. [46], is the use of a general set of non-Hermitian modes as a basis: the constant-flux (CF) modes.In a multidimensional domain R, the CF modes {φ m (r)} are solutions to the generalized eigenproblem (in normalized units) ∇ 2 φ m (r) = −n 2 (r)ω 2 m φ m (r), with complex eigenfrequencies {ω m }.Crucially, the modes obey outgoing boundary conditions past the boundary of the spatial domain, ∂R, correctly accounting for losses due to the fields leaving the laser cavity.The non-Hermitian CF basis has successfully provided the foundation for Steady-state Ab Initio Lasing Theory (SALT) [45].
While the CFTD method has been developed using this most general form of the CF basis, it can be greatly simplified for the ring laser geometry we consider here.Assuming further a complex but uniform refractive index n(r) → n, the eigenproblem can be simplified to the effectively onedimensional case: where x denotes the angular variable along the ring cavity, and the modes satisfy periodic boundary conditions.The modes satisfy the orthogonality relationship: where the integral is defined over the spatial domain of the ring cavity.The modes obtained by solving the resulting Eq. ( 6) and imposing periodic boundary conditions then take the simple form of propagating plane waves: where k m = 2mπ, m ∈ Z is the wavevector for the mode indexed by integer m.The complex eigenfrequencies ω m can be found exactly, and are parameterized in terms of their real parts ν m describing the mode frequencies and imaginary parts κ m describing losses: The cold cavity mode spacing ∆ (or free spectral range (FSR)) is constant and is given by ∆ = 2π nR .It is then clear that the orthogonality relationship of Eq. ( 7) is simply that of the complex Fourier basis.We note that while we consider here the special case of modes of a ring cavity with only internal losses, the time-dependent theory we discuss holds for non-Hermitian modes defined by Eq. ( 6) for arbitrary geometries in multiple dimensions, including random lasers [48], and incorporating losses via open boundary conditions [45].
Having defined our spatial basis modes, we will now return to the lasing cavity case and expand the coupled slowlyvarying envelopes of the electric field and polarization using the following ansätze: The spatial complexity of the laser cavity including boundary conditions is entirely captured by the cavity modes {φ m (x)}, while the nontrivial time dynamics are encoded in the expansion coefficients {E m (t), P m (t)}.We then substitute Eqs.(10a), (10b) into the MBEs (Eqs.(1a)-(1c)), and make use of the orthogonality relationship to integrate out the spatial degrees of freedom over the ring cavity domain.Doing so also projects the inversion onto a set of the basis modes via: For the special case of a ring cavity, using the explicit form of the spatial basis modes, it is clear that D nm (t) are simply the spatial Fourier components at the wavevector difference k n − k m of the inversion field at a time t; however, the projection and resulting dynamical equations are more general and hold for basis modes that are not simply complex exponentials.
Leaving details of the spatial projection for Appendix D, we present the final equations of motion for the variables {E m , P m , D nm } below: where we have introduced the dimensionless mode overlap tensor A nmrs : The second line above specializes the tensor to the case of multimode ring lasers.Finally, we note that D 0 nm defines the projected matrix elements for the incoherent pump, and is defined analogously to Eq. ( 11), with D(x, t) → D 0 .This definition can equivalently be used for spatially inhomogeneous pump profiles.
Eqs. (12a)-(12c) thus define the CFTD description of multimode dynamics for ring lasers: a time-domain description derived from the underlying spatiotemporal MBEs via a specialized spatial projection, most generally using the CF basis.This leads to a close connection between CFTD and SALT, which we expand on below.

A. Dynamics across regimes of stationary and non-stationary inversion
The assumption of stationary inversion, Ḋ(t) = 0, forms the basis of successful steady-state spectral descriptions of lasing, such as SALT.However, this assumption places strong constraints on lasing phenomena described by the MBEs.To see this, note that if the inversion is stationary Eqs.(1a)-(1b) for the electric and polarization fields form a set with no nonlinear mixing of time-dependent terms.While the inversion does depend nonlinearly on E and P via Eq.(1c), it is assumed to be time-independent and hence does not lead to any frequency mixing in Eqs.(1a)-(1b).The same observation holds for the subsequently-derived CFTD equations.Stationary inversion thus precludes the generation of new frequencies via intermodulation products of the lasing field and the gain medium population dynamics (which are strongly suppressed).New frequencies can arise in this regime, but via standard multimode lasing: with increasing pump power, a new spatial mode φ m can lase if a corresponding pump threshold is crossed.These pump thresholds provide a natural truncation scheme within SALT, determining when a specific mode begins to lase and hence must be included in the spectral description.
In contrast, CFTD is not bound by the stationary inversion approximation, and is able to simulate dynamics across regions of both stationary and non-stationary inversion.For aforementioned regimes where SALT is valid, the CFTD ansatz can capture lasing modes φ m with nonzero coefficients E m (t) containing only a single frequency in the long-time limit, defining the lasing frequency for that mode m.However, by allowing E m (t) to have more a general time dependence, CFTD can also capture transient dynamics where the inversion evolves with time, to the final steady-state lasing regime where SALT operates; we discuss this dynamics in Sec.IV.We note that correspondence between CFTD and SALT in regimes where the stationary inversion approximation is briefly violated were also investigated in Ref. [46].
Most importantly, in this paper we analyze regimes where the dynamics of the inversion give rise to lasing phenomena that are qualitatively different from multimode lasing: instabilities that lead to simultaneous generation of multiple coherent frequencies from a single lasing frequency background.
The coefficients E m (t) associated with distinct spatial modes now have much more general, coupled dynamics: several coefficients can become nonzero simultaneously past a single threshold, instead of sequentially with multiple thresholds.Each E m (t) can even possess multiple distinct frequency components, a feature increasingly prevalent in more complex laser geometries [49].This more general class of dynamics is the main focus of this paper, and the subject of Secs.V and VI.

B. Numerical verification
We investigate the accuracy and efficiency of the CFTD model by comparing against standard numerical approaches to simulating Eqs.(4a)-(4c).The most natural comparison of the CFTD method, which is geometry-agnostic, is against a completely general spatiotemporal FDTD scheme.For specific comparisons, however, we also employ the split-step Runge-Kutta method (SSRK) [50,51].This scheme is tailored to ring lasers, and hence can be expected to perform optimally for the present case, but lacks generalizability to other geometries.We begin with comparisons in both simple single-mode lasing regimes but with nontrivial initial spatial field distributions (i.e.spanning multiple cavity modes) in Sec.IV.When this single-mode lasing state becomes unstable, complex multimode lasing dynamics can ensue, including the formation of frequency combs which we analyze in Sec.V or even the emergence of chaotic dynamics, discussed in Sec.VI.

IV. SINGLE-MODE LASING DYNAMICS
We will begin our analysis of lasing dynamics with the simplest operating regime: single-mode lasing.Within our spectral approach, the single-mode lasing regime is characterized by restricting the expansion coefficients in Eqs.(10a), (10b) to the single-mode forms: so that the only mode with nonzero expansion coefficients is indexed by l, and is thus the solitary mode that lases.By virtue of how the inversion matrix elements are constructed -as projections onto the modes constituting the expansion of Eqs.(10a), (10b) -the resulting inversion matrix elements, Eqs.(11), similarly reduce to: As mentioned earlier, we consider a lasing cavity that experiences a spatially uniform pump gain and uniform loss profile.This is not a restriction of CFTD, but is a simplification we make for convenience of later comparisons.Under this assumption, it is easily found that the mode with lowest threshold pump power is one that is spectrally closest to the atomic transition frequency Ω (for details see Appendix E),; we index this mode by l = 0, so by definition |ν 0 − Ω| ≪ |ν l − Ω| ∀ l ̸ = 0. We emphasize, however, that this simple result holds provided the gain and loss distribution in the laser cavity is uniform.To simplify the analysis further, we consider a cavity such that ν 0 ≃ Ω, in which case the lasing occurs at the frequency Ω.This frequency is explicitly extracted in the CFTD ansätze, Eq. ( 10a), (10b), so that the lasing mode in this frame is at zero frequency.As a result, the steady-state fields E 0 = E ss , P 0 = P ss , D 00 = D ss in the single-mode lasing regime are all stationary, and are given by: Details of the derivation of the above expression can be found in Appendix E.
In this single-mode regime, it would appear that a spectral approach retaining only the mode that eventually lases would suffice, namely restricting Eqs.(10a)-( 11) to m = n = l = 0.However, our more general ansätze allows us to quantitatively capture dynamics that require modes beyond the lasing mode, for example pulsed initial conditions or spatially nonuniform pumping schemes that excite modes other than the lasing mode, even if such modes decay away in the long-time limit.
To simulate this nontrivial regime, we explore the multimode dynamics of a ring laser pumped above the single-mode lasing threshold, with an initial intensity distribution within the laser cavity that is described by a Gaussian profile.We consider pumping the laser at a power five times above the single-mode lasing threshold, p = D 0 00 D th = 5, with decay parameters γ ⊥ = 5.0, γ ∥ = 0.5, κ = 0.1, cold cavity refractive index set to n R = 1.96, and then simulate the dynamics of the electric and polarization fields at a fixed position x = 0 of the ring cavity as a function of time, using both FDTD simulations of Eqs.(4a)-(4c), and integration of the CFTD Eqs.(12a)-(12c).The CFTD approach here takes into consideration modes m ∈ −5, . . ., 5 for a total of N = 11 modes.A time-step ∆t = 2.45 × 10 −4 in units of roundtrips is used for FDTD simulations, while CFTD employs an ODE solver with an adaptive time-step.
The results are plotted in Fig. 2, in red for the FDTD simulations and blue for the CFTD simulations.We find excellent quantitative agreement between the two approaches; the lower panel zooms in on a length of time equal to 50 roundtrips, plotting both CFTD and FDTD results to highlight the agreement.Under stationary inversion, only a single-mode solution exists (see Appendix E).However, the nontrivial time evolution indicates multimode dynamics in a transient period of ∼ 200 roundtrips after which a single mode lases in the long-time limit.The CFTD dynamics are substantially more efficient to simulate than the FDTD, requiring simulation times that are about two orders-of-magnitude shorter than the FDTD for the same timestep in regimes captured by several cavity modes (a more detailed comparison of CFTD simulation times versus the number of retained modes in more complex dynamical regimes is shown in Sec.VII).
A unique aspect of the CFTD approach is the use of a projected set of inversion matrix elements D nm (t), Eq. ( 11), to capture the dynamics of the inversion field.Formally, this projection cannot be straightforwardly inverted to extract the spatio-temporal inversion field D(x, t).For the case of a ring laser, these matrix elements have a simple physical interpretation.For a spatially-uniform inversion field, D(x, t) = D(t), all off-diagonal inversion matrix elements vanish, as is seen from Eq. ( 11) and the orthogonality relationship, Eq. ( 7).This case represents the single-mode lasing regime for ring lasers.However, when multiple spatial modes with distinct wavevectors are active, a spatial inversion grating is established in the laser cavity [26]: D nm (t) then precisely represent the timeevolution of spatial Fourier components of the inversion field at the difference of wavevectors k n − k m .For more general cavities where the projecting modes {φ m (x)} are not simply Fourier basis elements, D nm (t) for n ̸ = m still represents the amplitude of this inversion grating.
The evolution of the inversion matrix elements can be seen in Fig. 3, where D nm (t) are plotted in the 2-D plots at indicated time values 1 through 4 scaled by the single-mode lasing steady-state inversion, D ss .For the first ≈ 200 roundtrips the multimode initial pulse is travelling through and in fact building up in the ring laser, and this dynamics leads to the generation of off-diagonal inversion matrix elements due to the inversion grating.However, with increasing time, additional modes decay away and the laser returns to a single-mode lasing state; concurrently, the off-diagonal inversion matrix elements are suppressed.
We thus see that the CFTD can accurately capture spatiotemporally-complex initial conditions in the regime of single-mode lasing.With increasing pump power, it is natural to ask whether the single-mode lasing regime analyzed in this section breaks down, and if other cavity modes can be coherently excited in the long-time limit, unlike the dynamics observed in Fig. 2.This is the subject of the next section.

V. COHERENT MULTIMODE LASING DYNAMICS: RNGH INSTABILITY
It can be shown that uniformly incoherently-pumped ring lasers experience gain clamping above the single-mode lasing threshold.The inversion field attains a stationary, spatiallyhomogeneous value, and prevents additional modes from crossing the lasing threshold while the first mode is still lasing.However, this does not mean that multi-mode dynamics cannot ensue: other types of instabilities do exist, where the single lasing mode itself becomes unstable, leading to the generation of new frequencies mediated by the gain medium.Such a regime is marked by the inversion field becoming timedependent, leading to so-called population pulsations [52].
The LH and RNGH instabilities [27,28] represent prototypical examples of such a strongly-pumped, nonlinear multimode lasing regime.The salient features of these instabilities were originally predicted via a spatiotemporal linear stability analysis (ST LSA) of the MBEs [27] and then verified numerically; these features are summarized in Fig. 4 and discussed below.Above a critical (normalized) pump power that we refer to as the instability threshold p th (the so-called second threshold in early literature), the single-mode lasing solution becomes unstable (see red dashed red line in Fig. 4), giving rise to symmetric sidebands at the frequencies Ω±ω inst .At higher pump powers, a continuous range of frequencies as marked by the shaded region are predicted to be unstable; the width of this region in frequency is given by ∆ inst .
Here, we note the distinction that is often made in the literature, according to the ratio of the gain linewidth γ ⊥ to the FSR ∆.In particular, for γ ⊥ /∆ < 1, a single cavity modethe lasing mode -falls under the gain curve.The resulting instability is referred to as the single-mode or LH instability.In the opposite regime where γ ⊥ /∆ > 1, several cavity modes fall under the gain curve; consequently, the emerging instability is labelled the multi-mode or RNGH instability.While dynamics in both regimes can be distinct (as we will see), both are described by the same ST LSA, with different parameters.
In specific parameter regimes, these ST LSA results can be cast in simplified analytic forms.In particular, in the RNGH case and further considering γ ⊥ ≫ ∆, γ ∥ , κ, the minimum unstable frequency and the range of unstable frequencies close to the instability threshold p th take the forms [27], Note that both ω inst and ∆ inst generally grow with γ ⊥ .While the magnitude of this linewidth in comparison to the FSR (γ ⊥ /∆) is heuristically understood to increase the participation of cavity modes in lasing, here we see its role in increasing the participation of modes in multimode instabilities.As we will see via numerical simulations, this increase in ∆ inst via γ ⊥ /γ ∥ will lead to a growing complexity of multimode dynamics.
Having described the general multimode instability, we now discuss how the CFTD approach can also predict this instability using an efficient modal linear stability analysis (modal LSA).

A. Modal LSA: instability via non-stationary inversion
Within the CFTD approach, instability of the single-mode lasing solution arises naturally via the growth of unstable sidebands for m ̸ = 0. We analyze the dynamics of these sidebands in a linearized approximation: and retaining only terms to linear order in {δE m , δP m , δD nm }.
The multimode LSA then proceeds by first obtaining linearized dynamical equations for the sideband fluctuations and inversion matrix elements.Leaving details of the derivation for Appendix F, the resulting equations can be conveniently written in matrix form as: where the component vectors δ⃗ v m include fluctuation variables for mode m: Immediately, we note that a closed set of linearized equations can only be obtained when considering the coupled dynamics of pairs of modes {m, −m}.Here J m is a 6-by-6 dynamical matrix of the fluctuations of mode m alone, including those of the inversion matrix elements.Importantly, the dynamics due to J m alone lead to a decay of the mth mode fluctuations in the long-time limit, so that the single-mode lasing solution (m = 0) remains stable.However, the matrix C m,−m couples fluctuations of different modes, in particular the symmetric mode pair {m, −m}.As J m does not lead to instability, any multimode instabilities arise due to this coupling matrix.Crucially, as we show in Appendix F, the coupling is mediated by nonzero off-diagonal inversion matrix elements D 0m , namely ±m ̸ = 0.For a ring laser with stationary inversion, the inversion is also spatially homogeneous, and from Eq. ( 11) will thus only lead to nonzero diagonal inversion matrix elements D mm .As a result, the requirement of D 0m ̸ = 0 for nonzero m in a ring laser necessitates non-stationary inversion, accounted for within the framework of CFTD.
We now present results of the CFTD stability analysis enabled by Eqs.(20), which evaluates the stability of each mode pair {m, −m} via the eigenvalues of the dynamical matrix M m .As an example, we plot the largest real part of the eigenvalue spectrum, Re Λ m , for a selection of mode pairs in Fig. 4, as a function of pump power.We see that with increasing pump power, the eigenvalue for a specific mode pair can cross the instability threshold (its real part becomes positive).These shaded regions correspond to mode pairs experiencing gain and growing around the single-mode lasing solution, as depicted in the lower panel.However, note that the eigenvalue flow can be highly monotonic, so that a specific mode pair can become sub-threshold again with increasing pump power.At any pump power, the predicted unstable mode pairs are marked with solid lines in Fig. 4(a), where they coincide perfectly with the RNGH instability region for the parameters considered.
Importantly, pump powers also exist where the CFTD modal LSA predicts no unstable mode pairs.For such pump powers, the RNGH stability analysis predicts an unstable frequency, but this frequency does not coincide with that of a cold cavity mode (dotted vertical lines).The CFTD analysis indicates that for such pump powers, there will be no multimode instability; we will verify these predictions numerically in the following section.
The modal LSA enabled by CFTD therefore agrees well with the original RNGH analysis for the considered parameters, but is based on a more general expansion that could be applied to studying multimode instabilities in alternative laser geometries.

B. Good cavity regime and slow polarization relaxation
We are now in a position to simulate the dynamics of ring lasers in parameter and pump regimes where both the approximate ST stability analysis and its modal counterpart enabled by CFTD predict multimode instabilities.We will do so by integrating the exact CFTD equations, Eqs.(12), and compare the results against integration of the MBEs under the slowlyvarying envelope approximation, Eqs.(4a)-(4c).To begin, we focus on the "good cavity" regime where κ < γ ⊥ +γ ∥ , considering γ ⊥ = 1.0, γ ∥ = 0.5, κ = 0.1 and n R = 1.96; dynamics in the "bad cavity" limit are analyzed in Sec.VI.Note that γ ⊥ < ∆ = 2π nR , so the emergent instability would typically be characterized as of the LH type.Finally, as the carrier and polarization decay rates are similar (hence, their ratio is relatively small), this regime may be representative of so-called class A lasers, depending on the choice of κ.This class of lasers includes dye lasers, HeNe lasers and quantum cascade lasers.
We simulate the CFTD equations to obtain {E m (t), P m (t), D nm (t)} for m ∈ −5, . . ., 5 (i.e.N = 11 spatial modes).This allows us to reconstruct the electric field using Eq.(10a), namely E(x, t) = m E m (t)φ m (x), as a function of normalized pump powers p over a large range past the single mode lasing threshold, evolving for t = 1000 roundtrips at each pump power.This detailed study is greatly aided by the numerical efficiency of the CFTD equations, which provide a significant simulation speedup in comparison to spatiotemporal schemes like FDTD (see Sec. VII).We first extract the resulting normalized frequency spectrum of the electric field at x = 0, F{E(x = 0, t)/E ss | 2 }, which is shown in Fig. 5(a), with amplitudes depicted according to the listed color scale.
Until the RNGH threshold p th is reached, a single-mode lasing regime is clearly observed.Past this threshold, the predicted instability of single mode lasing dynamics is observed, giving way to a stable waveform exhibiting a frequency comb in its electric field frequency spectrum.Importantly, the dominant unstable sidebands occur at Ω ± ∆, consistent with the prediction of both LSA methods.This instability propagates to a broad frequency comb with a spacing of one FSR.A more detailed look at the observed comb is presented in Fig. 5(b) for p = 16, which also shows the comparison between CFTD and FDTD, finding excellent agreement.Fig. 5(c) shows the CFTD and FDTD simulations in the time domain, demonstrating the ability of CFTD to capture not only the steady-state frequency comb but also the transient approach to this state.
Finally, the right panel of Fig. 5 (c) shows the spatial electric field profile at t = 450 (roundtrips), reconstructed from the CFTD solution.In this particular regime, each E m (t) associated with the spatial mode φ m evolves at a single dominant frequency; hence the multiple frequency components observed in Fig. 5 (b) indicate a waveform that consists of several spatial modes of the cavity being coherently excited.This yields a localized propagating pulse in the laser cavity, sometimes referred to as a Turing roll in the spatiotemporal pattern formation literature [53].We verify these CFTD results against the FDTD and SSRK (not shown), finding excellent agreement once more.
With increasing pump power, the frequency comb suddenly collapses back to stable single mode lasing operation.Here, even though the ST LSA predicts a continuous instability band, no specific mode falls into this band.The modal linear stability approach, on the other hand, predicts discrete pairs of unstable modes, and here finds no such pairs to be unstable, consistent with the simulated dynamics.At a representative pump power p = 30, the observed frequency spectrum in Fig. 5(b) using both CFTD and FDTD shows no comb formation, and the time traces in Fig. 5(c) exhibit no oscillations, instead settling into the single mode steady state solution E ss , with a uniform spatial distribution as viewed in the right panel.
With further increase in pump power, the frequency comb emerges again, now with a 2-FSR spacing.Here, the unstable sidebands at Ω ± 2∆ fall in the instability region, and are also predicted to be unstable by the modal LSA.Again at the representative value of p = 50, CFTD and FDTD show excellent agreement in capturing the extent of the frequency spectrum of the comb in Fig. 5(b), oscillations in the time domain in Fig. 5(c), and the spatial profile in the rightmost panel, which now displays two peaks and half the propagation period, consistent with the higher spacing of peaks in the frequency spectrum.
In summary, in this regime of slower polarization relaxation, we find excellent quantitative agreement between FDTD simulations of the MBE and our CFTD approach, as well as between the ST LSA and its modal counterpart derived from the CFTD equations.We also note that in this regime, both methods agree with the SSRK method specific to ring lasers (not shown).

C. Good cavity regime and slow polarization relaxation
The ratio γ ⊥ /γ ∥ plays an important role in the features of the ST LSA as explained by Eq. (18b).Thus far we have explored regimes for which γ ⊥ /γ ∥ ≤ 2, and have found the CFTD to provide excellent quantitative agreement with exact integration using FDTD, while providing the speedup advantages afforded by a temporal scheme as opposed to a spatiotemporal one.As the ratio γ ⊥ /γ ∥ increases, both the threshold frequency at which instability occurs and the width of the instability band increase, leading to more involved comb formation dynamics at high pump powers.In these more complex regimes, we find that the CFTD still provides very good qualitative agreement with more numericallyexpensive methods.For a concrete demonstration, we consider the ratio γ ⊥ /γ ∥ = 20, with decay parameters γ ⊥ = 10.0, γ ∥ = 0.5, κ = 0.1 (for supplementary simulations with different parameters, see Appendix G).Now, γ ⊥ > ∆, so the emergent instability would generally be characterized as the RNGH instability.We note also that this case where the polarization decay rate is significantly larger than the population decay rate typically describes lasers in class B, depending on the value of κ [54].This class includes lasers of high technological value such as semiconductor lasers and solid state lasers.
As the ratio of decay parameters increases, so does the number of relevant modes that must be included in the simulation.In this case, the number of modes increases up to ∼ 21 (corresponding to Ω ± 10∆).We then once again calculate the frequency spectrum of the electric field using CFTD, with the results shown in the left panel of Fig. 6.
With increasing pump power beyond the RNGH threshold p th , a stable frequency comb emerges from the single-mode lasing state, as before.However the observed spacing is now equal to 3 FSRs, as the increased ratio of γ ⊥ /γ ∥ increases the unstable mode frequency ω inst via Eq.(B1).Importantly, the unstable mode pair is exactly that predicted by both LSA schemes.With further increase in pump power, the predicted unstable mode pair moves outwards relative to the singlemode lasing frequency.Importantly, the CFTD simulations show an increase in the comb FSR with pump power, consistent with both LSA methods.Furthermore, the observed combs can even possess multiple, non-commensurate dominant frequency components, unlike the simpler dynamics at lower values of γ ⊥ /γ ∥ .
In contrast, in these regimes our benchmark FDTD method performs rather poorly.In Fig. 6(b), we compare simulated time traces using CFTD and FDTD for 700 roundtrips for three selected pump powers p = 14, p = 25 and p = 50.While CFTD demostrates sustained oscillations as necessitated by the comb spectra observed Fig. 6(a), the FDTD simulation exhibits a suppressed amplitude and quickly settles to E ss .The difference is further highlighted by a detailed look at the frequency spectrum of the electric field at p = 14 shown in Fig. 6(c) for both CFTD and FDTD.The FDTD simulation clearly does not yield a broad frequency comb like the CFTD approach; rather the central mode and only two sidebands are the main active frequencies.The "washing out" of higher frequency components is a known issue of FDTD methods, due to the increased rate of accumulation of phase errors of higher frequency components, whose the phase evolves more rapidly.By explicitly operating in the frequency domain, the CFTD appears more robust to such effects.
Deeper in this regime of fast polarization relaxation where increasing disagreements are observed between the CFTD and FDTD, the SSRK scheme provides a more exact method that can be used to understand the discrepancies.However, we recall at the outset that the SSRK is subject to a number of limitations.First, it is tailored to ring geometries only, whereas the CFTD can accommodate arbitrary geometries.Second, simulation times for the SSRK are similar to the FDTD; hence the CFTD remains a considerably faster method.Finally, the SSRK cannot be easily extended to simulate photonic molecule or other coupled laser arrangements.With these caveats, we present comparisons of the CFTD against the SSRK in regimes of fast polarization relaxation, where the FDTD performs poorly.For γ ⊥ /γ ∥ = 20, the SSRK generally shows good agreement with the CFTD.An example of the electric field frequency spectrum in the unstable regime is shown in orange in the lower panel of Fig. 6 revealing a 3-FSR comb for p = 14; the prominent frequency peaks show good agreement with the corresponding spectrum using CFTD shown in blue in the left panel, and strongly disagrees with the FDTD in red in the right panel.
For even larger γ ⊥ /γ ∥ ratios, we expect the CFTD to require the retention of even more spatial modes for simulation, leading to diminishing returns in terms of the relative simulation time advantage in comparison to the SSRK, and also possibly introducing numerical errors.We illustrate one such regime where γ ⊥ /γ ∥ = 50 in Fig. 7, showing the electric field frequency spectrum obtained using the SSRK method in orange and the CFTD in blue.We note that the SSRK predicts more than 100 distinct frequency components to be excited.In the CFTD simulation we include 35 spatial modes, which leads to a reasonable simulation time (shorter than for the SSRK) and no obvious numerical artefacts.While CFTD is not restricted to only exhibiting as many frequency components as included spatial modes (See Ref. [46,49] and discussion in Sec.III A), here we see the CFTD spectrum is signifi- cantly restricted in comparison to the SSRK.However, while CFTD may not be ideally suited to simulating such extremely broadband dynamics in ring lasers where the SSRK exists as an alternative, the generalizability to arbitrary geometries still makes CFTD an attractive proposition for analyzing multimode instabilities in more general laser systems.

VI. CHAOTIC DYNAMICS
Not all multimode instabilities yield desired stable frequency combs.In this section, we will show that the CFTD can also efficiently capture a broader class of multimode instabilities, namely chaotic lasing dynamics.To this end, we now consider a different parameter regime to our prior simulations: we maintain γ ∥ = 0.5, but set γ ∥ /γ ⊥ = 1, and importantly now choose κ = 1.1, which is just large enough for the laser system to operate in the "bad cavity" regime, where κ > γ ⊥ + γ ∥ .Note that we again have γ ⊥ < ∆, the regime of the LH instability, which in conjunction with the bad cavity limit is known to feature chaotic emission at high pump powers [44,55].We now demonstrate that this is indeed the case, and analyze CFTD performance in such regimes.
In Fig. 8 we show the spectrum of the electric field as a function of pump power using CFTD.For pump powers past p th up to around p ∼ 30, the single-mode lasing regime persists.This is consistent with the modal LSA, which predicts no unstable mode pairs in this pump range, and the ST LSA as no cold mode frequencies fall in the predicted instability region.For higher powers up to p ∼ 55, a stable frequency comb solution emerges, initially with an FSR predicted by both LSAs, but exhibiting complex frequency pulling at higher pump powers.
Such dynamics is similar to that observed in simulations in the good cavity regime in Sec.V.However, beyond p ∼ 55, the spectrum of the electric field suddenly broadens into several finely spaced, but still coherent, peaks.At even higher pump powers around p ∼ 70, however, any sharp peaks decohere into an elevated and extremely broad noise floor, features typical of chaotic regimes.This qualitative change in lasing dynamics coincides with regions where the instability bands of the ST LSA, which are symmetric with respect to the central frequency, start to overlap.Such an overlap requires a system operating in the bad cavity regime (see Appendix B), as we consider here.At a specific pump power p = 80 in this regime, a more detailed look at the frequency spectrum of the electric field is shown in Fig. 8(c), using both CFTD and FDTD schemes.Both show clear signatures of chaos: a broadband, noisy spectrum with several active frequencies and no significant coherent peaks.Time traces of the electric field using both CFTD and FDTD are also shown in Fig. 8(b).For p = 80, both methods are in good agreement, showing noisy traces with no recurrent oscillations.We also see that for specific pump powers in the other two aforemention regimes, CFTD and FDTD also show very good agreement.
We therefore clearly see that the CFTD is able to capture chaotic dynamics, as well as transitions from stable comb formation regimes to chaotic regimes.This makes the CFTD approach promising to identify -and avoid -regions of parameter space that do exhibit multimode instabilities, but do not allow useful, stable frequency comb formation.FIG. 9. Simulation time comparison between the CFTD method for various numbers of modes included in the simulation (circles) and the FDTD method (triangle).The CFTD is at least two orders of magnitude faster, even when a relatively large number of modes is included in the simulation.

VII. SIMULATION TIME COMPARISON
Having explored the numerical fidelity of the CFTD in comparison to more standard simulation schemes, we now quantify the speedup in simulation time provided by the CFTD.In Fig. 9, we show the simulation time for the CFTD as a function of the number of included modes N , and the simulation time for an equivalent FDTD integration for comparison.The results shown are for the laser system defined by parameters γ ⊥ = 5.0, γ ∥ = 0.5, κ = 0.1, n R = 1.96, and a fixed pump power p = 25; however they are representative of the speedup typically observed using the CFTD.We note that to ensure a fair comparison, the number of roundtrips simulated and the spatial discretization used is kept the same between both methods.
We consider various mode numbers from N = 7 to N = 19 for the CFTD, even though N = 15 is sufficient for the CFTD to accurate results in this particular example.The simulation time for the CFTD increases with the number of modes, but even for the highest mode numbers it remains about two orders of magnitude faster than the FDTD.This highlights one of the main advantages of our spectral CFTD approach: its computational speed, attained by the ability to reduce the integration of multidimensional PDEs to a set of ODEs in time only.The stark difference in the one-dimensional ring laser case explored here will only become larger for simulations of larger laser systems in two-or three-dimensions.

VIII. CONCLUSIONS AND OUTLOOK
In this paper, we have demonstrated a spectral approach, CFTD, as a viable and efficient method to simulate one of the most complex time-dependent phenomena in lasers: the emergence and resultant dynamics of coherent multimode instabilities.Our approach projects the standard MBEs of laser dynamics onto a suitably chosen spatial basis accounting for the lasing cavity geometry and losses, obtaining a set of timedependent coefficients of the electric field, polarization, as well as inversion which is captured via a projected set of matrix elements.Expanding far beyond our previous work [46], we have considered regimes where non-stationary inversion plays a crucial role in the emergence of coherent multimode instabilities.By benchmarking the CFTD method with an FDTD approach, as well as an SSRK scheme specific to ring lasers, we have found excellent qualitative agreement across a very wide parameter space, and quantitative agreement in regimes with up to ≈ 30 relevant modes.We reveal the ability to capture various waveforms and complex phenomena in our analysis, starting from single mode behavior to frequency combs of different FSR spacings to broadband chaotic spectra.
Our method provides not only an increase in simulation speed of at least two orders of magnitude over standard spatiotemporal integration schemes, but also more analytic insight into the system under study, via a linear stability analysis based on discrete included modes.While we have analyzed here a single ring laser, the CFTD method we introduce is a fast and efficient tool that can easily be extended to study complex time-dependent phenomena in much more general laser systems, a research area of significant contemporary interest.Such systems include complex cavity geometries without specific spatial symmetries, as well as coupled laser arrays.The complexity of analyzing such systems with conventional numerical schemes is substantial even in standard multimode lasing regimes (i.e. with stationary inversion), let alone coherent instabilities with nontrivial population dynamics.This only highlights the utility of the CFTD method.
The above Helmholtz equation can be solved by assuming ansätze for the mode functions φ m (x) of the form: where k m is a dimensionless wavevector, related to the dimensionful wavevector k m via k m = Lk m , and where N is a normalization constant.Note that for Eq.(C3) to satisfy the Helmholtz equation for the ring cavity, it must satisfy periodic boundary conditions, so that: Substituting the above and simplifying, we obtain the relation: The above immediately yields the dispersion relation: Recalling that n is complex-valued, the eigenfrequencies ω m must also be complex, and describe the decay of field intensity for the lossy ring cavity modes (in the absence of pumping and subsequent lasing).
Recalling that we have chosen parameters such that ν 0 = Ω, we find: Finally, imposing the orthogonality relationship in Eq. ( 7), we find that N = 1 √ L , so that the final expression for lossy ring cavity modes takes the form: with complex eigenfrequencies given by Eq. (C7).

Appendix D: Derivation of CFTD equations
We will derive here Eq. ( 12) starting with the Maxwell-Bloch equations for the scalar electric field amplitude E(x, t), polarization P (x, t) and inversion density D(x, t), Eqs.(1a)-(1c).The procedure for obtaining the CFTD equations can be summarized as follows: we substitute the CF ansätze, Eqs.(10a), (10b) into each of the Maxwell-Bloch equations, integrate out the spatial component by taking advantage of properties of the CF states, and apply the slowly varying envelope approximation to eliminate second-order time derivatives.
Beginning by substituting Eqs.(10a), (10b) into Eq.(1a), we obtain: Using Eq. (C1) the above equation simplifies to: Multiplying through by φ n (x) and integrating over the spatial domain of the cavity, we find: This spatial projection enables us to make use of the orthogonality of the cavity modes {φ m (x)}, Eq. ( 7), which collapses the sum over modes, finally yielding: We now also perform a slowly-varying envelope approximation to neglect second-order time derivatives as before, namely Ëm ≪ Ω Ėm and Pm ≪ Ω Ṗm ≪ Ω 2 P m .Also making use of the explicit form of scaling factors E c , P c and rearranging terms, we arrive at: We will now introduce some replacements and a final set of approximations.Using , and scaling time and frequency scales according to Eq. ( 5), we arrive at the dynamical equation for the electric field amplitude of the nth mode is given by: Which is Eq.(12a) from the main text.Similarly substituting Eqs.(10a), (10b), into the dynamical equation for the polarization field, Eq. (1b), we obtain: Multiplying through by φ n (x) and integrating over the spatial domain of the cavity, we find once more: where we have also used the explicit forms of E c , P c , and D c .The term in brackets on the right defines the inversion matrix elements introduced in Eq. (11).Introducing these matrix elements, the equation of motion for the polarization field expansion coefficients is given by: The equations of motion for the inversion matrix elements D nm are given by: Finally, substituting Eqs.(10a), (10b), into the dynamical equation for the inversion field, Eq. (1b), we obtain: Multiplying through by the scaling factor D c , and re-ordering the sum by a redefinition of labels, we obtain: We can project the inversion field onto the spatial basis.We multiply through by φ n (x)φ * m (x) and integrate over the spatial domain of the laser cavity, finally arriving at: where we have used the explicit forms of the scaling factors, moved to dimensionless time and frequency variables, and introduced the dimensionless mode overlap tensor of the main text, Eq. ( 13): Appendix E: Derivation of the single mode threshold In this appendix section, we derive the pump threshold and steady-state lasing solution in the regime of single-mode lasing.We allow for a general lasing frequency Ω l : where the inversion is stationary.Under this ansätz, the modal equations then simplify to: For steady-state, we require Ėss = Ṗss = Ḋss = 0.In this regime, Eq. (E2a) reduces to: Similarly, Eq. (E2b) yields for the same quantity: Comparing the two equations, we immediately find for D ss : Or, rationalizing the left-hand side: Comparing both sides, the imaginary parts immediately yield: while the real parts yield: Substituting the expression for D ss into the above, we immediately find: which finally yields: or: Finally, using Eq.(E2c) we find: which simplifies to: Using Eq. (E7), we can solve for E ss : which can be rewritten in the form: The requirement of |E ss | 2 > 0 yields: which yields the single-mode lasing threshold D th : Clearly minimizing this threshold for fixed damping parameters requires minimizing κ l and Ω l .Assuming all modes have equal damping rates κ l = κ ∀ l, the lasing threshold is minimized for mode l that minimizes Ω l .From Eq. (E11), Ω l = 0 if Ω 2 − ν 2 l + κ 2 l = 0, which requires: as κ l ≪ Ω. Namely, the mode l closest to the atomic transition frequency will have the lowest single-mode lasing threshold.Under this linearization, the equation of motion for the electric field amplitudes, which is already linear, retains its form and is given by Eq. (12a): The coupling matrix has a sparse form, and clearly emphasizes that sideband coupling arises via the inversion field sector of the dynamics (the final two rows).
solutions reach steady state, and the steady-state oscillation amplitude.
The intermediate pump power result, p = 55, describes a special case: here, two mode pairs are simultaneously predicted to be unstable by the ST LSA (Ω ± 4∆ and Ω ± 5∆ fall simultaneously in the instability band).As seen in Fig. 5, for simpler instabilities with a single unstable seed mode pair, this mode number determines the FSR spacing of the emergent multimode comb.With two adjacent unstable mode pairs, the nature of the comb waveform is not a priori obvious, as the mode numbers are incompatible with a fixed FSR comb.The time domain electric field traces shown in Fig. 10 using both simulation methods reveal that the dynamics is more complex than the single mode pair case.

The role of the retained number of modes
As a modal approach, one of the main parameters in the CFTD simulations is the number of modes.Generally, a larger number of modes leads to greater accuracy as defined by the discrepancy with the FDTD solution; however, if the number of modes is too large, it will lead to a greater simulation time without any benefit in accuracy and may even give rise to arti-facts.In our analysis in the main text, the choice of number of modes has been guided by two factors.First, in order to produce an initial pulse that closely resembles a Gaussian pulse, at least 11 modes have been included in each simulation regardless of the relevant number of modes in the system.For some of the regimes we have considered, this is a sufficient number of modes.The second factor we have consulted is the LSA by including at least enough modes to cover the instability bands.To highlight the importance of the number of modes, we have shown in Fig. 11 a comparison of two simulations where the difference lies only in the total number of modes included.The top panel shows the CFTD and FDTD time traces as well as spectra for a greater number of modes, in this case 15.We do see a discrepancy between the two methods in the time domain, where the FDTD is slower to reach the steady state oscillations.As we have analyzed in the main text, this parameter space belongs to a relatively complex regime, where the agreement between the two methods is more qualitative.However, the agreement over the modes present in each simulation is very good.In the bottom panel we show the same comparison between the two methods but now the number of modes included in the CFTD is only 9. The results are now qualitatively very different between the CFTD and FDTD where the CFTD settles to steady state and thus shows fewer modes to be present.In fact, the few side-bands that are present in the CFTD spectrum are due to the transient dynamics in the beginning of the time trace.

FIG. 1 .
FIG. 1.(a) Schematic of multimode ring laser under incoherent pumping.(b) Longitudinal mode structure of the ring cavity and gain medium response curve (grey).Here N = 20 cold cavity modes are depicted.

FIG. 2 .
FIG. 2. Multimode dynamics below the instability threshold (the pump power normalized by the single mode lasing threshold is p = 5).Top left panel shows the initial Gaussian electric field profile, E(x, t = 0) = 0.1e −100x 2 , scaled by Ess.Projection of this initial profile on the CFTD spatial basis is shown in the top right plot.Center panel: comparison of dynamics of the electric field intensity |E(x = 0, t)| 2 and polarization |P (x = 0, t)| 2 using CFTD (blue) and FDTD (red).Lower panel shows |E(x = 0, t)| 2 scaled by Ess for the first few roundtrips to highlight the agreement.Loss parameters are γ ⊥ = 5.0, γ ∥ = 0.5, κ = 0.1, while nR = 1.96.The number of modes included is N = 11, the minimum for reproducing an accurate Gaussian profile as the initial condition.

FIG. 3 .
FIG. 3. Top panel: Evolution of the electric field intensity in the single-mode lasing regime for 500 roundtrips.Bottom panel: Dynamics of inversion matrix elements Dnm(t) in CFTD, evaluated at labelled time points, below the single-mode instability threshold (p = 5) for an initial Gaussian electric field profile corresponding to Fig. 2. Parameters are the same as that figure.Off-diagonal elements of the matrix are suppressed as the solution reaches the single-mode steady state.

FIG. 4 .
FIG. 4. (a) ST LSA results based on the original formulation of Risken et al. [27], and corresponding modal LSA results.Beyond a specific pump threshold p th [See Eq. (17)], the single-mode lasing solution becomes unstable, at a frequency given by ωinst.Dashed lines indicate the results of the modal LSA based on linearization of the CFTD equations.(b) Real part of the eigenvalues Λm of the dynamical matrix in Eq. (20) as a function of pump power.A positive real part indicates an instability.Specific ring laser parameters considered here are γ ⊥ = 1.0, γ ∥ = 0.5, κ = 0.1, and nR = 1.96.

FIG. 5 .
FIG. 5. Emergence of frequency combs via the LH instability in ring lasers for loss parameters γ ⊥ = 1.0, γ ∥ = 0.5, κ = 0.1, and nR = 1.96.(a) Frequency spectrum of the electric field as a function of pump power obtained using CFTD.Shaded regions mark the continuous frequencies predicted to be unstable by the ST LSA, while yellow boxes mark the discrete mode pairs predicted to be unstable by the modal LSA.Cold cavity mode frequencies are shown in dashed black, while the horizontal dashed line is the instability pump threshold p th .(b) Normalized electric field spectra at specific pump powers, and (c) their corresponding time traces |E(x = 0, t)| 2 scaled by Ess, computed using both CFTD (blue) and FDTD (red) for a direct comparison.Right panel: Spatial distribution of the steady-state pulsing electric field in the ring cavity at a fixed time, |E(x, t = 450)/Ess| 2 (marked by (1) in (c)), computed using CFTD across the three pump powers simulated in the previous panels.

FIG. 6 .
FIG. 6. RNGH instability and comb formation for decay parameters: γ ⊥ = 10.0, γ ∥ = 0.5, κ = 0.1, and nR = 1.96.(a) Frequency spectrum of electric field as a function of pump power obtained using CFTD.Shaded regions mark the continuous frequencies predicted to be unstable by the ST LSA, while yellow boxes mark the discrete mode pairs predicted to be unstable by the modal LSA.The horizontal dashed line marks the RNGH instability threshold.(b) Side by side comparison of |E(x = 0, t)| 2 scaled by Ess for the first 700 roundtrips of the CFTD (blue) and FDTD (red) results.The number of modes included in the CFTD simulations is N = 19.(c) Frequency spectrum of the electric field for p = 14 for CFTD (blue) and SSRK (orange) in the left panel, and FDTD (red) in the right panel, calculated for 1100 roundtrips.Note the qualitatively good agreement between the CFTD and the SSRK, and the significant deviations of FDTD.

FIG. 7 .
FIG. 7. SSRK and CFTD comparison for decay parameters: γ ⊥ = 25.0,γ ∥ = 0.5, κ = 0.1, and nR = 1.96.The selected pump power is p = 30.Shown on the left is the frequency spectrum of the electric field as computed Fourier Transform from 1100 roundtrips of the SSRK result showing more than 100 modes.Shown on the right is the frequency spectrum of the electric field as obtained from 1000 roundtrips of the CFTD with N = 35 modes included in the simulation.

FIG. 8 .
FIG. 8. Emergence of chaotic behavior for decay parameters γ ⊥ = 0.5, γ ∥ = 0.5, κ = 1.1, and nR = 1.96.(a) Frequency spectrum of the electric field as a function of pump power obtained using CFTD with N = 17 retained modes.Shaded regions mark the continuous frequencies predicted to be unstable by the ST LSA, while yellow boxes mark the discrete mode pairs predicted to be unstable by the modal LSA.The horizontal dashed line indicates the instability threshold.(b) Side by side comparison of |E(x = 0, t)| 2 scaled by Ess for the first 2000 roundtrips of the CFTD (blue) and FDTD (red) results.(c) Fourier Transform of |E(x = 0, t)| 2 for p = 80 where we find chaotic behaviour in both the CFTD (blue) and FDTD (red) simulation.

Appendix F :
Details of the multimode LSA To analyze the linear stability of single-mode lasing dynamics described by the CFTD equations, we substitute the ansätze of Eqs.(19c) into Eqs.(12a)-(12c), and retain only terms linear in the fluctuation variables (δE m , δP m , δD nm ).
F1)while the equation of motion for polarization field, Eq. (12b), for sidebands m ̸ = 0 take the form:δ Ṗm = −γ ⊥ δP m − iγ ⊥   E 0 δD m0 + l̸ =0 δE l D 00 δ ml   (F2)where we have dropped terms involving the product of two fluctuation terms ∼ O(δE m δD ml ), m ̸ = 0 in the second line.It is clear that the nontrivial inversion matrix elements that couple to the emerging sidebands are off-diagonal elements namely D m0 .Using Eq. (12c), we can write equations of motion for these inversion matrix elements:δ Ḋm0 = −γ ∥ (δD m0 − D 0 m0 ) + iγ ∥ 2 [E 0 δP * m + δE −m P * 0 − δE * m P 0 − E * 0 δP −m ] (F3)where we have now dropped second-order terms such as O(δE r δP m+r ), r ̸ = 0, −m.Eqs.(F1)-(F3) can be compactly written in matrix form, yielding the matrix system in Eq. (20) of the main text.The component matrices are given by: mode index m arises only via the electric field sector.The coupling matrices C m,−m are symmetric under m → −m are given by: