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) waveguides1,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 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.

FIG. 1.

We study the stability of a single Bloch-periodic lasing mode under aperiodic perturbations. The stability eigenproblem can be solved using Bloch theorem by writing perturbations as a general Bloch wave. The lasing mode is stable when real parts of the eigenvalues σ ( q ) are negative for all wavevectors q.

FIG. 1.

We study the stability of a single Bloch-periodic lasing mode under aperiodic perturbations. The stability eigenproblem can be solved using Bloch theorem by writing perturbations as a general Bloch wave. The lasing mode is stable when real parts of the eigenvalues σ ( q ) are negative for all wavevectors q.

Close modal

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 

(1)

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 γ ), D is the population inversion (with relaxation rate γ ), D0 is the pump strength, ϵc 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.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 · x ω t ) satisfying the stationary ( D ̇ = 0 ) SALT equation,

(2)

where Γ ( ω ) = γ / ( ω ω a + i γ ) , P k = Γ ( ω k ) D k E k , D k = D 0 / ( 1 + | Γ ( ω k ) E k | 2 ) , and Θ k = e i k · x × × e i k · 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 + δ X , for X { 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 theorem32 to decompose the perturbations themselves into Bloch-wave modes δ E q , solving a separate linear-stability eigenproblem for each wavevector q.

The well-known linear-stability analysis19 of the Maxwell–Bloch equations (1) proceeds as follows. Linearization of (1) in δX gives

(3)

where d ω = ( d d t i ω ) . 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 ( δ E ) , Im ( δ E ) , Re ( δ P ) , Im ( δ P ) , δ 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 σ t ) , which leads to a quadratic eigenproblem,

(4)

The sign of Re ( σ ) 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 · x , where U q is periodic. The eigenvalues σ ( q , D 0 ) then determine the stability: if there exists a wavevector q so that Re ( σ ( q , D 0 ) ) > 0 , then the single-mode solution is unstable at the pump rate D0, with exponential growth at the wavevector k ± q . Since (A, B, C) are real, we also have σ ( q , D 0 ) = σ ( 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 σ c = 0.001 ω a and a two-level gain medium with ω a a / 2 π c = 0.31 and γ a / 2 π 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 Dt, defined as the lowest pump rate D0 that compensates the loss and leads to a real eigenfrequency ωk in (2). As expected, the smallest Dt is obtained at the band edge k = π / a of the first band, which we, therefore, take to be the first lasing mode. However, as discussed earlier, Dt varies continuously with k and other modes are expected to reach threshold for arbitrary close values of the pump in the linear model.

FIG. 2.

The cold cavity is a 1D photonic crystal with uniform conductivity loss σ c = 0.001 ω a . The two-level gain medium is characterized by ω a a / 2 π c = 0.31 and γ a / 2 π c = 0.008 . The frequency (dots) and pump (dashed lines) at the lasing threshold are computed for modes of the first band. The minimum pump at threshold is obtained at the band edge k a = π . In the absence of gain, the decay rate for the band edge mode is equal to κ 5.8 × 10 5 ( 2 π c / a ) .

FIG. 2.

The cold cavity is a 1D photonic crystal with uniform conductivity loss σ c = 0.001 ω a . The two-level gain medium is characterized by ω a a / 2 π c = 0.31 and γ a / 2 π c = 0.008 . The frequency (dots) and pump (dashed lines) at the lasing threshold are computed for modes of the first band. The minimum pump at threshold is obtained at the band edge k a = π . In the absence of gain, the decay rate for the band edge mode is equal to κ 5.8 × 10 5 ( 2 π c / a ) .

Close modal

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, γ for P, γ 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 D0 for the onset of instability (in this case, γ κ ).

FIG. 3.

(a) Stability region obtained from Maxwell–Bloch stability eigenproblem as a function of γ and pump strength D0. The inset shows the pump threshold of the second lasing mode using multimode SALT (assuming that one first mode at k a = π is lasing). This represents the limit γ 0 of the stability eigenproblem. (b) Detailed stability map for γ a / 2 π c = 10 4 as a function of q. We compare results to FDTD simulations using a finite supercell with periodic boundary conditions (unstable in shaded regions), initialized with the SALT solution plus 1 % noise and checking stability after 10 5 optical periods. Stars show the allowed q due to the finite supercell ( 2 π / a N cells ). (c) Modal intensity of lasing modes with FDTD ( N cells = 50 ) and multimode SALT (assuming second lasing mode at q = 4 π / 50 a ).

FIG. 3.

(a) Stability region obtained from Maxwell–Bloch stability eigenproblem as a function of γ and pump strength D0. The inset shows the pump threshold of the second lasing mode using multimode SALT (assuming that one first mode at k a = π is lasing). This represents the limit γ 0 of the stability eigenproblem. (b) Detailed stability map for γ a / 2 π c = 10 4 as a function of q. We compare results to FDTD simulations using a finite supercell with periodic boundary conditions (unstable in shaded regions), initialized with the SALT solution plus 1 % noise and checking stability after 10 5 optical periods. Stars show the allowed q due to the finite supercell ( 2 π / a N cells ). (c) Modal intensity of lasing modes with FDTD ( N cells = 50 ) and multimode SALT (assuming second lasing mode at q = 4 π / 50 a ).

Close modal

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 γ 0 . Solving (3) for γ 0 is indeed equivalent to having δ D 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 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 γ , 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 σ ( 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 ω d k 2 ) and the laser detuning ( ω t ω 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) simulations29,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 π / a N cells for = 0 , , N cells 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 D0 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 γ 0 ). In particular, the inset of Fig. 3(a) shows that the threshold of the multimode SALT (for q = 4 π / 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 γ 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 γ . The two-mode regime here also exhibits a chaotic behavior, typical in certain classes of lasers.33 

We next consider a 2D (Ez-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 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 Dt at different wavevectors k and find the lasing mode corresponding to the smallest Dt. In this example, the first lasing mode corresponds to the BiC at k a = 0.4 π , with D t 7 × 10 3 and a lasing frequency ω t a / 2 π c 0.65 . The results of the stability analysis are shown in Fig. 5(a) for γ a / 2 π c = 5 × 10 3 . We first note that the lasing mode is stable near threshold and that instability occurs at a higher pump value D0 [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 π (BiC) and k a = π (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 ωa 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 q).

FIG. 4.

The inset shows a 2D array of cylindrical rods with diameter = 0.7 a , ϵ c = 2.58 , σ c = 0.001 ω a , and a separation L = 1.078 a to a perfect mirror. Gain inside the rods is characterized by ω a a / 2 π c = 0.625 and γ a / 2 π c = 0.01 . Three BiCs are shown at k a = 0 , 0.4 π , 0.8 π . The minimum pump at threshold Dt is obtained at k a = 0.4 π , which is the first lasing mode. In the absence of gain, the decay rate for this mode is equal to κ 8 × 10 5 ( 2 π c / a ) . The top inset shows a positive band curvature at threshold.

FIG. 4.

The inset shows a 2D array of cylindrical rods with diameter = 0.7 a , ϵ c = 2.58 , σ c = 0.001 ω a , and a separation L = 1.078 a to a perfect mirror. Gain inside the rods is characterized by ω a a / 2 π c = 0.625 and γ a / 2 π c = 0.01 . Three BiCs are shown at k a = 0 , 0.4 π , 0.8 π . The minimum pump at threshold Dt is obtained at k a = 0.4 π , which is the first lasing mode. In the absence of gain, the decay rate for this mode is equal to κ 8 × 10 5 ( 2 π c / a ) . The top inset shows a positive band curvature at threshold.

Close modal
FIG. 5.

(a) Result from stability eigenproblem. The shaded region indicates instability. The inset shows FDTD results using a supercell with 20 unit cells and periodic boundary conditions. Plots show the Fourier transform of the electric field at a point near a rod. Small insets show the eigenvectors obtained from (4) along with their frequencies ω a / 2 π c . They do match modes obtained in the linear regime (below threshold) at k a = 0.8 π and k a = π . (b) Re ( σ ) as a function of q and D 0 / D t for different transition frequencies ω a a / 2 π c ( = 0.625 , 0.675 ) . The threshold lasing frequency ω t a / 2 π c is maintained at 0.65 . The system is unstable near threshold when the laser detuning ( ω t ω a ) has the opposite sign to the band curvature. The black solid line corresponds to Re ( σ ) = 0 .

FIG. 5.

(a) Result from stability eigenproblem. The shaded region indicates instability. The inset shows FDTD results using a supercell with 20 unit cells and periodic boundary conditions. Plots show the Fourier transform of the electric field at a point near a rod. Small insets show the eigenvectors obtained from (4) along with their frequencies ω a / 2 π c . They do match modes obtained in the linear regime (below threshold) at k a = 0.8 π and k a = π . (b) Re ( σ ) as a function of q and D 0 / D t for different transition frequencies ω a a / 2 π c ( = 0.625 , 0.675 ) . The threshold lasing frequency ω t a / 2 π c is maintained at 0.65 . The system is unstable near threshold when the laser detuning ( ω t ω a ) has the opposite sign to the band curvature. The black solid line corresponds to Re ( σ ) = 0 .

Close modal

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.

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

1.
H.
Kogelnik
and
C.
Shank
, “
Stimulated emission in a periodic structure
,”
Appl. Phys. Lett.
18
,
152
154
(
1971
).
2.
J. E.
Carroll
,
J.
Whiteaway
,
D.
Plumb
, and
R.
Plumb
,
Distributed Feedback Semiconductor Lasers
(
IET
,
1998
), Vol.
10
.
3.
M.
Imada
,
S.
Noda
,
A.
Chutinan
,
T.
Tokuda
,
M.
Murata
, and
G.
Sasaki
, “
Coherent two-dimensional lasing action in surface-emitting laser with triangular-lattice photonic crystal structure
,”
Appl. Phys. Lett.
75
,
316
318
(
1999
).
4.
M.
Meier
,
A.
Mekis
,
A.
Dodabalapur
,
A.
Timko
,
R.
Slusher
,
J.
Joannopoulos
, and
O.
Nalamasu
, “
Laser action from two-dimensional distributed feedback in photonic crystals
,”
Appl. Phys. Lett.
74
,
7
9
(
1999
).
5.
S.
Noda
,
M.
Yokoyama
,
M.
Imada
,
A.
Chutinan
, and
M.
Mochizuki
, “
Polarization mode control of two-dimensional photonic crystal laser by unit cell structure design
,”
Science
293
,
1123
1125
(
2001
).
6.
Y.
Kurosaka
,
S.
Iwahashi
,
Y.
Liang
,
K.
Sakai
,
E.
Miyai
,
W.
Kunishi
,
D.
Ohnishi
, and
S.
Noda
, “
On-chip beam-steering photonic-crystal lasers
,”
Nat. Photonics
4
,
447
(
2010
).
7.
W.
Zhou
,
M.
Dridi
,
J. Y.
Suh
,
C. H.
Kim
,
D. T.
Co
,
M. R.
Wasielewski
,
G. C.
Schatz
,
T. W.
Odom
 et al, “
Lasing action in strongly coupled plasmonic nanocavity arrays
,”
Nat. Nanotechnol.
8
,
506
(
2013
).
8.
K.
Hirose
,
Y.
Liang
,
Y.
Kurosaka
,
A.
Watanabe
,
T.
Sugiyama
, and
S.
Noda
, “
Watt-class high-power, high-beam-quality photonic-crystal lasers
,”
Nat. Photonics
8
,
406
411
(
2014
).
9.
D.
Zhao
,
S.
Liu
,
H.
Yang
,
Z.
Ma
,
C.
Reuterskiöld-Hedlund
,
M.
Hammar
, and
W.
Zhou
, “
Printed large-area single-mode photonic crystal bandedge surface-emitting lasers on silicon
,”
Sci. Rep.
6
,
18860
(
2016
).
10.
A.
Kodigala
,
T.
Lepetit
,
Q.
Gu
,
B.
Bahari
,
Y.
Fainman
, and
B.
Kanté
, “
Lasing action from photonic bound states in continuum
,”
Nature
541
,
196
(
2017
).
11.
S. T.
Ha
,
Y. H.
Fu
,
N. K.
Emani
,
Z.
Pan
,
R. M.
Bakker
,
R.
Paniagua-Domínguez
, and
A. I.
Kuznetsov
, “
Directional lasing in resonant semiconductor nanoantenna arrays
,”
Nat. Nanotechnol.
13
,
1042
1047
(
2018
).
12.
S.-L.
Chua
,
Y.
Chong
,
A. D.
Stone
,
M.
Soljačić
, and
J.
Bravo-Abad
, “
Low-threshold lasing action in photonic crystal slabs enabled by Fano resonances
,”
Opt. Express
19
,
1539
1562
(
2011
).
13.
S.
Wuestner
,
A.
Pusch
,
K. L.
Tsakmakidis
,
J. M.
Hamm
, and
O.
Hess
, “
Overcoming losses with gain in a negative refractive index metamaterial
,”
Phys. Rev. Lett.
105
,
127401
(
2010
).
14.
R.
Marani
,
A.
D'Orazio
,
V.
Petruzzelli
,
S. G.
Rodrigo
,
L.
Martín-Moreno
,
F. J.
García-Vidal
, and
J.
Bravo-Abad
, “
Gain-assisted extraordinary optical transmission through periodic arrays of subwavelength apertures
,”
New J. Phys.
14
,
013020
(
2012
).
15.
M.
Dridi
and
G. C.
Schatz
, “
Model for describing plasmon-enhanced lasers that combines rate equations with finite-difference time-domain
,”
JOSA B
30
,
2791
2797
(
2013
).
16.
J.
Cuerda
,
F.
Rüting
,
F.
García-Vidal
, and
J.
Bravo-Abad
, “
Theory of lasing action in plasmonic crystals
,”
Phys. Rev. B
91
,
041118
(
2015
).
17.
S.
Droulias
,
A.
Jain
,
T.
Koschny
, and
C. M.
Soukoulis
, “
Novel lasers based on resonant dark states
,”
Phys. Rev. Lett.
118
,
073901
(
2017
).
18.
P.
Glendinning
,
Stability, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations
(
Cambridge University Press
,
1994
), Vol.
11
.
19.
S.
Burkhardt
,
M.
Liertzer
,
D. O.
Krimer
, and
S.
Rotter
, “
Steady-state ab initio laser theory for fully or nearly degenerate cavity modes
,”
Phys. Rev. A
92
,
013847
(
2015
).
20.
D.
Liu
,
B.
Zhen
,
L.
Ge
,
F.
Hernandez
,
A.
Pick
,
S.
Burkhardt
,
M.
Liertzer
,
S.
Rotter
, and
S. G.
Johnson
, “
Symmetry, stability, and computation of degenerate lasing modes
,”
Phys. Rev. A
95
,
023835
(
2017
).
21.
H. E.
Türeci
,
A. D.
Stone
, and
B.
Collier
, “
Self-consistent multimode lasing theory for complex or random lasing media
,”
Phys. Rev. A
74
,
043822
(
2006
).
22.
L.
Ge
,
R. J.
Tandy
,
A. D.
Stone
, and
H. E.
Türeci
, “
Quantitative verification of ab initio self-consistent laser theory
,”
Opt. Express
16
,
16895
16902
(
2008
).
23.
S.
Esterhazy
,
D.
Liu
,
M.
Liertzer
,
A.
Cerjan
,
L.
Ge
,
K.
Makris
,
A.
Stone
,
J.
Melenk
,
S.
Johnson
, and
S.
Rotter
, “
Scalable numerical approach for the steady-state ab initio laser theory
,”
Phys. Rev. A
90
,
023816
(
2014
).
24.
Y.
Liang
,
T.
Okino
,
K.
Kitamura
,
C.
Peng
,
K.
Ishizaki
, and
S.
Noda
, “
Mode stability in photonic-crystal surface-emitting lasers with large κ 1 D L
,”
Appl. Phys. Lett.
104
,
021102
(
2014
).
25.
K.
Iga
and
H.
Li
,
Vertical-Cavity Surface-Emitting Laser Devices
(
Springer
,
2003
).
26.
M.
San Miguel
,
Q.
Feng
, and
J. V.
Moloney
, “
Light-polarization dynamics in surface-emitting semiconductor lasers
,”
Phys. Rev. A
52
,
1728
(
1995
).
27.
P.
Mandel
and
M.
Tlidi
, “
Transverse dynamics in cavity nonlinear optics (2000–2003)
,”
J. Opt. B
6
,
R60
(
2004
).
28.
C. W.
Hsu
,
B.
Zhen
,
S.-L.
Chua
,
S. G.
Johnson
,
J. D.
Joannopoulos
, and
M.
Soljačić
, “
Bloch surface eigenstates within the radiation continuum
,”
Light
2
,
e84
(
2013
).
29.
A. F.
Oskooi
,
D.
Roundy
,
M.
Ibanescu
,
P.
Bermel
,
J. D.
Joannopoulos
, and
S. G.
Johnson
, “
MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method
,”
Comput. Phys. Commun.
181
,
687
702
(
2010
).
30.
A.
Cerjan
,
A.
Oskooi
,
S.-L.
Chua
, and
S. G.
Johnson
, “
Modeling lasers and saturable absorbers via multilevel atomic media in the Meep FDTD software: Theory and implementation
,” arXiv:2007.09329 (
2020
).
31.
H.
Haken
,
Laser Light Dynamics
(
North-Holland, Amsterdam
,
1986
), Vol.
2
.
32.
M.
Tinkham
,
Group Theory and Quantum Mechanics
(
Courier Corporation
,
2003
).
33.
J.
Ohtsubo
,
Semiconductor Lasers: Stability, Instability and Chaos
(
Springer
,
2012
), Vol.
111
.

Supplementary Material