Frequency combs in a microring optical parametric oscillator
A. Villois,1 N. Kondratiev,2	I. Breunig,3,4 D. N. Puzyrev,1 and D. V. Skryabin1,2,*
1Department of Physics, University of Bath, Bath BA2 7AY, UK
2Russian Quantum Centre, Skolkovo 143025, Russia
3Department of Microsystems Engineering—IMTEK, University of Freiburg, 79110 Freiburg, Germany
4Fraunhofer Institute for Physical Measurement Techniques, 79110 Freiburg, Germany
*Corresponding author: d.v.skryabin@bath.ac.uk
Received 29 July 2019; accepted 11 August 2019; posted 13 August 2019 (Doc. ID 373956); published 3 September 2019

We report the soliton frequency comb generation in microring optical parametric oscillators operating in the downconversion regime and with the simultaneous presence of the χ ^2^ and Kerr nonlinearities. The combs are studied considering a typical geometry of a bulk LiNbO3 toroidal resonator with the normal group velocity dispersion spanning an interval between the pump and the down- converted signal. We have identified critical power signaling a transition between the relatively low pump power predominantly χ ^2^ combs and the high pump power ones shaped by the competition between the χ ^2^ and Kerr nonlinearities. © 2019 Optical Society of America
https://doi.org/10.1364/OL.44.004443
Frequency combs are a key concept and a tool in the fields of short pulse generation, optical metrology, and spectroscopy [1]. The relatively recent development of the high-quality factor whispering gallery mode microring resonators opened up doors for the chip and micro-scale frequency comb (microcomb) devices [2]. The association of the microcombs with the soliton solutions of the Lugiato–Lefever (LL) equation [3] has led to the emergence of a burgeoning research area of the microring solitons supported by the interplay of the positive Kerr nonlinearity and anomalous group velocity dispersion (GVD) [4]. Reducing power thresholds for the soliton comb and extending their range of existence to the normal dispersion are often desirable for practical applications and, therefore, quadratic nonlinear response, χ^2^ , of the noncentrosymmetric materials is a viable exploration direction in the microcomb area that has been recently studied by a number of authors [5–12]. Importantly, the comb solitons embracing the second harmonic are naturally octave wide and, therefore, suitable for selfreferencing [1,2]. Among the possible candidate materials to study χ^2^ soliton combs, the lithium-niobate (LN) results are particularly promising, given the large χ^2^ coefficient and its transparency over a wide spectral window. Note that the area of quadratic comb solitons connects to the prior research on the quadratic solitons in bulk [13] and parametric conversion in microresonators [14].
While the soliton combs associated with the second- harmonic generation process in microring resonators have received some preliminary attention in the recent publications [5,6,11,12], we demonstrate below that the parametric downconversion is equally promising for the soliton comb generation. We account for both χ ^2^ and χ ^3^ nonlinearities and demonstrate how their competition may alter the comb generation.
Recent experiments with a microring pumped with the external comb signal with the matched repetition rate have reported the line by line conversion into multiple higher and lower frequency bands [8]. In addition, recent publications on the meter scale bow-tie cavities loaded with a χ^2^ crystal and operating in the downconversion regime have reported low-coherence non-soliton combs with very few sidebands due to GVD-induced instabilities [16] and the so-called simul- ton-solitons supported solely by the group velocity mismatch, where GVD plays no role [17]. While the coherent multisideband microring parametric soliton combs reported below rely on the normal GVD and, thus, are the first of its kind.
We start from the averaged along the transverse dimensions scalar wave equation ∂z2E - c-2∂t2 -∞∞ g^t - t0^ E^t0^dt0 ^ μ0∂t2hP, f i∕hf , f i. Here z ^ Rθ is the coordinate along the ring, R is the ring radius, θ is the azimuthal angle, fE is the electric field, and c-2 ^ μ0ε0. f is the dimensionless transverse modal profile with the modal area S ≃ 30 μm2 . h·, ·i means integration over the transverse dimensions of the modal profile. g is the resonator response function accounting for a combination of the material and geometrical dispersions, such that -∞∞ g ^t^e-iωt dt ^ ne2 ^ω^, where ne ^ω^ is the effective refractive index. P ^ ε0χ2f 2E2 ^ ε0χ3f 3E3 is the nonlinear polarization combining quadratic, χ2 ≃ 5 · 10-12 m∕V, and cubic, χ3 ≃ 10-21 m2 ∕V2 , nonlinear susceptibilities of LN. f and χ2,3 are assumed non-dispersive for the sake of brevity and transparency of notations. For P ^ 0, E ∼ eimθ-iωmt, the linear resonator spectrum is recovered from c2 m2 ^ R2ω2mne2 ^ωm^. Taking E ^ bΣmAmeimθ-iωmt , P ^ Pm Pmeimθ-iωmt, assuming that the Fourier amplitudes vary slowly in time, and dispersion and nonlinearity are relatively small, we find ib∂t Am ^ - 2εω0mn2 hPm, f i∕hf ,f i, where n ≃ 2.2. Here

0146-9592/19/184443-04 Journal © 2019 Optical Society of America

b2 ^ 2∕^ε0cnS^ ≃ 1013 V2 m-2 W-1 and, thereby, jAmj2 is measured in watts. The above procedure renders explicit expressions for the nonlinear coefficients in the LL model; see γ2,3 below.
The derivation should be applied for both extraordinary (polarized along the resonator axis) and ordinary (polarized in the resonator plane) waves, which appropriately couple through the nonlinear polarization tensor. We now assume that the extraordinary pump (p) and ordinary signal (s) waves have their own modal profiles fp,s , and their electric fields are given by bp,s f p,s ψ p,s ^t, θ^eimp,sθ-iωp,st ^ c:c:, where ψp,s are the intracavity pump and signal amplitudes spectrally capturing a sufficiently large number of the resonator lines, and bp2,s ^ 2∕^ε0 cnSp,s^ are the respective power scaling coefficients. ωp is the cavity resonance nearest to the pump frequency Ω, and ωs is the one nearest to Ω∕2. Then transforming from the modal expansion into the real space as per, e.g., Refs. [15,22], including the pump field with frequency Ω nearest to ωp and losses, we obtain the LL-like system
i∂tψs ^
δs
iD1s∂θ
- 12 D2s ∂2θ ψs
- iκsψs
-		γ2sψs^ψpei2εt^i^mp-2ms^θ - γ3s^jψsj2 ^ 2asjψpj2^ψs,
(1)
1
i∂tψp ^ δp - iD1p∂θ - D2p ∂2θ ψp - iκpψp ^ h
2
-		γ2pψs2e-i2εt-i^mp-2ms^θ - γ3p^jψpj2 ^ 2apjψsj2^ψp:
(2)
Here δs ^ ωp∕2 - Ω∕2 and δp ^ ωp - Ω are the detuning parameters. In the following, we assume that χ2 varies periodically in θ, so that a mismatch of the modal indices, mp - 2ms, is compensated for in the leading order, while the higher- order corrections are disregarded. ε ^ ωs - ωp ∕2 is the frequency mismatch parameter. D1p,1s are the local free spectral ranges (FSRs) and D2p,2s set the GVD associated FSR dispersion. γ 0s are the quadratic and cubic nonlinear coefficients. In the following, we assume for brevity γ2s,2p ≃ ωpχ2N2pbp	ωpχ3N3pbp2	R j^1
γ2 ^	2n2 , γ3s,3p ≃ γ3 ^	2n2 , where Njp ^ f p
dS^ f p2dS j^-1 . For the parameters specified aboveγ2 ≃ 9 · 109 W-1∕2 s-1 and γ3 ≃ 106 W-1 s-1 . h ^ pffiηffiffiκffiffipffiffiDffiffiffi1ffiffipffiffiPffiffiffi∕ffiffi^ffiffi2ffiffiffiπffiffi^ffi is the pump parameter, where P is the pump power [15,22]. The quality factor of LN microresonators strongly depends on the chosen platform. The chip integrated samples are known to reach Q ≃ 107 [18]. However, the bulk LN resonators can have a better Q, e.g., the ones used in Ref. [19] have had extinction coefficients r ≃ 5 · 10-4 cm-1 , which corresponds to the intrinsic quality factor 2πn∕^rλ^ ≃ 2 · 108 in the required wavelength range, λ. Total Q is half of the above for the critical coupling (η ^ 1∕2) and, thus, for the moderate Q ^ 107 ,we have 2κp,s ≃ ωp,s ∕Q ∼ 108 s-1 .
To determine practical FSRs and GVDs, we consider a spheroidal LN microresonator (major radius R ∼ 500 μm, minor radius r ∼ 250 μm) and use analytical expressions for the resonance frequencies [20]. By varying r and R, we found a pair of resonances ωp (extraordinarily in the proximity of 0.8 μm) and ωs (ordinary in the proximity of 1.6 μm) such that the
FSR parameters D1p and D1s are approximately matched and, simultaneously, the frequency mismatch parameter ε is reduced to nearly zero. Close matching of FSRs is also known as an important condition to achieve soliton combs via second-harmonic generation [6]. Choosing ωp∕^2π^^ 381.749381 THz (azimuthal mode number mp ^ 8668) gives ωs∕^2π^ ^ 190.874707 THz (azimuthal mode number ms ^ 4401). Then ε ≃ 0, D1s ≃ D1p ^ 2π · 42.17 GHz, and D2s ^ -2π · 166.94 kHz, D2p ^ -2π · 548.30 kHz. Non-zero values of ε and 1 - D1s ∕D1p do not qualitatively alter the results described below.
The nonlinear resonances of Eqs. (1) and (2) are characterized by their continuous-wave (cw) solutions (∂t ^ ∂θ ^ 0). In the downconversion setting, there exists a cw state with the nontrivial pump, ψp ≠ 0, and the zero signal, ψs ^ 0. This pump-only state exists both above and below the threshold, where the parametric state, ψp,s ≠ 0, splits away from ψs ^ 0. The threshold condition, h2 γ22 ^ ^δp2 ^ κp2^^δs2 ^ κs2^, can be found explicitly for γ3 ^ 0. The blue and red lines in Fig. 1 mark the pump-only and the parametric cw states, respectively. The resonance line of the pump-only state tilts towards positive values of δp through the Kerr effect-induced bistability. By approximately solving the Kerr bistability equation, h2 ^ ^κp2 ^^δp - γ3jψpj2^2^jψpj2, in the limit of large detunings, we find that the extent of the tilt in δp is approximated as δpKerr ^ γ3h2∕κp2^1 ^ O^hγ13∕2∕jδpj3∕2^^.
In the presence of the parametric effect, the linear cavity resonance splits into two resonances oppositely tilted towards both negative and positive detunings; see Fig. 1, which is similar to what was previously described for the second- harmonic generation in microrings [5,6]. Thus, the parametric effect acts similarly to the positive/negative Kerr nonlinearity for δp > 0∕δp < 0. For γ3 ^ 0, the extent of the parametric tilt can be found from the exact condition h2γ22 ^ ^κs δp ^ κpδs^2 , so that for ε ^ 0 and κp ^ κs, we get δpprm ^ ^2γ2h∕^3κp^. While the width of the parametric resonance in ψs field is proportional to κs , the resonance line in the
Fig. 1. (a), (b) Pump and (c), (d) signal amplitudes of the cw states versus detuning. The parametric effect is dominant in (a) and (c) (Q ^ 107, P ^ 10 mW ≪ Pcr, χ2∕χ3 ^ 5· 109 V∕m). The Kerr and parametric effects compete in (b) and (d) (Q ^ 108 , P ^ 2.5W∼ Pcr, χ2 ∕χ3 ^ 5 · 107 V∕m). The blue/red lines show the pump-only/parametric cw states. The black lines denote the maxima of the intensities of the bright parametric solitons. The zooms for the rectangular areas in (c) and (d) are shown in the insets. The thinner lines mark unstable cw states and solitons.

(3)

Fig. 2. (a) MI growth rate Re Λ versus δp ^ 2δs and the modal index μ for the high amplitude parametric cw state. (b) Spatial distribution of the signal amplitude after a half-million round trips versus δp . θ is shown between -π∕2 and -π for better visibility. (c), (d) Log-scale modal spectrum of the pump/signal after a half-million round trips versus δp. All parameters are as in Figs. 1(a) and 1(c), i.e., P ≪ Pcr.

Fig. 3. Time-space evolution of the signal amplitudes when the cavity is initialized with an unstable cw state for (a) P ≪ Pcr and (b) P ∼ Pcr . The parameters in (a) and (b) are as in Figs. 1(a) and 1(b) with δp∕^2π^ ^ -1.97∕-0.79 GHz. θ is shown only between 0 and π for better visibility.

ψ component has a much smaller, even vanishingly small p
(for γ3 ^ 0) width; see Figs. 1(a) and 1(c).
Thus, the nonlinear shift of the resonance frequency is proportional to h2 for the Kerr effect and to h for the parametric one. Hence, the latter dominates for the relatively low pump powers, when jδpprm j ≫ δpKerr , while the two processes start to compete when jδpprm j ∼ δpKerr, i.e., at h ∼ hcr ^ 2γ2κp∕^3γ3^, which corresponds to P ∼ Pcr , where
4πωp γ22	4π	χ22 N 22p
Pcr	9ηQD1pγ23	9ηF · b2χ23 ·N23p, where F ^ QD1p ∕ωp is the cavity finesse. Using parameters as in the above text and Q ^ 107 we find Pcr ≃ 70 kW. By taking a higher, but still realistic, Q ^ 108 and a material with χ2 ∕χ3 ^ 5 · 107 V∕m, instead of previous 5 · 109 V∕m, we find a very practical Pcr ≃ 1W. The oppositely tilted parametric resonances are expected to lose their symmetry for P approaching Pcr; see Figs. 1(b) and 1(d). For the positive detunings, the parametric and Kerr effects impact the resonance in the same way by extending its range towards larger positive detunings. However, for the negative detunings nonlinearities, competing nonlinearities reshape the negatively tilted resonance by reversing the tilt direction. The signal part of the folded back solution hits ψs ^ 0 at some positive δp , while the corresponding ψp naturally merges with the pump-only cw state.
To study stability of the cw solutions with respect to the sideband excitation (modulational instability [MI]), we make an ansatz ψp,s ^θ, t^ ^ ψp,s ^ ^up,s ^t, θ^ ^ ivp,s ^t, θ^^, where up,s , vp,s are small perturbations. Solving the linearized problem for the latter numerically, we find the sideband growth rates Re Λ as a function of their modal indexes μ: up,s , vp,s ∼ eiμθ^Λt . The results of this analysis are shown in Fig. 1, using thick/thin lines for the stable/unstable cw states. Figure 2(a) shows the MI growth rate of the parametric P ≪ Pcr cw state with the highest amplitude as a function of δp ^ 2δs and of the modal index μ. The MI structure remains qualitatively similar for P ∼ Pcr.
To check the outcomes of the MI development, we have performed an extensive series of dynamical simulations. Figure 2(b) shows the spatial distributions of the signal field intensity, jψs j2 , after a half-million cavity round trips for a broad interval of the detunings and for P ≪ Pcr. Small amplitude and spectrally narrow waves are generated in the system for the large positive δp ; see Figs. 2(b)–2(d). More modes and broad, albeit non-solitonic, frequency combs are excited and the system exhibits complex spatio-temporal dynamics characterized by the high amplitude oscillations for δp approaching zero. When δp ∼ 0, the cw solutions becomes stable in a narrow interval of δp . For the negative values of δp the system enters the soliton regime; see the well localized spots in Fig. 2(b). Note that the initial noise realisations applied vary as δp is scanned. Figure 3 shows the spatio-temporal dynamics of the signal field component leading to the formation of multiple solitons from the cw input. For P ∼ Pcr; see Fig. 3(b), the solitary filaments exhibit a strong breathing dynamics associated with the repulsion and attraction of the steep wave fronts connecting the low and high amplitude cw states. With more round trips, the breathers are gradually converging to the multisoliton state.
To find stationary soliton solutions, we have numerically integrated Eqs. (1) and (2) with ∂t ^ 0, using a Newton– Raphson method in combination with a bi-conjugate gradient method. Fourier space with 2048 modes was used so that the periodic boundary conditions were accounted for automatically. The maxima of the soliton amplitudes are plotted in Fig. 1 using solid black lines. Soliton families were found to have the resonance shapes, similar to the cw ones. The linear stability of the solitons was also analyzed numerically. The higher amplitude soliton branch is typically stable, and the low amplitude one is typically unstable. Both nonlinearities provide a positive shift of the cavity resonance for δp ^ 2δs > 0, which does not favor bright solitons for the normal dispersions, D2p,2s < 0. If, however, δp ^ 2δs < 0, then the negative shift of the resonance by the quadratic nonlinearity prevails for P ≪ Pcr and creates conditions for stable bright solitons. It is instructive to present an approximate analytical soliton solution, which can be found assuming γ3 ^ D2p ^ 0, D1p ^ D1s, κp∕δp ≪ 1. Then, for ψpδp ≃γ2ψs2 - h,

Fig. 4. Parametric (a), (b) soliton profiles and (c), (d) and the corresponding comb spectra for P ≪ Pcr. (a), (c)/(b), (d) correspond to the points 1/2 in Fig. 1(c). The blue and red colors correspond to the pump and signal, respectively.

Fig. 6. Comparison between the pulse amplitudes of the pump field obtained in the dynamical modeling (blue circles) with the stationary soliton amplitudes (black lines) for the (a) P ≪ Pcr and (b) P ∼ Pcr cases. The parameters are as in Figs. 1 and 2.

x ^ θ - D1pt, we find that ψs ^ ψ eiα, where ψ is a solution of 12 D2s∂x2ψ ^ q^ψ - γ22ψ3∕δp and sin 2α ^ -κsδp ∕^γ2 h^ [21]. Thereby,
ψs ^ eiα qffi2ffiffiδffiffipffiffiqffiffi^ffiffiffi · sech x 2q^, γ2	D2s
γ22h2
q^ ^ δs ^	2 - κs2:
δp

(4)

Equation (4) provides a good qualitative understanding of the parametric soliton combs. It explicitly predicts that the bright comb solitons exist for normal GVD, D2s < 0, only for the negative detunings, δp ^ 2δs < 0, and that the high amplitude q- branch merges with the always unstable q^ branch at some critical detuning. Examples of the numerically found soliton profiles and their associated combs for P ≪ Pcr are shown in Fig. 4. The low pump solitons have typical sech profiles. They are becoming narrower, and the associated combs are becoming broader with the increased detuning, as expected from Eq. (4). Equation (4) is not applicable for P ∼ Pcr, since it completely disregards Kerr nonlinearity. In this case, the soliton branches show a tendency to fold back, together with the parametric cw’s; see Figs. 1(b) and 1(d). The two soliton branches do not join together as in the P ≪ Pcr case, but rather terminate at a cw state; see the inset in Fig. 1(d). For the relatively low intensities, these solitons resemble the ones for P ≪ Pcr ; see Fig. 5(b). However, when their intensities are high, the pulse profiles become the square-like, connecting the low and high amplitude cw states via steep fronts; see Fig. 5(a).

Fig. 5. Same as Fig. 4, but for P ∼ Pcr. (a), (c)/(b), (d) correspond to the points 1/2 in Fig. 1(d).

The associated comb develops the non-sech tails typical for the square pulses. In order to confirm that pulses emerging from the cw instabilities are indeed the soliton pulses, we have compared their amplitudes with stationary soliton amplitudes for both P ≪ Pcr and P ∼ Pcr cases and found an impressive agreement; see Fig. 6.
Funding. Horizon 2020 Framework Programme (812818); Russian Foundation for Basic Research (17-0200081); Leverhulme Trust (RPG-2015-456).
REFERENCES
1.		S. T. Cundiff and J. Ye, Rev. Mod. Phys.75, 325 (2003).
2.		A. Pasquazi, M. Peccianti, L. Razzari, D. J. Moss, S. Coen, M. Erkintalo, Y. K. Chembo, T. Hansson, S. Wabnitz, P. Del’Haye, X. X. Xue, A. M. Weiner, and R. Morandotti, Phys. Rep.729, 1 (2018).
3.		T. Herr, V. Brasch, J. D. Jost, C. Y. Wang, N. M. Kondratiev, M. L. Gorodetsky, and T. J. Kippenberg, Nat. Photonics8, 145 (2014).
4.		T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, Science361, eaan8083 (2018).
5.		T. Hansson, P. Parra-Rivas, M. Bernard, F. Leo, L. Gelens, and S. Wabnitz, Opt. Lett.43, 6033 (2018).
6.		A. Villois and D. V. Skryabin, Opt. Express27, 7098 (2019).
7.		P. Parra-Rivas, L. Gelens, T. Hansson, S. Wabnitz, and F. Leo, Opt. Lett.44, 2004 (2019).
8.		S. J. Herr, V. Brasch, J. Szabados, E. Obrzud, Y. Jia, S. Lecomte, K. Buse, I. Breunig, and T. Herr, Opt. Lett.43, 5745 (2018).
9.		D. J. Wilson, K. Schneider, S. Hoenl, M. Anderson, T. J. Kippenberg, and P. Seidler, “Integrated gallium phosphide nonlinear photonics,” arXiv:1808.03554 (2018).
10.		Y. He, Q.-F. Yang, J. Ling, R. Luo, H. Liang, M. Li, B. Shen, H. Wang, K. Vahala, and Q. Lin, “Self-starting bi-chromatic LiNbO3 soliton Microcomb,” arxiv:1812.09610 (2018).
11.		X. Xue, X. Zheng, and B. Zhou, Photonics Res.6, 948 (2018).
12.		R. Ikuta, M. Asano, R. Tani, T. Yamamoto, and N. Imoto, Opt. Express26, 15551 (2018).
13.		A. V. Buryak, P. Di Trapani, D. V. Skryabin, and S. Trillo, Phys. Rep.370, 63 (2002).
14.		I. Breunig, Laser Photonics Rev.10, 569 (2016).
15.		Y. K. Chembo, I. S. Grudinin, and N. Yu, Phys. Rev. A92, 043818 (2015).
16.		S. Mosca, M. Parisi, I. Ricciardi, F. Leo, T. Hansson, M. Erkintalo, P. Maddaloni, P. De Natale, S. Wabnitz, and M. De Rosa, Phys. Rev. Lett.121, 093903 (2018).
17.		M. Jankowski, A. Marandi, C. R. Phillips, R. Hamerly, K. A. Ingold, R. L. Byer, and M. M. Fejer, Phys. Rev. Lett.120, 053904 (2018).
18.		M. Zhang, C. Wang, R. Cheng, A. Shams-Ansari, and M. Loncar, Optica4, 1536 (2017).
19.		M. Leidinger, S. Fieberg, N. Waasem, F. Kuhnemann, K. Buse, and I. Breunig, Opt. Express23, 21690 (2015).
20.		M. L. Gorodetsky and A. E. Fomin, IEEE J. Sel. Top. Quantum Electron.12, 33 (2006).
21.		I. V. Barashenkov, M. M. Bogdan, and V. I. Korobov, Europhys. Lett.15, 113 (1991).
22.		Y. K. Chembo and N. Yu, Phys. Rev. A82, 033801 (2010).