Modeling of dual frequency combs and bistable solitons in third-harmonic generation

Phase-matching of the third-harmonic generation process can be used to extend the emission of radiation from Kerr microresonators into new spectral regions far from the pump wavelength. Here, we present a theoretical mean-field model for optical frequency combs in a dissipative and nonlinear χ(3)-based cavity system with parametric coupling between fundamental and third-harmonic waves. We investigate temporally dispersive dual-comb generation of phase-matched combs with broad bandwidth and anomalous dispersion of the fundamental field, individuating conditions for accessing a multistable regime that simultaneously supports two types of coupled bright cavity solitons. These bistable cavity solitons coexist for the same pump power and frequency detuning, while featuring dissimilar amplitudes of their individual field components. Third-harmonic generation frequency combs grant telecom pump laser sources a simultaneous and direct access to both the near-infrared and the visible regions, which may prove advantageous for the development of optical clocks and sensing applications. Third-harmonic generation frequency combs grant telecom pump laser sources the direct and simultaneous access to both the near infrared and the visible spectral regions. The authors model the broadband and temporally dispersive dual-comb generation, and identify conditions for accessing a regime supporting two distinct and coexisting cavity solitons.


INTRODUCTION
Optical frequency combs (OFCs) utilize the nonlinear polarization response of a cavity-enclosed dielectric medium, in order to convert an externally applied pump field to multiple new frequencies.Acting as broadband and coherent optical sources, OFCs are a key technology for enabling a diverse range of applications such as frequency metrology, optical communications and spectroscopy [1][2][3].However, conventional Kerr comb synthesizers only emit radiation in a spectral range that is centered around the pump laser frequency, and generally require anomalous group-velocity dispersion for the experimentally accessible formation of phase-locked states [4,5].This makes it challenging to form combs in wavelength ranges which lack suitable pump laser sources, and in spectral regions that exhibit an effective waveguide and material dependent normal dispersion.
One way of overcoming these limitations is to exploit the third-harmonic generation (THG) process of the χ (3) nonlinearity, in order to couple pump field excitations to parametric waves at thrice the fundamental frequency (FF), 3ω 1 .The THG process is inherent in all transparent nonlinear media that display a strong Kerr effect, but in practice it is hampered by the requirement of maintaining a fixed phase relation, which is necessary for efficient frequency conversion [6].While many experimental observations of THG in microcomb devices have relied on refractive index matching between the fundamental and higher-order whispering gallery modes, it is also feasible to accomplish phase-matching through birefringence, periodic poling and other quasi-phase-matching techniques [7].
In this work we consider a centrosymmetric nonlin-ear Kerr resonator system that is engineered to phasematch the third-harmonic process in order to enable resonant dual-comb generation around both the FF and the third-harmonic (TH) frequency, when the dissipative cavity is driven by a continuous-wave (CW) pump source at the fundamental frequency ω 1 .We go beyond previous studies of cavity THG that have been restricted to the non-dispersive case with only two interacting frequencies [8,9], by considering the mutual coupling between sidebands around each carrier wave and the simultaneous interaction of all frequency modes.This system shares similarities with non-phase-matched Kerr microresonators, that can be modelled by the Lugiato-Lefever equation (LLE) [10,11] or driven-and-damped nonlinear Schrödinger equation [12]; it is also analogous to OFCs in quadratically nonlinear resonators, which exploit cascaded processes, found in χ (2) -nonlinear media without inversion symmetry, in order to enable coupling and dual-comb generation around both the FF and the second-harmonic frequency [13,14].We note that a similar model of THG-assisted four-wave mixing was published in Ref. [15] during the final preparation of this manuscript, but with the inclusion of simplifying assumptions that limit its applicability to a perturbative THG regime, and exclude the possibility of generating bistable cavity solitons.
Previous experimental work has demonstrated the direct emission of visible light by THG from an infrared pump using both high-Q whispering-gallery-mode and integrated microresonators [16][17][18][19][20].The generation of OFCs by THG acting together with Raman-assisted spectral broadening in silica based microcavities was also reported [21].Theoretical studies of spatially diffractive beam propagation in conservative, cavityless systems have also shown the possibility of generating both bright and dark coupled solitary wave structures in the presence of THG; that exhibit properties such as a power threshold and bistability [22,23].Given that the dynamics of our system is governed solely by the fundamental field in the limit of vanishing parametric coupling, one may expect to find a similar homotopic extension of the phaselocked bright cavity solitons (CSs) that are supported by the LLE in the case of anomalous group-velocity dispersion [24,25].Surprisingly, we find that coupled FF and TH combs can also support an additional type of CS with a partially overlapping range of existence.These dual, two component, solitons share a common trapping refractive index potential through self-and cross-phase modulation (SPM/XPM).Moreover, there is the intriguing possibility that mutual XPM coupling can be used to overcome the group-velocity mismatch, in order to create various types of synchronized dual-frequency combs with locked repetition rates around both the FF and the TH frequency.
In the following sections, we develop a theoretical mean-field model for a doubly-resonant, cavity-enhanced, dispersive and nonlinear system, phase-matched for THG.In particular we show that, as the TH field grows larger, it may not be simply considered as an upconverted replica of the fundamental comb, but it will reciprocally influence the latter through both parametric coupling and XPM.We investigate the multistability properties of the homogeneous solution and consider the importance of modulational instability (MI) in generating various types of multi-frequency combs.Additionally, we make a detailed numerical study of the multistable regime, where we demonstrate the occurrence of bistable cavity solitons that coexist, when both the FF and the TH frequency lie in the anomalous dispersion regime.

Model
We consider OFC generation in a dispersive χ (3) -based resonator system with coupling between fundamental and third-harmonic fields, as schematically illustrated in Fig. 1.The system is assumed to be resonant around both the driving frequency of the fundamental field (FF) ω 1 as well as the frequency of the third-harmonic (TH) field ω 2 = 3ω 1 .For simplicity, we assume an isotropic nonlinear polarization response PNL = 0 χ (3) Ē3 and a linearly polarized electric field propagating along the zaxis, viz.
where 0 is the vacuum permittivity, χ (3) is the thirdorder nonlinear susceptibility, ê is a unit vector in the 1. Schematic of the THG resonator system.The χ (3) -based nonlinear microresonator is phase-matched for third-harmonic generation.The resonator is driven by a CW pump field at the fundamental frequency ω1 and generates simultaneous frequency combs around both ω1 and ω2 = 3ω1.
polarization direction and c.c. denotes complex conjugate.The OFC generation dynamics is modelled by means of a scalar Ikeda map [26] for the evolution of the temporal field during each roundtrip, together with appropriate boundary conditions for the injection of the external pump and the coupling of the fields from one roundtrip to the next.Starting from Maxwell's equations and applying the slowly-varying envelope approximation, the envelopes of the co-propagating electric field of the fundamental A m and the third-harmonic B m at the mth roundtrip are found to obey the following coupled nonlinear equations: where z is the longitudinal coordinate and t is time.
The dispersive properties of the medium that are associated with the non-equidistant resonance spacing are described by the Taylor series expansions Here k 1,2 are inverse group-velocities, k 1,2 are group-velocity dispersion coefficients and ∆k = 3k 1 −k 2 is a wave-vector mismatch.Moreover, we have that Q ij are modal overlap integrals, α c1,2 are the FF/TH absorption losses, c is the speed of light in vacuum and n 2 (ω) = 3χ (3) (ω)/8n(ω) is the nonlinear coefficient, with n(ω) the linear refractive index.In the case of natural phase-matching we have ∆k = 0 which requires matching of the FF/TH refractive indices n(ω 1 ) = n(ω 2 ) and implies that ∆ 2 = 3∆ 1 , which is assumed in the following.The transverse overlap integrals that are needed to account for the variation in spatial mode profiles between families of different mode orders are given by where Q 21 = Q 12 , Q * 23 = Q 13 and dS = dxdy.These definitions reduce to the familiar Kerr coefficient γ = 11 in the absence of any TH field.
The fields at the beginning of the (m + 1)th roundtrip are assumed to be related to the fields at the end of the mth roundtrip through the boundary conditions that model a generic optical coupling, such as the evanescent field overlap from a nearby waveguide or tapered fiber, that partially transmits the external pump field A in while reflecting the intracavity fields from the previous roundtrip.Here, L is the length of the cavity, while θ 1,2 are the power transmission coefficients and δ 1,2 are the phase detunings of the FF/TH fields from the nearest cavity resonance.We note that the complementary case of a down-converting, 3ω 1 → ω 1 , optical parametric oscillator can be modelled by simply moving the pump term to Eq. ( 6) for the TH field.The Ikeda map constitutes a complete model for the temporal and spectral dynamics of a THG cavity based OFC generation system for general resonance and phasematching conditions.But analytical and numerical investigations can be simplified in the doubly-resonant case (θ 1,2 ≈ 1) by averaging the above map over the roundtrip length into a mean-field model, similar to the LLE.In the following, we truncate the dispersion to the second-order; assume the phase-matching to be almost perfect, so that the coherence length is longer than the cavity length; and shift to a retarded reference frame moving with the group-velocity of the driving field (k 1 ) −1 .Following a derivation, whose details are presented in the Methods section, we obtain our main system of normalized meanfield evolution equations for the FF and TH fields A and B as where t and τ are normalized slow-and fast-time variables, respectively.α = α 2 /α 1 is the ratio of the roundtrip losses, ∆ j = δ j /α 1 is the normalized detuning, d = 2L/(|k 1 |α 1 )∆k is the walk-off parameter that depends on the group-velocity mismatch ∆k = is the ratio of the group-velocity dispersion coefficients and f = θ 1 ω 1 n 2 (ω 1 )LQ 11 /(cα 3 1 )A in is the normalized input pump field.The nonlinear interaction is governed by the four dimensionless parameters that can be assumed to be close to unity, unless the phase-matching is significantly multimodal.It is interesting to note that Eqs.(7)(8) are formally similar to models describing phase-matched doubly-resonant secondharmonic generation (SHG) in quadratic nonlinear media with a simultaneous Kerr nonlinearity [14,27].The two systems differ mainly in the magnitude of the terms and in the exponents of the FF that appears in the parametric coupling terms: these read as BA * and A 2 in the case of second-harmonic generation.

Homogeneous solutions
The response of the system for pump powers below the threshold for parametric comb generation is characterized by CW emission at both FF and TH frequencies.We find a set of stationary mixed-mode homogeneous solutions by setting the derivatives in Eqs.(7)(8) to zero.Eliminating terms that are linear in B one finds that where we have defined the power dependent functions Eqs. (10)(11) are written in a resonance form, where a maximum occurs for a detuning that makes the imaginary part zero.The stationary fields are related through the power conservation law f (A + A * ) = 2|A| 2 + 2(α/ρ)|B| 2 , and the two FF/TH intracavity powers |A| 2 and |B| 2 are seen to satisfy a closed set of real nonlinear equations, viz.
where as before ∆ 2 = 3∆ 1 in the case of natural phasematching.A detailed bistability analysis is complicated, owing to the high-order of the coupled Eqs.(13)(14).However, the equations can be solved numerically, in order to determine their number of solutions, as shown in Fig. 2.
Here, we observe the presence of multistability, with an odd number of simultaneous solutions (1, 3, 5, or 7) in different ranges.We find no separate bistability of the TH: only a single TH solution corresponds to each value of the FF.In fact, because of pump depletion, as well as self-and cross-phase modulation, the solution for the TH field is not simply proportional to A 3 at high powers, but it can be expressed as an explicit function of the FF through the relation In Fig. 3 we show an example that illustrates the TH modification of the homogeneous solution as a function of the FF detuning ∆ 1 for f = 3.As can be seen, the FF exhibits a Kerr tilted resonance shape, similar to that of the LLE model, but with a secondary peak associated with the resonance of the TH at higher values of the detuning.In particular, for a detuning around ∆ 1 = 6 we observe a multistable range with 5 different solutions, with three separate branches that are found to be stable to homogeneous perturbations (solid curves in Fig. 3).

Modulational instability analysis
The stability of the homogeneous solutions against periodic perturbations is important for determining the onset of comb generation and the accessibility of different solution branches.A positive MI gain causes the spontaneous growth of signal and idler sidebands that seed a cascade of phase-dependent four-wave mixing processes, eventually leading to broadband comb formation through the interplay of nonlinearity and dispersion.We analyze the MI gain by linearizing Eqs. where ) and ᾱ = α + idΩ, and the off-diagonal coupling elements are given by p 2 , p 5 = 3ρp 3 and p 6 = i3ρµB 2 0 .The characteristic equation for the eigenvalues λ is somewhat unwieldy, but it can be written in a factorized form that allows for an explicit solution in certain limiting cases, see Methods.It is clear that the MI is independent of the absolute phase through the invariance A 0 → A 0 e iθ0 and B 0 → B 0 e i3θ0 .However, the MI does depend on the relative phase of the fields, because of their coherent coupling.
The presence of TH coupling together with SPM/XPM also causes the appearance of new complex modulational instabilities, that have no counterpart in the LLE model, c.f. Ref. [15].To illustrate this, we calculated the MI growth rates for the two homogeneously stable branches with the highest power in the multistable solution of Fig. 3. Figure 4 shows the dependence of the maximum MI growth rate, max{Re[λ(Ω)]}, that is attainable for a finite perturbation frequency Ω, versus the normalized group-velocity dispersion η 2 and group-velocity mismatch d of the TH.It is seen that both branches generally are unstable to periodic perturbations when the FF dispersion is anomalous, except for the middle branch that has a small stability window for normal dispersion around η 2 = 0.2 and d = 0. We therefore fulfill the conditions for spontaneous sideband amplification that are a prerequisite for comb formation, possibly leading to separate phase-locked patterns and solitons that are associated with each branch.We note that the upper branch is dominated by the FF and only has a weak dependence on the TH, while the FF/TH fields are comparable on the middle branch, c.f. Fig. 3.This is reflected in the stability properties and the sensitivity to the sign of the TH dispersion in the latter case.

Bistable solitons
Next, let us investigate the fascinating possibility of finding bistable cavity soliton solutions in the multistable regime.Soliton bistability refers to the possibility of generating two (or more) localized structures with different temporal and spectral profiles, that coexist for the same values of pump power and cavity detuning.Such bistable (super) cavity solitons have previously been predicted in the framework of an Ikeda map for a Kerr medium without parametric coupling [28], and have subsequently been experimentally confirmed to exist in fiber-ring resonators [29].There it was established that bistable solitons form in regions of multistability, owing to the overlapping tilt of adjacent resonances at high pump powers.The coexistence of bistable vector solitons with different polarization states has also been found to occur in birefringent fiber-ring resonators [30,31].In general, we expect to find CSs in the vicinity of a detuning range where the homogeneous solution is bistable.Bright CSs are typically found in media with anomalous dispersion, and sit on a finite background that constitutes the lowest branch of the homogeneous solution.This background should be stable, while the upper branch may be modulationally unstable in favor of switching the stability to a periodic orbit corresponding to a stationary Turing pattern.A CS can then form by the locking of fronts that connect the homogeneous background and a cycle of the patterned state [32].
We made a numerical search for CSs by performing a series of detuning sweeps over the resonance shown in Fig. 3, for different values and signs of the TH dispersion parameter.We assumed d = 0, since walk-off leads to a separation of the pulse components, and is generally detrimental to soliton formation.The corresponding results are shown in Fig. 5(a): here we plot the total comb power as a function of the FF cavity detuning ∆ 1 and the TH group-velocity dispersion η 2 .We observe a dynamical sequence of evolving comb states with stable and chaotic MI regions followed by a series of steps that are characteristic of multi-soliton states [25].The length of the steps is found to vary with the sign and magnitude of the dispersion parameter, and is seen to be significantly longer in the case of anomalous dispersion (η 2 < 0).The termination point of the steps shows that only the upper branch of the homogeneous solution permits the formation of solitons for normal dispersion with η 2 > 0, while a sharp secondary step is observed at a detuning ∆ 1 ≈ 8.5 for negative TH dispersion with η 2 < −1.This secondary step extends to a detuning ∆ 1 ≈ 12.5, which lies beyond the endpoint of the middle branch in Fig. 3.In Fig. 5(b) we show an example of the temporal evolution of the FF intracavity power for the case η 2 = −2, when both the FF and TH fields experience anomalous group-velocity dispersion.Here, we can identify a broad multi-soliton region that is split between the detuning ranges ∆ 1 ≈ 5 − 8.5 and ∆ 1 ≈ 8.5 − 12.5; and where in the first part we find traces of coexisting localized structures with two different amplitudes, i.e., bistable solitons.Isolated solitons corresponding to a detuning ∆ 1 = 6 are shown in Fig. 6.The two localized solitons may be associated with fronts that connect the homogeneous back-ground with stationary states on the upper or middle branch, respectively.It is seen that the relative amplitudes of the FF/TH soliton components correspond to roughly twice the power of the homogeneous solutions.We observe that the power of the cavity soliton pair in panel (a) is dominated by the FF, with the TH being almost an order of magnitude smaller.Whereas the soliton pair in panel (b) has a much larger TH contribution: the total comb power is more equally split between the FF and TH components.It is notable that the background solution for the TH is very small in both cases, and has a spectral line magnitude that is comparable to that of the sidebands.
To characterize the bifurcation structure undergone by the bistable solitons, we apply a path-continuation method based on a Newton-Raphson solver [33].The results of these computations are depicted in Fig. 7.  6, respectively.In both panels, solid (dashed) lines correspond to stable (unstable) states.The bistability region between CS1 and CS2 is marked using a shadowed box.The different symbols correspond to the temporal profiles shown below.FIG. 8. Characterization of soliton width and peak power.Individual variation of root-mean-square (RMS) width (dashed lines, left axis) and peak power (fully drawn lines, right axis) for the soliton components on the CS1 branch (a) and branch (b).Blue and green lines correspond to the FF while red and purple colors correspond to the TH.Only the stable part of the detuning interval is shown.
, as a function of the detuning ∆ 1 for f = 3.In Fig. 7(b) we show the same result after removing the homogeneous background field, which allows to better illustrate the organization of the solution branches.To do so we plot ||∂ τ A|| 2 against ∆ 1 .The bifurcation diagram associated with CS2 (in blue) originates from the saddle-node bifurcation of the CW solution SN h .In contrast, the one corresponding to CS1 (in orange) forms an isola (i.e., a loop).Both curves exhibit a number of folding branches, with sets of both stable (fully drawn) and unstable (dashed) solutions.Some examples of temporal intensity profiles found along these curves are shown in the bottom row of Fig. 7, and correspond to the locations marked in Fig. 7(b).The stable soliton branches have turning points at the detuning values that correspond to the endpoints of the soliton steps, and the CS1 branch undergoes a Hopf-instability at the point Hopf CS1 .The detuning range where CS1 and CS2 coexist (i.e., the soliton bistability region) is shown using a shadowed box.We further characterize the bistable solitons in Fig. 8, by plotting the variation of the root-mean-square (RMS) width of the temporal duration (dashed) and the peak power (fully drawn), of the soliton components for the fundamental and TH fields along the stable part of each branch.While we can observe a large difference in the peak power of the CS1 soliton components in Fig. 8(a), we see that their RMS time widths remain nearly the same.TH power of the CS2 soliton in Fig. 8(b) makes it preferable as an operating state, in order to achieve optimal conversion efficiency for the TH comb.We observe similar trends of increasing peak power and simultaneously decreasing RMS duration, showing that the total power increases and the solitons become more energetic as the FF detuning grows larger.The soliton pair in panel (a) clearly has a smaller existence range, with the leftmost endpoint corresponding to the Hopf CS1 bifurcation, where the stable soliton transitions into a periodic breather state.
We emphasize that the observed bistability of dissipative CSs is due to the existence of two separate soliton attractors, and is different from the previously reported bistability mechanism of conservative spatial solitons in cavityless THG, that requires the presence of a phase-mismatch [23].A similar bistability of vectorial dissipative solitons has recently been shown to occur for the physically distinct situation of nonlinear polarization mode coupling [30,31].Bistability of conservative soli-tons the same power can also occur due to the appearance of two propagation constants when the nonlinear polarization has a certain functional dependence on the intensity, see Ref. [34].
Admittedly, it may be challenging to find suitable nonlinear cavities where the simultaneous assumptions of zero walk-off and anomalous group-velocity dispersion for both FF and TH fields that were used to find the bistable CSs hold.Nevertheless, it has been previously demonstrated that both criteria can be satisfied under realistic experimental conditions.The temporal walkoff can, e.g., be made to vanish by considering different mode families, or sets of suitably chosen FF/TH frequencies on opposite sides of the zero-dispersion wavelength [35]; whereas anomalous TH dispersion may be obtained through dispersion engineering, or by exploiting avoided mode-crossings [36,37].It is also likely that there may exist a similar bistability of dark dissipative solitons in the same parameter regime when the FF and TH fields both experience normal dispersion.However, the investigation of this case is beyond the scope of this work, and will be the subject of future studies.
Finally, we note that coupled dual solitary wave structures can also be found in the more common situation of normal TH group-velocity dispersion, η 2 > 0.Here the coupling is mainly perturbative, and we observe only short soliton steps that correspond to the upper branch of the homogeneous solution, see Fig. 5(a).These solitons have a weak TH with a broad profile, as shown in Fig. 9(a).They are similar to previously reported quadratic cavity solitons [35], and may persist also in the presence of a large walk-off.The mixed-dispersion condition is clearly non-ideal, but it can be compensated by adjusting the TH detuning ∆ 2 : specifically, supposing that ∆ 2 is not required to satisfy the condition of natural phase-matching (∆ 2 = 3∆ 1 ), but instead can be individually adjusted in order to change its sign.It is then possible to find coupled bright soliton pairs analogous to Fig. 6(a) also for normal TH dispersion, see for example the case of Fig. 9(b).

CONCLUSIONS
In conclusion, we have presented a theoretical model for optical Kerr frequency combs in a doubly-resonant and dispersive cavity system that is phase-matched for third-harmonic generation.We have reported conditions for achieving simultaneous dual-comb generation, and investigated a multistable regime that supports two types of bistable cavity solitons for anomalous dispersion.The parametric coupling between fundamental and thirdharmonic waves allows the formation of simultaneous combs around multiple wavelengths, and is expected to be important for future applications of frequency combs in the visible and ultraviolet spectral range.

Derivation of the mean-field model
To derive the mean-field model, we start by expanding the intracavity fields in a power series as m where is a small parameter, c.f. Ref. [38].The fields remain unchanged to the lowest order and the solution of Eqs.
m (0).Inserting these solutions on the right hand side of the first-order equations and assuming that all terms are small, we can immediately carry out an integration to find that A (1)  m (L) =A (1)  m (0) + − where ∆k = k 2 − k 1 is the group-velocity mismatch and κ = e i∆kL/2 sinc(∆kL/2).The boundary condition Eq. (5-6) similarly become A (0) m (L) to the lowest order, while the first order relations are given by By combining the above expressions, and introducing a slow time variable t we obtain a continuation of the map by setting A m+1 (0) − A m (0) → t R ∂A/∂t and B m+1 (0) − B m (0) → t R ∂B/∂t where t R is the roundtrip time.This leads to the following coupled system of meanfield evolution equations for the fundamental and thirdharmonic fields where α j = (α cj L + θ j )/2 is the total roundtrip loss.We note that, besides the form of the nonlinearity, the resulting mean-field equations are analogous to models which have previously been obtained for doubly-resonant SHG frequency combs in χ (2) resonators, see Ref. [39] for a derivation using an alternative approach.
The above equations are normalized with respect to the nonlinear coefficient and the timescales for the losses and dispersion of the FF by rescaling A → ω 1 n 2 (ω 1 )LQ 11 /(cα 1 )A, B → ω 1 n 2 (ω 1 )LQ 11 /(cα 1 )B, t → (α 1 /t R )t and τ → 2α 1 /(|k 1 |L)τ .This choice of normalization results in the mean-field Eqs.(7-8) of the main text.We note that the standard normalization for the LLE is obtained in the case that B = 0, which permits easy comparisons with previous Kerr comb results.

Numerical methods
The mean-field Eqs.(7-8) have been numerically solved using a split-step Fourier method.The two equations are solved simultaneously by integrating them in a series of short alternating linear and nonlinear propagation steps with different basis.The dispersive linear step is performed in the frequency domain with N = 2048 modes and the forward and transforms are implemented using Fast-Fourier Transforms (FFTs).The nonlinear step is meanwhile performed in the time-domain using a 4th-order Runge-Kutta algorithm.
The soliton bifurcation diagrams shown in Fig. 7 have been obtained through a path-continuation algorithm, by computing the stationary solutions of Eqs. and This allows us to compute soliton states as a boundary value problem, imposing Neumann boundary conditions at 0 and t R /2, and utilizing the open distribution software package AUTO-07p [40].The linear stability of these states is ascertained through computation of the eigenvalues of the Jacobian matrix associated with the system (28).

FIG. 2 .
FIG. 2. Multistability of homogeneous solutions.The plot shows colored parameter regions with 1,3,5 or 7 simultaneous homogeneous solutions that coexist for the same normalized detuning ∆1 and pump power f .

FIG. 3 .
FIG.3.Multistable resonance response.The plots show intracavity power of the homogeneous solution as a function of detuning when the cavity is driven by an external pump with normalized power f = 3. Panels (a) and (b) show the intracavity power of the FF and TH fields, respectively, where blue and red colors denote branches that are stable/unstable to homogeneous perturbations.The dotted vertical line at ∆1 = 6 indicates the location of the bistable cavity solitons considered below.

FIG. 4 .
FIG. 4. Instability diagrams for the homogeneous solution.Maximum MI growth rates for the multistable homogeneous solution in Fig. 3 as a function of TH dispersion η2 and frequency Ω, panels (a),(b), and walk-off d, panels (c),(d).The panels (a),(c) and (b),(d) show MI growth rates for the CW stable upper and middle branches, respectively, while the lowest branch is unconditionally stable (not shown).Case of ∆1 = 6 and anomalous dispersion for the FF, η1 = −1.

FIG. 5 .FIG. 6 .
FIG. 5. Numerical simulation of detuning sweeps.(a) Variation of comb power for different values of TH group-velocity dispersion η2.The dispersion of the FF is assumed to be anomalous, η1 = −1.(b) Temporal evolution of the FF power for the case of η2 = −2.Characteristic signatures of a multisoliton regime with two different types of bistable cavity solitons indicated by CS1 and CS2 are seen for ∆1 > 5.

FIG. 7 .
FIG. 7. Bifurcation structure of bistable CS solutions.(a) Bifurcation diagram showing the modification of ||A|| 2 against ∆1 for f = 3.(b) shows the same bifurcation curves as (a) but plotting ||∂τ A|| 2 vs. ∆1.Orange and blue curves are associated with the solitons CS1 and CS2 of Fig.6, respectively.In both panels, solid (dashed) lines correspond to stable (unstable) states.The bistability region between CS1 and CS2 is marked using a shadowed box.The different symbols correspond to the temporal profiles shown below.