In this Letter, we present a rigorous method to study the stability of periodic lasing systems. In a linear model, the presence of a continuum of modes (with arbitrarily close lasing thresholds) gives the impression that stable single-mode lasing cannot be maintained in the limit of an infinite system. However, we show that nonlinear effects of the Maxwell–Bloch equations can lead to stable systems near threshold given a simple stability condition on the sign of the laser detuning compared to the band curvature. We examine band edge (1D) and bound-in-continuum (2D) lasing modes and validate our stability results against time-domain simulations.

Many lasers rely on resonances in periodic systems, ranging from band edge modes of grated distributed-feedback (DFB) waveguides^{1,2} or photonic-crystal surface-emitting lasers (PCSELs)^{3–9} to more exotic bound-in-continuum (BiC) states.^{10,11} In this Letter, we address a fundamental question for periodic lasers: does stable single-mode lasing exist in an infinite periodic structure or does it inherently require the boundaries of a finite structure to stabilize? A number of theoretical works have studied lasing with periodic boundary conditions as in Fig. 1(left) and found lasing modes,^{12–17} but neglected a key concern: even if the structure and the lasing mode are periodic, stable lasing requires that arbitrary aperiodic electromagnetic perturbations [as in Fig. 1(right)] must decay rather than grow.^{18–20} At first glance, such stability may seem unlikely: any resonance in a periodic system is part of a continuum of resonances at different Bloch wavevectors with arbitrarily close lasing thresholds, and this seems to violate typical assumptions for stable lasing.^{21–23} A finite-size structure discretizes the resonance spectrum and hence may suppress this problem, but instabilities have been observed in large enough finite periodic lasers where the resonances become very closely spaced.^{24} Analogous transverse instabilities are known to occur in translation-invariant ( $ period \u2192 0 $) lasers such as vertical-cavity surface-emitting lasers (VCSELs),^{25} for which stability analysis has been performed with various assumptions.^{26,27} In fact, however, we show that single-mode lasing is possible even in infinite periodic structures for a range of powers above threshold, by applying a Bloch adaptation of linear-stability analysis to the full Maxwell–Bloch equations.^{19,20} (Instabilities can still arise if our criteria are violated or from effects such as disorder not considered in this work.) We consider examples for both 1D DFB-like lasers and 2D BiC-based lasing,^{10,11,28} and validate our result against brute-force time-domain simulations.^{29,30} Using perturbation theory (in the supplementary material), we also obtain a simple condition for stability near threshold of low-loss resonances and confirm it numerically: the sign of the laser detuning from the gain frequency should match the sign of the band curvature at threshold.

We consider lasing systems described by the semi-classical Maxwell–Bloch equations (with the rotating-wave approximation), which fully include nonlinear mode-competition effects (such as spatial hole-burning),^{31}

where $ E + $ is the positive-frequency component of the electric field (the physical field being given by $ 2 Re [ E + ] $), $ P + $ is the positive-frequency polarization describing the transition between two atomic energy levels (with frequency *ω _{a}* and linewidth $ \gamma \u22a5 $),

*D*is the population inversion (with relaxation rate $ \gamma \u2225 $),

*D*

_{0}is the pump strength,

*ϵ*is the cold-cavity real permittivity, and

_{c}*σ*is a cold-cavity conductivity loss. Here, we are assuming that the orientation of the atomic transition is parallel to the electric field and have written all three fields in their natural units.

_{c}^{19}

A steady-state solution of these equations can be obtained via steady-state *ab initio* lasing theory (SALT), which is exact for single-mode lasing and approximate for multi-mode lasing with well-separated modes.^{21–23} For a periodic system, we consider a Bloch-mode steady-state solution $ E + = E k e i ( k \xb7 x \u2212 \omega t ) $ satisfying the stationary ( $ D \u0307 = 0 $) SALT equation,

where $ \Gamma ( \omega ) = \gamma \u22a5 / ( \omega \u2212 \omega a + i \gamma \u22a5 ) , \u2009 P k = \Gamma ( \omega k ) D k E k , \u2009 D k = D 0 / ( 1 + | \Gamma ( \omega k ) E k | 2 ) $, and $ \Theta k = e \u2212 i k \xb7 x \u2207 \xd7 \u2207 \xd7 e i k \xb7 x $ is a periodic operator.

Given this steady-state solution, one can then apply linear-stability analysis to the full Maxwell–Bloch equations, linearizing arbitrary aperiodic perturbations $ X = X k + \delta X $, for $ X \u2208 { E , P , D} $, to determine whether perturbations *δX* exponentially grow (unstable) or shrink (stable).^{18–20} Here, our key point is that, because the linearized equations for the perturbations *δX* are periodic (for a Bloch-mode steady state), we can apply Bloch's theorem^{32} to decompose the perturbations themselves into Bloch-wave modes $ \delta E q $, solving a separate linear-stability eigenproblem for each wavevector **q**.

The well-known linear-stability analysis^{19} of the Maxwell–Bloch equations (1) proceeds as follows. Linearization of (1) in *δX* gives

where $ d \omega = ( d d t \u2212 i \omega ) $. Splitting complex variables into real and imaginary parts yields a set of linear equations $ ( C d 2 d t 2 + B d d t + A ) u ( x , t ) = 0 $,^{19} where $ u = ( Re ( \delta E ) , Im ( \delta E ) , Re ( \delta P ) , Im ( \delta P ) , \delta D ) $, and *A*, *B*, and *C* are operator matrices readily obtained from (3). Stability analysis consists of looking for solutions of the form $ u = Re ( U e \sigma t ) $, which leads to a quadratic eigenproblem,

The sign of $ Re ( \sigma ) $ determines the stability of the single-mode solution.^{19}

Since the operators *A*, *B*, and *C* are periodic in our case, however, we can use Bloch's theorem to further simplify the problem: the eigenfunctions can be chosen in the Bloch form $ U = U q e i q \xb7 x $, where $ U q $ is periodic. The eigenvalues $ \sigma ( q , D 0 ) $ then determine the stability: if there exists a wavevector **q** so that $ Re ( \sigma ( q , D 0 ) ) > 0 $, then the single-mode solution is unstable at the pump rate *D*_{0}, with exponential growth at the wavevector $ k \xb1 q $. Since (*A*, *B*, *C*) are real, we also have $ \sigma ( q , D 0 ) = \sigma ( \u2212 q , D 0 ) * $, so we need only consider one side of **q** within the Brillouin zone.

We can now use this method to study a simplified model for a DFB laser formed by a 1D photonic crystal with alternating layers of equal thickness and dielectric constants equal to 1 and 3 (Fig. 2). We assume a uniform conductivity loss $ \sigma c = 0.001 \omega a $ and a two-level gain medium with $ \omega a a / 2 \pi c = 0.31 $ and $ \gamma \u22a5 a / 2 \pi c = 0.008 $. Figure 2 shows part of the band diagram, with *ω _{a}* chosen near the first band edge. For every wavevector

*k*of the first band, we compute the pump threshold

*D*, defined as the lowest pump rate

_{t}*D*

_{0}that compensates the loss and leads to a real eigenfrequency

*ω*in (2). As expected, the smallest

_{k}*D*is obtained at the band edge $ k = \pi / a $ of the first band, which we, therefore, take to be the first lasing mode. However, as discussed earlier,

_{t}*D*varies continuously with

_{t}*k*and other modes are expected to reach threshold for arbitrary close values of the pump in the linear model.

In order to study the stability of the lasing band edge mode, we first solve the steady-state nonlinear equation (2) at higher pump values with a Newton–Raphson solver as described in Ref. 23. We then use the obtained steady-state solution to solve the stability eigenproblem (4) for different pump values. The results are summarized in Fig. 3. First, note that the single mode solution is stable close to threshold, unlike a linear model (Fig. 2). This can be attributed to the nonlinear gain saturation, which prevents arbitrary close modes from reaching threshold. In general, the stability of the laser depends on the relationship between the decay rates of the three fields, $ \gamma \u22a5 $ for **P**, $ \gamma \u2225 $ for *D*, and *κ* for **E**, the decay rate of the cavity in the absence of gain.^{33} When two (or more) of these decay rates become similar, we notice a sharp reduction of *D*_{0} for the onset of instability (in this case, $ \gamma \u2225 \u223c \kappa $).

Stability can also be studied using a multimode SALT by including the first lasing mode in the gain saturation and computing the pump threshold for a second lasing mode as a function of *k* [inset of Fig. 3(a)]. In particular, this coincides with the results from the stability eigenproblem in the limit $ \gamma \u2225 \u2192 0 $. Solving (3) for $ \gamma \u2225 \u2192 0 $ is indeed equivalent to having $ \delta D \u2192 0 $ and *δX* being a solution to SALT equation. As can be seen in the inset of Fig. 3(a), the nonlinear gain saturation pushes the threshold of the arbitrary close modes ( $ q \u2192 0 $) to a higher pump value compared to what is expected from a linear model. However, this multimode SALT predicts a second lasing mode that is arbitrary close to the first lasing mode, which is outside the domain of validity of SALT. Furthermore, the instability onset depends rather strongly on $ \gamma \u2225 $, emphasizing the need for a full Maxwell–Bloch stability analysis.

In order to check the stability of the lasing mode close to threshold for a general system, we use perturbation theory to compute $ \sigma ( q , D 0 ) $ near $ ( 0 , D t ) $. Analytical details are shown in the supplementary material, using methods similar to those developed in Ref. 20. In the case of small loss, we obtain a simple approximate condition for stability near threshold: the band curvature $ Re ( d 2 \omega d k 2 ) $ and the laser detuning ( $ \omega t \u2212 \omega a $) should have the same sign at threshold. When lasing at the band edge, this is equivalent to requiring *ω _{a}* to lie inside the bandgap.

We now validate the results of stability analysis against finite-difference time-domain (FDTD) simulations^{29,30} with a finite supercell and periodic boundary conditions. We initialize the simulation fields with the SALT solution plus additional noise and analyze whether the system remains in the same steady-state at later times. Note that for a supercell with $ N cells $ periods, only a finite set of values for *q* is allowed ( $ = 2 \pi \u2113 / a N cells $ for $ \u2113 = 0 , \u2026 , N cells \u2212 1 $). Figure 3(b) shows a perfect match between the two computations. In particular, the instability onset for the FDTD simulations corresponds to the value of the pump *D*_{0} for which at least one allowed *q* reaches the instability region obtained from the stability eigenproblem (4). Once instability is reached, a second lasing mode starts. This second lasing mode corresponds to the first *q* that hits the instability region. However, the new lasing solution is not accurately described by two-mode SALT [Fig. 3(c)] because the small frequency difference violates the SALT assumptions (exact in the limit $ \gamma \u2225 \u2192 0 $). In particular, the inset of Fig. 3(a) shows that the threshold of the multimode SALT (for $ q = 4 \pi / 50 a $) does not match the actual threshold for the stability eigenproblem. As $ N cells $ increases, the second lasing frequency becomes arbitrary close to the first mode, requiring an ever-smaller $ \gamma \u2225 $ for the multimode SALT approach to be viable. On the other hand, for a fixed $ N cells $, the multimode SALT approach becomes increasingly accurate for smaller $ \gamma \u2225 $. The two-mode regime here also exhibits a chaotic behavior, typical in certain classes of lasers.^{33}

We next consider a 2D (*E _{z}*-polarized) example to study the stability of a BiC lasing mode. The structure is a periodic line of surface rods placed at a distance

*L*from a perfect-metal boundary (Fig. 4, inset), which is known to have multiple BiCs.

^{28}BiCs are characterized by a quality factor $ Q \u2192 \u221e $ in the absence of an external pump and absorption loss, as seen in the inset. As in the previous 1D example, we compute the pump threshold

*D*at different wavevectors

_{t}*k*and find the lasing mode corresponding to the smallest

*D*. In this example, the first lasing mode corresponds to the BiC at $ k a = 0.4 \pi $, with $ D t \u2248 7 \xd7 10 \u2212 3 $ and a lasing frequency $ \omega t a / 2 \pi c \u2248 0.65 $. The results of the stability analysis are shown in Fig. 5(a) for $ \gamma \u2225 a / 2 \pi c = 5 \xd7 10 \u2212 3 $. We first note that the lasing mode is stable near threshold and that instability occurs at a higher pump value

_{t}*D*

_{0}[Fig. 5(b, left)]. This matches our condition for stability near threshold (positive band curvature and laser detuning). As clear from the corresponding

*q*and eigenfrequencies, instabilities at higher pump correspond to modes that become active at $ k a = 0.8 \pi $ (BiC) and $ k a = \pi $ (guided mode). A comparison between our stability results and FDTD simulations is shown in Fig. 5(a, inset), where we plot the Fourier transform of the electric field at a given point outside a rod for different pump values. The number and frequencies of lasing modes match our stability computations. Finally, in order to confirm our simple stability condition, we study the same system with a larger

*ω*corresponding to a negative laser detuning. As shown in Fig. 5(b, right), the lasing system is indeed not stable for any value of pump above threshold. Such instabilities may arise in very large systems (small

_{a}*q*).

The method presented in this Letter gives a rigorous answer to the fundamental question of stable lasing in infinite periodic systems and provides practical guidance in the form of theoretical criterion for stability. If these criteria are satisfied, the main theoretical challenges for future work are to analyze the effects of boundaries (which we expect are negligible for sufficiently large systems) and manufacturing disorder (which must eventually limit single-mode lasing).

See the supplementary material for analytical details of perturbation theory.

This work was supported in part by the U.S. Army Research Office through the Institute for Soldier Nanotechnologies under Award No. W911NF-18-2-0048.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.