**2. Intracavity Second Harmonic Generation**

The first system that we investigated for the generation of quadratic OFCs was a cw-pumped, cavity enhanced SHG system. The system, shown in Figure 1a, was based on a 15-mm-long periodically poled LiNbO3 crystal, placed inside a traveling-wave optical resonator (free spectral range FSR = 493 MHz, quality factor *Q* = 108), resonating at the fundamental laser frequency *ω*0. Mirror reflectivities were chosen in order to facilitate the onset of an internally pumped OPO [48]. The crystal was pumped by a narrow-line, 1064-nm-wavelength Nd:YAG laser, amplified by a Yb-doped fiber amplifier. Frequency locking of a cavity resonance to the laser was achieved by the Pound–Drever–Hall technique [56].

**Figure 1.** Singly resonant cavity second harmonic generation (SHG). (**a**) Experimental setup: periodically poled lithium niobate crystal (PPLN), piezoelectric actuator (PZT), photodiode (PD), dichroic mirror (DM). The output beams are detected and processed by radio-frequency (RF) analyzers, while optical spectral analysis is performed by an optical spectrum analyzer in the infrared range and a confocal Fabry-Pèrot interferometer (CFP) in the visible range. (**b**) Schematic representation of the first steps leading to the formation of a dual optical frequency comb in cavity-enhanced second-harmonic generation: (left) second-harmonic generation with cascaded nondegenerate optical parametric oscillator (OPO) gives rise to two subharmonic sidebands, which in turn (right) lead to successive, multiple second-harmonic and sum-frequency generations. Adapted with permission from [48]. Copyrighted by the American Physical Society.

The phase-matching condition for SHG was achieved by properly adjusting the crystal temperature. Under this condition, we observed a first regime of pure harmonic generation, where the harmonic power increased with the input pump power. As shown in Figure 2a, when the input power exceeded the threshold for internally pumped OPO, the second harmonic power ceased to grow, and two parametric waves started to oscillate at frequencies *ω*<sup>0</sup> ± Δ*ω*, symmetrically placed around the fundamental frequency (FF). As the power was further increased, additional sidebands appeared, displaced by multiples of Δ*ω*, leading to a multiple-FSR-spaced frequency comb, as sketched in Figure 2b. Finally, when the input power exceeded 5 W, secondary combs appeared around each of the primary comb lines, shown in Figure 2c. These secondary combs were spaced by 1 cavity FSR, as confirmed by the intermodal beat notes detected by fast photodetectors, both in the IR and in the visible spectral regions and processed by a radio frequency (RF) spectrum analyzer.

**Figure 2.** Optical spectral power around the fundamental mode for (**a**) 170 mW, (**b**) 2 W, and (**c**) 9 W of input powers. Adapted with permission from [48]. Copyrighted by the American Physical Society.

Subsequently, wave vector mismatch Δ*k* was changed to finite values by varying the crystal temperature. Figure 3 shows infrared spectra observed for different values of the mismatch vector. For a positive mismatch, Δ*k* > 0, the spectra (a)–(d) show widely separated sidebands, similar to the spectra observed at Δ*k* = 0 (see Figure 2b). The spacing between sidebands, as well as the pump power threshold for cascaded optical parametric oscillation, rapidly increases with the mismatch. For Δ*k* < 0, the spectra (e)–(h) consist of closely spaced (1 FSR) comb lines, and the spectral bandwidth increases with the magnitude of the mismatch. Larger negative phase mismatches are precluded by the limited accessible temperature range. Figure 3i,j show the beat notes corresponding to the comb in Figure 2c and the comb in Figure 3g, respectively. The broad feature of the beat note (i) reveals a strong intermodal phase noise and, as a consequence, a low degree of coherence between the comb teeth. This feature is consistent with a scenario where comb modes are weakly coupled with each other, as they originate independently from each other. On the contrary, the beat note (j) is extremely

narrow, being limited by the detection resolution bandwidth and indicates a low intermodal phase noise and thus a strong phase coupling between all the comb teeth.

It is worth noting that the nonlinear resonator exhibits a noticeable thermal effect, mainly due to light absorption in the nonlinear crystal, which generates heat and leads to an increase of the cavity optical path, via thermal expansion and thermo-optic effect [57]. The photothermal effect introduces an additional nonlinear dynamical mechanism, with a temporal scale determined by the thermal diffusion time over the typical optical beam size [58]. In our case, the photothermal effect was helpful in thermally locking a cavity resonance to the laser frequency [59] when, especially at higher power, the PDH locking scheme was less effective. However, a better comprehension of the effect of thermal dynamics on comb formation requires further investigations.

**Figure 3.** Optical spectra for phase-mismatched singly resonant cavity SHG. (**a**–**d**) Positive phase mismatch; (**e**–**h**) negative phase mismatch. Intermodal beat notes corresponding to (**i**) comb spectrum of Figure 2c, (**j**) comb spectrum in panel (**g**).

As anticipated in the introduction, the onset of internally pumped OPO marks the beginning of a cascade of second-order nonlinear processes, which eventually produces a comb of equally spaced frequencies. As depicted in Figure 1b, once generated, each parametric mode can generate new field modes through second harmonic, (*ω* + Δ*ω*)+(*ω* + Δ*ω*) → 2*ω* + 2Δ*ω* and sum frequency with the fundamental wave *ω* + (*ω* + Δ*ω*) → 2*ω* + Δ*ω*, processes, respectively. All these processes have been considered for the derivation of a simple system of coupled mode equations for the three intracavity subharmonic electric field amplitudes, the fundamental *A*0, at *ω*0, and the parametric intracavity fields *A<sup>μ</sup>* and *Aμ*¯, at *ωμ* = *ω*<sup>0</sup> + Δ*ω* and *ωμ*¯ = *ω*<sup>0</sup> − Δ*ω*, respectively, which read [48]

$$A\_0 = -\left(\gamma\_0 + i\Delta\_0\right)A\_0 - 2\g\,\eta\_{00\mu\emptyset}A\_0^\*A\_\mu A\_\mu - \left(\eta\_{0000}|A\_0|^2 + 2\eta\_{0\mu 0\mu}|A\_\mu|^2 + 2\eta\_{0\mu 0\mu}|A\_\mu|^2\right)A\_0 + F\_{\text{in}} \tag{1}$$

$$A\_{\mu} = -\left(\gamma\_{\mu} + i\Lambda\_{\mu}\right)A\_{\mu} - \lg \eta\_{\mu\emptyset0} \left.A\_0^2 A\_{\mu}^\* - \lg \left(2\eta\_{\mu 00\mu}|A\_0|^2 + \eta\_{\mu\mu\mu\mu}|A\_{\mu}|^2 + 2\eta\_{\mu\emptyset\mu\emptyset}|A\_{\mu}|^2\right)A\_{\mu} \tag{2}$$

$$\dot{A}\_{\varPhi} = -\left(\gamma\mu + i\Lambda\mu\right)A\_{\varPhi} - \lg\eta\mu\alpha\alpha\ A\_0^2 A\_{\varmu}^\* - \lg\left(2\eta\mu\alpha\mu\left|A\alpha\right|^2 + 2\eta\mu\mu\mu\left|A\_{\varmu}\right|^2 + \eta\mu\mu\mu\left|A\_{\varmu}\right|^2\right)A\_{\varmu}.\tag{3}$$

Here, *F*in = -2*γ*0/*t*R*A*in is the cavity coupled amplitude of the constant input driving field *A*in, at frequency *ω*0; the *γ*'s are the cavity decay constants; the Δ's are the cavity detunings of the respective modes; *g* = (*κL*)2/2*t*<sup>R</sup> is a gain factor depending on the crystal length *L* (hereafter we consider the cavity length equal to the crystal length); *<sup>t</sup>*<sup>R</sup> is the cavity round trip time; *<sup>κ</sup>* <sup>=</sup> <sup>√</sup>8*ω*0*χ*(2) eff / *c*3*n*<sup>2</sup> <sup>1</sup>*n*2<sup>0</sup> is the second-order coupling strength. The latter is normalized so that the square modulus of the field amplitudes is measured in watts, with *χ*(2) eff the effective second-order susceptibility, *c* the speed of light, *n*1,2 the refractive indices, and <sup>0</sup> the vacuum permittivity. The integer mode number *μ* denotes the *μ*th cavity mode, starting from the central mode at *ω*0, and overline stands for negative (lower frequencies). The *η*'s are complex nonlinear coupling constants, depending on the wave-vector mismatches associated with a pair of cascaded second-order processes,

$$\eta\_{\mu\nu\rho\upsilon} = \frac{2}{L^2} \int\_0^L \int\_0^z \exp\left[-i(\zeta\_{\mu\nu}z - \zeta\_{\rho\nu}z')\right] \,\mathrm{d}z' \,\mathrm{d}z \tag{4}$$

where *ξjk* = *kω<sup>j</sup>* + *kω<sup>k</sup>* − *kωj*+*ω<sup>k</sup>* .

A linear stability analysis of Equations (1)–(3) predicts the conditions for which a *μ*-pair of parametric fields starts to oscillate. By calculating the eigenvalues corresponding to Equations (1)–(3) linearized around the cw steady state solution, one obtains [49]

$$\lambda \pm = -\gamma - \lg(\eta\_{\mu 00\mu} + \eta\_{\mu 00\bar{\mu}}^{\*})|A\_0|^2 \pm \sqrt{g^2 |\eta\_{\mu \bar{\mu} 00}|^2 |A\_0|^4 - \left[\Delta\_0 - D\_2 \mu^2 - ig(\eta\_{\mu 00\mu} - \eta\_{\mu 00\bar{\mu}}^{\*})|A\_0|^2\right]^2},\tag{5}$$

where *<sup>D</sup>*<sup>2</sup> −2*π*2*c*3*β*/*L*<sup>2</sup> *<sup>n</sup>*<sup>3</sup> <sup>0</sup> = −(*c*/2*n*0)*D*<sup>2</sup> <sup>1</sup>*β* accounts for the group velocity dispersion at *ω*0, with *β* = *<sup>d</sup>*2*<sup>k</sup> dω*<sup>2</sup> *ω*0 , and *n*<sup>0</sup> = *n*(*ω*0) the refractive index at *ω*0. Side modes start to oscillate, i.e., the zero solution for the parametric fields becomes unstable when the real part of an eigenvalue goes from negative to positive values. The coupling constants which appear in Equation (5) are: *ημμ*¯00, which is the parametric gain related to cascaded SHG and OPO, whereby two photons at frequency *ω*<sup>0</sup> annihilate and two parametric photons at *ωμ* and *ωμ*¯ are created, mediated by a SH photon; and *ημ*00*<sup>μ</sup>* (*ημ*¯00*<sup>μ</sup>*¯), which is related to the sum frequency process between a parametric photon at *ωμ* (*ωμ*¯) and the pump. The latter process is the most relevant nonlinear loss at the threshold (second term of r.h.s of Equation (5)), and provides a nonlinear phase shift (last term in the square brackets of r.h.s of Equation (5)). The lowest threshold occurs for a pair of parametric fields which starts to grow close to the minima of the sum frequency generation (SFG) efficiency.

A general expression for the dynamic equations for any number of interacting fields can be derived heuristically [49], yielding for each field *Aμ*, nearly resonant with the *μ*-th cavity mode,

$$\dot{A}\_{\mu} = -\left(\gamma\_{\mu} + i\Delta\_{\mu}\right)A\_{\mu} - \underset{\rho,\sigma}{\text{g}} \sum\_{\substack{\rho,\sigma\\\nu=\rho+\sigma-\mu}} \eta\_{\mu\nu\rho\sigma} A\_{\nu}^{\*} A\_{\rho} A\_{\sigma} + F\_{\text{in }\prime} \tag{6}$$

where the summation over the indices *ρ* and *σ* goes over all the cavity resonant modes. The complex coupling constants are given by Equation (4), while the constraint over *ν* accounts for energy conservation. The coupled mode Equation (6) is formally analogous to the modal expansion for Kerr combs [60,61] and describes the whole comb dynamics. It is worth noting that the information provided by the linear stability analysis only holds for the very beginning of comb formation. Very quickly, a large number of cavity modes under the gain curve grow from noise. At the same time, they interact with each other through multiple nonlinear processes. These processes are not considered in the linear stability analysis, which intrinsically considers only three interacting modes. The long-term spectral configuration is thus the result of a complex interaction between many modes, over thousands of cavity round trips [52].
