**3. Time-Domain Model for Quadratic Combs**

An alternative description of quadratic comb dynamics can be given in terms of time evolution of the slowly varying intracavity field envelopes. Let us define the envelopes *A*(*z*, *τ*) for the fundamental and *B*(*z*, *τ*) for the second harmonic electric fields in a resonator. Field dynamics can be described by an infinite dimensional map (Ikeda map) for the field amplitudes [50,51], which describes the evolution of cavity fields over the *m*th round trip, along with the boundary condition for the fields at the end of each round trip. The propagation equations for the fields *Am*(*z*, *τ*) and *Bm*(*z*, *τ*) read as

$$\frac{\partial A\_m}{\partial z} = \left[ -\frac{a\_{c1}}{2} - i\frac{k\_1^{\prime\prime}}{2} \frac{\partial^2}{\partial \mathbf{r}^2} \right] A\_m + i\kappa B\_m A\_m^\* e^{-i\Lambda kz},\tag{7}$$

$$\frac{\partial B\_m}{\partial z} = \left[ -\frac{a\_{c2}}{2} - \Delta k' \frac{\partial}{\partial \tau} - i \frac{k\_2''}{2} \frac{\partial^2}{\partial \tau^2} \right] B\_m + i\kappa A\_m^2 e^{i\Delta kz},\tag{8}$$

where *z* ∈ [0, *L*] is the position along the cavity round trip path; *αc*1,2 are propagation losses (hereafter, subscripts 1 and 2 denote fields at *ω*<sup>0</sup> and 2*ω*0, respectively); *k* 1,2 = <sup>d</sup>2*k*/d*ω*2|*ω*0,2*ω*<sup>0</sup> are the group velocity dispersion coefficients; Δ*k* = d*k*/d*ω*|2*ω*<sup>0</sup> − d*k*/d*ω*|*ω*<sup>0</sup> is the corresponding group-velocity mismatch or temporal walk-off. The "fast-time" variable *τ* describes the temporal profiles of the fields in a reference frame moving with the group velocity of light at *ω*0.

For the case of intracavity SHG, the fields at the beginning of the (*m* + 1)th round trip are related to the fields at the end of the previous *m*th round trip according to the following cavity boundary conditions,

$$A\_{m+1}(0,\tau) = \sqrt{1-\theta\_1} \, A\_m(L,\tau) \, e^{-i\delta\_1} + \sqrt{\theta\_1} A\_{\text{in}} \tag{9}$$

$$B\_{m+1}(0,\tau) = \sqrt{1-\theta\_2} \, B\_m(L,\tau) \, e^{-i\delta\_2},\tag{10}$$

where *θ*1,2 are power transmission coefficients at the coupling mirror, *δ*<sup>1</sup> (*ω*<sup>0</sup> − *ωc*1)*t*<sup>R</sup> and *δ*<sup>2</sup> (2*ω*<sup>0</sup> − *ωc*2)*t*<sup>R</sup> are the round trip phase detunings for the fields at *ω*<sup>0</sup> and 2*ω*0, respectively, with *ωc*<sup>1</sup> and *ωc*<sup>2</sup> the frequencies of the respective nearest cavity resonance, and *A*in is the external, constant driving field amplitude. It is worth noting that the Ikeda map of Equations (7)–(10) can describe different nonlinear systems (SHG or OPO, either singly or doubly resonant), by suitably choosing the boundary conditions. For a singly resonant cavity SHG, *θ*<sup>2</sup> = 1, and the SH field resets at the beginning of each round trip, i.e., *Bm*+1(0, *τ*) = 0.

For a relatively high-finesse resonator, the fundamental field evolves slowly during each round trip, and the infinite dimensional map may be averaged over one round trip length *L*. This averaging procedure yields a single mean field equation for the fundamental field amplitude [50],

$$\mathcal{A}\_{\rm R} \frac{\partial A(t, \tau)}{\partial t} = \left[ -\mathfrak{a}\_{1} - i\delta\_{1} - iL\frac{k\_{1}^{\prime\prime}}{2} \frac{\partial^{2}}{\partial \tau^{2}} \right] A - \rho A^{\*} \left[ A^{2}(t, \tau) \otimes I(\tau) \right] + \sqrt{\theta\_{1}} A\_{\rm in} \tag{11}$$

where *t* is a "slow time" variable, linked to the round trip index as *A*(*t* = *mt*R, *τ*) = *Am*(*z* = 0, *<sup>τ</sup>*) [62–65], *<sup>α</sup>*<sup>1</sup> = (*αc*1*<sup>L</sup>* + *<sup>θ</sup>*1)/2, *<sup>ρ</sup>* = (*κL*)2, ⊗ denotes convolution and the nonlinear response function *I*(*τ*) = F−1[ˆ*I*(Ω)], with ˆ*I*(Ω) = (<sup>1</sup> − *<sup>e</sup>*−*ix* − *ix*)/*x*<sup>2</sup> , *x*(Ω) = Δ*k* + *i*ˆ *k*(Ω) *L*, and ˆ *k*(Ω) = −*αc*,2/2 + *i* Δ*k* Ω + (*k* 2/2)Ω<sup>2</sup> . Here, we define the direct and inverse Fourier transform operator as <sup>F</sup> [·] <sup>=</sup> <sup>∞</sup> <sup>−</sup><sup>∞</sup> ·*ei*Ω*<sup>τ</sup>* <sup>d</sup>*<sup>τ</sup>* and <sup>F</sup>−<sup>1</sup> [·] = (2*π*)−<sup>1</sup> <sup>∞</sup> <sup>−</sup><sup>∞</sup> ·*e*−*i*Ω*<sup>τ</sup>* <sup>d</sup>Ω, respectively.

Similarly to the coupled mode equations in frequency domain, also the mean field Equation (11) exhibits an effective cubic nonlinearity, with a noninstantaneous response analogous to the delayed Raman response of cubic nonlinear media and other generalized nonlinear Schrödinger models.

Linear stability analysis of the cw solution of Equation (11) leads to the following expression for the eigenvalues [50],

$$\lambda\_{\pm} = -\left(\mathfrak{a}\_1 + \rho P\_0[\mathfrak{l}(\Omega) + \mathfrak{l}^\*( - \Omega)]\right) \pm \sqrt{|\mathfrak{l}(0)|^2 \rho^2 P\_0^2 - \left(\delta\_1 - \frac{k\_1'' L}{2} \Omega^2 - i\rho P\_0[\mathfrak{l}(\Omega) - \mathfrak{l}^\*( - \Omega)]\right)^2},\tag{12}$$

which, baring the notation, is substantially equivalent to Equation (5). Figure 4a shows the MI gain, [*λ*+] profile as a function of the walk-off parameter Δ*k* . Clearly, there is no MI for zero walk-off, and MI appears for sufficiently large values of walk-off, revealing the fundamental role of group-velocity mismatch for the formation of quadratic optical frequency combs and related dissipative temporal patterns.

**Figure 4.** Modulation instability gain profiles as a function of temporal walk-off. (**a**) Singly resonant cavity SHG. (**b**) Doubly resonant cavity SHG (parameters are normalized according to Ref. [51]). Adapted with permission from [50,51]. Copyrighted by the American Physical Society.

Hansson et al. [52] demonstrated that the general system of coupled mode Equation (6) can be derived from the map of Equations (7)–(10). However, frequency domain coupled mode equations are not exactly equivalent to the time domain mean field Equation (11): the two approaches differ in the way the dispersion is averaged, although they provide almost equal results for the system of Ref. [48].

Theoretical models, in addition to providing useful insight into the physics of quadratic combs, can be a practical tool for simulating the comb dynamics, giving access to information not always available from the experiment. Both the frequency and time domain formalisms here described lend themselves to the numerical simulation of comb dynamics. Coupled mode Equation (6) is in general more time consuming than time domain approaches, unless it can be cast in a way where fast Fourier transform (FFT) algorithms can effectively reduce the computation time [66]. Numerical integration of the Ikeda map or the derived mean-field equation usually relies on split-step Fourier methods [67,68]. According to this method, propagation along each integration step is carried out in two steps. In a first step, the nonlinear and driving terms are propagated by means of a 4th-order Runge–Kutta method. The dispersive and absorption terms are propagated in a second step, where their propagation operator is evaluated in the Fourier domain, using an FFT algorithm. The simulation initiates by assuming a constant amplitude, input driving field *A*in that describes the resonant pump laser. More importantly, in the first step a numerical white-noise background of one photon per mode must be added in order to seed the nonlinear processes which lead to the comb. Whereas the numerical integration of Ikeda map requires a spatial step size smaller than the cavity round trip length, the mean-field equation can be numerically integrated with temporal step sizes of the order of the round trip time, for the benefit of the computation time.

Figure 5 shows two spectra, (a) and (b), and the respective temporal patterns, (c) and (d), obtained by numerically integrating Equation (11). The simulations have been performed using the parameters from Ref. [48], in the case of quasi-phase matched SHG, for a constant input power of 2 and 7 W, respectively, and a small positive detuning. The simulated spectra are in good agreement with the experimental spectra shown in Figure 2b,c. For the moment, we cannot determine the temporal profile corresponding to a comb spectra. Hence, numerical simulations provide insights on the temporal feature of comb dynamics. We notice that the temporal pattern (c) associated to spectrum (a) has a stable periodic structure (also called Turing or roll pattern), which entails a strong phase coupling between the spectral modes, i.e., a mode-locked regime. Instead, the spectrum of Figure 5b, with secondary combs around the primary sidebands, corresponds to an irregular temporal pattern with no evidence of intermodal phase coupling. Moreover, it does not appear to reach a stationary regime. In both cases, the emission is not purely pulsed, as typically occurs for combs generated in femtosecond, mode-locked lasers, but the temporal patterns coexist with a flat background. The coexistence of a temporal pattern with a flat background is frequent for Kerr combs [69], as well as for combs generated in quantum cascade lasers [70,71]. In fact, in femtosecond laser combs the emission of short pulses is due to a particular phase relation between laser mode—i.e., all the modes have equal phases. However, in a wider sense, mode-locking only requires that a stable phase relation holds between all the mode fields. Finally, numerical simulations also reveal a slow drift of the temporal patterns (both at the fundamental and the SH fields) in the reference frame moving with the group velocity of the FF.

**Figure 5.** Numerical simulation of Equation (11), using the parameters of the system in Ref. [48]. (**a**) Input power 2 W, *δ*<sup>1</sup> = 0.001. (**b**) Input power 7 W, *δ*<sup>1</sup> = 0.01. (**c**,**d**) Details of the respective temporal patterns.

When *θ*<sup>2</sup> < 1, the infinite dimensional map of Equations (7)–(10) describes the case of a doubly resonant optical cavity, where also second harmonic fields may resonate. Leo et al. theoretically analyzed this system [51] and derived a couple of two mean-field equations, which accurately model comb dynamics. These equations read, assuming phase-matched SHG,

$$t\_{\rm R} \frac{\partial A}{\partial t} = \left[ -\kappa\_1 - i\delta\_1 - i \frac{k\_1^{\prime\prime} L}{2} \frac{\partial^2}{\partial \tau^2} \right] A + i\kappa LBA^\* + \sqrt{\theta\_1} A\_{\rm in\prime} \tag{13}$$

$$t\_R \frac{\partial B}{\partial t} = \left[ -\kappa\_2 - i\delta\_2 - \Delta k' L \frac{\partial}{\partial \tau} - i \frac{k\_2' L}{2} \frac{\partial^2}{\partial \tau^2} \right] B + i\kappa L A^2 \,,\tag{14}$$

where *α*<sup>2</sup> is the cavity loss of the SH field.

Under realistic conditions, the two mean-field Equations (13) and (14) can be reduced to a single mean-field equation, analogously to Equation (11) for singly resonant cavity SHG. One obtains

$$t\_R \frac{\partial A}{\partial t} = \left[ -a\_1 - i\delta\_1 - i \frac{k\_1'' L}{2} \frac{\partial^2}{\partial \tau^2} \right] A - \rho A^\* \left[ A^2 \otimes I \right] + \sqrt{\theta\_1} A\_{\text{in} \nu} \tag{15}$$

where the Fourier transform of the kernel function *J* is

$$f(\Omega) = \frac{1}{n\_{2} + i\delta\_{2} - i\Delta k'L\Omega - i\frac{k\_{2}''L}{2}\Omega^{2}}.\tag{16}$$

A linear stability analysis of the cw solution (for both the Ikeda map and the mean-field approximations) reveals the significant role of temporal walk-off in enabling comb formation. However, in this case, MI gain may also occur for zero or relatively small values of the walk-off (Figure 4b).

#### **4. Combs in Optical Parametric Oscillators**

Degenerate optical parametric oscillation is the inverse process of cavity SHG, when the pump field *A*in at the FF *ω*<sup>0</sup> is replaced by a pump field *B*in at the SH frequency 2*ω*0. Its dynamics can be described by an infinite dimensional map as well, where, in addition to Equations (7) and (8), the following boundary conditions hold for the fields at the beginning of each round trip,

$$A\_{m+1}(0,\tau) = \sqrt{1-\theta\_1} \, A\_m(L,\tau) \, e^{-i\delta\_1} \tag{17}$$

$$B\_{m+1}(0,\tau) = B\_{\text{in}}.\tag{18}$$

Here, we consider an OPO cavity where only the parametric field resonates. It is straightforward to extend the analysis to the case when the harmonic pump field also resonates.

Following the approach of Ref. [50], the infinite dimensional map can be combined into a single mean-field equation for the parametric field *A*, which reads, assuming Δ*k* = 0 [54],

$$t\_{\mathbb{R}}\frac{\partial A(t,\tau)}{\partial t} = \left[-a\_1 - i\delta\_1 - i\frac{Lk\_1''}{2}\frac{\partial^2}{\partial \tau^2}\right]A(t,\tau) - \mu^2 A^\*(t,\tau)\left[A^2(t,\tau)\otimes I(\tau)\right] + i\mu B\_{\text{in}}A^\*(t,\tau), \tag{19}$$

where all the physical parameters and the kernel function *I* are the same as in Equation (11). We note that Equation (19) is similar to the corresponding mean-field equation for comb dynamics in cavity SHG, except for the parametric driving force (last term on the r.h.s.). Equation (19) has a trivial zero solution, *<sup>A</sup>*<sup>0</sup> <sup>=</sup> 0, and a nontrivial time independent solution, *<sup>A</sup>*<sup>0</sup> <sup>=</sup> <sup>|</sup>*A*0|*eiφ*. From a linear stability analysis of the constant solution, we derived the following expression for the eigenvalues [54],

$$\lambda\_{\pm} = -\left[a\_1 + \mu^2 |A\_0|^2 \mathcal{Z}\_+(\Omega)\right] \pm \sqrt{(a\_1^2 + \delta\_1^2) - \left[\delta\_1 - D\_2 \Omega^2 - i\mu^2 |A\_0|^2 \mathcal{Z}\_-(\Omega)\right]^2},\tag{20}$$

where |*A*0| <sup>2</sup> = [−*α*<sup>1</sup> ± *μ*2*B*<sup>2</sup> in − *<sup>δ</sup>*<sup>2</sup> <sup>1</sup> ]/*μ*<sup>2</sup> <sup>ˆ</sup>*I*(0) is the squared modulus of the nontrivial solution and I±(Ω) = <sup>ˆ</sup>*I*(Ω) <sup>±</sup> <sup>ˆ</sup>*I*∗(−Ω). Similarly, for the zero solution the eigenvalues of the linearized system are

$$
\lambda \pm = -\pi\_1 \pm \sqrt{\mu^2 B\_{\text{in}}^2 - \left(\delta\_1 - D\_2 \Omega^2\right)^2}.\tag{21}
$$

Both solutions exhibit MI gain for Re[*λ*+] > 0, which is shown in Figure 6a,b as a function of the cavity detuning. From Equation (20) it clearly appears that MI gain for the nontrivial solution depends both on walk-off Δ*k* , through I±(Ω), and GVD. As for singly resonant cavity SHG, MI only manifests itself for relatively high walk-off values, while it is absent for zero walk-off, as shown in Figure 6c. The instability of the zero solution, which is not expected in the usual dispersionless analysis of the OPO, does not depend on the walk-off, but it is rather induced by GVD. Actually, GVD is responsible for the unequal spacing between cavity resonances, so that they are asymmetrically displaced with respect to the degeneracy frequency *ω*0, when the latter is perfectly resonant. Thus, GVD effectively favors parametric oscillations close to the degeneracy frequency. For normal dispersion, a positive detuning between the degeneracy frequency and the nearest cavity resonance can make symmetric an initially asymmetric pair of distant resonances, which now can more favourably oscillate than the degeneracy frequency *ω*0. The larger the detuning, the more distant the symmetric resonances are. For small negative detunings no resonance pair can be symmetrically displaced around *ω*0, and MI gain is maximum at the degeneracy frequency, decreasing as a function of the detuning amplitude. The same occurs in the case of anomalous dispersion, provided that the detuning sign is reversed.

**Figure 6.** Optical frequency combs (OFC) in a degenerate OPO. (**a**,**b**) show the MI gain as a function of the normalized cavity detuning Δ = *δ*1/*α*1, for the constant solution and the zero solution, respectively. (**c**) MI gain profiles as a function of the temporal walk-off. Adapted with permission from [54]. Copyrighted by the American Physical Society.

Frequency comb generation in an OPO has been demonstrated by using a nearly degenerate OPO pumped by a frequency doubled cw Nd:YAG laser (Figure 7). The OPO was based on a 15-mm-long periodically-poled 5%-MgO-doped lithium niobate crystal, with a grating period of Λ = 6.92 μm, enclosed in a bow-tie cavity resonating for the parametric wavelengths around 1064 nm, similar to that used for cavity SHG. The nonlinear crystal was located between two high-reflectivity spherical mirrors (with radius of curvature = 100 mm), while a flat high-reflectivity mirror was mounted on a piezoelectric actuator for cavity length control. A fourth, partially reflective flat mirror (R = 98%) allowed us to couple out the generated parametric radiation. The SH beam entered the OPO cavity from a first spherical mirror, passed through the nonlinear crystal, and left the cavity at the second spherical mirror. The FSR of the cavity was 505 MHz. We observed combs for pump powers higher than 85 mW (about three times the OPO threshold of 30 mW) and studied the effect of small cavity detunings on the comb spectra. Figure 8a–c show the experimental comb spectra recorded for Δ = −0.30, 0.00, 0.30, respectively, with 300 mW of pump power. We found a good agreement with the corresponding spectra, shown in Figure 8d–f, calculated by numerically integrating the mean-field Equation (19). Experimental spectra for negative and zero detunings are very similar, displaying 1 FSR line spacing, whereas for the positive detuning the experimental spectrum consists of two pairs of widely spaced symmetric lines.

**Figure 7.** OFC in a degenerate OPO. Scheme of the experimental setup: beam splitter (BS), electro-optic phase modulator (EOM), periodically poled lithium niobate crystal (PPLN), piezoelectric actuator (PZT), photodiode (PD). Adapted with permission from [54]. Copyrighted by the American Physical Society.

**Figure 8.** (**a**–**c**) Experimental OPO optical spectra for detunings Δ = −0.30, 0.00, 0.30, respectively. (**d**–**f**) Corresponding numerically calculated spectra. From [54]. Copyrighted by the American Physical Society.

## **5. Single Envelope Equation**

Models based on the two field envelopes, i.e., Equations (7)–(10) and their approximations hold as long as there is a single dominant nonlinear process and the combs are confined around two carrier frequencies. When the combs start to overlap, or multiple nonlinear processes play a prominent role, frequency comb generation may be studied by means of a more general model, based on a single-envelope equation combined with the boundary conditions that relate the fields between successive round trips and the input pump field [55],

$$\mathcal{F}[A^{m+1}(\tau,0)] = \sqrt{\hat{\theta}(\Omega)}\mathcal{F}[A\_{\text{in}}] + \sqrt{1-\hat{\theta}(\Omega)}\,\epsilon^{i\phi\_0}\mathcal{F}[A^m(\tau,L)] \tag{22}$$

$$\left[\partial\_{\overline{z}} - D\left(i\frac{\partial}{\partial \tau}\right) + \frac{a\_d}{2}\right] A^m(\tau, z) = i\rho\_0 \left(1 + i\tau\_{\text{sh}}\frac{\partial}{\partial \tau}\right) p\_{\text{NL}}(\tau, z, A^m) \,. \tag{23}$$

The boundary condition, Equation (22), is written in the Fourier domain, in order to account for the frequency dependence of the transmission coefficient *θ* at the input port of the resonator. It determines the intra-cavity field *Am*+1(*τ*, *z* = 0) at the beginning of (*m* + 1)th round trip in terms of the field at the end of the previous round trip *Am*(*τ*, *z* = *L*) and the pump field *A*in. Equation (23) is written in a reference frame moving at the group velocity at *ω*0: *p*NL is the broadband envelope of the nonlinear polarization *<sup>P</sup>*NL <sup>=</sup> *<sup>P</sup>*(2) NL <sup>+</sup> *<sup>P</sup>*(3) NL <sup>+</sup> ... <sup>=</sup> (*χ*(2)*E*<sup>2</sup> <sup>+</sup> *<sup>χ</sup>*(3)*E*<sup>3</sup> <sup>+</sup> ...); *<sup>ρ</sup>*<sup>0</sup> <sup>=</sup> *<sup>ω</sup>*0/2*n*0*c*0; *<sup>τ</sup>*sh <sup>=</sup> 1/*ω*<sup>0</sup> is the shock coefficient that describes the frequency dependence of the nonlinearity, and *α<sup>d</sup>* is the distributed linear loss coefficient. Dispersion to all orders is included by the operator *D*,

$$D\left(i\frac{\partial}{\partial \tau}\right) = \sum\_{l\geq 2} i\frac{\beta\_l}{l!} \left(i\frac{\partial}{\partial \tau}\right)^l,\tag{24}$$

where *β<sup>l</sup>* = (*d<sup>l</sup> β*/*dω<sup>l</sup>* )*ω*=*ω*<sup>0</sup> are expansion coefficients of the propagation constant *β*(*ω*).

Figure 9 shows a spectrum obtained from the numerical simulation of Equations (22) and (23), when SHG and nondegenerate optical parametric oscillation are simultaneously quasi-phase matched in a radially poled, lithium niobate microresonator, pumped at 1850 nm (162 THz). In this case, the quasi-phase matching period for SHG (25.56 μm) simultaneously quasi-phase matches a nondegenerate OPO with idler (signal) at 56 THz (106 THz). The broadband power spectral density shows a generation of a multi-comb array, extending from the mid-infrared into the ultraviolet with a spacing of a single FSR (around 92 GHz). In addition to combs at the FF, SH, and third-harmonic (TH), two additional combs are generated around signal and idler frequencies. Moreover, several secondary combs appear between the FF and the SH and between the SH and the TH, respectively. These combs are generated by sum-frequency generation and difference frequency generation processes. For instance, the comb SC1 centered at 218 THz results from SFG between the idler and the FF, while SC3 (around 380 THz) results from SFG between the idler and the SH. On the other hand, DFG between the SH (TH) and the idler leads to a secondary comb SC2 (SC4) centered at 268 THz (430 THz).

**Figure 9.** Numerical simulation of the single-envelope map when SHG and OPO processes are simultaneously phase matched in a lithium niobate microresonator pumped by 100 mW of cw power at 1850 nm. Reprinted with permission from [55] c The Optical Society.

#### **6. Perspectives**

Because of the intrinsically higher strength of the quadratic nonlinearity with respect to the third-order one, quadratic comb generation can be less demanding in terms of power density and cavity quality factor. Although quadratic combs have been generated in bulk cavities with moderate pump powers, their performance could increase if implemented in miniaturized devices, thus further extending and stimulating new applications [72,73]. Scaling the resonator to micrometric size may be beneficial for quadratic combs, allowing for a dramatic reduction of threshold power and a flexible management of the dispersion through a geometric design, allowing for a broader comb emission. As a matter of fact, direct generation of quadratic frequency combs has been very recently observed in chip-scale lithium niobate devices, such as periodically poled linear waveguide resonators [74,75], or exploiting naturally phase-matched SHG in whispering-gallery-mode resonators [76,77].

Several materials with second-order nonlinearity are suitable to be shaped into low-loss small-footprint resonators. Most of them have been used to generate Kerr combs [78–82], and, in some cases, secondary quadratic effects have been reported [78,79] or explicitly considered [83]. In contrast to Kerr combs, quadratic combs usually require more stringent conditions on phase matching and group velocity mismatch between different spectral components. Natural [84,85], cyclic [86], and quasi-[87,88] phase matching have been used in crystalline whispering-gallery-mode resonators. More recently, significant progress has been made in the fabrication of integrated, high-Q, lithium niobate microresonators for *χ*(2) processes [89–92]. III-V materials provide an interesting photonic platform for second-order nonlinear optics, and different techniques have been devised to

achieve phase matching [93], in particular for resonant structures [94–96]. It is worth noting that the well developed silicon platform can also be exploited for second-order nonlinear interaction. In fact, Timurdogan et al. demonstrated that a large "dressed" *χ*(2) nonlinearity can be induced by breaking the crystalline center-symmetry of silicon when a direct-current field is applied across p-i-n junctions in ridge waveguides [97], enabling the implementation of quasi-phase matching schemes.

Unlike optical frequency combs in mode-locked lasers, parametrically generated combs do not usually correspond to a stable pulsed emission in the time domain. Different temporal regimes are possible, from chaotic to perfectly coherent states. The formation of temporal cavity solitons in a cw-pumped nonlinear resonator has attracted a particular interest in connection with parametrically generated combs [98]. Combs associated to a cavity soliton are broadband and highly coherent, which makes them ideal for low noise and metrological applications. In fact, cavity solitons are robust states which circulate indefinitely in a cavity, thanks to the double compensation between nonlinearity and chromatic dispersion and between cavity losses and cw driving. Recent theoretical works aim at identifying the dynamical regimes that exhibit soliton states or localized solutions in cavity SHG systems [99–102] or OPOs [103,104].

Finally, optical frequency combs are attracting a growing interest as sources of complex quantum states of light for high-dimensional quantum computation [26,105,106]. Second-order nonlinear optical systems are efficiently used for generation of quantum states of light: the classical correlations that establish in three-wave-mixing processes hold at the quantum level as well, leading, for instance, to generation of squeezed light or bipartite entanglement in an OPO. Tripartite, or quadripartite multicolor entanglement has been predicted in second-order nonlinear devices [107], in particular when multiple cascaded second-order nonlinear interactions occur, in traveling-wave or intracavity processes [108–110]. Interestingly, a recent study based on the three-wave model of Equations (1)–(3) predicts five-partite entanglement between one-octave-distant modes [111]. This result suggests that quadratic combs could exhibit multipartite entanglement between frequency modes, which are essential for scalable measurement-based quantum computing [112]. To fully explore these features, a general and complete analysis of the quantum dynamics of quadratic combs is needed [113].

**Author Contributions:** Designed and performed the experiments, I.R., S.M., M.P., P.M., and M.D.R.; theoretical modeling I.R, F.L., T.H., M.E., S.W., and M.D.R.; writing–original draft preparation, I.R. and M.D.R.; funding acquisition, M.E., P.D.N., S.W., and M.D.R.; all authors analyzed the data, discussed the results, read and edited the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR), PRIN 2015KEZNYM (NEMO); Ministero degli Affari Esteri e della Cooperazione Internazionale, project NOICE Joint Laboratory; European Union's Horizon 2020 research and innovation programme (Qombs Project, FET Flagship on Quantum Technologies grant no. 820419); the Rutherford Discovery Fellowships of the Royal Society of New Zealand and the Marsden Fund of the Royal Society of New Zealand. The work of S. W. is supported by the Ministry of Education and Science of the Russian Federation (Minobrnauka) (14.Y26.31.0017). T.H. acknowledges funding from the Swedish Research Council (Grant No. 2017-05309).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
