## Abstract

Single solitons are a special limit of more general waveforms commonly referred to as cnoidal waves or Turing rolls. We theoretically and computationally investigate the stability and accessibility of cnoidal waves in microresonators. We show that they are robust and, in contrast to single solitons, can be easily and deterministically accessed in most cases. Their bandwidth can be comparable to single solitons, in which limit they are effectively a periodic train of solitons and correspond to a frequency comb with increased power. We comprehensively explore the three-dimensional parameter space that consists of detuning, pump amplitude, and mode circumference in order to determine where stable solutions exist. To carry out this task, we use a unique set of computational tools based on dynamical system theory that allow us to rapidly and accurately determine the stable region for each cnoidal wave periodicity and to find the instability mechanisms and their time scales. Finally, we focus on the soliton limit, and we show that the stable region for single solitons almost completely overlaps the stable region for both continuous waves and several higher-periodicity cnoidal waves that are effectively multiple soliton trains. This result explains in part why it is difficult to access single solitons deterministically.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

It has long been known that single solitons are a special limit of a more general family of stationary nonlinear waves that are almost universally referred to as cnoidal waves in the plasma and fluid physics communities [1–6]. When loss and gain can be neglected, many fluid, plasma, and optical systems can be described at the lowest order by the nonlinear Schrödinger equation (NLSE). Examples in optics include light propagation in optical fibers [7], passively modelocked lasers [8], and microresonators [9–11]. In the NLSE approximation, solitons can be analytically expressed in terms of hyperbolic-secant functions [$\mathrm{sech}(x)$], while more generally cnoidal waves can be expressed in terms of Jacobi elliptic functions [$\mathrm{cn}(x)$, $\mathrm{sn}(x)$, $\mathrm{dn}(x)$] [12,13]. When loss and gain are present—as is the case in all experimental fluid and optical systems—then exact analytical solutions to the equations that describe these systems almost never exist. Nonetheless, it is common to refer to pulse and periodic solutions as “solitons” and “cnoidal waves,” respectively. In the past 15 years, it has become increasingly common to refer to real-world solitons as “dissipative solitons” [14], and by analogy it seems reasonable to refer to real-world stationary periodic solutions as “dissipative cnoidal waves”.

The literature on cnoidal waves in fluids is extensive, but that is not the case in optics, except in the soliton limit. (See, however, [15–17] for theoretical discussions.) Solitons or cnoidal waves that are described at the lowest order by the NLSE are created by a balance between nonlinearity and dispersion. The pulses in optical fiber propagation experiments are typically separated by many times their duration in order to avoid instabilities [18]. Even in experiments in optical fiber ring resonators [19] or passively modelocked lasers [20] with many pulses in the cavity, it is useful to model the pulses as a periodic train of weakly interacting solitons, rather than as a special limit of cnoidal waves. However, the small dimensions of microresonators make it easier to observe more general cnoidal waves than is possible in other optical systems, and they have frequently been observed [21–28], although often remarked upon only as a stepping stone on the path to single solitons [21–24].

In recent years, a large experimental effort has focused on obtaining broadband microresonator combs using single solitons, and this effort has achieved remarkable experimental successes, including the demonstration of octave-spanning combs [29,30], a $2f-3f$ self-referenced comb [31], and an optical frequency synthesizer [32]. However, these combs cannot be obtained in a straightforward way. In contrast to many modelocked lasers, the microresonators do not “self-start,” i.e., when the microresonator is turned on, single solitons do not simply appear. Approaches to generate single solitons include pump tuning, thermal control, and engineered spatial mode interactions [33–36]. However, there is often a random element in the generation of single solitons [35,36], and the search for robust, deterministic paths to obtain broadband combs continues. Recent work indicates that is possible to generate single solitons deterministically by using trigger pulses [37], by using two pumps [38–41], or by adding an additional continuous wave component [42]. However, these approaches make the system significantly more complex. Another drawback to the use of single bright solitons to generate combs is that they make inefficient use of the pump. Typically, less than 1% of the pump power is transferred to the comb [24].

In this paper, we will show that more general cnoidal waves can potentially solve these problems. Cnoidal waves can be broadband, while at the same time they can be straightforwardly accessed and use the pump more efficiently. They are sufficiently easy to obtain that they have frequently been observed, although, as previously noted, they have often been passed over. Their robustness has, however, been previously noted [25,26], and narrowband combs have potential applications to optical communication systems [26].

An essential difficulty in studying the potential of cnoidal waves is that there is a large parameter space to explore. As previously noted, single solitons correspond to period-1 cnoidal waves, in which the azimuthal repetition period of the light in the microresonator cavity is equal to the entire azimuthal period, so that there is only a single pulse in the cavity. By contrast, cnoidal waves can have any periodicity, and a period-${N}_{\mathrm{per}}$ cnoidal wave has an optical amplitude that repeats ${N}_{\mathrm{per}}$ times in one azimuthal period. Moreover, the cnoidal waves with a given periodicity ${N}_{\mathrm{per}}$ can vary greatly from waveforms in which the light intensity is never close to zero to waveforms that are effectively a periodic train of solitons or a soliton crystal. Indeed, this variability is their great advantage. We will show that it is possible to continuously tune the microresonator parameters to effectively obtain a train of solitons without ever passing through a chaotic regime, as is required to obtain single solitons.

Given the size of the parameter space to be explored, more theoretical guidance is needed, and direct solution of the evolution equations is inadequate for reasons that we will discuss shortly.

To address this issue, we have adapted a unique set of computational tools that were previously developed to study lasers [43–46] and that allow us to rapidly determine where in the parameter space cnoidal wave solutions exist and to find their frequency spectra. These tools are based on dynamical systems theory [47–49] and complement widely available theoretical tools that are based on direct solutions of the evolution equations. The basic idea is that once a cnoidal wave solution of a given periodicity has been found, one varies parameters, finding the cnoidal wave solutions as parameters vary by solving a nonlinear root-finding problem, instead of solving the evolution equations. In parallel, one determines the stability of the solution by finding the eigenvalues of the linearized equations. Once a stability boundary is encountered, where the cnoidal waves become unstable or cease to exist, we track the boundary [43]. This approach is computationally rapid, and it allows us to identify all the stable cnoidal wave solutions that exist for a given set of parameters, in contrast to standard evolutionary simulations that will typically only identify a single solution, depending randomly on the initial conditions that are used. This approach also allows us to find unstable as well as stable cnoidal waves. That is important because unstable cnoidal waves can persist for a long time, and this approach allows us to calculate that time scale. When the lifetime of an unstable cnoidal wave is long, standard evolutionary simulations can incorrectly conclude that it will be stable.

The basic ideas that we are using are old, dating back at least to Maxwell’s study of the stability of Saturn’s rings [47]. Soliton perturbation theory is based on these ideas, and they have been adapted by Haus [48,49] and others [50] to study passively modelocked lasers. They have been used by Matsko and Maleki [51] to study soliton timing jitter and by Godey *et al.* [52] to study the stability of continuous waves in microresonators. Barashenkov *et al.* [53,54] have applied these ideas to the study of multiple soliton solutions of the driven-damped NLSE [55], which in the microresonator community is typically referred to as the Lugiato–Lefever equation (LLE) [56]. However, almost all of this work is based on analytical approximations of the stationary solutions. (The work in Refs. [53,54] is an exception.) In our approach, we calculate computationally exact stationary solutions, which allows us to accurately study the solutions in any parameter regime. That is particularly important when studying cnoidal waves in the anomalous dispersion regime, which is the focus of this paper. While analytical solutions exist for the cnoidal wave solutions of the LLE when damping is neglected [57], these solutions are disconnected from the stable cnoidal wave solutions in the anomalous dispersion regime except when the periodicity is small, as discussed in Sections 4 and 1 of Supplement 1. As a consequence, the utility of a perturbation theory based on these analytical solutions appears likely to have limited value.

The microresonator community has developed a nomenclature that differs in some respects from the nomenclature in common use in the nonlinear waves community. The use of the LLE instead of the “driven-damped NLSE” is one example. The periodic solutions of the LLE have been referred to as “Turing rolls,” “primary combs,” and “soliton crystals,” as well as “cnoidal waves” or “dissipative cnoidal waves.” The choice “Turing roll” connects the cnoidal wave solutions with earlier work on pattern formation [58]; the choice “soliton crystal” emphasizes the particle-like nature of solitons and connects the periodic solutions to single soliton solutions. As noted previously, our choice of “cnoidal waves” connects these solutions to the large body of work on periodic solutions in fluids, plasmas, and nonlinear waves.

The remainder of this paper is organized as follows: in Section 2, we present our basic equations, discuss our normalizations and how they relate to the experimental literature, and review our dynamical approach. In Section 3, we present our results on the stability of the cnoidal wave solutions—showing where in the parameter space stable solutions exist. We discuss their accessibility and how to optimize their parameters, including their bandwidth. In Section 4, we discuss the single soliton solutions. Our approach sheds new light on the difficulties in accessing them. Finally, Section 5 contains the conclusions.

## 2. BASIC EQUATIONS, NORMALIZATIONS, AND DYNAMICAL METHODS

Light propagation in microresonators is described at the lowest order by the LLE [9–11], which may be written as

Equation (1) can be written in terms of three normalized parameters. One standard way to normalize Eq. (1) is to let $t=l\tau /2{T}_{R}$, $\psi ={(2\gamma /l)}^{1/2}A$, $\beta =2{\beta}_{2}/l$, $\alpha =2\sigma /l$, and $F=i{(8\gamma /{l}^{3})}^{1/2}\sqrt{{P}_{\mathrm{in}}}$, in which case Eq. (1) becomes

An alternative is to normalize Eq. (1) with respect to the detuning, assumed positive, so that $t=\sigma \tau /{T}_{R}$, $x={(2\sigma /|{\beta}_{2}|)}^{1/2}\theta $, $\psi ={(\gamma /2\sigma )}^{1/2}A$, $\delta =l/2\sigma $, and $h={(\gamma /2{\sigma}^{3})}^{1/2}\sqrt{{P}_{\mathrm{in}}}$. In this case, Eq. (1) becomes

In experimental work, typical values of $F$ range from 0 to 4, while typical values of $\alpha $ range from $-5$ to 10. The corresponding device parameters vary greatly, depending on the material from which the device is made and its size. Table 1 shows the experimental parameters that correspond to $\alpha =1$ and $F=1.5$ for several different experimental devices. The normalized mode circumference $L$ ranges between 30 and 200.

The standard evolutionary methods for solving the LLE are typically based on the split-step method and its variants [7]. We use a variant that we recently demonstrated is significantly faster computationally than most other variants when large loss and gain are present [60]. In a standard evolutionary simulation, one starts from an initial condition, and one steps the LLE forward in time, sometimes allowing the equation’s parameters to also change as a function of time. When a stationary solution is observed at the end of the simulation, one has learned where in the parameter space that particular stationary solution can be accessed, starting from the initial conditions and along the path in the parameter space that was used in the simulation. By contrast, dynamical methods for solving the LLE reveal where in the parameter space the stable solutions exist, but do not reveal how to access them. The dynamical methods function much like an x-ray machine that tells a surgeon where a bone is broken, but not how best to repair it. Hence, the evolutionary and dynamical methods are complementary.

While evolutionary methods are widely used in the microresonator community, dynamical methods are not. Our implementation is an adaptation of computational methods that were developed to study passively modelocked lasers [43–46]. Our starting point is a cnoidal wave solution that we find using evolutionary methods with equation parameters for which the solution is highly stable. This solution is stationary and satisfies Eq. (3) with its time derivative set equal to zero,

We now solve Eq. (5), while allowing $F$ or $\alpha $ to vary, which becomes a nonlinear root-finding problem. We use the Levenberg–Marquart algorithm [61], which is a variant of the well-known Newton’s method.In parallel with solving Eq. (5) to find the cnoidal wave solutions, we determine their stability. We write

*et al.*[43] and as shown schematically in Fig. 1. As we track the boundary, we find that transitions occur among these three cases.

The ultimate fate of the system when the parameters pass through a stability boundary must be determined by solving the evolution equations using evolutionary methods. We have found that case (1) leads to a transition from a cnoidal wave solution to a continuous wave solution, case (2) leads to a transition to another cnoidal wave solution, and case (3) leads to a transition to a breather or chaos. These three cases are referred to in the mathematics literature as saddle-node, transcritical, and Hopf bifurcations, respectively [62,63]. In the mathematics literature, the set of eigenvalues $\{{\lambda}_{j}\}$ are referred to as the “spectrum” of the operator $\mathcal{L}$ [64,65]. Here, we refer to the $\{{\lambda}_{j}\}$ as the “dynamical spectrum” in order to avoid confusion with the frequency spectrum. The dynamical spectrum of the LLE has a large amount of symmetry. It is symmetric about the real $\lambda $ axis, since if ${\lambda}_{j}$ is an eigenvalue, then so is ${\lambda}_{j}^{*}$. Additionally, it is symmetric about the axis $\mathrm{Re}(\lambda )=-1$. We discuss some features of the dynamical spectrum in Section 3 and in Section 2 of Supplement 1. While interesting, the details of this structure are of little importance for this study, since almost all of the eigenvalues correspond to perturbations that rapidly damp out. Our focus is almost exclusively on the eigenvalues whose real parts become equal to zero as the microresonator parameters change.

Once the stable region of a cnoidal wave with a given periodicity has been found, its properties can be rapidly found throughout the region by solving Eq. (5) as $\alpha $ and $F$ vary. We have used this technique to find the fraction of the pump power that goes into the cnoidal wave frequency comb, as well as the comb bandwidth. By solving the evolution equations, we have verified that any solution within the stability boundary can be experimentally accessed in principle by changing the pump power and the detuning.

The dynamical methods that we are using have significant advantages over just using evolutionary methods to identify stable operating regions in the parameter space and paths through the parameter space to access them. First, directly solving the evolution equations for given initial conditions and system parameters can miss stable solutions that exist for those parameters. As shown schematically in Fig. 2, a stable waveform has a basin of attraction in the phase space of all possible waveforms. If the initial conditions are within this basin, then the evolution will eventually converge to this waveform, but, if it does not, then the evolution will not converge to this solution, and it is possible to miss it completely. This issue is particularly important for the cnoidal wave solutions, where several different periodicities can be stable at the same point in the parameter space. Second, evolutionary methods yield ambiguous results near a stability boundary, since the time for a perturbation to either decay or grow tends toward infinity. Examples of this ambiguity appear in Section 3 and in Section 1 of Supplement 1. Third, the dynamical methods are computationally rapid, making it possible to determine the stability boundaries for all the cnoidal wave solutions within a broad parameter range. We have typically found that this approach is 3–5 orders of magnitude faster than evolutionary methods, depending on how well one wants to resolve the stability boundary.

These same issues arise in any driven-damped system, in which stationary solutions such as single solitons, molecules, and periodic trains of solitons appear. Akhmediev *et al.* [66], for example, have discussed these issues in the context of the complex Ginzburg–Landau equation. This equation approximately describes a variety of optical and non-optical systems [14].

In most computational studies of the dynamics in microresonators, the instabilities are seeded by random numerical noise due to roundoff. This seeding implies a random noise amplitude on the order of ${10}^{-15}$ in modern-day 64 bit computers. This noise amplitude is many orders of magnitude lower than what is expected in an experimental system, which is largely set by environmental or “technical” noise [67,68]. By contrast, the randomly varying amplitude due to quantum noise is in the range of ${10}^{-5}\u2013{10}^{-4}$ for typical parameters, as discussed in Section 3 of Supplement 1. The growth of the amplitude of an unstable eigenmode is given by $A(t)={A}_{0}\text{\hspace{0.17em}}\mathrm{exp}({\lambda}_{R}t)$, where ${\lambda}_{R}$ is the real part of the eigenvalue. As long as ${\lambda}_{R}^{-1}$ is small compared to the time scale of the simulation by factors larger than about $\mathrm{ln}\text{\hspace{0.17em}}{10}^{-15}\simeq -35$, then the small level of seeding from roundoff noise does not matter. When making a transition out of a region in the parameter range that is chaotic or where continuous waves are stable into a region where cnoidal waves are stable, the growth of the cnoidal waves is sufficiently rapid that it is not unreasonable to seed the growth from roundoff noise. We will see, however, that the time scale for instabilities to manifest themselves when making transitions from one cnoidal wave to another can be a sizable fraction of a second in real time, and these instabilities are computationally challenging to simulate. For that reason, we use a quantum noise seed, as described in Section 3 of Supplement 1. While the quantum noise power is lower than the technical noise power in microresonators by one to three orders of magnitude [68], it is large enough for instabilities to manifest themselves on physically realistic time scales.

## 3. STABILITY, ACCESSIBILITY, AND OPTIMIZATION OF CNOIDAL WAVES

In Fig. 3, we show the regions where the cnoidal waves are stable for different periodicities ${N}_{\mathrm{per}}$ in the range $-2<\alpha <6$ and for $L=10$, 25, and 50. For clarity, we only show some of the stable regions. However, we have found all the stable regions in the parameter range that we show, and a complete map of the stable regions for $L=50$ is given in Section 4 of Supplement 1. When $L\to \infty $, we find that $L/{N}_{\mathrm{per}}$ tends to a constant for a fixed value of $\alpha $ and $F$, as we show in Fig. 4. The value in the limit $L\to \infty $ is calculated using an approach like that of Godey *et al.* [52]. More details are in Section 5 of Supplement 1. As a consequence, the three-dimensional parameter space in which all solutions of the LLE exist effectively reduces to a two-dimensional parameter space when $L\ge 50$.

In Fig. 3, we see that the region where cnoidal waves are stable form a U-shaped belt (red-dashed curve) in the $\alpha -F$ plane. Below this belt, only continuous waves (flat solutions) are stable and are always obtained. Above this belt, there are no stationary solutions; only breathers or chaotic solutions are obtained, depending on the parameters. An analytical expression for this curve is derived in [52] and is given by ${F}^{2}=1+{(1-\alpha )}^{2}$. An important transition occurs at $\alpha =41/30$. Below this value of $\alpha $, stable continuous waves and cnoidal waves do not both exist at the same parameter values. Hence, cnoidal waves can be easily accessed by simply raising the pump power. Above this value of $\alpha $, stable cnoidal waves and continuous waves can both exist at the same parameter values. In particular, continuous waves are stable whenever single solitons are stable. Hence, single solitons cannot be accessed by simply raising the pump power. Additionally, the stable region for single solitons overlaps almost completely the stable region for several higher-periodicity cnoidal waves when $L\ge 25$. The overlap of the stable region for single solitons with other stationary solutions explains in part why they are difficult to access.

In Fig. 5, we show the optical waveform, its frequency spectrum, and its dynamical spectrum for the ${N}_{\mathrm{per}}=8$ solution with $L=50$, $\alpha =4.4$, and $F=2.3$, corresponding to the red dot in Fig. 3(c). Figure 5(a) shows the stationary waveform’s intensity, ${|\psi (x)|}^{2}$. Figure 5(b) shows ${P}_{n}$, defined as ${P}_{n}={|{\tilde{\psi}}_{n}|}^{2}/{|{\tilde{\psi}}_{0}|}^{2}$, where ${\tilde{\psi}}_{n}=(1/L){\int}_{-L/2}^{L/2}\psi (x)\mathrm{exp}(-i2\pi n{N}_{\mathrm{per}}x/L)\mathrm{d}x$. Each comb line is spaced eight times the FSR apart. This ${N}_{\mathrm{per}}=8$ cnoidal wave is effectively a set of eight periodically spaced solitons—a periodic train of solitons or a soliton crystal—that is stabilized by sitting on a small pedestal. The power ${P}_{n}$ of the comb lines falls off exponentially from the central line. This clean, exponential falloff is characteristic of all cnoidal waves and distinguishes them from a random assortment of solitons or other pulses. In Fig. 5(c), we show the dynamical spectrum of the waveform, and in Fig. 5(d), we show an expanded view of a small region near $\lambda =(0,0)$. The dynamical spectrum has a fourfold symmetry about the real axis and the axis $\mathrm{Re}(\lambda )=-1$. There is always an eigenvalue at $\lambda =(0,0)$, corresponding to azimuthal invariance of the solution, and, hence, there is an eigenvalue at $\lambda =(-2,0)$. Except for these two points, we have found that all of the eigenvalues on the real axis are degenerate. We see that a group of real eigenvalues cluster near zero. These eigenvalues correspond roughly to modes that shift the amplitudes and phases of the pulses or solitons that make up the cnoidal wave, while eigenvalues along the axis $\mathrm{Re}(\lambda )=-1$ correspond roughly to modes that perturb the pedestal. We explain some of the principal features of this spectrum in Section 2 of Supplement 1. However, most of these features are of little relevance for this study. Our focus is almost entirely on eigenvalues that reach the imaginary axis as $\alpha $ and $F$ change, implying that the stationary solution either becomes unstable or disappears.

In Figs. 5(e) and 5(f), we show the optical waveform and its frequency spectrum for $\alpha =1$ and $F=1.15$, corresponding to the blue dot in Fig. 3(c). In this case, the ${N}_{\mathrm{per}}=8$ cnoidal wave appears to be a modulated continuous wave rather than a periodic train of solitons. Starting with this $F$ and $\alpha $ and continuously varying the parameters to $\alpha =4.4$ and $F=2.3$ makes it possible to lock in place the periodic soliton train in Fig. 5(a).

In Fig. 6, we show examples of the transitions that occur when a stability boundary is crossed, starting in the region where the ${N}_{\mathrm{per}}=8$ solutions are stable. In Fig. 6(a), we show the stable regions for ${N}_{\mathrm{per}}=8$ and ${N}_{\mathrm{per}}=9$, as well as four trajectories in the parameter space. We show the transitions in the dynamical spectrum for all four trajectories, focusing on the eigenvalues that hit or cross the imaginary axis. For clarity, we have omitted other eigenvalues. We show the initial eigenvalues as red dots and the final eigenvalues as blue circles. We also show the stationary solution as a function of the variable that we are varying, $\alpha $ or $F$. In Figs. 6(b) and 6(c), we show the dynamical spectrum and the stationary pulse intensities for trajectory A, for which we start at $\alpha =2$, $F=1.29$ and decrease $F$ to 1.27. Slightly below $F=1.28$, a saddle-node bifurcation occurs, at which point the cnoidal wave solution and, hence, its dynamical spectrum cease to exist. In this case, the final point in the dynamical spectrum is at $\lambda =(0,0)$. The solution evolves from a cnoidal wave into a continuous wave. In Figs. 6(d) and 6(e), we show the dynamical spectrum and stationary pulse intensities for trajectory B, for which we start at $\alpha =4$, $F=2.3$, and we increase $F$. A Hopf bifurcation occurs; the cnoidal wave continues to exist, but it is unstable. It evolves into a breather solution, also referred to as a secondary comb. In Fig. 6(e), we show the region beyond the Hopf bifurcation in gray since no stationary solution exists. In Figs. 6(f) and 6(g), we show the dynamical spectrum and pulse intensities for trajectory C, for which we start at $\alpha =2$, $F=1.83$, and increase $F$. A Hopf bifurcation occurs, and the ${N}_{\mathrm{per}}=8$ cnoidal wave evolves into a stable ${N}_{\mathrm{per}}=9$ cnoidal wave. Finally, in Figs. 6(h) and 6(i), we show the dynamical spectrum and pulse intensities for trajectory D, for which we start at $\alpha =0.85$, $F=1.12$, and decrease $\alpha $. A transcritical bifurcation occurs, and the stable ${N}_{\mathrm{per}}=8$ cnoidal wave evolves into a stable ${N}_{\mathrm{per}}=9$ cnoidal wave. However, the time scale on which this evolution occurs is noticeably longer than is the case for trajectories A, B, and C. That is a common feature of transcritical bifurcations and underscores the importance of using a realistic background noise level and sufficiently long integration times when using evolutionary methods.

Comparing Figs. 6(d) and 6(f) to Fig. 6(h), we see that ${\lambda}_{R}$, the maximum real value that appears in the dynamical spectrum, is about 25 times smaller for a transcritical bifurcation than it is after a Hopf bifurcation for a comparable excursion in the parameter space. As a consequence, the time scale for the instability to manifest itself is considerably longer for a transcritical bifurcation. In Fig. 7, we show the evolution as a function of normalized time and for physical time, using the parameters of Wang *et al.* [24], starting from the unstable stationary cnoidal wave solutions at the final points of Figs. 6(d), 6(f), and 6(h). As expected the time scale in Fig. 7(c), corresponding to Fig. 6(h), is about 25 times larger than the time scale in Figs. 7(a) and 7(b), corresponding to Figs. 6(d) and 6(f). The dynamical spectrum is particularly useful when studying transcritical bifurcations since, without its guidance, one can easily integrate for too short a time to accurately assess the stability of the solution.

In Fig. 8, we compare the stability regions for $L=50$ with evolutionary simulations for $L\simeq 46$. We show a similar comparison in Section 1 of Supplement 1 for the normalization of Eq. (4). In Fig. 8(b), we show the number of peaks in the final solution after jumping from an initial cnoidal wave at $\alpha =0$ and $F=6.3$ to another point in the parameter space and remaining at that point for $\tau =1.8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu s}$ [70]. While the results correspond roughly to the stability curves, differences should be noted. The dwell time is not always sufficient to determine the stability, and the location of the stability boundaries is ambiguous. Additionally, the passage time through the parameter space affects the final state. We will discuss this issue in more detail in Section 4. By moving along the trajectory in the phase space, shown as the solid black curves in Fig. 8(a) or the red-dashed curve in Fig. 8(b), it is possible to stay within the region in the parameter space where the cnoidal waves are stable and deterministically control, with occasional back-tracking, the periodicity of the cnoidal wave that is asymptotically obtained.

In Fig. 9, we show an optimization study for the ${N}_{\mathrm{per}}=8$ cnoidal wave when $L=50$. The point at $\alpha =4.4$, $F=2.3$, shown as a red dot in Figs. 3 and 9, corresponds to the frequency spectrum that we show in Fig. 5(b). The point at $F=1.15$ and $\alpha =1$ corresponds to the blue dot. The frequency spectrum always decreases exponentially away from the central peak. In Fig. 9(a), we show a contour plot of the rate of exponential decrease, ${P}_{n}/{P}_{n+1}$ (dB) for $n\ge 3$. In Fig. 9(b), we show a contour plot for the ratio $\rho $ of the power in the comb lines to the input pump power, $\rho =(1/{F}^{2}){\int}_{-L/2}^{L/2}{|\psi |}^{2}(\mathrm{d}x/L)$. The corresponding change in the physical ratio of the power in the comb to the power in the pump is given by ${\rho}_{\mathrm{phys}}=(4/{l}^{2})\rho $, which for the parameters of [24] and a single soliton is given by ${\rho}_{\mathrm{phys}}=73.6\%$, which is close to the measured value at the drop port. There is a tradeoff between the efficiency with which the pump is used and the bandwidth of the frequency spectrum. At the point in the parameter space shown as a red dot, we find that $\gamma <6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{dB}$, while $\rho \simeq 0.3$. The comb lines in the frequency spectrum are separated by eight times the FSR. The lines at $n=\pm 6$ are down by 30 dB from the central line [See Fig. 5(b)]. For an FSR of 230 GHz, corresponding to the experiments in [24], the corresponding bandwidth is 21 THz, which is a large fraction of the operating frequency of 200 THz. This bandwidth is consistent with the bandwidth of single solitons in microresonators, whose bandwidth has not been enhanced through the use of dispersive radiation. It should be possible to apply techniques that have been used to enhance the bandwidth of single solitons [30,31] to cnoidal waves with larger periodicity. While the comb line spacing (1 THz) is large, it is consistent with the line spacing of the single soliton comb with large line spacing in the dual-microcomb experiment of [32].

To characterize the increase in efficiency of pump utilization as ${N}_{\mathrm{per}}$ increases, we computed $\rho $ for periodicities ${N}_{\mathrm{per}}=1-8$ with $\alpha =4.4$ and $F=2.3$. We found that $\rho $ increases almost linearly, so that $\rho =0.038$ when ${N}_{\mathrm{per}}=1$ and $\rho =0.287$ when ${N}_{\mathrm{per}}=8$. This linear increase is not surprising since at these values of $\alpha $ and $F$, the cnoidal wave is effectively a periodic train of ${N}_{\mathrm{per}}$ solitons. Since the total power is nearly proportional to ${N}_{\mathrm{per}}$, and the comb lines are spaced ${N}_{\mathrm{per}}\times \mathrm{FSR}$ apart, the power in each comb line is nearly proportional to ${N}_{\mathrm{per}}^{2}$. That can be important for applications where the power in a single comb line must be made large.

## 4. SOLITON LIMIT

Single solitons are a special limit of cnoidal waves for which ${N}_{\mathrm{per}}=1$. They share many of the properties of cnoidal waves for which ${N}_{\mathrm{per}}>1$. In particular, their frequency spectrum falls off exponentially away from the central peak. Single solitons are also a special limit of another class of nonlinear stationary waveforms, referred to as multi-bound solitons or soliton molecules [54,71]. Except in special limits where they are single solitons or higher-periodicity cnoidal waves, soliton molecules are characterized by a complicated frequency spectrum.

Single solitons in microresonators are not easy to access, and a large experimental effort has been aimed at finding ways to do that [21–27,29–36,43]. Here, we focus on understanding the issues with accessing them from a dynamical perspective, and we discuss the reasons that there is often a random element in their generation [35,36]. We extend an earlier discussion that proposed a deterministic path through the parameter space to obtain single solitons [70], and we describe the critical role that is played by the time scale on which one moves through that path.

In Fig. 10, we show the stable regions for the ${N}_{\mathrm{per}}=1$, 3, and 5 cnoidal waves when $L=50$. For clarity, we do not show the regions of stability for ${N}_{\mathrm{per}}=2$, 4, or 6. However, they may be seen in the map of the stable regions in Section 4 of Supplement 1. While the stable region for single solitons (${N}_{\mathrm{per}}=1$ cnoidal waves) is large, it almost entirely overlaps the regions of stability for several higher-periodicity cnoidal waves in the parameter range that we show, as well as with continuous waves. Hence, it is not surprising that multiple solitons are often obtained in experiments rather than single solitons. However, the multiple solitons are usually randomly spaced, leading to a complicated frequency spectrum, which is usually undesirable in applications. While these states are unlikely to be stationary states, simulations indicate that they can be long-lived, with lifetimes that far exceed standard simulation runs.

We can obtain insight into this behavior by examining how single solitons are generated. The most common path to obtain single solitons is to start in the blue-detuned region, where either continuous waves or high-periodicity cnoidal waves are stable, at the left side of the U-shaped region shown in Fig. 3 [21–26]. Solitons are then obtained by red detuning the pump laser through the chaotic region above the U-shaped region, where cnoidal waves are stable into the region where single solitons, other low-periodicity cnoidal waves that are effectively a periodic train of solitons, and continuous waves are stable. This path is labeled A in Fig. 11. An alternative path, labeled B in Fig. 11, is to raise the pump power above the line where continuous waves are stable and then to lower the power until single or multiple solitons are stable [23,29,34]. A third alternative, labeled C in Fig. 11, is to tune backward after entering the region where single solitons are stable [35]. Since the soliton duration is inversely proportional to the detuning, backward tuning enhances the interaction between randomly created multiple solitons and thus hastens their collapse into a single soliton, a continuous wave, or another stationary state. In all these paths, the region where single solitons are stable is entered from the chaotic region, in which both temporal and spatial power fluctuations are comparable to the average power and are several orders of magnitude larger than the underlying noise level in the device [67]. As a consequence, the spacing of the solitons when they form will be random, and both their number and spacing will vary from shot to shot. Their ultimate evolution will depend on the details of their formation. Since continuous waves and a number of higher-periodicity cnoidal waves are also stable, nothing in principle prevents their formation.

The stable region for single solitons that we show in Fig. 11 is bounded by curves at which either a saddle-node bifurcation occurs (cyan curves) or a Hopf bifurcation occurs (blue curve). In Fig. 12, we show the dynamical spectra for the stable soliton for $L=50$, $\alpha =5.5$, and $F=2.11$, 2.45, and 2.90, shown as the red dots in Fig. 11. We also show the corresponding waveforms and frequency spectra. It is evident that the dynamical spectra differ significantly from the dynamical spectrum that appears in soliton perturbation theory [72] or the Haus modelocking equation [48,49]. In that case, the only discrete eigenvalues correspond to amplitude, frequency, time, and phase shifts. There are also continuous eigenvalues that correspond to background radiation. In this case, near the value at which the Hopf bifurcation occurs ($F=2.9$), the least-damped eigenvalues correspond to damped-periodic oscillations of the single soliton. Near the value at which the saddle-node bifurcation occurs ($F=2.11$), the least-damped eigenvalue corresponds to a damped amplitude and phase modulation. The cluster of eigenvalues near $\lambda =(0,0)$ that was visible in the case of ${N}_{\mathrm{per}}=8$, shown in Fig. 5, and that correspond to damped interactions between the solitons that make up the cnoidal wave, are not present.

As a consequence, we might expect single solitons to be more robust in the presence of large noise or environmental fluctuations. We have verified, doing a long simulation run for which $t={10}^{8}$, that the ${N}_{\mathrm{per}}=8$ cnoidal wave illustrated in Fig. 5 remains stable. These simulations are computationally taxing, requiring many hours on a standard desktop computer. However, this normalized time $t$ only corresponds to about 1 s of physical time for a resonator quality factor $Q=1.67\times {10}^{6}$, corresponding to the parameters of [70]. Determining their nonlinear stability in the presence of noise fluctuations over the lifetime of an experiment remains to be done.

A promising approach to deterministically obtain low-periodicity, large-bandwidth cnoidal waves, including single solitons, is to avoid the chaotic region [70]. Rather than just changing the detuning or the power, it is advantageous to change both together so that the system moves along one of the paths, shown as the black curves in Fig. 8(a) or the red-dashed curve in Fig. 8(b). In this case, transcritical bifurcations occur along the trajectory through the parameter space when cnoidal waves become unstable. The time scale of these instabilities, which is typically milliseconds, is long compared to most simulations, although short compared to most experiments. It is possible to take advantage of this long time scale to deterministically create single solitons [70].

In Fig. 13, we show the evolution when the system parameters are shifted from $\alpha =0$ to $\alpha =6$ along the red-dashed curve in Fig. 8(b) at different speeds, and the system is then allowed to evolve to a stationary final state. We set $L=50$ and use $Q=2.7\times {10}^{6}$ to convert from normalized time to physical time, ${t}_{\mathrm{phys}}$. The expression for the red-dashed curve is

In Fig. 13(a), we show the evolution when the initial evolution from $\alpha =0$ to $\alpha =6$ occurs rapidly, analogous to the case considered in [70]. In this case, the parameters are initially fixed at $\alpha =0$, $F=6$ and remain there until $t=100$ ($\tau =0.46\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu s}$). By that time, a stable ${N}_{\mathrm{per}}=13$ cnoidal wave has emerged. The detuning is then increased linearly up to $\alpha =6$, arriving at $t=250$ ($\tau =1.15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu s}$). Finally, we wait until $t=600$ ($\tau =2.76\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu s}$). The original ${N}_{\mathrm{per}}=13$ cnoidal wave that forms at $\alpha =0$ does not have time to evolve into anything else before the system arrives at $\alpha =6$. Subsequently, the unstable waveform evolves into a single soliton. This behavior consistently occurs with different noise realizations. In Figs. 13(b) and 13(c), we show the evolution when the transition between $\alpha =0$ and $\alpha =6$ occurs between $t=100$ and $t=1.201\times {10}^{5}$ ($\tau =0.55\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ms}$). We then allow the solution to evolve up to $t=2\times {10}^{6}$ ($\tau =9.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ms}$). In this case, the evolution is sufficiently slow for several transcritical bifurcations to take place as $\alpha $ increases. The ${N}_{\mathrm{per}}=8$ cnoidal wave that appears at the end of the initial evolution is only weakly unstable, i.e., $\mathrm{max}({\lambda}_{R})\simeq 6.7\times {10}^{-5}$. The solution that we show ultimately collapses to a continuous wave at around 8 ms, as shown in Fig. 13(c). We have run this simulation with four other noise realizations. In all cases, the solution collapses to continuous waves. The computational integration time is quite long in this case, although the physical time is less than a second. Again, the use of the dynamical approach to point to the presence of an instability was critical in integrating sufficiently long enough to detect it.We find that the final state of the system depends on how quickly the system moves through the parameter space, as well as the trajectory.

## 5. CONCLUSION

Single solitons are a special case (${N}_{\mathrm{per}}=1$) of cnoidal waves, more commonly referred to in the microresonator community as Turing rolls. We have determined the parameter ranges within which different periodicities are stable by solving the LLE for parameters that are relevant for microresonators. We have also described methods for accessing these different periodicities. We have demonstrated that cnoidal waves with ${N}_{\mathrm{per}}>1$ can have a broad bandwidth, comparable to single solitons, while at the same time they are easier to access and use the pump more efficiently. In this limit, they are effectively a periodic train of solitons or a soliton crystal.

In order to determine the ranges of stable operation in the parameter space and to optimize the cnoidal wave parameters, we used a set of software algorithms that are based on dynamical systems theory. While the basic ideas are old, we are the only research group in optical sciences of which we are aware that has implemented them in software in a form that makes it possible to rapidly and accurately determine where in the experimentally adjustable parameter space stable waveforms exist. As part of these algorithms, we calculate the dynamical spectrum, which is the set of eigenvalues for the linearized operator about a stationary waveform. In addition to its utility in determining the stability of the stationary waveform, the dynamical spectrum also allows us to calculate the time scale on which an instability will manifest itself.

We focus, in particular, on the ${N}_{\mathrm{per}}=8$ cnoidal waves for a normalized mode circumference $L=50$, which corresponds approximately to the experimental parameters of [70]. We describe the mechanisms by which the cnoidal wave can become unstable or cease to exist.

The cnoidal waves for a fixed device length occupy an approximately U-shaped band in the pump power detuning parameter space. Below this band, only continuous waves are stable. Above this band, only breathers or chaotic solutions exist. For a normalized detuning $\alpha <41/30$, cnoidal waves are easily accessed by simply raising the pump power. That is the case for the ${N}_{\mathrm{per}}=8$cnoidal wave. Once a cnoidal wave has been accessed, its parameters can be changed by moving inside its stable region in the parameter space. Typically, cnoidal waves with several different periodicities can exist for the same system parameters. We have found that they all have comparable bandwidths. In particular, that is the case for cnoidal waves that exist for the same set of parameters as single solitons.

The most common way to obtain single solitons is to start with negative detunings, where continuous waves or high-periodicity cnoidal waves exist. The system is then red-detuned through the high-periodicity cnoidal wave region into the chaotic region. After further red detuning, the system moves into the region where low-periodicity cnoidal waves, including single solitons, exist. This path through the chaotic region makes it possible for multiple solitons to appear, whose number and spacing vary randomly from shot to shot. Since the region where single solitons are stable nearly overlaps with the regions where continuous waves and several other low-periodicity cnoidal waves exist, it is hard to deterministically ensure that only a single soliton will appear.

We have investigated an alternative approach, in which the system moves through a U-shaped trajectory in the parameter space where cnoidal waves are stable. We showed that the stationary cnoidal wave that is obtained depends on the time scale at which the system moves along a trajectory through this space, as well as the trajectory itself.

All of the theoretical work presented here is based on the LLE. We have not discussed noise issues in any detail, and we have not discussed higher-order dispersion or thermal effects. It is reasonable to suppose that techniques that use dispersive waves to increase the bandwidth of single solitons would be useful for cnoidal waves and could perhaps stabilize them against the effects of noise, in combination with acoustic effects, as is the case in some fiber lasers with multiple pulses in the cavity [20].

It is already known that thermal effects have a profound effect on the region in the parameter space where single solitons are stable [30,35,73,74]. Thermal effects can differ significantly depending on the material system, geometry, and operating temperature. For silicon-nitride microresonators operating at or near room temperature, preliminary results indicate that thermal effects distort the U-shaped region where stable cnoidal waves exist, but do not change the basic results. Recently, Moille *et al.* [75] have studied microresonators at cryogenic temperatures where thermal effects can be almost eliminated, and our theoretical results apply directly. Regardless of the material system, geometry, and operating temperature, our results are a necessary starting point for further study.

## Funding

National Science Foundation (ECCS-1807272, ECCS-1809784); Air Force Office of Scientific Research (FA9550-15-1-0211); Aviation and Missile Research, Development, and Engineering Center (AMRDEC); Defense Advanced Research Projects Agency (W31P4Q-14-1-0002).

## Acknowledgment

We thank Yanne Chembo and Omri Gat for useful discussions. We thank Prem Kumar for encouraging and supporting this work. A portion of our computing work was carried out at the UMBC High Performance Computing Center (https://hpcf.umbc.edu/).

See Supplement 1 for supporting content.

## REFERENCES

**1. **D. J. Korteweg and G. de Vries, “On the change of form of long waves in a rectangular canal, and on a new type of long stationary waves,” Philos. Mag. **39**(240), 422–443 (1895). [CrossRef]

**2. **G. B. Whitham, *Linear and Nonlinear Waves* (Wiley, 1974).

**3. **H. Schamel, “Analytical BGK modes and their modulational instability,” J. Plasma Phys. **13**, 139–145 (1975). [CrossRef]

**4. **A. E. Walstead, “A study of nonlinear waves described by the cubic nonlinear Schrödinger equation,” Report UCRL-52920 (1980).

**5. **A. R. Osbourne, *Nonlinear Ocean Waves and the Inverse Scattering Transform* (Academic, 2010).

**6. **T. Kaladze and S. Mahmood, “Ion-acoustic cnoidal waves in plasmas with warm ions and kappa distributed electrons and positrons,” Phys. Plasmas **21**, 032306 (2014). [CrossRef]

**7. **G. P. Agrawal, *Nonlinear Fiber Optics* (Academic, 2001).

**8. **S. T. Cundiff, “Soliton dynamics in mode-locked lasers,” in *Dissipative Solitons*, N. Akhmediev and A. Ankiewicz, eds. (Springer, 2005).

**9. **A. B. Matsko, A. A. Savchenko, W. Liang, V. S. Ilchenko, D. Seidel, and L. Maleki, “Mode-locked Kerr frequency combs,” Opt. Lett. **36**, 2845–2847 (2011). [CrossRef]

**10. **S. Coen, H. G. Randle, T. Sylvestre, and M. Erkintalo, “Modeling of octave-spanning Kerr frequency combs using a generalized mean-field Lugiato–Lefever equation,” Opt. Lett. **38**, 37–39 (2013). [CrossRef]

**11. **Y. K. Chembo and C. R. Menyuk, “Spatiotemporal Lugiato–Lefever formalism for Kerr-comb generation in whispering-gallery-mode resonators,” Phys. Rev. A **87**, 053852 (2013). [CrossRef]

**12. **P. F. Byrd and M. D. Friedman, *Handbook of Elliptic Integrals for Engineers and Scientists* (Springer, 1971).

**13. **L. M. Milne-Thomson, “Jacobian elliptic functions and theta functions,” in *Handbook of Mathematical Functions*, M. Abramowitz and I. A. Stegun, eds. (National Bureau of Standards, 1972).

**14. **N. Akhmediev and A. Ankiewicz, “Dissipative solitons in the complex Ginzburg–Landau and Swift–Hohenberg equations,” in *Dissipative Solitons*, N. Akhmediev and A. Ankiewicz, eds. (Springer, 2005).

**15. **V. Aleshkevich, Y. Kartashov, and V. Vysloukh, “Cnoidal wave compression by means of multisoliton effect,” Opt. Commun. **185**, 305–314 (2000). [CrossRef]

**16. **Y. V. Kartashov, V. A. Vysloukh, and L. Torner, “Cnoidal wave patterns in quadratic nonlinear media,” Phys. Rev. E **67**, 066612 (2003). [CrossRef]

**17. **D. J. Kedziora, A. Ankiewicz, and N. Akhmediev, “Rogue waves and solitons on a cnoidal background,” Eur. Phys. J. Spec. Top. **223**, 43–62 (2014). [CrossRef]

**18. **J. P. Gordon, “Interaction forces among optical solitons,” Opt. Lett. **8**, 596–598 (1983). [CrossRef]

**19. **F. Leo, S. Coen, P. Kockaert, S.-P. Gorza, P. Emplit, and M. Haelterman, “Temporal cavity solitons in one-dimensional Kerr media as bits in an all-optical buffer,” Nat. Photonics **4**, 471–476 (2010). [CrossRef]

**20. **M. Pang, X. Jiang, W. He, G. K. L. Wong, G. Onishchukov, N. Y. Joly, G. Ahmed, C. R. Menyuk, and P. St.J. Russell, “Stable subpicosecond soliton fiber laser passively mode-locked by gigahertz acoustic resonance in photonic crystal fiber core,” Optica **2**, 339–342 (2015). [CrossRef]

**21. **M. R. E. Lamont, Y. Okawachi, and A. L. Gaeta, “Route to stabilized ultrabroadband microresonator-based frequency combs,” Opt. Lett. **38**, 3478–3481 (2013). [CrossRef]

**22. **T. Herr, V. Brasch, J. D. Jost, N. M. Kondratiev, M. L. Gorodetsky, and T. J. Kippenberg, “Temporal solitons in microresonators,” Nat. Photonics **8**, 145–152 (2014). [CrossRef]

**23. **X. Yi, Q.-F. Yang, K. Y. Yang, M.-G. Suh, and K. Vahala, “Soliton frequency comb at microwave rates in a high-Q silica microresonator,” Optica **2**, 1078–1085 (2015). [CrossRef]

**24. **P.-H. Wang, J. A. Jaramillo-Villegas, X. Yuan, X. Xue, C. Bao, D. E. Leaird, M. Qi, and A. M. Weiner, “Intracavity characterization of micro-comb generation in the single-soliton regime,” Opt. Express **24**, 10890–10897 (2016). [CrossRef]

**25. **A. Coillet, I. Balakireva, R. Henriet, K. Saleh, L. Larger, J. M. Dudley, C. R. Menyuk, and Y. K. Chembo, “Azimuthal Turing patterns, bright and dark solitons in Kerr combs generated with whispering-gallery-mode resonators,” IEEE Photon. J. **5**, 6100409 (2013). [CrossRef]

**26. **J. Pfeifle, A. Coillet, R. Henriet, K. Saleh, P. Schindler, C. Weimann, W. Freude, I. V. Balakireva, L. Larger, C. Koos, and Y. Chembo, “Optimally coherent Kerr combs generated with crystalline whispering gallery mode resonators for ultrahigh capacity fiber communications,” Phys. Rev. Lett. **114**, 093902 (2015). [CrossRef]

**27. **D. C. Cole, E. S. Lamb, P. Del’Haye, S. A. Diddams, and S. B. Papp, “Soliton crystals in Kerr resonators,” Nat. Photonics **11**, 671–676 (2017). [CrossRef]

**28. **C. Joshi, A. Klenner, Y. Okawachi, M. Yu, K. Luke, X. Ji, M. Lipson, and A. L. Gaeta, “Counter-rotating cavity solitons in a silicon nitride microresonator,” Opt. Lett. **43**, 547–550 (2018). [CrossRef]

**29. **V. Brasch, M. Geiselmann, T. Herr, G. Lihachev, M. H. P. Pfeiffer, M. L. Gorodetsky, and T. J. Kippenberg, “Photonic chip-based optical frequency comb using soliton Cherenkov radiation,” Science **351**, 357–360 (2016). [CrossRef]

**30. **Q. Li, T. C. Briles, D. A. Westly, T. E. Drake, J. R. Stone, B. R. Ilic, S. A. Diddams, S. B. Papp, and K. Srinivasan, “Stably accessing octave-spanning microresonator frequency combs in the soliton regime,” Optica **4**, 193–203 (2017). [CrossRef]

**31. **V. Brasch, E. Lucas, J. D. Jost, M. Geiselmann, and T. J. Kippenberg, “Self-referenced photonic chip soliton Kerr frequency comb,” Light Sci. Appl. **6**, e16202 (2017). [CrossRef]

**32. **D. T. Spencer, T. Drake, T. C. Briles, J. Stone, L. C. Sinclair, C. Frederick, Q. Li, D. Westly, B. R. Ilic, A. Bluestone, N. Volet, T. Komljenovic, L. Chang, S. H. Lee, D. Y. Oh, M.-G. Suh, K. Y. Yang, M. H. P. Pfeiffer, T. J. Kippenberg, E. Norberg, L. Theogarajan, K. Vahala, N. R. Newbury, K. Srinivasan, J. E. Bowers, S. A. Diddams, and S. Papp, “An optical-frequency synthesizer using integrated photonics,” Nature **557**, 81–85 (2018). [CrossRef]

**33. **C. Joshi, J. K. Jang, L. Luke, X. Ji, S. A. Miller, A. Klenner, Y. Okawachi, M. Lipson, and A. L. Gaeta, “Thermally controlled comb generation and soliton modelocking in microresonators,” Opt. Lett. **41**, 2565–2568 (2016). [CrossRef]

**34. **X. Yi, Q.-F. Yang, K. Y. Yang, and K. Vahala, “Active capture and stabilization of temporal solitons in microresonators,” Opt. Lett. **41**, 2037–2040 (2016). [CrossRef]

**35. **H. Guo, M. Karpov, E. Lucas, A. Kordts, M. H. P. Pfeiffer, V. Brasch, G. Lihachev, V. E. Lobanov, M. L. Gorodetsky, and T. J. Kippenberg, “Universal dynamics and deterministic switching of dissipative Kerr solitons in optical microresonators,” Nat. Phys. **13**, 94–102 (2017). [CrossRef]

**36. **C. Bao, Y. Xuan, S. Wabnitz, M. Qi, and A. M. Weiner, “Spatial mode-interaction induced single soliton generation in microresonators,” Optica **4**, 1011–1015 (2017). [CrossRef]

**37. **Z. Kang, F. Li, J. Yuan, K. Nakkeeran, J. N. Kutz, Q. Wu, C. Yu, and P. K. A. Wai, “Deterministic generation of single soliton Kerr frequency comb in microresonators by a single shot pulsed trigger,” Opt. Express **26**, 18563–18577 (2018). [CrossRef]

**38. **H. Taheri, A. A. Eftekhar, K. Wiesenfeld, and A. Adibi, “Soliton formation in whispering-gallery-mode resonators via input phase modulation,” IEEE Photon. J. **7**, 2200309 (2015). [CrossRef]

**39. **G. D’Aguanno and C. R. Menyuk, “Nonlinear mode coupling in whispering-gallery-mode resonators,” Phys. Rev. A **93**, 043820 (2016). [CrossRef]

**40. **R. J. Weiblen and I. Vurgaftman, “Bichromatic pumping in mid-infrared microresonator frequency combs with higher-order dispersion,” Opt. Express **27**, 4238–4260 (2019). [CrossRef]

**41. **D. C. Cole, J. R. Stone, M. Erkintalo, K. Y. Yang, X. Yi, K. J. Vahala, and S. B. Papp, “Kerr-microresonator solitons from a chirped background,” Optica **5**, 1304–1310 (2018). [CrossRef]

**42. **S. Zhang, M. Silver, L. Del Bino, F. Copie, M. T. M. Woodley, G. N. Ghalanos, A. Ø. Svela, N. Moroney, and P. Del’Haye, “Sub-milliwatt-level microresonator solitons with extended access range using an auxiliary laser,” Optica **6**, 206–212 (2019). [CrossRef]

**43. **S. Wang, A. Docherty, B. S. Marks, and C. R. Menyuk, “Boundary tracking algorithms for determining stability of mode-locked pulses,” J. Opt. Soc. Am. B **31**, 2914–2930 (2014). [CrossRef]

**44. **S. Wang, B. S. Marks, and C. R. Menyuk, “Nonlinear stabilization of high-energy and ultrashort pulses in passively modelocked lasers,” J. Opt. Soc. Am. B **33**, 2596–2601 (2016). [CrossRef]

**45. **Y. Shen, J. Zweck, S. Wang, and C. R. Menyuk, “Spectra of short pulse solutions of the cubic-quintic Ginzburg–Landau equation near zero dispersion,” Stud. Appl. Math. **137**, 238–255 (2016). [CrossRef]

**46. **C. R. Menyuk and S. Wang, “Spectral methods for determining the stability and noise performance of passively modelocked lasers,” Nanophotonics **5**, 332–350 (2016). [CrossRef]

**47. **J. C. Maxwell, *On the Stability of the Motion of Saturn’s Rings* (MacMillan, 1859).

**48. **H. A. Haus, “Parameter ranges for CW passive modelocking,” IEEE J. Quantum Electron. **12**, 169–176 (1976). [CrossRef]

**49. **H. A. Haus, “Mode-locking of lasers,” IEEE J. Sel. Top. Quantum Electron. **6**, 1173–1185 (2000). [CrossRef]

**50. **J. N. Kutz, “Mode-locked soliton lasers,” SIAM Rev. **48**, 629–678 (2006). [CrossRef]

**51. **A. B. Matsko and L. Maleki, “On timing jitter of mode locked Kerr frequency combs,” Opt. Express **21**, 28862–28876 (2013). [CrossRef]

**52. **C. Godey, I. V. Balakireva, A. Coillet, and Y. K. Chembo, “Stability analysis of the spatiotemporal Lugtiato–Lefever model for Kerr optical frequency combs in the anomalous and normal dispersion regimes,” Phys. Rev. A **89**, 063814 (2014). [CrossRef]

**53. **I. V. Barashenkov and Y. S. Smirnov, “Existence and stability for the ac-driven, damped nonlinear Schrödinger solitons,” Phys. Rev. E **54**, 5707–5725 (1996). [CrossRef]

**54. **I. V. Barashenkov, Y. S. Smirnov, and N. V. Alexeeva, “Bifurcation to multisoliton complexes in the ac-driven, damped nonlinear Schrödinger equation,” Phys. Rev. E **57**, 2350–2364 (1998). [CrossRef]

**55. **M. Haelterman, S. Trillo, and S. Wabnitz, “Dissipative modulation instability in a nonlinear dispersive ring cavity,” Opt. Commun. **91**, 401–407 (1992). [CrossRef]

**56. **L. A. Lugiato and R. Lefever, “Spatial dissipative structures in passive optical systems,” Phys. Rev. Lett. **58**, 2209–2211 (1987). [CrossRef]

**57. **Z. Qi, G. D’Aguanno, and C. R. Menyuk, “Nonlinear frequency combs generated by cnoidal waves in microring resonators,” J. Opt. Soc. Am. B **34**, 785–794 (2017). [CrossRef]

**58. **A. M. Turing, “The chemical basis of morphogenesis,” Philos. Trans. R. Soc. B **237**, 37–72 (1952). [CrossRef]

**59. **Y. Okawachi, K. Saha, J. S. Levy, Y. H. Wen, M. Lipson, and A. L. Gaeta, “Octave-spanning frequency comb generation in a silicon nitride chip,” Opt. Lett. **36**, 3398–3400 (2011). [CrossRef]

**60. **S. Wang, A. Docherty, B. S. Marks, and C. R. Menyuk, “Comparison of numerical methods for modeling laser mode locking with saturable gain,” J. Opt. Soc. Am. B **30**, 3064–3074 (2013). [CrossRef]

**61. **J. Nocedal and S. Wright, *Numerical Optimization* (Springer, 2006).

**62. **S. H. Strogatz, *Nonlinear Dynamics and Chaos: With Applications to Biology, Chemistry, and Engineering* (Westview, 1994).

**63. **M. W. Hirsch, S. Smale, and R. L. Devaney, *Dynamical Systems and an Introduction to Chaos* (Academic, 2013).

**64. **I. Stakgold, *Boundary Value Problems of Mathematical Physics* (MacMillan, 1967), Vol. 1, Chap. 2.

**65. **G. Helmberg, *Introduction to Spectral Theory in Hilbert Space* (Dover, 1997).

**66. **N. N. Akhmediev, A. Ankiewicz, and J. M. Soto-Crespo, “Multisoliton solutions of the complex Ginzburg–Landau equation,” Phys. Rev. Lett. **79**, 4047–4051 (1997). [CrossRef]

**67. **A. B. Matsko, A. A. Savchenkov, N. Yu, and L. Maleki, “Whispering-gallery-mode resonators as frequency references. I. Fundamental limitations,” J. Opt. Soc. Am. B **24**, 1324–1335 (2007). [CrossRef]

**68. **W. Liang, D. Eliyahu, V. S. Ilchenko, A. B. Matsko, D. Seidel, and L. Maleki, “High spectral purity Kerr frequency comb radio frequency photonic oscillator,” Nat. Commun. **6**, 7957 (2015). [CrossRef]

**69. **Z. Qi, S. Wang, J. A. Jaramillo-Villegas, M. Qi, A. M. Weiner, G. D’Aguanno, and C. R. Menyuk, “Stability of cnoidal wave frequency combs in microresonators,” in *Conference on Lasers and Electro-Optics* (Optical Society of America, 2018), paper SF2A.6.

**70. **J. A. Jaramillo-Villegas, X. Xue, P.-H. Wang, D. E. Leaird, and A. M. Weiner, “Deterministic single soliton generation and compression in microring resonators avoiding the chaotic region,” Opt. Express **23**, 9618–9626 (2015). [CrossRef]

**71. **P. Parra-Rivas, D. Gomila, M. A. Matías, S. Coen, and L. Gelens, “Dynamics of localized and patterned structures in the Lugiato–Lefever equation determine the stability and shape of optical frequency combs,” Phys. Rev. A **89**, 043813 (2014). [CrossRef]

**72. **D. J. Kaup, “Perturbation theory for solitons in optical fibers,” Phys. Rev. A **42**, 5689–5694 (1990). [CrossRef]

**73. **C. Bao, Y. Xuan, J. A. Jaramillo-Villegas, D. E. Leaird, M. Qi, and A. M. Weiner, “Direct soliton generation in microresonators,” Opt. Lett. **42**, 2519–2522 (2017). [CrossRef]

**74. **J. A. Jaramillo-Villegas, C. Wang, P.-H. Wang, C. Bao, Y. Xuan, K. Han, D. E. Leaird, M. Qi, and A. M. Weiner, “Experimental characterization of pump power and detuning in microresonator frequency combs,” in *Latin America Optics & Photonics Conference* (2016), paper LTh3B.2.

**75. **G. Moille, X. Lu, A. Rao, Q. Li, D. A. Westly, L. Ranzani, S. B. Papp, M. Soltani, and K. Srinivasan, “Kerr microresonator soliton frequency combs at cryogenic temperatures,” arXiv 1906.0654 (2019).