Parametrically driven pure-Kerr temporal solitons in a chip-integrated microcavity

The injection of monochromatic continuous wave (CW) laser light into dispersive optical resonators with purely Kerr-type $\chi^{(3)}$ nonlinearity can lead to the generation of localized structures known as dissipative Kerr cavity solitons (CSs) leo_temporal_2010 ; kippenberg_dissipative_2018 . These CSs correspond to ultrashort pulses of light that can persist within the resonator [Fig. 1(a)], indefinitely maintaining constant shape and energy wabnitz_suppression_1993 . While first observed in macroscopic optical fiber ring resonators leo_temporal_2010 , CSs have attracted particular attention in the context of monolithic Kerr microcavities kippenberg_dissipative_2018 , where they underpin the generation of coherent and broadband optical frequency combs herr_temporal_2014 ; brasch_photonic_2016 ; pasquazi_micro-combs_2018 . By offering a route to coherent frequency comb generation in chip-integrated, foundry-ready platforms, CSs have enabled ground breaking advances in applications including telecommunications marin-palomo_microresonator-based_2017 ; corcoran_ultra-dense_2020 , artificial intelligence xu_11_2021 ; feldmann_parallel_2021 , astronomy suh_searching_2019 ; obrzud_microphotonic_2019 , frequency synthesis spencer_optical-frequency_2018 , microwave generation lucas_ultralow-noise_2020 ; kwon_ultrastable_2022 , and distance measurements suh_soliton_2018 ; riemensberger_massively_2020 .

The conventional CSs that manifest themselves in resonators with pure Kerr nonlinearity sit atop a CW background, and they gain their energy through four-wave-mixing (FWM) interactions with that background leo_temporal_2010 . In the frequency domain, the solitons are (to first order) centred around the frequency of the external CW laser that drives the resonator [Fig. 1(a)]. They are (barring some special exceptions nielsen_coexistence_2019 ; anderson_coexistence_2017 ; hansson_frequency_2015 ; xu_spontaneous_2021 ; lucas_spatial_2018 ) unique attracting states: except for trivial time translations, all the CSs that exist for given system parameters are identical. These features can be disadvantageous or altogether prohibitive for selected applications: noise on the external CW laser can degrade the coherence of nearby comb lines, removal of the CW background may require careful spectral filtering, whilst applications that require coexistence of distinguishable binary elements takesue_10_2016 ; okawachi_quantum_2016 ; okawachi_dynamic_2021 ; inagaki_Large-scale_2016 ; mohseni_ising_2022 are fundamentally beyond reach. Interestingly, recent experiments reveal that qualitatively different types of CSs can exist in resonators that display a quadratic $\chi^{(2)}$ in addition to a cubic $\chi^{(3)}$ nonlinearity [Fig. 1(b)]; in particular, degenerate optical parametric oscillators driven at $2\omega_{0}$ can support CSs at $\omega_{0}$ englebert_parametrically_2021 ; bruch_pockels_2021 . In this configuration, the solitons are parametrically driven through the quadratic down-conversion of the externally-injected field, which endows them with fundamental differences compared to the conventional CSs emerging in monochromatically-driven, pure-Kerr resonators. Specifically, parametrically driven cavity solitons (PDCSs) are spectrally separated from the driving frequency (e.g. $\omega_{0}$ versus $2\omega_{0}$), and they come in two binary forms with opposite phase. These traits render PDCSs of interest for an altogether new range of applications.

Optical PDCSs have so far been generated only via the quadratic $\chi^{(2)}$ nonlinearity, which is not intrinsically available in integrated (foundry-ready) resonator platforms, such as silicon thomson_roadmap_2016 or silicon nitride liu_high-yield_2021 ; luSHG2020 ; NitissNat.Photon.2022 . However, it is well-known that phase-sensitive amplification analogous to $\chi^{(2)}$ parametric down-conversion can also be realised in pure Kerr resonators when driven with two lasers with different carrier frequencies mecozzi_long-term_1994 ; agrawal_nonlinear_nodate ; radic_two-pump_2003 ; okawachi_dual-pumped_2015 ; Andrekson_fiber-based_2020 , allowing e.g. for novel random number generators takesue_10_2016 ; okawachi_quantum_2016 ; okawachi_dynamic_2021 and coherent optical Ising machines inagaki_Large-scale_2016 ; mohseni_ising_2022 . A natural question that arises is: is it possible to generate PDCSs in foundry-ready, pure-Kerr resonators with bichromatic driving? Whilst a related question has been theoretically explored in the context of diffractive Kerr-only resonators de_valcarcel_phase-bistable_2013 , the presence of dispersion substantially changes the physics of the problem. The impact of bichromatic driving in the dynamics of conventional Kerr CSs has also been considered hansson_bichromatically_2014 ; ceoldo_multiple_2016 ; zhang_spectral_2020 ; moille_ultra-broadband_2021 ; qureshi_soliton_2021 ; taheri_all-optical_2022 , but the possibility of using the scheme to generate temporal PDCSs remains unexplored.

Here, we theoretically predict and experimentally demonstrate that a dispersive resonator with pure Kerr nonlinearity can support PDCSs in the presence of bichromatic driving [Fig. 1(c)]. We reveal that, under appropriate conditions, a signal field with carrier frequency in between two spectrally-separated driving fields obeys the damped, parametrically driven nonlinear Schrödinger equation (PDNLSE) that admits PDCS solutions, and we unveil the system requirements for the practical excitation of such solutions. Our experiments are performed in a 23 $\mu$m-radius, chip-integrated silicon nitride microring resonator whose dispersion is judiciously engineered to facilitate PDCS generation at 253 THz (1185 nm) when bichromatically pumping at 314 THz (955 nm) and 192 THz (1560 nm). We observe PDCS frequency comb spectra that are in good agreement with numerical simulations, as well as clear signatures of the anticipated $\mathbb{Z}_{2}$ symmetry, i.e., coexistence of two PDCSs with opposite phase. By revealing a fundamentally new pathway for the generation of coherent PDCS frequency combs far from any pump frequency, in a platform that has direct compatibility with foundry-ready fabrication, our work paves the way for integrated, low-noise frequency comb generation in new spectral regions, as well as photonic integration of applications requiring combs with a binary degree of freedom.

We first summarise the main points that lead to the prediction of PDCSs in bichromatically-driven Kerr resonators [for full details, see Methods]. To this end, we consider a resonator made out of a dispersive, $\chi^{(3)}$ nonlinear waveguide that is driven with two coherent CW fields with angular frequencies $\omega_{\pm}$ [see Fig. 1(c)]. The dispersion of the resonator is described by the integrated dispersion pasquazi_micro-combs_2018 at the cavity resonance $\omega^{\prime}_{0}$ (apostrophes highlight resonance frequencies throughout the article) closest to the frequency $\omega_{0}=(\omega_{+}+\omega_{-})/2$:

$$D_{\mathrm{int}}(\mu)=\omega^{\prime}_{\mu}-\omega^{\prime}_{0}-\mu D_{1}=\sum_{k\geq 2}\frac{D_{k}}{k!}\mu^{k}.$$

(1)

Here, $\mu$ is a relative mode number with respect to the resonance $\omega_{0}^{\prime}$ and $D_{1}/(2\pi)$ is the cavity free-spectral range (FSR) at $\omega^{\prime}_{0}$. The terms $D_{k}$ with $k>1$ account for deviations of the resonance frequencies $\omega^{\prime}_{\mu}$ from an equidistant grid defined by $\omega^{\prime}_{0}+\mu D_{1}$.

Under particular conditions [see Methods], the evolution of the slowly-varying electric field envelope centred at $\omega_{0}$ can be shown to be (approximately) governed by the PDNLSE, with the parametric driving ensuing from non-degenerate FWM driven by the intracavity fields at the pump frequencies [$\omega_{+}+\omega_{-}\rightarrow\omega_{\mu}+\omega_{-\mu}$, see Fig. 1(d)]. (Note: in stark contrast to standard Kerr CSs, for which only one comb line is externally driven, all of the components of a PDCS frequency comb are separately driven via non-degenerate FWM.) Because the PDNLSE is well-known to admit PDCS solutions englebert_parametrically_2021 ; miles_parametrically_1984 ; barashenkov_stability_1991 , it follows that the system may support such solitons with a carrier frequency $\omega_{0}$ in between the two driving frequencies, provided however that the system parameters – particularly resonator dispersion – are conducive for soliton existence.

The resonator dispersion must meet three key conditions for PDCS excitation to be viable [Methods]. First, for solitons to exist, the dispersion around the degenerate FWM frequency $\omega_{0}$ must be anomalous, i.e., $D_{2}>0$ in Eq. (1). Second, the effective detuning [see Methods] between the degenerate FWM frequency ($\omega_{0}$) and the closest cavity resonance ($\omega^{\prime}_{0}$) must be within the range of soliton existence, essentially requiring that the degenerate FWM process $\omega_{+}+\omega_{-}\rightarrow 2\omega_{0}$ (approximately) satisfies linear phase-matching [Fig. 1(e)]. This second condition can be written as $\delta\omega=(\omega^{\prime}_{+}+\omega^{\prime}_{-})/2-\omega_{0}^{\prime}=[D_{\mathrm{int}}(p)+D_{\mathrm{int}}(-p)]/2\approx 0$, where $\pm p$ correspond to the modes excited by the driving lasers at $\omega_{\pm}$. Given that $D_{2}>0$, this requires at least one higher-even-order dispersion coefficient (e.g. $D_{4}$) to be negative. Third, the intracavity field amplitudes at the driving frequencies, $|E_{\pm}|$, must remain (approximately) homogeneous and stationary to ensure a constant parametric driving strength for the PDCS field $E_{0}$ centred at $\omega_{0}$ [Fig. 1(f)]. This final condition can be met by ensuring dispersion at the driving frequencies is (i) normal (or driving amplitudes small), such that the corresponding intracavity fields do not undergo pattern forming (modulation) instabilities coen_universal_2013 , and (ii) such that the temporal walk-off between the driving frequencies $\omega_{\pm}$ and the signal frequency $\omega_{0}$ is sufficiently large so as to mitigate pump depletion in the vicinity of the soliton that would otherwise break the homogeneity of the fields at $\omega_{\pm}$ [Fig. 1(f)]. As will be demonstrated below, all of these conditions can be met through judicious dispersion engineering that is within the reach of contemporary microphotonic fabrication.

Figure 2(b) shows the evolution of the numerically simulated intracavity intensity profile with an initial condition consisting of two hyperbolic secant pulses with opposite phases. As can be seen, after a short transient, the field reaches a steady-state that is indicative of two pulses circulating around the resonator. The pulses sit atop a rapidly oscillating background that is due to the beating between the quasi-homogeneous fields at the pump frequencies [Fig. 2(c)]. Correspondingly, the spectrum of the simulation output [Fig. 2(d)] shows clearly the presence of a hyperbolic secant-shaped feature that sits in between the strong quasi-monochromatic components at the pump frequencies. In accordance with PDCS theory [see Methods], there is no significant CW peak at the parametric signal frequency $\omega_{0}$ at which the solitons are spectrally centred. To highlight the phase disparity of the steady-state pulses, we apply a numerical filter to remove the quasi-monochromatic intracavity components around the pump frequencies, and plot in Fig. 2(e) the real part of the complex intracavity electric field envelope. The simulation results in Fig. 2(e) are compared against the real parts of the exact, analytical PDCS solutions [Methods], and we clearly observe excellent agreement.

The results in Fig. 2(a)–(e) corroborate the fundamental viability of our scheme. However, they were obtained assuming a completely symmetric dispersion profile with no odd-order terms, which may be difficult to realise even with state-of-the-art microphotonic fabrication (including the resonators considered in our experiments). We find, however, that PDCSs can exist even in the presence of odd-order-dispersion, albeit in a perturbed form. This point is highlighted in Figs. 2(f)–(h), which show results from simulations with all parameters as in Fig. 2(a)–(e) except an additional non-zero third-order dispersion term $D_{3}=-2\pi\times 58\leavevmode\nobreak\ \mathrm{Hz}$. As for conventional (externally-driven) Kerr CSs coen_modeling_2013 ; jang_observation_2014 ; milian_soliton_2014 ; brasch_photonic_2016 , we find that third-order dispersion causes the solitons to emit dispersive radiation at a spectral position determined by the phase-matching condition $D_{\mathrm{int}}(\mu_{\mathrm{DW}})\approx(\omega_{0}-\omega_{0}^{\prime})$ [Fig. 2(g)]. This emission results in the solitons experiencing constant drift in the temporal domain, and endows them with oscillatory tails [Fig. 2(h)]. Yet, as can clearly be seen, the PDCSs continue to exist in two distinct forms with near-opposite phase. It is worth noting that, for the parameters considered in Fig. 2(f)–(h), the low-frequency driving field experiences anomalous group-velocity dispersion; however, the intracavity intensity at that frequency is below the modulation instability threshold coen_universal_2013 , thus allowing the corresponding field to remain quasi-homogeneous (the modulation on the total intensity profile arises solely from the linear beating between the different fields).

The orange curve in Fig. 3(b) depicts an estimate of the resonator’s integrated dispersion around a cavity mode at 253 THz, obtained through a combination of finite-element-modelling and fitting to our experimental observations [see Methods]. This data is consistent with experimentally measured resonance frequencies (blue circles), yet we caution that our inability to probe the resonances around 253 THz prevents unequivocal evaluation of the dispersion at that frequency. The estimated dispersion can be seen to be such that the requisite phase-matching for generating a PDCS at 253 THz ($\delta\omega\approx 0$) can be satisfied, provided that the pump lasers are configured to drive cavity modes at 314 THz and 192 THz [Fig. 3(c)].

In our experiments, we set the on-chip driving power for both driving fields to be about 150 mW and tune the high-frequency pump to the cavity mode at 314 THz. We then progressively tune the low-frequency pump to the cavity mode at 192 THz (from blue to red), maintaining the high-frequency pump at a fixed frequency. As the low-frequency pump tunes into resonance, we initially observe non-degenerate parametric oscillation characterised by the generation of two CW components symmetrically detuned about 253 THz. These CW components progressively shift closer to each other as the pump tunes into the resonance, concomitant with the formation of a frequency comb around the degenerate FWM frequency $\omega_{0}$ [see Fig. 3(c)]. To characterise the comb noise, we performed a heterodyne beat measurement using a helper laser at 230 THz within the vicinity of a single comb line. Initially, no beat note is observed, which is characteristic of an unstable, non-solitonic state within the resonator. Remarkably, as the 192 THz driving field is tuned further into resonance, we observe that the parametric signals reach degeneracy, concomitant with the emergence of a broadband comb state with smooth spectral envelope [Fig. 3(e)] and a heterodyne beat note (comparable with the helper laser linewidth of 250 kHz) that is considerably narrower than the 300 MHz microcavity linewidth [Fig. 3(f) and Supplementary Figure 1].

The emergence of the smooth comb state [Fig. 3(e)] is associated with an abrupt drop in the photodetector signal recorded around 253 THz, giving rise to a noticeable step-like feature [Fig. 3(f)]. Similar steps are well-known signatures of conventional CSs in monochromatically-driven Kerr resonators herr_temporal_2014 . Moreover, as shown in Fig. 3(e), the smooth spectral envelope observed in the step-region is in very good agreement with the spectrum of a 24 fs (full-width at half-maximum) PDCS derived from numerical modelling that use estimated experimental parameters [see Supplementary Figure 2]. The simulations faithfully reproduce the main features of the experimentally observed spectrum, including a strong dispersive wave peak at about 210 THz. We note that the prominent dip at about 275 THz arises due to the frequency-dependence of the pulley coupler moille_broadband_2019 , which was taken into account ad hoc when estimating the spectrum of the out-coupled PDCS [shown as blue curve in Fig. 3 – see also Methods and Supplementary Figure 2].

It is interesting to note that, in addition to the frequency comb around the degenerate FWM frequency 253 THz, frequency combs arise also around both of the pump frequencies. These combs originate from FWM interactions between the pump fields and the comb lines around 253 THz, in a manner similar to spectral extension zhang_spectral_2020 ; moille_ultra-broadband_2021 ; qureshi_soliton_2021 and two-dimensional frequency comb MoillearXiv2023 schemes studied in the context of conventional Kerr CSs. The combs around the pump frequencies share the line spacing with the comb around 253 THz, but there is a constant offset between the pump and PDCS combs. In our experiments, this comb offset is directly observable in the optical spectrum [inset of Fig. 3(e)] and found to be about 50 GHz $\pm$ 2 GHz (uncertainty defined by the optical spectrum analyzer resolution), which is in good agreement with the value of 49 GHz predicted by our modelling [see Methods]. All in all, given the considerable uncertainties in key experimental parameters (particularly dispersion and detunings), we find the level of agreement between the simulations and experiments remarkable.

The results shown in Fig. 3 are strongly indicative of PDCS generation in our experiments. Further confirmation is provided by observations of low-noise combs with complex spectral structures that afford a straightforward interpretation in terms of multi-PDCS states [Fig. 4]. Specifically, whilst a single PDCS circulating in the resonator is expected to yield a smooth spectral envelope, the presence of two (or more) PDCSs results in a spectral interference pattern whose details depend upon the soliton’s relative temporal delay and – importantly – phase.

Figures 4(a) and (b) show selected examples of multi-soliton comb spectra measured in our experiments. Also shown as solid curves are spectral envelopes corresponding to fields with two linearly superposed, temporally delayed PDCSs [Figs. 4(c) and (d) and Methods]. We draw particular attention to the fact that, in the measured data shown in Fig. 4(b), the comb component at the degenerate FWM frequency $\omega_{0}$ is suppressed by about 40 dB compared to neighbouring lines, which is in stark contrast with results in Fig. 4(a), in which the degenerate FWM component is dominant. This suppression is indicative of a relative phase shift of $0.992\pi$ between the two solitons [Fig. 4(d)] – a clear signature of PDCSs.

We have shown theoretically, numerically, and experimentally that dispersive resonators with a purely Kerr-type $\chi^{(3)}$ nonlinearity can support parametrically driven cavity solitons under conditions of bichromatic driving. Our theoretical analysis has revealed the salient conditions that the system dispersion must meet to allow for PDCS persistence, with approximation-free numerical simulations confirming the fundamental viability of the scheme. Experimentally, we realise suitable dispersion conditions in a chip-integrated silicon nitride microresonator, observing low-noise frequency comb states that evidence PDCS generation. Significantly, our measurements show spectral interference patterns that indicate the co-existence of two localized structures with opposite phase – a defining feature of PDCSs.

Our work fundamentally predicts and demonstrates that dispersive Kerr resonators can support a new type of dissipative structure – the PDCS – in addition to conventional Kerr CSs. We envisage that studying the rich nonlinear dynamics leo_dynamics_2013 ; yu_breather_2017 ; lucas_breathing_2017 , interactions jang_ultraweak_2013 ; cole_soliton_2017 ; wang_universal_2017 , and characteristics (including quantum chembo_quantum_2016 ; bao_quantum_2021 ; yang2021squeezed ; guidry_quantum_2022 ) of pure-Kerr PDCSs will draw substantial future research interest, echoing the extensive exploration of conventional Kerr CS dynamics over the past decade leo_dynamics_2013 ; yu_breather_2017 ; lucas_breathing_2017 ; jang_ultraweak_2013 ; cole_soliton_2017 ; wang_universal_2017 ; chembo_quantum_2016 ; bao_quantum_2021 ; yang2021squeezed ; guidry_quantum_2022 . In this context, to the best of our knowledge, the results reported in our work represent the first prediction and observation of dispersive wave emission by PDCSs in any physical system.

From a practical vantage, our scheme offers a route to generate PDCS frequency combs in foundry-ready, chip-integrated platforms with characteristics that are fundamentally different from those associated with conventional Kerr CSs. For example, forming between the two input frequencies, PDCSs could permit comb generation at spectral regions where direct pump lasers may not be available. Moreover, the lack of a dominating CW component at the PDCS carrier frequency alleviates the need for careful spectral shaping, and could result in fundamental advantages to noise characteristics. We emphasise that PDCSs are underpinned by phase-sensitive amplification mecozzi_long-term_1994 , which can theoretically offer a sub-quantum-limited (squeezed) noise figure McKinstrie_phase-sensitive_2004 ; slavik_all-optical_2010 ; tong_towards_2011 ; Zhang2008 ; Ye2021 . Finally, the fact that PDCSs come in two forms with opposite phase opens the doors to a new range of applications that require a binary degree of freedom, including all-optical random number generation and realisations of coherent optical Ising machines. Whilst the potential of PDCSs for such applications has been noted earlier englebert_parametrically_2021 , our work provides for the first time a route for chip-integrated realisations with potential to CMOS-compatible mass manufacturing.

We present here a heuristic derivation of the Ikeda-like map used in our numerical simulations [Eq. (5) of the main manuscript]. Including the rapid temporal oscillations at $\omega_{0}=(\omega_{+}+\omega_{-})/2$, where $\omega_{\pm}$ are the input frequencies, the intracavity electric field during the $m^{\mathrm{th}}$ cavity transit is written as $E^{(m)}(z,\tau)\exp[-i\omega_{0}T]$, where $\tau$ is time in a co-moving reference frame defined as $\tau=T-z/v_{\mathrm{g}}$ with $T$ absolute laboratory time, $z$ the coordinate along the waveguide that forms the resonator, $v_{\mathrm{g}}$ the group velocity of light at $\omega_{0}$, and $E^{(m)}(z,\tau)$ is the slowly-varying electric field envelope that follows the generalized nonlinear Schrödinger equation [Eq. (2) of the main manuscipt]. The boundary equation for the full electric field can then be written as

	$\displaystyle E^{(m+1)}(0,\tau)e^{-i\omega_{0}T}$	$\displaystyle=\sqrt{1-2\alpha}E^{(m)}(L,\tau)e^{-i\delta_{0}-i\omega_{0}T}$
		$\displaystyle+\sqrt{\theta_{+}}E_{\mathrm{in,+}}e^{-i\omega_{+}T}$
		$\displaystyle+\sqrt{\theta_{-}}E_{\mathrm{in,-}}e^{-i\omega_{-}T},$		(16)

where $\delta_{0}=2\pi k-\beta(\omega_{0})L$ is the linear phase detuning of the reference frequency $\omega_{0}$ from the closest cavity resonance (with order $k$), and $\omega_{\pm}$ are the frequencies of the pump fields. Multiplying each side with $\exp[i\omega_{0}T]$ and replacing $T=\tau+mL/v_{\mathrm{g}}=\tau+mt_{\mathrm{R}}$ where $t_{\mathrm{R}}$ is the round trip time, yields

	$\displaystyle E^{(m+1)}(0,\tau)$	$\displaystyle=\sqrt{1-2\alpha}E^{(m)}(L,\tau)e^{-i\delta_{0}}$
		$\displaystyle+\sqrt{\theta_{+}}E_{\mathrm{in,+}}e^{-i\Omega_{\mathrm{p}}\tau+i(\omega_{0}-\omega_{+})mt_{\mathrm{R}}}$
		$\displaystyle+\sqrt{\theta_{-}}E_{\mathrm{in,-}}e^{i\Omega_{\mathrm{p}}\tau+i(\omega_{0}-\omega_{-})mt_{\mathrm{R}}}.$		(17)

Note that, in the above formulation, the co-moving time variable $\tau$ should be understood as the “fast time” that describes the envelope of the intracavity electric field over a single round trip, i.e., the distribution of the envelope within the resonator pasquazi_micro-combs_2018 . As such, the $\tau$ variable spans a single round trip time of the resonator, and the intracavity envelope must obey periodic boundaries within that range. These conditions stipulate that the frequency variable $\Omega_{\mathrm{p}}=2\pi p\times\text{FSR}$, where $p$ is a positive integer and $\text{FSR}=t_{\mathrm{R}}^{-1}$ is the free-spectral range of the resonator. The fact that the frequency difference $(\omega_{0}-\omega_{\pm})/(2\pi)$ may not, in general, be an integer multiple of the FSR is captured by the additional phase shifts accumulated by the driving fields with respect to the intracavity field from round trip to round trip.

To link the frequency differences $\omega_{0}-\omega_{\pm}$ to the respective phase detunings, we first recall that the phase detuning of a driving field with frequency $\omega$ from a cavity resonance at $\omega^{\prime}$ obeys $\delta\approx(\omega^{\prime}-\omega)t_{\mathrm{R}}$. We can thus write

$$\omega_{0}-\omega_{\pm}\approx\omega^{\prime}_{0}-\frac{\delta_{0}}{t_{\mathrm{R}}}-\omega^{\prime}_{\pm}+\frac{\delta_{\pm}}{t_{\mathrm{R}}},$$

(18)

where the frequency variables with (without) apostrophes refer to resonance (pump) frequencies. We next write the resonance frequencies as $\omega^{\prime}_{q}=\omega^{\prime}_{0}+qD_{1}+\hat{D}_{\mathrm{int}}(q)$, where $D_{1}=2\pi\text{FSR}$ with $\text{FSR}=t^{-1}_{\mathrm{R}}$ the free-spectral range of the cavity (at $\omega^{\prime}_{0}$), $q$ an integer that represents the mode index (with $\omega_{0}$ corresponding to $q=0$), and the integrated dispersion

$$D_{\mathrm{int}}(q)=\sum_{k\geq 2}\frac{D_{k}}{k!}q^{k},$$

(19)

where $D_{k}$ are the expansion coefficients. Assuming that the resonance frequencies $\omega^{\prime}_{\pm}$ are associated with indices $\pm p$ (with $p>0$), respectively, we can write Eq. (18) as

$$\omega_{0}-\omega_{\pm}\approx-\frac{\delta_{0}}{t_{\mathrm{R}}}\mp pD_{1}-D_{\mathrm{int}}(\pm p)+\frac{\delta_{\pm}}{t_{\mathrm{R}}}.$$

(20)

The second term on the right-hand-side of Eq. (20) can be ignored, as it yields an integer multiple of $2\pi$ when used in Eq. (17). Next, we use the fact pasquazi_micro-combs_2018 that the coefficients $D_{k}$ with $k\geq 2$ can be linked to the Taylor series expansion coefficients of the propagation constant $\beta(\omega)$ viz. $D_{k}\approx-D_{1}^{k}L\beta_{k}/t_{\mathrm{R}}$. This allows us to write the integrated dispersion corresponding to the resonance frequencies $\omega^{\prime}_{\pm}$ as

	$\displaystyle D_{\mathrm{int}}(\pm p)$	$\displaystyle\approx\sum_{k\geq 2}-\frac{D_{1}^{k}L\beta_{k}}{t_{\mathrm{R}}k!}(\pm p)^{k},$		(21)
		$\displaystyle=-\frac{L}{t_{\mathrm{R}}}\sum_{k\geq 2}\frac{\beta_{k}}{k!}(\pm D_{1}p)^{k}.$		(22)

Using Eq. (22) in Eq. (20) and substituting the latter into Eq. (17) yields the Ikeda-like map described by Eq. (5) of the main manuscript with coefficients $b_{\pm}$ as defined in Eq. (6) of the manuscript, and the pump frequency shift $\Omega_{\mathrm{p}}=D_{1}p=2\pi p\times\text{FSR}$. We also note that the map can be straightforwardly extended to include arbitrarily many driving fields following the procedure above.

To derive the relationship between the parametric signal detuning $\delta_{0}$ and the pump detunings $\delta_{\pm}$ (i.e., Eq. (7) of the main manuscript), we write out the parametric signal frequency as

$$\omega_{0}=\frac{\omega_{+}+\omega_{-}}{2}=\frac{\omega_{+}^{\prime}-\Delta\omega_{+}+\omega_{-}^{\prime}-\Delta\omega_{-}}{2},$$

(23)

where $\omega_{\pm}^{\prime}$ are the resonance frequencies closest to the pump frequencies and $\Delta\omega_{\pm}$ are the angular frequency detuning of the pump frequencies from those resonances. Substituting $\Delta\omega_{\pm}=\delta_{\pm}/t_{\mathrm{R}}$ and expanding the pump resonance frequencies as $\omega_{\pm}^{\prime}=\omega_{0}^{\prime}\pm pD_{1}+\hat{D}_{\mathrm{int}}(\pm p)$ yields

$$\omega_{0}=\omega_{0}^{\prime}-\frac{\delta_{+}+\delta_{-}}{2t_{\mathrm{R}}}+\frac{\hat{D}_{\mathrm{int}}(p)+\hat{D}_{\mathrm{int}}(-p)}{2}.$$

(24)

Then using Eq. (22) and rearranging, we obtain

	$\displaystyle\delta_{0}$	$\displaystyle=(\omega_{0}^{\prime}-\omega_{0})t_{\mathrm{R}}$
		$\displaystyle=\frac{\delta_{+}+\delta_{-}+L[\hat{D}_{\mathrm{S}}(\Omega_{\mathrm{p}})+\hat{D}_{\mathrm{S}}(-\Omega_{\mathrm{p}})]}{2}.$		(25)

This is Eq. (7) of the main manuscript.

The boundary equation of the Ikeda-like map [Eq. (5) of the main manuscript] shows that the driving fields experience an additional relative phase shift per round trip determined by the coefficient

$$b=\frac{\delta_{+}-\delta_{-}+L[\hat{\beta}_{\mathrm{S}}(\Omega_{\mathrm{p}})-\hat{\beta}_{\mathrm{S}}(-\Omega_{\mathrm{p}})]}{2}.\vspace{2pt}$$

(26)

Using $\delta_{\pm}=(\omega^{\prime}_{\pm}-\omega_{\pm})t_{\mathrm{R}}$ and $L\hat{\beta}_{\mathrm{S}}(\pm\Omega_{\mathrm{p}})=-t_{\mathrm{R}}D_{\mathrm{int}}(\pm p)$, we obtain

$\displaystyle b=t_{\mathrm{R}}\left[\frac{\omega^{\prime}_{+}-\omega^{\prime}_{-}}{2}-\frac{\omega_{+}-\omega_{-}}{2}-\frac{D_{\mathrm{int}}(p)-D_{\mathrm{int}}(-p)}{2}\right].$

(27)

By then using $D_{\mathrm{int}}(p)-D_{\mathrm{int}}(-p)=\omega^{\prime}_{+}-\omega^{\prime}_{-}-2pD_{1}$, we obtain

$\displaystyle b=t_{\mathrm{R}}\left[pD_{1}-\frac{\omega_{+}-\omega_{-}}{2}\right].$

(28)

Given that $\omega_{0}=(\omega_{+}+\omega_{-})/2$, we recognise $(\omega_{+}-\omega_{-})/2$ as the angular frequency shift between the pump at $\omega_{+}$ and the degenerate FWM signal at $\omega_{0}$. Moreover, because the integer $p$ corresponds to mode number of the driven mode $\omega^{\prime}_{+}$ relative to $\omega^{\prime}_{0}$, we may write

$$b=t_{\mathrm{R}}\Delta\omega,$$

(29)

where $\Delta\omega$ is the angular frequency difference between the pump at $\omega_{+}$ and the closest component of the frequency comb (with spacing $D_{1}/(2\pi)$) that forms around $\omega_{0}$. Defining the ordinary comb offset as $\Delta f=|\Delta\omega/(2\pi)|$, and using the fact that the combs that form around $\omega_{0}$ and $\omega_{\pm}$ have the same spacing, we obtain

$$\Delta f=\frac{|b|}{2\pi t_{\mathrm{R}}}=\frac{|b|D_{1}}{(2\pi)^{2}},$$

(30)

where we used $D_{1}=2\pi/t_{\mathrm{R}}$.

Under the assumption that the intracavity envelope $E^{(m)}(z,\tau)$ evolves slowly over a single round trip (i.e., the cavity has a high finesse, and the linear and nonlinear phase shifts are all small), the Ikeda-like map described by Eqs. (2) and (5) of the main manuscript can be averaged into a single mean-field equation. The derivation is well-known haelterman_dissipative_1992 , proceeding by integrating Eq. (2) using a single step of the forward Euler method to obtain $E^{(m)}(L,\tau)$, which is then substituted into Eq. (5). After linearizing with respect to $\delta_{0}$ and $\alpha$ and introducing the slow time variable $t=mt_{\mathrm{R}}$ (such that the round trip index $m=t/t_{\mathrm{R}}$), one obtains:

	$\displaystyle t_{\mathrm{R}}\frac{\partial E(t,\tau)}{\partial t}$	$\displaystyle=\left[-\alpha+i(\gamma L|E|^{2}-\delta_{0})+iL\hat{\beta}_{\mathrm{S}}\left(i\frac{\partial}{\partial\tau}\right)\right]E$
		$\displaystyle+\sqrt{\theta_{+}}E_{\mathrm{in,+}}e^{-i\Omega_{\mathrm{p}}\tau+ib_{+}t/t_{\mathrm{R}}}$
		$\displaystyle+\sqrt{\theta_{-}}E_{\mathrm{in,-}}e^{i\Omega_{\mathrm{p}}\tau+ib_{-}t/t_{\mathrm{R}}}.$		(31)

To obtain the normalized Eq. (10) of the main manuscript, we first introduce the variable transformations $\tau\rightarrow\tau\sqrt{2\alpha/(|\beta_{2}|L)}$, $t\rightarrow\alpha t/t_{\mathrm{R}}$, $\Omega_{\mathrm{p}}\rightarrow\Omega_{\mathrm{p}}\sqrt{|\beta_{2}|L/(2\alpha)}$ and $E\rightarrow E\sqrt{\gamma L/\alpha}$, yielding

	$\displaystyle\frac{\partial E(t,\tau)}{\partial t}$	$\displaystyle=\left[-1+i(|E|^{2}-\Delta_{0})+i\hat{\beta}\left(i\frac{\partial}{\partial\tau}\right)\right]E$		(32)
		$\displaystyle+S_{+}e^{-i\Omega_{\mathrm{p}}\tau+ia_{+}t}+S_{-}e^{i\Omega_{\mathrm{p}}\tau+ia_{-}t},$

where $S_{\pm}=E_{\mathrm{in,\pm}}\sqrt{\gamma L\theta_{\pm}/\alpha^{3}}$, and $\Delta_{0}=\delta_{0}/\alpha$. The normalized dispersion operator $\hat{\beta}$ is defined as

$$\hat{\beta}\left(i\frac{\partial}{\partial\tau}\right)=\sum_{k\geq 2}\frac{d_{k}}{k!}\left(i\frac{\partial}{\partial\tau}\right)^{k},$$

(33)

where the normalized dispersion coefficients are given by

$$d_{k}=\frac{\beta_{k}L}{\alpha}\left(\frac{2\alpha}{|\beta_{2}|L}\right)^{k/2}.$$

(34)

Finally, the coefficients $a_{\pm}=\Delta_{\pm}-\Delta_{0}+\hat{\beta}(\pm\Omega_{\mathrm{p}})$, where $\Delta_{\pm}=\delta_{\pm}/\alpha$ are the normalized detunings of the external driving fields. Note that, for the particular configuration considered in our work, where the signal frequency $\omega_{0}$ is strictly linked to the pump frequencies $\omega_{\pm}$ via $\omega_{0}=(\omega_{+}+\omega_{-})/2$, the coefficients $a_{\pm}=\pm a$, where $a$ is defined by Eq. (11) of the main manuscript.