Multistability and coexisting soliton combs in ring resonators: the Lugiato-Lefever approach

We are reporting that the Lugiato-Lefever equation describing the frequency comb generation in ring resonators with the localized pump and loss terms also describes the simultaneous nonlinear resonances leading to the multistability of nonlinear modes and coexisting solitons that are associated with the spectrally distinct frequency combs.


Introduction
Generation of frequency combs spanning an optical octave [1] and the soliton formation [2][3][4] are the primary nonlinear effects in the microring and in the passive fiber-loop ring resonators currently attracting a growing attention. Both of these effects and both of these systems are closely related and are connected to a variety of applications, such as e.g., all-optical signal processing [5] and high precision spectroscopy [6]. Solitons in the ring resonators circulate indefinitely providing losses are compensated by the continuous wave pumping. On the contrary, the solitons freely propagating in optical fibers and other waveguides are, and in their essence, quasi-solitons since the unbalanced losses and higher order effects lead to the energy radiation [7].
A classic Lugiato-Lefever equation (LLE) [8] is widely used to model the frequency combs and solitons in ring resonators, see, e.g., [2,9]. If one sweeps the pump frequency over the range of several resonances, then the LLE successfully describes the nonlinear tilt and the associated bistability effect only for a single resonance frequency, see, e.g., [2,9,10]. This is an obvious shortcoming since doing such a sweep in an experiment produces a bistable response at every resonance. An alternative approach is to apply the Ikeda map model, which was first used for the solitons in ring resonators a few decades ago [11]. The Ikeda approach was extended to the nonlinear fiber rings [12] and more recently to the microring comb generation, where it predicted several nonlinear resonances, multistability and supersolitons [13]. The title of Ref. [13] explicitly says that these effects are beyond what the LLE can describe. Disadvantages of the Ikeda approach is that the periodic boundary conditions are imposed in the evolution coordinate, which complicates mathematical methods and physical interpretation of even the simplest stability analysis of a homogeneous state [12,13], which is standard and routine in the LLE, see, e.g., [2,10].
The LLE typically utilizes the space independent pump term to describe frequency combs and associated solitons, see, e.g., [2,[8][9][10], while the Ikeda approach intrinsically relies on the pump that is localized at one point [11][12][13]. The salient feature of the latter is that its expansion into the resonator modes has an infinite spectrum of the Dirac-delta function. Thus, in the LLE methods used so far the pump term has nonzero projection on only one mode, while the Ikeda method implies equal pump strength for the entire spectrum of the resonator modes. Clearly, both of these limits are mathematical idealizations, which have their regions of validity. In this work, we demonstrate how the LLE approach can be used to describe the multiple nonlinear resonances leading to the multistability effects and coexisting solitons with different widths and amplitudes that are associated with the spectrally distinct frequency combs.

Lugiato-Lefever equation with the localised pump and loss
The LLE model, well established in the context of the ring resonators [2,3,9,10], reads as It describes the evolution of the dimensionless electric field amplitude ψ in time T and along the spatial coordinate Z. All the parameters and variables in the above equation, apart from Ψ itself, are measured in physical units. The dimensionless amplitude Ψ obeys an obvious and important periodicity condition Ψ(T, Z − L/2) = Ψ(T, Z + L/2), where L is the resonator length. Here κ is the loss rate, P and ω p are the pump rate and the pump frequency, T −1 r |Ψ| 2 is the nonlinear shift of the resonances. The coefficient before |Ψ| 2 has the dimension of the inverse time, but its value can be scaled to any convenient number, through the scaling of P and Ψ only. Our choice was to make it equal to the inverse round trip time, T r = L/ω .
Disregarding the pump, loss and nonlinearity and assuming Ψ = e imβZ−iω m T , where |m| = 0, 1, 2, 3, . . . , one finds for the frequency of the m'th mode ω m = ω 0 +ω mβ+ω m 2 β 2 /2. Here ω is the group velocity and ω is the dispersion coefficient at the frequency of the reference m = 0 mode. ω 0 makes the cavity spectrum slightly non-equidistant. In our case the dispersion is anomalous ω > 0, see, e.g., Refs. [2,10] for the equations connecting the waveguide and the resonator dispersion coefficients. The condition βL = 2π comes from the periodicity of Ψ. β = 2π/L has the meaning of the intermodal distance in the Fourier space reciprocal to the coordinate z, while ω β is the free spectral range in the frequency space. In a typical setup, pump P = P(Z) and losses κ = κ(Z) are the functions of Z that have a sharp maximum in the vicinity of a point where the resonator is coupled to the bus waveguide, which is used to couple in the pump and to couple out the signal. The effects of multistability and co-existence of different solitons studied below are introduced into the LLE model through the spatial inhomogeneity of the pump P and are the primary subjects of this Express communication.
In order to transform our system into a dimensionless form we assume Z = θ L/(2π). Here θ is a new dimensionless coordinate, θ ∈ (−π, +π), which in the case of a ring shaped resonator is exactly the polar angle.
Here h and γ are the maximal pump and loss values at θ = 0, θ 0 2π is the angular width of the pump spot. The θ-dependent separation between the straight pump waveguide and the ring shaped resonator waveguide is well approximated by l + Lθ 2 /(2π), where l is the minimal separation at θ = 0. The coupling strength between the two waveguides is known to decay exponentially with the separation, which gives the Gaussian profile for H(θ). Note, that if one transforms Eq. (2) into a reference frame making one round trip per unit time, then the ∂ θ -term cancels out, but the pump and loss terms become periodic in time, providing their spatial inhomogeneity was accounted for: The dimensionless resonator dispersion parameter d is the ratio between the round trip time and the dispersion time of a pulse having the width of the resonator length: d = 2π 2 L/ω L 2 /ω . The soliton width scales as (∆/d) 1/2 , which for d = 10 −5 allows to fit ∼ 1000 solitons over the cavity length. Here ∆ is the relative detuning from a chosen resonance, 0 < ∆ 2π. Our choice of the relatively high resonator dispersion is convenient to reduce computational requirements, since the higher dispersion increases the width of the solitons relative to the resonator length, but it does not alter general applicability of our results to the resonators with a broad range of ω values, since we can use L as a fitting parameter. The resonator finesse F can be defined as the ratio between the free spectral range (FSR) and the cavity loss rate: F = ω β min κ = 2π γ . We choose γ = 0.05 in what follows, which gives the relatively low F 125. Relatively low F achieved primarily through the small FSR are critical for observing multistability and coexisting comb solitons discussed below. While our aim here is to introduce a proof of principle theoretical results, the explicit scaling we presented above allows to relate our approach to a number of geometries and material choices. A practical implementation of the regimes considered here is likely to be possible, e.g., in the fiber loop resonators [14] and in the low finesse relatively long microrings [15].

Linear and nonlinear modes, multistability and instabilities
We first consider stationary, ∂ τ = 0, nonlinear modes of Eq. (2). We expand ψ(θ, τ) =ψ(θ), Γ(θ) and H(θ) into Fourier series {ψ(θ), Γ(θ), H(θ)} = ∞ m=−∞ {ψ m , γ m , h m }e imθ and obtain the following algebraic system: When losses are uniform and nonlinearity is absent the single mode solutions are given by Plotting ψ m as a function of δ one finds a set resonances separated roughly by 2π, with a small offset due to d 0. Since the amplitude of the m-th resonance is proportional to h m , the resonance peaks are located under an envelope with the spectral width ∼ θ −2 0 . If the pump is uniform, then only the m = 0 resonance has a nonzero amplitude, simply because the pump projection on all other modes is zero. The positive Kerr nonlinearity tilts this resonance towards large δ introducing a bistability loop.
Taking as an example θ 0 = π/50 (variation of θ 0 within the realistic interval of few degrees produces quantitatively similar results) and accounting for 2000 modes, we solved Eqs. (4) numerically using Newton method. Figs. 1(b) and 1(c) show the dependencies of the amplitude of the nonlinear modeψ(θ) = m ψ m e imθ defined as max θ |ψ(θ)| on the detuning parameter δ. One can see the multiple tilted resonances, that start overlapping for the sufficiently strong pump and form what we call the multistability or tristability effect, see Fig. 1(c). Spectra of the pump H and of the nonlinear solutionψ are shown in Fig. 1(a).ψ is generally asymmetric with respect to the θ = 0 point in the real space and with respect to m = 0 in the Fourier space due to the simultaneous presence of the ∂ θ −term in Eq. (1) and of the non-uniformities of H and Γ. At the lowest branches of the nonlinear resonance curves, this asymmetry is almost negligible, but it becomes noticeable for the high amplitude solutions, see Fig. 1(a).
By taking the real-space Eq. (2) we have linearised it aroundψ using the substitution ψ(θ, τ) =ψ(θ) + u(θ)e λτ + v * (θ)e λ * τ with |ψ| |u|, |v|: The intervals of the instabilities calculated numerically from Eqs. (5) are shown in Figs. 1(b) and 1(c) with the red lines. Similarly to the case of the uniform pump and losses, the upper branches of every resonance are unstable [8]. For every δ there is a finite set of the unstable u, v pairs [8]. The spatial profiles of these pairs are dominated by the linear modes e imθ with the modal numbers m spanning a finite interval. Representative spectra of u(θ) = m u m e iθm and v(θ) = m v m e iθm corresponding to the maximal instability growth rate found for a given δ are shown in Fig. 1(d).
The instability in this particular example is dominated by the growth of the m = 549th mode.

Coexisting solitons and frequency combs
It is well known by now that the instability of the upper branch of the nonlinear resonance leads to the generation of the frequency combs through the cascaded four-wave mixing process, see, e.g., [9]. As δ approaches its value corresponding to the rightmost point of the tilted resonance the generally complex space-time comb signals are replaced by sets of spontaneously generated solitons, see, e.g., [2,10]. We observed this dynamics in our system too. However, when the multistability is present, then we have an option of starting a simulation, for a given δ, from either of the two coexisting upper branch solutions. Note, that the stable lowest amplitude solution can also coexist with the two unstable upper branches. The former provides a stable background for the soliton pulses with the peak amplitudes close to either the highest or the intermediate nonlinear resonances. When we performed simulations of the comb generation using Eqs. (2) with 2 16 modes and initializing them with the highest amplitude nonlinear mode found from the time independent Eqs. (4), we have typically observed the dynamics leading to the comb generation through the formation of multiple solitons, see Fig. 2(a). The associated frequency comb is shown in Fig. 2(b). A close examination reveals that the solitons constituting this comb can have two different amplitudes and widths. Extensive numerical simulations show that output field distributions always contain a considerable number of the low-amplitude solitons and a few of the high-amplitude ones. Fig.  3(a) shows a zoomed interval of θ's from Fig. 2(a) capturing the two different solitons stably propagating one next to the other. The taller one is also narrower since it corresponds to the larger effective detuning from the associated linear resonance. Isolating either of the two solitons from the rest of the signal, we have observed that both of them persist in the resonator practically indefinitely. The comb associated with the high amplitude soliton is much broader, cf. Figs. 3(c) and 3(d). If one compares spectra of the individual solitons with the spectrum of the whole signal in Fig. 2(b), one can clearly see the two different scales in tails of the latter. Thus the central part of the spontaneously generated comb is associated with the low-amplitude solitons, while the outer wings are determined by the high-amplitude narrow-width ones.
The comb solitons reported here, both inside and outside the multistability regime, reside on top of the oscillating in space and time background wave and demonstrate pronounced breathing. The breathing soliton amplitude is shown in Fig. 3(b). The oscillation period matches the resonator round trip time and hence is associated with the localizations of the pump and loss terms. However, the maxima and minima of these amplitude oscillations do not generally coincide with the moments when the soliton overlaps with the maximum of H. Also, the amplitude of the oscillations varies slightly even for the solitons of the same type, but located at different points. This temporal dynamics is associated with the modulation of the comb soliton spectra in Figs. 3(c) and 3(d), which is similar to the spectral modulation reported in Ref. [16].

Conclusions
We demonstrated how the LLE model describing the frequency comb generation can be generalized to predict multiple simultaneous nonlinear resonances, multistability and coexistence of the small and high amplitude solitons and of the frequency combs with different spectral widths.