Aurélien Coillet

Rémi Henriet

Patrice Salzenstein

Kien Phan Huy

Laurent Larger

Yanne K Chembo
email: yanne.chembo@femto-st.fr

* A Coillet

Time-domain Dynamics and Stability Analysis of Optoelectronic Oscillators based on Whispering-Gallery Mode Resonators

Optoelectronic oscillators (OEOs) are microwave photonics systems intended to generate ultrastable radio-frequency signals with unprecedented phase noise performance for aerospace and communication engineering applications. They had originally been introduced in a configuration where the energy storage element was a fiber delay line. However, recent research in view of size and power consumption optimization has led to novel configurations where this fiber delay line is replaced by a ultrahigh Q whisperinggallery mode (WGM) resonator. So far, there has been no theoretical framework enabling to understand the dynamical behavior of these new architectures of OEOs. In this paper, we propose for the first time a deterministic time-domain model to investigate the dynamics of these OEOs based on WGM resonators. This model enables us to perform the stability analysis of the microwave oscillations, and to determine rigorously their range of stability as the loop gain is varied. After building the model, we perform a full stability analysis of the various stationary solutions for the microwave output. We then perform extensive numerical simulations, which are in complete agreement with the stability analysis. The theoretical analysis is also found to be in excellent agreement with our experimental measurements.

Introduction

The optoelectronic oscillator (OEO) is nowadays considered as one of the most promising ultra-stable microwave generator for applications in time-frequency metrology, frequency synthesis, detection and navigation systems [START_REF] Maleki | The optoelectronic oscillator[END_REF]. The first architecture of OEO, proposed by Yao and Maleki [START_REF] Yao | High frequency optical subcarrier generator[END_REF][START_REF] Yao | Optoelectronic microwave oscillator[END_REF][START_REF] Yao | Optoelectronic oscillator for photonic systems[END_REF], performed energy storage in the feedback loop by using a few kilometer-long fiber-delay line instead of a high-finesse radio-frequency (RF) filter. The idea to store laser light energy instead of microwave energy was a conceptual breakthrough which provided a technological pathway towards improved stability for RF signal generators. Another noteworthy feature of OEOs is their nearly absolute frequency versatility, since the microwave can be arbitrarily set at any frequency belonging to the band 0.1-100 GHz. The upper limit of this frequency band is in fact imposed by currently available optoelectronics components, and nothing theoretically prevents OEOs to generate millimeter waves as well.

Despite their excellent stability performances, the original fiber-based OEOs unfortunately have the disadvantage to be bulky because of the temperature-stabilized box containing the optical fiber delay line. Effectively, these fiber-based OEOs are not very transportable, and their weight might affect negatively the payload of spacecrafts. Their size (few dm 3 ) also raises problems related to temperature stabilization, which becomes overly energy-greedy in this case. In addition, despite active stabilization, the long fiber delay line also induces an unavoidable phase drift that deteriorates the longterm stability of the oscillator. The ring-cavity modes induced by the fiber delay line are also arising as very strong (even though narrow) parasite peaks close to the carrier in the phase noise spectrum. These spurious peaks are indeed very detrimental in most applications.

In order to improve the features of this oscillator, many novel architectures have been proposed in recent years. Some of them involve multiple loops in order to suppress the spurious peaks [START_REF] Yao | Dual microwave and optical oscillator[END_REF][START_REF] Okusaga | Spurious mode reduction in dual injection-locked optoelectronic oscillators[END_REF]. Others lock the oscillator to atomic resonances [START_REF] Strekalov | Stabilizing an Optoelectronic Microwave Oscillator with Photonic Filters[END_REF], or enhance the functionalities of the optical branch to improve the phase noise figure [START_REF] Devgan | Improvement in the Phase Noise of a 10 GHz Optoelectronic Oscillator Using All-Photonic Gain[END_REF], the tunability of the oscillator [START_REF] Li | An Optically Tunable Optoelectronic Oscillator[END_REF], or to mode-lock the optical modes [START_REF] Lasri | Self-starting optoelectronic oscillator for generating ultra-low jitter high-rate (10 GHz or higher) optical pulses[END_REF][START_REF] Yu | Ultralow-noise mode-locked laser with coupled optoelectronic oscillator configuration[END_REF] for ultra-low jitter pulse generation.

However, one of the most interesting architecture to overcome the shortcomings of fiber-based OEOs is undoubtedly the configuration which replaces the fiber delay-line by an ultra-high Q WGM resonator (see, e. g., refs. [START_REF] Matsko | Whispering gallery mode based optoelectronic microwave oscillator[END_REF][START_REF] Savchenkov | Whispering-gallery mode based opto-electronic oscillators[END_REF][START_REF] Ilchenko | Optical Resonators With Whispering-Gallery Modes-Part II: Applications[END_REF][START_REF] Merrer | Characterization technique of optical whispering gallery mode resonators in the microwave frequency domain for optoelectronic oscillators[END_REF]). Whispering gallery mode resonators are low-loss dielectric disks or rings that perform optical energy storage through trapping photons by total internal reflection [START_REF] Matsko | Optical Resonators With Whispering-Gallery Modes-Part I: Basics[END_REF][START_REF] Chiasera | Spherical whispering-gallery-mode microresonators[END_REF]. The optical Q factor of these resonators can be defined as Q = ω 0 /∆ω, where ω 0 is the central angular frequency of the mode of interest, and ∆ω is the corresponding linewidth. When the resonators are almost perfectly shaped (with sub-nanometer surface roughness) with a ultra-low loss bulk material (fused silica, calcium or magnesium fluoride crystals, etc.), they can achieve a quality factor that is typically above 10 8 at 1550 nm; they can even exceptionally reach record values higher than 10 11 [START_REF] Savchenkov | Optical resonators with ten million finesse[END_REF]. The linewidth ∆ω = ω 0 /Q of these WGMs is typically of the order of 100 kHz for Q ∼ 10 9 at 1550 nm. These linewidths are at least two orders of magnitude narrower than typical RF filters. Al- The "through" port of the coupling fibers is used to monitor the optical spectra with the OSA, in both the transient and stationary regimes. In the electric branch, a fast oscilloscope enables to resolve the temporal dynamics of the microwave V(t), whose complex envelope is V (t) (or A(t) in the dimensionless form). An ESA also enables to monitor the corresponding RF spectrum.

In the open-loop configuration (no oscillations), the "through" port is also used to perform the cavity-ring down measurement using an oscilloscope, thereby enabling the determination of the intrinsic and in-coupling quality factors.

ternatively, the optical storage properties of WGM resonators can also be understood in terms of photon lifetime τ ph = Q/ω 0 , which is of the order of 1 µs for Q ∼ 10 9 at 1550 nm.

The modal structure of these resonators is such that their eigenmodes are grouped within families where WGMs are nearly equidistant when dispersion is neglected. For WGMs belonging to the same family, the intermodal frequency (which is sometimes referred to as free spectral range, or FSR) is a free parameter that depends on the resonator's principal radius. It may vary from GHz (millimeter-size disks) to THz (micrometer-size) frequencies. Therefore, in the case of OEOs based on WGM resonators, we can obtain a microwave oscillation by extracting the intermodal frequency of an optically pumped WGM resonator, while the energy storage would be performed by trapping photons in the long-lifetime WGM cavity. This architecture solves almost all the problems raised by fiber-based OEOs, and yields an oscillator that is versatile, compact, energy efficient, and free from parasite spurious peaks (see, e. g., ref. [START_REF] Volyanskiy | Compact optoelectronic microwave oscillators using ultra-high Q whispering gallery mode disk-resonators and phase modulation[END_REF]). Most importantly, it is expected that these WGM-based architectures would enable us to downsize the OEO from a shoe-box (fiber based OEO) to a smartphone (WGM-based OEO), and ultimately, to a match-box (integrated WGM-based OEO), without deteriorating the phase noise performance in the most favorable case.

Even though it appears clearly that WGM-based OEOs will play an increasingly im-portant role in microwave photonics, there is currently no dynamical framework to study their dynamics. This lack of analytical insight into the dynamical properties of this oscillator does not enable us to optimize its properties. In particular, it is impossible to know under which conditions the oscillator is merely stable, and it is worth reminding that this stability problem is far from being trivial. Our aim in this paper is, therefore, to propose a deterministic model to understand and analyze the time-domain dynamics of WGM-based OEOs. Using this model, we will study the stability of the various solutions, and perform numerical simulations to confirm the stability analysis study. We have also built a new architecture of WGMbased OEO, and have performed experimental measurements to check the validity of the model. Both the theoretical and experimental studies agree with excellent precision, thereby confirming the validity of the nonlinear dynamics framework of analysis.

This paper is organized as follows. The next section is devoted to the description of the experimental system under study, which is an architecture of WGM-based OEO with amplitude modulation and add-drop coupling. In Section 3, we present the main lines of our theoretical approach for OEOs, which is based on the nonlinear dynamics of the slowly-varying microwave envelope, while Section 4 is focused on the construction of the theoretical model for this WGM-based OEO. Then, we determine the stationary states of the oscillator in Section 5, while their stability is investigated in Sections 6 and 7. Hence, we perform numerical simulations that are compared to experimental measurements in Section 8. The last section concludes this paper with a resume of our work and perspectives of future research.

The experimental system

The WGM-based OEO under study is displayed in Fig. 1. The various elements of this single-loop architecture are as follows.

• A continuous-wave (CW) distributed feedback (DFB) semiconductor laser of optical power P L and whose central wavelength is λ L = 1552.2 nm, corresponding to an angular frequency ω L = 2πc/λ L , where c is the velocity of light in vacuum.

• An erbium-doped fiber amplifier (EDFA) delivering a maximal optical power of 30 dBm, and of optical gain G o when driven by the input laser.

• A wideband integrated optics LiNbO 3 Mach-Zehnder (MZ) intensity modulator characterized by the half-wave voltages V π DC = 4 V and V π RF = 4.7 V.

• A polarization controller to tune the polarization at the input of the MZ intensity modulator and the WGM resonator.

• A crystalline calcium fluoride (CaF 2 ) WGM resonator coupled in the add-drop configuration. The intrinsic, excitation and drop Q factors are respectively Q i , Q e and Q d . They respectively correspond to photon lifetimes Q i,e,d /ω L = τ i,e,d /2, where ω L is the angular frequency of the laser while the τ parameters are the laser field decay times. The internal quality factor Q i is fixed and equal to 4 × 10 8 , while the coupling quality factors Q e and Q d can be varied by changing the position of the coupling tapered fibers relatively to the resonator; in all case, these coupling Q-factors are of the order of 10 8 . The WGM disk has a refractive index equal to n g = 1.43 and a principal radius a = 3.2 mm. In consequence, its free spectral range is a microwave frequency equal to Ω M = c/an g = 2π × 10.4 GHz.

• A fast photodiode with a conversion factor S = 50 Ω × 0.75 A/W = 37.5 V/W, and bandwidth 0-12 GHz.

• A narrow-band microwave RF filter of central frequency at 10.5 GHz and a bandwidth of 1 GHz, intended to reject the RF noise outside the frequency band of interest.

• Cascaded microwave amplifiers with overall gain G e are used to close the loop and drive the intensity modulator.

The principle of operation of this WGM-based OEO is, therefore, the following: noise in the microwave branch of the loop modulates a laser light beam, and generates a broadband optical spectrum at the output of the MZ modulator. This broadband, initially noisy optical spectrum is narrowly filtered by the WGM resonator. At the output, the optical spectrum is now a set of equidistant spectral lines, separated by the FSR of the resonator (that is, 10 GHz). The photodiode detects the various beating frequencies kΩ M (k being an integer), but because of its limited bandwidth, it filters them out except the fundamental frequency Ω M and the dc component of the field. The RF filter rejects this dc component. Then, the microwave signal at kΩ M is finally amplified and sent as a driving signal to the MZ in order to close the feedback loop and the same cycle is started again.

It is known from oscillator engineering theory that the oscillation might be sustained when the Barkhausen conditions are fulfilled, that is, when 1) the loop gain overcomes the loop losses and 2) the round-trip phase of the microwave signal is null modulo 2π. However, the Berkhausen theory can not describe what might occur above threshold (amplitude of the oscillations, multi-stability, hysteresis, higher-order bifurcations, chaos, etc.), nor does it enable the stability study of the oscillating solutions. The purpose of the nonlinear dynamics approach is to shed the light on all these blind spots, as it provides a complete understanding of the oscillator behavior below and above threshold. We will explain in the next section how this nonlinear dynamics is used in the context of OEOs in general, and WGM-based OEOs in particular.

Modelling OEOs: a microwave envelope approach

A microwave oscillator whose output angular frequency is around Ω 0 is always expected to have an output of the form

V(t) = A(t) cos[Ω 0 t + ψ(t)], (1)
where A(t) and ψ(t) are, respectively, the amplitude and the phase of the microwave. In the case where A and ψ are constant, the microwave output V(t) is perfectly sinusoidal. The microwave is still sinusoidal as well if ψ = ψ 0 + σt, even though the oscillation frequency is shifted to Ω 0 +σ. In the general case, both A and ψ will be time-dependent, but their time variation will be significantly slower than the period of the oscillation. More concretely, A(t) and ψ(t) will vary at a slow time-scale comparable to the inverse RF bandwidth 1/∆Ω of the oscillation loop, while the full microwave V(t) will vary at a fast time-scale of 1/Ω 0 . Therefore, the slow and fast dynamics are split in by a factor Q RF = Ω 0 /∆Ω corresponding to the RF quality factor of the oscillator. This is why A(t) and ψ(t) are referred to as slowly varying amplitude and slowly-varying phase, respectively.

In general, the only variables of interest are A and ψ since Ω 0 is known. Equation ( 1) can be rewritten under the useful form

V(t) = 1 2 V (t)e iΩ 0 t + 1 2 V * (t)e -iΩ 0 t , ( 2
)
where V (t) = A(t)e iψ(t) , and the star denotes the complex conjugation. Here, the complex valued variable V (t) synthetically gathers all the information about the amplitude and the phase of the microwave: it is referred to as the complex slowly varying envelope of the microwave. Hence, it is the idoneous variable to investigate the dynamics of the RF output of a microwave oscillator. It is interesting to note that complex slowly varying envelopes are routinely used in optics and laser theory, where the nominal frequencies (from the lasers, optical cavities, etc.) are generally known, while the optical amplitudes and phases are the variables of interest. It is also worth noting that the complex envelope takes into account eventual frequency detunings σ from the nominal frequency Ω 0 through ψ(t).

The microwave envelope approach has already been used with great success in fiberbased OEOs. In fact, fiber-based OEOs belong to the large family of electro-optic systems with delayed feedback, which can be described by Ikeda-like delay differential equations where the main variable is the real-valued voltage at the input of the modulator (see [START_REF] Larger | Complexity in electro-optic delay dynamics: modeling, design, and applications[END_REF]). In the particular case of fiber-based OEOs, the narrowband filtering around the frequency Ω 0 of interest enables us to rewrite this voltage under the form of (2). A nonlinear delay differential equation ruling the dynamics of the complex microwave envelope had, therefore, been obtained, thereby enabling to demonstrate many fundamental results. Just to name a few, this nonlinear dynamics approach enables us to prove that fiber-based OEOs can turn unstable if the feedback gain exceeds a precise bifurcation value [START_REF] Chembo | Dynamic instabilities of microwaves generated with optoelectronic oscillators[END_REF][START_REF] Chembo | Nonlinear dynamics and spectral stability of optoelectronic microwave oscillators[END_REF]. The same approach also enabled to show that under certain conditions, the abrupt switch-on behavior of the fiber-based OEOs leads to robust multimode oscillations instead of an ultrastable single-tone microwave [START_REF] Chembo | Effects of gain and bandwidth on the multimode behavior of optoelectronic microwave oscillators[END_REF]. The time-domain deterministic model was also an essential prerequisite that enabled to perform a phase noise analysis based on stochastic differential equations (or Langevin equations) and which enabled to predict phase noise characteristics with remarkable precision [START_REF] Chembo | Determination of phase noise spectra in optoelectronic microwave oscillators: a Langevin approach[END_REF][START_REF] Volyanskiy | Contribution of Laser Frequency and Power Fluctuations to the Microwave Phase Noise of Optoelectronic Oscillators[END_REF]. The same formalism also enabled us to analyze more complex OEO architectures, like dual-loop OEOs [START_REF] Nguimdo | On the Phase Noise Performance of Nonlinear Double-Loop Optoelectronic Microwave Oscillators[END_REF], or hybrid configurations whose outputs are an ultra-stable microwave in the RF domain and an ultralow jitter picosecond pulse train in the optical domain [START_REF] Chembo | Generation of Ultralow Jitter Optical Pulses Using Optoelectronic Oscillators With Time-Lens Soliton-Assisted Compression[END_REF].

We show in the next section that a complex envelope formalism can be developed as well for WGM-based OEOs.

Model

The dynamics of this system is essentially defined by two variables. The first variable is the laser electric field E(t) at the input of the resonator and the second variable is the input microwave voltage V(t) of the integrated modulator.

Instead of working directly with the real-valued E(t) and V(t), we use their complex slowly varying amplitudes

E(t) = |E(t)|e iϕ(t) and V (t) = |V (t)|e iψ(t) defined through E(t) = 1 2 E(t)e iω L t + 1 2 E * (t)e -iω L t V(t) = 1 2 V (t)e iΩ M t + 1 2 V * (t)e -iΩ M t , (3)
where ω L and Ω M are the angular frequencies associated with the 1550 nm infrared laser beam, and to the 10 GHz microwave signal, respectively. The slowly varying amplitude of the optical beam at the input of the MZ modulator simply reads E cw = √ P 0 , where P 0 = G o P L is the optical power at the output of the EDFA. This pumping field sets the optical phase reference, and as a consequence is real (its phase is null). This beam is amplitude-modulated with a driving RF signal

V(t) = |V (t)| cos[Ω M t + ψ(t)]
, so that the slowly varying amplitude of the optical field at the input of the optical fiber is

E(t) = E cw cos πV(t) 2V π RF + πV B (t) 2V π DC (4) = P 0 cos π|V (t)| 2V π RF cos[Ω M t + ψ(t)] + πV B 2V π DC .
The Jacobi-Anger expansion gives

e iz cos α = +∞ n=-∞ i n J n (z)e inα (5)
where J n is the n-th order Bessel function of the first kind. Therefore, we have

E(t) = P 0 +∞ n=-∞ E n (t)e inΩ M t (6)
where

E n (t) = n (φ) J n π|V (t)| 2V π RF e inψ(t) , (7)
and

n (φ) = 1 2 [e iφ + (-1) n e -iφ ]i n (8) = (-1) n 2 cos φ if n is even (-1) n+1 2 sin φ if n is odd .
The parameter

φ = πV B 2V π DC (9)
is the offset phase due to the bias voltage of the integrated MZ modulator. The intracavity field F (t) inside the resonator can also be spectrally decomposed as

F (t) = +∞ n=-∞ F n (t)e inΩ M t
. According to the Haus formalism [START_REF] Haus | Waves and fields in optoelectronics[END_REF], the dimensionless components F n obey

dF n dt = - 1 τ F n -iσF n + 2 τ e E n , ( 10
)
where τ is defined as

1 τ = 1 τ i + 1 τ e + 1 τ d , (11)
and stands for the overall loss-induced decay time for the electric fields inside the resonator, while σ = ω Lω 0 is the laser detuning relatively to the central frequency of the pumped mode. On the other hand, the dimensionless components G n of the output field G(t) = +∞ n=-∞ G n (t)e inΩ M t can simply be recovered as

G n = 2 τ d F n . ( 12
)
The spectra of the three optical fields E(t), F (t), and G(t) have been schematically represented in Fig. 1.

The optical power at the input of the photodiode is equal to

P 0 |G(t)| 2 = P 0 +∞ n=-∞ G n (t)e inΩ M t 2 (13) = 1 2 C 0 (t) + +∞ k=1 1 2 C k (t)e ikΩ M t + c.c. ,
where the slowly varying Fourier coefficients C k express the multifrequency nature of the input optical power. The photodiode has an inbuilt filter that rejects the harmonics at kΩ M with k ≥ 2. On the other hand, the dc component is rejected by the RF bandpass filter. Hence, only the spectral component at frequency Ω M (that is, C ±1 (t)) is allowed to pass through. The slowly varying amplitude V of the voltage at the output of the photodiode is, therefore, SC 1 (t), where S stands for the photodiode sensitivity. We can multiply this output by the overall gain G = G e G o and the overall losses κ to obtain the slowly-varying voltage at the input of the modulator as

V (t) = κGSe iς C 1 (t) = κGSP 0 +∞ n=-∞ 2G n+1 G * n , ( 14
)
where the phase factor e iς accounts for the effect of the microwave round-trip phase shift ς. If necessary, it can be tuned to any desired value (modulo 2π) using a RF phase shifter.

We define the dimensionless microwave voltage as

A(t) = π 2V π RF V (t), (15)
and the optoelectronic gain as

β = πκGSP 0 2V π RF . ( 16
)
It is very important to note that since the loss parameter κ does not consider the losses associated with the WGM resonator, the gain coefficient β is not a loop-gain parameter as it is the case for fiber-based OEO studies, where κ also takes into account the losses induced by the fiber. We therefore have the following three-step model for numerical simulation:

E n = n (φ) Jc n [|A|] A n (17) Ġn = - 1 τ + iσ G n + 2 √ τ e τ d E n ( 18
)
A = 2βe iς +∞ n=-∞ G n+1 G * n , ( 19
)
where the overdot stands for time derivative. Note that in [START_REF] Chiasera | Spherical whispering-gallery-mode microresonators[END_REF], we introduced the Bessel cardinal functions to simplify the expression. Their definition and properties are further developed in the Appendix. The model describes the following phenomenology. Small noise in A generates the fields E n , which are fed in the WGM resonator and yield the output fields G n . The photodiode extracts the intermodal frequency and the bandpass RF filter outputs a microwave of complex envelope A, which is plugged back to the modulator to generate new fields E n , and the previous sequence of events takes place again.

It is interesting to note that fiber-based OEOs have an optical delay line that performs the optical energy storage, and a RF filter to select the microwave oscillation frequency. On the other hand, WGM-based OEOs rather have a WGM disk that performs at the same time the optical storage and filtering functions. The physics is intrinsically different, and so are the corresponding models: this explains why the former models built

A st2 A st3 A st4 A max1 st A st J 1 (2A st ) and A st /Γ Figure 2:
Geometrical interpretation of the stationary states A st , which are given by the intersection of the functions J 1 (2A st ) and A st /Γ . When Γ = 0.8, both curves only intersect for A st = 0, which is therefore the unique fixed point. This will be the case whenever Γ < 1. For Γ > 1, both functions will intersect for other points than 0, thereby generating non-trivial equilibria. For Γ 2 = 1.4, the non-trivial solution is A st2 = 0.79, and it will be equal to A st3 = 1.46 for Γ 3 = 4. When the gain is increased to Γ 4 = 15.52, the same branch yields a non-trivial solution A st4 = 1.77, while a pair of new solutions are created around 4.21. An infinity of other pairs of solutions are sequentially created as Γ increases to infinity.

for fiber-based OEOs are not valid anymore in this context. Fiber-based OEO models relied on delay-differential equations, with the de-lay being induced by the fiber delay line, and the dynamics (i. e., the derivative term) was on the microwave variable A(t).

For WGM-based OEOs, we rather have a modal expansion model, in the sense that we have one equation per optical mode, and the dynamics is on the optical modes G n (t).

It is also important to note that this model is nonlinear and continuous. Hence, we can analytically and numerically determine the various dynamical behaviors of the system, and investigate their stability.

Stationary states

The equilibria (or fixed points) of autonomous oscillators are obtained by setting all the derivatives to zero. In the complex envelope formalism, a trivial equilibrium corresponds to the absence of oscillation, while nontrivial equilibria correspond to a steadystate oscillation (because of the amplitudes are constant and not null). The aim of this section is to determine the fixed points of the OEO.

The stationarity conditions Ġn ≡ 0 yield

G n = T E n , ( 20
)
where

T = 2 √ τ e τ d 1 τ + iσ (21)
is the transmission coefficient of the WGM resonator at the drop port. It is interesting to note that we always have |T | ≤ 1, and ideal transmission (T = 1) only occurs when a resonant laser radiation (σ = 0) is coupled to a lossless resonator (τ i → +∞) with coupling photon lifetimes that are matched (τ e = τ d ).

According to [START_REF] Volyanskiy | Compact optoelectronic microwave oscillators using ultra-high Q whispering gallery mode disk-resonators and phase modulation[END_REF] and ( 20), we have

A st = 2βe iς +∞ n=-∞ G n+1 G * n = 2βe iς |T | 2 +∞ n=-∞ E n+1 E * n . ( 22
)
Hence, using [START_REF] Chiasera | Spherical whispering-gallery-mode microresonators[END_REF], the stationary amplitude A st = |A st | obeys

A st = -β sin(2φ) |T | 2 e iς +∞ n=-∞ (-1) n J n+1 (A st )J n (A st ), (23)
and using successively the Bessel equalities

J -n (x) = (-1) n J n (x) (24)
J m (x + y) = +∞ n=-∞ J n (x)J m-n (y), (25)
we are finally led to the transcendental equation

A st = Γ J 1 (2A st ), (26)
where

Γ = -β sin(2φ) |T | 2 e iς (27)
is the overall loop gain, which is essentially the product of the optoelectronic gain and the power transmission factor of the coupled WGM resonator. Existence of a stationary state requires the phase factor e iς to be real, and equal to ±1 such that Γ is real and positive. Effectively, this round-trip phase matching condition corresponds to the necessity of a constructive interference between successive round-trip replicas of the microwave (Barkhausen condition for the phase).

The possible solutions of the transcendental (26) are therefore 26) as Γ is increased. Only the trivial solution exists for Γ < 1. Then, the primary branch is the unique non-trivial solution for 1 < Γ < 15.52. This primary branch of solutions, which exists for all Γ > 1, increases monotonously but has an horizontal asymptote A max 1 st = 1.91. New non-trivial solutions (secondary branch) emerge for Γ > 15.52. The bifurcation analysis shows that the microwave envelope A undergoes a pitchfork bifurcation at Γ = 1, which corresponds to a Hopf bifurcation for the microwave voltage variable V. The same analysis also shows that A undergoes a saddle-node bifurcation at Γ = 15.52 (emergence of two new fixed points, one being stable and the other one unstable). As Γ is further increased, the system undergoes saddle-node bifurcations any time a new branch of paired solutions emerges.

A st = A tr = 0 valid for all Γ A osc = 1 2 Jc -1 1 1 2Γ valid for all Γ ≥ 1 , (28)
where Jc -1 1 is the inverse Bessel-cardinal function. The trivial equilibrium is therefore a solution for all values of Γ , while nontrivial (oscillatory) solutions A osc can only exist for Γ > 1.

As explained in Fig. 2, the stationary solutions of the transcendental equation ( 26) are given by the intersection of the functions J 1 (2A st ) and A st /Γ . When Γ < 1, there is only one solution (the trivial one), while for Γ > 1, both functions will intersect for other points than 0, thereby generating nontrivial equilibria. When 1 < Γ < 15.52, there is only one oscillatory solution. However, at Γ = 15.52, a new pair of nontrivial solutions emerges and coexists with the previous oscillatory state. As Γ increases to infinity, an infinity of branches generating paired solutions are created; it appears very clearly that all of them are converging to the zeros of the Bessel function J 1 (2x) when Γ → +∞.

It is important to note that realistic values for the normalized gain Γ are generally not as high as 15. That is why in all the experiments of OEOs, the oscillator is generally operated in the first branch of non-trivial solutions, whose maximal value is the first zero of J 1 (2A st ), yielding A max 1 st = 1.91. This asymptotic saturation can be observed in Fig. 3. In this paper, we will refer to this branch of nontrivial solution as the primary branch of oscillatory solutions. Figure 3 also displays the emergence of the secondary branch of paired solutions above Γ = 15.52.

After the determination of the various stationary states of the system, we will perform in the next two sections their stability analysis. 51). This function is inferior to 1/4 for the primary branch of nontrivial solutions (starting at Γ = 1). For the secondary branch (starting at Γ = 15.52), one solution is stable (R < 1/4) while the other one is unstable (R > 1/4).

Stability of the trivial equilibrium

The onset of oscillations generally occurs when the trivial fixed point A tr = 0 loses its stability. In this section, we perform the stability analysis of this fixed point in order to determine the conditions leading to oscillations. Let us consider a perturbation δA around the trivial equilibrium A tr = 0. If that perturbation decreases with time, the trivial equilibrium is stable; otherwise, if it increases, the rest point is unstable and oscillations are triggered. The optical spectral components E n excited by the perturbation δA explicitly read

E n n (φ) ×        1 2 n n! (δA) n if n ≥ 0 (-1) -n 2 -n (-n)! (δA -n ) * if n < 0 . ( 29
)
It appears that whenever |n| > 1, the fields E n are of higher order of perturbation, and can therefore be neglected in a linear stability analysis. On the other hand, the induced fields at orders n = 0, ± 1 are explicitly given by

E 1 = - 1 2 sin(φ) δA (30)
E 0 = cos φ (31)

E -1 = - 1 2 sin(φ) δA * , ( 32
)
In other words, the mode n = 0 is of zeroth order and is not influenced by the microwave perturbation δA in the linear approximation, while the modes n = ±1 are of first order and are directly proportional to |δA|. For the sake of mathematical clarity, we will then rewrite the fields E ±1 as δE ±1 since they are first order perturbations. It straightforwardly appears that the modal output variables related to G will have the same order of magnitude as their input counterpart E. Hence, we will have to consider only the variables G 0 and δG ±1 , and neglect all the remaining ones. The thick gray curve is the experimental optical power at the output of the "through" port of the resonator when the input wavelength is rapidly swept. The optical signal at resonance decays exponentially and interferes with the next wavelengths coming from the input laser. The theoretical fit that enables to extract the intrinsic and excitation photon lifetimes is derived in [START_REF] Dumeige | Determination of coupling regime of high-Q resonators and optical gain of highly selective amplifiers[END_REF]. Here, the intrinsic and excitation photon lifetimes are measured at 0.65 µs and 2.9 µs, respectively.

Using [START_REF] Volyanskiy | Compact optoelectronic microwave oscillators using ultra-high Q whispering gallery mode disk-resonators and phase modulation[END_REF] the microwave perturbation can be rewritten as

δA = 2βe iς [G 0 δG * -1 + G * 0 δG 1 ], (33)
while according to [START_REF] Savchenkov | Optical resonators with ten million finesse[END_REF], the output field perturbations G ±1 obey

δ Ġ1 = - 1 τ + iσ δG 1 + 2 √ τ e τ d δE 1 (34) δ Ġ-1 = - 1 τ + iσ δG -1 + 2 √ τ e τ d δE -1 . ( 35
)
In the matrix form, the above equation can be rewritten as

δ Ġ1 δ Ġ * -1 = [S tr ] δG 1 δG * -1 , ( 36
)
where

[S tr ] = 2 √ τ e τ d       1 T Γ 2 -1 Γ 2T * Γ 2T 1 T * Γ 2 -1       ( 37
)
is a 2 × 2 matrix ruling the dynamics of the perturbation flow (Jacobian). The trivial fixed point will be stable whenever the eigenvalues of this Jacobian matrix have strictly negative real parts. The determination of these eigenvalues is straightforward and it is found that the trivial equilibrium is stable only when Γ < 1, and unstable otherwise. We will show in the next section that crossing that threshold value triggers unconditionally stable microwave oscillations in the primary branch of nontrivial equilibria, while for the higher order branches of paired solutions, some solutions are stable while others are not.

Stability of the oscillatory solution

So far, the theoretical analysis has shown that there are two types of stationary solutions. The trivial equilibrium exists for all Γ but is stable only for Γ < 1. The oscillatory solutions only exist for Γ > 1, and the purpose of this section is to demonstrate that they might be stable or unstable. Once again, to demonstrate that this solution is stable, we have to show that any perturbation δA of the oscillatory solution A osc of interest exponentially decays to zero. Otherwise, the oscillatory solution is unstable.

According to [START_REF] Volyanskiy | Compact optoelectronic microwave oscillators using ultra-high Q whispering gallery mode disk-resonators and phase modulation[END_REF], the perturbation of the steady state microwave solution A osc obeys

δA = 2βe iς +∞ n=-∞ G n+1 δG * n + G * n δG n+1 . ( 38
)
On the other hand, the steady-state input and output electric fields obey

E n = (φ) J n (A osc ) (39)
G n = T E n . ( 40
)
In the demonstration we have made to investigate the stability of the trivial fixed point, it appeared that the Jacobian matrix [S tr ] had to be expressed relatively to the The microwave noise from the amplifier is visible as a plateau on each sidemode, on top of which the oscillation at the FSR frequency stands. It is apparent on this figure that the amplitudes of the sidemodes decreases very rapidly, so that only few of them are necessary to describe the behavior of the oscillator in the numerical simulations.

variables δG 1 and δG * -1 . In the case of the nontrivial solutions, (38) shows that we have an infinity of perturbations to consider, but however, we will decompose by analogy the microwave perturbation as

δA = δA + + δA - (41)
where

δA + = 2βe iς +∞ n=-∞ T * E * n-1 δG n (42) δA -= 2βe iς +∞ n=-∞ T E 1-n δG * -n (43)
are global variables associated to the output field perturbations δG n and their counterparts δG * -n (complex conjugate, opposite sidemode), respectively. We will hereafter use these variables to obtain a Jacobian matrix whose eigenvalues will decide the stability of the nontrivial stationary states.

In the feedback loop, the perturbation δA will first induce perturbations δE n . The first-order Taylor expansion of a perturbation of the amplitude can be determined as

|A osc + δA| A osc + 1 2 [δA + δA * ]. (44)
Hence, using (39), the input field perturbations can be calculated as

δE n = 1 2 n (φ)J' n (A osc ) [δA + δA * ], (45)
where the prime denotes the derivative of the Bessel function relatively to its argument A osc . The input field δE n do induce output field perturbations δG n , which obey

δ Ġn = - 1 τ + iσ δG n + 2 √ τ e τ d δE n (46)
By multiplying the aforementioned equation by 2βe iς T * E * n-1 and summing over all modal indices n, we are led to the following equation for δA + :

δ Ȧ+ = - 1 τ + iσ δA + + 4βe iς T * √ τ e τ d +∞ n=-∞ E * n-1 δE n , (47)
and analogously, it can be found that δA -obeys [ns]

δ Ȧ-= - 1 τ -iσ δA -+ 4βe iς T √ τ e τ d +∞ n=-∞ E 1-n δE * -n . ( 48
)
Figure 8: Experimental time-domain dynamics. The (light gray) experimental trace of the microwave signal is obtained with a ultra-fast oscilloscope, just after abruptly switching-on the laser. At the very beginning, the signal consists of noise from the RF amplifier. The RF oscillation at the FSR frequency rapidly grows above the noise level, and its amplitude reaches the stationary state with a time-scale in the µs range. The inset is a zoom-in presenting the fastscale dynamics of the same curve. It shows that the microwave frequency is indeed around 10 GHz. The apparent amplitude noise is an artefact originating from the fact that albeit ultra-fast, the 30 GHz sampling rate only provides 6 points for every period of our 10 GHz signal. The experimental microwave envelope signal (thick blue line) is in very good agreement with the envelopes A(t) obtained numerically in Fig. 9.

We demonstrate in the Appendix that the perturbation (47) and (48) can be rewritten as

δ Ȧ+ = - 1 τ + iσ δA + -R[δA + δA * ] (49) δ Ȧ-= - 1 τ -iσ δA --R[δA + δA * ] , ( 50
)
where R is a function of the gain Γ :

R = 1 4 Γ [J 0 (2A osc ) -J 2 (2A osc )] . ( 51
)
Since δA = δA + + δA -, we can finally rewrite (49) and (50) under the form of a four-dimensional autonomous flow

             δ Ȧ+ δ Ȧ- δ Ȧ * + δ Ȧ * -              = [S osc ]              δA + δA - δA * + δA * -              , ( 52
)
where the Jacobian matrix [S osc ] can be written under the synthetic block matrix form

[S osc ] = [U] [V] [V * ] [U * ] (53) with [U] = 2 √ τ e τ d R-1 T R T R T * R-1 T * (54) [V] = 2 √ τ e τ d R T R T R T * R T * . ( 55
)
The oscillatory solution A osc is therefore stable if all the eigenvalues of the constant and complex-valued matrix [S osc ] have a strictly negative real part. A straightforward method would be to actually compute the eigenvalues of this four-dimensional matrix and evaluate the sign of their respective real parts. However, we can circumvent that tedious task by noting due to the particular structure of [S osc ], the set of its eigenvalues is the union of the eigenvalues of the two-dimensional matrices [U + V] and [U -V]. Hence, the stability of the oscillations is guaranteed as long as the eigenvalues of the 2 × 2 matrices [U ± V] have strictly negative real parts. This stability condition is trivially satisfied for [U -V] which is a diagonal matrix whose diagonal elements (i. e., eigenvalues) are -1/τ ± iσ. On the other hand, as far as the matrix [U + V] is concerned, the above conditions respectively yield R < 1 4 . According to Fig. 4, this inequality is indeed satisfied whenever Γ > 1 for the branch of primary solutions, since R has an absolute maximum value equal to 1 4 for Γ = 1, and decreases monotonously afterwards. On the other hand, the higher-order branches of paired solutions have the typical stability pattern of saddle-node fixed points, as one solution remains stable while the other is unstable.

We recall again that experimentally, it is extremely difficult to reach Γ values of the order of 15. Most experimental studies can in fact hardly achieve gain values superior to 3. Therefore, the solution of practical interest belongs to the primary branch: our analysis has demonstrated that this oscillating solution is unique in the range 1 < Γ < 15.52, and is always exponentially stable for any Γ > 1. It is interesting to note that in the case of fiber-based OEOs, we proved that the microwave oscillation was stable only at up to Γ = 2.31, in agreement with experiments [START_REF] Chembo | Dynamic instabilities of microwaves generated with optoelectronic oscillators[END_REF][START_REF] Chembo | Nonlinear dynamics and spectral stability of optoelectronic microwave oscillators[END_REF]. Here, we prove that the oscillator is unconditionally stable at up to a much higher value (15.52), thereby demonstrating that this WGM-based OEO is significantly more stable than its fiberbased counterpart.

Comparison between numerical and experimental results

The experimental setup is presented in detail in Fig. 1. A preliminary measurement is the evaluation of the intrinsic and coupling Q factors of the cavity. This measure, presented in Fig. 5 is performed in the open loop configuration and is mathematically explained in ref. [START_REF] Dumeige | Determination of coupling regime of high-Q resonators and optical gain of highly selective amplifiers[END_REF]. It enables to confirm that all our quality factors are of the order of 10 8 , yielding characteristic time scales of the order of 1 µs. When the oscillation loop is closed, we can monitor different optical and microwave variables of interest.

The optical output of the coupled resonator's through-port is used to control the detuning σ between the input laser and the resonance. Once the detuning is set and the gain is above threshold, oscillations are sustained in the optical branch. Figure 6 shows the optical spectrum taken at the "through" output of the resonator (as drawn in Fig. 1). The spectral line due to the laser is in the center at f L = 193 276.6 GHz, and two pairs of sidemodes are visible in this case, separated by the FSR f M = 10.4 GHz. As we can see on this 75 dB dynamic-range figure, the amplitude of these sidemodes decreases very rapidly decreases with their order |n|.

In the RF branch, the spectrum of the microwave oscillation is measured using a 22 GHz RF spectrum analyzer. A typical RF spectrum is presented in Fig. 7 and demonstrates that the oscillation arises at a fixed frequency given by the FSR f M of the resonator. A plateau is visible below the oscillation spectral line and is due to the noise of the RF amplifier filtered between 10 and 11 GHz.

Simultaneously, the temporal dynamics of the microwave signal is monitored using a 40 GHz-bandwidth oscilloscope. To investigate this time-domain dynamics, another MZ modulator driven by a square signal is used to switch the input laser light on and off. Therefore, the gain Γ of the oscillator is abruptly changed from 0 to a value that is higher than 1, and the transient dynamics can be monitored using a photodiode and a fast oscilloscope. The resulting signal is displayed in Fig. 8. The actual signal V(t) is shown in light gray while the envelope V (t) is the thick blue curve. The inset of this figure is provided to show the oscillation occurring at a faster time scale (in the range of 100 ps, which corresponds to the period T M = 1/f M ), while the envelope time scale is in the µs range. It is worth noting that the apparent amplitude noise in this inset is an artefact due to the barely sufficient sampling rate of the oscilloscope (30 GHz) compared to the oscillation frequency (10 GHz).

The envelope curve of Fig. 8 can be compared to the numerical simulations presented on Fig. 9. These simulations were obtained using a fourth-order Runge-Kutta algorithm for the [START_REF] Chiasera | Spherical whispering-gallery-mode microresonators[END_REF], [START_REF] Savchenkov | Optical resonators with ten million finesse[END_REF] and [START_REF] Volyanskiy | Compact optoelectronic microwave oscillators using ultra-high Q whispering gallery mode disk-resonators and phase modulation[END_REF]. The initial conditions where random complex value for A, taken in a gaussian distribution centered around 0.01 with standard deviation of 0.01. Both numerical and experimental curves feature the same characteristics, and the transition from the initial state to the stationary solution occurs in similar duration, of the order of 1 µs. This very good agreement between the simulation and the experimental result validates the experimental interest of this model.

It is interesting to note that from a purely theoretical point of view, the model is infinite dimensional because there is an infinity of fields G n to consider. However, only a few of them are necessary to yield accurate results, as foreshadowed by the experimental spectrum of Fig. 6. In practice, our simulations were performed with 10 pairs of sidemodes, even though considering 5 pairs would have already been very accurate.

Conclusion

In this article, we have proposed a nonlinear dynamics approach to study WGM-Based OEOs. We have used a complex microwave envelope variable to investigate the timedomain behavior and stability properties of this oscillator. Our study has enabled to determine the various stationary states and their stability. It was shown that above threshold, the principal branch of the oscillations, which is the only one experimentally accessible so far, is always exponentially stable regardless of the gain. However, the analysis has also evidenced higher order branches of solutions whose stability properties are more complex, with some states being stable while the others are not. Both the analytical and numerical analysis have been confirmed by the experimental measurements.

Future work will consist in using the model to optimize the metrics of the oscillator. In particular, we aim to investigate the phase noise performance of this WGM-based OEO by adding calibrated noise terms in our dynamical equations. We would then obtain stochastic differential equations that would enable to predict the phase noise the spectra and the Allan deviation of the oscillator.

It is already known that very high microwave frequencies, at the edge of the millimeterwaves band ( 100 GHz), can be generated by selecting a higher harmonic of the beatnote signal detected by the photodetector. Our model enables to analyze the microwave envelope of such harmonics kΩ M by summing the quadratic terms G n+k G * n in [START_REF] Ilchenko | Optical Resonators With Whispering-Gallery Modes-Part II: Applications[END_REF]. Another prospective work would be to extend this formalism in order to account for nonlinear phenomena in the WGM resonator. Some research works (see e. g. [START_REF] Saleh | Optical scattering noise in high Q fiber ring resonators and its effect on optoelectronic oscillator phase noise[END_REF]) have already demonstrated that scattering has a measurable effect on the phase noise performance of OEOs, particularly when the optical cavity has a high Q factor. Various nonlinear phenomena (Kerr, Raman, Brillouin) have an effect on the phase noise performance, as well as the chromatic dispersion which converts laser frequency noise into phase fluctuations. We expect this theoretical approach to be able to give both a quantitative and qualitative insight into all into these phenomenologies. Since -n (φ) = (-1) n n (φ), it can also be shown that in the non-trivial stationary states, E * -n = E n and δE * -n = δE n , so that the infinite sums in the right-hand sides of (47) and (48) are identical. Then, using [START_REF] Chembo | Generation of Ultralow Jitter Optical Pulses Using Optoelectronic Oscillators With Time-Lens Soliton-Assisted Compression[END_REF] and (61) the demonstration of (49) and (50) is straightforward.

Figure 1 :

 1 PhotodetectorBandpass filter

4 Γ 4 =

 44 15.52

Figure 3 :

 3 Figure3: Variation of the non-trivial solutions of (26) as Γ is increased. Only the trivial solution exists for Γ < 1. Then, the primary branch is the unique non-trivial solution for 1 < Γ < 15.52. This primary branch of solutions, which exists for all Γ > 1, increases monotonously but has an horizontal asymptote A max 1

Figure 4 :

 4 Figure 4: Variations of the function R expressed in (51). This function is inferior to 1/4 for the primary branch of nontrivial solutions (starting at Γ = 1). For the secondary branch (starting at Γ = 15.52), one solution is stable (R < 1/4) while the other one is unstable (R > 1/4).

Figure 5 :

 5 Figure5: Cavity ring down measurement. The thick gray curve is the experimental optical power at the output of the "through" port of the resonator when the input wavelength is rapidly swept. The optical signal at resonance decays exponentially and interferes with the next wavelengths coming from the input laser. The theoretical fit that enables to extract the intrinsic and excitation photon lifetimes is derived in[START_REF] Dumeige | Determination of coupling regime of high-Q resonators and optical gain of highly selective amplifiers[END_REF]. Here, the intrinsic and excitation photon lifetimes are measured at 0.65 µs and 2.9 µs, respectively.

Figure 6 :

 6 Figure 6: Optical spectrum of the signal at the through-port of the resonator. The central frequency f L = 193 276.6 GHz is the optical frequency of the input laser. The sidemodes created by the MZ intensity modulator driven by the microwave oscillation are clearly visible on each side of the input frequency.The microwave noise from the amplifier is visible as a plateau on each sidemode, on top of which the oscillation at the FSR frequency stands. It is apparent on this figure that the amplitudes of the sidemodes decreases very rapidly, so that only few of them are necessary to describe the behavior of the oscillator in the numerical simulations.

Figure 7 :

 7 Figure 7: RF spectrum of the generated microwave. The 10.41 GHz oscillation corresponding to the FSR of the resonator. This microwave signal is strong, and stands 50 dB above the filtered noise of the RF amplifier.

Figure 9 :

 9 Figure 9: Numerical simulation of the microwave envelope dynamics, and comparison with experiments. The continuous black (down) and dashed blue (top) curves are numerical simulations obtained for Γ = 1.18 and 1.55, respectively. The other parameters are τ i = 0.65 µs, τ e = 2.9 µs (same as the experimental values), τ d = 0.1 µs, σ = 0 (at resonance), and φ = π/4. The initial condition is a random value taken in a gaussian distribution with average value and standard deviation of 0.01. The envelopes of experimental time traces are plotted in thick lighter lines on top of the continuous black and dashed blue curves (the photodiode and artefact noise have been substracted).

  ) n-1 (φ) J n-1 (A osc ) × [J n-1 (A osc ) -J n+1 (A osc )]However, we have n (φ) n-1 (φ) = (-1) n sin φ cos φ, while (24) and (25) yield the following Bessel relationships:n J n-1 (x)J n+1 (x) = -+∞ n=-∞ J 1-n (x)J n+1 (x) = -J 2 (2x). φ (δA + δA * ) × [J 0 (2A osc ) -J 2 (2A osc )].(61)

Acknowledgements

The authors would like to thank Prof. Patrice Féron for insightful discussions about the theoretical and experimental characterization of WGM resonators.

† The authors acknowledge financial support from the Centre National d'Etudes Spatiales through the project SHYRO (Action R & T-R-S10/LN-0001-004/DA:10076201), from the ANR project ORA (BLAN 031202), and from the Région de Franche-Comté. They also acknowledge financial support from the European Research Council through the project NextPhase (ERC StG 278616).

Appendix A. Bessel-cardinal functions

We define the Bessel-cardinal function of order n as

where J n is the n-th order Bessel function of the first kind. From a qualitative point of view, Bessel-cardinal functions look like the sine-cardinal function sinc(x) = sin x/x when n > 0., with an absolute maximum centered at x = 0, and an oscillatory behavior converging to zero as x → ±∞. On the other hand, the Bessel-cardinal function diverges to infinity as x → ±∞ with an oscillatory behavior when n < 0. Since J n (r)e inθ can be rewitten under the analytical form z n Jc n (|z|) with z = re iθ , the Bessel-cardinal formalism is very useful to carry out some of the mathematical calculations.

B. Demonstration of (49) and (50)

This demonstration relies on the explicit calculation of the infinite sums in the righthand side of (47) and (48). Let's fist calculate the sum in (49). Using (39) and (44) and the recurrence relationship