Bi-stability and period-doubling cascade of frequency combs in exceptional-point lasers

1 Introduction

A frequency comb is an optical phenomenon where a system produces a series of equally-spaced spectral lines. Optical frequency combs (OFCs) are essential to optical-frequency synthesizers [1] and precision metrology [2], [3], [4]. In the past decade, OFCs have also been applied in optical communications [5] and quantum computation [6]. Conventionally, OFCs are generated by mode locked lasers [7], optical cavities with nonlinearity [8], [9], [10], and quantum cascade lasers [11], [12], [13]. However, these traditional comb generation methods all require the repetition rate to match the free spectral range (FSR) of the laser cavity. Consequently, large cavity sizes are required to generate radio-frequency OFCs.

Recently, it was discovered that exceptional point (EP) optical cavities can develop into frequency combs with repetition rates independent of the cavity’s FSR [14], [15], [16]. EPs are points of degeneracy in the phase space where two or more eigenmodes become identical [17], [18], [19], [20], [21]. When two cavity modes are sufficiently close to an EP, there exists a gain threshold above which any perturbation to the system induces a long-lived periodic modulation to the carrier populations in the gain medium. These dynamic populations subsequently modulate any active lasing modes in the system, generating equally spaced comb lines in the output spectrum. As such, an EP comb is self-generated without any external modulation, and its repetition rate ω _d, equal to the self-modulation rate of the inversion, can be significantly smaller than the cavity’s free spectral range as it is approximately set by the frequency spacing of the EP modes. In principal, arbitrarily small repetition rates can be achieved as the laser system is tuned sufficiently close to the EP. However, in practice, it is technically difficult to reach an exact EP in experiments. Moreover, the robustness of the comb will be compromised due to the enhanced sensitivity near EP [16], [21], [22]. Therefore, the minimum achievable repetition rate in EP combs is thought to be limited.

In this work, we show that one can significantly reduce the repetition rate of EP combs without pushing the system closer to the EP, thus realizing robust radio-frequency OFCs in small laser cavities. Specifically, we demonstrate that an EP comb can halve its repetition rate multiple times through a period-doubling cascade, which was typically observed in more complicated laser systems with light injection or an external modulation [23], [24]. In particular, we carry out a stability analysis with the perturbation method on the Maxwell-Bloch equations [25], [26]. Driven by the periodic population inversion, an EP laser becomes a Floquet system [27], where any infinitesimal perturbation can be decomposed into Floquet eigenmodes, with each having a complex Floquet frequency. A Floquet frequency with a positive (negative) imaginary part leads to the corresponding Floquet mode growing (decaying) over time. The stability of an EP comb is thus determined by the sign of Im(ω _F), with ω _F being the Floquet frequency with the largest imaginary part. Given an existing EP comb with a line spacing of ω _d, by solving for ω _F, we find a series of pumping thresholds at which a Floquet mode turns on, with Im(ω _F) = 0 and Re(ω _F) = 0.5ω _d. Through gain saturation, this Floquet mode induces an additional modulation to the population inversion, which has twice the period of the inversion fluctuation from the original EP comb. The re-modulated inversion then doubles the period of the lasing field’s envelope, hence inserting extra lines into the original EP comb. The period doubling occurs through each of these thresholds cascade, eventually leading to arbitrarily small repetition rate. With the perturbation method, we also find a bistability zone in EP lasers, where two different EP combs exist at the same pumping strength. Thus, the laser state depends on initialization, which is potentially applicable in optical signal-processing devices and all-optical computer systems [28].

2 Stability analysis on frequency combs

Period doubling is a special transition from one frequency comb to another. For the transition to occur, the former comb must first become unstable. To predict the stability of a comb, we apply the first-order perturbation method to the fundamental equations governing lasing materials. We begin by deriving generalized perturbation equations for arbitrary lasing states. Then, we extend the analysis to limit-cycle lasing states, proving that the stability of a frequency comb is associated with a complex Floquet frequency ω _F, which can be determined by solving a linear eigenmode equation.

Lasers can be described rigorously by the Maxwell–Bloch (MB) equations [16], [25], [26], a semi-classical model depicting the relations among the population inversion D(r, t) of gain media, the electric field E(r, t) and the polarization density P(r, t). To simplify the notation, we focus on a one-dimensional (1D) laser cavity with E ( r , t ) = E ( x , t ) z ̂ , P ( r , t ) = P ( x , t ) z ̂ and D(r, t) = D(x, t) (Our method can be readily extended to three-dimensional systems, as shown in Supplementary section I). In particular, the MB equations are

D, E, and P here have been normalized by R ²/(ɛ ₀ ℏγ _⊥), 2 R / ( ℏ γ ⊥ γ ‖ ) , and 2 R / ( ε 0 ℏ γ ⊥ γ ‖ ) , respectively, with R being the amplitude of the atomic dipole moment, ɛ ₀ the vacuum permittivity, ℏ the Planck constant, γ _‖ the population relaxation rate and γ _⊥ the dephasing rate of the gain-induced polarization (i.e., the bandwidth of the gain). D _p(r) is the normalized net pumping strength and profile, ω _ba is the frequency gap between the two atomic levels, ɛ _c (x) is the relative permittivity profile of the cold cavity, σ(x) is a conductivity profile that produces linear absorption, and c is the vacuum speed of light.

Consider a fixed-point or limit-cycle solution to the MB equations (1)–(3), D = D _s(x, t), E = E _s(x, t) and P = P _s(x, t). To determine the stability of the solution, we add a small perturbation, such that D = D(x,t)_s + Δd(x, t), E = E _s(x, t) + Δϵ(x, t), and P = P _s(x, t) + Δp(x, t), where Δ is a real infinitesimal number. The perturbation equations are derived by substituting the perturbed D, P and E into Eqs. (1)–(3) and then extracting the linear terms of Δ,

Notably, for small D _p, the laser is off, hence E _s = 0, P _s = 0 and D _s = D _p. At such a trivial state, equations (4)–(6) yield linear-cavity wave equations [29],

where Γ ⊥ ( ω ) = γ ⊥ ω − ω b a + i γ ⊥ . Equation (7) determines the linear cavity’s resonant modes ϵ _m and the related resonant frequencies ω ̃ m . An EP is approached when two resonant modes merge.

Above the first lasing threshold, a single mode turns on. The solution is a non-trivial fixed point, E s = E 0 ( x ) e − i ω 0 t , P s = P 0 ( x ) e − i ω 0 t and D _s = D ₀(x). Fixed-point stability analysis has been applied to the study of the sensitivity and signal-to-noise ratio of nonlinear EP sensors [30]. Under the stationary-inversion approximation (SIA), Eqs. (4)–(6) yield active-cavity modes and determine the second lasing threshold [29].

Beyond the SIA, Eqs. (4)–(6) determine the comb threshold at which single-mode lasing becomes unstable and the system transitions to a limit cycle [16],

where the repetition rate ω _d, spectral center ω ₀ and the Fourier components {D _m , E _m , P _m } can all be determined by “periodic-inversion ab initio laser theory” (PALT) [16]. Eqs. (8)–(10) act as a periodic temporal modulation in Eqs. (4)–(6). By Floquet theory [27], the solution to Eqs. (4)–(6) should include Floquet modes, f k ( x , t ) e − i ω F , k t , where f _k (x, t) has a period of 2π/ω _d. Due to the difference-frequency generation terms involving ϵ* and p* in Eq. (4), each harmonic oscillation e − i ω F , k t generates a complex-conjugate term e i ω F , k * t . Hence, we postulate the following form of solutions,

Here, f _ak,bk with f = ϵ, p, d are functions of x and t, with a time period of 2π/ω _d, f _ak,bk = f _ak,bk (x, t) = f _ak,bk (x, t + 2π/ω _d). To determine ω _F and the space-time profile of ϵ _a,b (x, t), we substitute Eqs. (8)–(10) and Eqs. (11)–(13) into Eqs. (4)–(6). By expanding the periodic functions in Fourier series, we derive the following wave equations,

in which ϵ ̄ a k , b k ( x ) is a column vector that includes the Fourier components of ϵ _ak,bk (x, t). Ω(ω _F,k ) and X(ω _F,k ) are matrices determined by the limit cycle, {E _s, P _s, D _s}. (d ²/dx ²) acts element-wise on the column vector. d 2 ϵ ̄ a k / d x 2 isssss associated with outgoing boundaries while d 2 ϵ ̄ b k / d x 2 is associated with incoming boundaries. The derivation of Eq. (14) and the expressions for Ω(ω _F,k ) and X(x, ω _F,k ) are provided in the supplementary section I.

Floquet modes can be solved from Eq. (14). Their Floquet frequencies ω _F,k are determined such that the eigenvectors ϵ ̄ a k , b k ( x ) are non-trivial. For each ω _F,k , the expressions of Eq. (11)–(13) suggest that ω F , k ′ ≡ − ω F , k * and ω _F,k ″ ≡ ω _F,k + mω _d with m ∈ Z are degenerate solutions to the same Floquet mode. Therefore, all the Floquet frequencies can be mapped into the “Floquet zone”, defined by Re(ω) ∈ [0, 0.5ω _d ] on the complex plane. In the Floquet zone, we define the primary Floquet frequency ω _F as the ω _F,k with the largest imaginary part, then denote the related mode profile in the expression of Eq. (11) as ϵ _a,b (x, t). A comb is stable if and only if all perturbations decay over time, which requires Im(ω _F,k ) < 0 for all k, equivalently Im(ω _F) < 0.

3 Bistability and period doubling cascade of EP combs

While the imaginary part of the primary Floquet frequency ω _F generally determines the stability of a frequency comb, the real part implies how the comb evolves after it becomes unstable. Specifically, period doubling occurs when (i) Im(ω _F) = 0 and (ii) Re(ω _F) = 0.5ω _d. Under condition (i), a random perturbation reduces to a single Floquet mode over time, ϵ ( x , t ) → e − i ( ω 0 + ω F ) t g ( x , t ) , where g ( x , t ) = ϵ a + ϵ b * e 2 i ω F t . A stronger pumping will lead to the nonlinear effect of wave mixing between the growing Floquet mode and the original comb. If Re(ω _F)/ω _d is irrational, such wave mixing process will generate an overcomplicated lasing spectrum spreading over all frequencies, illustrated in Figure 1a. However, condition (ii) indicates that g(x, t + 2π/ω _d) = g(x, t), hence ϵ(x, t) consists of Fourier components at ω ₀ + (m + 0.5)ω _d, right in the middle of the existing comb lines. Therefore, wave mixing only generates frequencies of ω ₀ + 0.5mω _d , forming another frequency comb with half the repetition rate as before, as shown by Figure 1b. The fact that ω _F reaches exactly 0.5ω _d in condition (ii) is not a coincidence. As discussed in Section 2, both ω _F and − ω F * + ω d are eigenvalues of Eq. (14). These two frequencies are necessarily symmetric about the Floquet zone boundary at 0.5ω _d, and are therefore constrained to appear as a pair. Consequently, when ω _F reaches the zone boundary, it must collide with its partner. As the pumping strength increases further, this degeneracy splits in one of two possible ways: horizontally or vertically. In the case of horizontal splitting, the two solutions acquire different real parts while maintaining the same imaginary part. They remain paired and generically move away from the zone boundary, yielding irregular dynamics if they reach the real axis as illustrated in Figure 1a. In contrast, vertical splitting keeps both solutions pinned to the zone boundary with equal real parts but different imaginary parts, such that each solution becomes its own partner. We denote these vertically separated solutions as ω _F and ω _F,↓ in Figure 1b, where ω _F is the one that eventually crosses the real axis. Any deviation from this behavior – for example, if one branch were to leave the zone boundary – would require the appearance of additional solutions such as − ω F * + ω d or − ω F , ↓ * + ω d , violating the conservation of eigenvalue count ω _F,k as defined in Eq. (14).

We demonstrate such period doubling mechanism in a one-dimensional EP laser cavity shown in Figure 2a. We adopt a smooth pumping profile D _p(x) = 0.5D _max[1 − cos(2πx/L)] to improve the accuracy of the finite-difference time-domain (FDTD) simulations. Figure 2b shows the trajectories of two eigen frequencies solved from Eq. (7); they are tuned near an EP at the first threshold D max = D 1 th , where ω ̃ 0 reaches the real axis, as shown in Figure 2b. Without nonlinear gain saturation, ω ̃ 0 would move quickly upward as D _max increases above D 1 th , while ω ̃ 1 would almost stay still, due to the counteraction between being pumped and being repelled by ω ̃ 0 .

Figure 1:

The Floquet frequencies ω _F and the spectra resulting from ω _F > 0. The original comb lines are located at ω _m = ω ₀ + mω _d for m = 0, ±1, …, with ω ₀ being the comb center and ω _d the repetition rate. In the left column, pink arrows are the Floquet frequency ω _F and its partner eigenvalues in other zones, namely − ω F * + k ω d and (ω _F + kω _d), evolving with the pumping strength. (a) ω _F crosses the real axis at an irrational ratio to ω _d, resulting in a randomly fluctuating spectrum. (b) ω _F reaches the zone boundary (0.5ω _d) simultaneously with − ω F * + ω d below the real axis, then splits into two non-degenerated solutions, which we label as ω _F and ω _{F, ↓} , respectively. ω _F crosses the real axis along the zone boundary, while ω _{F, ↓} is repelled down. Such behavior leads to another comb with half the repetition rate as the original comb.

Figure 2:

The EP laser with gain-loss coupled cavity. (a) A 1D laser cavity with gain on the left and loss on the right. The length of the gain cavity is L and the length of the loss cavity is L ₂ = 1.2L. The absorption rate inside the loss cavity is σ/ɛ ₀ = 7.6(c/L). The passive refractive index ε c is 3.4 in the gain cavity and 3.67 in the passive cavity, typical values of commonly used semiconductor lasing materials. x ₀ is the location where we plot the field intensity in Figure 3 and the phase portraits in Figure 5. The values of all passive parameters, including the thickness and the index of each layer, are elaborated in supplementary section II. (b) the trajectories of two near-EP resonant frequencies solved from Eq. (7). The first lasing threshold is D 1 th = 1.2 . Dashed lines show the would-be above-threshold trajectories in the absence of gain saturation.

For this system, the frequency gap | Re ( ω ̃ 1 − ω ̃ 0 ) | = 0.03 ( c / L ) at the first threshold is an estimation to the repetition rate of the EP comb near the comb threshold. It is approximately 30 times smaller than the gain cavity’s FSR, ω FSR = π c / ε c ≃ 0.9 ( c / L ) . Without tuning the system closer to the EP, we now show how the repetition can be significantly reduced through period-doubling cascade. Figure 3a shows the PALT calculation of the laser states depending on the pumping strength. The upper panel shows a continuous transition from the single lasing mode D 1 th < D max ≤ D c th to the major comb line E ₀(x ₀) across the comb threshold D c th . The lower panel shows the corresponding repetition rate of the EP combs above D c th . Different comb solutions are labeled as “C”-branches with different colors. Figure 3b shows the solutions of primary Floquet frequencies ω _F for the EP combs in Figure 3a. Stable combs are associated with Im(ω _F) < 0. After an ω _F crosses the real axis, the corresponding comb transits from one branch to another.

Figure 3:

Pump dependence of EP combs. (a) PALT calculation of EP combs. (Upper panel) The intensity of central comb line and (lower panel) the repetition rate of EP combs. (Solid lines) Stable solutions and (dashed lines) unstable solutions. (Blue diamonds) the boundaries of bistability. (circles) Period doubling, where ω _d drops by half. (b) The imaginary part (upper panel) and real part (lower panel) of primary Floquet frequencies ω _F solved from Eq. (14) for the combs in a. The inset is the trajectory of C2’s Floquet frequency ω _F in the complex plane, where arrows indicate the increasing of D _max. The result corresponds the case of vertical splitting illustrated in Figure 1b.

The filled circles in Figure 3a mark the thresholds of period doubling, corresponding to the points where ω _F = 0.5ω _d in Figure 3b. Above each of these thresholds, the repetition rate reduces by half, shown in the lower panel of Figure 3a. Figure 4 shows the PALT calculation and FDTD simulation of the lasing spectra at several different pumping strengths selected from Figure 3. The extra comb lines induced by the Floquet mode can be identified on the comb spectrum at D _max = 4.8, compared to the spectrum at D _max = 4.0. The two spectra also show similar comb bandwidth, implying that the emergence of new comb lines does not compromise the intensity of the existing lines.

As the number of comb lines increases, solving PALT and Eq. (14) becomes increasingly time-consuming due to the large number of Fourier components required to maintain accuracy. As a result, although both PALT and the Floquet-theorem-based stability analysis are theoretically sound, they become impractical beyond C4. In this regime, FDTD simulations reveal a continuous spectrum at D _max = 8.0, which can arise from one of two possible evolution paths: either an infinite cascade of period doublings or a random crossing of the real axis by ω _F, as illustrated in Figure 1a. Additional simulation results provided in the Supplementary Section III demonstrate two more period doublings (i.e., C5–C6). Beyond this point, the repetition rate becomes too small to resolve within the finite simulation times considered here. A future theoretical framework is needed to ultimately distinguish between these two evolution paths.

In addition to period doubling cascade, our stability analysis also predicts bistable EP combs within the pumping range of 2.8 < D _max < 5.6. In this region, both C1 and the C2–C3–C4 chain are stable. The two stable branches are accessed by different initial conditions. The accessibility is demonstrated by FDTD simulations in Figure 4. First, we set the pumping strength to be D max = 2.0 > D c th and initialize the simulation with a random pulse of the electrical field inside the gain cavity. After the lasing state converges onto the limit cycle, we increase the pumping strength by a small step, and then continue the simulation. We iterate such simulation process until D _max = 8.0. The simulated comb stays on C1 branch (lower row of the spectra in Figure 4), then suddenly jumps up onto C4 (indicated by the right arrow) as D _max crosses the right end of C1 at 5.6. Second, we start the simulation at D _max = 8.0, then regularly reduce the pumping strength. The simulated comb changes backwardly along C4 → C3 → C2, then suddenly jumps down onto C1 as D _max crosses the left end of C2 at 2.8 (indicated by the left arrow). Thus, the simulation results demonstrate both PALT and the stability analysis on EP combs.

Figure 4:

PALT and FDTD simulation of EP-comb spectra at different pumping strengths. Upper row shows the period-doubling cascade. The bi-stable region is highlighted in orange. Orange arrows show two different paths of how an EP comb evolves. When the pumping D _max increases regularly from 2.0 to 8.0, the simulation converges to the spectra on the lower row for D _max = 4.0 and 4.8, then jumps to C4 spectrum on the upper row. When the pumping D _max reduces regularly from 8.0 to 2.0, the simulation converges to the upper row for D _max = 4.0 and 4.8, then jumps down to C1.

Finally, we summarize the period doubling and bistability phenomena using the phase trajectories of EP combs. The laser’s phase space is defined as a manifold with dimensions of D as well as the real and imaginary parts of E and P from Maxwell–Bloch equations Eqs. (1)–(3). A frequency comb is then recognized as a limit cycle in the phase space [31]. For the solutions in Figures 3a and 5 plots the projections of their phase portraits on the D(x ₀, t)–Re[E(x ₀, t)] plane. The red dot in the upper left plot is the single-mode lasing state, known as a fixed point. The first row shows how the fixed point opens and continuously extends into a stable limit cycle. The second row shows the limit cycle shifting from an attractor (right) to a repeller (middle), then reverting to an attractor (left) again. This gives rise to bi-stability. The third row shows period-doubling cascade, where the limit cycle with a basic period doubles its orbit twice.

Figure 5:

The projection of phase portraits on the plane of population inversion versus electrical field at x = x ₀ . The horizontal axis is the population inversion D _s (x ₀, t) and vertical axis is the envelope of the electrical field E _s . Line styles and colors are consistent with Figure 3.

4 Conclusions

In this work, we develop a limit-cycle stability analysis based on Floquet theory. The analysis predicts the novel phenomena of bi-stability and period-doubling cascade of frequency combs in EP lasers. Period doubling cascade occurs when the Floquet frequency crosses real axis at half of the comb’s repetition rate. It reduces the repetition rate without shrinking the comb bandwidth, hence allowing for the design of extremely compact OFC generators. The theoretical results are confirmed by FDTD simulations. The four-wave mixing phenomenon observed in a laser diode coupled to a high-Q resonator [32] may represent an experimental demonstration of such period-doubling cascade in EP combs. This system possesses all the essential ingredients of an EP: two modes with similar frequencies (one from the laser diode and the other from the high-Q micro-ring cavity), field coupling via backscattering, and gain-loss compensation (gain from the lasing material and loss from the passive cavity).

Our perturbation analysis can be extended to solve scattering problems of periodically-driven nonlinear systems. By including a term of incident wave in Eq. (11), one can re-derive Eq. (14) with an extra inhomogeneous source term. The Flouqet frequency will then be recognized as ω _in − ω ₀, where ω _in is the frequency of the incident wave. Future work can study the scattering spectrum of EP lasers by solving such inhomogeneous perturbation equation.