Ladder of Eckhaus instabilities and parametric conversion in chi(2) microresonators

Low loss microresonators have revolutionised nonlinear and quantum optics over the past decade. In particular, microresonators with the second order, chi(2), nonlinearity have the advantages of broad spectral tunability and low power frequency conversion. Recent observations have highlighted that the parametric frequency conversion in chi(2) microresonators is accompanied by stepwise changes in the signal and idler frequencies. Therefore, a better understanding of the mechanisms and development of the theory underpinning this behaviour is timely. Here, we report that the stepwise frequency conversion originates from the discrete sequence of the so-called Eckhaus instabilities. After discovering these instabilities in fluid dynamics in the 1960s, they have become a broadly spread interdisciplinary concept. Now, we demonstrate that the Eckhaus mechanism also underpins the ladder-like structure of the frequency tuning curves in chi(2) microresonators. A series of discrete transitions between the different mode pairs in the output of the microresonator OPOs has hinted at a connection with Eckhaus instabilities originally discovered in fluid dynamics. The theory of this effect is developed to describe the OPO signal tuning by the pump laser frequency adjustment in the proximity of one resonance.

N onlinear and quantum optics in the ring microresonators have attracted a great deal of attention over the past few years 1,2 . Its applications include classical and quantum information processing [3][4][5] , precession spectroscopy 1 and exoplanet discovery 6 . While the third-order nonlinear response (Kerr effect) of microresonators remains the workhorse of research in this area, the second-order, χ (2) , nonlinearity is attracting an increased attention 2,7,8 . χ (2) response is generally stronger than Kerr and offers greater flexibility with the spectral tunability and access to a range of quantum effects 5,[9][10][11] . The respective material and fabrication issues have matured significantly 12 , and the phase-matching by the self-induced gratings in Si 3 N 4 has been demonstrated [13][14][15] .
With all the innovations in mind, the temperature, T, control of the refractive index to provide the signal-idler-pump energy conservation, ℏω s + ℏω i = ℏω p , remains an essential method of the frequency tuning of the down-converted signal in microresonators since their very early days 16,17 and till now 18,19 . For the total (material+geometry) resonator dispersion being well approximated by the second-order expansion, the index matching curves, T vs ω s,i , take the parabolic shape. It was demonstrated already in Ref. 17 , and then has become a staple for many results on frequency conversion in χ (2) microresonators, see, e.g., [19][20][21] and references therein. Recent parametric down-conversion results in the thin-film quasi-phase-matched 20 and whispering gallery 22 LiNbO 3 microresonators have also used the pump laser frequency as a control parameter leading to the parabolic ω p vs ω s,i tuning curves 20,22 . In what follows, we develop a theory underpinning the pump-frequency-tuning approach.
A distinct feature of the tuning process in microresonators is that the signal and idler frequencies vary discretely, i.e., with a jump, when one signal-idler pair switches to the next in sequence 19,20,22 . This discreteness raises some interesting questions. For example, is there a bifurcation scenario leading to the stepwise mode-number switch, or this is a bifurcation free parameter-drag effect, and what could be the locking interval in the parameter space providing the device operation in the desired pair of the resonator modes? Our theory answers the above questions and reveals that the instability scenario, or rather a ladder-like sequence of the instabilities, behind the tuning of the parametric frequency conversion in microresonators is known in the contexts of fluid dynamics and pattern formation as the Eckhaus instability [23][24][25][26] . The term of the Eckhaus, or so-called secondary, instability is used to describe the bifurcation that changes the period of a nonlinear wave [23][24][25][26][27][28][29][30] . We also find an explicit condition on the power vs loss balance which explains the cutting of the optical parametric oscillator (OPO) tuning parabola at one or both of its ends and limits the achievable spectral separation between the signal and idler photons.

Results and discussion
Model and classification of OPO regimes. Frequencies of the modes in high-finesse resonators are rigidly linked to the quantised set of the wavenumbers. Therefore, the step in frequency also implies the abrupt change of the spatial period, as it happens in the Eckhaus instability scenario. Accounting for the quantisation of the spectrum calls for the model formulation and methodology different from the continual models of transverse nonlinear optics 27,[31][32][33][34][35] , fibre-loop 28,29,36 and bow-tie 37-39 cavities, cf. with the soliton and frequency conversion theories in the high-finesse microresonators [40][41][42][43][44][45][46][47][48] .
We now set the model in a transparent and rigorous manner. We assume that the pump laser frequency, ω p , is tuned around the frequency, ω 0b , of the resonator mode with the number 2M. 2M equals the number of wavelengths fitting along the ring circumference. We also express the multi-mode intra-resonator electric field centred around ω p and its half-harmonic, ω p /2, as e iMϑÀi 1 2 ω p t ∑ μ a μ e iμθ þ c:c:; respectively. Here, ϑ ¼ 0; 2π ð is the angular coordinate, and θ is its transformation to the reference frame rotating with the rate D 1 /2π. 'a' and 'b' mark the half-harmonic, i.e., signal, and pump fields, respectively. μ = 0, ±1, ±2, … is the relative mode number, and the resonator frequencies are where, D 1ζ /2π are the free spectral ranges, FSRs, and D 2ζ are dispersions. In what follows, we choose D 1 = D 1a . The frequency, i.e., phase, matching parameter for the non-degenerate parametric process is defined as Here, n m is the effective refractive index taken for the frequencies of the modes with the absolute numbers m = M ± μ (signal and idler) and 2M (pump), c is the vacuum speed of light and R is the resonator radius. For example, ε 0 = 2ω 0a − ω 0b = 0 corresponds to the exact matching for the degenerate parametric conversion, n M = n 2M .
Simple algebra reveals that all modal detunings in the halfharmonic signal can be expressed via δ 0b and the respective Fig. 1 Generic and staggered frequency combs. a, b is an illustration of the generic comb with the spatial period 2π/ν. c, d is the staggered comb with the period π/ν. Red and green colours mark the signal and pump combs, respectively. The resonator parameters 7 are: linewidths: κ a /2π = 1 MHz, κ b / 2π = 2 MHz; repetition rates: D 1a /2π = 21 GHz, D 1b /2π = 20 GHz; dispersions: phase-matching parameters, While δ μa and ε μ do not depend on the repetition-rate difference, δ μb does, Depending on the pump power, the classical half-harmonic signal can be either zero or not. This is reflected in the structure of Eq. (4) and its solutions. Three types of solutions we should highlight are (i) no-OPO state: (ii) degenerate OPO state: and a family of (iii) non-degenerate OPO states: Though Eq. (4) does not have a closed analytical solution for the non-degenerate OPO, the experimental data demonstrate that the states with the |b 0 | 2 and |a ±ν | 2 powers strongly dominating across the whole spectrum both exist and can be tuned to change from one ν to the other, and therefore, they represent the practically desirable regimes of the microresonator operation [18][19][20][21] . The explicit expressions for the non-zero modes in Eqs. (8) and (9) are introduced later.
In the fluid dynamics context [23][24][25] , the transition from the no-OPO to an OPO state would correspond to the Benjamin-Feir instability (emergence of the signal), while moving from the OPO state operating in the ±ν pair of modes to the ±(ν + 1) pair would be the Eckhaus instability, see, e.g., Ref. 23 that elucidates the difference between the two. The general form of Eq. (4) is not tractable analytically regarding the issue of the interplay between the different OPO states, and, therefore, in the next section, we derive the reduced model that allows such analysis to be carried out in a transparent form.
Apart from the OPO regimes listed above, Eq. (4) allows for the multi-mode frequency comb solutions that could be either stationary or time-dependent. The left column of Fig. 1 schematically illustrates a solution of Eq. (4) corresponding to the generic frequency comb with the spatial period 2π/ν. The structure of Eq. (4) also admits a family of the spectrally staggered combs, see the right column in Fig. 1, and we will show below that the non-degenerate OPO states are, in fact, approximations of the staggered combs.
Method of slowly varying amplitudes. Here we assume a condition that is quite common for the frequency conversion experiments in χ (2) microresonators. If a resonator is made to operate close to the μ = 0 phase-matching, i.e., |ε 0 |~κ ζ , then the simultaneous control of the repetition-rate difference between the pump and signal, D 1b − D 1a , is hard to achieve. Therefore, μ(D 1b − D 1a ) easily becomes the dominant frequency scale in Eq. (4), i.e., μ|D 1a − D 1b | ≫ |δ 0b |, κ ζ , |ε 0 |, γ a b 0 . For example, for a bulk-cut (R~1 mm) 22 and integrated (R~100 μm) 19 LiNbO 3 resonators, (D 1b − D 1a )/κ ζ~1 0 3 and~10, respectively. Thus, the above conditions work very well for the former starting from |μ| = 1 and for the latter from |μ|~10. Now, the natural methodological step is to separate the fast and slow time scales in the pump sidebands, where B μ are the slowly varying amplitudes. Then, the μ = 0 part of Eq. (4) becomes and the μ ≠ 0 part is Here, Δ μζ are the auxiliary dimensionless detuning parameters, which include the losses and, hence, are complex-valued. We note that Δ μζ are free from D 1b − D 1a which has been absorbed by the fast oscillating exponents, see Eq. (10).
Integrating the ∂ t B μ equation, while assuming that a μ is a slow function of time, we express the pump sidebands via the signal ones, Substituting Eq. (14) into Eqs. (11) and (12) would make up the Kerr-like nonlinear terms. These terms represent the so-called cascaded Kerr nonlinearity 49 , which is, however, negligible in the leading order, because it scales inversely with μ(D 1b − D 1a ). Hence, Eqs. (11) and (12), and the whole of the master system, Eq. (4), simplify to Thus, the pump sidebands, b μ≠0 , play no significant role in the frequency conversion when the repetition-rate difference, μ(D 1b − D 1a ), is large. The latter simply is not featured in Eq. (15). The pumped mode, b 0 , is driven by the sum-frequency processes of the ±μ signal sidebands, which feeds back to the equations for the signal sidebands a ±μ via the b 0 -terms. The absence, in the leading order, of the sum-frequency interaction between the sidebands with |μ 1 | ≠ |μ 2 | hints that the generation of the isolated sideband pairs should be a preferential regime over the broadband frequency combs.
Our approach is different from, e.g., the method when the whole of the high-frequency field is adiabatically eliminated by one way or the other so that the low-frequency field becomes driven by the cascaded Kerr effect, see, e.g., 28,38,49 . The transition from Eq. (4) to Eq. (15) reduces the phase-space dimensionality of the pump field to one, but does not eliminate it entirely, and retains the leading order quadratic nonlinearity.
Solutions and thresholds. The reduced model, Eq. (15), allows examining in detail the properties of the OPO states (this section) and studying their instabilities with respect to each other (next section). In what follows we keep using μ as the running sideband index and ν designates a specific OPO state. Fixing ∂ t a ±ν = ∂ t b 0 = 0, and after some engaging algebra with Eq. (15), the explicit solutions for the non-degenerate OPO states are found in the following form, where, H is the dimensionless pump parameter, and q ν = 2 for ν ≠ 0. Eq. (16) with ν = 0 also covers for the degenerate case, but q 0 = 1. Two possible solutions for the sideband amplitudes are found by taking modulus squared of the equation for e iϕ ν , The recent microresonator theories 19,45,50 reported the doublevalued solutions for the degenerate case. The non-degenerate case was also considered in 19 , but only for the zero detunings and exact phase-matching, while Eq. (16) and stability analysis that follow depend critically on the full account and ability to vary both of those flexibly.
The bifurcation points from the no-OPO to OPO regimes, i.e., Benjamin-Feir instabilities, are conditioned by the zeros of the sideband powers, The laser power W ν at the OPO threshold is then calculated from W ν ¼ 8πκ 2 a H 2 ν =F γ 2 a . The left column in Fig. 2 shows W ¼ W ν vs δ 0b , for the negative, zero, and positive phase-mismatch ε 0 . The minimum of W near δ 0b = 0 exists in all three cases. For ε 0 ≤ 0, this minimum is provided by the ν = 0 mode, i.e., the no-OPO state losses its stability to the degenerate OPO first. The second minimum of W at δ 0b = − ε 0 , i.e., δ 0a = 0 (see Eq. (6)), also happens for ν = 0.
The top row of Fig. 3 shows the results of numerical simulations of Eq. (4) across the regions of the unstable no-OPO state for the relatively small input powers marked with the black dashed lines in the first column of Fig. 2. The ε 0 = 0 case in Fig. 3a demonstrates how the degenerate OPO is first excited from the no-OPO state and then switches, in a cascaded manner, to the non-degenerate OPOs. Figure 3b shows how this scenario can happen twice for ε 0 < 0. We note that the two cascades in Fig. 3b occur in the reverse order. The right-to-left cascade involves oscillations of the μ = 0 and other modes (breathing comb), while all the left-to-right cascades in Fig. 3a-c appear as the direct transitions between the neighbouring non-degenerate OPO states, as it is expected to happen in the Eckhaus instability scenario. We recall that we consider the near-phase-matching between the μ = 0 modes in the pump and signal fields. It leads to For ε 0 > 0, the minimum of W is also found at δ 0b ≈ 0, but now it happens for ν ≠ 0, see Fig. 2b. Thus, here, the no-OPO state transits directly to the non-degenerate regime. Fixing δ 0b = 0 in Eq. (17) gives 16H 2 ν ¼ ε 2 ν =κ 2 a þ 1, where ε ν = ε 0 + ν 2 D 2a . Hence, the minimum threshold power, H 2 b ν ¼ 1=16, for the Benjamin-Feir instabilities is achieved at the exact phase-matching, Figure 3c demonstrates how the below threshold OPO switches directly to the non-degenerate state and passes through the cascade of νs.
Solving H 2 ¼ H 2 ν , Eqs. (18), one could find either two or four real values of δ 0b where the ±ν sidebands bifurcate from zero, cf., Fig. 3d, f with Fig. 3e. In the left column of Fig. 2, the different colours mark the parts of the H 2 ¼ H 2 ν thresholds where either ja þ ν j 2 ¼ 0 (red) or ja À ν j 2 ¼ 0 (blue). The points of transition between the two colours are found by setting Re ðΔ νa Δ 0b Þ ¼ 0: Tuning the pump laser across the red boundaries and moving into the instability interval leads to ja þ ν j 2 gradually increasing from zero, which corresponds to the soft-excitation regime (supercritical bifurcation), see Fig. 3d. Entering the instability tongue across the blue boundary leads to ja þ ν j 2 popping out stepwise and ja À ν j 2 bifurcating from zero sub-critically (hard excitation), see Fig. 3e, f.
The right column in Fig. 2 shows the existence and stability ranges of the degenerate, ν = 0, OPO states. They bifurcate from the no-OPO state super-critically along the red line and sub-critically from the blue line. The parameter range where the ja þ 0 j 2 and ja À 0 j 2 solutions coexist is located between the blue and dashed-grey lines. Generalising for arbitrary ν, the dashed-grey boundary is found from Eckhaus instabilities and OPO tuning. While the top row in Fig. 3 shows the sequential switching between the modes with different numbers and intervals of the pump detuning that select a particular sideband pair, the bottom row shows the changes of the sideband amplitudes computed from Eq. (4) and maps them on the analytical solutions for the OPO states. Apart from the oscillatory instabilities around δ 0b /κ b ≈ 0.5 in Fig. 3b, e, all the insatiabilities of a given OPO state converge to the nearby stable one. In other words, these instabilities lead to the ν → ν + 1 swaps and, hence, to the change of the spatial period and frequency of the waveform in the resonator, i.e., these are the discrete spectrum Eckhaus instabilities. By taking the OPO state with an arbitrary ν ⩾ 0, adding small perturbationsâ μ ðtÞ, where |μ| ≠ ν, and linearising the slowly varying amplitude model (15) we find where b 0 is a function of ν in accordance with Eq. (16). Solving Eq. (21) withâ μ ðtÞ ¼ã μ expftλ ν;μ g andâ Ã Àμ ðtÞ ¼ã Àμ expftλ ν;μ g  17): ja þ ν j 2 þ ja þ Àν j 2 (red) and ja À ν j 2 þ ja À Àν j 2 (blue). Black circles show the numerically computed sideband powers corresponding to the data in the top row. The detuning range around δ 0b ≈ 0.5κ b in the middle column corresponds to the breather states.
yields a set of the sideband-pair growth rates, Thus, Eq. (22) describe the growth rates of the Eckhaus instabilities in the microresonator OPO, i.e., destabilization of the non-degenerate OPO state corresponding to the ±ν sideband pair through the excitation of any other pair ±μ. The instability threshold is reached when the curly brackets become zero, while the generation of the ±ν sideband pair is stable if the pump frequency is tuned to provide δ 0b such that δ 2 νa ⩽ δ 2 μa . The oscillatory bifurcations, i.e., the birth of breathers highlighted by the white lines in Fig. 2, correspond to ν = μ, and, therefore, are exempt from the above theory.
To find conditions of the switching from one sideband pair to the next, we set μ = ν ± 1. For D 2a < 0 (normal dispersion), the interval of stable generation of the ±ν sideband pair is found as The instability boundaries, λ ν,ν±1 = 0, as given by the above condition, are shown in Figs. 3 and 4a by the white dashed vertical lines and they match perfectly with the transitions found from the numerical modelling of Eq. (4).
The detuning corresponding to the midpoint of every step on the Eckhaus ladder, i.e., the left plus right limits divided by two, is δ 0b ¼ Àε 0 À D 2a ðν 2 þ 1 2 Þ. Thus, the latter reproduces the characteristic parabolic shape of δ 0b vs ν in Fig. 4a. Recalling Eqs. (2) and (3), one can show that the same parabola is described by Thus, the analysis done so far has demonstrated that a sequence of the OPO regimes achieved by the scan of the pump frequency follows the steps of the Eckhaus instabilities ladder and the phasematching condition, cf., Eqs. (3) and (24). It also reveals that the size of the Eckhaus ladder's steps and the tuning curve's discreetness are controlled by the second-order dispersion. We shell note that the applicability of Eq. (22) extends beyond the second-order expansion applied for ω μ , Eq. (2), while Eq. (24) relies on it. Therefore, the microresonators with dispersion dominated by the higher orders could be an interesting case to consider in future. Figure 4a shows a sequence of the OPO transitions for the power which is much higher than in Fig. 3 (cf. black and white lines in Fig. 2b). The Eckhaus ladder in Fig. 4a is shown up to ν = 30, but extends up to ν = ν max = 120. The data in the top row of Fig. 3 also demonstrate that, for ε 0 > 0, the parabolic tuning curve is also cut at its minimum. The power reduction, from the levels in Fig. 4a to the ones in Fig. 3, does not change the shape of the parabola but reduces its extent in ν and δ 0b . Since the value of ν max sets the practical limit for the OPO tunability, it is now mandatory to apply our theory to elaborate this point further.
First, we should recall that the data in Fig. 3 have confirmed that Eqs. (16) and (17) provide an excellent approximation for the OPO states. Figure 4d includes the much weaker, practically negligible, part of the spectrum and, thereby, explicitly reveals the degree of accuracy of Eqs. (16) and (17). It also uncovers that the non-degenerate OPO states belong to the family of staggered combs. It follows from Eq. (16) that the power of the signal sidebands, a ±ν , depends on the laser power, H 2 $ W, and the power of b 0 does not, and only its phase does. Indeed, This is a reason why H 2 does not explicitly enter the Eckhaus instability rate, see Eq. (22). However, the limits of existence of ja ± ν j 2 ⩾ 0, and, hence, of b ± 0 in Eq. (21) do critically depend on H 2 as is given by Eq. (20).
Substituting the detunings δ 0b = −ε ν (corresponding to the numbered tips in the left column of Fig. 2 From here one can work out the power dependencies of the minimal, ν min vs W, and maximal, ν max vs W, mode numbers corresponding to the ±ν OPO states. Recalling that the dispersion is normal, D 2a < 0, the maximal number is The ν min is conditioned by the sign of ε 0 . It is zero for ε 0 ⩽ 0 and ν 2 min % ðÀε 0 þ 2κ b HÞ=D 2a if ε 0 > 0, see Fig. 3. The plots of ν max , ν min vs the laser power are shown in Fig. 5.  4), spectra of the optical parametric oscillator (OPO) signal vs pump detuning, δ 0b , for the pump laser power W ¼ 3 mW (4.8d Bm) and ε 0 = 5κ b , see the horizontal white dashed line in Fig. 2a. The white vertical dashed lines map the sequence of the Eckhaus instability thresholds, λ ν,ν±1 = 0, see Eqs. (22) and (23). The colour scale shows the sideband power in dBm. b-d show the details of the spectra of the pump (green) and OPO signal (red) corresponding to the dashed green lines in (a). Spectra in (c) and (d) are the staggered combs. d corresponds to the non-degenerate OPO state that is very well approximated by Eq. (17). The vanishing power levels have been included in (d) in order to demonstrate that the non-degenerate OPO states belong to the staggered comb family.
For the whispering gallery and integrated LiNbO 3 resonators, the mode number change by one corresponds to the wavelength step~0.1 nm 22 and 2 nm 20 , respectively. Hence ν max = 50 corresponds to the wavelength difference between the signal and idler 10 and 200 nm, respectively. The pump-frequency tuning data in Ref. 20 report up to 200-nm signal-idler separations. The experimental reports of the parabolic pumpfrequency tuning curves in the microresonator [18][19][20] , fibre-loop 36 , and bow-tie 38 OPOs do not contain the data sufficient for making a comparison with the reported here power-dependent parabola cut-off at ν max and ν min . Therefore, further research into this problem is necessary.

Conclusions
We have demonstrated that the degenerate OPO transiting to the ladder of the non-degenerate OPO states via a sequence of Eckhaus instabilities is a universal OPO tuning scenario if the repetition-rate difference between the signal and pump fields, |D 1b − D 1a |, dominates over all other frequency scales. The emerging combs are most typically staggered, with the two dominant sidebands in the signal field, which corresponds to the non-degenerate OPO regime, see Figs. 1, 3 and 4. For the highfinesse whispering gallery resonator used as an example in our work, this scenario covers the span of the pump laser powers from nano to milli Watts. The more complex frequency combs start to appear for the pump powers above 3mW, see δ 0b ≈ −20κ b in Fig. 4a, b. Apart from increasing the pump power, the ways to facilitate the microresonator OPOs to operate in the nonstaggered comb could be, e.g., using resonators with the bettermatched repetition rates taken relative to the losses, i.e., having smaller |D 1b − D 1a |/κ a,b , and using pumping into the modes with the odd numbers 8,45,49,50 .
The slowly varying amplitude model, Eq. (15), that we derived and employed here was the key step that has allowed us to solve the problem of switching between the OPO operating in the ±ν sideband pair to any other mode pair ±μ, where ν ≠ μ. This model relies on the aforementioned dominance of the repetition-rate difference, |D 1b − D 1a |, which drives the fast oscillations of the pump sidebands, see Eq. (10). It also reveals that the growth rate of the Eckhaus instability does not depend on |D 1b − D 1a | in the leading order, see Eq. (22). Regarding the quadratic nonlinearity, the slowly varying amplitude model highlights the dominant role of the terms responsible for the nonlinear mixing of the central mode of the pump field with all the signal sidebands.
Our methodology has also let us derive the equation for the maximally allowed sideband order generated by the pumpfrequency tuning method in the microresonator OPOs, which is easy to apply and use as a practical guideline for the parametric down-conversion experiments. All our analytical results are in excellent agreement with the numerical modelling of the full coupled-mode system, see Eq. (4).

Methods
Numerical integration of the coupled-mode model, see Eq. (4), has been performed using a fourth-order Runge-Kutta method. The procedure was optimised by applying the fast Fourier transform algorithm to the nonlinear terms in the real space to compute the sums in Eq. (4), see Ref. 44 for details. The data shown in Figs. 3 and 4 were obtained for 256 optical modes in the 'a' and 'b' fields, i.e., for μ = −127, … 128. The initial conditions were the no-OPO state in the μ = 0 modes and random perturbations in the rest of the spectrum.

Data availability
The data supporting the findings of this study are available from the authors upon reasonable request.

Code availability
The codes used for this study are available from the first author (dp710@bath.ac.uk) upon reasonable request. Fig. 5 Maximal and minimal optical parametric oscillator (OPO) sideband orders achievable by the pump-frequency scan. The pump power dependencies of the maximal, ν max , and minimal, ν min , mode numbers corresponding to the ±ν OPO states that can be generated by the pumpfrequency tuning method. ν min is different from zero only for ε 0 > 0. b ν % ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Àε 0 =D 2a p corresponds to the lowest possible excitation threshold for ε 0 > 0, see text between Eqs. (18) and (19). Dispersion is normal, D 2a < 0, see Fig. 1 for other parameters.