Two-colour dissipative solitons and breathers in microresonator second-harmonic generation

Abstract

Frequency conversion of dissipative solitons associated with the generation of broadband optical frequency combs having a tooth spacing of hundreds of giga-hertz is a topical challenge holding the key to practical applications in precision spectroscopy and data processing. The work in this direction is underpinned by fundamental problems in nonlinear and quantum optics. Here, we present the dissipative two-colour bright-bright and dark-dark solitons in a quasi-phase-matched microresonator pumped for the second-harmonic generation in the near-infrared spectral range. We also found the breather states associated with the pulse front motion and collisions. The soliton regime is found to be typical in slightly phase-mismatched resonators, while the phase-matched ones reveal broader but incoherent spectra and higher-order harmonic generation. Soliton and breather effects reported here exist for the negative tilt of the resonance line, which is possible only via the dominant contribution of second-order nonlinearity.

Introduction

Optical frequency comb generation in microresonators has attracted significant attention in recent years^1,2. The key results in this area are the demonstration of temporal dissipative Kerr solitons³ and octave-spanning combs suitable for self-referencing^4,5,6. These developments have enabled a broad range of applications, such as optical clocks, coherent optical communication, exoplanet detection, and many others, see e.g.,^7,8,9. Beyond becoming an outstanding test-bed for dissipative solitons, nonlinear and quantum effects in microresonators have made a profound impact in such interdisciplinary areas as pattern formation¹⁰, synchronization of oscillators^11,12, light crystals, and topological physics in space and time^{13,14,15,16,17}.

Dissipative bright solitons and associated frequency combs in microresonators possessing Kerr nonlinearity require anomalous group-velocity dispersion (GVD) at the pump wavelength^1,3, and the dark ones are observed in the normal GVD regime¹⁸. Simultaneous bright and dark Kerr soliton pairs spectrally located on different sides of the zero GVD wavelength have also been recently observed¹⁹. Normal-dispersion Kerr resonators with the modulated circumference of the inner ring, which couples forward and backward waves, have been used to demonstrate the continuum of bright and dark dissipative Kerr soliton states²⁰.

The need to expand the family of microresonator combs in the visible and mid-infrared (IR) ranges stimulates interest in modelocking involving simultaneous harmonic generation, e.g., using χ⁽²⁾, i.e., second-order nonlinearity. χ⁽²⁾ effect allows generating combs at twice or half of the pump frequency at the comparatively low input powers^{21,22,23,24,25,26,27}. An attractive feature of microresonator χ⁽²⁾ combs is their immediate octave width, which offers a pathway to compact self-referencing arrangements in the integrated setups⁶.

Reliable generation of χ⁽²⁾ solitons in microresonators remains a challenge. The existence of such solitons assumes, as a necessary condition, mutual modelocking of the two groups of modes located around the pump frequency and either half- or second-harmonic²⁸. One of the obstacles is, therefore, the large accumulated dispersion across the octave bandwidth. As a result, the group-velocity difference between the two modal groups dominates the nonlinear frequency shifts, which complicates the generation of solitons. Also, the spectral non-equidistance of neighboring mode pairs in microresonators is very substantial if compared to fiber-loop, open multi-mirror, and other types of low-repetition-rate and low-finesse resonators, where modelocking using χ⁽²⁾-effects have been demonstrated^29,30,31,32. Ref. ³³ provides an overview of theoretical studies of χ⁽²⁾ solitons in resonators between the 1990s and our days.

Lithium niobate (LN) is one of the favored χ⁽²⁾ materials to use in nano-fabrication for nonlinear and quantum optics applications³⁴. However, LN and the other χ⁽²⁾ materials possess appreciable χ⁽³⁾, i.e., Kerr, nonlinearity. In particular, the prior soliton demonstrations in the thin-film LiNbO₃ resonators^35,36 have been attributed to the Kerr effect. Also, recent experiments with comb generation in AlN microresonators have revealed the strong competition between χ⁽²⁾ and Kerr effects^37,38. The AlN microresonator used for the parametric down-conversion has allowed observation of bright solitons in the IR signal accompanied by the non-localized modelocked waveform in the near-IR pump³⁹. Thus, the challenge of observing the two-color bright-bright or dark-dark frequency-comb solitons in the χ⁽²⁾-mediated high-repetition-rate microresonator frequency conversion has so far remained unresolved.

Our present work demonstrates the dissipative two-color solitons and breathers in a quasi-phase-matched microresonator pumped for second-harmonic generation in the near-IR spectral range. The soliton regime is found to be typical for phase-mismatched resonators, while phase-matched ones reveal broader but incoherent spectra and higher-order harmonic generation. Positive phase-mismatching by less than one free spectral range induces tilting of the resonance line towards negative detunings, which is possible only via the dominant contribution of second-order nonlinearity.

Results

Here, we study the second-harmonic generation from the IR (1550 nm) to near-IR (780 nm) spectrum in the periodically polled thin-film LiNbO₃ microresonator. Our experimental setup is illustrated in Fig. 1. Our resonator radius is 70 μm, which provides a high repetition rate ≃290 GHz. We use the quasi-phase-matching grating to provide a large controllable positive phase mismatch, which interplays with the χ⁽²⁾ nonlinearity, making the resonance peak tilt towards negative detunings, see Fig. 2a.

Fig. 1: Experimental setup and resonator dispersion.

a Scanning-electron-microscopy image of the lithium niobate microresonator. The radius of the microresonator is 70 μm, corresponding to a 286 GHz repetition rate. b Photograph of the chip. c Measured (blue circles) and computed (blue and red lines) integrated dispersion, ω_μζ − ω_0ζ − D_1ζμ, vs. μ. Blue marks the infrared pump, ζ = a, and red marks the near-infrared second harmonic, ζ = b. d Measurement setup: EDFA erbium-doped fiber amplifier, WDM wavelength-division multiplexer, ESA electrical spectrum analyzer, OSA optical spectrum analyzer, FBG fiber Bragg grating.

Fig. 2: Nonlinear single-mode solutions and soliton families.

a Magenta lines show the single-mode (continuous-wave (CW)) solutions numerically computed for the on-chip powers ${{{{{{{\mathcal{W}}}}}}}}=60$ mW and 300 mW. Blue-yellow colors show the CW instability boundaries relative to the generation of the ±μ sideband pairs. μ values are indicated where possible. b A family of the soliton states (red is unstable and black is stable) and CW state (black is stable and magenta is unstable) vs. detuning: ε/2π = 95 GHz, pump wavelength is 1552 nm.

Below, we present measurements of the optical and radio-frequency (RF) spectra, which we interpret by the existence of bright-bright and dark-dark two-color soliton pairs and breathers. We demonstrate that the bright and dark solitons merge into a single family continuously on the variation of the system parameters. The merging becomes possible through angular periodicity and small ring sizes. The blurred difference between the bright and dark solitons manifests itself in the measured and modeled periodic expansion and shrinking of the solitons. Our experiment deals with the case when the pump experiences large anomalous GVD, and the second-harmonic is in the large normal GVD range. Despite this, the pump and second-harmonic solitons have the same type, e.g., if one is bright, the other is too, which appears to be the case not yet met in the resonator and modelocking contexts. The detailed numerical analysis guides our interpretations of the data.

The width of the resonator ridge is 1.8 μm, and the vertical dimensions are 410 and 590 nm to the air–LiNbO₃ and LiNbO₃–SiO₂ interfaces, respectively. A bus waveguide is specifically designed for simultaneous telecom and near-visible light coupling. Following the design principle elaborated in⁴⁰, the waveguide width, wrap-around angle, and resonator-waveguide coupling gap are optimized to be 1.8 μm, 60°, and 400 nm. The resonator spectrum near the pump, ζ = a, and second-harmonic, ζ = b, is approximated by ω_μζ = ω_0ζ + ∑_nD_nζμⁿ/n!. Here, ω_0a and ω_0b are the resonator frequencies with the mode numbers M and 2M + Q, respectively, where 2πR/Q is the poling period. Q equals 150 in the resonator sample used in the experiments described below. μ = 0, ± 1, ± 2, …  is the relative mode number. The phase mismatch is characterized by parameter ε³³,

$$\varepsilon=2{\omega }_{0a}-{\omega }_{0b},$$

(1)

which is determined by Q and temperature tuning. The resonator repetition rates are D_1a/2π = 286.24 GHz, D_1b/2π = 289.24 GHz, and make the 3 GHz difference. Second-order dispersion is large anomalous near the pump, D_2a/2π = 14 MHz, and large normal near the second harmonic, D_2b/2π = −18 MHz. Linewidths are rounded to κ_a/2π = 600 MHz and κ_b/2π = 1.2 GHz. The other parameters, as well as the coupled-mode equations, are described in Methods.

The microresonator chip is placed on a piezo-positioning stage with a standard laser temperature controller set at 130 °C, limiting the variations to less than 0.01 °C and providing ε/2π = 95 GHz for the pump around 1550 nm. The pump is coupled into the microresonator using an aspheric lens (numerical aperture = 0.6) and out of waveguide with a butt-coupled fiber. The pump is amplified in an erbium-doped fiber amplifier, whose output power is set at 27 dBm. The coupling loss for the 1550 and 780 nm light is estimated to be 7–8 dB/facet and 12–13 dB/facet, respectively. The output spectra are recorded using two grating-based optical spectrum analyzers covering 350–1750 and 1500–3400 nm. The generated dual-band frequency combs are spectrally separated using a wavelength-division multiplexer. The comb’s optical and RF spectra are characterized using the grating-based optical spectrum analyzer and electrical spectrum analyzer, respectively. The resolution bandwidth of 100 kHz is utilized for the RF noise measurement. Detuning, δ = ω_0a − ω_p, between the laser frequency, ω_p, and the resonance at ω_0a, is scanned from its negative (blue detuned) to the positive (red detuned) values.

The comb generation occurs when the pump frequency moves towards the resonance and makes the intra-resonator power exceed the modulational instability threshold⁴¹. It triggers the simultaneous growth of sidebands around the pump and its second harmonic, which further develops into the dual-band comb; see the experimental and numerical spectra in Fig. 3. To compute the regions of instabilities of the single-mode, i.e., continuous-wave (CW), operation relative to the generation of the ±μ pairs, we apply the approximation-free part of the formalism developed in ref. ⁴¹, see Fig. 2a. The parameter space in Fig. 2a is spanned by the pump detuning δ and intra-resonator pump power in the μ = 0 mode, ∣a₀∣². The CW states are shown with the magenta lines. When the CW crosses into an instability tongue, it becomes unstable relative to the respective ±μ mode pair. Two crossing points of the yellow (μ = 0 instability) and magenta lines limit the range of the CW bistability for a given on-chip laser power, ${{{{{{{\mathcal{W}}}}}}}}$.

Fig. 3: Experimental and numerical data for soliton states.

Panels a–i show four pairs of the low-noise soliton spectra experimentally measured with the increase of the pump detuning, δ. Panels e, j show the experimental RF spectra with the −90 dBm noise levels representative of the simultaneous soliton formation at 1550 and 780 nm. Black lines in e, j mark background electrical noise. Spectral envelopes shown with black lines in a–d, f–i are computed numerically and correspond to points (i)–(iv) in Fig. 2b. Panels k–t show numerical soliton spectra, and u–y show the respective pulse shapes at the points (i)–(v) marked in Fig. 2b and between the third and fourth columns here. The transition from the dark to bright solitons is clear from (u–y), which follows the collap4”n sed snaking trajectory in Fig. 2b. Left and right vertical axes in the u–y mark power for the 1550 and 780 nm pulses, respectively.

The negative direction of the tilt of the resonance curve, see Fig. 2, is determined by the χ⁽²⁾ effect and ε > 0. Hence, solitons and breathers reported below for negative detunings, δ < 0, are attributed to the χ⁽²⁾ interaction. Kerr nonlinearity is accounted for in all our simulations, see Methods, and plays only a corrective role since the dominant χ⁽³⁾ effect would cause the opposite, i.e., the positive tilt of the resonance well known from the theory and observations of Kerr solitons^1,3.

The power conversion efficiency of the IR (1550 nm) frequency comb is defined as ${{{{{{{{\mathcal{W}}}}}}}}}_{{{{{{{{\rm{IR}}}}}}}}}/{{{{{{{\mathcal{W}}}}}}}}$, where ${{{{{{{{\mathcal{W}}}}}}}}}_{{{{{{{{\rm{IR}}}}}}}}}$ is the total power in all IR comb lines excluding the central one. The measured conversion increases from ~5% to ~25% as the pump detuning increases, which could be further improved by optimizing the extraction efficiency in the IR. The improvement in conversion with growing δ is evident from the measured and simulated spectra, where the central peak first dominates over the IR comb, see Fig. 3f, and then blends with it, see Fig. 3g, h.

By solving the coupled-mode equations, we have found a family of the stationary modelocked pulses, i.e., solitons, associated with the observed spectra, see Figs. 2b and 3, and determined the corresponding soliton repetition rate ${\widetilde{D}}_{1}\,\ne \,{D}_{1}$. The soliton branch splits from the unstable high-power CW state, follows the snaking trajectory, and ends on the lower CW state. The snake line in Fig. 2b starts and ends at the points where the CW magenta line in Fig. 2a becomes unstable relative to the generation of the μ = ±1 sidebands. Pulse profiles near the upper and lower CW states correspond to the two-color dark–dark and bright-bright solitons, respectively. The pulse profiles in the IR and NIR are practically the same, while the IR power is around one order of magnitude smaller.

The snaking soliton line in Fig. 2b corresponds to the definition of the collapsed snaking used to describe a sequence of bifurcations of dark localized structures in the Kerr models with normal dispersion⁴², see also the earlier results in, e.g., ref. ⁴³. A feature of our resonator sample is that the pulse size is comparable with the ring circumference. Therefore, periodic boundary conditions make the dark soliton transform into the bright one after several turns of the snake, see Fig. 3u–y. The above-mentioned solitons with the high conversion into the NIR comb correspond to the nearly vertical region of the snake trajectory. The duty cycles of the respective pulses in Fig. 3u–y match the measured conversion efficiency. The low (−90 dBm) levels of noise in RF spectra are characteristic of the high degree of coherence typical for the dissipative multi-mode solitons, see Fig. 3e, j.

While doing the detuning scan and before entering the soliton regime, we first observed the characteristic three-peak spectra, see Fig. 4. The corresponding RF spectra are characterized by several well-defined peaks, see Fig. 4c–f and similar RF spectra of breathers in Kerr microresonators^44,45. Numerical analysis reveals that such regimes correspond to the two wavefronts moving in the ring resonator. These fronts eventually meet to create a pulse, which then starts expanding again, and the cycle repeats periodically; see numerical spectra and the space–time dynamics in Fig. 4g–i and Fig. 5, respectively.

Fig. 4: Experimental and numerical data for breathers.

a, b, d, e Experimentally measured infrared and NIR spectra found on the approach to the soliton range. The corresponding RF spectra in c, f reveal the onset of modelocking via the formation of the soliton breather. Optical and RF spectra corresponding to the numerically found breathers are shown in (g–i). δ/κ_a = − 1.64, which is left from the snake turns in Fig. 2b and soliton data in Fig. 3.

Fig. 5: Numerically computed breathers.

Space–time evolution of the two-color breather state, which spectra are shown in Fig. 4g–i. a The 1550 nm field intensity, i.e., ∣A∣² vs. t and θ, and b The 780 nm field, ∣B∣², see Eq. (2) in Methods.

Our experiments have also revealed the trade-off between the soliton regimes for the resonator samples with large ε and the generation of the high-bandwidth incoherent combs and higher order harmonics for the phase-matched resonators with ε close to zero achieved by tuning the temperature to T = 30 °C. The low-noise RF noise operation becomes inaccessible for ε close to zero. Experimentally recorded combs in the ε = 0 resonator feature the broad bandwidth incoherent spectra centered around 1560 nm (pump), 780 nm (2nd harmonic), 520 nm (3rd harmonic), and 390 nm (4th harmonic), which are plotted in blue, orange, green, and purple in Fig. 6. The measured on-chip pump-to-second-harmonic-comb conversion efficiencies is around 20.7%, which is mainly limited by the bus waveguide-microring coupling condition at both near-IR and near-visible wavelength bands. The inset shows a visible light emission from the resonator captured using a CCD camera.

Fig. 6: Broadband incoherent spectra in a phase-matched microresonator.

Measurements of the broadband incoherent frequency comb generation-spanning across four octaves in the phase-matched microresonator, ε = 0. The on-chip pump power is ${{{{{{{\mathcal{W}}}}}}}}=100$ mW.

Discussion

Our observations demonstrate two-colour dissipative solitons in the thin-film periodically polled LiNbO₃ microresonator with the ~300 GHz pulse repetition rate. A short resonator circumference limits the soliton number to one, makes possible the merging of bright and dark soliton families, and plays a role in the front-motion instability triggering breather dynamics. Solitons, breathers, and frequency comb generation reported here happen for the negative tilt of the resonance line, which is only possible via the dominant contribution of second-order nonlinearity. Future research directions include, e.g., implementing resonance tilting towards positive detunings, engineering different combinations of dispersion signs, and generation of shorter solitons and χ⁽²⁾ soliton crystals.

Methods

Device fabrication

Devices are fabricated from a commercial LN on an insulator wafer (supplied by NANOLN), in which a 590 nm thick Z-cut LN layer sits on top of 2 μm silicon dioxide on a silicon handle. The pattern is first defined using electron beam lithography (EBL) with a negative FOx-16 resist and subsequently transferred onto the LN layer using an optimized inductively coupled plasma reactive ion etching process with Ar⁺ plasma. A thin layer of hafnium oxide is deposited on top of the fabricated photonic device using the atomic layer deposition technique, which serves as a protection layer from metal contamination induced during the poling process and also aids in confining the electric field for high-fidelity poling as a high-k material. The radial nickel electrodes are patterned on top of the LN microring concentrically via the EBL with alignment and following a bi-layer lift-off process. An optimized poling sequence was then applied to create the desired polling pattern. Afterward, the electrodes and oxide interface were sequentially removed by wet etching. Finally, the chip is cleaved to expose the waveguide facets for fiber-to-chip coupling.

Numerical simulation

Multimode intra-resonator pump field (1552 nm, TM polarized) and it is second-harmonic (TM polarized) are expressed via their mode expansions as³³

$$ A{e}^{iM\vartheta -i{\omega }_{p}t}+c.c.,\,\,A=\mathop{\sum}\limits_{\mu }{a}_{\mu }(t){e}^{i\mu \theta },\\ B{e}^{i(2M+Q)\vartheta -i2{\omega }_{p}t}+c.c.,\,\,B=\mathop{\sum}\limits_{\mu }{b}_{\mu }(t){e}^{i\mu \theta },\\ \theta=\vartheta -{\widetilde{D}}_{1}t.$$

(2)

Here, ϑ = (0, 2π] is the angular coordinate in the laboratory frame, θ is the coordinate in the reference frame rotating with the rate ${\widetilde{D}}_{1}$, Q = 150 is the number of the poling periods, M = 515 and μ = 0, ± 1, ± 2, …  is the relative mode number. a_μ, b_μ are the amplitudes of the pump and second-harmonic modes.

The resonator spectrum around the pump and second harmonic is approximated as

$$\begin{array}{l}{\omega }_{\mu \zeta }={\omega }_{0\zeta }+\mu {D}_{1\zeta }+\frac{1}{2!}{\mu }^{2}{D}_{2\zeta }+\frac{1}{3!}{\mu }^{3}{D}_{3\zeta }+\frac{1}{4!}{\mu }^{4}{D}_{4\zeta },\\ \zeta=a,b,\end{array}$$

(3)

where, D_1ζ are the linear repetition rates, D_2ζ are the second order dispersions, and D_3ζ, D_4ζ are the third and fourth order ones. D_1b − D_1a is the walk-off parameter, i.e., the repetition rate difference. The values of the dispersion coefficients are specified in Fig. 1. The pump laser, ω_p, is tuned around the 1552.3 nm wavelength and targets the ω_0a resonance. The pump detuning is defined as

$$\delta={\omega }_{0a}-{\omega }_{p}.$$

(4)

Coupled-mode equations governing the evolution of a_μ(t), b_μ(t) include both χ⁽²⁾ and χ⁽³⁾ nonlinearities^33,37. The equations have been derived under the assumption that the e^iQϑ Fourier component of the quasi-phase-matching grating, χ⁽²⁾G(ϑ) = χ⁽²⁾G(ϑ + 2π/Q), provides the required phase-matching³³,

$$i{\partial }_{t}{a}_{\mu }= {\delta }_{\mu a}{a}_{\mu }-\frac{i{\kappa }_{a}}{2}\left({a}_{\mu }-{\widehat{\delta }}_{\mu,0}{{{{{{{\mathcal{H}}}}}}}}\right)\\ -{\gamma }_{2a}\mathop{\sum}\limits_{{\mu }_{1}{\mu }_{2}}{\widehat{\delta }}_{\mu,{\mu }_{1}-{\mu }_{2}}{b}_{{\mu }_{1}}{a}_{{\mu }_{2}}^{*}\\ -{\gamma }_{3a}\mathop{\sum}\limits_{{\mu }_{1}{\mu }_{2}{\mu }_{3}}{\widehat{\delta }}_{\mu,{\mu }_{1}+{\mu }_{2}-{\mu }_{3}}{a}_{{\mu }_{1}}{a}_{{\mu }_{2}}{a}_{{\mu }_{3}}^{*}\\ -2{\gamma }_{3a}\mathop{\sum}\limits_{{\mu }_{1}{\mu }_{2}{\mu }_{3}}{\widehat{\delta }}_{\mu,{\mu }_{1}+{\mu }_{2}-{\mu }_{3}}{a}_{{\mu }_{1}}{b}_{{\mu }_{2}}{b}_{{\mu }_{3}}^{*},\\ i{\partial }_{t}{b}_{\mu }= {\delta }_{\mu b}{b}_{\mu }-\frac{i{\kappa }_{b}}{2}{b}_{\mu }\\ -{\gamma }_{2b}\mathop{\sum}\limits_{{\mu }_{1}{\mu }_{2}}{\widehat{\delta }}_{\mu,{\mu }_{1}+{\mu }_{2}}{a}_{{\mu }_{1}}{a}_{{\mu }_{2}}\\ -{\gamma }_{3b}\mathop{\sum}\limits_{{\mu }_{1}{\mu }_{2}{\mu }_{3}}{\widehat{\delta }}_{\mu,{\mu }_{1}+{\mu }_{2}-{\mu }_{3}}{b}_{{\mu }_{1}}{b}_{{\mu }_{2}}{b}_{{\mu }_{3}}^{*}\\ -2{\gamma }_{3b}\mathop{\sum}\limits_{{\mu }_{1}{\mu }_{2}{\mu }_{3}}{\widehat{\delta }}_{\mu,{\mu }_{1}+{\mu }_{2}-{\mu }_{3}}{b}_{{\mu }_{1}}{a}_{{\mu }_{2}}{a}_{{\mu }_{3}}^{*}.$$

(5)

Here, ${\widehat{\delta }}_{\mu,{\mu }^{{\prime} }}=1$ for $\mu={\mu }^{{\prime} }$ and is zero otherwise. ${{{{{{{\mathcal{H}}}}}}}}$ is the pump parameter, ${{{{{{{{\mathcal{H}}}}}}}}}^{2}=\eta {{{{{{{\mathcal{F}}}}}}}}{{{{{{{\mathcal{W}}}}}}}}/2\pi$, where ${{{{{{{\mathcal{W}}}}}}}}$ is the laser power, and ${{{{{{{\mathcal{F}}}}}}}}={D}_{1a}/{\kappa }_{a}$ is finesse³³. η is the coupling coefficient, which was used as the fitting parameter, η = 0.33333. δ_μζ are the modal detuning parameters in the rotating reference frame,

$${\delta }_{\mu a}= ({\omega }_{\mu a}-{\omega }_{p})-\mu {\widetilde{D}}_{1},\\ {\delta }_{\mu b}= ({\omega }_{\mu b}-2{\omega }_{p})-\mu {\widetilde{D}}_{1},$$

(6)

where δ_0a = δ, δ_0b = 2δ − ε and ε is the phase mismatch parameter³³,

$$\varepsilon=2{\omega }_{0a}-{\omega }_{0b}=2\frac{cM}{R{n}_{M}}-\frac{c(2M+Q)}{R{n}_{2M+Q}}.$$

(7)

Here, c is the vacuum speed of light, R is the resonator radius, n_M is the effective refractive index experienced by the resonator mode with the number M. The value of ε can be controlled by the temperature and pump wavelength. T = 130^oC yields ε/2π ≈ 95GHz and was set to get the soliton and breather generation in Figs. 2–4. T = 30 °C gave ε ≈ 0 and was used to generate the incoherent multi-octave spectra in Fig. 5.

γ_2ζ and γ_3ζ are parameters specifying the strength of the second and third-order nonlinear effects³³. Using a simplifying assumption that the effective mode area, S, does not disperse, we estimate γ_2ζ as

$${\gamma }_{2\zeta }=\frac{d{\omega }_{0\zeta }q}{3{n}_{0}^{2}},\,q=\sqrt{\frac{2{{{{{{{{\mathcal{Z}}}}}}}}}_{vac}}{S{n}_{0}}},\,\zeta=a,b,$$

(8)

where n₀ = 2.2 is the linear refractive index, ω_0a/2π = 193 THz, ω_0b/2π = 386 THz, ${{{{{{{{\mathcal{Z}}}}}}}}}_{vac}=1/{\epsilon }_{vac}c=377$ V²/W is the free space impedance, and the averaged effective area S ≈ 1.5 μm². These values yield q ≈ 15 × 10⁶ W^−1/2 V/m d ~χ⁽²⁾ is the relevant element of the reduced χ⁽²⁾ tensor, d ≈ 20 pm/V and γ_2a/2π ≈ 4GHz/$\sqrt{{{{{{{{\rm{W}}}}}}}}},\,{\gamma }_{2b}/2\pi \,\approx \,8$GHz/$\sqrt{{{{{{{{\rm{W}}}}}}}}}$. Kerr parameters γ_3ζ are estimated using the results derived in ref. ⁴⁶

$${\gamma }_{3\zeta }=\frac{{\omega }_{0\zeta }{n}_{2}}{2S{n}_{0}},$$

(9)

where n₂ = 9 × 10⁻²⁰ m²/W² (χ⁽³⁾ = 1.6 × 10⁻²¹ m²/V²) is the Kerr coefficient of LiNbO₃ and γ_3a/2π ≈ 2.5 MHz/W, γ_3b/2π ≈ 5 MHz/W.

According to the estimates based on the comparison of the nonlinear resonance shifts induced by the χ⁽²⁾ and χ⁽³⁾ terms⁴¹, the latter is expected to play a notable role in ∣ε∣ becoming close to and exceeding ε_cr,

$${\varepsilon }_{{{{{{{{\rm{cr}}}}}}}}}={\gamma }_{2a}{\gamma }_{2b}/{\gamma }_{3a}\,\approx \,2\pi \times 2{{{{{{{\rm{THz}}}}}}}},$$

(10)

which corresponds to the mismatch by about six modes in a resonator with a 300 GHz repetition rate and is much larger than ε/2π ≃ 0.1 THz in our resonator.

Typical time-dependent simulations of Eq. (5) were performed using the fourth-order Runge–Kutta method applying ${\widetilde{D}}_{1}={D}_{1a}$. Stationary soliton profiles were found using the Newton method after ∂_t was set to zero and ${\widetilde{D}}_{1}$ assumed as one of the unknowns. The value of ${\widetilde{D}}_{1}$ after calculations was close, but not equal, to D_1a. A typical number of modes around ω_p and 2ω_p used in the modeling was 256, i.e., μ = − 127, …, 0, …, 128.

To differentiate between stable and unstable solitons, we have analyzed the linear stability of the soliton family. We perturbed the time-independent soliton amplitudes, ${\hat{a}}_{\mu }$, ${\hat{b}}_{\mu }$ (${\partial }_{t}{\hat{a}}_{\mu }={\partial }_{t}{\hat{b}}_{\mu }=0$) with small perturbations, ${\varepsilon }_{a\mu }(t)={x}_{a\mu }(t)+{y}_{a\mu }^{*}(t)$ and ${\varepsilon }_{b\mu }(t)={x}_{b\mu }(t)+{y}_{b\mu }^{*}(t)$, i.e.,

$${a}_{\mu }(t) ={\hat{a}}_{\mu }+{x}_{a\mu }(t)+{y}_{a\mu }^{*}(t),\\ {b}_{\mu }(t) ={\hat{b}}_{\mu }+{x}_{b\mu }(t)+{y}_{b\mu }^{*}(t).$$

(11)

and then linearized Eq. (5). By substituting $\{{x}_{a\mu }(t),{y}_{a\mu }(t),{x}_{b\mu }(t),{y}_{b\mu }(t)\}=\{{X}_{a\mu },{Y}_{a\mu },{X}_{b\mu },{Y}_{b\mu }\}{e}^{\lambda t}$, we have reduced the linearized differential equations to the algebraic (4 × 256) × (4 × 256) eigenvalue problem, which was solved numerically⁴¹. The soliton is stable if all Reλ < 0. The stability of the CW state, a_μ≠0 = 0, b_μ≠0 = 0, was found from the simpler four-by-four eigenvalue problem⁴¹. In this case, each eigenvalue λ is attributed to a particular pair of ±μ sidebands producing its own instability boundary, which are all plotted in Fig. 2b.