1. 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 $\chi ^{(2)}$ nonlinearity [7] and bi-chromatic pump [8–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–11]. For the sake of brevity, we refer below to the carrier–envelope offset frequencies as the offset frequencies or offsets. 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. Soliton instabilities happen via the oscillatory scenario leading to the formation of soliton breathers [13–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.

2. Fast and slow solitons

We assume that the electric field, ${{\cal E}}$, inside the resonator is

Soliton frequency combs realise a particular case of Eq. (2) such that

 

Fig. 1. (a) Soliton peak power as a function of the frequency detuning $\delta /\kappa$ for the fixed laser power ${{\cal W}} = 29.5$mW. (b) Difference between the soliton repetition rate, $\tilde D_1$, and the bare resonator rate, $D_1$. Full black lines in (a), (b) show stable solitons. Dashed-black and dashed-blue lines show unstable solitons. The insets in (a) show the stationary (unstable) soliton profiles vs $\theta /\pi$ at points (2) and (3). Panels (c)-(h) show the soliton breathers and their radio frequency (RF) spectra at points (1), (2), and (3). $2\pi /\Omega$ is the breather period. $\Omega$ is given in units of $\kappa$.

Download Full Size | PDF

The relatively large angular soliton width, $\approx 2\pi /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 $\tilde D_1- D_1$ vs $\delta$ plot in Fig. 1(b). Sufficiently rapid change of $\tilde D_1- D_1$ around $\delta /\kappa =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 $\simeq 20$mW.

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 $D_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. 2. Space-time trajectories of the slow (a) and fast (b) solitons plotted in the reference frame rotating with the rate of the bare resonator, $D_1$. The laser power $\mathcal {W}=29.5$mW and $\delta /\kappa =4.6$. (c) Random switching between the fast (longer time intervals) and slow (shorter intervals) for the two solitons simultaneously present in the resonator.

Download Full Size | PDF

3. Soliton breathers and combs with multiple offsets

The unstable soliton branches marked with the blue dashed lines represent the so-called low-branch solitons. If the model is initialised with these solitons, it collapses to the single mode state, $\mu =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 $\delta$’s) or extend along the entire resonator circumference (large $\delta$’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), see Figs. 1(a),(b), is shown in Figs. 1(c),(d) and (e). The reference frame is chosen to rotate with the breather rate, which is generally different from the rate of the unstable stationary solitons taken for the same parameters. Figures 1(f)-(h) show the RF (power) spectra computed in the rotating frame, $\int _0^{\tau } \sum _\mu | \psi _\mu |^{2} e^{-i\Delta t}dt/\tau$. Points (2) and (3) are nearest to the stable solitons, and the RF spectra are near harmonic. Here, the breather frequency $\Omega$ practically coincides with the imaginary part $\beta$ of the eigenvalues driving instabilities of the stationary solitons, $e^{\lambda t-i\beta t}$, 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 $\Omega =\beta /2$. The primary oscillations are still happening with $\beta$, 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,

 

Fig. 3. Illustration of the breather spectrum. The breathing period is $2\pi /\Omega$. $m=0$ is the primary comb and $m=\pm 1$ are the sub-combs. The repetition rate, $\tilde D_1$, is the same for the primary and all sub-combs, while the offset frequencies differ by $m\Omega$.

Download Full Size | PDF

The colour map in Fig. 4(a) shows the dB levels of the $|\hat \psi _{\mu,m}|^{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 $\Omega$. The upward, $\beta _\mu ^{+}$, and downward, $\beta _\mu ^{-}$, parabolas, and the sparsely separated circles for larger $\mu$ show the spectrum of the stable small amplitude excitations, $e^{(\lambda _\mu ^{\pm }-i\beta _\mu ^{\pm })t}$, $\lambda _\mu <0$, existing on top of the single mode, $\mu =0$, state,

 

Fig. 4. Panel (a) shows the two-dimensional map representing the power spectrum of the soliton breather, $|\hat \psi _{\mu,m}|^{2}$, computed at the point (2) in Fig. 1, see Eq. (5). The horizontal lines are separated by the breather frequency $\Omega /\kappa =3.19$. The small panels below show the full (b), the primary (c), $|\hat \psi _{\mu,0}|^{2}$, and the higher order, $|\hat \psi _{\mu,|m|>0}|^{2}$, sub-combs; see Eq. (7). White circles in (a) show the spectrum of excitations around the single mode solution, $\mu =0$, see Eq. (8). The colour map interpolation obscures the discreteness of the spectrum in (a).

Download Full Size | PDF

Data in Fig. 4(a) bring us to an immediate conclusion that, in this instance, the breather frequency $\Omega$ equals $\beta _0^{+}=\sqrt {3(g-\delta )(g-\delta /3)}$. Hence, the primary physical reason for breathing comes from the resonance excitation of the group of modes around $\mu =0$. This is further corroborated by the spectral shapes of the $m=\pm 1$ sub-combs, where the central group of modes is either dominant, see Fig. 4(d), or more pronounced than the $\mu =-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 $|\hat \psi _{\mu,m}|^{2}$ in Fig. 5(a) shows that the breather frequency $\Omega$ is now much larger and is given by $\beta _{52}^{+}\approx \beta _{20}^{+}\approx \beta _{-42}^{+}$. The shape of the $m=\pm 1$ sub-combs in Figs. 5(d), (e) is unambiguous that spectral groups around $\mu =-42,20$, and $52$ are now dominant in the breather structure, and the spectrum around $\mu =0$ is practically absent.

 

Fig. 5. Panel (a) shows the two-dimensional map representing the power spectrum of the soliton breather, $|\hat \psi _{\mu,m}|^{2}$, computed at the point (3) in Fig. 1, see Eq. (5). The horizontal lines are separated by the breather frequency $\Omega /\kappa =28.76$. The small panels below show the full (b), the primary (c), $|\hat \psi _{\mu,0}|^{2}$, and the higher order, $|\hat \psi _{\mu,|m|>0}|^{2}$, sub-combs; see Eq. (7). White circles in (a) show the spectrum of excitations around the single mode solution, $\mu =0$, see Eq. (8). The colour map interpolation obscures the discreteness of the spectrum in (a). Note, the difference of the y-axis scale in (a) here and in Fig. 4.

Download Full Size | PDF

The higher-order sub-combs with $|m|\ge 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 $m=0,\pm 1,\pm 2,\pm 3$, and $\pm 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 spectrum in Fig. 1(f). Figures 4–6 show the mode-number spectra, while the optical frequencies, as measured on a spectrum analyser, must be found using Eq. (7). Figure 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 $\Omega$.

 

Fig. 6. Full (a), primary (b) and higher-order sub-combs computed away from the Hopf instability boundary at point (1) in Fig. 1.

Download Full Size | PDF

 

Fig. 7. Selected spectral lines (a-j) and a whole of the numerically computed frequency spectrum (k) corresponding to the mode-number data in Fig. 5. Zooming on the individual spectral lines reconstructs the schematic illustration in Fig. 3. The red lines in (a-j) correspond to the primary $m=0$ comb. The green and blue lines are the $m=1$ and $m=-1$ sub-combs: $\Omega /2\pi =3.6$GHz.

Download Full Size | PDF

4. 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–11,15–20]. What was particularly interesting for us to find and report here is the significant power and bandwidth of the breather sub-combs with different offset frequencies. 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.

Appendix

The starting point of the derivation of our model is the Maxwell equations in the ring-resonator geometry and with Kerr nonlinearity [21,24,25]. The modal amplitudes $\psi _\mu$ and the field envelope $\psi$, see Eq. (2), comply with the coupled-mode equations that can be cast to the following form [25],

(9)$$i\partial_t\psi_\mu=\Big(\omega_\mu-\omega_p-\mu \tilde D_1\Big)\psi_\mu -\frac{i\kappa}{2}\left(\psi_\mu-\hat\delta_{0,\mu}{{\cal H}}\right) -\gamma \int_0^{2\pi}\psi |\psi|^{2} e^{{-}i\mu\theta}\frac{d\theta}{2\pi}.$$

Eq. (9) contains several parameters in addition to the ones defined in the main text. ${{\cal H}}$ is the pump parameter linked to the laser power, ${{\cal W}}$, as ${{\cal H}}^{2}=\eta {{\cal W}}D_1/2\pi \kappa$, where $\eta$ is the dimensionless coupling parameter. $\hat \delta _{\mu ',\mu }=1$ for $\mu =\mu '$ and $\hat \delta _{\mu ',\mu }=0$ otherwise, and $\gamma$ is the nonlinear parameter, $\gamma /2\pi =20$MHz/W. The $\psi _\mu$ representation of the nonlinear term yields [25],

(10)$$i\partial_t\psi_\mu=\Big(\omega_\mu-\omega_p-\mu \tilde D_1\Big)\psi_\mu -\frac{i\kappa}{2}\left(\psi_\mu-\hat\delta_{0,\mu}{{\cal H}}\right) -\gamma \sum_{\mu_1\mu_2\mu_3} \hat\delta_{\mu_1+\mu_2-\mu_3,\mu} \psi_{\mu_1}\psi_{\mu_2}\psi_{\mu_3}^{*}.$$

The stationary soliton solutions, $\psi _\mu =\hat \psi _\mu$, 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, $\hat \psi _\mu$, with small perturbations, $\epsilon _\mu (t)=x_\mu (t)+y_\mu ^{*}(t)$,

(11)$$\psi_\mu(t) = \hat\psi _\mu + x_\mu(t) +y^{*}_\mu( t),$$

and then linearised Eq. (10),

(12)$$\begin{aligned} i(\partial_t+\tfrac{1}{2}\kappa)x_\mu&= \Big(\omega_\mu-\omega_p-\mu \tilde D_1\Big)x_\mu\\ &-\gamma \sum_{\mu_1\mu_2\mu_3} \hat\delta_{\mu_1+\mu_2-\mu_3,\mu} \left( \psi_{\mu_2}\psi_{\mu_3}^{*}x_{\mu_1}+ \psi_{\mu_1}\psi_{\mu_3}^{*}x_{\mu_2}+ \psi_{\mu_1}\psi_{\mu_2}y_{\mu_3} \right),\\ i(\partial_t+\tfrac{1}{2}\kappa)y_\mu&={-}\Big(\omega_\mu-\omega_p-\mu \tilde D_1\Big)y_\mu\\ &+\gamma \sum_{\mu_1\mu_2\mu_3} \hat\delta_{\mu_1+\mu_2-\mu_3,\mu} \left( \psi_{\mu_1}^{*}\psi_{\mu_2}^{*}x_{\mu_3}+ \psi_{\mu_1}^{*}\psi_{\mu_3}y_{\mu_2}+ \psi_{\mu_2}^{*}\psi_{\mu_3}y_{\mu_1} \right). \end{aligned}$$

By substituting $\left \{x_\mu (t),y_\mu (t)\right \}= \left \{X_\mu,Y_\mu \right \}e^{\lambda t-i\beta t}$, we have reduced Eq. (12) to the algebraic $2N\times 2N$ eigenvalue problem, which was solved numerically. The soliton is stable if all $\lambda <0$ and becomes Hopf unstable when a pair of identical $\lambda$’s with $\pm \beta \ne 0$ crosses zero. Near the Hopf threshold, $\beta$ closely matches the breather frequency $\Omega$, see Figs. 1(g), (h). Components of the eigenvectors in physical space are expressed as $X(\theta )=\sum _\mu X_\mu e^{i\mu \theta }$ and $Y(\theta )=\sum _\mu Y_\mu e^{i\mu \theta }$.

Numerically computed eigenvalues of the unstable solitons at points (1), (2) and (3) in Fig. 1(a) are shown in Figs. 8(a)-(c). As one can see, all three cases correspond to the Hopf instability with the pairs of complex conjugate eigenvalues. For points (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 $\mu =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 $\mu =-42,20$, and $52$. The eigenvector analysis supports the conclusions from decomposing the dynamical data into sub-combs described in the main text.

 

Fig. 8. Eigenvalues and eigenvectors of the perturbations of the solitons unstable relative to the breather excitations. The first, second and third rows correspond to the points (1), (2) and (3) in Fig. 1, respectively. Each panel in the first column has a distinct pair of unstable eigenvalues, $\lambda >0$. The second and third columns show the corresponding eigenvectors in the physical and Fourier spaces, respectively.

Download Full Size | PDF

Frequencies of the excitations supported by the single mode, $\psi _{\mu \ne 0}=0$, state, see Eq. (8), are found from the simple two-by-two eigenvalue problem,

(13)$$i\left(\lambda_\mu^{{\pm}}-i\beta_\mu^{{\pm}}+\tfrac{1}{2}\kappa\right)\begin{bmatrix} x_\mu\\ y_{-\mu} \end{bmatrix}=\begin{bmatrix}\displaystyle (\omega_{\mu}-\omega_p-\mu\tilde D_1)-2{g} & \displaystyle -{g} e^{i\phi}\\ \displaystyle{g} e^{{-}i\phi} & \displaystyle -(\omega_{-\mu}-\omega_p+\mu\tilde D_1)+2{g} \end{bmatrix} \begin{bmatrix} x_\mu\\ y_{-\mu} \end{bmatrix}.$$

Here, ${g}=\gamma |\psi _0|^{2}$ and ${g} e^{i\phi }=\gamma \psi _0^{2}$. In this case, each eigenvalue is attributed to a particular pair of sidebands, $\pm \mu$. The ’$\pm$’ superscripts in the left side of Eq. (13) are used to distinguish between the two eigenvalues of the matrix, see Eq. (8).

Funding

Engineering and Physical Sciences Research Council (2119373); H2020 Marie Skłodowska-Curie Actions (812818).

Acknowledgments

Authors thank S.B. Papp for sharing and discussing experimental data.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data presented in this study are openly available at [26].

References

1. S. A. Diddams, K. Vahala, and T. Udem, “Optical frequency combs: Coherently uniting the electromagnetic spectrum,” Science 369(6501), eaay3676 (2020). [CrossRef]  

2. A. L. Gaeta, M. Lipson, and T. J. Kippenberg, “Photonic-chip-based frequency combs,” Nat. Photonics 13(3), 158–169 (2019). [CrossRef]  

3. Q. Li, T. C. Briles, D. A. Westly, T. E. Drake, J. R. Stone, R. B. Ilic, S. A. Diddams, S. B. Papp, and K. Srinivasan, “Stably accessing octave-spanning microresonator frequency combs in the soliton regime,” Optica 4(2), 193–203 (2017). [CrossRef]  

4. M. H. P. Pfeiffer, C. Herkommer, J. Liu, H. Guo, M. Karpov, E. Lucas, M. Zervas, and T. J. Kippenberg, “Octave-spanning dissipative Kerr soliton frequency combs in Si3N4 microresonators,” Optica 4(7), 684 (2017). [CrossRef]  

5. C. Milian and D. V. Skryabin, “Soliton families and resonant radiation in a micro-ring resonator near zero group-velocity dispersion,” Opt. Express 22(3), 3732 (2014). [CrossRef]  

6. V. Brasch, M. Geiselmann, T. Herr, G. Lihachev, M. H. P. Pfeiffer, M. L. Gorodetsky, and T. J. Kippenberg, “Photonic chip-based optical frequency comb using soliton Cherenkov radiation,” Science 351(6271), 357–360 (2016). [CrossRef]  

7. X. Liu, Z. Gong, A. W. Bruch, J. B. Surya, J. Lu, and H. X. Tang, “Aluminum nitride nanophotonics for beyond-octave soliton microcomb generation and self-referencing,” Nat. Commun. 12(1), 5428 (2021). [CrossRef]  

8. S. Zhang, J. M. Silver, T. Bi, and P. Del’Haye, “Spectral extension and synchronization of microcombs in a single microresonator,” Nat. Commun. 11(1), 6384 (2020). [CrossRef]  

9. G. Moille, E. F. Perez, J. R. Stone, A. Rao, X. Lu, T. S. Rahman, Y .K. Chembo, and K. Srinivasan, “Ultra-broadband Kerr microcomb through soliton spectral translation,” Nat. Commun. 12(1), 7275 (2021). [CrossRef]  

10. P. C. Qureshi, V. Ng, F. Azeem, L. S. Trainor, H. G. L. Schwefel, S. Coen, M. Erkintalo, and S. G. Murdoch, “Soliton linear-wave scattering in a Kerr microresonator,” Commun. Phys. 5(1), 123 (2022). [CrossRef]  

11. S. B. Papp (private communication).

12. J. R. Stone and S. B. Papp, “Harnessing Dispersion in Soliton Microcombs to Mitigate Thermal Noise,” Phys. Rev. Lett. 125(15), 153901 (2020). [CrossRef]  

13. K. Nozaki and N. Bekki, “Solitons as attractors of a forced dissipative nonlinear Schrödinger equation,” Phys. Lett. A 102(9), 383–386 (1984). [CrossRef]  

14. D. V. Skryabin, “Energy of internal modes of nonlinear waves and complex frequencies due to symmetry breaking,” Phys. Rev. E 64(5), 055601 (2001). [CrossRef]  

15. C. Milian, A. V. Gorbach, M. Taki, A. V. Yulin, and D. V. Skryabin, “Solitons and frequency combs in silica microring resonators: Interplay of the Raman and higher-order dispersion effects,” Phys. Rev. A 92(3), 033851 (2015). [CrossRef]  

16. E. Lucas, M. Karpov, H. M. Guo, L. Gorodetsky, and T. J. Kippenberg, “Breathing dissipative solitons in optical microresonators,” Nat. Commun. 8(1), 736 (2017). [CrossRef]  

17. M. Yu, J. Jang, Y. Okawachi, A. G. Griffith, K. Luke, S. A. Miller, X. Ji, M. Lipson, and A. L. Gaeta, “Breather soliton dynamics in microresonators,” Nat. Commun. 8(1), 14569 (2017). [CrossRef]  

18. M. Johansson, V. E. Lobanov, and D. V. Skryabin, “Stability analysis of numerically exact time-periodic breathers in the Lugiato-Lefever equation: Discrete vs continuum,” Phys. Rev. Res. 1(3), 033196 (2019). [CrossRef]  

19. C. Bao, B. Shen, M. G. Suh, H. Wang, K. Şafak, A. Dai, A. B. Matsko, F. X. Kartner, and K. Vahala, “Oscillatory motion of a counterpropagating Kerr soliton dimer,” Phys. Rev. A 103(1), L011501 (2021). [CrossRef]  

20. A. A. Afridi, H. Weng, J. Li, J. Liu, M. McDermott, Q. Lu, W. Guo, and J. F. Donegan, “Breather solitons in AlN microresonators,” Opt. Continuum 1(1), 42 (2022). [CrossRef]  

21. T. Herr, V. Brasch, J. D. Jost, C. Y. Wang, N. M. Kondratiev, M. L. Gorodetsky, and T. J. Kippenberg, “Temporal solitons in optical microresonators,” Nat. Photonics 8(2), 145–152 (2014). [CrossRef]  

22. D. J. Kaup and A. C. Newell, “Theory of nonlinear oscillating dipolar excitations in one-dimensional condensates,” Phys. Rev. B 18(10), 5162–5167 (1978). [CrossRef]  

23. I. V. Barashenkov and Y. S. Smirnov, “Existence and stability chart for the ac-driven, damped nonlinear Schrödinger solitons,” Phys. Rev. E 54(5), 5707–5725 (1996). [CrossRef]  

24. Y. K. Chembo and N. Yu, “Modal expansion approach to optical-frequency-comb generation with monolithic whispering-gallery-mode resonators,” Phys. Rev. A 82(3), 033801 (2010). [CrossRef]  

25. D. V. Skryabin, “Hierarchy of coupled mode and envelope models for bi-directional microresonators with Kerr nonlinearity,” OSA Continuum 3(6), 1364 (2020). [CrossRef]  

26. D. N. Puzyrev and D. V. Skryabin, “Dataset related to the article "Frequency combs with multiple offsets in THz-rate microresonators",” University of Bath Research Data Archive (2022), https://doi.org/10.15125/BATH-01206