Optical Frequency Combs in Quadratically Nonlinear Resonators

Abstract

Optical frequency combs are one of the most remarkable inventions in recent decades. Originally conceived as the spectral counterpart of the train of short pulses emitted by mode-locked lasers, frequency combs have also been subsequently generated in continuously pumped microresonators, through third-order parametric processes. Quite recently, direct generation of optical frequency combs has been demonstrated in continuous-wave laser-pumped optical resonators with a second-order nonlinear medium inside. Here, we present a concise introduction to such quadratic combs and the physical mechanism that underlies their formation. We mainly review our recent experimental and theoretical work on formation and dynamics of quadratic frequency combs. We experimentally demonstrated comb generation in two configurations: a cavity for second harmonic generation, where combs are generated both around the pump frequency and its second harmonic and a degenerate optical parametric oscillator, where combs are generated around the pump frequency and its subharmonic. The experiments have been supported by a thorough theoretical analysis, aimed at modelling the dynamics of quadratic combs, both in frequency and time domains, providing useful insights into the physics of this new class of optical frequency comb synthesizers. Quadratic combs establish a new class of efficient frequency comb synthesizers, with unique features, which could enable straightforward access to new spectral regions and stimulate novel applications.

1. Introduction

Twenty years ago, optical frequency combs (OFCs) were established as powerful tools for accurate measurements of optical frequencies and timekeeping [1,2], a result of a long-standing effort, which was recognized with the Nobel Prize in Physics in 2005 [3,4]. The discrete ensemble of equally spaced laser frequencies that distinguish OFCs from other light sources is the spectral counterpart of the regular train of short pulses emitted by mode-locked lasers, which were initially used for comb generation. OFCs have become a critical component in many scientific and technological applications [5], from highly accurate optical frequency measurements for fundamental tests of physics [6,7,8,9] to exoplanet exploration [10,11,12] from air pollution detection [13,14,15,16,17,18] to telecommunication systems [19,20,21], while a growing interest has arisen in the quantum properties of OFCs [22,23,24,25,26].

Thereafter, comb emission was also demonstrated in continuous-wave (cw) laser-pumped resonators through cascaded third-order ${\chi }^{\left(3\right)}$ parametric processes [27]. In such Kerr resonators, a first pair of sidebands is generated around the pump frequency by cavity modulation instability or degenerate four-wave mixing (FWM); subsequently, cascaded four-wave mixing processes lead to the formation, around the pump frequency, of a uniform frequency comb, where self- and cross-phase modulation act to compensate for the unequal cavity mode spacing induced by the group velocity dispersion (GVD) [28,29]. Because of the relatively low strength of third-order nonlinearity, generation of Kerr combs requires small interaction volumes and high-Q resonators. For these reasons, small size resonators are particularly suited to reach broadband comb generation with quite moderate pump power [30]. Moreover, when the mode size is comparable with the light wavelength, a careful design of the resonator geometry can effectively modify the GVD of the resonator, leading to a broader spectral emission.

While ${\chi }^{\left(2\right)}$ three-wave mixing processes have been widely used for spectral conversion of femtosecond laser combs since their inception [31,32,33,34,35,36], only in recent years it was demonstrated that quadratic ${\chi }^{\left(2\right)}$ processes can lead to direct generation of optical frequency combs in cw-pumped quadratic nonlinear resonators. Actually, in 1999 Diddams at al. generated an OFC in a second-order nonlinear system, by actively inducing intracavity phase modulation inside a cw-pumped nearly degenerate optical parametric oscillator (OPO) [37], following a long development of phase modulation in lithium niobate for comb generation [38]. According to this scheme, besides the nonlinear crystal for parametric amplification, a phase modulator was placed inside the OPO cavity and driven at a modulation frequency equal to the cavity free spectral range. The modulator thus generated a family of phase-coupled sidebands, around the nearly degenerate signal and idler waves, which coincided with the resonator mode frequencies. Unlike other works presented in the following, where combs arise through purely ${\chi }^{\left(2\right)}$ optical processes, in that work combs were initially seeded by the sidebands generated in the intracavity modulator. Optical parametric amplification further increased the number of resonant sidebands, eventually leading to an 18-nm wide comb of equally spaced, mode-locked lines around the degenerate OPO frequency, only limited by the dispersive shift of the cavity modes, where mode-locking is imposed by phase modulation.

More recently, an optical frequency comb was produced by adding a second nonlinear crystal in a nondegenerate OPO [39]. The authors observed comb formation around the signal wavelength when the second crystal was phase mismatched for second harmonic generation (SHG) of the signal wave. Subsequent investigations of the same system reported experimental evidence of a comb around the second harmonic of the signal wave, whereas the comb around the signal was simultaneously transferred to the idler spectral range by parametric amplification [40]. In this case, the phase mismatched crystal behaves like a Kerr medium, producing a phase shift of the fundamental wave, which is proportional to the field intensity [41,42,43]. This phase shift can be explained as the consequence of cascaded quadratic processes which occur in the crystal when SHG is not phase matched. Indeed, when the fundamental pump wave, at frequency $\omega /2\pi$ enters a nonlinear crystal, a second harmonic field is generated, $\omega +\omega \to 2\omega$. If the process is not phase matched, the second harmonic (SH) field travels at a different phase velocity and, after half a coherence length, down-converts back to the fundamental frequency, $2\omega -\omega \to \omega$, with a different phase from that of the unconverted pump field.

As we will see later, a different cascaded three-wave-mixing process is decisive for the onset of frequency combs in phase-matched intracavity SHG—namely, internally pumped optical parametric oscillation [44,45]. In fact, degenerate optical parametric oscillation and SHG are mutually inverse processes, which satisfy the same phase matching condition, $\Delta k=2k_1-k_2=0$, between wave vectors $k_1=k\left(\omega \right)$ of the fundamental field and $k_2=k\left(2\omega \right)$ of the second harmonic field, respectively. Therefore, a properly phase-matched crystal placed inside an optical resonator can work either for SHG or parametric oscillation, depending on whether it is pumped at the fundamental or second harmonic frequency, respectively. However, the harmonic field generated in the first case can act as a pump for a nondegenerate cascaded OPO, and a pair of parametric fields start to oscillate with frequencies symmetrically placed around the fundamental pump. Although internally pumped OPO was observed and investigated for a long time, before the importance of OFCs was established [44,45,46,47], the observation of frequency combs in quadratic nonlinear media was postponed to recent years.

Here, we present a concise introduction to the physical mechanism that underlies quadratic comb formation, as well as an extended theoretical framework that has been developed so far. We particularly focus on our recent activity in this field, discussing our experimental and theoretical work on direct generation of quadratic combs. As a whole, it represents a systematic and coherent, although not exhaustive, approach to this new field. After the work of Ref. [39], Ricciardi et al. experimentally demonstrated direct frequency comb generation in an optical resonator with a single nonlinear crystal inside, originally conceived for cavity-enhanced SHG. OFCs were observed in the case of both phase-matched and phase-mismatched SHG. Moreover, the authors presented a simple theoretical model, which explained comb generation as the result of cascaded ${\chi }^{\left(2\right)}:{\chi }^{\left(2\right)}$ processes [48,49]. A more general theoretical description of comb generation in SHG cavity was successively developed by Leo et al., who modeled the dynamics of the cavity field in the time domain [50,51,52], and described comb formation in the framework of a modulation instability (MI), i.e., the growth of sidebands around a carrier frequency by amplification of small modulations on the carrier wave [53]. A similar theoretical description was adopted to describe the dynamics of quadratic combs observed in a degenerate OPO [54]. Finally, the most general approach, based on a single-envelope equation, has been also developed in order to study multi-octave, quadratic comb formation [55].

Quadratically nonlinear resonators thus emerge as the basis of an entirely new class of highly efficient synthesizers of OFCs, with unique features, such as the simultaneous generation of frequency combs in spectral regions far from the pump frequency and the role of phase matching in mitigating the effect of dispersion. Compared to Kerr combs, quadratic combs exploit the intrinsically higher efficiency of second-order nonlinearity, reducing the requirement in terms of pump power. Quadratic combs are still at an early stage but they are attracting the interest of an increasing number of research groups. More recent works are briefly reviewed in Section 6, where we conclude by giving an overview of promising developments of quadratic combs in terms of material platforms for chip-scale devices, steady low-noise dynamical regimes, and their potential interest for quantum optics.

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 LiNbO${}_3$ crystal, placed inside a traveling-wave optical resonator (free spectral range FSR = 493 MHz, quality factor $Q={10}^8$), resonating at the fundamental laser frequency ${\omega }_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].

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 ${\omega }_0\pm \Delta \omega$, symmetrically placed around the fundamental frequency (FF). As the power was further increased, additional sidebands appeared, displaced by multiples of $\Delta \omega$, 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.

Subsequently, wave vector mismatch $\Delta 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, $\Delta k>0$, the spectra (a)–(d) show widely separated sidebands, similar to the spectra observed at $\Delta 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 $\Delta 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.

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, $(\omega +\Delta \omega )+(\omega +\Delta \omega )\to 2\omega +2\Delta \omega$ and sum frequency with the fundamental wave $\omega +(\omega +\Delta \omega )\to 2\omega +\Delta \omega$, 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 ${\omega }_0$, and the parametric intracavity fields $A_{\mu }$ and $A_{\overline{\mu }}$, at ${\omega }_{\mu }={\omega }_0+\Delta \omega$ and ${\omega }_{\overline{\mu }}={\omega }_0-\Delta \omega$, respectively, which read [48]

$$\begin{array}{cc}\hfill {\dot{A}}_0=& -({\gamma }_0+i{\Delta }_0) A_0-2g {\eta }_{00\mu \overline{\mu }}{A}_0^{*}{A}_{\mu }{A}_{\overline{\mu }}-g({\eta }_{0000}|A_0{|}^2+2{\eta }_{0\mu 0\mu }|A_{\mu }{|}^2+2{\eta }_{0\overline{\mu }0\overline{\mu }}|A_{\overline{\mu }}{|}^2) A_0+F_{\mathrm{in}}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\dot{A}}_{\mu }=& -({\gamma }_{\mu }+i{\Delta }_{\mu }) A_{\mu }-g {\eta }_{\mu \overline{\mu }00} A_0^2{A}_{\overline{\mu }}^{*}-g(2{\eta }_{\mu 00\mu }|A_0{|}^2+{\eta }_{\mu \mu \mu \mu }|A_{\mu }{|}^2+2{\eta }_{\mu \overline{\mu }\mu \overline{\mu }}|A_{\overline{\mu }}{|}^2) A_{\mu }\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\dot{A}}_{\overline{\mu }}=& -({\gamma }_{\overline{\mu }}+i{\Delta }_{\overline{\mu }}) A_{\overline{\mu }}-g {\eta }_{\overline{\mu }\mu 00} A_0^2{A}_{\mu }^{*}-g(2{\eta }_{\overline{\mu }00\overline{\mu }}|A_0{|}^2+2{\eta }_{\overline{\mu }\mu \overline{\mu }\mu }|A_{\mu }{|}^2+{\eta }_{\overline{\mu }\overline{\mu }\overline{\mu }\overline{\mu }}|A_{\overline{\mu }}{|}^2) A_{\overline{\mu }} .\hfill \end{array}$$

(3)

Here, $F_{\mathrm{in}}=\sqrt{2{\gamma }_0/t_{\mathrm{R}}}{A}_{\mathrm{in}}$ is the cavity coupled amplitude of the constant input driving field $A_{\mathrm{in}}$, at frequency ${\omega }_0$; the $\gamma$’s are the cavity decay constants; the $\Delta$’s are the cavity detunings of the respective modes; $g={\left(\kappa L\right)}^2/2t_{\mathrm{R}}$ is a gain factor depending on the crystal length L (hereafter we consider the cavity length equal to the crystal length); $t_{\mathrm{R}}$ is the cavity round trip time; $\kappa =\sqrt{8}{\omega }_0{\chi }_{\mathrm{eff}}^{\left(2\right)}/\sqrt{c^3{n}_1^2{n}_2{ϵ}_0}$ is the second-order coupling strength. The latter is normalized so that the square modulus of the field amplitudes is measured in watts, with ${\chi }_{\mathrm{eff}}^{\left(2\right)}$ the effective second-order susceptibility, c the speed of light, $n_{1,2}$ the refractive indices, and ${ϵ}_0$ the vacuum permittivity. The integer mode number $\mu$ denotes the $\mu$th cavity mode, starting from the central mode at ${\omega }_0$, and overline stands for negative (lower frequencies). The $\eta$’s are complex nonlinear coupling constants, depending on the wave-vector mismatches associated with a pair of cascaded second-order processes,

$${\eta }_{\mu \sigma \rho \nu }=\frac{2}{L^2}{\int }_0^L{\int }_0^{z}exp[-i({\xi }_{\mu \sigma }z-{\xi }_{\rho \nu }{z}^{\prime })] \mathrm{d}{z}^{\prime } \mathrm{d}z$$

(4)

where ${\xi }_{jk}=k_{{\omega }_j}+k_{{\omega }_k}-k_{{\omega }_j+{\omega }_k}$.

A linear stability analysis of Equations (1)–(3) predicts the conditions for which a $\mu$-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 -g({\eta }_{\mu 00\mu }+{\eta }_{\overline{\mu }00\overline{\mu }}^{*}){\left|A_0\right|}^2\pm \sqrt{g^2|{\eta }_{\mu \overline{\mu }00}{|}^2 {\left|A_0\right|}^4-{\left[{\Delta }_0-D_2{\mu }^2-ig({\eta }_{\mu 00\mu }-{\eta }_{\overline{\mu }00\overline{\mu }}^{*}){\left|A_0\right|}^2\right]}^2},$$

(5)

where $D_2\simeq -2{\pi }^2{c}^3{\beta }^{″}/L^2 n_0^3=-(c/2n_0)D_1^2{\beta }^{″}$ accounts for the group velocity dispersion at ${\omega }_0$, with ${\beta }^{″}={\left.\frac{d^{2}k}{d{\omega }^2}\right|}_{{\omega }_0}$, and $n_0=n\left({\omega }_0\right)$ the refractive index at ${\omega }_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: ${\eta }_{\mu \overline{\mu }00}$, which is the parametric gain related to cascaded SHG and OPO, whereby two photons at frequency ${\omega }_0$ annihilate and two parametric photons at ${\omega }_{\mu }$ and ${\omega }_{\overline{\mu }}$ are created, mediated by a SH photon; and ${\eta }_{\mu 00\mu }$ (${\eta }_{\overline{\mu }00\overline{\mu }}$), which is related to the sum frequency process between a parametric photon at ${\omega }_{\mu }$ (${\omega }_{\overline{\mu }}$) 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_{\mu }$, nearly resonant with the $\mu$-th cavity mode,

$${\dot{A}}_{\mu }=-({\gamma }_{\mu }+i{\Delta }_{\mu }) A_{\mu }-g\sum _{\begin{array}{c}\rho ,\sigma \\ \nu =\rho +\sigma -\mu \end{array}}{\eta }_{\mu \nu \rho \sigma } A_{\nu }^{*}{A}_{\rho }{A}_{\sigma }+F_{\mathrm{in}} ,$$

(6)

where the summation over the indices $\rho$ and $\sigma$ goes over all the cavity resonant modes. The complex coupling constants are given by Equation (4), while the constraint over $\nu$ 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].

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,\tau )$ for the fundamental and $B(z,\tau )$ 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 mth round trip, along with the boundary condition for the fields at the end of each round trip. The propagation equations for the fields $A_m(z,\tau )$ and $B_m(z,\tau )$ read as

$$\begin{array}{cc}\hfill \frac{\partial A_m}{\partial z}=& \left[-\frac{{\alpha }_{c1}}{2}-i\frac{k_1^{″}}{2}\frac{{\partial }^2}{\partial {\tau }^2}\right]A_m+i\kappa B_m{A}_m^{*}{e}^{-i\Delta kz},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{\partial B_m}{\partial z}=& \left[-\frac{{\alpha }_{c2}}{2}-\Delta k^{\prime }\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},\hfill \end{array}$$

(8)

where $z\in [0,L]$ is the position along the cavity round trip path; ${\alpha }_{c1,2}$ are propagation losses (hereafter, subscripts 1 and 2 denote fields at ${\omega }_0$ and $2{\omega }_0$, respectively); $k_{1,2}^{″}={\mathrm{d}}^{2}k/\mathrm{d}{\omega }^2{|}_{{\omega }_0,2{\omega }_0}$ are the group velocity dispersion coefficients; $\Delta k^{\prime }=\mathrm{d}k/\mathrm{d}\omega {{|}_{2{\omega }_0}-\mathrm{d}k/\mathrm{d}\omega |}_{{\omega }_0}$ is the corresponding group-velocity mismatch or temporal walk-off. The “fast-time” variable $\tau$ describes the temporal profiles of the fields in a reference frame moving with the group velocity of light at ${\omega }_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 mth round trip according to the following cavity boundary conditions,

$$\begin{array}{cc}\hfill A_{m+1}(0,\tau )& =\sqrt{1-{\theta }_1} A_m(L,\tau ) e^{-i{\delta }_1}+\sqrt{{\theta }_1}{A}_{\mathrm{in}}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill B_{m+1}(0,\tau )& =\sqrt{1-{\theta }_2} B_m(L,\tau ) e^{-i{\delta }_2},\hfill \end{array}$$

(10)

where ${\theta }_{1,2}$ are power transmission coefficients at the coupling mirror, ${\delta }_1\simeq ({\omega }_0-{\omega }_{c1})t_{\mathrm{R}}$ and ${\delta }_2\simeq (2{\omega }_0-{\omega }_{c2})t_{\mathrm{R}}$ are the round trip phase detunings for the fields at ${\omega }_0$ and $2{\omega }_0$, respectively, with ${\omega }_{c1}$ and ${\omega }_{c2}$ the frequencies of the respective nearest cavity resonance, and $A_{\mathrm{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, ${\theta }_2=1$, and the SH field resets at the beginning of each round trip, i.e., $B_{m+1}(0,\tau )=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],

$$\begin{array}{cc}\hfill t_{\mathrm{R}}\frac{\partial A(t,\tau )}{\partial t}=& \left[-{\alpha }_1-i{\delta }_1-iL\frac{k_1^{″}}{2}\frac{{\partial }^2}{\partial {\tau }^2}\right]A-\rho A^{*}\left[A^2(t,\tau )\otimes I\left(\tau \right)\right]+\sqrt{{\theta }_1} A_{\mathrm{in}},\hfill \end{array}$$

(11)

where t is a “slow time” variable, linked to the round trip index as $A(t=m{t}_{\mathrm{R}},\tau )=A_m(z=0,\tau )$ [62,63,64,65], ${\alpha }_1=({\alpha }_{c1}L+{\theta }_1)/2$, $\rho ={\left(\kappa L\right)}^2$, ⊗ denotes convolution and the nonlinear response function $I\left(\tau \right)={\mathcal{F}}^{-1}\left[\widehat{I}\left(\mathsf{\Omega }\right)\right]$, with $\widehat{I}\left(\mathsf{\Omega }\right)=\left[(1-e^{-ix}-ix)/x^2\right]$, $x\left(\mathsf{\Omega }\right)=\left[\Delta k+i\widehat{k}\left(\mathsf{\Omega }\right)\right]L$, and $\widehat{k}\left(\mathsf{\Omega }\right)=-{\alpha }_{c,2}/2+i\left[\Delta k^{\prime }\mathsf{\Omega }+(k_2^{″}/2){\mathsf{\Omega }}^2\right]$. Here, we define the direct and inverse Fourier transform operator as $\mathcal{F}\left[·\right]={\int }_{-\infty }^{\infty }· e^{i\mathsf{\Omega }\tau } \mathrm{d}\tau$ and ${\mathcal{F}}^{-1}\left[·\right]={\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }· e^{-i\mathsf{\Omega }\tau } \mathrm{d}\mathsf{\Omega }$, 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({\alpha }_1+\rho P_0[\widehat{I}\left(\mathsf{\Omega }\right)+{\widehat{I}}^{*}(-\mathsf{\Omega })]\right)\pm \sqrt{|\widehat{I}{\left(0\right)|}^2{\rho }^2{P}_0^2-{\left({\delta }_1-\frac{k_1^{″}L}{2}{\mathsf{\Omega }}^2-i\rho P_0[\widehat{I}\left(\mathsf{\Omega }\right)-{\widehat{I}}^{*}(-\mathsf{\Omega })]\right)}^2},$$

(12)

which, baring the notation, is substantially equivalent to Equation (5). Figure 4a shows the MI gain, $\Re \left[{\lambda }_{+}\right]$ profile as a function of the walk-off parameter $\Delta k^{\prime }$. 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.

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_{\mathrm{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.

When ${\theta }_2<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,

$$\begin{array}{cc}\hfill t_{\mathrm{R}}\frac{\partial A}{\partial t}& =\left[-{\alpha }_1-i{\delta }_1-i\frac{k_1^{″}L}{2}\frac{{\partial }^2}{\partial {\tau }^2}\right]A+i\kappa LB{A}^{*}+\sqrt{{\theta }_1}{A}_{\mathrm{in}},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill t_{\mathrm{R}}\frac{\partial B}{\partial t}& =\left[-{\alpha }_2-i{\delta }_2-\Delta k^{\prime }L\frac{\partial }{\partial \tau }-i\frac{k_2^{″}L}{2}\frac{{\partial }^2}{\partial {\tau }^2}\right]B+i\kappa L{A}^2,\hfill \end{array}$$

(14)

where ${\alpha }_2$ 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_{\mathrm{R}}\frac{\partial A}{\partial t}=\left[-{\alpha }_1-i{\delta }_1-i\frac{k_1^{″}L}{2}\frac{{\partial }^2}{\partial {\tau }^2}\right]A-\rho A^{*}\left[A^2\otimes J\right]+\sqrt{{\theta }_1}{A}_{\mathrm{in}},$$

(15)

where the Fourier transform of the kernel function J is

$$\widehat{J}\left(\mathsf{\Omega }\right)=\frac{1}{{\alpha }_2+i{\delta }_2-i\Delta k^{\prime }L\mathsf{\Omega }-i\frac{k_2^{″}L}{2}{\mathsf{\Omega }}^2}.$$

(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_{\mathrm{in}}$ at the FF ${\omega }_0$ is replaced by a pump field $B_{\mathrm{in}}$ at the SH frequency $2{\omega }_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,

$$\begin{array}{cc}\hfill A_{m+1}(0,\tau )& =\sqrt{1-{\theta }_1} A_m(L,\tau ) e^{-i{\delta }_1}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill B_{m+1}(0,\tau )& =B_{\mathrm{in}}.\hfill \end{array}$$

(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 $\Delta k=0$ [54],

$$\begin{array}{cc}\hfill t_{\mathrm{R}}& \frac{\partial A(t,\tau )}{\partial t}=\left[-{\alpha }_1-i{\delta }_1-i\frac{L{k}_1^{″}}{2}\frac{{\partial }^2}{\partial {\tau }^2}\right] A(t,\tau )-{\mu }^2{A}^{*}(t,\tau )\left[A^2(t,\tau )\otimes I\left(\tau \right)\right]+i\mu B_{\mathrm{in}}{A}^{*}(t,\tau ),\hfill \end{array}$$

(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, $A_0=0$, and a nontrivial time independent solution, $A_0=\left|A_0\right|e^{i\varphi }$. From a linear stability analysis of the constant solution, we derived the following expression for the eigenvalues [54],

$${\lambda }_{\pm }=-\left[{\alpha }_1+{\mu }^2{\left|A_0\right|}^2{\mathcal{I}}_{+}\left(\mathsf{\Omega }\right)\right]\pm \sqrt{({\alpha }_1^2+{\delta }_1^2)-{\left[{\delta }_1-D_2{\mathsf{\Omega }}^2-i{\mu }^2{\left|A_0\right|}^2 {\mathcal{I}}_{-}\left(\mathsf{\Omega }\right)\right]}^2} ,$$

(20)

where $|A_0{|}^2=[-{\alpha }_1\pm \sqrt{{\mu }^2{B}_{\mathrm{in}}^2-{\delta }_1^2}]/{\mu }^2 \widehat{I}\left(0\right)$ is the squared modulus of the nontrivial solution and ${\mathcal{I}}_{\pm }\left(\mathsf{\Omega }\right)=\widehat{I}\left(\mathsf{\Omega }\right)\pm {\widehat{I}}^{*}(-\mathsf{\Omega })$. Similarly, for the zero solution the eigenvalues of the linearized system are

$$\begin{array}{cc}\hfill {\lambda }_{\pm }=& -{\alpha }_1\pm \sqrt{{\mu }^2{B}_{\mathrm{in}}^2-{\left({\delta }_1-D_2{\mathsf{\Omega }}^2\right)}^2}.\hfill \end{array}$$

(21)

Both solutions exhibit MI gain for $\mathrm{Re}\left[{\lambda }_{+}\right]>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 $\Delta k^{\prime }$, through ${\mathcal{I}}_{\pm }\left(\mathsf{\Omega }\right)$, 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 ${\omega }_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 ${\omega }_0$. The larger the detuning, the more distant the symmetric resonances are. For small negative detunings no resonance pair can be symmetrically displaced around ${\omega }_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.

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 $\mathsf{\Lambda }=6.92 \mathsf{\mu }\mathrm{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 $\Delta =-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.

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}\left[A^{m+1}(\tau ,0)\right]=\sqrt{\widehat{\theta }\left(\mathsf{\Omega }\right)}\mathcal{F}\left[A_{\mathrm{in}}\right]+\sqrt{1-\widehat{\theta }\left(\mathsf{\Omega }\right)} e^{i{\varphi }_0}\mathcal{F}\left[A^m(\tau ,L)\right]$$

(22)

$$\left[{\partial }_z-D\left(i\frac{\partial }{\partial \tau }\right)+\frac{{\alpha }_d}{2}\right]A^m(\tau ,z)=i{\rho }_0\left(1+i{\tau }_{\mathrm{sh}}\frac{\partial }{\partial \tau }\right)p_{\mathrm{NL}}(\tau ,z,A^m).$$

(23)

The boundary condition, Equation (22), is written in the Fourier domain, in order to account for the frequency dependence of the transmission coefficient $\theta$ at the input port of the resonator. It determines the intra-cavity field $A_{m+1}(\tau ,z=0)$ at the beginning of $(m+1)$th round trip in terms of the field at the end of the previous round trip $A_m(\tau ,z=L)$ and the pump field $A_{\mathrm{in}}$. Equation (23) is written in a reference frame moving at the group velocity at ${\omega }_0$: $p_{\mathrm{NL}}$ is the broadband envelope of the nonlinear polarization $P_{\mathrm{NL}}=P_{\mathrm{NL}}^{\left(2\right)}+P_{\mathrm{NL}}^{\left(3\right)}+\dots =ϵ({\chi }^{\left(2\right)}{E}^2+{\chi }^{\left(3\right)}{E}^3+\dots )$; ${\rho }_0={\omega }_0/2n_{0}c{ϵ}_0$; ${\tau }_{\mathrm{sh}}=1/{\omega }_0$ is the shock coefficient that describes the frequency dependence of the nonlinearity, and ${\alpha }_d$ 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\ge 2}i\frac{{\beta }_l}{l!}{\left(i\frac{\partial }{\partial \tau }\right)}^l,$$

(24)

where ${\beta }_l={(d^l\beta /d{\omega }^l)}_{\omega ={\omega }_0}$ are expansion coefficients of the propagation constant $\beta \left(\omega \right)$.

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 $\mathsf{\mu }$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).

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,79,80,81,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 ${\chi }^{\left(2\right)}$ processes [89,90,91,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,95,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” ${\chi }^{\left(2\right)}$ 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,100,101,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,109,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].