Slow-light frequency combs and dissipative Kerr solitons in coupled-cavity waveguides

TL;DR

This paper investigates Kerr frequency combs and dissipative Kerr solitons in silicon photonic crystal coupled-cavity waveguides, highlighting how dispersion engineering and slow-light effects influence comb generation thresholds.

Contribution

It introduces a detailed analysis of comb generation thresholds in CCWs considering mode volume, quality factor, and group index, advancing understanding of slow-light enhanced nonlinear optics.

Findings

Threshold depends on mode volume, quality factor, and group index.

Dispersion engineering enables optimized comb generation.

Slow-light effects significantly enhance Kerr comb formation.

Abstract

We study Kerr frequency combs and dissipative Kerr solitons in silicon photonic crystal coupled-cavity waveguides (CCW) with globally optimized dispersion at telecom wavelengths. The corresponding threshold for comb generation is found to explicitly depend on the main CCW figures of merit, namely, mode volume, normal mode quality factor and slow-light group index. Our analysis is carried out by solving the non-linear dynamics of the CCW Bloch modes in presence of Kerr non-linearity and two-photon absorption. Our results open the way to CCW comb generation via dispersion engineering and slow-light enhancement.

Equations

$${\left|{\mathcal{A}}_{{\alpha }_0}\right|}_{\mathrm{th}}^2=\frac{ϵV_c}{2l{n}_{g,{\alpha }_0}{n}_2{Q}_{{\alpha }_0}}f(\kappa )=138\text{mbox}mW,$$

$$\nabla \times \mathbf{E}(\mathbf{r},t)+{\mu }_0\frac{\partial \mathbf{H}(\mathbf{r},t)}{\partial t}$$

$$\nabla \times \mathbf{H}(\mathbf{r},t)-{\epsilon }_0ϵ(\mathbf{r},|\mathbf{E}{|}^2)\frac{\partial \mathbf{E}(\mathbf{r},t)}{\partial t}$$

$$\left[\hat{\mathcal{L}}+\frac{ϵ(\mathbf{r},|\mathbf{E}{|}^2)}{c^2}\frac{{\partial }^2}{\partial t^2}\right]\mathbf{E}(\mathbf{r},t)=0,$$

$$\mathbf{E}(\mathbf{r},t)=\sum _{\mu }{\mathcal{A}}_{\mu }^{{}^{\prime }}(t)e^{-i{\omega }_{\mu }t}{\mathbf{E}}_{\mu }(\mathbf{r})+E_{\mathrm{ext}}{e}^{-i{\Omega }_{0}t}{\hat{e}}_0,$$

$${\mathbf{E}}_{\mu }(\mathbf{r})=\sqrt{N_{\mu }}{\vec{\mathcal{E}}}_{\mu }(\mathbf{r}),$$

$$\hat{\mathcal{L}}{\vec{\mathcal{E}}}_{\mu }(\mathbf{r})=\frac{{\tilde{\omega }}_{\mu }^2}{c^2}ϵ(\mathbf{r}){\vec{\mathcal{E}}}_{\mu }(\mathbf{r}),$$

$${\int }_lϵ(\mathbf{r}){\vec{\mathcal{E}}}_{\alpha }^{\ast }(\mathbf{r})\cdot {\vec{\mathcal{E}}}_{\mu }(\mathbf{r})dV={\delta }_{\alpha ,\mu },$$

$${\tilde{\omega }}_{\mu }={\omega }_{\mu }-i\frac{{\gamma }_{\mu }}{2},$$

$$v_g(\mu )=\frac{{\int }_l\hat{\mu }\cdot {\mathbf{S}}_{\mu }dV}{{\mathcal{U}}_{{\mathbf{E}}_{\mu }}+{\mathcal{U}}_{{\mathbf{H}}_{\mu }}},$$

$$v_g(\mu )=\frac{{\int }_{l}dy\left[{\int }_{\infty }dxdz{S}_{\mu ,y}\right]}{2{\mathcal{U}}_{{\mathbf{E}}_{\mu }}},$$

$$v_g(\mu )=\frac{l{P}_{\mu }}{2{\mathcal{U}}_{{\mathbf{E}}_{\mu }}}.$$

$$v_g(\mu )=\frac{l{P}_{\mu }}{\frac{1}{2}{\int }_l{\epsilon }_0ϵ(\mathbf{r})|{\mathbf{E}}_{\mu }(\mathbf{r}){|}^{2}dV},$$

$${\int }_lϵ(\mathbf{r})|{\mathbf{E}}_{\mu }(\mathbf{r}){|}^{2}dV=\frac{2l{P}_{\mu }}{{\epsilon }_0{v}_g(\mu )}.$$

$$N_{\mu }=\frac{2l{P}_{\mu }}{{\epsilon }_0{v}_g(\mu )},$$

$$\mathbf{E}(\mathbf{r},t)=\sqrt{\frac{2l}{{\epsilon }_0}}\sum _{\mu }\sqrt{P_{\mu }}{\mathcal{A}}_{\mu }^{{}^{\prime }}(t)e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_g(\mu )}}{\vec{\mathcal{E}}}_{\mu }(\mathbf{r})+E_{\mathrm{ext}}{e}^{-i{\Omega }_{0}t}{\hat{e}}_0.$$

$$\mathbf{E}(\mathbf{r},t)=\sqrt{\frac{2l}{{\epsilon }_0}}\sum _{\mu }{\mathcal{A}}_{\mu }(t)e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}{\vec{\mathcal{E}}}_{\mu }(\mathbf{r})+E_{\mathrm{ext}}{e}^{-i{\Omega }_{0}t}{\hat{e}}_0,$$

$${\int }_Lϵ(\mathbf{r}){\vec{\mathcal{E}}}_{\alpha }^{\ast }(\mathbf{r})\cdot {\vec{\mathcal{E}}}_{\mu }(\mathbf{r})dV=M{\delta }_{\alpha ,\mu },$$

$$\frac{{\partial }^2\mathbf{E}(\mathbf{r},t)}{\partial t^2}=\sqrt{\frac{2l}{{\epsilon }_0}}\sum _{\mu }\left[{\ddot{\mathcal{A}}}_{\mu }(t)-2i{\omega }_{\mu }{\dot{\mathcal{A}}}_{\mu }(t)-{\omega }_{\mu }^2{\mathcal{A}}_{\mu }(t)\right]e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}{\vec{\mathcal{E}}}_{\mu }(\mathbf{r})-{\Omega }_0^2{E}_{\mathrm{ext}}{e}^{-i{\Omega }_{0}t}{\hat{e}}_0,$$

$$\hat{\mathcal{L}}\mathbf{E}(\mathbf{r},t)=\sqrt{\frac{2l}{{\epsilon }_0}}\frac{1}{c^2}\sum _{\mu }{\tilde{\omega }}_{\mu }^2{\mathcal{A}}_{\mu }(t)e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}ϵ(\mathbf{r}){\vec{\mathcal{E}}}_{\mu }(\mathbf{r}),$$

$$\sum _{\mu }{\tilde{\omega }}_{\mu }^2{\mathcal{A}}_{\mu }(t)e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}ϵ(\mathbf{r}){\vec{\mathcal{E}}}_{\mu }(\mathbf{r})$$

$$+ϵ(\mathbf{r},|\mathbf{E}{|}^2)\left\{\sum _{\mu }\left[{\ddot{\mathcal{A}}}_{\mu }(t)-2i{\omega }_{\mu }{\dot{\mathcal{A}}}_{\mu }(t)-{\omega }_{\mu }^2{\mathcal{A}}_{\mu }(t)\right]e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}{\vec{\mathcal{E}}}_{\mu }(\mathbf{r})-\sqrt{\frac{{\epsilon }_0}{2l}}{\Omega }_0^2{E}_{\mathrm{ext}}{e}^{-i{\Omega }_{0}t}{\hat{e}}_0\right\}=0.$$

$$\sum _{\mu }{\tilde{\omega }}_{\mu }^2{\mathcal{A}}_{\mu }(t)e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}ϵ(\mathbf{r}){\vec{\mathcal{E}}}_{\mu }(\mathbf{r})$$

$$+ϵ(\mathbf{r},|\mathbf{E}{|}^2)\sum _{\mu }\left[{\ddot{\mathcal{A}}}_{\mu }(t)-2i{\omega }_{\mu }{\dot{\mathcal{A}}}_{\mu }(t)-{\omega }_{\mu }^2{\mathcal{A}}_{\mu }(t)\right]e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}{\vec{\mathcal{E}}}_{\mu }(\mathbf{r})-ϵ(\mathbf{r})\sqrt{\frac{{\epsilon }_0}{2l}}{\Omega }_0^2{E}_{\mathrm{ext}}{e}^{-i{\Omega }_{0}t}{\hat{e}}_0=0.$$

$$ϵ(\mathbf{r},|\mathbf{E}{|}^2)={\left[n(\mathbf{r})+n_2(\mathbf{r})I(\mathbf{r},|\mathbf{E}{|}^2)+i\frac{c}{2{\omega }_0}{\beta }_{\mathrm{TPA}}(\mathbf{r})I(\mathbf{r},|\mathbf{E}{|}^2)\right]}^2,$$

$$I(\mathbf{r},|\mathbf{E}{|}^2)=\frac{{\epsilon }_{0}c}{2}n(\mathbf{r})|\mathbf{E}{|}^2,$$

$$ϵ(\mathbf{r},|\mathbf{E}{|}^2)\backsimeq ϵ(\mathbf{r})+2n(\mathbf{r})\Delta n(\mathbf{r},|\mathbf{E}{|}^2),$$

$$\Delta n(\mathbf{r},|\mathbf{E}{|}^2)=\frac{{\epsilon }_{0}c}{2}n(\mathbf{r})\left[n_2(\mathbf{r})+i\frac{c}{2{\omega }_0}{\beta }_{\mathrm{TPA}}(\mathbf{r})\right]|\mathbf{E}{|}^2.$$

$$\sum _{\mu }({\tilde{\omega }}_{\mu }^2-{\omega }_{\mu }^2){\mathcal{A}}_{\mu }(t)e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}ϵ(\mathbf{r}){\vec{\mathcal{E}}}_{\mu }(\mathbf{r})+\sum _{\mu }\left[{\ddot{\mathcal{A}}}_{\mu }(t)-2i{\omega }_{\mu }{\dot{\mathcal{A}}}_{\mu }(t)\right]e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}ϵ(\mathbf{r}){\vec{\mathcal{E}}}_{\mu }(\mathbf{r})$$

$$-2\sum _{\mu }{\omega }_{\mu }^2{\mathcal{A}}_{\mu }(t)e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}n(\mathbf{r})\Delta n(\mathbf{r},|\mathbf{E}{|}^2){\vec{\mathcal{E}}}_{\mu }(\mathbf{r})$$

$$+2\sum _{\mu }\left[{\ddot{\mathcal{A}}}_{\mu }(t)-2i{\omega }_{\mu }{\dot{\mathcal{A}}}_{\mu }(t)\right]e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}n(\mathbf{r})\Delta n(\mathbf{r},|\mathbf{E}{|}^2){\vec{\mathcal{E}}}_{\mu }(\mathbf{r})-ϵ(\mathbf{r})\sqrt{\frac{{\epsilon }_0}{2l}}{\Omega }_0^2{E}_{\mathrm{ext}}{e}^{-i{\Omega }_{0}t}{\hat{e}}_0=0.$$

$$({\tilde{\omega }}_{\alpha }^2-{\omega }_{\alpha }^2){\mathcal{A}}_{\alpha }(t)e^{-i{\omega }_{\alpha }t}\frac{1}{\sqrt{v_{\alpha }}}+\left[{\ddot{\mathcal{A}}}_{\alpha }(t)-2i{\omega }_{\alpha }{\dot{\mathcal{A}}}_{\alpha }(t)\right]e^{-i{\omega }_{\alpha }t}\frac{1}{\sqrt{v_{\alpha }}}$$

$$-\frac{2}{M}\sum _{\mu }{\omega }_{\mu }^2{\mathcal{A}}_{\mu }(t)e^{-i{\omega }_{\mu }t}\frac{1}{\sqrt{v_{\mu }}}{\int }_{L}n(\mathbf{r})\Delta n(\mathbf{r},|\mathbf{E}{|}^2){\vec{\mathcal{E}}}_{\alpha }^{\ast }(\mathbf{r})\cdot {\vec{\mathcal{E}}}_{\mu }(\mathbf{r})dV$$

Full text

I Introduction

Kerr frequency combs have revolutionized several fields in optical sciences since they were first proposed in monolithic microresonators Del’Haye et al. (2007). To date, they have been successfully applied to a vast variety of state-of-the-art technologies and highly sophisticated measurement techniques as ultra-high resolution spectroscopy Coddington et al. (2016); Picqué and Hänsch (2019), massively parallel coherent optical communications Marin-Palomo et al. (2017), light detection and ranging (LIDAR) Trocha et al. (2018), optical-frequency synthesizers Spencer et al. (2018) and microphotonic astrocombs Benedick et al. (2010); Obrzud et al. (2019). Frequency combs in optical microresonators are generated by the interaction of either, a continuous-wave (cw) pump laser or pumped optical pulses Obrzud et al. (2017), with the resonator modes via parametric four-wave mixing (FWM), assisted by the Kerr non-linearity of the material. This parametric process fulfills the energy conservation and is enhanced when the side bands created by FWM coincide with the resonator frequencies. Since the intensity-dependent refractive index induces a frequency shift on the modes, an increasing free spectral range (FSR) with frequency (anomalous dispersion) is required to compensate this effect and effectively produce equidistant spectral lines suitable to support the cascaded FWM Kippenberg et al. (2011). When a high number of mode-locked modes are involved in this non-linear interaction and dissipation cannot be neglected, the complex non-linear dynamics may give rise to disipative Kerr solitons (DKS), which arise as a double balance between non-linearity and dispersion (preserving their shape), as well as dissipation and parametric gain (preserving their amplitude) Kippenberg et al. (2018). Moreover, multiple-soliton formations have become very interesting in the context of localized-state interactions Vladimirov et al. (2018); Parra-Rivas et al. (2017); Milián et al. (2018), breather solitons Yu et al. (2017); Lucas et al. (2017); Guo et al. (2017) and soliton crystals Cole et al. (2017); Karpov et al. (2019). DKS are notably relevant because their corresponding spectra exhibit highly coherent frequency combs with perfect single-FSR spacing between side bands Herr et al. (2014); Gaeta et al. (2019), and dispersive waves in presence of high-order dispersion Coen et al. (2013); Milián and Skryabin (2014); Parra-Rivas et al. (2014); Brasch et al. (2016). The latter is particularly important because it allows a further spanning of the frequency combs, highlighting the important role of dispersion engineering on parametric gain, and consequently, on comb generation Foster et al. (2006). Such task has been successfully achieved in non-linear microresonators by varying the ring geometrical parameters in order to tailor the dispersion for specific comb functionalities Zhang et al. (2013); Okawachi et al. (2014); Grudinin and Yu (2015). Nevertheless, this geometry has a limited parameter space, thus restricting the choice of materials, operation wavelength and size of the final device. Recently, Fabry-Perot resonators have been proposed as appealing candidates for frequency comb generation since different methods to reshape the cavity dispersion can be applied Cole et al. (2018); Yu et al. (2019), however, advanced dispersion engineering still remains challenging.

In this work, we propose a coupled-cavity waveguide (CCW, also known as coupled-resonator optical waveguide or CROW) as a new candidate to efficiently generate low-threshold Kerr frequency combs and DKS. Similar to ring resonators, CCWs define a discrete spectrum of propagating modes which may trigger cascaded parametric FWM in the anomalous dispersion regime, and lead to a highly coherent comb if the waveguide dispersion is conveniently optimized. In fact, since their appearance two decades ago Yariv et al. (1999), CCWs have shown to be extremely flexible for advanced dispersion engineering and enhancement of non-linear phenomena via slow-light Soljačić et al. (2002); Chen and Blair (2004); O’Brien et al. (2007); Notomi et al. (2008); Matsuda et al. (2011); Minkov and Savona (2015). We show in Fig. 1(a) a schematic representation of the CCW studied in this work. The periodic waveguide is formed by a staggered distribution of cavities where first (dashed red) and second (dashed green) neighbor cavity-cavity coupling is assumed. The system is pumped with a cw source and a DKS is propagating along the waveguide direction $y$.

II Non-linear photonic crystal CCW

The system of Fig. 1(a) is realized with the waveguide of coupled L3 photonic crystal slab cavities shown in Fig. 1(b), where we have plotted the intensity profile of the displacement field [$\mathbf{D}(\mathbf{r})=\epsilon(\mathbf{r})\vec{{\cal E}}(\mathbf{r})$] at the boundary of the Brillouin zone. The photonic crystal is formed by a hexagonal lattice of holes, pitch $a=400$ nm and hole radii $r=0.25a$, etched in a silicon membrane of thickness $d=0.55a$, and the L3 cavity is introduced by removing three holes along the $\Gamma K$ direction of the lattice. This CCW configuration allows us to tune first and second neighbor coupling between the L3 cavities by varying the size of the red and green holes, respectively, and optimize the out-of-plane losses by varying the position and size of the white ones. The waveguide of Fig. 1(b) has been previously proposed as a compact CCW with outstanding figures of merit obtained via automated global optimization Minkov and Savona (2015) and successfully measured experimentally Lai et al. (2018); Mohamed et al. (2018). Here, we adopt the design with largest averaged group index and small second order dispersion reported in Ref. 38 with parameters $(\Delta r_{1},\Delta r_{2},\Delta r_{3},\Delta x)=(-0.0049,-0.0340,-0.1016,0.2204)a$, where $r_{1}=r+\Delta r_{1}$, $r_{2}=r+\Delta r_{2}$ and $r_{3}=r+\Delta r_{3}$ are the radii of the red, green and white holes, respectively, and $\Delta x$ is the outward displacement of the latter ones. The dispersion relation $\nu_{\alpha}$ of the waveguide and decay rates $\gamma_{\alpha}$ of the Bloch modes of momentum $\alpha$ are calculated with the guided mode expansion method (GME) Andreani and Gerace (2006) for a system length of $M=400$ cavities (400 normal modes within the Brillouin zone of the CCW) with period $l=\sqrt{3}a$ (waveguide length $L=Ml$), while the non-linear mode volume is computed with a commercial FDTD solver lum (2019) and found to be $V_{c}=0.345~{}\mu$m3. This sample length represents a good compromise between total number of normal modes [avoiding finite-size effects in the coupled mode equations Eq. (56)] and robustness to fabrication disorder Song et al. (2010); Mohamed et al. (2018). Silicon parameters at telecom frequencies are considered for the material, namely, dielectric constant $\epsilon=12.04$, Kerr coeficient $n_{2}=5.52\times 10^{-18}$ m${}^{2}/$W and two-photon absorption (TPA) coefficient $\beta_{\rm TPA}=1\times 10^{-11}$ m$/$W Hon et al. (2011). The GME dispersion relation, second order dispersion and group index are shown in Fig. 2(a)-2(c), respectively, for which red and blue segments correspondingly highlight the regions of normal and anomalous dispersion. Notice that these dispersion curves are valid for a straight waveguide with periodic boundary conditions, nevertheless, we expect them to describe very well a close loop of coupled cavities as bending losses is usually very small, with respect to out-of-plane losses, in photonic crystal geometries Zhang and Li (2006). The system is pumped at $a\alpha_{0}=1.2742$, where $\nu_{\alpha_{0}}=193.39$ THz, $(ac)^{-1}(d^{2}\nu/\alpha^{2})_{\alpha_{0}}=8.63\times 10^{-4}$, $n_{g,\alpha_{0}}=119.34$ and the normal mode quality factor is $Q_{\alpha_{0}}=5.72\times 10^{4}$. This sets an internal mode threshold given by

[TABLE]

where ${\cal A}_{\alpha_{0}}$ is the slowly-varying amplitude of the Bloch mode $\alpha_{0}$, and $f(\kappa)=(\sqrt{1+\kappa^{2}}+2\kappa)/(1-3\kappa^{2})$, with $\kappa=c\beta_{\rm TPA}/(2n_{2}\omega_{\alpha_{0}})$, is a function of the material parameters only which gives $f(\kappa=0.2236)=1.73$ for silicon at $\omega_{\alpha_{0}}/2\pi=\nu_{\alpha_{0}}=193.39$ THz, and reduces to $f(\kappa=0)=1$ for $\beta_{\rm TPA}=0$. This result is particularly remarkable because, even by considering that TPA increases the threshold by a factor of 1.73, $|{\cal A}_{\alpha_{0}}|_{\rm th}^{2}=138~{}$mW is still one order of magnitude smaller than in typical mm-size crystalline non-linear ring resonators Chembo (2016). Equation (1) undoubtedly strengths the potential capabilities of CCW structures for low-threshold comb generation. While the factor $V/Q$ also enters in the expression for the internal mode threshold of ring resonators Chembo et al. (2010); Chembo and Yu (2010); Chembo and Menyuk (2013), CCWs can drastically decrease this minimal power in the slow-light regime, i.e., for $n_{g,\alpha_{0}}\gg\sqrt{\epsilon}$. Such structural effect allows to access the rich physics of DKS in CCWs at much lower powers than in their monolithic counterparts even in presence of strong TPA, which is the main source of non-linear losses in silicon structures at telecom wavelengths Lau et al. (2015); Sahoo et al. (2017). It is worth emphasizing that, since we are not specifying any particular coupling architecture with the external pump, we are referring to the internal mode threshold which scales as $1/Q$, and not to the external power threshold which is found to scale as $1/Q^{2}$ for ring resonators at critical coupling with the external source Herr et al. (2012). While the latter is an important parameter to set the excitation power in an experimental setup, the former quantity, also called the minimum intracavity power for comb generation Chembo (2016), gives the corresponding effective power into the specific normal mode of the system. We plot in Fig. 2(d) the driven intensity $|{\cal F}_{\alpha_{0}}|^{2}$, in $|{\cal A}_{\alpha_{0}}|_{\rm th}^{2}$ units, evaluated at threshold (continuous red lines) and at the bistability boundaries (continuous blue curves), as a function of the laser detuning $\sigma_{\alpha_{0}}=\Omega_{0}-\omega_{\alpha_{0}}$ in units of $\gamma_{\alpha_{0}}$, where the dashed lines are the corresponding curves for $\beta_{\rm TPA}=0$. The effects of TPA on the threshold are clearly seen in this plot, moreover, the hysteresis region (where the system has bistable states) is red-shifted and the area between the bistability boundaries is significantly reduced. Hence, the system needs to be pumped stronger and the laser frequency has to be further decreased to see Kerr frequency combs and DKS in presence of TPA.

III Frequency comb and DKS solutions

In order to find the possible steady state solutions determined by the non-linear CCW dynamics, we carry out a frequency scan at $|{\cal F}_{\alpha_{0}}|^{2}=6|{\cal A}_{\alpha_{0}}|_{\rm th}^{2}$ in Fig. 2(d) within the $\sigma_{\alpha_{0}}/\gamma_{\alpha_{0}}$ interval $[-3,1]$. The dynamical non-linear coupled-mode equations Eq. (56) are thus propagated in time for each $\sigma_{\alpha_{0}}$ value along this trajectory using an explicit Runge-Kutta integrator and fast Fourier transform (see Ref. 53) until the steady state is reached. Our simulation takes into account the momentum-dependent decay rate $\gamma_{\alpha}$ of the Bloch modes, and the initial condition of the integrator is a sharp Gaussian pulse $\psi(y,0)=\exp[-0.5(y/l)^{2}]$, where $\psi(y,t)=\sum_{\alpha}{\cal A}_{\alpha}(t)e^{i\sigma_{\alpha}t}e^{-i(\alpha-\alpha_{0})y}$ is the envelope function along the CCW direction. Results of this analysis are presented in Fig. 2(e) where we show the averaged waveguide power in the steady state $\langle|\psi^{s}|^{2}\rangle=\int_{L}|\psi^{s}(y)|^{2}dy$, in units of threshold, as a function of the laser detuning. When the laser frequency is red-shifted the intra-waveguide power start to increase until a clear series of discrete steps (see inset), decreasing the average power, are seen in the yellow area, which corresponds to the region where there is hysteresis in the system and the soliton solutions are expected to appear Kippenberg et al. (2018); Godey et al. (2014). In fact, these steps have been previously measured experimentally in the transmission of non-linear ring resonators, and they are associated to the formation of different soliton states within the system, in excellent agreement with the predictions of the non-linear coupled mode equations and the Lugiato-Lefever formalism Herr et al. (2014). We show in Figs. 2(f)-2(i) the respective frequency combs at the four representative points marked in Fig. 2(e). Figures 2(f) and 2(g) displays a sequence of primary combs separated by 40-FSR and 41-FSR, respectively, which correspond to supercritical Turing patterns as they are excited above threshold. In Figs. 2(h)-2(i), both combs are subcritical (excited right below threshold) with a single FSR spacing, and they are the signature of soliton complexes. Notably, because of the presence of non-trivial high order terms in the dispersion relation of the photonic crystal CCW, these frequency combs have signatures of dispersive wave formation or soliton Cherenkov radiation at the linear phase matching condition $\omega_{\alpha}-[\omega_{\alpha_{0}}+\zeta_{1}(\alpha-\alpha_{0})]=0$, where $\zeta_{1}$ represents the group velocity Brasch et al. (2016). Particularly, the dispersive peak clossest to the zero comb line is predicted around $a\alpha\simeq 1.63$, which is in perfect agreement with our numerical simulations. The corresponding steady state intensity of the envelope functions are shown in Fig. 3. For the I (40 FSR repetition rate) and II (41 FSR repetition rate) states, we see 40 and 41 Turing rolls, respectively, in agreement with the Lugiato-Lefever theory of Kerr frequency combs in the anomalous dispersion regime Godey et al. (2014). Moreover, for the III and IV states we clearly identify, respectively, a soliton molecule composed of two pulses and a single soliton pulse, propagating along the waveguide direction while keeping their shape and amplitude. It must be said that because these soliton structures are subcritical (pumped below threshold), they are extremely sensitive to the initial conditions, and therefore, the system may follow different trajectories in Fig 2(e) when slightly changing the Bloch mode amplitudes at $t=0$ Godey et al. (2014); Herr et al. (2014). Furthermore, the characteristic width of the soliton pulse is much smaller than the waveguide length thus ruling out finite-size effects in the basis expansion. While results presented in Fig. 2 are strictly valid for a CCW system with periodic boundary conditions, we still expect highly coherent combs in finite waveguides where non-negligible corrections to the non-linear interaction are predicted by the theory of frequency combs in Fabry-Perot resonators Cole et al. (2018). Nevertheless, a rigorous quantitative analysis is out of the scope of this paper, and will be focus of a future work.

Figure 4 shows the frequency spectrum of the DKS in Fig. 3(d), where the unloaded modes (normal modes of the CCW) are represented by the vertical dashed lines. As expected from energy conservation, the spectral lines are equally spaced with a repetition rate of 9.1 GHz. Such rates have been previously achieved in non-linear ring resonators of 7 mm diameter Suh and Vahala (2018), in contrast to our system length of $\sim 277~{}\mu$m. Notice that the minimal external power for comb generation may significantly vary depending on the material, geometry and more importantly, the coupling architecture with a bus waveguide Chembo (2016), however, this important result notably highlights the potential capabilities of $\mu$m-scale low-threshold CCW systems to generate low-rate frequency combs, which are desirable in high-precision comb applications Kippenberg et al. (2018).

It is important to stress that the specific choice of photonic crystal CCWs is mainly motivated by their characteristic diffraction-limited mode volumes Kristensen and Hughes (2014) which play an important role on decreasing the threshold power and final device size. Nevertheless, the results presented in this section only rely on the dispersion relation of the CCW and are therefore expected for any CCW system, e.g. coupled microrings, displaying anomalous dispersion with, more generally, $n$-th neighbor coupling.

IV Conclusions

In conclusion, we have presented results on Kerr frequency combs and DKS in CCWs where the main source of non-linear loses is given by TPA. We found that the internal mode threshold for comb generation depends on the main CCW figures of merit, namely, cavity volume, normal mode quality factor and group index. While TPA losses has the main effect of increasing the threshold power and inducing a red shift in the optimal laser detuning, structural slow-light plays an important role on reducing the minimal power to trigger FWM phenomena between the CCW normal modes, evidencing the capabilities of CCW systems for low-threshold frequency combs. Specifically, we have demonstrated the possibility of DKS in a realistic dispersion-engineered silicon photonic crystal CCW at telecom wavelengths, where highly-complex combs were seen with signatures of soliton Cherenkov radiation. Repetition rates of few gigahertz in a $\sim 277~{}\mu$m-length photonic crystal CCW have been obtained into the soliton regime, which are commonly achieved in $\sim 7$ mm size ring resonators, thus demonstrating the potential of CCW systems for high-precision comb applications in ultra-compact devices. Finally, although the spectral spanning of the combs generated in CCW is limited to the waveguide bandwidth, our results open the way to Kerr frequency comb and DKS generation via advanced dispersion engineering and slow-light non-linear enhancement.

Appendix A Electromagnetic wave equation

We start from Maxwell’s curl equations

[TABLE]

where the material is assumed to be non-dispersive, non-linear and isotropic. If we apply the curl operator to Eq. 2 and then use Eq. 3 we get the wave equation for $\mathbf{\mathbf{E}(\mathbf{r},t)}$

[TABLE]

where $\hat{\mathcal{L}}=\nabla\times\nabla\times$ and $c=1/\sqrt{\varepsilon_{0}\mu_{0}}$.

Appendix B Electric field expansion

Because the absolute intensity of the electric field determines the effective strength of the Kerr non-linearity, it is very important to write $\mathbf{E}(\mathbf{r},t)$ in proper units. We specifically adopt the SI system, where the electric field has units of [V/m], and expand it in terms of the of the normal modes of the coupled-cavities waveguide (CCW):

[TABLE]

where ${\cal A}^{{}^{\prime}}_{\mu}(t)$ is the slowly-varying amplitude of the normal mode $\mathbf{E}_{\mu}(\mathbf{r})$ with eigenfrequency $\omega_{\mu}$ and momentum $\mu$ (along the waveguide direction $\hat{\mu}$), and $E_{\rm ext}$ is the electric field amplitude of the driving continuous-wave (cw) source with frequency $\Omega_{0}$ and vector direction $\hat{e}_{0}$. Equation 5 is the so called slowly-varying amplitude expansion, as the time and spatial variations have been separated, and the fast-varying temporal dependence is explicitly represented by the oscillatory term $e^{-i\omega_{\mu}t}$. Notice that the field $\mathbf{E}_{\mu}(\mathbf{r})$ in Eq. 5 also has units of [V/m], then ${\cal A}^{{}^{\prime}}$ is dimensionless. In order to express the CCW normal modes in a more convenient form we write

[TABLE]

where $N_{\mu}$ is a normalization factor and $\vec{{\cal E}}_{\mu}(\mathbf{r})$ fulfills the generalized eigenvalue equation

[TABLE]

subject to the normalization condition

[TABLE]

with $l$ representing the period of the waveguide and $\epsilon(\mathbf{r})$ its dielectric function. Equation 8 sets the units of $\vec{{\cal E}}_{\mu}(\mathbf{r})$ to [m*-3/2*]. In Eq. 7 we have considered the complex frequency

[TABLE]

where $\gamma_{\mu}$ is the overall loss rate of the Bloch mode $\vec{{\cal E}}_{\mu}(\mathbf{r})$. For the particular case of a photonic crystal slab, $\gamma_{\mu}$ represents the out-of-plane losses due to the coupling with the leaky modes of the homogeneous system Andreani and Gerace (2006). Notice that the normalization condition of Eq. 8, satisfied by normal modes, is not rigorously valid when the eigenvalue of Eq. 7 is complex, i.e., when dealing with quasinormal modes Kristensen and Hughes (2014). Nevertheless it still represents a very good approximation in the low loss limit $\gamma_{\mu}\ll\omega_{\mu}$, which is the case of interest in the present formulation. The normalization factor $N_{\mu}$ can be easily found by recalling the fundamental relation between the group velocity $v_{g}(\mu)$, averaged energy flux and averaged energy density of the electromagnetic field, which holds in periodic waveguides Joannopoulos et al. (2008)

[TABLE]

where $\mathbf{S}_{\mu}(\mathbf{r})=\mbox{Re}[\mathbf{E}^{\ast}_{\mu}(\mathbf{r})\times\mathbf{H}^{\ast}_{\mu}(\mathbf{r})]/2$ is the time-averaged Poynting vector (energy flux), ${\cal U}_{\mathbf{E}_{\mu}}=\int_{l}\varepsilon_{0}\epsilon(\mathbf{r})|\mathbf{E}_{\mu}(\mathbf{r})|^{2}dV/4$ is the physical electric field energy (integrated density energy) and ${\cal U}_{\mathbf{H}_{\mu}}=\int_{l}\mu_{0}|\mathbf{H}_{\mu}(\mathbf{r})|^{2}dV/4$ the magnetic one in the waveguide period. If we set the waveguide direction along the $y$ axis, i.e., $\hat{\mu}=\hat{y}$, and use the result that for harmonic modes the electric and magnetic field energies are equal Joannopoulos et al. (2008), i.e., ${\cal U}_{\mathbf{E}_{\mu}}={\cal U}_{\mathbf{H}_{\mu}}$, Eq. 10 takes the following form

[TABLE]

where the last integral is carried out over all space and it is the power flux $P_{\mu}=\int_{\infty}dxdzS_{\mu,y}$ through the transverse cross section of the waveguide, which has been shown to be independent of $y$ Song et al. (2010), therefore

[TABLE]

Using the definition of ${\cal U}_{\mathbf{E}_{\mu}}$, Eq. 12 turns into

[TABLE]

which means that the electric field $\mathbf{E}_{\mu}(\mathbf{r})$ must fulfill

[TABLE]

The expression for $N_{\mu}$ is readily found from Eqs. 6, 8 and 14

[TABLE]

and the field expansion in Eq. 5 is thus expressed in terms of the electric field eigenmodes $\vec{{\cal E}}_{\mu}(\mathbf{r})$ as follows

[TABLE]

We now define the new amplitude ${\cal A}_{\mu}(t)=\sqrt{P_{\mu}}{\cal A}^{{}^{\prime}}_{\mu}(t)$ and simplify the notation for the group velocity of the mode $\mu$ as $v_{g}(\mu)\rightarrow v_{\mu}$. The electric field expansion takes the following final form

[TABLE]

where $|{\cal A}_{\mu}(t)|^{2}$ is the instantaneous power, in [W] units, of the mode $\mu$ propagating in the waveguide direction.

Appendix C Non-linear coupled mode equations

Let’s now consider a CCW of $M$ cavities where the total electric field of the system is represented by Eq. 17 and the normalization condition of Eq. 8 is written as

[TABLE]

when integrating over the total waveguide length $L=Ml$. The second time-derivative of Eq. 17 reads

[TABLE]

while the operator $\hat{\mathcal{L}}$ applied to Eq. 17 gives

[TABLE]

where we have used Eq. 7. Because $\mathbf{E}(\mathbf{r},t)$ must satisfy the wave equation in Eq. 4, we employ Eqs. 19 and 20 to get

[TABLE]

The intensity dependent dielectric function in Eq. C can be set to the reference value $\epsilon(\mathbf{r})=\epsilon(\mathbf{r},0)$ in the pump contribution as the perturbation induced by $|\mathbf{E}|^{2}$ can be neglected when driving the system at $\epsilon(\mathbf{r})$. Therefore, Eq. C turns into

[TABLE]

At this point we explicitly write the intensity dependent dielectric function as Boyd (2007)

[TABLE]

where $n(\mathbf{r})=\sqrt{\epsilon(\mathbf{r})}$ is the linear refractive index, $n_{2}(\mathbf{r})$ and $\beta_{\rm TPA}(\mathbf{r})$ are the spatial dependent Kerr and two-photon absorption (TPA) coefficients, respectively, of the material at $\omega_{0}$ and

[TABLE]

is the instantaneous electric field intensity. As the non-linear shift induced on the dielectric function is expected to be small, Eq. 23 can be approximated to

[TABLE]

with

[TABLE]

Replacing Eq. 25 in Eq. C gives (after reorganizing terms)

[TABLE]

In order to arrive to a system of differential equations for the slowly-varying amplitudes, Eq. C is multiplied by $\vec{{\cal E}}^{\ast}_{\alpha}(\mathbf{r})$ and integrated over the whole waveguide. After the orthogonality condition of Eq. 18 is explicitly applied, we get

[TABLE]

We now apply the slowly-varying assumption where the second derivative terms are assumed to be small, i.e., $|\ddot{{\cal A}}_{\alpha}(t)|\ll|2\omega_{\alpha}\dot{{\cal A}}_{\alpha}(t)|$, and consequently they can be neglected, as well as the first derivative term multiplying the non-linear perturbation $\Delta n(\mathbf{r},|\mathbf{E}|^{2})$ in the second integral of Eq. C. We also assume low-loss CCWs, for which we can approximate $\tilde{\omega}_{\alpha}^{2}-\omega_{\alpha}^{2}\backsimeq-i\omega_{\alpha}\gamma_{\alpha}$ in the first term of the equation. Therefore, Eq. C is rewritten as

[TABLE]

where we have defined the driving amplitude

[TABLE]

Let’s focus on the non-linear term of Eq. 29. From the field expansion in Eq. 17, the intra-waveguide modulus squared field is given by

[TABLE]

therefore, the non-linear perturbation in Eq. 26 is written as

[TABLE]

and the integral of Eq. 29 becomes

[TABLE]

where we have defined the complex Kerr coefficient

[TABLE]

In order to write the last integral term of Eq. 33 in a more convenient form, we assume that the normal modes of the CCW can be expanded in the cavity mode basis

[TABLE]

which is a very good approximation as long as the modes are strongly confined in the cavity region 111Rigorously, one should employ the set of Wannier functions to carry out the expansion of Eq. 35, as the cavity modes do not define a complete orthogonal set of functions. Nevertheless, for strongly localized cavity modes this expansion describes very well the physics of a CCW system.. In Eq. 35, $\vec{\Upsilon}_{m}(\mathbf{r})$ is the eigenmode of the $m$-th cavity in units of [m*-3/2*] and normalized as

[TABLE]

with $\epsilon_{c}(\mathbf{r})$ being the dielectric function of the single cavity system. By replacing Eq. 35 into Eq. 33, the non-linear integral can be rewritten as

[TABLE]

where the non-linear overlapping elements $\Pi_{lmhp}$ are given by

[TABLE]

Since the mode $\vec{\Upsilon}_{m}(\mathbf{r})$ is strongly confined in the cavity region, the Kerr shift on cavity $m$ induced by the field of neighbor cavities ($m\pm 1$, $m\pm 2$, $\cdots$), is expected to be very small in comparison to the shift induced by its own field strength, therefore, we can safely neglect the cross-Kerr terms in Eq. 38 and consider only the diagonal one in the sum of Eq. 37, i.e.,

[TABLE]

with

[TABLE]

Furthermore, the coefficient $\tilde{n}_{2}(\mathbf{r})$ is non-zero only within the non-linear material, which allow us to write Eq. 40 as

[TABLE]

where $\tilde{n}_{2}$ is the same of Eq. 34 but with no spatial dependence, $\epsilon$ is the dielectric constant of the non-linear material, and the function $\Theta(\mathbf{r})$ is defined as

[TABLE]

In order to write Eq. 41 in a more intuitive way, it is an excellent approximation to express the integral of $|\vec{\Upsilon}_{m}(\mathbf{r})|^{4}$ in terms of the cavity non-linear mode volume $V_{m}$ as follows

[TABLE]

Therefore

[TABLE]

Using Eqs. 39 and 44 the system of Eq. 29 can be written as

[TABLE]

where we have introduced the frequency detuning of the mode $\mu$ with respect to the cw frequency $\Omega_{0}$

[TABLE]

and explicitly written the resonance condition for the driven term, which excites only the normal mode $\alpha_{0}$ whose eigenfrequency is closest to the laser frequency $\Omega_{0}$. The explicit time dependence of Eq. C can be removed by carrying out the following transformation

[TABLE]

which turns Eq. C into

[TABLE]

Equation C dictates the non-linear dynamics of the slowly-varying normal mode amplitudes in a CCW system with different coupled cavities, and no particular choice on the boundary conditions. This equation can be highly simplified if we consider identical single mode cavities and periodic boundary conditions, in which the expansion coefficients in the cavity basis are analytical Yariv et al. (1999)

[TABLE]

allowing to evaluate the second sum in Eq. C

[TABLE]

where $V_{c}$ represents the cavity non-linear mode volume and the last approximation is valid for large $M$. Physically, the $\delta$ term in Eq. 50 represents the momentum conservation which is naturally expected from the periodic boundary condition assumption. Using this result, Eq. C takes the following form

[TABLE]

Equation 51 can be further simplified by assuming that the strength of the non-linearity mainly depends on the frequency at which the non-linear effect takes place, $\omega_{\alpha}=\omega_{\alpha_{0}}$, and not on the frequencies contributing to the non-linear dynamics. This allows us to set the frequencies in the non-linear term to $\omega_{\alpha_{0}}$ (frequency of the driven mode) as well as the group velocity terms to $v_{\alpha_{0}}$. When the different frequencies are not independent and fulfill the energy conservation $\omega_{\alpha}+\omega_{\eta}=\omega_{\mu}+\omega_{\alpha+\eta-\mu}$, which is our case of interest, such approximation is expected to describe the non-linear phenomena accurately Colman (2015). Equation 51 is thus turned into the following simpler expression

[TABLE]

where $n_{g,\alpha_{0}}=c/v_{\alpha_{0}}$ is the group index at $\omega_{\alpha_{0}}$. If we define the complex gain as

[TABLE]

with

[TABLE]

being the real gain comming from the non-linear Kerr effect, and

[TABLE]

being the imaginary gain, representing the non-linear losses comming from TPA, we finally get the non-linear system of equations for the mode slowly-varying normal mode amplitudes in a CCW, with identical single mode cavities and periodic boundary conditions

[TABLE]

Equation 56 is equivalent to the one previously derived by Chembo and Yu in non-linear ring resonators Chembo and Yu (2010).

Appendix D Stationary solution and stability analysis

In order to simplify the analysis of the stationary solutions of Eq. 56 and their stability we introduce shifted momentum indices $\{\alpha^{\prime},\mu^{\prime},\eta^{\prime}\}\rightarrow\{\alpha,\mu,\eta\}-\alpha_{0}$, where the driven mode has zero index and the side modes are symmetrically distributed around the driven one. In terms of this new indexes Eq. 56 is written as

[TABLE]

D.1 Stationary solution

The stationary solution of Eq. 57 is found by setting the derivative term to zero

[TABLE]

D.2 System below threshold

When the system is below threshold only the $\alpha^{\prime}=0$ mode is excited, then, Eq. 58 becomes

[TABLE]

and we get the very well known cubic relation between $|{\cal F}_{0}|^{2}$ and $|{\cal B}_{0}^{s}|^{2}$ which leads to the hysteresis phenomenon and bistability of the stationary solution

[TABLE]

The bistability boundaries $B_{\pm}$ correspond to the local extrema of Eq. 60, which are the solutions of the quadratic equation $\partial|{\cal F}_{0}|^{2}/\partial|{\cal B}_{0}^{s}|^{2}=0$ on $|{\cal B}_{0}^{s}|^{2}$. Thus, we readily obtain

[TABLE]

Hysteresis exists if the boundaries $B_{\pm}$ are real and positive, implying the following condition on the laser detuning

[TABLE]

where

[TABLE]

and $\rho(\kappa)$ is a function which depends on the material properties only

[TABLE]

with the dimensionless parameter $\kappa$ defined as

[TABLE]

Notice that by setting $g_{0}^{\rm TPA}=0$ in Eqs. 60, 61 and 62, we recover exactly the same expressions for the ring resonator Chembo and Yu (2010).

D.3 System at threshold

The onset of the oscillations of a pair of side modes ${\cal B}_{\pm\alpha^{\prime}}$ can be determined by studying the linear stability around the trivial equilibrium, i.e., ${\cal B}_{\pm\alpha^{\prime}}^{s}=0$. This technique consists in adding a small time-dependent perturbation $\delta{\cal B}_{\pm\alpha^{\prime}}(t)$ to the stationary equilibrium, and the set of parameters leading to the divergence of $\delta{\cal B}_{\pm\alpha^{\prime}}(t)$ defines the threshold for side mode oscillations. Hence, the mode amplitude ${\cal B}_{\alpha^{\prime}}(t)$ is written as

[TABLE]

and replaced in Eq. 58, which gives

[TABLE]

where we have written explicitly the time derivative operator for the sake of clarity. By keeping only terms linear in $\delta{\cal B}$, and using the steady state solution of the central mode in Eq. 59, we get

[TABLE]

For studying the dynamical behavior of the time-dependent perturbations, Eq. 68 can be rewritten in a more convenient form by making the transformations

[TABLE]

with

[TABLE]

Moreover, in order to simplify the stability analysis, we replace the loss rates of the side modes $\gamma_{\pm\alpha^{\prime}}$ by $\gamma_{0}$, which is a good approximation when $\gamma$ does not have strong fluctuations within the interval $[-\alpha^{\prime},+\alpha^{\prime}]$. Equation 68 is finally turned into

[TABLE]

Notice that the coefficients of Eq. 71 are invariant to the change $\alpha^{\prime}\rightarrow-\alpha^{\prime}$ because $\overline{\omega}_{\alpha^{\prime}}=\overline{\omega}_{-\alpha^{\prime}}$, therefore, the dynamics of the side mode perturbations are described by the following two linearized equations

[TABLE]

where

[TABLE]

The stability of the trivial equilibrium is then assessed by analyzing if the perturbation decays to zero, which means that such equilibrium is stable and there are no side mode oscillations, or if the perturbation diverges, meaning that the trivial equilibrium is unstable thus leading to the onset of side mode oscillations. This is carried out by writing

[TABLE]

which together with Eq. D.3 sets the following $2\times 2$ eigenvalue problem

[TABLE]

The solutions for $\lambda$ in Eq. 75 are readily found

[TABLE]

where we have defined

[TABLE]

The divergence of Eq. D.3, i.e., the onset of the side mode oscillations, is clearly obtained when ${\rm Re}\{\lambda_{+}\}>0$, which in terms of the system parameters reads

[TABLE]

The zeros of Eq. 78 determine the boundaries in which the trivial equilibrium becomes unstable. Therefore, by defining $\tilde{B}=|{\cal C}_{0}^{s}|^{2}$, these boundaries are the solutions of the quadratic equation ${\cal S}_{\alpha^{\prime}}(\tilde{B})=0$, i.e.,

[TABLE]

Notice that Eqs. 78 and 79 reduce exactly to the same corresponding expressions reported in Ref. Chembo and Yu (2010) for $g_{0}^{\rm TPA}=0$. The minimum mode power leading to side mode oscillation is determined by the equation

[TABLE]

In terms of the fundamental CCW parameters and the actual momentum index (no shifted), this modal threshold power for comb generation is found to be

[TABLE]

where $Q_{\alpha_{0}}=\omega_{\alpha_{0}}/\gamma_{\alpha_{0}}$ is the quality factor of the driven CCW normal mode $\alpha_{0}$, and $f(\kappa)$ is a function of $\kappa$ only

[TABLE]

From Eqs. 81 and 82 we clearly see that, as far as the threshold is concerned, TPA increases the minimum power required to start comb generation. The corresponding driven amplitude threshold $|{\cal F}_{\alpha_{0}}|_{\rm th}^{2}(\sigma_{\alpha_{0}})$, which depends on the frequency detuning between the external laser and driven mode, is obtained by replacing Eq. 81 on Eq. 60 and recalling that $|{\cal A}|=|{\cal B}|$, thus

[TABLE]

The optimal detuning that leads to the absolute minimal driven amplitude, required by comb generation, is easily found by solving the equation $\partial|{\cal F}_{\alpha_{0}}|_{\rm th}^{2}/\partial\sigma_{\alpha_{0}}=0$, and it is given by

[TABLE]

Equation 84 shows that in presence of TPA, the optimal laser detuning to start frequency combs is red-shifted.

A closer look to Eq. 81 evidences that $\left|{\cal A}_{\alpha_{0}}\right|_{\rm th}^{2}$ depends actually on the effective area of the cavity mode

[TABLE]

where the mode area $A_{\rm eff}=V_{c}/l$ is widely employed in waveguide physics. Moreover, if we consider the limit of no TPA, i.e., $f(\kappa=0)=1$, and the limit of no internal nanostrcuture (homogeneous slab for photonic crystals) where the group index is very close to the actual refractive index of the dielectric $n$, i.e., $\epsilon=n^{2}\backsimeq n_{g,\alpha_{0}}^{2}$ , Eq. 85 becomes

[TABLE]

which is exactly the same threshold expression (in Watts) for the internal mode power in a non-linear ring resonator Chembo (2016).

D.4 Role of dispersion in frequency comb generation

In order to the boundaries of Eq. 79 be real and positive, the discriminant of the square root has to be positive and the term $\left(2\Delta_{\alpha^{\prime}}g_{0}+\gamma_{0}g_{0}^{\rm TPA}\right)$ has to be negative. Both conditions are satisfied if

[TABLE]

where the critical detuning $\sigma_{\rm cr}$ is defined as

[TABLE]

If the dispersion relation $\omega_{\alpha^{\prime}}$ is expanded in Taylor series around $\omega_{0}$ up to third order, we have

[TABLE]

with $\zeta_{1}$, $\zeta_{2}$ and $\zeta_{3}$ representing the group velocity, group velocity dispersion and third order dispersion, respectively. Replacing Eq. 89 in Eq. 70 gives

[TABLE]

turning the condition of Eq. 87 into

[TABLE]

For normal dispersion $\zeta_{2}<0$ and Eq. 91 defines the following condition for $\alpha^{\prime}$

[TABLE]

while for anomalous dispersion $\zeta_{2}>0$ the condition determined by Eq. 91 reads

[TABLE]

Therefore, from Eqs. 92 and 93, we clearly see that if the system is pumped where dispersion is normal, there is an upper bound on the side mode momentum, thus constraining the spanning of the frequency comb. Nevertheless, if the dispersion is anomalous, there is no upper bound on the momentum for side mode oscillations, and the comb generation regime is easily achieved.

Appendix E Lugiato-Lefever equation

The corresponding spatio-temporal equation for the non-linear formalism presented in Sec. C can be obtained by following the derivation of Ref. Chembo and Menyuk (2013). We start from the system of coupled mode equations in Eq. 56, which can be rewritten as

[TABLE]

where we have explicitly considered the delta factor in the non-linear term and used Eqs. 53, 54, 55 and 65 to write $G_{\alpha_{0}}=g_{\alpha_{0}}(1+i\kappa)$. We define the following spatio-temporal slowly-varying envelope along the waveguide direction

[TABLE]

from which we have

[TABLE]

Equation 94 is now replaced in Eq. 96 leading to

[TABLE]

and because

[TABLE]

we readily get from Eq. 98

[TABLE]

where we have simplified the notation $\psi(y,t)\rightarrow\psi$. By writing the mode detuning as

[TABLE]

and expanding the last term in Taylor series around $\omega_{\alpha_{0}}$

[TABLE]

Eq. 100 is turned into

[TABLE]

If we use Eqs. 95 and 97, and assume an overall losses given by the loss rate at the pump frequency, i.e., $\gamma_{\alpha}=\gamma_{\alpha_{0}}$ 222We are ultimately assuming that the loss rates do not depend on the frequency or, equivalently, on the momentum. Such losses may be given by $\gamma_{\alpha_{0}}$ or by any reliable constant value within the Brillouin zone of the CCW., Eq. E takes the following form

[TABLE]

where the group velocity $\zeta_{1}$ and the group velocity dispersion $\zeta_{2}$ have been written explicitly. In order to get rid of the first spatial-derivative term we change to a moving reference frame with velocity $\zeta_{1}$, i.e.,

[TABLE]

which gives

[TABLE]

Replacing Eq. E in Eq. 104, and introducing the normalized detuning $\varsigma=-2\sigma_{\alpha_{0}}/\gamma_{\alpha_{0}}$, we arrive to

[TABLE]

By defining the quantities $\Psi=\sqrt{2g_{\alpha_{0}}/\gamma_{\alpha_{0}}}\psi$, $F=\sqrt{2g_{\alpha_{0}}/\gamma_{\alpha_{0}}}{\cal F}_{\alpha_{0}}$, $\varrho_{m}=-2\zeta_{m}/\gamma_{\alpha_{0}}$ and $\tau=(\gamma_{\alpha_{0}}/2)t^{\prime}$ we finally get the normalized Lugiato-Lefever equation with non-linear losses (TPA) and high-order dispersion

[TABLE]

Appendix F Soliton dynamics

We report in Figs. 5 and 6 the spectral dynamics of soliton formation for the two relevant cases III and IV, respectively, presented in Figs. 2 and 3 in the manuscript. Panels correspond to different times (in ns units) during the transient dynamics and the associated envelope function along the waveguide direction is shown in the inset.