Introduction

Commonly described as "Rulers for Light", optical frequency combs1,2,3,4 and their generation are topics with wide-reaching research interest. This interest is especially owed to their diverse application in metrology, optical communications5, and novel photonic devices3. They are used, for example, in high precision spectroscopy, optical atomic clocks, and navigation systems. One particular method for producing an optical frequency comb utilizes dissipative temporal cavity solitons (TCS) in optical micro-resonators.

We concern ourselves here not only with the generation of TCS but also with combining them with the phenomenon known as spontaneous symmetry breaking (SSB). Broadly defined, SSB describes the situation when two or more properties of a system suddenly change from displaying an equal property or state (symmetric) to having these states becoming unequal (asymmetric) following an infinitely small change to some control parameter.

The SSB of light in Kerr resonators has seen a flurry of study and interest in the past half-decade, being predicted and observed for systems with counter-propagating light components6,7,8,9,10,11,12,13,14,15, for orthogonally polarized light components16,17,18,19,20,21,22,23, and also very recently for a single system combining both counter-propagating light and two orthogonal polarization components together24, and systems containing multiple resonators25,26.

Although recent interest about TCS in Kerr resonators mainly focused on ring geometries, frequency combs based on TCS have also been found in Fabry-Pérot resonators (linear cavities)27,28,29. Figure 1 shows a basic schematic of this type of resonator, and this is the configuration of interest in this paper. One can see that it is comprised primarily of two highly reflective mirrors, which bounce light from the linearly polarized input light back and forth between them. Here we consider that the space between the mirrors is filled with a Kerr, or χ3, nonlinear medium.

Fig. 1: Schematic of a Kerr Fabry-Pérot cavity.
figure 1

The cavity is bounded by two highly reflective mirrors, with a Kerr nonlinear medium filling the space between them. Linearly polarized laser input is introduced via mirror 1, allowing the light to enter the cavity. Once inside, the light undergoes many reflections between the mirrors before eventually exiting. The intracavity field is analyzed by decomposing it into its left- and right-circularly polarized components, E±.

Full size image

Vector solitons, opposed to scalar solitons, refer to solitons involving two or more field components30. There has been great success in observing the SSB of a pair of vector TCS in Kerr ring resonators18,19, where the TCS have orthogonal polarizations, but similar phenomena in linear Fabry-Pérot cavities, Fig. 1, has remained unreported. This is despite recent separate experimental observations in Fabry-Pérot cavities of both the SSB of flat solutions16 and of scalar TCS28. Here we outline the SSB of a pair of vector TCS in Fabry-Pérot cavities.

Results and discussions

To model the intra-cavity, and slower, temporal dynamics of a field propagating in a Fabry-Pérot resonator with consideration for its polarization, Fig. 1, we derive, see the Methods section, a system of two coupled, generalized, Lugiato-Lefever Equations (LLEs)31 with fast-time averaged terms. Here we rewrite the final integro-partial differential equations from the Methods section for the complex amplitudes E± of the circularly polarized components

$$\frac{\partial {E}_{\pm }}{\partial t}=\, {E}_{{{{{{{{\rm{in}}}}}}}}}-{E}_{\pm }-{{{{{{{\rm{i}}}}}}}}\theta {E}_{\pm }-{{{{{{{\rm{i}}}}}}}}\eta \frac{{\partial }^{2}{E}_{\pm }}{\partial {\tau }^{2}}\\ +{{{{{{{\rm{i}}}}}}}}\left(A| {E}_{\pm }{| }^{2}{E}_{\pm }+B| {E}_{\mp }{| }^{2}{E}_{\pm }+2A\langle | {E}_{\pm }{| }^{2}\rangle {E}_{\pm }+B\langle | {E}_{\mp }{| }^{2}\rangle {E}_{\pm }+B\langle {E}_{\pm }{E}_{\mp }^{* }\rangle {E}_{\mp }\right),$$
(1)

where Ein is the input pump in the circular basis32, θ the cavity detuning, η controls the type of dispersion (with its value ± 1 referring to normal or anomalous dispersion, respectively), t and τ are both temporal variables but on relative slow and fast time-scales, respectively33, and A and B control the strengths of self- and cross-phase modulation effects11,24. The terms with angled brackets, i.e., 〈f〉, represent the temporal averages of the encapsulated functions over a single resonator round-trip, and are defined by the following:

$$\langle f\rangle =\frac{1}{{t}_{{{{{{{{\rm{R}}}}}}}}}}\int\nolimits_{-{t}_{{{{{{{{\rm{R}}}}}}}}}/2}^{{t}_{{{{{{{{\rm{R}}}}}}}}}/2}f\left(\tau \right)\,d\tau .$$
(2)

Comparing it with the equations of Cole et al.27,29, we see that the effects of the polarization components cause not only additional cross-phase modulation terms (and their averages) but also a further additional term, the final term of Eq. (1), which is caused by an energy exchange between the two circular components of each beam24,34. We note that, if so desired, one may transition Eq. (1) to the linear polarization basis through its manipulation around the substitution of \({E}_{\pm }=\frac{1}{\sqrt{2}}\left({E}_{x}\pm i{E}_{y}\right)\). Note also that the main symmetry of Eq. (1) is the invariance under the exchange of the plus and minus indexes. We recognize the potential for further generalization and study of this system when taking into account higher order dispersion effects35,36,37,38,39, which could result in additional localized structures to those discussed here.

Homogeneous stationary states

The homogeneous (meaning here unchanging over the τ domain) and stationary (meaning here unchanging with t) solutions (HSS) to Eq. (1) are obtained by setting ∂E±/∂t and ∂2E±/∂τ2 both to zero. We may further make use of the fact that when a function f(τ) is homogeneous over τ its average over the τ domain is equal to its value at a single point—that is to say Eq. (2) becomes trivially f. Hence, under homogeneous and stationary conditions, and following suitable algebraic manipulation, Eq. (1) becomes

$$| {E}_{\pm }{| }^{2}=\frac{| {E}_{{{{{{{{\rm{in}}}}}}}}}{| }^{2}}{1+{\left(3A| {E}_{\pm }{| }^{2}+3B| {E}_{\mp }{| }^{2}-\theta \right)}^{2}}.$$
(3)

Equation (3) is identical in its mathematical form, although not physical meaning, to the homogeneous stationary states of two linearly polarized fields counter-propagating a Kerr ring resonator6,7,10,40, or of two orthogonally polarized fields co-propagating a Kerr ring resonator11,17,32. Due to its mathematical analogies to other such systems, Eq. (3) has been studied extensively, and we shall not repeat that analysis here.

We stress that SSB emerges as a property of Eq. (3). In this context, SSB describes the phenomena where the two circulating fields E± go from having symmetric intensities, ∣E+∣2 = ∣E−∣2, to having asymmetric intensities, ∣E+∣2 ≠ ∣E−∣2, upon an infinitely small change to input conditions (input conditions such as the input intensity ∣Ein∣2 or the cavity detuning, θ, for example). With Eq. (3) and its SSB largely explored elsewhere, we display, for the later benefit of the reader, in Fig. 2(a) & (b), examples of SSB in input intensity and cavity detuning scans for Eq. (3), respectively, for self- and cross-phase modulation values A = 2/3 and B = 4/3 – the values for silica glass fibers. These values give a very general B/A ratio11.

Fig. 2: Spontaneous symmetry breaking of homogeneous solutions.
figure 2

Variations in the circulating fields' intensities, ∣E±∣2 as (a) input intensity, and (b) cavity detuning, are scanned. For this figure we use A = 2/3 and B = 4/3, with θ = 2 for (a) and ∣Ein∣2 = 0.9 for (b). Observe how in both panels "red'', asymmetric, solutions emerge spontaneously from "black", symmetric, solutions at points where SSB occur via pitchfork bifurcations. Note regions where the symmetric solution line is unstable to perturbations in t (dashed line) on the middle branches of optical bistability, and also between SSB points.

Full size image

In the next section, we proceed to describe the stability of the homogeneous and stationary states and hence their susceptibility to oscillations following fast- and/or slow-time perturbations.

Linear stability analysis of homogeneous stationary states

To assess the system’s susceptibility to temporal perturbations, both on t and τ, we performed a linear stability analysis27 on the modal expansion of Eq. (1), see methods section. General mathematical methods to find instability thresholds in single polarization Fabry-Pérot resonators have recently been established41. Here, we linearized the modal equations around a homogeneous stationary solution E±HS with Un = E±HSδn,0 + δUn, which, without loss of generality, had its phases adjusted such that it is real. The index n is the modal number of the fast-time solution. We find that the four linear stability eigenvalues have the form

$$\lambda =-1\pm \frac{\sqrt{-{\alpha }_{+}{\beta }_{+}-{\alpha }_{-}{\beta }_{-}-\frac{1}{6}\left(1-{\delta }_{n,0}\right)C\pm S}}{\sqrt{2}}$$
(4)

with

$$S=\sqrt{{({\alpha }_{+}{\beta }_{+}-{\alpha }_{-}{\beta }_{-})}^{2}+{\alpha }_{+}{\alpha }_{-}C+\left(1-{\delta }_{n,0}\right)\frac{C}{3}\left({\alpha }_{+}{\beta }_{+}+\frac{1}{3}{\beta }_{+}{\beta }_{-}+{\alpha }_{-}{\beta }_{-}\right)},$$
(5)

where

$${\alpha }_{\pm } =\, -\theta +{\eta }_{n}+3A{E}_{\pm {{{{{{{\rm{HS}}}}}}}}}^{2}+B\left(2+{\delta }_{n,0}\right){E}_{\mp {{{{{{{\rm{HS}}}}}}}}}^{2},\\ {\beta }_{\pm } =\, -\theta +{\eta }_{n}+A\left(5+4{\delta }_{n,0}\right){E}_{\pm {{{{{{{\rm{HS}}}}}}}}}^{2}+B\left(2+{\delta }_{n,0}\right){E}_{\mp {{{{{{{\rm{HS}}}}}}}}}^{2},\\ C = \, 36{B}^{2}\left(1+3{\delta }_{n,0}\right){E}_{+{{{{{{{\rm{HS}}}}}}}}}^{2}{E}_{-{{{{{{{\rm{HS}}}}}}}}}^{2},$$
(6)

where

$${\eta }_{n}=\eta {\left(n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\right)}^{2}=\eta {k}^{2}.$$
(7)

In Fig. 3, we display the HSS, Eq. (3), and their corresponding linear stability eigenvalues, Eq. (4) from the Methods section, over a cavity detuning scan of Eq. (3) for anomalous dispersion, k = 1.9, ∣Ein∣2 = 4, and A = 2/3 and B = 4/3.

Fig. 3: Homogeneous stationary state stability analysis.
figure 3

Regions of stability on t and τ, (a), are displayed for a cavity detuning scan of (b) the homogeneous and stationary circulating field intensities (see Eq. (3)), (c) the Real part of the associated stability eigenvalues see Eq. (4)), and (d) the Imaginary part of the associated stability eigenvalues (see Eq. (4)), for anomalous dispersion, k = 1.9, |Ein|2 = 4, and A = 2/3 and B = 4/3. Eigenvalues in (c) and (d) use the color corresponding to the different branches of HSS presented in (b). When at least one real part of the stability eigenvalues associated with a solution is greater than zero, \({{{{{{{\rm{Re}}}}}}}}(\lambda) > 0\), then the solution is unstable to perturbations. If the imaginary part of said eigenvalue is not equal to zero, \({{{{{{{\rm{Im}}}}}}}}(\lambda)\ne 0\) then the solution is susceptible to oscillations following perturbations, with the oscillations occurring over slow-time, t, fast-time, τ, (Turing patterns), or both. a This summarizes the stability of the various colored HSS solutions of (b) based upon the stability eigenvalues of (c) and (d).

Full size image

Regarding solution stability, if the eigenvalues of Eq. (4) all have real parts less than zero, i.e., Re(λ) < 0, then the solution in question is stable, whereas if at least one eigenvalue has a real part greater than zero, then the solution is unstable. If the stability eigenvalues with real parts greater than zero and imaginary part different from zero occur when n = 0, then the solution will begin to oscillate over slow-time, if they occur when n ≠ 0 the solution will oscillate over fast-time (Turing patterns).

We discuss first the eigenvalues associated with the symmetric solution line of this plot. They confirm the instability of the middle branch of the tilted symmetric Lorentzian curve, and also that the upper branch of the same curve is unstable between the SSB bifurcations points as expected from previous work. Note, however, that the instability of the symmetric line extends even beyond the higher detuning value corresponding to the reverse SSB bifurcation owed to the unstable n ≠ 0 eigenvalues, resulting in symmetric Turing patterns. Similarly, the instability of the asymmetric solution line also extends beyond slow-time instabilities, again owed to the unstable n ≠ 0 eigenvalues resulting in asymmetric Turing patterns, since they remain unstable for detuning values larger than those where the n = 0 eigenvalues have stabilized.

Inhomogeneous solutions: patterns, temporal cavity solitons

Due to Eq. (1)’s susceptibility to instabilities on the fast-time, we further explored the inhomogeneous solutions of the system. In this section, we always use A = 2/3 and B = 4/3, and here analyze situations only for anomalous dispersion, η = −1, –situations related to normal dispersion, η = +1, will be discussed elsewhere. Figure 4 shows in panels (a)–(d) the variation in the field intensity profiles across the cavity as the cavity detuning is scanned, for a set input pump of ∣Ein∣2 = 4. We observe eight distinct regions of characteristic behaviors labeled in the figure as regions (1) – (8). The detuning scans begin with the field intensity profiles following the symmetric HSS line (region 1), before forming symmetric Turing patterns across τ (region 2), in agreement with the linear stability analysis. A characteristic example of these symmetric patterns is presented in Fig. 4e, its peaks grow in intensity as the cavity detuning is further increased. While these pattern states may be initially symmetric, this property later breaks spontaneously, resulting in a region of asymmetric patterns (region 3), with characteristics similar to that of Fig. 4f. Progressing further still across the cavity detuning scan, the pattern’s intensity profiles soon become unstable to chaotic oscillations (region 4), see the example in Fig. 4g, before finally reaching a detuning value where the pattern states naturally form vector soliton-pair structures (region 5), see Fig. 4h. One will note that initially the solitons in each of the field intensity profiles ∣E±∣2 are both asymmetric and breathing (region 5). In Fig. 4i, however, these asymmetric soliton pairs stop breathing and stabilize (region 6). The asymmetric soliton pairs eventually converge to symmetric profiles at the point starting region 7, as the example shows in Fig. 4j. Finally, in region 8, the soliton pairs die, and the field intensity profiles return to the HSS line.

Fig. 4: Vector soliton pair generation.
figure 4

Cavity detuning is scanned for an input intensity of ∣Ein∣2 = 4. a shows the maximum value achieved across the cavity by the intensity profiles of the two fields, ∣E±∣2 in red and blue, respectively. b and c respectively show the full field intensity profiles, ∣E±∣2, over the cavity, while (d) shows the difference between them. We note eight distinct regions in (a–d), which we label (1)–(8). These regions are discussed in the main text. e–j We show characteristic intensity profiles reflecting those existing in regions (2)–(7), respectively, with regions (1) and (8) being simply homogeneous stationary profiles.

Full size image

Long-range talking of soliton pairs and their polarization conformity

In this section, we show that the angle-bracketed terms of Eq. (1) have fascinating repercussions on the system’s evolution. Although one of their effects has already been felt, the requirement to analyze separately the n = 0 stability eigenvalues, and their wider effects have remained somewhat hidden until this section. We have stated that, mathematically, these angled terms amount to an average of the field intensity profiles across the cavity, but owed to this they cause a global coupling across the cavity40. Related to the TCS pairs, these global coupling terms effectively mean that all the soliton pairs feel the impact of, and in turn influence, any others in the cavity.

One effect of this global communication, which was noted by Cole et al.27 for the single polarization case, is that the maximum number of solitons existing within the cavity at any given time can be controllably limited by one’s choice of system parameters. This ability could be extremely useful, for example, for applications that require robust temporal delays between subsequent soliton generations, such as for frequency combs. We show that this control on the maximum number of solitons present in the cavity is maintained for our vector soliton pairs.

Normally solitons in LLE systems, of the nature similar to those discussed in this manuscript, “sit” upon the bottom branch of the symmetric HSS, and hence naturally lose viability when the HSS bottom branch ends. However, as displayed in Fig. 5, and owed to the 〈f〉 terms, these soliton bases are in fact elevated above the HSS in our system. One can also see that these 〈f〉 terms further have the effect of raising the cavity detuning value at which the solitons gain viability. This is because the 〈f〉 terms cause an effective detuning in the system27,40, i.e., they actively alter the system’s resonance mismatch and hence the value of the light intensity in the soliton background. Further still, one notes that the more solitons present in the cavity, the more by which this detuning limit is raised. This increase in the cavity detuning value required to support additional solitons effectively limits the maximum soliton pair number supported at any one time, based upon the current cavity detuning value. In relation to Fig. 5, for example, at a cavity detuning of θ = 5.45, the system can only support a maximum of a single pair of vector solitons, with the system being unable to support even a single additional pair until well after θ = 5.5. This effect amounts to a system which can guarantee the generation of a maximum of a single vector soliton pair if so desired. Looking at the intensity profile of a single asymmetric soliton pair in the frequency domain, Fig. 6, we see the types of symmetry-broken frequency combs that can be produced utilizing this phenomenon.

Fig. 5: Existence of a viable soliton base.
figure 5

a Variation in the soliton-base intensity as the cavity detuning is varied, and also as the number of solitons within the cavity varies, for ∣Ein∣2 = 4. The solid black line displays the homogeneous stationary state (HSS) line, the green line tracks the soliton-base when a single soliton pair exists within the cavity, the blue line tracks when two soliton pairs are present, and the red line tracks when three soliton pairs are present. Note that due to the asymmetry between the solitons existing in ∣E±∣2 respectively there are two lines for each color tracking the soliton-base. b we show a cavity containing three symmetry broken vector solitons in the field profiles of ∣E±∣2 respectively (blue and red curves) and the position where we track the soliton-base for our scans in (a). The results of (a) come from numerical simulations with initial conditions containing the n-soliton pairs and tracks the detuning from high to low in order to find the different detuning limits of the n-soliton pair viabilities.

Full size image
Fig. 6: Examples of symmetry-broken TCSs and dual frequency combs.
figure 6

Field intensity profiles, ∣E±∣2, in red and blue, respectively, across the cavity containing a single asymmetric soliton pair, (a), and the corresponding frequency combs in (b). Here ∣Ein∣2 = 4 and θ = 6.2.

Full size image

We now move to the final effect of the angle-bracketed terms that will be discussed in this manuscript. The end result of this effect was already apparent in panels (h) and (i) of Fig. 4. One will note that the soliton pairs in these panels, both those breathing, and those that are stable, all share the same dominant polarization although generated autonomously. This result is particularly intriguing when compared with that for soliton pairs produced in a Kerr ring18,19, as opposed to a FP cavity, where there each asymmetric soliton pair’s dominant polarization was independent of that held by any other pair. This meant that provided those soliton pairs were sufficiently spaced apart from each other, they could behave and be addressed independently. The fact that we always observe all simultaneous asymmetric soliton pairs in our cavities sharing the same dominant polarization is due to the global coupling terms that facilitate long-range interactions between all soliton pairs existing within the cavity, terms not present in the Kerr ring system, in particular the final term of Eq. (1) caused by the energy exchange between the two circular components of each beam.

We find that the conformity of the soliton pairs to a globally dominant polarization is a strong and robust effect. We demonstrate this in Fig. 7. The main panel of the figure shows the slow-time evolution of a tracked point in the fast-time field intensity profiles, a point chosen to line up with the intensity maxima of one of the cavity soliton pairs. In panel (a) one sees the initial stable cavity condition made up of a train of four asymmetric soliton pairs, all sharing the same dominant polarization. This configuration is stable in the slow time. At time t = 25 we attempt to force non-conformity on the system. We do this by splicing the field intensity profiles and swapping the field roles of one of the soliton pairs such that in the swapped pair (third) the opposing polarization is now dominant. Allowing the system to continue to evolve after this attempt to force non-conformity shows how the system evolves back to a state where soliton pair polarization conformity appears once again. This behavior can have advantages in, for example, protecting communicated binary data encoded in the polarization of a chain of soliton pairs in noisy systems. The ability of the system to return to its original polarization dominance at conformity is limited in the way one would likely expect - if a minority of soliton pairs have the polarizations reversed, then the system conforms to the original dominant polarization. If, however, the majority of soliton pairs have their polarization reversed, then the system evolves to the new dominant-on-average polarization until that polarization is conformed to across the system. If exactly half of the soliton pairs have their polarization inverted, then the eventual conformed polarization is determined by random perturbations.

Fig. 7: Long-range “talking” and conforming solitons.
figure 7

For ∣Ein∣2 = 4, θ = 9, (a) of the figure tracks the paired soliton peaks in, ∣E±∣2, in red and blue, respectively, at a set position within the cavity as they evolve over time. The stable initial cavity condition is formed by a train of four pairs of asymmetric solitons and is displayed in the smaller panel (b). At a slow-time value of t = 25, the third soliton pair is forced to switch its polarization profiles, small panel (c). The system then evolves until it becomes stable once again, small panel (d). One should note that all four soliton pairs in (d) have conformed again to share the same dominant polarization.

Full size image

Conclusions

In conclusion, we have demonstrated the spontaneous symmetry breaking of a pair of vector temporal cavity solitons in Fabry-Pérot cavities with possible experimental realizations in fiber loop set-ups similar to those of ref. 18 but in a Fabry-Pérot configuration. Our investigation not only revealed the stability and dynamical behavior of these vector solitons in various parameter regimes but also uncovered useful phenomena arising from the global coupling terms in our derived model. We showed that the maximum number of soliton pairs in the cavity can be self-limited with important implications for frequency comb applications. Furthermore, we discovered that asymmetric soliton pairs in our system exhibit a conformity in their dominant polarization states, highlighting the effects on soliton pair behavior caused by the global coupling terms. These findings significantly expand our understanding of temporal cavity solitons and spontaneous symmetry breaking in Kerr cavities and may have extensions into studies around frequency combs from high order dispersion, thus opening up new avenues for further exploration and potential applications in optical metrology, pulse coding, optical communications and beyond.

Methods

The nonlinear interaction terms in our model are based on those of Pitois et al.34 where counter-propagation of polarization components is considered in an isotropic fiber. In this work, we add dispersion and consider the reflective boundary conditions, the input pump, and the cavity detuning and other losses, all of which are inherent to a Fabry-Pérot resonator27,32. Then we generalize the equations for various self- (A) and cross-phase (B) modulation strengths,11,24, resulting in the following Eqs. (8)–(11)).

Right circular polarization

Forward propagating, u(t, τ):

$$\frac{\partial u}{\partial t}+{v}_{{{{{{{{\rm{g}}}}}}}}}\frac{\partial u}{\partial \tau }=\, {E}_{{{{{{{{\rm{in}}}}}}}}}-u-{{{{{{{\rm{i}}}}}}}}\theta u-{{{{{{{\rm{i}}}}}}}}\eta \frac{{\partial }^{2}u}{\partial {\tau }^{2}}\\ +{{{{{{{\rm{i}}}}}}}}\left(A| u{| }^{2}u+B| v{| }^{2}u+2A| \bar{u}{| }^{2}u+B| \bar{v}{| }^{2}u+B\bar{u}{\bar{v}}^{* }v\right),$$
(8)

Backward propagating, \(\bar{u}(t,\tau )\):

$$\frac{\partial \bar{u}}{\partial t}-{v}_{{{{{{{{\rm{g}}}}}}}}}\frac{\partial \bar{u}}{\partial \tau }=\, {E}_{{{{{{{{\rm{in}}}}}}}}}-\bar{u}-{{{{{{{\rm{i}}}}}}}}\theta \bar{u}-{{{{{{{\rm{i}}}}}}}}\eta \frac{{\partial }^{2}\bar{u}}{\partial {\tau }^{2}}\\ +{{{{{{{\rm{i}}}}}}}}\left(A| \bar{u}{| }^{2}\bar{u}+B| \bar{v}{| }^{2}\bar{u}+2A| u{| }^{2}\bar{u}+B| v{| }^{2}\bar{u}+Bu{v}^{* }\bar{v}\right),$$
(9)

Left circular polarization

Forward propagating, v(t, τ):

$$\frac{\partial v}{\partial t}+{v}_{{{{{{{{\rm{g}}}}}}}}}\frac{\partial v}{\partial \tau }=\, {E}_{{{{{{{{\rm{in}}}}}}}}}-v-{{{{{{{\rm{i}}}}}}}}\theta v-{{{{{{{\rm{i}}}}}}}}\eta \frac{{\partial }^{2}v}{\partial {\tau }^{2}}\\ +{{{{{{{\rm{i}}}}}}}}\left(A| v{| }^{2}v+B| u{| }^{2}v+2A| \bar{v}{| }^{2}v+B| \bar{u}{| }^{2}v+B\bar{v}{\bar{u}}^{* }u\right),$$
(10)

Backward propagating, \(\bar{v}(t,\tau )\):

$$\frac{\partial \bar{v}}{\partial t}-{v}_{{{{{{{{\rm{g}}}}}}}}}\frac{\partial \bar{v}}{\partial \tau }=\, {E}_{{{{{{{{\rm{in}}}}}}}}}-\bar{v}-{{{{{{{\rm{i}}}}}}}}\theta \bar{v}-{{{{{{{\rm{i}}}}}}}}\eta \frac{{\partial }^{2}\bar{v}}{\partial {\tau }^{2}}\\ +{{{{{{{\rm{i}}}}}}}}\left(A| \bar{v}{| }^{2}\bar{v}+B| \bar{u}{| }^{2}\bar{v}+2A| v{| }^{2}\bar{v}+B| u{| }^{2}\bar{v}+Bv{u}^{* }\bar{u}\right),$$
(11)

where u and v are the forward propagating field components with right- and left-circular polarizations, respectively, with \(\bar{u}\) and \(\bar{v}\) being the backwards propagating variants (representing the reflections of fields u and v respectively), Ein accounts for the input pump, θ is the cavity detuning (the difference between the frequency of the input pump laser and the closest cavity resonance frequency), η controls the type of dispersion (with its value ± 1 referring to normal or anomalous dispersion, respectively), and t and τ are both temporal variables but on relative slow and fast time-scales, respectively33. We also set the following boundary conditions

$$\begin{array}{ll}u(t,0)=\bar{u}(t,0),&u(t,{t}_{{{{{{{{\rm{R}}}}}}}}}/2)=\bar{u}(t,{t}_{{{{{{{{\rm{R}}}}}}}}}/2),\\ v(t,0)=\bar{v}(t,0),&v(t,{t}_{{{{{{{{\rm{R}}}}}}}}}/2)=\bar{v}(t,{t}_{{{{{{{{\rm{R}}}}}}}}}/2).\end{array}$$
(12)

where tR is the round trip time of the resonator.

Drawing inspiration from Cole et al.27, we seek to combine Eq. (8) and (9) and Eq. (10) and (11), respectively, together, such that Eqs. ((8)–(11)) reduce to a system with only two equations. To achieve this, we introduce the following modal expansions in terms of the modal amplitudes \({\tilde{u}}_{n},{\tilde{v}}_{n},{\tilde{E}}_{in}\):

$$u =\, \mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{u}}_{n}(t){e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }, \bar{u}=\mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{u}}_{n}(t){e}^{-{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau },\\ v =\, \mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{v}}_{n}(t){e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau },\bar{v}=\mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{v}}_{n}(t){e}^{-{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau },\\ {E}_{{{{{{{{\rm{in}}}}}}}}} =\, \mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{E}}_{{{{{{{{\rm{in}}}}}}}}}{\delta }_{n,0}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau },$$
(13)

where δn,0 is the Kronecker delta function, and the modal amplitudes are given by

$${\tilde{u}}_{n}(t) =\frac{1}{{t}_{{{{{{{{\rm{R}}}}}}}}}}\int\nolimits_{-{t}_{{{{{{{{\rm{R}}}}}}}}}/2}^{{t}_{{{{{{{{\rm{R}}}}}}}}}/2}\,d\tau \,u\left(t,\tau \right){e}^{-{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }=\frac{1}{{t}_{{{{{{{{\rm{R}}}}}}}}}}\int\nolimits_{-{t}_{{{{{{{{\rm{R}}}}}}}}}/2}^{{t}_{{{{{{{{\rm{R}}}}}}}}}/2}\,d\tau \,\bar{u}\left(t,\tau \right){e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau },\\ {\tilde{v}}_{n}(t) =\frac{1}{{t}_{{{{{{{{\rm{R}}}}}}}}}}\int\nolimits_{-{t}_{{{{{{{{\rm{R}}}}}}}}}/2}^{{t}_{{{{{{{{\rm{R}}}}}}}}}/2}\,d\tau \,v\left(t,\tau \right){e}^{-{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }=\frac{1}{{t}_{{{{{{{{\rm{R}}}}}}}}}}\int\nolimits_{-{t}_{{{{{{{{\rm{R}}}}}}}}}/2}^{{t}_{{{{{{{{\rm{R}}}}}}}}}/2}\,d\tau \,\bar{v}\left(t,\tau \right){e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }.$$
(14)

These allow us to extend our field equations over a full round trip since

$$u(t,\tau ) = \bar{u}(t,-\tau ),\bar{u}(t,\tau )=u(t,-\tau )\\ v(t,\tau ) = \bar{v}(t,-\tau ),\bar{v}(t,\tau )=v(t,-\tau ).$$
(15)

Focusing momentarily alone on Eq. (8), we insert the modal expansions of Eqs. (13) and (14) to obtain

$$ \mathop{\sum }\limits_{n=-\infty }^{\infty }\frac{\partial {\tilde{u}}_{n}}{\partial t}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }+{v}_{{{{{{{{\rm{g}}}}}}}}}\mathop{\sum }\limits_{n=-\infty }^{\infty }{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}{\tilde{u}}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }\\ =\mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{E}}_{in}{\delta }_{n,0}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }-\mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{u}}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }\\ -i\theta \mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{u}}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }-{{{{{{{\rm{i}}}}}}}}\eta \mathop{\sum }\limits_{n=-\infty }^{\infty }{\left({{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\right)}^{2}{\tilde{u}}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }\\ +{{{{{{{\rm{i}}}}}}}}\mathop{\sum }\limits_{n=-\infty }^{\infty }\mathop{\sum }\limits_{{n}^{{\prime} }=-\infty }^{\infty }\mathop{\sum }\limits_{{n}^{{\prime\prime} }=-\infty }^{\infty }\left(A{\tilde{u}}_{n+{n}^{{\prime} }-{n}^{{\prime\prime} }}{\tilde{u}}_{{n}^{{\prime} }}^{* }{\tilde{u}}_{{n}^{{\prime\prime} }}+B{\tilde{v}}_{n+{n}^{{\prime} }-{n}^{{\prime\prime} }}{\tilde{v}}_{{n}^{{\prime} }}^{* }{\tilde{u}}_{{n}^{{\prime\prime} }}\right.\\ +\left.2A{\tilde{u}}_{-n+{n}^{{\prime} }+{n}^{{\prime\prime} }}{\tilde{u}}_{{n}^{{\prime} }}^{* }{\tilde{u}}_{{n}^{{\prime\prime} }}+B{\tilde{v}}_{-n+{n}^{{\prime} }+{n}^{{\prime\prime} }}{\tilde{v}}_{{n}^{{\prime} }}^{* }{\tilde{u}}_{{n}^{{\prime\prime} }}+B{\tilde{u}}_{-n+{n}^{{\prime} }+{n}^{{\prime\prime} }}{\tilde{v}}_{{n}^{{\prime} }}^{* }{\tilde{v}}_{{n}^{{\prime\prime} }}\right){e}^{in\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau },$$
(16)

with a similar equation focusing on \({\tilde{v}}_{n}\). We can then decompose \(\tilde{u}\) and \(\tilde{v}\) into the product of two functions on distinct time scales by setting

$${\tilde{u}}_{n}(t)={U}_{n}(t){e}^{-{{{{{{{\rm{i}}}}}}}}{v}_{{{{{{{{\rm{g}}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}t},\quad {\tilde{v}}_{n}(t)={V}_{n}(t){e}^{-{{{{{{{\rm{i}}}}}}}}{v}_{{{{{{{{\rm{g}}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}t}$$
(17)

If we take an average of the decomposed Eq. (16) over an extended slow-time, t, interval then the exponential terms on the average vanish for n″ ≠ n. This turns Eq. (16) to

$$\mathop{\sum }\limits_{n=-\infty }^{\infty }\frac{\partial {U}_{n}}{\partial t}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }= \mathop{\sum }\limits_{n=-\infty }^{\infty }{\tilde{E}}_{in}{\delta }_{n,0}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }-\mathop{\sum }\limits_{n=-\infty }^{\infty }{U}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }\\ -{{{{{{{\rm{i}}}}}}}}\theta \mathop{\sum }\limits_{n=-\infty }^{\infty }{U}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }-{{{{{{{\rm{i}}}}}}}}\eta \mathop{\sum }\limits_{n=-\infty }^{\infty }{\left(in\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\right)}^{2}{U}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }\\ +{{{{{{{\rm{i}}}}}}}}\mathop{\sum }\limits_{n=-\infty }^{\infty }\mathop{\sum }\limits_{{n}^{{\prime} }=-\infty }^{\infty }\mathop{\sum }\limits_{{n}^{{\prime\prime} }=-\infty }^{\infty }(A{U}_{n+{n}^{{\prime} }-{n}^{{\prime\prime} }}{U}_{{n}^{{\prime} }}^{* }{U}_{{n}^{{\prime\prime} }}+B{V}_{n+{n}^{{\prime} }-{n}^{{\prime\prime} }}{V}_{{n}^{{\prime} }}^{* }{U}_{{n}^{{\prime\prime} }}){e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }\\ +{{{{{{{\rm{i}}}}}}}}\mathop{\sum }\limits_{n=-\infty }^{\infty }{U}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }\mathop{\sum }\limits_{{n}^{{\prime} }=-\infty }^{\infty }(2A{U}_{n{\prime} }{U}_{{n}^{{\prime} }}^{* }+B{V}_{{n}^{{\prime} }}{V}_{{n}^{{\prime} }}^{* })\\ +{{{{{{{\rm{i}}}}}}}}\mathop{\sum }\limits_{n=-\infty }^{\infty }{V}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau }\mathop{\sum }\limits_{{n}^{{\prime} }=-\infty }^{\infty }(B{U}_{{n}^{{\prime} }}{V}_{{n}^{{\prime} }}^{* }),$$
(18)

with again a similar equation focusing on \({\tilde{v}}_{n}\).

Finally, if we collapse the modal expansions by defining:

$${E}_{+}=\mathop{\sum }\limits_{n=-\infty }^{\infty }{U}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau },\quad {E}_{-}=\mathop{\sum }\limits_{n=-\infty }^{\infty }{V}_{n}{e}^{{{{{{{{\rm{i}}}}}}}}n\frac{2\pi }{{t}_{{{{{{{{\rm{R}}}}}}}}}}\tau },$$
(19)

one obtains the used model, Eq. (1).

$$ \frac{\partial {E}_{\pm }}{\partial t}={E}_{{{{{{{{\rm{in}}}}}}}}}-{E}_{\pm }-{{{{{{{\rm{i}}}}}}}}\theta {E}_{\pm }-{{{{{{{\rm{i}}}}}}}}\eta \frac{{\partial }^{2}{E}_{\pm }}{\partial {\tau }^{2}}\\ +{{{{{{{\rm{i}}}}}}}}\left(A| {E}_{\pm }{| }^{2}{E}_{\pm }+B| {E}_{\mp }{| }^{2}{E}_{\pm }+2A\langle | {E}_{\pm }{| }^{2}\rangle {E}_{\pm }+B\langle | {E}_{\mp }{| }^{2}\rangle {E}_{\pm }+B\langle {E}_{\pm }{E}_{\mp }^{* }\rangle {E}_{\mp }\right),$$
(20)

where the terms with angled brackets correspond to \(\langle f\rangle =(1/{t}_{{{{{{{{\rm{R}}}}}}}}})\int\nolimits_{-{t}_{{{{{{{{\rm{R}}}}}}}}}/2}^{{t}_{{{{{{{{\rm{R}}}}}}}}}/2}f\left(\tau \right)\,d\tau\), i.e., the temporal averages over the round-trip time tR. The presence of this kind of averaged terms in Fabry-Pérot configurations was first suggested by Firth in ref. 42 when describing the phase shift due to the counter-propagating fields in the Kerr medium.

Linear stability analysis

To derive Eq. (4), we focus on the nth term of Eq. (18) and linearized it around a homogeneous stationary state Un,HS where we defined Un(τ) = δn,0UHS + δUn(τ). Details of this method are discussed by Cole et al.27.