Hierarchy of coupled mode and envelope models for bi-directional microresonators with Kerr nonlinearity

We consider interaction of counter-propagating waves in a bi-directionally pumped ring microresonator with Kerr nonlinearity. We introduce a hierarchy of the mode expansions and envelope functions evolving on different time scales set by the cavity linewidth and nonlinearity, dispersion, and repetition rate, and provide a detailed derivation of the corresponding hierarchy of the coupled mode and of the Lugiato-Lefever-like equations. An effect of the washout of the repetition rate frequencies from the equations governing dynamics of the counter-propagating waves is elaborated in details.

A variety of models has been reported in the context of experiments dealing with a single mode operation in each direction [29][30][31][32][33][34]. We note, here that studies into single mode bidirectional lasers, laser gyroscopes and symmetry breaking in them have history going back to 1980's, see, e.g., [35][36][37]. To interpret recent soliton experiments, [15,18] have used models without nonlinear cross-coupling, while [16] has accounted for it. As we will see below, neglecting by the nonlinear cross-coupling was probably a better approach to analyse the experimental measurements under the circumstances, when modelling in neither of [15,16,18] included the effect of opposing group velocities, i.e., opposite signs of the resonator repetition rates for counter-propagating waves.
Due to complexity of the problem and diversity of equations both met in literature and the ones that are encountered during first principle analysis of the problem, it appears beneficial to have a detailed reference derivation that can be followed and tailored by a reader. Such mathematically transparent and physically motivated derivation that can be readily mapped onto a variety of experimental setups is present below. Focus of our work is to identify a hierarchy of the mode expansions and envelope functions evolving on different time scales set by the cavity linewidth, nonlinearity and 2nd order dispersion (slow time scales), and by the repetition rate (fast time scale), which can be used to derive a hierarchy of the coupled mode and envelope equations. We pay particular attention to comprehensive explanations of our derivation steps and interpretation of the results.

Hierarchies of mode expansions and envelope functions
This Section introduces physical system and discusses a hierarchy of mode amplitudes and envelope functions accounting for different time-scales. It also outlines plan of work for the rest of the paper.
Maxwell equations written for the electric field components E α using Einstein's notations read as Here ε is the dielectric response function varying in space and time. It is assumed to be scalar for the sake of brevity. θ ∈ [0, 2π) is the azimuthal coordinate varying along the ring circumference. z axis is perpendicular to the ring, while r = x 2 + y 2 measures distance from the ring centre. N α is the nonlinear part of the material polarization and c is the vacuum speed of light. We assume 3rd order nonlinearity, so that where α 1,2,3 and α represent either of the three Cartesian projections, x, y or z, of a physical quantity they are used with. An implicit summation is assumed over any repeated α ′ s. χ (3) is a 4th rank tensor of the third order nonlinear susceptibility, which is taken to be nondispersive (Kleinman condition, i.e., interchangeability of all four indices). Electric field vector E α inside a ring resonator is expressed as a superposition of its linear modes F αj (r, z)e ijθ±iω j t , which are solutions of Eq. (1) with N α = 0: Here j > 0 is an azimuthal mode number, or angular momentum, and ω j > 0 is the corresponding mode frequency. B ± j are the amplitudes of the clockwise (CW) and counter-clockwise (CCW) modes. For typical microresonators geometries, either bulk crystalline or chip integrated, the transverse mode profiles F αj can be divided into quasi-TE quasi-radial modes (|F x j,y j | ≫ |F z j |) and quasi-TM (|F z j | ≫ |F r j,θ j |). For many practical purposes, which is in our case calculation of the overlap integrals in the nonlinear terms, it often suffices to neglect the smaller components of F αj . We also assume that the dominant components of F αj and ω j are real, so that F αj = F * αj , ω j = ω * j . Thus for TE modes F x j ≈ cos θF r , F y j ≈ sin θF r , F z j ≈ 0 and for TM modes F x j,y j ≈ 0.
In order to cut notational complexity and drop the α index, we consider TM family, so that from now on F z j → F j , and Results of our derivations would be the same for TE modes, F r j → F j , and therefore, what we are loosing is only formal consideration of the nonlinear coupling between the TE and TM families.
We assume that inhomogeneities of the resonator surfaces result in scattering in general and in backscattering, in particular, and hence lead to the linear coupling between the modes. We account for these effects assuming Here ε id is the dispersive dielectric function of the ideal (no backscattering) geometry, that does not depend on θ, while relatively small ε in (θ) accounts for inhomogeneities along the ring. Mode profiles F j (r, z) are calculated for ε in = 0. E is measured in V/m, hence normalising linear modes as max r,z |F j | = 1 makes units of b j B ± j to be V/m. Real field amplitude of a CCW mode is 2b j |B + j |, so that its intensity is where r 0 is the distance between the z axis and a point of maximum of |F j |. We define scaling factors b j as so that the |B ± j | 2 are measured in Watts. ǫ vac is the vacuum susceptibility. We assume that the resonator is pumped into its j p mode and introduce the mode index offset µ = j − j p . The real field expression, Eq. (3), is then where j µ = j p + µ = j and µ = −| j p − j min |, . . . , 0, . . . , | j max − j p |. Here E(t, θ) is the real electric field measured in W 1/2 , which dependence on the transverse coordinates has been factored out. Introducing pump laser frequency Ω, we define mode detunings where δ 0 is the detuning for mode j p , and where a desired number of the dispersion orders can be included to approximate ω µ over a required spectral range. D 1 is the resonator repetition rate (or free spectral range (FSR)) and D 2 is its group velocity dispersion. D 2 > 0 implies anomalous and D 2 < 0 normal dispersion. For example, the work [15] deals with a bi-directionally pumped silica ring with radius 1.5mm and it has D 1 = 2π × 22GHz, D 2 = 2π × 16kHz. The linewidth of this resonator is κ = 2π × 1.5MHz, and hence the corresponding finesse F = D 1 /κ ≃ 13000. The mode area estimate is S j p ≃ 30µm 2 , which gives b 2 j ≃ 4 × 10 12 V 2 W −1 m −2 . Pump laser wavelength was ≃ 1550nm (ω 0 ≃ 2π × 193THz) and the comb spectra observed there were relatively narrow and span over ∼ 20nm bandwidth, corresponding to about 300 modes, and the momentum of a mode nearest to the pump is estimated as j p = 8700.
In order to introduce a new set of mode amplitudes important in what follows, we transform Eq. (7) further: Newly introduced mode amplitudes Q ± µ are defined as and as we can see they absorb frequency scales associated with both D 1 and D 2 . The corresponding CW and CCW envelope functions are Inclusion of the backscattering effects to the envelope equations, see Section 5, requires introducing of the envelope functions with reflections of their spatial coordinate, The above definition of the space reflected functions follows a text-book list of properties of Fourier transforms, where an equivalent transformation is typically introduced in time domain and could be called as either time reflection or time inversion transformation. Differential equations involving functions with reflections of their arguments also attracted some attention from a more general mathematics prospective, see, e.g., [38], while our system reveals their role in nonlinear photonics. In order to take control of D 1 in our future calculations, we define yet another set of slow amplitudes Here D 1 is moved away from the exponential factors defining our third and final set of amplitudes A ± µ . Instead, exponents with D 1 appear explicitly in the total field equation that uses A ± µ , We also use A ± µ to define the corresponding envelope functions and their reflections Though the envelopes A ± can not be used themselves to define the electric field E (only their mode amplitudes can), cf. Eq. (15) and (16), they play a pivotal role in the transition from the coupled mode to the partial differential equations, see Section 5.
To summarize this section: B ± µ amplitudes absorb only the slowest time scales associated with the nonlinear effects and resonator losses. A ± µ absorb time scales associated with the second and higher order dispersions, in addition to the ones already inside B ± µ . Q ± µ amplitudes evolve with the highest in our hierarchy frequency determined by the resonator repetition rate. To see how these different time scales and mode amplitudes are used to express the total field, E, one should compare Eqs. (10a), (10b) and (15).
The rest of this work is structured as follows: In section 3, we first derive a system of equations for B ± µ and perform its exact reduction to the equations for Q ± µ . In Section 4, we come back to the equations for B ± µ , make the D 1 role explicit, eliminate the associated fast oscillations and derive a simpler system for A ± µ . Corresponding mean-field equations for the envelope functions Q ± and A ± and their counter parts with the reflected spatial coordinates are derived in Section 5.

Separating equations for CW and CCW amplitudes
. Neglecting all the 2nd and higher order time derivatives of B ± µ we find that Eq. (1) transforms to where s µ = n 2 µ + 1 2 ω µ ∂ ω n 2 µ ≃ n 2 µ . We then expand nonlinear polarization N in Fourier series In order to carry out separation of the CW and CCW equation we also need to define CW and CCW components of nonlinear polarization, P ± j µ e ∓iω jµ t , such that Explicit expressions for P ± j µ are given by Eqs. (29) below. We now multiply the left and right hand-sides of Eq. (17) by b j p F j p exp −ij µ ′ θ , integrate in r, z and θ, and approximate ω µ ≃ ω 0 = ω j p , n µ ≃ n 0 inside all the pre-factors, but not in the powers of the exponents. The resulting model, see Eqs. (22), makes use of the two scattering matrices having dimensions of angular frequencies. One characterises scattering induced coupling between the co-propagating modes and the other one describes backscattering induced mode coupling, The projected equation itself is where V p = 2π ∬ F 2 j p rdrdz is the mode volume for j = j p . Eq. (22) can now be split, as per rotating wave approximation, into the parts proportional to e ±iω jµ t exponents, so that we have two equations defined on the slow, D 2 related, time scales: where we have also swapped µ and µ ′ . In order to be used to describe laboratory experiments with microresonators, Eqs. (23) have to be amended with the single mode pump term and losses accounting for the finite linewidth. We take, for the laser frequency at the exact cavity resonance Ω = ω µ=0 = ω j=j p and for the low pump levels, i.e., linear regime, the intracavity powers of CW and CCW waves to be |H ± | 2 = |B ± µ | 2 . This is achieved via a phenomenological substitution Here, Kronecker delta is defined as δ µ,µ 1 = 1 for µ = µ 1 and is 0 otherwise. If pump is absent, then the field power would decay with the rate κ (full width of the resonance). An expression linking H ± with the laser powers W ± is where W ± are the laser powers pumping, respectively, CW and CCW waves. η < 1 is the coupling efficiency via, e.g., a prism or a waveguide, into a resonator mode. η = κ c /κ, where κ c is the coupling pump rate (equals coupling loss rate). F /π is the cavity induced power enhancement. Detailed theoretical and experimental studies of the power enhancement effect and coupling in and out considerations for ring cavities can be found in, e.g., [39,40]. R µµ ∼ 2π × 4 kHz in Ref. [15]. In this regime, it is safe to assume that κ dominates over Γ and R terms. Using this we disregard Γ µ ′ µ in what follows, and retain only the dominant diagonal terms in R µ ′ µ , i.e., R µ ′ µ µ ′ ≈ 0. Dispersion of the diagonal terms is also disregarded, R µµ ≃ R 00 = R. Accounting for all of the above and complex conjugating second of Eqs. (23) we conclude this subsection with

Opening up nonlinearity
Using Eqs. (4), (7) we have and Comparing Eqs. (27), (28) and Eqs. (18), (19), one can define explicit expressions for P ± j µ . Assuming spectrally narrow combs, and therefore omitting all terms with exponential factors oscillating in space with multiples of j p and in time with multiples of Ω, we find where nonlinear coefficient is Total refractive index, n, for a single mode operation is n = Assuming that the j p mode is well confined within the resonator material, the mode shape can be approximated by a Gaussian function (allowing for different widths along z and x), and rdr ≈ r 0 dx (see text before Eq. (6)), gives 2π ∬ Eq. (32) and Eq. (33) have been compared using mode profiles calculated with Comsol and it was found that the latter provides a very practical approximation. For ω 0 = 2π × 193THz, n 0 = 1.47 and n 2 ≃ 3.2 × 10 −20 m 2 /W 2 (silica glass), and mode area S j p ≈ 30µm 2 we have γ ≃ 2π × 70kHz/W. S j p is an order of magnitude smaller and n 2 is an order of magnitude larger in integrated Si 3 N 4 microresonators, and their combined effect boosts γ up by two orders of magnitude.
Using Eqs. (11) to express amplitudes B ± µ via Q ± µ we find that all the time dependent exponents cancel out and the resulting coupled mode equations for Q ± µ amplitudes are where Q ± envelopes are given by Eqs. (12).

Washout of the repetition rate timescales from the coupled mode equations
Systems of Eqs. (30), (11), (12) on one side, and Eqs. (34), (12) on the other, are mathematically and physically equivalent. However, there are important observations to be made here. If one could assume that |Q + | 2 + 2|Q − | 2 under the integrals in the right hand sides of Eqs. (30) and Eqs. (34) is a slow function of time, then these integrals would be approximately equal to Q ± µ e −iδ µ t , see Eqs. (11). Balancing these with the e iδ µ t exponents before the integrals in Eqs. (30), one would end up with equations involving time scales determined only by the linewidths, pump detuning and nonlinear resonance shifts, which are all order of MHz. MHz frequencies would be far simpler to resolve numerically, compare to GHz-THz frequencies associated with D 1 , that are directly implicated inside δ µ in the linear parts of Eqs. (34).
In this Section, we demonstrate that there are both slow and fast time scales inside the nonlinear terms in Eqs. (30), and that the latter can be eliminated resulting in a simpler and better balanced system of equations for the A ± µ amplitudes, see Eqs. (14), (39). We proceed by taking Eqs. (30a), express Q ± via B ± µ , see Eqs. (11), (12), calculate integrals in the nonlinear terms, see Eq. (35a), and perform the two step transformation, see Eqs. (35b), (35c), The four-wave mixing momentum matching conditions are reflected in the Kronecker delta's in front of the nonlinear terms in the second line of the above and directly follow from taking the integrals in θ. Swapping of µ 2 and µ 3 inside the nonlinear cross-coupling is a critical step that a reader should pay attention to, see Eq. (35b). This operation equals the Kronecker delta's, but it re-orders the amplitudes and respective frequency detunings in the second nonlinear term. After inserting explicit expressions for δ µ , see Eqs. (8), (9), and using the momentum matching condition, we find that D 1 frequencies cancel out inside the nonlinear self-action terms, but remain in the cross-action ones providing µ 2 µ 3 , see Eq. (35c). Thus if D 1 oscillations are much faster than dynamics associated with the other time scales left in the equations, i.e., κ, R and nonlinear frequency shifts, then the fast oscillating components can be disregarded [19,41]. This leaves us only with µ 2 = µ 3 components in the cross-action terms, so that Nonlinear terms in Eqs. (37) are now grouped into the phase insensitive pure cross-Kerr term, that contains nonlinear shift of the CW resonance frequencies due to CCW wave, and into the term that mixes both phase sensitive and phase insensitive four-wave mixing CW-CW nonlinearities. The phase sensitive effects come only from the CW-CW interaction, because all the phase sensitive CW-CCW dynamics develops with the 2D 1 frequencies and is washed out by the high repetition rates. This can be called the washout effect of high repetition rates on nonlinear frequency mixing of the counter-propagating waves in a ring resonator. Using A ± µ amplitudes and detunings δ ′ µ , which are both D 1 free, see Eqs. (14), allows to hide e iD 2 µ 2 t/2 exponents in Eqs. (37). Adding the CCW equation, we have The difference of the above nonlinear terms with the ones in the equations for Q ± µ , see Eq. (34), that include un-averaged D 1 oscillations, becomes more obvious, if the sums in Eqs. (38) are replaced with the integrals, see also Eq. (16a), The last terms in Eqs. (39) follow from the Parseval's theorem. Thus Eqs. (37) include only effects of the second and higher order dispersions in both linear and nonlinear terms, that in microresonators are associated with the kHz to MHz time scales. Hence solving Eqs. (37) is expected to provide significant computational advantages over all other versions of the coupled mode equations.

Envelope models
Connection of the coupled mode equations to the wave dynamics becomes more intuitive, if one now derives the envelope, Lugiato-Lefever like, equations. First, we take the Q ± µ model, see Eqs. (34), and multiply Eq. (34a) with e iµθ and (34b) with e −iµθ . We then sum up each of the equations in µ and use Eqs. (12), (13) connecting the envelopes Q ± and the reflected envelopes Q (r) ± to their mode amplitudes. This procedure is free from approximations and it leads to a system of partial differential equations for Q ± and Q (r) ± , To form a closed system, the above pair of equations should be supplemented with two more equations for the Q (r) ± , see Eqs. (13) defining θ reflection. Starting from the equations for A ± µ , see Eqs. (39), we follow a modified procedure. Namely, we multiply both CW and CCW equations by the same exponent e iµθ , use the envelope definitions in Eqs. (16), observe that ∫ 2π , sum up in µ, and derive the following envelope equations The above equations do not only exclude the D 1 dynamics, but also form a closed system of two equations for the CW A + envelope and for the reflected CCW A (r) − envelope, see Eqs. (16). They can also be supplemented with equations for A (r) + , A − , but this time those are left as an independent pair. Again numerical modelling of Eqs. (41) is expected to have great advantages relative to working with Eqs. (40). Similar to ours procedure to remove the D 1 linked time scales has been developed for the Kerr Fabry-Perot cavities supporting a single family of standing waves and hence yielding a one-component Lugiato-Lefever model [41]. The respective ring geometry model in [19] mixes all four envelope functions, i.e., A ± , A (r) ± , and is limited by the second order dispersion.
Eqs. (40) (not Eqs. (41)) could in fact, be written without a rigorous derivation, by simply relying on common knowledge, let aside reflected envelopes in the backscattering terms. These equations include traditional cross-phase modulation, and also repetition rates terms and other odd order dispersion terms with the opposite signs. Contrary, Eqs. (41) have no repetition rate terms, i.e., D 1 -terms, and the remaining odd dispersions, i.e., D 3 , D 5 , etc., come with the same signs. Simultaneously, phase sensitive nonlinear wave mixing effects induced by CW-CCW interaction have been washed out. The only nonlinear cross-interaction left comes from the integrated power, which merely shifts the detuning parameters. Thus, in the absence of backscattering a nonlinear bi-directional resonator operates as a uni-directional one, but with the detuning parameter altered by the total power of the counter-propagating wave.

Summary
We have derived coupled mode equations describing nonlinear wave mixing processes in Kerr microresonator with counter-propagating waves. Features of the first two coupled mode for-mulations given by Eqs. (30) and Eqs. (34) are that they fully account for the repetition rate effects and that nonlinear terms are taken in the real space, and can be evaluated via Fourier transforms, see also [42]. We then proceeded to present simplified multi-mode equations that neglect the repetition rate dynamics driving the phase sensitive terms responsible for nonlinear interaction between the counter-propagating fields (washout effect, section 4), and again deal with the nonlinearity in the real space, see Eqs. (39).
Finally, we demonstrated that coupled mode equations (34) and Eqs. (39) are equivalent to two different, Lugiato-Lefever-like, envelope models. The one that involves the repetition rate dynamics, see Eqs. (40), links two usual envelopes for the CW and CCW fields, with two of their space reflections. While the one with the repetition rate averaged out, see Eqs. (41), makes a closed system already for two envelopes, A ± , one of which is reflected. We note, that Q ± can be used directly to reconstruct total electric field, see Eq. (10c), while A ± can not, but their respective mode amplitudes can, see Eqs. (15), (16a).
We have taken care to reveal all mathematical transformations, that allow a reader to verify our derivation steps and apply modifications if required. Opportunities for future theoretical and numerical studies offered by the models presented here are numerous, as well as their potential to guide and interpret experimental work.