Self-Starting Soliton–Comb Regimes in χ ( 2 ) Microresonators

: The discovery of stable and broad frequency combs in monochromatically pumped high-Q optical Kerr microresonators caused by the generation of temporal solitons can be regarded as one of the major breakthroughs in nonlinear optics during the last two decades. The transfer of the soliton– comb concept to χ ( 2 ) microresonators promises lowering of the pump power, new operation regimes, and entering of new spectral ranges; scientiﬁcally, it is a big challenge. Here we represent an overview of stable and accessible soliton–comb regimes in monochromatically pumped χ ( 2 ) microresonators discovered during the last several years. The main stress is made on lithium niobate-based resonators. This overview pretends to be rather simple, complete, and comprehensive: it incorporates the main factors affecting the soliton–comb generation, such as the choice of the pumping scheme (pumping to the ﬁrst or second harmonic), the choice of the phase matching scheme (natural or artiﬁcial), the effects of the temporal walk off and dispersion coefﬁcients, and also the inﬂuence of frequency detunings and Q-factors. Most of the discovered nonlinear regimes are self-starting—they can be accessed from noise upon a not very abrupt increase in the pump power. The soliton–comb generation scenarios are not universal—they can be realized only under proper combinations of the above-mentioned factors. We indicate what kind of restrictions on the experimental conditions have to be imposed to obtain the soliton–comb generation.


Introduction
The generation of broad frequency combs, i.e., of long sequences of equidistant coherent light lines, in monochromatically pumped high-Q Kerr microresonators [1][2][3][4][5] can be regarded as one of the main discoveries in nonlinear optics during the last two decades. Remarkably, the spectral positions of these lines do not generally coincide with (or are not close to) the discrete resonator modes because of the dispersion of the latter. This is why the interpretation of the combs as a sequential parametric excitation of new resonator modes is rather fruitless. The right physical picture of the frequency comb in question is proven to be the generation of a spatially narrow soliton propagating along the resonator rim with a constant velocity v 0 [5]. Fourier expansion of the light field employing its periodicity over the circumference ensures the equidistance of the light lines with the frequency spacing v 0 /R, where R is the major resonator radius. Thus, the comb and soliton are two physical terms expressing different properties of the same nonlinear light state. Most probably, the soliton nature of broad frequency combs has no alternative.
The soliton-comb concept is nowadays generally accepted for χ (3) resonators. It is applicable not only to glass-based Kerr microresonators, but also to resonators made of χ (3) crystals [5][6][7][8][9][10]. One can speak of the corresponding comb technology. The amount of publications on the χ (3) combs can be estimated as ∼10 2 . Note that the dissipative solitons in question are essentially different from (and more complicated than) the so-called conservative solitons associated typically with the nonlinear Schrödinger equation [11,12]. The dissipative solitons balance not only the dispersion broadening and nonlinear compression but also the gain owing to external pumping and the resonator losses.  (2) comb generation; R and ϕ are the major resonator radius and the azimuth angle. The output comb spectra are due to a dual FH-SH soliton propagating along the resonator rim with velocity v 0 . The red spots indicate localization of the resonator modes near the rim.
Axially symmetric χ (2) microresonators are widespread [17,18,[20][21][22][23]. Their linear properties are largely common with those of the χ (3) resonators [17,[24][25][26]. Each resonator can be characterized by a major radius R, see Figure 1. The modes (whispering gallery modes) are localized near the resonator rim, and the rim shapes can be different. Each mode can be labeled by an integer azimuth number m and, additionally, by two transverse modal numbers. All modal functions are proportional to exp(imϕ), where ϕ is the azimuth angle. While exact analytical solutions exist only for spherical resonators, approximate methods to calculate the modal frequencies, functions, and nonlinear interactions are well developed [24,[27][28][29][30][31]. This is also relevant to the methods for coupling light into and out of the resonator. Owing to the presence of a coupler (tapered fiber, prism, etc.), the inverse total (loaded) quality factor Q −1 consists of internal and coupling contributions.
The most common lithium niobate (LN) based χ (2) resonators are discs with the major radius R ∼ 1 mm, the minor radius r R, and quality factors Q = (10 7 -10 8 ) [17,18]. The resonator modes can be viewed as quasiplane waves with frequencies ω ≡ 2πc/λ = k m c/n, discrete wavenumbers k m = m/R, and an effective refractive index n depending on the modal numbers, R, and r. For R > 1 mm, n is close to the bulk index n b (λ); the difference n − n b characterizes the effects of the geometric dispersion. In the optical range, the azimuth number m = 2πnR/λ is ∼10 4 . The frequency distance δω = c/nR between neighboring modes with δm = 1 is much larger than the line width. As LN is a uniaxial optical material, modes with ordinary (o) and extraordinary (e) polarizations have to be distinguished. In addition to lithium niobate, other polar crystals such as lithium tantalate, beta-lithiun borate, and aluminium nitride can be used to manufacture high-quality, Q ≈ (10 6 -10 8 ), microresonators for the visible and UV spectral ranges [18,32,33].
Investigations of the SH generation and parametric oscillation in χ (2) microresonators are also known [17,18]. In LN resonators, the PM condition for the SH generation, ω 2m = 2ω m , can be fulfilled naturally at the FH wavelength λ 1 1064 nm with the use of o-and e-polarized FH and SH modes [14]. The parametric oscillation for this natural (birefringent) PM corresponds to the SH wavelength λ 2 = λ 1 /2 532 nm. Here and later, the subscripts 1 and 2 refer to the FH and SH ranges, respectively. The employment of the radial poling in LN resonators [15][16][17][18] enables one to fulfill the quasi-PM conditions at practically any wavelength on demand via a controllable breaking of the axial symmetry. It is also necessary to mention that, because of very small resonator line widths, fine tuning means (such as the temperature or geometric tuning) have to be used to control the linear and nonlinear phenomena [17,18,25,34].
Both the χ (3) and χ (2) combs originate from parametric instabilities inside the resonators. This is why the continuous-wave monochromatic pump is able to generate temporal solitons. The threshold pump power for the instability Pth is largely controlled by the relevant modal Q factors. Roughly speaking, in the χ (3) and χ (2) cases, Pth is proportional to Q −2 and Q −3 , respectively. For Q 10 7 , the χ (2) nonlinearity is typically dominating, and Pth belongs to the µW range [17,18].
Initial experimental and theoretical attempts to obtain χ (2) combs dealt with bulk systems [35][36][37][38][39][40][41] that are quite specific in physics and relevant to much higher threshold pump powers. While some tens of comb lines were observed, neither solitons nor broad comb spectra were reported. Recently, the first experimental efforts to generate the χ (2) combs in microresonators have been made [42][43][44][45]. No soliton-comb regimes have been observed so far to the best of our knowledge. Furthermore, we can claim that this failure is in agreement with our theoretical results on the necessary conditions for the realization of such regimes. Therefore, we see the main tasks of this mini-review in (i) the determination of the conditions for the χ (2) soliton-comb generation and (ii) in the predictions of the main properties of the soliton-comb states. Additionally, we will indicate where and why the necessary conditions were violated in the known experiments. Lastly, we will compare our findings with the known early predictions for the χ (2) soliton-comb states [46][47][48][49]. As far as possible, the cases considered will be applied to experimentally realizable and important situations.
A considerable part of the literature on nonlinear phenomena in microresonators, see, e.g., [4,17,20,44], extensively employs the quantum-mechanical terminology and notation deeply within the classical domain, where the applicability of the classical Maxwell equations is beyond any doubt. Being well developed and attracting, the quantum-mechanical formalism here is out of place and confusing. We claim that it can easily be replaced by the classical Hamiltonian formalism [50,51]. Keeping all merits of the quantum formalism, it is fully adequate to the situation.

Periodic and Antiperiodic States
First, let the pump frequency ω p be close to a modal frequency and the FH-SH phase matching be natural, i.e., the azimuth symmetry takes place. Then, two variants have to be considered, namely pumping into a FH and into a SH mode. In the FH pumping case, the PM condition ω(2m 0 1 ) 2ω(m 0 1 ) can be satisfied with a high accuracy only for a single azimuth number m 0 1 . It is implied that the frequency detunings ω p − ω(m 0 1 ) and ω(2m 0 1 ) − 2ω(m 0 1 ) are much smaller in the values than the intermodal distance c/nR. The light electric field E(ϕ, t) can be represented as where E 1 (ϕ, t) and E 2 (ϕ, t) are slowly varying (and generally vectorial) FH and SH amplitudes (envelopes). Obviously, these amplitudes are 2π-periodic in ϕ, and the SH azimuth number m 0 2 = 2m 0 1 is even. In the SH pumping case, the situation is different, as illustrated by Figure 2. The above PM condition can be fulfilled only if the azimuth number of the pumped mode m 0 2 is even ( Figure 2a). If m 0 2 is odd (Figure 2b), two FH modes (with even and odd azimuth numbers) must be excited. Here, the simplest PM conditions are ω m 0 Let us use (regardless of the parity of m 0 2 ) the ansatz for the light field. Then, the first exponential factor is 2π periodic for an even m 0 2 and 2π antiperiodic for an odd m 0 2 . As the true light field E(ϕ, t) is 2π periodic, we must accept that the FH amplitude E 1 (ϕ, t) is antiperiodic for the odd m 0 2 ; i.e., E 1 (ϕ, t) = −E 1 (ϕ ± 2π, t). The SH amplitude is still 2π periodic, E 2 (ϕ, t) = E 2 (ϕ ± 2π, t). A dual FH-SH state with an antiperiodic amplitude E 1 and periodic amplitude E 2 is named the antiperiodic state [52,53]. The antiperiodic sates (solutions) are topologically different from the conventional periodic states. No temporal evolution can make the conversion between these two states. To the best of our knowledge, the separation of the light states into periodic and antiperiodic ones is original, and it was never implemented for bulk χ (2) systems. We stress that the antiperiodic states are relevant to pumping into an odd SH mode. Otherwise, we have a periodic dual state, such that E 1 (ϕ, t) = E 1 (ϕ ± 2π, t) and E 2 (ϕ, t) = E 2 (ϕ ± 2π, t). The excitation of periodic (a) and antiperiodic (b) states for SH pumping into even and odd resonator modes. The primary FH azimuth number m 0 1 = m 0 2 /2 is an integer in (a) and semi-integer in (b). Side harmonics arise automatically in (b) above the threshold. Using the above vectorial notation, we wanted to stress the generality of the periodic and antiperiodic light states for the χ (2) nonlinearity. While the FH and SH resonator modes can be different in polarization, one can always deal with scalar FH and SH amplitudes denoted in what follows as F and S.

Equations for the Modal Amplitudes
Now, we write down nonlinear equations for the complex modal amplitudes F m 1 and S m 2 relevant to the FH and SH frequency domains. It is assumed, as in all known theoretical studies, that the modal polarization and the transverse modal numbers stay the same within each of these domains. We start from the phase matching case relevant to an even azimuth number m 0 2 = 2m 0 1 in the absence of radial poling (full azimuth symmetry). This takes place for the FH pumping and is also applicable to the SH pumping into an even mode, see Figure 2a. It is useful to count the FH and SH frequencies from ω 1 and ω 2 . In accordance with Equations (1) and (2), ω 1 = ω p , ω 2 = 2ω p in the FH pumping case, while ω 1 = ω p /2, ω 2 = ω p in the case of SH pumping. Furthermore, instead of m 1 and m 2 , we introduce the more convenient FH and SH modal numbers j = m 1 − m 0 1 = 0, ±1, . . . and l = m 2 − m 0 2 = 0, ±1, . . .. With these preliminaries, we have Here, the dots indicate the differentiation in time; Z 1,j = ω j − ω 1 − iγ j and Z 2,l = ω l − ω 2 − iγ l , γ j,l = ω j,l /2Q j,l are the modal decay rates; µ and h 1,2 are the nonlinear coupling coefficient and the internal pump amplitudes, respectively; and δ j0,l0 are the Kronecker deltas. For the FH and SH pumping, we have h 2 = 0 and h 1 = 0, respectively. The modal amplitudes F j and S l are defined such that ω j |F j | 2 and ω l |S l | 2 are the modal energies. This definition corresponds to the classical Hamiltonian formalism [50]. The dimensionless quantities |F j | 2 /h 1 and |S l | 2 /h 1 are nothing but the number of modal FH and SH light quanta in the classical limit. The presence of a factor of 2 in set (3) reflects the Hamiltonian nature of the nonlinear interaction-the total energy is conserved for h 1,2 = 0 and γ j,l = 0.
The coefficients µ and h 1,2 can be determined in a quite general manner via modeling the modes and their coupling in and out of the resonator [29,31]: where d is the relevant component of the conventional susceptibility tensor, σ eff is the effective interaction cross-section, P 0 is the pump power,Q −1 is the coupling contribution to Q −1 , and n n o n e is a representative refraction index. Further details can be found in Section S1 of the Supplementary Materials. The rate coefficients Z 1,j and Z 2,l can be expressed by the group velocities v 1,2 and dispersion coefficients v 1,2 if we use the expansion ω(k) = ω(k 0 ) + vδk + v δk 2 /2 with δk = k − k 0 near the points k 0 = m 0 1 /R and m 0 2 /R: The parameters v and v are expressed by the refractive index n(λ): v = c/(n − λn ) and v = −cλ 3 n /2π(n − λn ) 2 , where the prime indicates the differentiation in λ. Knowledge of the dependences n o,e (λ), with or without the geometric dispersion, enables one to calculate v 1,2 and v 1,2 as functions of λ 2 and λ 1 = 2λ 2 . Further details can be found in Section S2 of the Supplementary Materials. Expressions for detunings ∆ 1,2 depend on the pumping scheme. For the FH pumping, Consider now the remaining case of SH pumping into a mode with an odd azimuth number m 0 2 when two FH modes with even and odd azimuth numbers enter the phase matching condition, see also Figure 2b. The SH modes can be numerated again with the integer l = m 2 − m 0 2 = 0, ±1, . . .. The only way to symmetrically numerate the FH modes is to introduce a semi-integer number j = ±1/2, ±3/2, . . . such that m 1 = m 0 2 /2 + j. The form of sets (3) and (5) remains the same, and the only difference is in the expression for detuning ∆ 1 :

The Impact of Radial Poling
The radial poling enables one to adjust the phase matching to practically any desirable spectral point with no effect on the linear properties. In polar χ (2) materials, such as LN, the largest components of the quadratic susceptibility tensor, d 13 and d 33 , change sign under the inversion of the spontaneous polarization. The periodic radial poling results in the ϕ periodic sign changing for any of these components, as illustrated schematically by If N is the number of poling periods, the nonzero components d s correspond to s = 0, ±N , ±2N , . . .. In the case of the ± symmetric radial domain structure (50% duty cycle ratio), which is the most suitable for quasiphase matching, only the Fourier harmonics with odd ratios s/N = ±1, ±3, . . . are nonzero. For these harmonics |d s | = 2d 0 N /π|s|, see also Figure 3b. The reduction factor |d s |/d 0 decreases with an increasing |s|/N but remains comparable with 1 for |s| = N . Typically N 1, and the neighboring peaks are well separated from each other.  The employment of the largest first Fourier harmonics of d(ϕ) for quasiphase matching corresponds to the PM condition 2ω m 0 1 = ω 2m 0 1 ±N for SH generation. The sign "plus" is relevant to the most typical case of decreasing wavelength dependence of the refraction index n(λ) [32]. We see that the SH azimuth number m 0 2 = 2m 0 1 ± N is even for an even alternation number N and odd for an odd number N . In the second case, SH pumping leads to the excitation of periodic states. In the case of the odd combination m 0 2 ∓ N , the above PM condition cannot be fulfilled. However, we can fulfill the PM condition ω m 0 1 + ω m 0 1 +1 = ω 2m 0 1 ±N , which leads to the excitation of the antiperiodic states. Thus, the quasi-PM by means of the radial poling results in a controllable shift of the FH azimuth number m 0 1 . The set (3) for the FH and SH amplitudes results in them staying unchanged. Note that microresonators made of commercial PPLN wafers [20] are not appropriate for the spectral adjustment in question. The spectrum |d s | includes in this case many significant peaks |d N ±1 |, |d N ±2 |, . . . in addition to |d N |, which leads to unwanted interference of many cascading χ (2) processes.
The possibility to adjust the phase matching to any spectral range raises the question about the specific goals and properties of the radial poling. Figure 4 shows the dependences v 1,2 (λ) and v 1,2 (λ) for the e-polarized modes in a LN-based resonator with R = 1.5 mm. The effects of geometric dispersion are very weak here. The zero walk-off point corresponds to λ c 2 1349 nm, where the dispersion coefficients v 1 and v 2 are opposite in signs. The vicinity of this point, as will be shown below, is especially attractive for the soliton-comb generation. The number of periods N necessary to achieve quasi-PM at this point can be estimated as 284. Further details can be found in Section S3 of the Supplementary Materials.

Fundamental z, t Representation
The primary set (3) of the ordinary differential equations is convenient for numerical simulations. However, many fundamental features of the soliton-comb states become more evident if we transform it to a set of partial differential equations. This can be performed in terms of the FH and SH amplitudes F(ϕ, t) and S(ϕ, t) defined by the Fourier expansions For integer values of j, l, this leads automatically to the periodic states, while for integer l and semi-integer j, we arrive at an antiperiodic state. Using Equations (3), (5) and (7) and introducing the rim coordinate z = Rϕ, one can make sure that the periodic and antiperiodic states obey the same set of nonlinear equations for any phase-matching scheme [52,53]: where Ω 1,2 = ∆ 1,2 − iγ 1,2 . At perfect resonant pumping and phase matching, we have ∆ 1,2 = 0. All parameters entering set (8) are known or controllable in experiment. For the χ (2) solitons, set (8) plays a role that is similar to the role of the nonlinear Lugiato-Lefever equation for the χ (3) soliton-comb states [4,5]. However, these two cases are fundamentally different. For the periodic states, equations analogous to set (8) can be found in [46][47][48]. Consider some important properties and useful transformations of set (8): circulating with a common velocity v 0 without shape changes. Solitons belong to this class of solutions. Velocity v 0 has to be determined simultaneously with the shape of the envelopes. -Set (8) is written for a static coordinate frame. It is practical to rewrite it for a coordinate frame moving with velocity v 1 . To conduct this, it is sufficient to drop the term v 1 ∂ z in the first equation and replace v 2 by −v 12 = v 2 − v 1 in the second one. This is assumed from now on. - In a dual steady state propagating with a constant velocity v 0 , the modal amplitudes F j and S l oscillate in time as exp(−ijv 0 t/R) and exp(−ilv 0 t/R). This elementary property is useful to control the establishment of dual steady states in numerical simulations; see Section 4. -Attempts to rewrite set (8) for the antiperiodic states using the replacements F(z, t) = F(z, t) exp(±iz/2R) lead to restoration of the periodicity but also to an explicit zdependence of the right-hand sides (nonautonomous system). This is inappropriate.
The structure of Equation (8) gives a hint about the normalization of parameters relevant to the linear properties of the system. Such a normalization enables one to decrease the number of variable parameters. In what follows, we will use the following normalized detunings, walk off, and dispersion parameters δ 1,2 , α, and β 1,2 : As a rule, the actual values of the detunings are restricted by the inequalities |δ 1,2 | 1. The actual values of α and β 1,2 substantially depend on the group velocity difference v 12 = v 1 − v 2 , the dispersion coefficients v 1,2 , and the decay rates γ 1,2 .
For LN resonators, typical values of the decay rates are ∼10 7 , such that |β 1,2 | 1. On the contrary, we typically have |α| 1. Exceptions to this rule are possible near the zero walk-off point λ 2 1349 nm. The decay rates γ 1 and γ 2 are typically not strongly different, such that the assumption γ 1,2 = γ can be used sometimes to simplify the modeling.
Generally speaking, the values of F, S, and h 1,2 also need some normalization. It is, however, different in the FH and SH pumping cases. We will provide the relevant information when necessary.
One more important property of set (8) is the presence of spatially uniform steadystate background solutionsF = 0,S = 0 [54]. In the steady state, F 0 =F and S 0 =S. Obviously, such a dual background is possible only for the periodic states. The dual background is closely related both to the properties of the solitons and to the thresholds of the instabilities leading to soliton-comb formation. Omitting all derivatives in (8), one can easily come to a closed algebraic equation for |F| 2 that is, however, different for the FH and SH pumping cases. We have |S| = µ|F| 2 /|Ω 2 | and |Ω 1 |/2µ for the FH and SH pumping cases, respectively. The dependence of |F| 2 on the pump amplitudes h 1,2 is presented graphically in Figure 5. In both pumping cases, the dependence of |F| 2 (h 1,2 ) is controlled by the phase Φ = arg(Ω 1 Ω 2 ), such that multivalued solutions exist for sufficiently small |Φ|. These dependences can be interpreted in the terms of the bifurcation theory [55]. This theory cannot, however, provide useful insights into the properties of light states well above the instability thresholds because of a very large number of the FH and SH light harmonics. Importantly, the dual background withF,S = 0 always exists in the FH pumping case, see Figure 5a. At the SH pumping, the situation is different. For sufficiently small values of h 2 , there is a single background state, such thatF = 0 and |S| = h 2 /|Ω 2 | = 0, see Figure 5b. This difference has dramatic consequences for the dual comb-soliton states, see Sections 4-6. Note also that Figure 5 gives a clear hint about the normalization of the amplitudes F, S and the pump amplitudes h 1,2 in the FH and SH pumping cases.
Lastly, we mention two limiting cases for set (8): - In the limit h 1,2 → 0, γ 1,2 → 0 (no gain, no modal decay), we proceed to the conservative case reviewed in [19]. A number of exact soliton solutions are known here, but they are far from the subject of our study. - In the dissipativeless driven limit γ 1,2 → 0 at ∆ 1,2 = 0 and h 1,2 = 0, one can also obtain exact soliton solutions [56]. It turns out, however, that they are all unstable, i.e., they cannot be realized.

Instability Thresholds
As is well known, the generation of second harmonic is thresholdless-it corresponds just to forced oscillation caused by the χ (2) nonlinearity. This thresholdless process becomes merely more efficient in the presence of FH-SH phase matching. On the contrary, the excitation of FH harmonics in the presence of SH modes is always a coherent threshold process. This explains a big difference in the properties of the parametric excitation for the FH and SH pumping schemes in the microresonators [54].
The simplest and most efficient scheme is the SH pumping scheme. At sufficiently small pump powers (small values of h 2 at h 1 = 0), we deal simply with the linear excitation of the pumped mode m 0 2 . Its steady-state amplitude is S 0 =S = h 2 /(γ 2 + i∆ 2 ), and detuning ∆ 2 characterizes just the distance to the linear resonance. To see the effect of S 0 on the FH harmonics, we consider the first of Equation (3) with S l = S 0 δ l,0 and h 1 = 0. The quantities F j and F * −j with an arbitrary j obviously obey two linear differential equations. Can they have a nonzero steady-state solution? Using the zero determinant solvability condition for linear algebraic equations and Equation (9) for the normalized parameters, we obtain the threshold condition in the form [52][53][54] where η 2 = 2µh 2 /γ 1 γ 2 is the normalized pump amplitude. The same threshold equation can be obtained using a more general approach; namely, one can set F j , F * −j ∝ exp(−iνt) to determine the complex characteristic exponent ν = ν + iν via a linear stability analysis. Then, one can make sure that ν > 0 for η 2 > η th 2,j , while ν < 0 for η 2 < η th 2,j . Thus, small amplitudes F ±j (t) decay below the threshold and grow exponentially above it. This fully justifies the use of the threshold Equation (10). Note that similar relations with different notations are known in the literature.
Consider in some detail the threshold properties given by Equation (10). Remarkably, it includes neither the walk-off parameter α nor the SH dispersion coefficient β 2 . The threshold equation refers to the FH pairs j and −j. The value of j can be an integer or semi-integer depending on the type of solution in question. The thresholds are generally different for different |j|. However, the minimal in δ 1 values of η th 2 are the same for all |j|. The value of δ 1 minimizing the threshold for a pair ±j is δ Using this property, one can excite selectively any j, −j pair on demand. The absolute threshold minimum, [η th 2 ] min = 1, can be achieved indeed at δ 2 = 0, i.e., at the exact resonant excitation of the pumped SH mode. It is clear now that the parameter η 2 can be expressed by the pump power as η 2 = (P /P th min ) 1/2 . Now we turn to the FH pumping scheme. It is relevant not only to the comb generation issue, but also to the SH excitation. This scheme is complex compared to the previous one for fairly simple reasons. The parametric instability requires sufficiently large values of the SH amplitude S 0 =S. They can be achieved, in turn, only for sufficiently large values of the amplitude of the pumped mode F 0 =F(h 1 ). Keeping this fact in mind, we consider what happens if we admit the presence of small perturbations δF(ϕ, t) and δS(ϕ, t), such that F =F + δF and S =S + δS. Substituting these expressions into Equation (8) with h 2 = 0, we see that the SH perturbation δS becomes linearly linked to the FH perturbation via the termFδF, which was absent earlier. As a result, we must deal with a quartet of instead of the former parametric pair F j , F * −j . The characteristic equation for ν becomes a cumbersome four-degree algebraic equation including all linear parameters δ 1,2 , β 1,2 , and α [54]. The threshold properties become complicated, and they cannot be covered by simple relations similar to Equation (10).
Nevertheless, a general overview and illustrations of the threshold properties, which do not pretend to be complete, are possible. We introduce first the normalized pump amplitude η 1 = √ 2µh 1 /γ 1 √ γ 1 γ 2 ; this normalization is different from that relevant to the SH pumping. The presence of parameter α, ranging from 0 to very large values, in the threshold equation raises questions about the effect of walk off. It turns out that the lowest threshold value of η 1 at a zero walk off (α = 0) corresponds to j = 0, i.e., to an internal instability of the dual background [54]. For detunings δ 1,2 = 0, this value is minimal, With increasing |α|, some quartets with j = 0 quickly become the easiest to excite. However, the optimal values of |j| depend significantly on the dispersion coefficients. This is illustrated by Figure 6, which was obtained with the use of numerical methods.  The effect of |α| is strongly different here, and it strongly depends on |j|. The reason is in the competition between the walk off and dispersion terms entering the rate coefficient Z 2,l given by Equation (5). Anyhow, the selective excitation of modes with |j| = 1, 2 is practically impossible. The effects of nonzero detunings of δ 1,2 can also be investigated numerically. They also do not facilitate the selective excitation of harmonics with small values of |j|.
With increasing |α|, the threshold values η th 1,j (α) slowly approach two from above. This range is, however, not favorable for the soliton-comb generation because of a strong suppression of the side SH harmonics with l = 0, see also Sections 5 and 6.

Numerical Methods
We simulated predominantly set (3) in the coordinate frame moving with velocity v 1 by using the conventional fourth-order Runge-Kutta method [57]. Initially, for verification purposes, we also solved set (8) of the partial differential equations by using the step-split Fourier method [58]. The results were practically equivalent. The number of Fourier harmonics taken into account (and the mesh points in ϕ) ranged from 32 to 1024 within each of the FH and SH frequency domains. The accuracy of the calculations was controlled by changing the time step and number of harmonics. Different values of the variable input dimensionless parameters of the system (η 1,2 , δ 1,2 , β 1,2 , α) were used. With harmonics F j (t), S l (t) calculated, one can investigate the temporal evolution of the spatial profiles F(ϕ), S(ϕ) and comb spectra |F j | 2 , |S l | 2 . Additionally, it is possible to verify the establishment of steady states moving with a constant velocity and, if applicable, determine this velocity v 0 .
Three main problems were considered numerically: -We numerically verified the analytical results relevant to the instability thresholds.
To do so, we used initial conditions with very small random complex values of F j (0) and S l (0) to see the initial exponential growth or decay of |F j |(t) for different values of η 1,2 . Coincidence of the analytical and numerical results with a high accuracy was always achieved. -We verified the stability of the dual soliton solutions obtained analytically in [52,56] for the known limiting cases. As the initial conditions, we used Fourier transforms of the relevant analytical expressions for F(ϕ) and S(ϕ). Numerical analyses has shown that all the analytical solutions are unstable: Temporal evolution inevitably leads to irregular spatial profiles and comb spectra. - We tried to generate stable steady-state soliton-comb solutions starting from very small random complex amplitudes F j (0) and S l (0). As a rule, the evolution of our nonlinear systems well above the threshold ultimately leads to highly irregular saturated behavior after a stage of exponential growth. When abruptly switching the pump on, the rough features of this irregular behavior depend also on the choice of the initial conditions. Nevertheless, for many well-defined input parameters, proper pumping schemes, and not very abrupt switching the pump on, temporal evolution leads reliably and uniquely to multiparametric families of comb-soliton states. Below we describe the corresponding adiabatic procedure of growing soliton-comb states as applied (for definiteness) to the SH pumping case.
Let FH harmonics ±j possess the lowest threshold. For η 2 < η th 2,j , we have then the steady state with S 0 = 0 and zero other harmonics. Switching to η 2 slightly above η th 2,j , we obtain a steady state including S 0 and only a few relatively small other essential FH and SH harmonics. After that, the steady-state harmonics are used as new initial conditions at a slightly larger value of η 2 . This procedure is repeated then many times to come to a steady state well above the threshold. The number of FH and SH harmonics taken into account has to be chosen properly. The final state is insensitive to the choice of the initial random amplitudes. The states achieved in this way can be qualified as self-starting steady states, and their stability is doubtless. Let us stress that the achievement of the soliton-comb states with the adiabatic procedure is not guaranteed. However, in many cases, it takes place.
The restriction on the rise time of the pump, τ p , that is necessary to achieve the steady state is not very obligating. It is sufficient to fulfill the inequality γ 1,2 τ p 1, which corresponds typically to τ p 1 µs. Note also that soliton-comb states can be achieved often in a few steps of step-wise increase in the pump power.
The establishment of soliton-comb steady states and the corresponding velocity v 0 can be quantified with a high accuracy by using a correlation method. We use the fact that the Fourier harmonic A j (t) (F j or S j ) known in the coordinate frame moving with velocity v 1 can be recalculated in the frame moving with an arbitrary velocity v * through multiplication by exp[ij(v 1 − v * )t/R]. Let us introduce the dimensionless error parameter where A j refers to a coordinate frame moving with velocity v * , t is the calculation time, and τ is a variable time shift. Obviously, ε(t, τ) turns to zero only when we deal with the steady state and, simultaneously, v * = v. The error parameter calculated for modestly large evolution times, γ 1,2 t 10 3 , and minimized over v * shows extremely small values (ε = 10 −14 -10 −15 ) caused by the numerical noise, see also below. For smaller t, i.e., during the transient stage, it is larger by many orders of magnitude. We thus have a numerical tool to control the establishment of steady states and to precisely determine the velocity v 0 . Figure 7 shows three examples of the temporal evolution of ε(t, τ) relevant to our adiabatic procedure applied to the antiperiodic case at γ 1,2 = γ = 10 7 s −1 and τ = 1/γ; the parameters β 1 = 0.02, β 2 = −0.01, α = 1 are representative for a close vicinity of the zero walk-off point in the LN-based resonators. The harmonics in question are F j . The normalized pump amplitude η 2 sequentially takes 17 values above the threshold: 1.01, 1.1, . . . , and 30. Figure 7 exhibits the evolution of ε after taking the first two and the last values at t 0 , t 1 , and t 17 . The most general feature of this evolution is always the same: it is an almost exponential decrease in ε from initial values of the order of 1 to extremely low values of 10 −14 -10 −15 . Note also some quantitative details: The longest evolution, shown in Figure 7a, starts from modal noise. It takes about (5-7)γ −1 for the harmonics to grow exponentially and reach the saturation. The corresponding evolution time can be identified with the rise time of the near-threshold optical parametric oscillation. The further evolution runs, which start from regular distributions of the FH and SH harmonics achieved during the previous runs, occur substantially faster. Our choice of taking 17 steps of evolution to reach the state with η = 30 is largely arbitrary. In many cases, the same final state can be achieved in a few steps. In the next sections, we will see in detail the spatial and spectral properties of the dual states grown with the use of our adiabatic procedure.

Soliton-Comb States for SH Pumping
Below we report the results on soliton-comb generation obtained by using our adiabatic procedure in the SH pumping case for LN resonators [53,59]. The ϕ-dependences of |F| 2 , |S| 2 , arg[F], and arg[S] and the spectra |F j | 2 and |S l | 2 are of our prime interest. Additionally, we are interested in the dependences of the soliton-comb characteristics on the pump amplitude η 2 , the walk-off parameter α, and other variable parameters of the system. Generally, the soliton position in the resonator depends on the choice of the initial random amplitudes. Since the azimuth angle ϕ can be counted from any point, we can chose its zero point as convenient.
We start our report with the case of quasi-PM via the radial poling. Here, the pump wavelength λ p can take practically any value. In particular, it can be made arbitrary close to the zero walk-off point λ c 2 1349 nm. This case is especially important. For sufficiently small values of δλ = λ p − λ c 2 , the group velocity difference is given by v 12 [10 5 cm/s] −5.9 × δλ [nm] in accordance with Figure 4a. Correspondingly, for |δλ| 10 nm and typical values of R and γ 2 , the walk-off parameter α = v 12 /γ 2 R ranges from 0 to large negative and positive values, while the dispersion coefficients β 1,2 (δλ) stay practically constant and opposite in signs, β 1 > 0, β 2 < 0.
It is relevant to the zero walk off (δλ = 0, α = 0); η 2 = 30; zero detunings (δ 1,2 = 0); R = 1.5 mm; and decay rates γ 1 = 10 7 , γ 2 = 10 8 s −1 (Q 1 3.3 × 10 7 , Q 2 0.7 × 10 7 ). The total number of harmonics taken into account is 1024. Figure 8a and Figure 8d show the central parts of the normalized FH and SH intensity profiles-the intensities stay constant far from the center. It is evident that we are dealing with a narrow dual dark-bright soliton. Far from the soliton core, the normalized values of |F(ϕ)| 2 and |S(ϕ)| 2 ( 29 and 1) nicely correspond to the dual background considered in Section 2. At the first sight, this feature is not consistent with the absence of the dual background for the antiperiodic states. The point is that the amplitude F(ϕ) tends to the opposite left and right constant values far from the center, and these values can be found from Equation (8). The facts that (i) the intensity profiles are even in ϕ and (ii) the FH intensity |F(ϕ)| 2 turns exactly to zero at the soliton center are exclusively due to our zero walk-off assumption.  The FH and SH comb spectra are shown in Figure 8c and Figure 8f. They are symmetric in j and l and very broad. The normalizations of F and S are slightly different and specific for the SH pumping case. On the other hand, they provide the maximum values of |F j | and |S l | of the order of one. The shapes of the FH and SH spectra are notably different. The central peak in Figure 8f is due to the pump. To quantify the comb spectra, we define the number of significant FH and SH lines (N 1 and N 2 ) whose normalized intensities are above the cut-off level 10 −4 . This characteristic is conditional and rather useful; the lines of this strength are typically well above the light noise. In our case, N 1 = 122 and N 2 = 161 correspond to the well developed dual combs.
With increasing η 2 , the dual soliton is progressively narrowing and the comb spectra are getting broader. No deterioration of the soliton-comb states was seen up to the maximal investigated values η max 2 ≈ 10 2 . Despite an apparent simplicity of the soliton profiles at α = 0, no analytical solutions for them is known so far.
It is not difficult to estimate the numbers of light quanta |F j | 2 /h and |S l | 2 /h in the FH and SH modes for the data of Figure 8c and Figure 8f. Using Equation (4) for µ and setting σ eff = 100 µm 2 , we have found that they are roughly of the order of 10 6 for |j| ∼ |l| ∼ 1. Thus, we are deeply within the classical range. This is valid also for the subsequent illustrations. At the same time, far tails of the comb spectra can experience quantum effects.
Next, we consider the effects of a nonzero walk-off. Nonzero values of α influence neither the thresholds nor the outcome of the adiabatic procedure-the system still evolves to a single-soliton steady state. The most evident effects are (i) a nonzero soliton velocity, v 01 = 0, and (ii) asymmetry of the soliton profiles and the comb spectra. The larger the |α|, the stronger is the asymmetry. This is illustrated by Figure 9. Compared to Figure 8, we changed δλ from 0 to −1.08 nm, which corresponds to α 0.05. The asymmetry of the soliton profiles and comb spectra is already evident. Furthermore, the soliton profiles acquire pronounced tails. The value of |F(ϕ)| 2 min is now far from zero, and |S(ϕ)| 2 max is about two times smaller than earlier. At the same time, the background values of the FH and SH intensities are the same, and the total number of the comb lines N = N 1 + N 2 is even slightly larger than it was at α = 0. Changing the sign of δλ results in the mirror reflection of the graphs of Figure 9 about the central vertical line and in the changed sign of the soliton velocity v 01 . To gain new insights, we consider color maps of the velocity parameter |v 01 |/2πR and of the total comb line number N = N 1 + N 2 on the δλ, η 2 plane for a few representative combinations of the decay rates γ 1 and γ 2 , see Figure 10. The first three subfigures of the upper row represent the maps of |v 01 |/2πR for the combinations γ 1,2 = γ (a1), γ 1 = γ, γ 2 = 10γ (b1), and γ 1,2 = 10γ (c1) with γ = 10 7 s −1 . Figure 10(a2,b2,c2) of the lower row represent the corresponding maps of the total comb line number N. Figure 10(d1,d2) give the additional cross-sections v 01 (δλ) and N(δλ) at η 2 = 20 for our combinations of γ 1 and γ 2 .
As evident from the maps, both |v 01 | and N grow monotonously with an increasing η 2 . Regarding the dependences of these parameters on |δλ| (or |α|), they are decreasing only outside a vertical central strip, whose width depends on γ 1,2 . The decrease of N(|δλ|) and |v 01 |(|δλ|) persists for |δλ| > 10 nm. In short, large values of the walk-off parameter, |α| 1, are harmful for the soliton-comb generation in question. The range of small |δλ| (or α), illustrated in some detail in the Figure 10(d1,d2), is worthy of attention. The dependences v 01 (δλ) and N(δλ) at large values of η 2 are not always continuous. This means that small variations in δλ can lead to notably different results when adiabatically increasing η 2 . Continuous changes in v 01 (δλ) and N(δλ) at η 2 = 20 in Figure 10 take place only for γ 1 = 0.1γ 2 = γ. The optimum value of |δλ| and the corresponding largest peak values of N take place for γ 1,2 = 10γ. In accordance with our definition of η 2 , they correspond to the largest h 2 (and the pump power). Anyhow, the increase in N via the optimization of δλ is not very large. Note lastly that the wavelength range of 20 nm relevant to Figure 10 strongly exceeds the intermodal wavelength distance of 0.08 nm so that the indicated spectral features are not available for the standard experimental fine-tuning means. Now we turn to the effects of the frequency detunings. These effects are of two kinds. One can consider the effects of δ 1,2 on the initiation of soliton-comb states via an adiabatic increase in η 2 . These strong and important effects are considered in Section 5.2. Alternatively, we can slowly vary δ 1,2 upon reaching the above-considered soliton states at large values of η 2 . This issue, which is largely (but not only) about the area of existence of the developed soliton-comb states, is considered below.
We have extended our adiabatic procedure for the slow independent changes in δ 1 and δ 2 . In fact, two similar scanning procedures were used. Within the first one, we initially determined the soliton velocity v 01 and the comb line number N by increasing and decreasing δ 1 at δ 2 = 0. After that, starting from point δ 1 , 0, we increased and decreased the detuning δ 2 . The achieved steady-state values of the amplitudes F j and S l were used as new initial conditions in each step of changing δ 1 or δ 2 . Within the second procedure, the scanning order was inverted: initially, we varied δ 2 at δ 1 = 0 and, after that, we increased and decreased δ 1 starting from points 0, δ 2 . We have verified whether the functions v 01 (δ 1 , δ 2 ) and N(δ 1 , δ 2 ) depend on the scanning order. Additionally, we have made sure that each steady state on the δ 1 , δ 2 plane is still relevant to a single-soliton antiperiodic state.
Consider first the case of zero walk-off, α = 0. In our coordinate frame moving with velocity v 1 , the positive and negative propagation directions are equivalent, such that the soliton velocity difference v 01 = v 0 − v 1 is expected to be zero. Surprisingly, we have found that this is not always the case. Figure 11a represents color map of |v 01 |/2πR for δλ = 0, γ 1 = 10 7 , and γ 2 = 10 8 s −1 (as in Figure 8). The uniformly dark-blue-colored lower and upper parts, which are separated from the rest by two slanted borders (bifurcation lines), correspond to v 01 /2πR = 0 within an accuracy of ∼ 10 −9 MHz. Between these borders we have a smooth distribution |v 01 |(δ 1 , δ 2 ) whose scale is comparable with that of Figure 10. This distribution practically does not depend on the scanning order, so what is the sign of v 01 ? We have found that it randomly depends on the fine particularities of our calculation procedure, such as the presence of very weak remnant velocity perturbations. In essence, we have a clear example of the spontaneous symmetry breaking when the state with a high symmetry (v 01 = 0) becomes unstable against the excitation of one of two equivalent states of lower symmetry with |v 01 | = 0. One of the most known examples of such a symmetry breaking is the ferroelectric second-order transition below the Curie temperature. Let us now take a look at the corresponding color map of N, Figure 11b. The distribution N(δ 1 , δ 2 ) is rather smooth. The lower bifurcation line of Figure 11a is barely seen here, while the upper line is pronounced. Remarkably, the increase in δ 2 at δ 1 ≈ 0 is favorable for the comb in spite of the increasing instability threshold, see Equation (10). The characteristic scale in δ 1 here is about one, which strongly exceeds the scale β 1 = 0.02 relevant to the changes in η th 2 (δ 1 ). This assertion is general for η 2 1; it is justified by plotting maps of N(δ 1 , δ 2 ) for different values of δλ and different combinations of γ 1 and γ 2 . In short, we have a vast family of single-soliton antiperiodic states continuously depending on the input parameters.
Our attempts to generate periodic soliton-comb solutions at δ 1,2 = 0 failed. The reason why and the possibilities to obtain new solutions are indicated in Section 5.2.

Selective Generation of Periodic and Antiperiodic Multisoliton States
The idea for the generation of periodic and antiperiodic multisoliton states is rooted in the threshold conditions of Section 3 and a simple observation. The point is that the FH profile of Figure 8 with a sharp minimum emerges with increasing η 2 as a continuous (without bifurcations) development of the near-threshold profile |F| 2 ∝ sin 2 (ϕ/2) corresponding to the lowest threshold for FH harmonics with j = ±1/2. One can thus suggest that the spatial structure of developed solitons can be determined by the threshold conditions.
To develop this idea, we consider the excitation diagrams of Figure 12. Diagrams (a) and (b) correspond to the antiperiodic and periodic cases. The cyan circles on the horizontal δ 1 /β 1 axis separate the regions where the lowest thresholds correspond to different values of the FH harmonic number |j|. Within a broad range in diagram (a), δ 1 > −β 1 /2, including the zero point, an adiabatic increase in η 2 is expected to lead to single-soliton antiperiodic states. This is in line with the above-considered results. For −9β 1 /4 < δ 1 < −β 1 /4, we can expect three-soliton antiperiodic states, etc.
The periodic case, diagram (b), is expected to be different. In the range δ 1 > −β 1 /2, including the zero point, the lowest threshold corresponds to j = 0, i.e., to the generation of the dual background. As this background is stable, no periodic solitons are expected. However, in the range −5β 1 /2 < δ 1 < −β 1 /2, where the lowest threshold corresponds to |j| = 1, a two-soliton periodic state is expected. Similarly, we can hope to generate four-soliton periodic states for δ 1 < −5β 1 /2. As concerned the detuning of δ 2 , it plays a passive role by nonselectively increasing all thresholds according to Equation (10).
For η 2 10, stable states grown via the adiabatic two-step (i-ii) procedure exist in a broad range of δ 1 /β 1 exceeding the shown one. Now we consider the results of the corresponding numerical experiments on the selective growth of multisoliton states [60]. Figure 13 shows the periodic two-soliton state grown via the adiabatic increase in η 2 for δλ = 0, γ 1,2 = γ = 10 7 s −1 , δ 2 = 0, and δ 1 = −β 1 . Remarkably, the phase arg[F](ϕ) shows two opposite π-steps when crossing each of the two narrow soliton areas. Thus, our periodic two-soliton solution consists of two equally spaced and mutually repelling antiperiodic solitons. This structure emerges every time spontaneously without imposing the antiperiodic conditions. The velocity of the state is v 0 = v 1 = v 2 . The FH and SH comb spectra, shown in Figure 13c and Figure 13f, are symmetric and structurally similar to those of Figure 8c and Figure 8f. Because of the π-periodicity in ϕ, the comb line spacing is 2v 1 /R, and it is doubled compared to the single-soliton case. Correspondingly, the number of significant comb lines (above 10 −4 ) is about two times smaller than it would be in the single-soliton antiperiodic case at η 2 = 20. Remarkably, the spatial structures analogous (or similar) to those presented in Figure 13a,b,d,e were obtained independently in [61][62][63] for the periodic states and were interpreted as opposite domain walls.
Changing δ 1 sequentially to δ (j) 1 = −4β 1 and −9β 1 , which corresponds to the lowest threshold for |j| = 2 and 3, we have grown periodic four-and six-soliton states. They all consist of pairs of the same antiperiodic solitons. The comb line spacing increases with |j| accordingly. A further increase in |j| does not lead, however, to new multisoliton states. Temporal development tends to break such states into the former two-, four-, and six-soliton solutions. This feature is not related to an insufficiently small time step and/or an insufficiently large number of harmonics taken into account. Additionally, we mention the case δ 1 = 0 relevant to the lowest threshold for j = 0. Our adiabatic procedure leads here, as expected, to a dual background steady stateF,S.
Generalization of these results to the case of nonzero walk-off (α = 0) and different decay rates γ 1 and γ 2 is rather trivial. We still have multisoliton periodic solutions consisting of pairs of the former antiperiodic solitons moving with velocity v 0 = v 1,2 .  Figure 13. Each of the three equally spaced solitons has a familiar structure with oscillating tails. Owing to a relatively large negative value of α, the tails are strongly pronounced and the spatial asymmetry is inverted as compared to Figure 9. The soliton velocity difference is v 01 2.2 × 10 5 cm/s. Interestingly, antiperiodic five-and seven-soliton states can be easily grown for |j| = 5/2 and 7/2 in spite of the strong tail overlaps. This contrasts with the multisoliton periodic states.

Natural PM: The Absence of Solitons
The natural phase matching that can be realized in LN-based resonators at λ 2 532 nm is of special interest. Avoiding the complicated procedure of the radial poling and the availability of tunable laser sources in this range make it very attractive. Furthermore, it was indicated [64] that the value of the walk-off parameter can be substantially decreased for sufficiently small microresonators (R < 1 mm) by using the effects of geometric dispersion. This implies the use of proper Sellmeier equations for n o,e (λ) [32,65] and relations for the modal frequencies including the size corrections [29], see Section S2 of the Supplementary Materials for details. It turns out, however, that the main problem with the use of natural PM for the soliton-comb generation is the dispersion issue [66]. In contrast to the aboveconsidered cases, the FH and SH dispersion parameters are both negative: v 1 −0.64 and v 2 −1.5 of 10 4 cm 2 /s.
Aiming to see the role of the signs of the dispersion parameters, we first conducted some test numerical experiments. For γ 1,2 = 10 7 s −1 , α = 0, δ 1,2 = 0, the absolute values |β 1 | = 0.02 and |β 2 | = 0.01, and four combinations of the signs of β 1,2 , we have launched our adiabatic procedure to grow the antiperiodic soliton states. In two cases with β 1 β 2 < 0, we obtained the same known single-soliton solution, see Figure 8. In two other cases with β 1 β 2 > 0, we obtained no soliton states, and the comb spectra remained undeveloped. Additionally, we varied the values of the input parameters α, γ 1,2 , and δ 1,2 and switched from the antiperiodic to periodic case. No soliton solutions were found as well. The sign of the product β 1 β 2 is thus crucial for the soliton formation.

Solutions for FH Pumping
The case of FH pumping is important for experiment because it corresponds to the conditions of SH generation, including the tuning issues. At the first sight, switching to the FH pumping case is trivial-it is sufficient to set h 2 = 0 and h 1 = 0 in sets (3) and (8). The normalized pump amplitude is here η 1 = √ 2µh 1 /γ 1 √ γ 1 γ 2 , and only the periodic states are allowed. While the threshold conditions are more complicated now, they can be evaluated numerically for any set of the input parameters. With these preliminaries, we consider the most important cases [66].
We start from the vicinity of the zero walk-off point 2λ c 2698 nm. The value of the walk-off parameter α is controlled here by the wavelength difference δλ = λ p − 2λ c , and the dispersion coefficients β 1,2 can be treated as constants. Setting δλ = 0 (α = 0) and δ 1,2 = 0 brings us to the situation where the lowest instability threshold, η th 1 = 3 √ 2 4.24, corresponds to j = 0, i.e., to auto-oscillations of the dual background [54]. Figure 15 illustrates what happens with increasing η 1 when applying our adiabatic procedure at γ 1,2 = γ = 10 7 s −1 . We do see auto-oscillations of |F 0 | 2 and |S 0 | 2 near the threshold, see Figure 15a. The oscillation amplitude and period increases and decreases with increasing η 1 , respectively, while the harmonics with j, l = 0 remain very small. For η 1 10, these spatial harmonics become notable. Their amplitudes grow gradually, and the number of significant comb lines reaches N ≈ 10 for η 1 = 11 with the comb line distance δj = δl = 1. The corresponding normalized FH and SH intensity profiles are shown in Figure 15b by lines 1 and 2. A further increase in η 1 results, however, in a sharp transition to an irregular spatial behavior, see Figure 15c and Figure 15d. This contrasts strongly with the case of SH pumping. Variations in the input parameters γ 1,2 and δ 1,2 do not change this unfavorable scenario. With the increase in |α|, the threshold values η th 1,j with j = 0 quickly decrease and become lower than the internal instability threshold. This is illustrated by Figure 6a relevant to our input parameters. For α 1, the modes with |j| 1 possess the lowest thresholds. The application of our adiabatic procedure to the case |α| = 1, δ 1,2 = 0 results, e.g., in the excitation of a π/3-periodic dual state (with six intensity peaks) near the threshold. With the increase in η 1 , it transforms via a bifurcation to a π/6 periodic state (at η 1 6.7). For η 1 16, the system shows again an irregular spatial behavior. The number of significant comb lines N does not exceed 25 for η 1 16. The employment of the nonzero detunings δ 1,2 does not help much.
For |α| 10, the threshold values η th 1,j (|α|) approach two from above, and the selective excitation of harmonics with |j| = |l| = 1, 2 becomes possible. However, the thresholds for |j| > 2 are only slightly higher. Consider representatively the case α = 11, δ 1 = 1/2, and δ 2 = 0 relevant to the lowest threshold for |j, l| = 1. Applying the adiabatic procedure, we have a π-periodic dual state near the threshold. For η 1 ≈ 7 and 18, it transforms sequentially via bifurcations to π/2-and π/4-periodic states. Correspondingly, the comb line distance δj = δl increases from 2 to 4 and lastly to 8. Figure 16a and Figure 16b show the complicated steady-state π/4-periodic asymmetric FH,SH intensity profiles for η max 1 = 20, while Figure 16c and Figure 16d exhibit the corresponding comb spectra. While the steady state found shows a regular spatial behavior, it cannot be regarded either as a soliton state or a pronounced comb state. For |α| = 50, the above scenario of the decreasing spatial period and increasing comb line distance with increasing η 1 remains the same. At η 1 = 20, we have a π/3 periodic dual state with a modest comb line number N = 20. The suppression of SH harmonics with l = 0 and the asymmetry of the intensity profiles and comb spectra are strongly pronounced.
Lastly, we touch on the case of the natural PM at λ p 1664 nm, when small and/or modest values of |α| are available owing to the geometric dispersion, see Section S2 of the Supplementary Materials for details. The values of the quality factors Q 1,2 are expected to be substantially larger than they are at 2698 nm [67], which leads to larger absolute values of the dispersion coefficients β 1,2 . As representative values, we chose R = 1 mm, β 1 = −0.04, and β 2 = −0.15, which corresponds to Q 1,2 ≈ 3 × 10 8 . The relevant dependences η th 1,j (α) are presented in Figure 6b. While they differ from the threshold dependences of Figure 6a, one important feature is the same: it is impossible to selectively excite the modes with |j| = 1, 2 for |α| 10.
A large number of our numerical experiments, conducted with different values of α and modest absolute values of δ 1,2 , persistently show the same scenario: increasing η 1 above the threshold leads (via bifurcations) to dual states with smaller and smaller spatial periods (larger and larger comb line spacing) and, ultimately, to an irregular spatial behavior. The larger α, the longer is the interval of η 1 with the regular behavior. On the other hand, the number of significant comb lines N within this interval remains modest; it does not exceed 15. No signs of soliton-comb states were detected. Possibly, such states were seen in numerical calculations of [68,69] for sufficiently large mismatches.

Discussion
Experimental and theoretical studies of dissipative χ (2) soliton-comb states in microresonators represent a vast, perspective, and almost unexplored research area. It is strongly different from the corresponding Kerr soliton research area. In turn, it is much broader than the studies of conservative Kerr solitons. The χ (2) nonlinearity is fundamentally different from Kerr's one, which results in a number of specific qualitative features including the duality of the nonlinear states and the necessity to employ the phase matching and minimize the temporal walk-off. The number of important variable parameters of the system is significantly larger than it is in the Kerr case. The deficit of analytical tools forces researches to rely on numerical methods.
Combining analytical and numerical tools and approaches, we predicted an unprecedentedly vast family of dual χ (2) soliton-comb steady states. When using the term "vast", we mean a continuous dependence of the parameters of the state on all the input parameters-the pump power, the modal quality factors, the frequency detunings, and the walk-off. The found states are self-starting, i.e., they are easily accessible via not very abrupt switching the pump on and changing the detunings. The soliton-comb states are divided into two topological classes-periodic and antiperiodic-depending on the parity of the azimuth number of the pumped SH mode. The steady state solutions can be singleand multisoliton ones. All our soliton-comb states are dissipative; they cannot be reduced to (obtained from) the known conservative soliton solutions.
Remarkably, the predicted self-starting soliton-comb states are allowed only for the SH pumping. Furthermore, they require opposite signs of the FH and SH dispersion coefficients. This condition is fulfilled near the zero walk-off point λ 2 1349 nm in LNbased resonators. However, the phase matching at this point requires a precise radial poling. Being available [15][16][17][18], it represents a serious technical problem. For natural phase matching in LN resonators, the necessary condition is not fulfilled, and soliton-comb states seem to be forbidden. Despite our strong efforts, no soliton-comb solutions were found for the FH pumping. This contradicts to the initial expectations about a similarity of the FH and SH pumping cases but has profound reasons behind it.
Both intuitive considerations and numerical modeling indicate that substantial differences in the FH-SH group velocities v 12 = v 1 − v 2 , which are typical of χ (2) resonators, are unfavorable for soliton-comb generation. Quantitatively, this means that the combination (|v 12 |/c) × (λ 2 /R) × Q 2 has to be not much larger than one. This condition is not generally fulfilled for nonstructured high-Q resonators. The minimization of the walk-off parameter can be accomplished via a proper radial poling. Thus, there is one more physical reason to employ the radial poling for the χ (2) comb generation.
Despite our massive failed efforts to find self-starting soliton-comb steady states for the FH pumping, we cannot claim that all sensible combinations of the input parameters (η 1 , α, β 1,2 , δ 1,2 , and γ 1,2 ) are exhausted. Additionally, we cannot claim that the initial assumption of the same transverse structure of the modes within each of the FH and SH frequency domains is fully justified. This problem still remains open for both the χ (2) and χ (3) resonators. It is likely that the suppression of the excitation of unwanted transverse modes requires further efforts for the mode management.
A distinctive feature of our approach is the employment of a simple classical equivalent of the Hamiltonian quantum formalism used in [4,5,17] and many other papers. While the quantum electrodynamic aspects of microresonators are important [17], the comb regimes belong to the classical range where the use of quantum language is redundant. As we have shown, the actual numbers of modal light quanta in the comb spectra are very large.
The effects of χ (3) nonlinearity in χ (2) microresonators is one more issue to consider. The possibility to neglect the χ (3) nonlinearity is largely controlled by the modal Q factors.
Let us consider for simplicity the SH pumping case. The pump power thresholds P (2) th and P (3) th for the parametric instability caused by the χ (2) and χ (3) nonlinearities are roughly proportional to 1/Q 3 and 1/Q 2 , such that for Q 1,2 10 6 , the effects of Kerr nonlinearity are expected to be negligible. For modest quality factors, Q 10 5 , which are typical, e.g., for ring waveguides, the situation can be more complicated. An interplay between the χ (2) and χ (3) nonlinear processes becomes possible [48].
The above-considered soliton-comb steady states in χ (2) microresonators are not the first predicted ones and are known in the literature. To the best of our knowledge, the other known stable dissipative soliton states [46][47][48][49] belong to the periodic ones. Furthermore, they are not self-starting, i.e., not easily accessible. The single-soliton numerical solutions of [46] are relevant to the FH pumping at the zero walk-off point. Both bright-bright and dark-dark locally stable solitons were exhibited. Similar numerical soliton solutions were found in [47,48]. Additionally, quasi-solitons at very large detunings, |δ 1,2 | 1, were predicted. Remarkably, the soliton-comb states of [46] are allowed when the dispersion coefficients have the same and opposite signs. The generation of flat-top solitons via an additional spatial modulation of the pump amplitude was found numerically in the SH pumping case when the dispersion coefficients had opposite signs [49]. A considerable amount of soliton-related models, see, e.g., Ref. [70], were proposed for the bulk parametric oscillators. Their detailed analysis is beyond the scope of this mini-review.
The first experimental attempts to obtain χ (2) combs in high-Q LN-based microresonators dealt with the FH pumping, the natural phase matching, and large walk-off [42,43]. Neither χ (2) solitons nor broad comb spectra similar to the Kerr soliton spectra were observed. This is in agreement with our theoretical results on the necessary conditions for the realization of such regimes. The authors of [44] claim for the discovery the Pockels soliton microcomb in an aluminium nitride microring resonator by employing SH pumping. The parameters used correspond to modest quality factors Q 1,2 < 10 6 , a large walk-off parameter |α| 1, and the threshold power P (2) th only slightly below P th . The observed broad FH comb spectrum was practically symmetric, which is the fingerprint of the absence of the dominating contribution of the χ (2) nonlinearity. The observed and numerically modeled combs are due to the Kerr nonlinearity. The quadratic nonlinearity is mostly involved in initiating the parametric instability, see also [71]. Thus, net χ (2) soliton-comb states are still waiting for experimental discovery and exploration.

Conclusions
A vast family of stable soliton-comb steady states is accessible in χ (2) microresonators by switching the monochromatic pump on not very abruptly. The necessary conditions for soliton-comb generation are (i) pumping into a SH resonator mode and (ii) the opposite signs of the FH and SH dispersion coefficients. Additionally, the minimization of the group velocity (free spectral range) difference for the phase-matched FH and SH modes is highly desirable. For LN-based resonators, the latter can be accomplished simultaneously with (ii) by using proper radial poling of the resonator. Pumping into FH resonator modes is not favorable for the realization of soliton-comb states. The employment of high-Q resonators (Q > 10 7 ) ensures low generation thresholds (within the µW range) and a negligibly small influence of the Kerr nonlinearity.