Dirac solitons in optical microresonators

Mode-coupling-induced dispersion has been used to engineer microresonators for soliton generation at the edge of the visible band. Here, we show that the optical soliton formed in this way is analogous to optical Bragg solitons and, more generally, to the Dirac soliton in quantum field theory. This optical Dirac soliton is studied theoretically, and a closed-form solution is derived in the corresponding conservative system. Both analytical and numerical solutions show unusual properties, such as polarization twisting and asymmetrical optical spectra. The closed-form solution is also used to study the repetition rate shift in the soliton. An observation of the asymmetrical spectrum is analysed using theory. The properties of Dirac optical solitons in microresonators are important at a fundamental level and provide a road map for soliton microcomb generation in the visible band. The intriguing properties of a new type of ultrashort light pulse in microresonators are described. Called a Dirac soliton, Heming Wang and colleagues at the California Institute of Technology say that understanding the properties of Dirac solitons is important at the fundamental level and also for providing a road map for generating ultrashort soliton pulses in the visible light range. The latter is relevant for applications such as miniature optical clocks and medical imaging techniques, for example. The solitons are generated in micron-scale devices called microresonators. And the soliton properties are controlled by the microresonator shape and choice of material. Wang and his colleagues also called attention to the importance of broadband nonlinear coupling in the microresonators for generation of the new Dirac soliton, which their equations show exhibits polarization twisting and asymmetrical optical spectra.


Introduction
Soliton mode locking in microresonators 1 provides a pathway for the miniaturization of frequency comb systems 2 . The dissipative solitons 3 formed in the resulting microcombs 4 are coherently pumped 5 and were first observed in optical fibre cavities 6 . In microresonators, such Kerr solitons (KSs) have been realized in a wide range of geometries and material systems [7][8][9][10][11][12][13][14] . Soliton microcomb devices have been tested in diverse system demonstrations, including spectroscopy [15][16][17] , coherent communications 18 , range detection [19][20][21] , optical frequency synthesis 22 , exoplanet studies 23,24 , and optical clocks 25 . Progress towards integration of the microcomb with pump and other control functions is also being made [26][27][28] . Modal coupling, wherein distinct mode families experience frequency degeneracy analogous to an energy level crossing 29 , is an important feature of soliton formation in microresonators. Such crossings impart structure to the soliton spectral envelope 30 and are responsible for an intriguing range of microcomb phenomena of both scientific and technical importance, including dispersive wave emissions 9,31 , dark soliton formation 32 , pump noise isolation 33 , improved pumping efficiency 34,35 , and dispersion engineering for near-visible emissions 36,37 .
Here, a new type of soliton in microresonators, termed Dirac solitons (DSs), is shown to result from broadband modal coupling. The name originates from the nonlinear Dirac equations, which govern the dynamics of these solitons and are discussed below. A similar soliton has been theoretically studied in fibre Bragg gratings 38,39 and later experimentally observed 40 . In these Bragg solitons, forward and backward propagating waves are coupled by a periodic structure, and a Dirac-like model has been applied to understand these systems 41,42 . As shown recently for broadband coupling in a dimer resonator system 43 , the co-existence of coupling and nonlinearity changes the solution behaviour qualitatively, and a full understanding requires a non-perturbative approach. We show that broadband nonlinear coupling results in a range of new phenomena in the Dirac soliton system, including polarization twisting along the soliton and asymmetrical soliton comb spectra. A closed-form expression for DSs is derived by solving the Lugiato-Lefever equation (LLE) 5,44 augmented with mode coupling. Curiously, the requisite coupling conditions for DS generation have been obtained experimentally for near-visible 36 and 1-μm-band 37 soliton generation. As shown here, these experiments were performed in a regime where Dirac solitons collapse into conventional Kerr solitons. New data from the 778-nmband measurement featuring asymmetrical spectra will be presented, showing initial deviation away from conventional Kerr soliton behaviour.

Polarization mode coupling and coupled LLEs
To illustrate the features of DSs, we consider the specific case of a circularly symmetric (whispering-gallery) resonator that has an initial reflection plane of symmetry in the plane of the resonator. The resulting geometry supports transverse-electric (TE) and transverse-magnetic (TM) mode families (Fig. 1a) that are symmetrical and antisymmetrical, respectively, with respect to the reflection plane. A pair of TE and TM modes can have accidental degeneracy for a particular wavenumber. By breaking the reflection geometry, it becomes possible to lift the degeneracy and create an avoided crossing 36 . In effect, the original modes are strongly coupled, and the eigenmodes of the non-symmetric system are hybrid modes, as shown in Fig. 1b. Loosely speaking, the resulting hyperbolic shape of the eigenfrequency dispersion creates an anomalous dispersion window that is suitable for soliton generation. However, we note that this dispersion is not associated with a single mode family across the entire avoided crossing. Indeed, the mode composition changes when the wavenumber increases, evolving from TM to TE mode for the upper branch or vice versa for the lower branch.
The standard form of the LLE for one transverse mode family (denoted by mode 1) describes the temporal soliton dynamics in a microresonator 44 : Here, E 1 is the slowly varying amplitude in a co-moving frame normalized to optical energy, defined via E 1 = E 1 A 1 , where E 1 is the electric field and A 1 is the normalized vector field distribution. The frequency detuning δω = ω c − ω p is the frequency difference between the cavity resonant frequency ω c and pump frequency ω p . The Kerr nonlinear coefficient is g 11 ¼ n ð2Þ ω c c=ðn 2 V 11 Þ, with speed of light in vacuum c, refractive index n, Kerr nonlinear index n (2) and mode volume V 11 ¼ ð R jA 1 j 2 dV Þ 2 = R jA 1 j 4 dV . κ 1 is the energy loss rate, and f 1 is the pumping term for mode 1. The linear dispersion operatorL 1 describes mode dispersion and can be Taylor expanded asL 1 % ÀiD 1;1 ∂ θ À D 2;1 ∂ 2 θ =2, where θ is the angular coordinate and D 1,1 /(2π) and D 2,1 / (2π) are the free spectral range (FSR) and second-order dispersion (proportional to the group velocity dispersion (GVD)), respectively, for mode 1. In the case of a conservative system (κ 1 = 0 and f 1 = 0) and D 1,1 = 0 (i.e., choosing a co-moving reference frame), the exact soliton solution to Eq. (1) is given by: which is also commonly used as an ansatz to describe a KS 7 .
To generalize the above LLE to the coupled pair of TE and TM modes (modes 1 and 2, respectively), we introduce mode coupling as well as cross-phase modulation into the equations. The following two-mode coupled LLE results: Here, g c > 0 gives the coupling rate between the two (originally uncoupled) modes, g ij ¼ n ð2Þ ω c c=ðn 2 V ij Þ (i, j = 1, 2), and the mode volumes V 11 and V 22 and cross mode volume V 12 are defined as  Fig. 1 Principle of mode hybridization and Dirac solitons. a For a symmetric resonator cross section (top-left insets), TE and TM modes within the resonator can become accidentally degenerate at the same wavenumber. The bottom (right) inset depicts the TE (TM) mode electric field directions. b For an asymmetric resonator cross section (top-left inset), degeneracy is lifted, and an avoided crossing is created. The left and right insets depict the hybridized mode electric field directions. c Schematic of balancing nonlinear and dispersion effects for a KS. d Schematic of balancing nonlinear, coupling and group velocity difference effects for one of the components in a DS frequencies are measured) is chosen as the degeneracy frequency. In the symmetric co-moving frame (moving at a speed corresponding to D 1 ¼ ðD 1;1 þ D 1;2 Þ=2), the group velocities of the two resulting hybrid modes become antisymmetric with δD 1 ¼ jD 1;1 À D 1;2 j=2. Here, second-and higher-order dispersions are ignored, as the couplinginduced dispersion of the eigenfrequency is typically one order of magnitude larger than the intrinsic mode dispersion. The nonlinear terms include self-phase and cross-phase modulation. Other four-wave mixing terms that induce nonlinear coupling, such as jE 1 j 2 E 1 E Ã 2 and E 2 1 ðE Ã 2 Þ 2 , have been dropped because these are either forbidden by reflection symmetry or strongly suppressed by the phase mismatching of the underlying modes.
The LLE Eq. (1) (without loss and pump terms) is known as the nonlinear Schrödinger equation in theoretical physics. Similarly, the coupled LLE Eq. (3) presented here are a generalization of the nonlinear Dirac equations. When only cross-phase modulation is considered (g 11 = g 22 = 0), Eq. (3) reduces to the massive Thirring model in quantum field theory 45 , which is known to support Dirac solitons 46,47 . With equal but nonzero self-phase modulations (g 11 = g 22 ), the Bragg soliton solution 38,39 is recovered, which has been realized in fibre Bragg grating systems 40 . On the other hand, when second-order dispersion is present and is much stronger than the effect induced by linear inter-mode coupling within the band being considered, Eq. (3) becomes the vector soliton model in birefringent systems, where soliton solutions are also known [48][49][50] . We note that the vector soliton relies on both modes having anomalous dispersion, while anomalous dispersion is not required in the DS model.
Before proceeding to solve Eq. (3), it is helpful to understand why a DS solution exists in the absence of second-order dispersion. The conventional KS is a delicate balance of the Kerr nonlinear effect, which creates chirping within a pulse, and the anomalous mode dispersion, which cancels the chirping effect (Fig. 1c). In the DS system, the hyperbolic-shaped upper-branch eigenfrequency spectrum (as in Fig. 1b) resembles a spectrum with anomalous dispersion. While this "dispersion" plays the same role as conventional dispersion and constitutes the foundation for the generation of the DS, this viewpoint only holds when the pumping frequency (and soliton spectrum) is close to this branch. Indeed, in this case, it will be shown that it reduces to the conventional KS. In general, and as noted in the introduction, dispersion is only locally well-defined in this spectrum, because the mode composition of the hybrid mode can change rapidly with respect to the wavenumber. Correspondingly, the dispersion interpretation fails for the general DS, and the coupling effect must be treated nonperturbatively. These rapid composition changes in the hybridized modes redistribute pulse energy in the frequency domain and produce a new contribution to chirping within the pulse, which manifests as phase differences between the two-mode components of the pulse. Coupling then makes the two components interfere differently at different positions and leads to both chirping and pulse shifting. These effects are delicately cancelled by nonlinear effects and group velocity differences, respectively (Fig. 1d), and maintain the DS pulse shape without anomalous dispersion. In the language of field theory, anomalous mode dispersion provides a positive "mass" for the KS field, which then becomes a well-defined non-relativistic field theory. The "mass" of the DS field is the inter-mode coupling, and the mode spectrum corresponds to a relativistic field theory.

Closed-form soliton solutions
To obtain an analytical solution for the coupled LLE, we consider a conservative system by setting κ 1 = κ 2 = 0 and f 1 = f 2 = 0. Equation (3) can then be solved by finding the invariants associated with the system (see 'Materials and methods'). The closed-form single bright soliton solution for Eq. (3), without periodic boundary conditions, can be obtained as: where E 1 (E 2 ) takes the upper (lower) sign in all instances of ± or ∓, v is the repetition rate shift in the symmetric comoving frame, G is the combined nonlinear coefficient, γ is a phase exponent related to v,ξ is the reduced detuning and θ is the reduced coordinate. While dark soliton solutions and bright soliton on a background solutions can also be found in the same conservative system (see Supplementary  Information), we will focus on this bright soliton solution and refer to it as the DS. In the following discussion of DS properties, we take the special case g 11 = g 22 (i.e., additional exchange symmetry between the modes) and v = 0 (i.e., the pulse is stationary in the symmetric co-moving frame), and the general solution simplifies to: Figure 2a shows the offset frequency for the hybrid modes, defined as ω off ¼ ω k À ω c À kD 1 where k is the relative wavenumber (the difference between the absolute wavenumber and the wavenumber at the degeneracy point) and ω k is the mode eigenfrequency at k. Due to the square roots in the special solution (Eq. (5)), the soliton detuning range can lie only in the band gap created by the avoided crossing. This phenomenon can be intuitively understood, as none of the comb lines can have the same frequency as the resonator modes, which would otherwise create infinite amplitudes on the modes due to perfect resonance with no loss. Geometrically, the resonance of the soliton can be described by the linear equation ω DS = −δω + vk, and the line cannot intersect the two hyperbolas of mode frequencies on the mode spectrum plot. The same also holds for the general solution (Eq. (4)) and is depicted in Fig. 2b (see 'Materials and methods'). As a result, the soliton cannot have a group velocity faster than the first mode or slower than the second mode. This argument does not apply to dark solitons and bright solitons on a background, where the resonance line always intersects the mode spectrum (see Supplementary  Information).
To understand the properties of the DS, we consider some special cases of detuning (marked in Fig. 2a). The first case is when δω approaches −g c , where the resonance line is close to the upper branch. By taking appropriate limits, Eq. (5) reduces to: which is exactly in the form of a conventional KS. A comparison of the exact DS near δω = −g c and the limiting KS is shown in Fig. 2c. The appearance of the sech-shaped KS here is not a coincidence. The effective nonlinear coefficient for a single component is g 11 + g 12 (using jE 1 j 2 ¼ jE 2 j 2 ), the effective detuning is g c + δω, and the curvature of the eigenfrequency (the hybrid-mode equivalent of the second-order dispersion) is δD 2 1 =g c , as derived from coupled mode theory. The reduction of a DS to a KS is straightforward when these quantities are substituted into the KS solution. Thus, if the DS is close to the resonance (which is usually the case when the hybridization coupling g c is large), the eigenfrequency spectrum is locally equivalent to a single mode in terms of dispersion, and mode composition differences for different wavenumbers can be ignored. This phenomenon also applies to the general solution (Eq. (4)) by explicitly reducing the coupled LLE to a single-mode LLE, and the truncation errors can also be estimated (see 'Materials and methods'). The second case is when δω = 0, where the resonance line passes through the degeneracy point. In this case, Eq. (5) simplifies to: and begins to deviate from a hyperbolic secant shape (Fig. 2d). With this analytical solution, it becomes apparent that each component of the wave packet has an overall phase shift when θ goes from −∞ to ∞ (π/2 in this special case). This phase twist within the pulse contributes to the chirping and shifting of the soliton pulse when they are coupled together, as discussed in the previous section. If pulse polarization is considered, it also twists from the start of the pulse to the end of the pulse (Fig. 2e).
The last case we consider is when δω approaches g c , where the resonance line moves towards the lower branch and maximum red detuning is approached. As this phenomenon occurs, the exponential tails of the soliton decay at an increasingly slower rate until finally in the limit δω → g c Eq. (5) becomes: showing that the solution decays polynomially rather than exponentially when θ → ∞ (Fig. 2f). The resulting polynomial tails can potentially enable long-range interactions of the DS. We now turn to the frequency domain profile for Eq. (5) by the Fourier transform, which can also be expressed in closed form using contour integration: where k is the relative wavenumber and Apart from the usual sechshaped envelope, the extra exponential factor causes the spectrum of each component to become asymmetrical around k = 0 (Fig. 2g). This phenomenon can be explained by different mode compositions on the different sides of the spectrum. We note that the spectrum for the total power is still symmetric, which is consistent with the soliton carrying no total momentum in the symmetric comoving frame. In the general v ≠ 0 case, the spectrum for the total power is expected to become asymmetric around the centre frequency of the soliton.

DS with dissipation and repetition rate shifts
In addition to being an exact solution to the conservative system, the solution (Eq. (4)) also serves as a DS ansatz for dissipative cases, similar to how a KS can be approximated by a sech-shaped soliton on a background 7 . As an example, we study the repetition rate shifts associated with the DS when externally pumped by continuous waves. Accordingly, we do not require g 11 = g 22 or v = 0 and return to work with the general solution (Eq. (4)).
The repetition rates of ideal KSs remain constant when pumped in the same mode while the detuning changes, as formulated by the standard LLE. In contrast, real-world KSs may experience additional nonlinear effects, including dispersive wave backactions 33 or Raman effects 51 , which lead to centre frequency changes and repetition rate shifts. The mode hybridization process is similar to a mode crossing in dispersive wave generation where two modes are strongly coupled, and therefore, the DS is also expected to experience repetition rate shifts. As the repetition rate shift parameter v is free in the conservative case, we need to find the conditions that determine v when dissipation is present.
By calculating the momentum integral of the solution, the following criterion is obtained (see 'Materials and methods'): where arg is the argument function and E 1,2 should be substituted by the DS solution. According to the above criterion, the phase twist of each component is essential in determining the repetition rate shift. Intuitively, this concept can be understood as follows: the pulse cannot carry any net momentum in the reference frame of the pumping, and any additional momentum will be damped out by the dissipation. All the above integrations can be carried out in closed form, leading to an equation in v, which can then be solved as the repetition rate shift. For a fixed detuning δω, the repetition rate shift v depends on the ratios of the nonlinear coefficients g 11 , g 22 and g 12 . Figure 3a plots the special case of δω = 0. When the nonlinearity on the second mode increases, the optical field will shift to the first mode to compensate, leading to an increased overall speed of the pulse, and vice versa. As the repetition rate shift results from the imbalance of selfphase modulations, increasing the proportion of g 12 leads to more stability in the repetition rate, while decreasing the proportion of g 12 allows more tunability. We note that, depending on the nonlinear nature of the resonator material and the mode overlap, the cross-phase modulation may be larger or smaller than the self-phase modulation 52,53 . Moreover, theoretical DS solutions exist for almost all combinations of nonlinear coefficients.
For the tunability of the DS repetition rate, Fig. 3b shows the repetition rate shift as the detuning changes. Near δω = −g c , the repetition rate shift approaches zero, which is consistent with the local KS equivalence argument in the previous section. With a more red-detuned δω, the effect of imbalance in the nonlinear coefficients is more apparent, leading to repetition rate changes in the corresponding direction. Simulations of the coupled LLE have also been performed and show that both the simulated pulse shape and the repetition rate shifts agree with the analytical solutions (Fig. 3b). A graphical representation of the repetition rate shift is also shown in Fig. 3c. As an aside, breather-like states 54 have also been observed in the simulations, but the origin of breathing and whether it behaves in the same way as for KS breathers [55][56][57][58] is not yet fully understood.
Although the discussion so far has focused on pumping at the central mode, it can be readily generalized to offcentre pumping by introducing additional detunings into each of the mode families and shifting the spectral centres of the solutions accordingly. The DS offers a novel and controllable way to tune the repetition rate of the frequency combs. Together with existing nonlinear processes for the resonator, the hybridization-induced shift can be tailored to enhance or suppress the overall repetition rate shift with respect to the pump detuning and may find application in optical frequency division or for pump noise isolation, such as what is performed using quiet point operation 33 .

Implementation of Dirac solitons
The wedge resonator 59 is used to induce mode hybridization and Dirac soliton formation. This resonator offers very high-quality factors 60 and independent control over key parameters during the fabrication process (Fig. 4a). A wedge is entirely characterized by three geometric parameters: the diameter D, which depends on the lithographic pattern; the thickness t, which depends on the oxidation growth time of the silicon wafer; and the wedge angle α, which depends on the adhesion between silica and the photoresist used for patterning. In the following, we will fix the resonator diameter as D = 3.2 mm (corresponding to a resonator FSR of 20 GHz at~1550 nm), but we note that this can be readily generalized to resonators of other sizes.
For a symmetrical wedge resonator (α = 90°), the typical simulated effective refractive index n eff versus wavelength is shown in Fig. 4b. At shorter wavelengths, TE1 and TM1 have the highest indices, followed by TE2, TE3 and other highorder modes. Since the electrical fields of the TM modes are along the thickness direction, their indices are more sensitive to changes in the wavelength scale, and the index of TM1 decreases faster than TE2 as the wavelength increases. Eventually, TM1 and TE2 cross, and their relative positions are interchanged at longer wavelengths. However, for α = 90°, no hybridization occurs, as the reflection symmetry prohibits interactions between modes of different parities. On the other hand, if we explicitly break the reflection symmetry of the resonator by decreasing the wedge angle (α < 90°), the original modes will see an asymmetric change in the refractive index profile, which causes them to couple. Such couplings lift the degeneracy, leading to avoided crossing. The two cases are compared in Fig. 4c, where n eff is first converted to the mode number m via m ¼ n eff Dω m =ð2cÞ, where ω m is the resonance (angular) frequency, and then plotted as offset (angular) frequencies ω off ¼ ω m À ω X À ðm À m X ÞD 1 versus the relative mode number m − m X , where the subscript X indicates the quantity at the degeneracy point. We note that the relative mode number has the same role as the relative wavenumber k in the theoretical analyses, except that it is restricted to integer values for periodic boundary conditions.
In view of perturbation theory 61 , the wedge part of the resonator perturbs the underlying symmetrical structure and induces polarization coupling similar to the coupling obtained in directional couplers 62 . Therefore, we expect that the centre wavelength of hybridization λ X is determined by the thickness t, while the wedge angle controls the coupling strength g c . A plot of λ X versus t is shown as the black curve in Fig. 4d. As t is the only geometry scale close to optical wavelengths in the system, we expect that λ X will scale linearly with t, which can be visually verified in the plot. This scaling allows for hybridization to occur at short wavelengths where the dispersion of the original modes (for example, the TE1 mode shown in the figure) is typically normal. A plot of g c versus α is shown in Fig. 4e. While only a particular wavelength (778 nm) is shown, g c depends on the wavelength very weakly, varying less than 5% from wavelengths of 400-1600 nm. The coupling strength scales linearly with α near α = 90°, which can also be independently verified by first-order perturbation theory (see 'Materials and methods'), but the coupling effect eventually saturates at shallow wedge angles because mode profiles cannot "squeeze" into the wedge tip as α decreases. The calculated GVD β 2 is shown in Fig.  4f, which is related to D 2 via β 2 ¼ ÀnD 2 =ðcD 2 1 Þ. Using suitably designed thicknesses and wedge angles greater than 30°, an anomalous dispersion window can be created all the way down to the blue side of the visible spectrum, where simple geometrical dispersion fails to compensate for normal material dispersion.

Demonstration of Dirac solitons
Guided by these design principles, devices that target 1550 nm and 778 nm as their hybridization wavelengths were fabricated. The mode spectra are measured for each device using a tuneable laser and a calibrated Mach-Zehnder interferometer 8 . As expected, each of the devices shows a pair of modes with hyperbolic dispersion and large curvatures (Fig. 5a, c). The local D 2 of the anomalous branch can be fit to give D 2 = 2π × 401 kHz (1550 nm) (Fig. 5b) and D 2 = 2π × 132 kHz (778 nm), corresponding to β 2 = −790 ps 2 km −1 and β 2 = −255 ps 2 km −1 , respectively, which are orders of magnitude larger than the mode intrinsic GVD without hybridization.
Finally, to demonstrate the existence of DSs in wedge resonators with hybridized modes, we generated solitons at 778 nm. The detailed experimental setup and measurement procedures can be found elsewhere 36 . The optical spectrum of the soliton is shown in Fig. 5. A direct sech 2 fit to the spectrum reveals that the frequency components are highly asymmetric around the spectrum centre, which results from the high-order dispersions of the mode and is a typical feature for DSs not being pumped in the crossing centre. Indeed, the analytical DS solution derived here provides a good fit to the measured spectrum, further confirming the above theoretical results and their applicability for soliton modelling in microresonators. The detuning relative to the upper branch can be inferred from the fitting and is approximately 90 MHz, which is less than 5% of the band gap. We did not attempt to increase the detuning further due to pump power limitations. As an aside, we believe that a similar 780-nm soliton generated To further validate the feasibility of DSs in the visible band, we also performed simulations of a DS near 532 nm. A resonator design is shown in Fig. 5e together with its simulated mode dispersions, and solitons can be found in the resonator with different pumping positions (Fig. 5f). The soliton at the crossing centre has a more symmetric spectrum, while the soliton far away from the crossing centre has a more skewed spectrum, as expected.

Discussion
Critical to the generation of solitons in microresonators is the control of mode dispersion. Bright soliton formation requires anomalous dispersion. The curvature of a circular resonator contributes to normal dispersion, shifting the zero-dispersion wavelength towards longer wavelengths as the resonator size decreases. Both this and material dispersion have been managed over a range of wavelengths by using geometrical dispersion introduced through optical waveguide confinement 1 . However, normal dispersion of bulk dielectric materials increases towards shorter wavelengths making visible soliton generation in dielectric microresonators difficult. And the offset of this normal dispersion through waveguide confinement increases unwanted scattering loss, which degrades the resonator Q factor and increases the pumping power (i.e., the comb threshold power varies inversely quadratically with the Q factor 8,63 ). While the use of intra-cavity nonlinear optical processes such as second-or third-harmonic generation and sumfrequency generation can provide a way to bypass normal dispersion for comb generation in the visible band 13,64-67 , managing dispersion by mode coupling provides an alternative approach that can avoid waveguide confinement loss 36,37,[68][69][70][71][72] . The practical management of dispersion for soliton formation at the edge of the visible band 36,37 has been possible using Dirac solitons, and they can provide a way to further extend operation well into the visible region.
In summary, we demonstrated the peculiar characteristics of the Dirac soliton in microresonators by solving the corresponding conservative coupled LLE non-perturbatively and using the exact solution as a soliton ansatz for the hybrid system. The balance between nonlinearity and coupling provides a new viewpoint on the soliton and opens up further directions for the study of soliton dynamics. From an experimental viewpoint, generating DSs in resonator platforms is straightforward, but it might be challenging to observe some of their more unusual properties at large detunings. Decreasing the linear coupling (and thus the band gap) is beneficial for pushing the soliton deeper into the band gap, and tuning the pumping polarization to match the soliton decreases the pumping power requirements. We believe the formalism described here can be readily generalized to other systems, including Dirac solitons in waveguides, solitons on ordinary (non-polarization) mode crossings 33 and counter-propagating solitons with moderate coupling 73 .

Solving the conservative coupled LLE
We copy the conservative coupled LLE here for convenience: We seek soliton solutions in the form of E 1,2 (θ − vt), where v is the repetition rate shift in the symmetric comoving frame, which reduces the partial differential equations to ordinary differential equations: Continuous symmetries of the system result in conservation laws 74 , which can reduce the dimensions of the system. As the system is conservative, we expect the equations will have a Hamiltonian structure. Indeed, the following quantity is conserved when θ is viewed as an evolution coordinate 46 : The conservation of H can be verified by rewriting Another quantity that is conserved is the photon number flow along the θ-axis: The conservation of N can be verified by observing that all the nonlinear terms do not change the individual numbers of particles, while the coupling terms do not change the total number of particles.
For soliton solutions, these two conserved quantities can be determined as H ¼ N ¼ 0 since the solution should vanish exponentially as θ → ∞ without periodic boundary conditions. This determination leads to the following amplitude-phase parametrization of the solutions: which automatically satisfies the N conservation (the negative sign is added for later convenience). The H conservation reads as: 2ðδD 1 ÀvÞ 2 þ g 22 2ðδD 1 þvÞ 2 þ g 12 from which the cosine of the phase difference χ 2 − χ 1 can be solved as: where for convenience, we defined a combined nonlinear coefficient: Turning back to the original equations of evolution along θ, we substitute E 1,2 with the parametrization and split the real and imaginary parts: For the differential equation for ψ 2 , expressing sin(χ 2 − χ 1 ) in terms of ψ 2 gives: which can be integrated (with the boundary condition ψ 2 → 0 as θ → ∞) in terms of elementary functions: where the pulse centre is chosen as θ = 0 without loss of generality and the reduced detuning and coordinate are defined as: As an aside, we also obtain that: The differential equation for χ 1,2 can be integrated after substitution of the above solution for ψ 2 and cos(χ 2 − χ 1 ).
Because the equation has global phase symmetry (E 1;2 ! e iϕ E 1;2 , where ϕ is an arbitrary constant phase), we can fix χ 1 (θ = 0) = 0, which also forces χ 2 = 0 through cosðχ 1 À χ 2 Þj θ¼0 ¼ 1. We obtain: With these results, the soliton solutions can be expressed as: where we introduce the phase exponent: Although we have not been very rigorous for multivalued functions encountered in the calculations, direct substitution shows that the E 1,2 obtained above is indeed a solution to the original conservative LLE when principal branches are used.

Resonance line and the band gap
The general bright soliton solution includes the square root of 1 Àξ 2 , which requires that jξj < 1. Expanded with resonator parameters, this gives: For a fixed v, the inequality gives the detuning range where the solution is well-defined. A quick plot of the range (Fig. 2b) shows that the boundaries are tangent to the mode spectrum curves. Indeed, using coupled mode theory, the frequencies can be described as: where ω + (ω − ) is the eigenfrequency for the lower symmetric branch (upper anti-symmetric branch) and k is the wavenumber. The tangent lines for the upper branch with slope v satisfy: Eliminating k recovers the previous boundaries. Thus, the soliton resonance lines can stay only in the band gap and cannot cut through the band curves. We note that in dissipative cases, v is not fixed but depends on the pumping details (as discussed in the main text), so this point should not be understood as a limitation on detuning when pumping the soliton. Instead, the detuning range should be determined from the momentum constraints imposed on the soliton at fixed k (the longitudinal mode being pumped).

Reduction of a DS to a KS
Following the above discussions on resonance lines, we focus on the case in whichξ ! À1 þ , where the resonance line is almost tangent to the upper branch of the mode spectrum. Taking the limits and expanding the reduced quantities results in: The hyperbolic secant form is now apparent, and to complete the reduction, we explicitly calculate the local quantities of the mode spectrum.
When the resonance line is tangent to the mode spectrum, the wavenumber k can be solved from the previous section: which matches the exponential term. The minimum detuning that can be achieved at this particular m also matches δω min . The local second-order dispersion is given by: where we have eliminated k using v and it matches the dispersion term. The mode composition can be found using coupled mode theory and is given by: which agrees with the prefactors in E 1,2 and is also consistent with the conservation of N. Finally, the effective nonlinear coefficient GðδD 2 1 À v 2 Þ=ð2δD 2 1 Þ can be calculated as a weighted average of the nonlinear coefficients, the weight being the power proportions on each mode derived above. It matches the nonlinear coefficient except for the extra factor ðδD 1 ± vÞ=ð2δD 1 Þ, which is the power ratio of each mode component to the total power and is a result of expressing the solution using components rather than the hybridized field.
To complete the discussion of reducing a DS to a KS, we also present a perturbative approach that is explicitly based on the hybridized field. We begin by defining the following auxiliary fields: We note that while the ψ + component is the normalized linear eigenstate of the lower branch at the wavenumber corresponding to v, the ψ − term defined here is, in general, not the eigenstate of the upper branch, and ψ − and ψ + are not orthogonal (although in the special case ψ + = 0, ψ − becomes proportional to the true field amplitude). Rewriting the conservative coupled LLE in terms of ψ ± results in: ffiffiffiffiffiffiffiffiffiffiffi ffi δD 2 1 Àv 2 p δD1þv δD1Àv g11 2 jψ þ þ ψ À j 2 ðψ þ þ ψ À Þ h þg 12 ðψ 2 þ À ψ 2 À Þψ Ã þ þ δD1Àv δD1þv g22 2 jψ þ À ψ À j 2 ðψ þ À ψ À Þ i where we have substituted δω and θ withξ andθ, respectively, for later convenience. Based on the structure of the above equation, we seek the following solution near ξ ! À1 þ : ψ À $ Oð1 þξÞ 1=2 ; ψ þ $ Oð1 þξÞ Keeping the lowest-order terms gives: Combining gives: Gjψ À j 2 ψ À ¼ 0 which is the steady-state single-mode LLE and its solution is the same as the limit of E 1,2 as derived above.
We thus calculate the first derivative of P with respect to t: where we have used integration by parts to move the spatial derivatives to the conjugated field. After plugging the equations of motion into the integral, all the conservative terms cancel each other out, and the pumping terms vanish by integration by parts. We are left with: Rewriting the above equation using arguments gives: To proceed further, we take the soliton ansatz as the exact solution of the DS derived earlier. In this case, the integration can be carried out analytically: þ γ ± 1 ðπ À arccosξÞξ þ ðγ ± 1Þ ffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Àξ 2 q ! All the quantities can be explicitly expressed in v, and the resulting equation can be solved numerically. In the special case of κ 1 = κ 2 and δω = 0, we haveξ ¼ 0 independent of v, and the criterion is greatly simplified: v δD 1 ¼ Àγ Expanding γ gives a cubic equation in v and is used in the plot of Fig. 3a.

First-order perturbation calculation of mode coupling in wedge resonators
Here, using first-order degeneracy perturbation theory and the integral form of the propagation constant, we derive the mode coupling in wedge resonators as an overlap integral of the unperturbed modes.
For a circular waveguide, the angular momentum number (angular propagation constant) of a mode can be expressed as 61 : where ω 0 is the angular frequency of the light and E r , E z , E θ (B r , B z , B θ ) are the mode electric field (magnetic flux density) components (the coordinate system in use is shown in Fig. 6). The linear propagation of a field (with fixed ω 0 ) is described by dE 1 =dθ ¼ imE 1 , where E 1 is the field amplitude at different angular positions. If the mode profile in another waveguide with a slightly different shape is nearly identical to the current waveguide, which is usually true up to the first-order of the geometry differences, then the same integral can be used to calculate the propagation constant using the known field profile and the perturbed refractive index profile. For a pair of nearly degenerate modes, the propagation constant generalizes into a matrix: