Frequency combs with multiple offsets in THz-rate microresonators

Octave-wide frequency combs in microresonators are essential for self-referencing. However, it is difficult for the small-size and high-repetition-rate microresonators to achieve perfect soliton modelocking over the broad frequency range due to the detrimental impact of dispersion. Here we examine the stability of the soliton states consisting of one hundred modes in silicon-nitride microresonators with the one-THz free spectral range. We report the coexistence of fast and slow solitons in a narrow detuning range, which is surrounded on either side by the breather states. We decompose the breather combs into a sequence of subcombs with different carrier envelope offset frequencies. The large detuning breathers have a high frequency of oscillations associated with the perturbation extending across the whole microresonator. The small detuning breathers create oscillations localised on the soliton core and can undergo the period-doubling bifurcation, which triggers a sequence of intense sub-combs.


Introduction
Broadband optical frequency combs find their applications in optical metrology, spectroscopy and optical information processing [1]. The high repetition rate, small footprint, and possibility of integrating the resonator and pump source on the same chip are the attractive features of the microresonator frequency comb technology [2]. Comb self-referencing requires achieving the octave-wide soliton spectral spans. Methods to generate the broad soliton combs rely on, e.g., reducing the microresonator radius [3,4], higher order dispersion [5,6], frequency doubling in a material with (2) nonlinearity [7] and bi-chromatic pump [8][9][10].
In the microresonators, the comb lines start making their equidistant spectral steps from the pump laser frequency. However, this ideal grid does not always extend across the entire comb. The spectral structure of the octave combs should be carefully examined on a case-by-case basis because it can be composed of the overlapping combs having different repetition rates and different carrier-envelope offset frequencies [8][9][10][11]. For the sake of brevity, we refer below to the carrier-envelope offset frequencies as the offset frequencies or offsets. The Significant changes in the dispersion and growth of specific modes triggering new four-wave mixing cascades can be the reasons behind this unwanted complexity. Therefore, developing more profound insights into the physical mechanisms behind the multiple offset frequencies and repetition rates is beneficial.
Our goal is to present a theoretical framework and associated numerical data using an example of the silicon nitride THz-rate resonator [3,4,12]. We report the coexisting stable soliton combs with different repetition rates, i.e., fast and slow solitons. The soliton instability happens via the oscillatory scenario leading to the formation of soliton breathers [13][14][15][16][17][18][19][20]. The range of the breather existence is broad and happens on either side from the pump detuning interval supporting stable solitons. For small detunings, the unstable waveforms are localised around the soliton core, while the large detuning instabilities are driven by the perturbations smeared along the resonator circumference. We demonstrate that the breather spectra are composed of several, sometimes very intense, overlapping combs having different offset frequencies and sharing the same repetition rates. Mathematics underpinning our formalism is described in Appendix, while the main text is focused on the data analysis.

Fast and slow solitons
We assume that the electric field, E, inside the resonator is where is the pump laser frequency, is time, ∈ [0, 2 ) is the angular coordinate in the laboratory frame, and is the mode number of the resonance nearest to . The envelope function is sought as the mode expansion in the rotating frame of reference, Here, are the complex modal amplitudes, are the relative mode numbers, and is the coordinate in the frame rotating with the repetition rate˜1. We assume here that the bare resonator spectrum, , is approximated by = 0 + 1 + 1 2! , where 0 /2 = 193THz. The dispersion parameters match the THz-rate Si 3 N 4 resonator used in Ref. [12], see Fig. 4 there: 1 /2 = 1THz, 2 /2 = 20MHz, 3 /2 = 0Hz, 4 /2 = −120kHz, Soliton frequency combs realise a particular case of Eq. (2) such that The interplay of nonlinearity and higher-order dispersion makes the soliton repetition rate˜1 to be different from the bare resonator rate value, 1 . The repetition rate difference,˜1 − 1 , depends on the resonator and pump parameters and has been determined numerically self-consistently with the modal amplitudesˆ. Numerical data assembled in Fig. 1 trace the soliton peak power, max | | 2 , and˜1 − 1 as functions of the pump detuning, = 0 − , for the fixed input power W = 29.5mW.
The relatively large angular soliton width, ≈ 2 /10, boosts the role played by the finite size effects. This makes the bifurcation structure of the THz-rate solitons to be more complex than the one in the case of low repetition rate resonators, see, e.g., Fig. 4 in Supplemental Materials of Ref. [21], which corresponds to the damped-driven Nonlinear Shrödinger equation (Lugiato-Lefever model) solved on the effectively infinite spatial interval [22,23]. The complex behaviour of the soliton peak power in Fig. 1(a) corresponds to the more transparent structure of the˜1 − 1 vs plot in Fig. 1(b). Sufficiently rapid change of˜1 − 1 around / = 5 separates the ranges of the slow and fast solitons. Numerical investigation of the power dependencies of the soliton families shows that the fast solitons emerge for powers above 20mW.
The stability of solitons relative to small perturbations is an important issue that we have analysed in detail using the formalism described in Appendix. The dashed and full lines in Figs. 1(a) and 1(b) correspond to the unstable and stable solitons, respectively. It should be noted that the stability intervals of the fast and slow solitons overlap, see the shading in Fig. 1(b). The space-time trajectories of the slow and fast solitons in the reference frame rotating with the bare resonator rate 1 are shown in Figs. 2(a) and 2(b), respectively. If the model is initialised simultaneously with the fast and slow solitons, e.g., located on the opposing sides of the resonator, see Fig. 2(c), then we have observed the convergence of the slow soliton to the fast one. This process could not, however, fully stabilise itself, and the pair was making random switches to the slow state and back to fast. Fig. 3. Illustration of the breather spectrum. The breathing period is 2 /Ω. = 0 is the primary comb and = ±1 are the sub-combs. The repetition rate,˜1, is the same for the primary and all sub-combs, while the offset frequencies differ by Ω.

Soliton breathers and combs with multiple offsets
The unstable soliton branches marked with the blue dashed lines represent the so-called lowbranch solitons. If the model is initialised with these solitons, it collapses to the single mode state, = 0. The solitons shown by the dashed black lines are unstable to the oscillatory (Hopf) instability. The waveforms driving the Hopf instability are either localised on the soliton core (small 's) or extend along the entire resonator circumference (large 's), see Appendix. Typical operation regimes emerging from the initial conditions on the black dashed lines correspond to the soliton breathers; see Figs. 1 (c)-(e).
The breather dynamics computed at points (1), (2) and (3) (2) and (3) are nearest to the stable solitons, and the RF spectra are near harmonic. Here, the breather frequency Ω practically coincides with the imaginary part of the eigenvalues driving instabilities of the stationary solitons, − , see Appendix for details. Point (1) is the furthest from the Hopf threshold and is beyond the period doubling threshold also. Here, the breather frequency becomes Ω = /2. The primary oscillations are still happening with , but the small up and down shifts of the neighbouring intense spots, see Fig. 1(c), double the period and half the breathing frequency.
Since the breather combs have periodic modal amplitudes, the time-domain Fourier expansion ofˆ( ) yields whereˆ, are the time-independent amplitudes of the breather spectrum. Using Eq. (1) and Eq. (5) and rearranging we find that the electric field of a breather is Thus, the breather is a superposition of the primary, = 0, and higher order, ≠ 0, sub-combs which all share the same repetition rate,˜1, but have different central frequencies + Ω, and, therefore, different offsets from zero. The multiple spectral grids of the breather combs are represented by , = + Ω +˜1, = 0, ±1, ±2, . . . , with |ˆ, | 2 being the sideband powers in the 's subcomb, see Fig. 3 for an illustration. The colour map in Fig. 4(a) shows the dB levels of the |ˆ, | 2 spectrum for the breather found at the point (2), i.e., on the left from the soliton stability interval in Fig. 1(a). The discreteness of the horizontal lines in Fig. 4(a) is smoothed over by the colour map interpolation. The vertical separation of these lines equals Ω. The upward, + , and downward, − , parabolas, and the sparsely separated circles for larger show the spectrum of the stable small amplitude excitations, ( ± − ± ) , < 0, existing on top of the single mode, = 0, state, Here, g ∼ | 0 | 2 is the nonlinear shift of the resonance; see Appendix for further details. Data in Fig. 4(a) bring us to an immediate conclusion that, in this instance, the breather frequency Ω equals + 0 = √︁ 3( − ) ( − /3). Hence, the primary physical reason for breathing comes from the resonance excitation of the group of modes around = 0. This is further corroborated by the spectral shapes of the = ±1 sub-combs, where the central group of modes is either dominant, see Fig. 4(d), or more pronounced than the = −40 group, see Fig. 4(e). Now, we are moving onto describing the sub-combs driving the breather dynamics for large detunings, i.e., on the right from the soliton stability interval in Fig. 1(a). The colour map of |ˆ, | 2 in Fig. 5(a) shows that the breather frequency Ω is now much larger and is given by The shape of the = ±1 sub-combs in Figs. 5(d), (e) is unambiguous that spectral groups around = −42, 20, and 52 are now dominant in the breather structure, and the spectrum around = 0 is practically absent.
The higher-order sub-combs with | | ≥ 2 have very small power in the examples presented in Figs. 4, 5. This is because points (2) and (3) are located not deep enough into the instability range, see Fig. 1(a). However, point (1) corresponds to well-developed instability, which has already crossed the period doubling threshold. The higher-order sub-combs become well developed here, see the = 0, ±1, ±2, ±3, and ±4 cases in Fig. 6. Predictably, the even order sub-combs are more powerful because they correspond to a subsequence of the stronger peaks in the RF Fig. 6. Full (a), primary (b) and higher-order sub-combs computed away from the Hopf instability boundary at point (1) in Fig. 1. spectrum in Fig. 1(f). Figs. 4-6 show the mode-number spectra, while the optical frequencies, as measured on a spectrum analyser, must be found using Eq. (7). Fig. 7, however, shows the selected spectral lines and a whole of the numerically computed frequency spectrum corresponding to the mode-number data in Fig. 5. Zooming on the individual lines reconstructs the schematic illustration in Fig. 3. The predicted spectral features are amenable to measurements on the high resolution spectrum analysers. Possible oscillatory instabilities due to interaction of different mode families may yield higher values of Ω.

Summary
Soliton breather states in the damped and driven nonlinear Schrödinger (Lugiato-Lefever) equation have been known for a long time [13,14]. The recent resurgence in studies of their properties has happened primarily through the context of the frequency comb generation in microresonators [8][9][10][11][15][16][17][18][19][20]. What was particularly interesting for us to find and report here is The modal composition of these sub-combs varies between the ones spectrally localised around the pump and others detuned away. These detuned sub-combs are associated with relatively high-frequency breathers. All this complexity happens through the interplay of the large second and higher order dispersions and a relatively small number of modes covering the significant bandwidth, which are all features of the small-size and high-repetition-rate microresonators. Stable solitons can also be found under this extreme spectral broadening. The bifurcation structure of the THz-rate solitons and breathers differs from the one found in the idealised Lugiato-Lefever equation solved on the infinite spatial domain. In addition to the breathers found at the small and large detunings, we have reported the coexisting stable solitons having two different repetition rates, i.e., slow and fast solitons. . (9) Eq. (9) contains several parameters in addition to the ones defined in the main text. H is the pump parameter linked to the laser power, W, as H 2 = W 1 /2 , where is the dimensionless coupling parameter.ˆ , = 1 for = andˆ , = 0 otherwise, and is the nonlinear parameter, /2 = 20MHz/W. The representation of the nonlinear term yields [25], The stationary soliton solutions, =ˆ, see Eq.
(3), have been found using the Newton method. To differentiate between stable and unstable solitons and identify the oscillatory (Hopf) instabilities we have analysed the linear stability of the soliton family. We perturbed the time-independent soliton amplitudes,ˆ, with small perturbations, ( ) = ( ) + * ( ), and then linearised Eq. (10), By substituting ( ), ( ) = , − , we have reduced Eq. (12) to the algebraic 2 × 2 eigenvalue problem, which was solved numerically. The soliton is stable if all < 0 and becomes Hopf unstable when a pair of identical 's with ± ≠ 0 crosses zero. Near the Hopf threshold, closely matches the breather frequency Ω, see Figs (1) and (2), i.e., left from the soliton stability interval, the eigenvectors are localised around the soliton core and are largely shaped by the group of modes around = 0, see Fig. 8. For point (3), i.e., right from the soliton stability interval, the eigenvectors are spread across the entire resonator circumference and spectrally centred around = −42, 20, and 52. The eigenvector analysis supports the conclusions from decomposing the dynamical data into sub-combs described in the main text.
Frequencies of the excitations supported by the single mode, ≠0 = 0, state, see Eq. (8), are found from the simple two-by-two eigenvalue problem, Here, g = | 0 | 2 and g = 2 0 . In this case, each eigenvalue is attributed to a particular pair of sidebands, ± . The '±' superscripts in the left side of Eq. (13) are used to distinguish between the two eigenvalues of the matrix, see Eq. (8).