ABSTRACT

Title of Dissertation: QUANTUM AND STOCHASTIC DYNAMICS OF
KERR MICROCOMBS

Fengyu Liu
Doctor of Philosophy, 2024

Dissertation Directed by: Professor Yanne K. Chembo
Department of Electrical and Computer Engineering

Kerr microcombs are sets of discrete, equidistant spectral lines and are typically generated

by pumping a high-quality factor optical resonator with a continuous- or pulse-wave resonant

laser. They have emerged as one of the most important research topics in photonics nowadays,

with applications related to spectroscopy, sensing, aerospace, and communication engineering.

A key characteristic of these microcombs is the threshold pump power. Below the threshold, two

pump photons are symmetrically up- and down-converted as twin photons via spontaneous four-

wave mixing, and they can be entangled across up to a hundred eigenmodes. These chipscale,

high-dimensional, and room-temperature systems are expected to play a major role in quantum

engineering. Above the threshold, the four-wave mixing process is stimulated, ultimately leading

to the formation of various types of patterns in the spatio-temporal domain, which can be ex-

tended (such as roll patterns), or localized (bright or dark solitons). The semiclassical dynamics

of Kerr microcombs have been studied extensively in the last ten years and the deterministic char-

acteristics are well understood. However, the quantum dynamics of the twin-photon generation



process, and the stochastic dynamics led by the noise-driven fluctuations, are still not so clear.

In the first part of our investigation, we introduce the theoretical framework to study the

semiclassical dynamics of the Kerr microcombs based on the slowly varying envelope of the

intracavity electrical fields. Two equivalent models – the coupled-mode model and the Lugiato-

Lefever model are used to analyze the spectro- and spatio-temporal dynamics, respectively. These

models can determine the impact of key parameters on the Kerr microcomb generation process,

such as detuning, losses, and pump power, as well as critical values of the system, such as thresh-

old power. Various types of patterns and combs can be observed through simulations that follow

experimental parameters. Furthermore, we show an eigenvalue analysis method to determine

the stability of the microcomb, and this method is applied to an unstable microcomb solution to

understand the generation of subcombs surrounding the primary comb.

In the second and third parts, we investigate a stochastic model where noise is added to the

coupled-mode equations governing the microcomb dynamics to monitor the influence of random

noise on the comb dynamics. We find the model with additive Gaussian white noise allows us to

characterize the noise-induced broadening of spectral lines and permits us to determine the phase

noise spectra of the microwaves generated via comb photodetection. Our analysis indicates that

the low-frequency part of the phase noise spectra is dominated by pattern drift while the high-

frequency part is dominated by pattern deformation. The dynamics of the Kerr microcomb with

multiplicative noises, including thermal and photothermal fluctuations, are also investigated in the

end. We propose that the dynamics of the noise can be included in the simulation of stochastic

dynamics equations, introduce the methods to solve the dynamics of the noise, and study a quiet

point method for phase noise reduction.

In the fourth part, we use canonical quantization to obtain the quantum dynamics for Kerr



microcombs generated by spontaneous four-wave mixing below the threshold and develop the

study of them using frequency-bin quantum states. We introduce a method to find the quantum

expansion of the output state and explore the properties of the eigenkets. A theoretical framework

is also developed to obtain explicit solutions for density operators of quantum microcombs, which

allows us to obtain their complete characterization, as well as for the analytical determination of

various performance metrics such as fidelity, purity, and entropy. Finally, we describe a quantum

Kerr microcomb generator with a pulse-wave laser and propose the time-bin entangled states

generated by it.



QUANTUM AND STOCHASTIC DYNAMICS OF KERR MICROCOMBS

by

Fengyu Liu

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy

2024

Advisory Committee:
Professor Yanne K. Chembo, Chair/Advisor
Professor Cheng Gong
Professor Mario Dagenais
Professor Avik Dutt
Professor Rajarshi Roy, Dean’s Representative



© Copyright by
Fengyu Liu

2024



Dedication

To my parents, my grandparents, and my fiancée.

ii



Acknowledgments

Completing this dissertation marks the end of an unforgettable graduate student journey

filled with challenges, growth, and valuable experiences. I express my profound gratitude to

those whose support and guidance made this endeavor possible and shaped it into a cherished

chapter in my life.

First and foremost, I extend my deepest appreciation to my advisor, Professor Yanne K.

Chembo, whose unwavering support and guidance have been instrumental throughout this jour-

ney. His willingness to share his expertise, offer insightful advice, and provide opportunities

for engaging in challenging and fascinating projects have been truly invaluable. Working along-

side such an extraordinary mentor has been an honor and a privilege. Without his guidance and

expertise, realizing this dream would have remained beyond reach.

My heartfelt thanks extend to Professor Cheng Gong, Professor Mario Dagenais, Profes-

sor Avik Dutt, and Professor Rajarshi Roy for graciously agreeing to serve on my dissertation

committee and for their invaluable feedback on the manuscript.

The enriching environment of the Photonics Systems Lab for Aerospace and Communica-

tion Engineering (PSACE Lab) was critical to my graduate experience. I am deeply grateful to

my colleagues, including Dr. Elham Heidari, Dr. Richard Brewster, Meenwook Ha, and Haoying

Dai, whose contributions and many valuable discussions have significantly enriched my research

journey. I also need to thank Dr. Helene Nguewou-Hyousse, Sadia Siraz, Benjamin Klimko,

iii



Jazzmin Robinson, Millicent Ayako, and Beck Saunders for their support and help in the group.

I am thankful for the fruitful interactions with my collaborators and co-authors, including

Professor Curtis R. Menyuk, Professor Guoping Lin, Professor Aurélien Coillet, and Dr. Damià

Gomila. I would also like to express my gratitude to all the faculty who have taught me at

the University of Maryland, College Park. Acknowledgment is also due to the dedicated staff

members in the Department of Electrical and Computer Engineering (ECE) and the Institute for

Research in Electronics and Applied Physics (IREAP) for their support during these five years.

Furthermore, I wish to express my appreciation to the Air Force Office of Scientific Re-

search (AFOSR) and the U.S. Department of Energy (DoE) for their financial support, which has

been instrumental in realizing the projects discussed in this dissertation.

My deepest thanks go to my family for their unconditional support and guidance throughout

my Ph.D. journey. Their steadfast encouragement has been a source of strength for me during

challenging times. I am immensely grateful to my parents, whose unwavering belief in me has

been a constant source of inspiration. I also extend my heartfelt gratitude to my fiancée Xinyuan

Chen, whose unwavering support has been a cornerstone of my journey.

Lastly, I extend my sincere gratitude to all for guiding me through this transformative

journey.

iv



Table of Contents

Dedication ii

Acknowledgements iii

Table of Contents v

List of Figures viii

List of Abbreviations xii

Chapter 1: Introduction 1
1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Kerr Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Whispering-Gallery Mode Resonators . . . . . . . . . . . . . . . . . . . . . . . 6

Chapter 2: Semiclassical Dynamics of Kerr Microcombs 9
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 System and Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Eigenfrequences and Eigenmodes . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4 Semiclassical Dynamics – Coupled-Mode Equations . . . . . . . . . . . . . . . 16

2.4.1 Coupled-Mode Equations with Modal Fields in the Normalization of Pho-
ton Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.4.2 Coupled-Mode Equations with Modal Fields in the Normalization of In-
tracavity Power . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.5 Semiclassical Dynamics – Lugiato-Lefever Equation . . . . . . . . . . . . . . . 18
2.5.1 Lugiato-Lefever Equation with Intracavity Field in the Normalization of

Photon Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.5.2 Lugiato-Lefever Equation with Intracavity Field in the Normalization of

Intracavity Power . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5.3 Lugiato-Lefever Equation with Intracavity Field in the Dimensionless

Time Normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.6 Numerical Simulation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.7 Spatio-Temporal and Spectro-Temporal Solutions . . . . . . . . . . . . . . . . . 22
2.8 Eigenvalue Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.8.1 Eigenvalue Analysis of Solitons . . . . . . . . . . . . . . . . . . . . . . 24
2.8.2 Eigenvalue Analysis of Rolls . . . . . . . . . . . . . . . . . . . . . . . . 26
2.8.3 Subharmonic Instabilities of Rolls . . . . . . . . . . . . . . . . . . . . . 27

v



2.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Chapter 3: Stochastic Dynamics and Phase Noise of Kerr Microcombs
– Additive Gaussian White Noise 31

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Stochastic Model for Additive Gaussian White Noise . . . . . . . . . . . . . . . 32
3.3 Optical Fluctuations Under Threshold . . . . . . . . . . . . . . . . . . . . . . . 34
3.4 Microcomb Above Threshold: Fluctuations from Pattern Drift and Pattern Defor-

mation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.4.1 Optical Stochastic Fluctuations of Solitons . . . . . . . . . . . . . . . . . 39
3.4.2 Optical Stochastic Fluctuations of Rolls . . . . . . . . . . . . . . . . . . 46

3.5 Optical Power Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.5.1 Contribution of Pattern Drift . . . . . . . . . . . . . . . . . . . . . . . . 54
3.5.2 Contribution of Pattern Deformation . . . . . . . . . . . . . . . . . . . . 56
3.5.3 Comparison Between Numerical and Analytical Results . . . . . . . . . 58

3.6 Microwave Phase Noise Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.6.1 Contribution of Pattern Drift . . . . . . . . . . . . . . . . . . . . . . . . 60
3.6.2 Contribution of Pattern Deformation . . . . . . . . . . . . . . . . . . . . 62
3.6.3 Comparison Between Numerical and Analytical Results . . . . . . . . . 63

3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

Chapter 4: Stochastic Dynamics and Phase Noise of Kerr Microcombs
– Multiplicative Noise 66

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.2 Stochastic Model for Multiplicative Noises . . . . . . . . . . . . . . . . . . . . . 66
4.3 Sources of Multiplicative Noises . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.3.1 Fluctuations of Refractive Index . . . . . . . . . . . . . . . . . . . . . . 68
4.3.2 Fluctuations of Radius . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.4 Thermal Fluctuations in Microresonators . . . . . . . . . . . . . . . . . . . . . . 70
4.4.1 Thermorefractive Fluctuations . . . . . . . . . . . . . . . . . . . . . . . 72
4.4.2 Thermoelastic and Thermal Expansion Fluctuations . . . . . . . . . . . . 75
4.4.3 Photothermal Fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . 78

4.5 Quiet Point Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

Chapter 5: Quantum Dynamics and Representations of Kerr Microcombs 82
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.2 Heisenberg Picture – Quantum Langevin Equations . . . . . . . . . . . . . . . . 82

5.2.1 Canonical Quantization of Coupled-Mode Equations . . . . . . . . . . . 83
5.2.2 Dynamics of Modal Fields with the Hamiltonian . . . . . . . . . . . . . 85
5.2.3 Canonical Quantization of Lugiato-Lefever Equations . . . . . . . . . . . 86

5.3 Quantum Kerr Microcombs Generated Under Threshold . . . . . . . . . . . . . . 87
5.3.1 Quantum Kerr Combs in Frequency Domain . . . . . . . . . . . . . . . . 90
5.3.2 Spectra of the Quantum Fields . . . . . . . . . . . . . . . . . . . . . . . 92
5.3.3 Output Photon Flux and Power . . . . . . . . . . . . . . . . . . . . . . . 94

vi



5.3.4 Second-order Correlations of Quantum Fields . . . . . . . . . . . . . . . 95
5.4 Schrödinger Picture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

5.4.1 The Lossless Configuration . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.4.2 Entanlged Biphoton State Kets . . . . . . . . . . . . . . . . . . . . . . . 104
5.4.3 The Lossy Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . 107

5.5 Time-Bin Entangled States . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

Chapter 6: Conclusions and Perspectives 116

Bibliography 119

vii



List of Figures

1.1 Schematic illustration of the experimental setup for Kerr microcombs genera-
tion based on a 12 mm magnesium fluoride (MgF2) disk resonator at 1550 nm.
WGMR: Whispering-gallery mode resonator; PC: Fiber polarization controller;
EDFA: Erbium-doped fiber amplifier; VOA: Variable fiber optical attenuator;
FC3dB: 1×2 50:50 fiber coupler; PD: Photodetector. Source: Ref. [8]. . . . . . . 2

1.2 (a) Ray optics viewpoint and (b) Wave optics viewpoint of the whispering-gallery
mode in a WGMR. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.3 Schematic diagrams of microresonators. From left to right: Whispering-gallery
mode-based bulk toroids, spheres, monolithic toroids, and integrated ring-resonators.
Source: Ref. [19]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.4 The coupling methods for the WGM resonators with optical waveguides: (a)
Prism coupling, (b) Waveguide side coupling, (c) Waveguide tip coupling. Source:
Ref. [40]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.1 Two main coupling configurations for Kerr comb generation. (a) Add-through
coupling. (b) Add-drop coupling. Source: Ref. [41]. . . . . . . . . . . . . . . . . 11

2.2 Distributions of electric fields of WGMs in a spherical WGM resonator. (a) Side
view (left) and top view (right) of the mode with n = 1, ℓ = 30, and ℓ − |m| =
0. (b) Examples of transverse field distributions for different radial and polar
eigennumbers n andm. Each transverse field distribution corresponds to a family
of modes with different ℓ. Source: Ref. [18]. . . . . . . . . . . . . . . . . . . . . 13

2.3 Schematic representation of stationary Kerr microcombs. The microcombs are in
space [(a)–(e)] and frequency [(f)–(j)], after simulation of Eq. (2.23). The upper
row displays the total intracavity field power |E(θ)|2, while the lower row displays
the corresponding stem plot for the modal intensities |Eout,l|2 in a logarithmic
scale (note that they are spectrotemporal snapshots and not Fourier spectra). (a)
and (f) Flat state (no comb). (b) and (g) Bright soliton (P

L
= 4.00 mW, σ = −2κ,

and ζ2 = 2π × 2.92 kHz). (c) and (h) Dark soliton (P
L
= 5.30 mW, σ = −2.5κ,

and ζ2 = −2π × 2.92 kHz). (d) and (i) Roll pattern of order L = 1 (P
L
= 3.09

mW, σ = −κ, and ζ2 = 2π × 874.7 kHz). (e) and (j) Roll pattern of order L = 4
(P

L
= 2.27 mW, σ = −κ, and ζ2 = 2π × 366.4 kHz). . . . . . . . . . . . . . . . 23

viii



2.4 [Left column] Five types of Kerr comb spectra snapshots of the output fields. The
quality factor of the through port is set to Qt = Q/4, and the nonlinear param-
eter is set to γ = 1.0 W−1km−1. From top to bottom, the detuning σ of each
case are 2κ, 2.3κ, 1.8κ, 1.83κ, and 0.55κ, respectively, and the 2nd order dis-
persion ζ2 are set to −0.0037κ/vgω

2
FSR, −0.0073κ/vgω

2
FSR, −0.0060κ/vgω

2
FSR,

−0.0015κ/vgω
2
FSR, and −0.0055κ/vgω

2
FSR. Pump powers PL are 75 mW, 90 mW,

100 mW, 100 mW, and 75 mW, corresponding to the thresholds: 48.4 mW, 57.5
mW, 42.8 mW, 43.6 mW, and 16.5 mW. [Center column] Eigenvalues of the solu-
tion in a complex plane. [Right column] The real part of the eigenvalues is plotted
as a function of the wavenumber q of the Bloch modes. The solution is unstable
when the real part of one of them becomes positive, and the eigenvalue with the
largest positive real part corresponds to the order of the subcomb. . . . . . . . . 29

3.1 Schematic illustration of microwave generation using a Kerr optical frequency
comb. WGMR: Whispering-gallery mode resonator; PD: Photodiode. . . . . . . 31

3.2 Schematic representation of stochastic Kerr microcombs. The combs are in space
[(a)–(e)] and frequency [(f)–(j)], and are at a given time t (snapshot), after simu-
lation of Eq. (3.1) with the same parameters as Fig. 2.3. The upper row displays
the total intracavity field power |E(θ)|2, while the lower row displays the corre-
sponding stem plot for the modal intensities |Eout,l|2 in a logarithmic scale (note
that they are spectrotemporal snapshots and not Fourier spectra). The noise in the
spatiotemporal patterns has been accentuated for the sake of visual clarity. (a)
and (f) Flat state (no comb). (b) and (g) Bright soliton. (c) and (h) Dark soliton.
(d) and (i) Roll pattern of order L = 1. (e) and (j) Roll pattern of order L = 4.
The numerical simulations will involve the four cases from Comb 1 to Comb 4. . 34

3.3 Several computation and simulation results of a Kerr microcomb below the thresh-
old. The parameters are set to: P

L
= 1.65 mW, σ = −κ, and ζ2 = 2π × 2.9 kHz.

(a) Snapshot of the frequency comb pattern under noise. (b) Spectral density of
the sidemode l = 18 obtained after computation of Eq. (3.16) (red line) and sim-
ulation of Eq. (3.1) (blue dots). (c) Average square amplitude of the fluctuations
in each mode: The red line shows the computation results from Eq. (3.17), while
the blue dots show the simulation results from Eq. (3.1). . . . . . . . . . . . . . . 37

3.4 Schematic illustration of pattern drift and pattern deformation. (a) Schematic il-
lustration of pattern drift for a soliton. The red curve represents the initial unper-
turbed soliton, while the blue curves represent the perturbed solitons with pattern
drift. (b) Schematic illustration of pattern deformation for a soliton. The red
curve represents the initial unperturbed soliton, while the green curves represent
the perturbed solitons with pattern deformation. Pattern drift corresponds to a
translational change of the soliton’s azimuthal position within the cavity. On the
other hand, pattern deformation corresponds to a change of the soliton’s shape.
Our analysis shows that pattern drift is the leading contribution to phase noise at
a low frequency offset from the carrier, while pattern deformation is the leading
contribution to phase noise at high frequency offset. . . . . . . . . . . . . . . . 40

ix



3.5 Variance of orbital deviation of the intra-cavity field. The figure shows the results
of (a) bright soliton (Comb 1), (b) dark soliton (Comb 2), (c) roll pattern of order
L = 1 (Comb 3), and (d) roll pattern of order L = 4 (Comb 4) displayed in
Fig. 3.2 as a function of the sidemode order l, respectively. The thick gray curves
correspond to analytical results obtained from Eqs. (3.41), (3.64), and (3.65),
while the thin red curves result from the simulations of Eq. (3.1). One can note
the excellent agreement between analytical and numerical results. . . . . . . . . . 46

3.6 Sidemode spectra of the output optical field. The figure corresponds to (a, e)
bright soliton (Comb 1), (b, f) dark soliton (Comb 2), (c, g) roll pattern of order
L = 1 (Comb 3) and (d, h) roll pattern of order L = 4 (Comb 4) displayed in
Fig. 3.2. (a-d) Spectra of excited sidemodes that are closest to the right of the
center mode. (e-f) Spectra of excited sidemodes that are fourth closest to the
right of the center mode. The blue lines correspond to the noise contribution of
pattern drift [the computation of Eq. (3.74)]. The green lines correspond to the
noise contribution of pattern deformation [the computation of Eq. (3.80)]. The
red lines correspond to the spectra obtained after the numerical simulation of the
stochastic modal equations. These simulations clearly indicate that pattern drift
noise dominates the spectra for low offset frequencies, while pattern deformation
noise dominates the spectra for high offset frequencies. . . . . . . . . . . . . . . 59

3.7 Phase noise spectra of the microwaves. The figure corresponds to the results ob-
tained via photodetection of (a, e) bright soliton (Comb 1), (b, f) dark soliton
(Comb 2), (c, g) roll pattern of order L = 1 (Comb 3) and (d, h) roll pattern of
order L = 4 (Comb 4) displayed in Fig. 3.2. (a-d) Phase noise spectra of the mi-
crowaves with the lowest frequency. (e-f) Phase noise spectra of the microwaves
with the fourth lowest frequency. The simulation results (redline) are obtained
from the simulation of Eq. (3.1). The blue lines correspond to pattern drift noise
[the computation of Eq. (3.89)] while the green lines correspond to pattern defor-
mation noise [the computation of Eq. (3.94)]. The numerical simulations show
that, as in the case of optical noise, pattern drift noise dominates the spectra for
low offset frequencies, while pattern deformation noise dominates the spectra for
high offset frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

5.1 Schematic of the system under study. SFWM creates twin-photons in a fountain-
like fashion via the photonic interaction 2ℏω0 → ℏωk + ℏω−k, where two pump
photons are annihilated in the mode k = 0 and create a biphoton in the sidemodes
of azimuthal eigennumbers ±k. Up to a hundred modes can be thereby entangled
as

∑
k ck| − k⟩|k⟩. However, the microcomb is partially incoherent because it is

also entangled to its environment through the loss coefficient 1− η. . . . . . . . 83
5.2 The spectral density (first column), real part (second column), imaginary part

(third column), and phase (fourth column) of the spectrum for different side-
modes. The parameters are set to g0|A0|2 = κ/10, ζ2 = κ/20, κe = 0.7κ
and σ = −κ. (a-d) Frequency-bin eigenket | − 1, 1⟩, which corresponds to a
single-peaked lineshape function. (e-h) Frequency-bin eigenket | − 8, 8⟩, which
corresponds to a doubled peaked lineshape function. . . . . . . . . . . . . . . . 94

x



5.3 The normalized second-order correlation functions g(2)−l,l(τ) for sidemodes l =
±9. The parameters are set to η = 0.9, ζ2 = κ/20 and g0|A0|2 = ηκ/50. Blue,
red, and green lines stand for the detuning σ = −κ, σ = 0, and σ = κ, respectively. 98

5.4 The spectral density Ssp(ω) (curves) and the coefficient |ck|2 (hollow circles) of
two types of combs. For the sake of visual clarity, we consider 10 pairs of the
side-modes and set Ω

R
= 20κ. The other parameters are set to g0|A0|2 = κ/10

and ζ2 = κ/20. (a) σ = κ. The output photon flux is double-peaked. (b) σ = −κ.
The output photon flux is single-peaked. . . . . . . . . . . . . . . . . . . . . . . 100

5.5 The comparison of the definition and approximate solutions. The parameters are
set to η = 0.9, σ = κ, and g0|A0|2 = ηκ/50. Red, green and blue dots stand for
Eq. (5.78), Eq. (5.79) and Eq. (5.80), respectively. (a) ζ2 = κ/100; (b) ζ2 = κ/20. 103

5.6 Color-coded representation of the density operator of Eq. (5.103) for (a) η = 0.1
and (b) η = 0.9. We consider K = 25 mode pairs and the parameters are set
to ζ2 = κ/20, σ = κ and g0|A0|2 = ηκ/50. For low η, the comb is mostly
incoherent with twinless photons, while for higher η, the comb is mostly coherent
with twin photons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

5.7 Schematic illustration of a time-bin entangled microcomb generator along with
the corresponding measurement system. WGMR: Whispering gallery mode res-
onator; BPF: Bandpass filter; PD: Single photon detector. . . . . . . . . . . . . . 113

xi



List of Abbreviations

BPF Bandpass filter

CME Coupled-mode equation
CW Continuous wave

EDFA Erbium-doped Fiber Amplifier

FFT Fast Fourier transform
FSR Free-spectral range
FWM Four-wave mixing

LLE Lugiato-Lefever equation

MgF2 Magnesium fluoride

PC Polarization controller
PD Photodetector

QED Quantum electrodynamics
QLE Quantum Langevin equation

RF Radio frequency

SFWM Spontaneous four-wave mixing

TIR Total internal reflection

VOA Variable optical attenuator

WDM Wavelength-division multiplexing
WGM Whispering-gallery mode
WGMR Whispering-gallery mode resonator

xii



Chapter 1: Introduction

1.1 Overview

In recent years, the nonlinear effects in optical resonators have become the focus of re-

searchers in photonics. Among these research, optical frequency combs based on high-quality

factor (high-Q) microresonators, which are also known as microcombs, are one of the most im-

portant topics as they provide a potential way to create efficient and compact multi-wavelength

sources for many applications. Microresonators can be fabricated by different materials and

to different shapes, including spheres, cylinders, and rings. Thanks to the new technologies,

researchers are able to fabricate high-Q optical resonators, such as whispering-gallery mode res-

onators (WGMR) and integrated ring resonators, with perimeters ranging from a few micrometers

to a few millimeters, corresponding to the free-spectral range (FSR) ranging from a few terahertz

to a few gigahertz. They are considered as the potential central elements of compact optical fil-

ters, optoelectronic oscillators, lasers, and modulators [1, 2], and have been found in numerous

applications such as coherent optical communications [3], spectroscopy [4], sensing [5, 6], and

ultra-pure microwave generation [7].

In the most basic configuration, when one eigenmode of the high-Q microresonator is

pumped by a continuous wave (CW) laser and the bulk medium of the cavity features a Kerr

nonlinearity, the pump photons can be frequency-converted via four-wave mixing (FWM) and

1



Figure 1.1: Schematic illustration of the experimental setup for Kerr microcombs generation
based on a 12 mm magnesium fluoride (MgF2) disk resonator at 1550 nm.
WGMR: Whispering-gallery mode resonator; PC: Fiber polarization controller; EDFA: Erbium-
doped fiber amplifier; VOA: Variable fiber optical attenuator; FC3dB: 1×2 50:50 fiber coupler;
PD: Photodetector. Source: Ref. [8].

populate adjacent cavity eigenmodes [9, 10]. This process can also be observed in nonlinear

optical fiber [11], while with a resonator of a small volume of confinement and a high-Q factor,

photons have a very long lifetime and the optical field stored in it is greatly enhanced, leading to

the same photon density at much lower pump power. Therefore, many nonlinear phenomena at

high-power optical fields can be observed in high-Q resonators with only a small power require-

ment. In amorphous media and centrosymmetric crystals, the Kerr effect is generally the leading

nonlinear effect compared the others, such as Raman and Brillouin effects, so Kerr microcombs

have received widespread interest. The pioneering breakthrough was made by two independent

groups at almost the same time in 2004, who demonstrated the Kerr-based degenerate parametric

oscillations in WGMRs [12, 13]. The broadband Kerr microcombs with non-degenerate four-

wave mixing were achieved in 2007 [14], and this idea was extended to integrated optical plat-

forms that are CMOS-compatible in 2009 by the other two groups [15, 16]. A typical schematic

of the experimental setup for Kerr microcombs generation is shown in Fig. 1.1.

The semiclassical dynamics of Kerr microcombs have already been well studied today by

2



the coupled-mode theory and the nonlinear Schrodinger equation, namely the Lugiato-Lefever

equation (LLE). In this prospect, four-wave mixing is considered to occur when the pump power

exceeds a certain threshold [17, 18, 19]. This four-wave mixing process is stimulated and the

comb lines exhibit strong phase correlation in the spectral domain, leading to the formation of

well-defined patterns in the spatio-temporal domain, which can be extended (such as roll pat-

terns), or localized (bright or dark solitons). These theoretical descriptions and the corresponding

modeling methods are explained in detail in Chapter 2.

However, there are still two main challenges for the semiclassical model under study – the

determination of the noise-driven phenomenon and the quantum description of the Kerr micro-

combs. Most applications of the comb above the threshold are based on the exceptional coher-

ence of the comb lines, which is ultimately limited by noise-driven fluctuations. Compared with

other frequency comb generation systems, Kerr comb generators are simpler and more compact

because of the small mode volume, high photon density, and long photon lifetime, but these ad-

vantages may introduce more noise. Researchers have studied various noises in the resonator

theoretically and experimentally [20, 21, 22], while the mechanism by which these noises affect

the fluctuations of the microcombs, and further the phase noise, remains unclear.

Chapter 3 introduces a stochastic approach to characterize the coherence of the combs in

the optical domain below and above the threshold. We focus on the case of additive Gaussian

white noise because, despite its conceptual simplicity, it already unveils a complex interaction

between nonlinearity and noise that leads to profound consequences in terms of metrological

performance. This analysis allows us to understand how the fluctuations of the optical comb

lines are converted to microwave phase noise after the photodetection of the comb.

We anticipate that our results in Chapter 3 can be extended to more complex resonators

3



with high-Q factors and more complex noise spectra. Thus, a more detailed description of the

noises in the system is given in Chapter 4, where the multiplicative noises induced by thermal

fluctuations and photothermal fluctuations are shown as examples. We also study a method to

reduce the phase noise generated by detuning fluctuations, which is called the quiet point effect.

On the other hand, the quantum description of microcombs requires the study of quan-

tum electrodynamics (QED). From the quantum viewpoint, entangled twin photons are gener-

ated by spontaneous four-wave mixing (SFWM) below the threshold and the quantum combs

may achieve entanglement over tens of cavity eigenmodes. Such entanglement properties make

quantum combs particularly suitable for wavelength-division multiplexing (WDM) schemes in

quantum communications, and the high dimensionality also helps to achieve robust distribution

of entanglement over multi-node networks [23, 24, 25, 26, 27, 28, 29, 30]. More recently, people

have started to focus on the quantum properties of the entangled light, such as the squeezing,

when the power of the pump is around or above the threshold [31, 32]. Therefore, deriving a

complete and comprehensive quantum model is necessary and important.

In order to solve the quantum properties of quantum microcombs, Chapter 5 demonstrates

a model to describe the biphoton quantum superposition in quantum Kerr microcombs below the

threshold. We propose a method to find the quantum expansion of the output state, as well as

the lineshape function, amplitude, and build-in phase. This frequency-bin approach also allows

us to determine an explicit steady-state density operator for quantum microcombs and a com-

plete characterization of these quantum systems in terms of photodetection statistics, which are

experimentally accessible measurements.

4



1.2 Kerr Effect

The Kerr effect, also known as the quadratic electro-optic effect, is a change in the refractive

index of a material due to an applied electric field. It is related to the third-order susceptibility

χ(3) and can be caused by an external electromagnetic field or the field of the incident light itself.

In our research, we mainly consider the second type, which is also known as the optical Kerr

effect.

In the optical Kerr effect, the change in refractive index is proportional to the local intensity

(irradiance) of light, which can be written as

n = n0 + n2I , (1.1)

with the nonlinear Kerr coefficient

n2 =
3

4

χ(3)

ε0cn2
0

. (1.2)

In this equation, the parameter n0 is the refractive index of the material, I is the optical intensity,

ε0 is the vacuum permittivity, and c is the velocity of light in vacuum. We can see the optical Kerr

effect will become significant in strong light (high photon density), such as in optical resonators

with small modal volumes and pumped by external laser beams.

5



Figure 1.2: (a) Ray optics viewpoint and (b) Wave optics viewpoint of the whispering-gallery
mode in a WGMR.

1.3 Whispering-Gallery Mode Resonators

The term ”whispering-gallery mode” was introduced and studied by Lord Raleigh in the

late 19th century to describe an acoustic mode observed at St. Paul’s Cathedral in London, which

allowed whispers to be heard distinctly from one side of the cathedral’s dome to the other. The

phenomena are firstly described as the total internal reflection (TIR) of acoustic waves along

the inner circumference of the dome and then analyzed using the wave’s point of view [33, 34].

For this reason, the optical modes along the inner circumferences of optical resonators are also

referred to by this name.

We focus on whispering-gallery mode resonators here because of their conceptual sim-

plicity and practical importance. As shown in Fig. 1.2, in the ray optics viewpoint of WGM

resonators, light only progrates along and is reflected by the edge. As a result, in wave optics, the

electrical field only focuses near the edge.

Various shapes, materials, and technologies have been used to fabricate WGM resonators,

mainly aiming at achieving high quality factors (or equivalently, low loss) in different platforms.

Nowadays, crystalline WGMRs made by fluoride crystals, such as calcium fluoride and magne-

6



Figure 1.3: Schematic diagrams of microresonators. From left to right: Whispering-gallery
mode-based bulk toroids, spheres, monolithic toroids, and integrated ring-resonators. Source:
Ref. [19].

sium fluoride, have been found to feature very low intra-cavity losses and surface scattering loss.

They are mostly millimeter-level dimensions and can achieve Q-factors on the scale of tens of

billions, or even higher [35, 36, 37].

On the other hand, integrated microresonators can be fabricated based on fused silica, sil-

icon nitride (Si3N4), silicon oxynitrides (SiOxNy), and Hydex glass, which achieve Q factors on

the scale of several millions or higher [14, 38, 39]. Although featuring more losses, the chip-

scale size of the integrated microresonators makes them compatible with large-scale integration.

Strictly speaking, waveguide-based integrated microresonators are not whispering-gallery mode

resonators because the working modes in them are waveguide modes rather than whispering-

gallery modes. However, there is no large difference between their theories, so we don’t make a

distinction here.

The optical fibers can be coupled to the WGMRs in different ways, as displayed in Fig. 1.4.

The prism coupling is generally very robust and can be used for WGM modes of the resonators

with very high refractive index, with a coupling efficiency of up to 70%, while waveguide cou-

7



Figure 1.4: The coupling methods for the WGM resonators with optical waveguides: (a) Prism
coupling, (b) Waveguide side coupling, (c) Waveguide tip coupling. Source: Ref. [40].

pling is very efficient, resulting in a coupling efficiency of more than 90%. In this dissertation,

we consider the waveguide side-coupling method as our default experimental method.

8



Chapter 2: Semiclassical Dynamics of Kerr Microcombs

2.1 Introduction

The theoretical study of Kerr microcombs starts from the research on their semiclassical

dynamics. In the previous chapter, we briefly introduced the basic principles of Kerr comb gen-

eration, including the background, the definition of the Kerr effect, and whispering-gallery mode

resonators. In this chapter, we aim to propose a time-domain mathematical model to describe the

formation process of Kerr microcombs.

2.2 System and Parameters

The main element of a Kerr microcomb generator is a high-Q resonator. In this chapter,

we consider a whispering-gallery mode resonator or integrated ring resonator with a closed-path

perimeter L = 2πa, pumped with a resonant and continuous wave laser of power P
L

and angular

frequency ω
L
. The size of the resonator defines its FSR as

F
R
=

1

T
R

=
Ω

R

2π
=
vg
L

=
c

ngL
, (2.1)

where Ω
R

is the angular FSR, c is the velocity of light in vacuum, ng is the group-velocity refrac-

tion index at the pump wavelength ω
L
, and vg = c/ng is the corresponding group velocity. The

9



parameter T
R

is the round-trip period of photons traveling along the closed-path perimeter of the

resonator.

The losses in the resonator are characterized by the half-linewidth

κ =
ω

L

2Q
, (2.2)

where Q is the total quality factor. In fact, the total quality factor can be calculated as

Q−1 = Q−1
i +Q−1

e , (2.3)

where Qi and Qe are the intrinsic and extrinsic quality factors, respectively. It leads to

κ = κi + κe , (2.4)

with κi and κe being the half-linewidth of the intrinsic and extrinsic coupling. Although the Q

factors and linewidths depend on their eigenmode order l, we consider they are degenerated and

identical across different modes for the sake of simplicity.

The extrinsic coupling comes from the input/output coupling system, which has two main

configurations, as shown in Fig. 2.1, depending on whether input and output are performed by

the same coupler. In the add-through coupling configuration, as there is only one coupler with a

loss κt, the extrinsic half-linewidth is equal to the loss of that coupler:

κe = κt . (2.5)

10



Figure 2.1: Two main coupling configurations for Kerr comb generation. (a) Add-through cou-
pling. (b) Add-drop coupling. Source: Ref. [41].

In the add-drop coupling configuration, there are two couplers with loss κt and κd, so the extrinsic

half-linewidth is the sum of them

κe = κt + κd . (2.6)

Generally speaking, since the coupling loss of the add-through configuration is smaller, the sys-

tem may have a smaller threshold power and therefore can run at a lower power. However, as

the second configuration can arbitrarily control the direction of the outgoing field, which makes

people potentially able to save space, and the corresponding output field is strictly proportional

to the intracavity field, it is still very useful for many applications.

11



2.3 Eigenfrequences and Eigenmodes

To model Kerr microcombs, the first step is solving the spatial mode distribution for the

resonator operating in the linear and lossless regime, that is, solving the Helmholtz equation

[
∇2 +

ω2
l ϵ(r, ωl)

c2

]
Υl(r) = 0 , (2.7)

where the relative permittivity ϵ(r, ωl) = n2(r, ωl) depends on the structure of the resonator. The

torus-like spatial mode profile Υl(r) is in units of m−3/2 and normalized to 1. The photons in the

resonator are trapped in the eigenmodes and the effective mode volume can be calculated as

Veff =

[∫
V

|Υl(r)|4dV
]−1

, (2.8)

which is generally much smaller than the bulk volume of the resonator. Sometimes, instead of the

effective volume, a parameter Aeff = Veff/L is used to define the effective region of the resonator

trapping the photons, which is called the effective mode area. Again, although this effective

region depends on the eigenmode order l, we assume they are degenerated and identical here to

simplify the problem.

We will see the solutions of Eq. (2.7) are the longitudinal eigenmodes in many transverse

mode families, as displayed in Fig. 2.2. In the same transverse mode family, the longitudinal

eigenmodes can be labeled with their azimuthal eigenumber ℓ. Assuming the laser is pumping

the eigenmode with the azimuthal eigenumber ℓ0, we can define a reduced azimuthal eigenumber

l = ℓ − ℓ0 such that l = 0 is the laser-pumped mode and the sidemodes are expanded as l =

12



Figure 2.2: Distributions of electric fields of WGMs in a spherical WGM resonator. (a) Side view
(left) and top view (right) of the mode with n = 1, ℓ = 30, and ℓ − |m| = 0. (b) Examples of
transverse field distributions for different radial and polar eigennumbers n andm. Each transverse
field distribution corresponds to a family of modes with different ℓ. Source: Ref. [18].

±1,±2,±3, · · · . With this reduced azimuthal eigenumber, the corresponding eigenfrequencies

can be Taylor-expanded as

ωl = ω0 +
+∞∑
n=1

ζnl
n

n!
, (2.9)

where ω0 is the eigenfrequency of the laser-pumped mode, ζ1 = Ω
R

, while ζn is the nth-order dis-

persion coefficient for n ≥ 2. It should be noted that these dispersion coefficients are determined

by both geometric dispersion and material dispersion.

The Kerr effect of the bulk material of this optical cavity is characterized by a Kerr coeffi-

cient n2, which can be rescaled to the nonlinear coefficient in units of W−1m−1

γ =
ω

L
n2

cAeff

, (2.10)

13



or to the FWM gain in units of rad/s

g0 =
n2cℏω2

L

n2
gVeff

, (2.11)

with ℏ being the reduced Planck constant.

With the above parameters, the semiclassical electric field (in units of V/m) inside the

resonator can be written as two equivalent expansions with different normalizations. The first

one follows

E(r, t) =

√
2
ℏω

L

ε0n2
g

∑
l

1

2
Al(t)e

iωltΥl(r) + c.c. , (2.12)

where Al(t) is the complex slowly-varying amplitude of the optical field for the mode l, nor-

malized such that |Al|2 is dimensionless and represents the number of intracavity photons in that

mode. The variable t stands for time, ε0 is the vacuum permittivity, and c.c. stands for the com-

plex conjugate of all the preceding terms. The other option is to write the total electrical field in

the resonator as

E(r, t) =

√
2
T

R

ε0n2
g

∑
l

1

2
El(t)e

iωltΥl(r) + c.c. , (2.13)

where El(t) is a complex slowly-varying amplitude for each mode l with |El|2 in units of watts.

In the quantum analysis, we are more concerned with the photon number and the output photon

flux, so we start from the first expansion, while in the semiclassical and stochastic analysis, the

power of the intracavity and outgoing fields is more important, so we prefer the second expansion

instead.

14



Although expansions based on eigenfrequencies are the most common and well-known, it

should be noted that for the process of four-wave mixing, this is not the optimal way. As we know,

the four-wave mixing process of photons needs to obey energy conservation law but because

of the dispersion, the eigenfrequencies are not equally spaced, which prevents photons from

jumping between two eigenfrequencies. Obviously, the center of each resonance in a stationary

Kerr microcomb is not located at the eigenfrequency, and instead, they are equally spaced by a

repetition rate Ωrep. Therefore, in the modeling, the best way is to expand the field based on this

repetition rate.

When we consider up to the second-order dispersion and ignore the other nonlinear ef-

fects, the repetition rate of a Kerr microcomb is exactly the angular FSR, so we can rewrite the

expansions in Eqs. (2.12) and (2.13) as

E(r, t) =

√
2
ℏω

L

ε0n2
L

∑
l

1

2
Al(t)e

iΩltΥl(r) + c.c. , (2.14)

and

E(r, t) =

√
2
T

R

ε0n2
L

∑
l

1

2
El(t)eiΩltΥl(r) + c.c. , (2.15)

with the complex slowly varying amplitudes defined as

Al(t) = Al(t) exp [(ωl − Ωl)t] , (2.16)

El(t) = El(t) exp [(ωl − Ωl)t] . (2.17)

15



The angular frequencies of the expansions become

Ωl = ω
L
+ Ω

R
l , (2.18)

which are symmetrical to the laser angular frequency and equally spaced.

2.4 Semiclassical Dynamics – Coupled-Mode Equations

The semiclassical dynamics of Kerr microcombs can be investigated using two equivalent

models, namely the coupled-mode model [42, 43] and the Lugiato-Lefever model [44, 45, 46, 47].

The coupled-mode model is the spectro-temporal model describing the interaction of the fields

of the sidemodes, while the Lugiato-Lefever model shows the spatio-temporal behavior of the

optical fields. We focus on the first one in this section.

2.4.1 Coupled-Mode Equations with Modal Fields in the Normalization of Pho-

ton Number

In the coupled mode equations, the fields of each mode are described by the envelope

functions Al from the expansion in Eq. (2.14). Therefore, the dynamics of these fields should be

a set of coupled equations of the loss, dispersion, nonlinearity, laser pump power, etc. Assuming

the eigenmodes under consideration are quasi-degenerate, the coupled-mode equations can be

16



written as

Ȧl = −κAl + i

[
σ −

+∞∑
k=2

ζk
k!
lk

]
Al + δ(l)

√
2κeAin

+ig0
∑
m,n,p

δ(m− n+ p− l)AmA∗
nAp , (2.19)

where Ȧl = dAl/dt is the time derivative and the parameter σ = ω
L
− ω0 is the detuning

between the laser angular frequency and the cold-cavity resonance frequency of the pumped

mode. The function δ(l) is the Kronecker delta function, indicating that the central mode l = 0

is the only mode being pumped by the laser. We consider the phase of the external pump field as

the reference and set it to zero, so the normalized pump field is a real number, following

Ain =
√
P

L
/ℏω

L
, (2.20)

with |Ain|2 being the pump photon flux (in units of s−1). The modal output fields can be calculated

as

Aout,l =
√
2κeAl −Ainδ(l) (2.21)

and

Aout,l =
√
2κdAl (2.22)

for the add-through and add-drop coupling configurations, respectively.

17



2.4.2 Coupled-Mode Equations with Modal Fields in the Normalization of In-

tracavity Power

The coupled mode equations with the field normalized to the intracavity field power in

Eq. (2.15) can be obtained from the previous result, following

Ėl = −κ El + i

[
σ −

+∞∑
k=2

ζk
k!
lk

]
El + δ(l)

√
2κe/TR

√
P

L

+ivgγ
∑
m,n,p

δ(m− n+ p− l) EmE∗
nEp , (2.23)

referring to the phase and frequency of the pump field as zero. The modal output fields for

add-through coupling and add-drop coupling configurations can be written respectively as

Eout,l =
√

2κeTR
El −

√
P

L
δ(l) (2.24)

and

Eout,l =
√
2κdTR

El , (2.25)

which are also in units of watts.

2.5 Semiclassical Dynamics – Lugiato-Lefever Equation

The LLE was initially introduced by Luigi Lugiato and René Lefever in 1987 to describe

spontaneous pattern formation in nonlinear optics. It has been shown that the spatiotemporal

18



dynamics of Kerr microcombs can be modeled by the LLE, and it’s exactly equivalent to the

CME approach we introduced above.

2.5.1 Lugiato-Lefever Equation with Intracavity Field in the Normalization of

Photon Number

With the mode expansion fields we introduced before, the intracavity field amplitude of the

resonator can be calculated as

A(θ, t) =
∑
l

Al(t)e
ilθ , (2.26)

where θ is the azimuthal angle of the resonator ranging from −π to π. Then, the LLE can be

written as

∂A
∂t

= −κA+ iσA+ i
+∞∑
k=2

ζk
k!

∂kA
∂θk

+ ig0 |A|2A

+
√
2κt Ain , (2.27)

with the output fields

Aout =
√
2κtA−Ain (2.28)

for the add-through coupling configuration, and

Aout =
√
2κdA . (2.29)

19



for the add-drop coupling configuration. We can see the output field predicted by the LLE relies

on two timescales θ and t, where θ serves as the fast timescale showing the periodic shape of the

mode in a round-trip period of photons, while t is the slow timescale describing the slow-varying

envelope of the output field.

2.5.2 Lugiato-Lefever Equation with Intracavity Field in the Normalization of

Intracavity Power

Equation (2.27) can also be rescaled to

∂E
∂t

= −κ E + iσE + i
+∞∑
k=2

ζk
k!

∂kE
∂θk

+ ivgγ |E|2E

+
√

2κt/TR

√
P

L
, (2.30)

where E(θ, t) = ∑
l El(t)eilθ is the total intracavity field, with the output fields (in units of Watts)

being written as

Eout =
√
2κtTR

E −
√
P

L
(2.31)

and

Eout =
√

2κtTR
E (2.32)

for the two configurations.

20



2.5.3 Lugiato-Lefever Equation with Intracavity Field in the Dimensionless

Time Normalization

The original LLE is defined with a dimensionless time by L. Lugiato and R. Lefever [48],

making it more convenient for mathematical analysis. In this way, the normalized intracavity

field is rescaled as

ψ(θ, t) = (g0/κ)
1/2A(θ, t) , (2.33)

while the LLE can also be rewritten as

∂ψ

∂τ
= −(1 + α)ψ − i

β

2
ψ + i|ψ|2ψ + F , (2.34)

where τ = κt is the dimensionless time, α = −σ/κ is the normalized detuning, and β = −ζ2/κ

is the second order dispersion. The dispersion is called normal dispersion when β is positive and

is called anomalous when it’s negative. We restrict the dispersion up to the second order and

ignore the higher-order dispersion here. The external pump power is therefore normalized to

F =

(
2g0κt
κ3

) 1
2

√
P

L

ℏω
L

, (2.35)

which is also real-valued.

21



2.6 Numerical Simulation Methods

LLEs can be simulated efficiently by the split-step Fourier algorithm, which is based on

the fast Fourier transform (FFT) algorithm and assumes a periodic boundary condition, while

the simulations of CMEs are used to be slower because of calculations of the nonlinear terms.

However, it has been shown later that FFT can also be applied to CMEs to optimize the simulation

of the nonlinear terms [49], making these two approaches equivalent in most respects of the

simulation. It’s interesting to note that the intracavity field is actually related to the modal fields

with a Fourier transform, so the optimization method to simulate CMEs is in theory the same as

the split-step Fourier algorithm we used before, only with a different starting point. More accurate

simulations can be achieved using Runge-Kutta methods of fourth order or higher, which is very

useful for the simulation of fast-changing fields, e.g. thermo-optical relaxation oscillations [50,

51].

2.7 Spatio-Temporal and Spectro-Temporal Solutions

Kerr microcomb solutions can be obtained from the simulations of CMEs or LLE. Below

the threshold, the central mode is solely excited by the pump laser, which is called a flat state (No

Kerr comb is generated). The Kerr combs are mainly of three types, namely rolls, solitons, and

chaos. Among them, two types of stationary solutions – rolls and solitons, are our focus in this

dissertation. Some examples of the stationary solutions are displayed in Fig. 2.3 considering only

the second-order dispersion. It should be noted that bright solitons are generated in resonators

with anomalous dispersions, while dark solitons are generated in them with normal dispersions.

22



(a)

a
0

−a

−a

0

a

0

20

40
|E
(θ
)|2

[a
.u
.]

(b)

a
0

−a

−a

0

a

0

20

40

|E
(θ
)|2

[a
.u
.]

(c)

a
0

−a

−a

0

a

0

20

40

|E
(θ
)|2

[a
.u
.]

(d)

a
0

−a

−a

0

a

0

20

40

|E
(θ
)|2

[a
.u
.]

(e)

a
0

−a

−a

0

a

0

20

40

|E
(θ
)|2

[a
.u
.]

(f)

−30 −15 0 15 30

0

−60

−120

l

|E
o
u
t,
l|2

[d
B
m
]

(g) Comb 1

−200 −100 0 100 200

0

−60

−120

l
|E

o
u
t,
l|2

[d
B
m
]

(h) Comb 2

−200 −100 0 100 200

0

−60

−120

l

|E
o
u
t,
l|2

[d
B
m
]

(i) Comb 3

−30 −15 0 15 30

0

−60

−120

l

|E
o
u
t,
l|2

[d
B
m
]

(j) Comb 4

−30 −15 0 15 30

0

−60

−120

l

|E
o
u
t,
l|2

[d
B
m
]

Figure 2.3: Schematic representation of stationary Kerr microcombs. The microcombs are in
space [(a)–(e)] and frequency [(f)–(j)], after simulation of Eq. (2.23). The upper row displays the
total intracavity field power |E(θ)|2, while the lower row displays the corresponding stem plot for
the modal intensities |Eout,l|2 in a logarithmic scale (note that they are spectrotemporal snapshots
and not Fourier spectra). (a) and (f) Flat state (no comb). (b) and (g) Bright soliton (P

L
= 4.00

mW, σ = −2κ, and ζ2 = 2π × 2.92 kHz). (c) and (h) Dark soliton (P
L
= 5.30 mW, σ = −2.5κ,

and ζ2 = −2π × 2.92 kHz). (d) and (i) Roll pattern of order L = 1 (P
L
= 3.09 mW, σ = −κ,

and ζ2 = 2π × 874.7 kHz). (e) and (j) Roll pattern of order L = 4 (P
L
= 2.27 mW, σ = −κ, and

ζ2 = 2π × 366.4 kHz).

2.8 Eigenvalue Analysis

The eigenvalue analysis can be applied to both CMEs and LLE to analyze the stability of

stationary solutions [52, 53, 54]. We show an example of applying the eigenvalue analysis to

CMEs with only the second-order dispersion in this section.

With only second-order dispersion, the coupled mode equations Eq. (2.23) can be simpli-

fied to

Ėl = −κ El + i

[
σ − 1

2
ζ2l

2

]
El + δ(l)

√
2κt/TR

√
PL (2.36)

+ivgγ
∑
m,n,p

δ(m− n+ p− l) EmE∗
nEp .

Considering 2N + 1 modes centered at l = 0 where N is a very large number so modes other

23



than these are small enough to be ignored, any stationary solution

E = (E−N , E−N+1, ..., EN−1, EN) (2.37)

can be found by setting Ėl = 0 and solving the equations

0 = −κEl + i

[
σ − 1

2
ζ2l

2

]
El + δ(l)

√
2κt/TR

√
PL (2.38)

+ivgγ
∑
m,n,p

δ(m− n+ p− l)EmE
∗
nEp .

Eigenvalue analysis is based on the Jacobian matrix of the coupled equations at this solution.

2.8.1 Eigenvalue Analysis of Solitons

For bright and dark solitons, all of the modes adjacent to the pump frequency are excited.

When we perturbed the stationary soliton solution, the modal field of mode l can be written as

El = El + δEl , (2.39)

where the perturbations form

δE = (δE−N , δE−N+1, ..., δEN−1, δEN) . (2.40)

24



Substrating Eq. (2.39) into the CMEs in Eq. (2.36) and ignoring the product of perturbations, we

can get the dynamics of the perturbations

δĖl = −κ δEl + i

[
σ −

+∞∑
k=2

ζk
k!
lk

]
δEl

+ivgγ
∑
m,n,p

δ(m− n+ p− l) [EmE∗
nδEp + E∗

nEpδEm + EmEpδE∗
n] , (2.41)

which can be rewritten in the form of a matrix, following

 δĖ

δĖ∗

 = J

 δE

δE∗

 , (2.42)

where the complex-valued Jacobian matrix

J =

 R S

S∗ R∗

 (2.43)

is explicitly defined as

Rlp =

[
−κ+ i(σ − ζ2

2
l2)

]
δ(p− l) (2.44)

+2ivgγ
∑
m,n

δ(m− n+ p− l) EmE∗
n ,

Slp = ivgγ
∑
m,n

δ(m+ n− p− l) EmEn . (2.45)

The stability of the stationary solution can be determined by calculating the real parts of the

eigenvalues. If any of them is positive, the soliton solution is unstable.

It is interesting to note that in fact, this Jacobian of stationary Kerr microcombs always has

25



an eigenvalue of 0, which means that the so-called “stable” Kerr comb solution is actually only

neutrally stable. The mode corresponding to this eigenvalue 0 is called the Goldstone mode and

we will discuss it in detail in Chapter 3.

2.8.2 Eigenvalue Analysis of Rolls

For roll patterns of orderL, only part of the sidemodes are excited and the spacings between

them are LΩ
R

. Assuming there are M excited modes in the sidemodes under consideration, the

Jacobian matrix will be block diagonal with L/2 + 1 boxes if L is even, and (L + 1)/2 boxes if

L is odd. For convenience, we rewrite the Jacobian matrix to make the q-th block Jq correspond

to the perturbations:

δE (q) = (δE−ML−q, ..., δEML−q, δE−ML+q, ..., δEML+q) , (2.46)

which is sometimes referred to as a “Bloch mode”. The block Jacobian Jq can be explicitly

written as

Jq =

 Rq Sq

S∗
q R∗

q

 . (2.47)

For the matrices with q ̸= 0 and q ̸= L/2, Rq and Sq can be calculated as

Rq =

 Uq 0

0 Vq

 and Sq =

 0 Wq

W∗
q 0

 , (2.48)

26



where Rq and Sq are block matrices of order (4M + 2) with complex-valued elements

U (q)
lp = RlL−q,pL−q (2.49)

V(q)
lp = RlL+q,pL+q (2.50)

W(q)
lp = SlL+q,pL−q = SlL−q,pL+q . (2.51)

For the cases q = 0 and q = L/2, we find that Uq and Vq are identical, so that Rq and Sq are

matrices of order (2M + 1) or (2M + 2), with

Rq = Uq and Sq = Wq . (2.52)

The stability of the stationary roll solutions can also be determined by calculating the real parts of

the eigenvalues. Any eigenvalues with positive real parts indicate the solution is unstable, while

the corresponding Bloch mode shows unstable sidemodes of the microcomb.

2.8.3 Subharmonic Instabilities of Rolls

The theory developed in the previous section can be used in explaining the phenomenon

of subharmonic excitations, where primary combs with subcomb lines have been observed in the

experiment.

The left column of Fig. 2.4 shows the simulation results of the subcombs generated by a

MgF2 disk resonator with a Q-factor of one billion at 1550 nm and a diameter of 12 mm. In the

experiment, the single FSR value is measured as 5.9 GHz, and the primary comb spacings for

these five combs are 50 FSR, 38 FSR, 42 FSR, 84 FSR, and 45 FSR, respectively, by exciting the

27



modes from different families of WGMs with different dispersions [55].

The results of this eigenvalue analysis are presented in the left column of Fig. 2.4, and they

confirm that the subcombs originate from the most unstable Bloch mode, i.e., the one with the

largest positive real part. We can see there is always an eigenvalue with real value 0 for the block

q = 0 (Data points in the upper left of the figures in the right column), corresponding to the

existence of the Goldstone mode.

In fact, the modulated patterns arising from this instability are typically unstable and can

be observed as transients. When a pattern becomes unstable, it gets modulated by the growth of

small-wavenumber perturbations. If the solution is decomposed as the modulus and phase, the

small-wavenumber modulations correspond to variations of the phase with a long spatial scale.

The dynamics of the modulus can be relaxational, i.e., it has negative eigenvalues and can be

adiabatically eliminated, while the dynamics of the phase are diffusive with significantly slower

dynamics. As a consequence, the transients persist for a long time before the small-wavenumber

modulations are smoothed out. However, the modulation can also be eventually stabilized by

other intracavity nonlinear effects.

2.9 Conclusion

In this chapter, we have introduced the fundamental parameters needed to understand our

Kerr microcomb model. We have presented two deterministic models that allow us to under-

stand the dynamics of Kerr microcombs generated in the whispering-gallery mode resonator in

both spectro-temporal and spatio-temporal regimes. In particular, we show an example using the

eigenfrequency analysis method of the deterministic model, which enables us to unveil that the

28



Figure 2.4: [Left column] Five types of Kerr comb spectra snapshots of the output fields. The
quality factor of the through port is set to Qt = Q/4, and the nonlinear parameter is set to
γ = 1.0 W−1km−1. From top to bottom, the detuning σ of each case are 2κ, 2.3κ, 1.8κ,
1.83κ, and 0.55κ, respectively, and the 2nd order dispersion ζ2 are set to −0.0037κ/vgω

2
FSR,

−0.0073κ/vgω
2
FSR, −0.0060κ/vgω

2
FSR, −0.0015κ/vgω

2
FSR, and −0.0055κ/vgω

2
FSR. Pump pow-

ers PL are 75 mW, 90 mW, 100 mW, 100 mW, and 75 mW, corresponding to the thresholds: 48.4
mW, 57.5 mW, 42.8 mW, 43.6 mW, and 16.5 mW. [Center column] Eigenvalues of the solution
in a complex plane. [Right column] The real part of the eigenvalues is plotted as a function of the
wavenumber q of the Bloch modes. The solution is unstable when the real part of one of them
becomes positive, and the eigenvalue with the largest positive real part corresponds to the order
of the subcomb.

29



origin of the additional subcomb lines is a subharmonic bifurcation leading to the excitation of a

Bloch mode. This method can also be used to explain the interaction of this instability with other

nonlinearities in the optical resonator [56].

30



Chapter 3: Stochastic Dynamics and Phase Noise of Kerr Microcombs

– Additive Gaussian White Noise

3.1 Introduction

Kerr microcombs are known for their stability and precision, so they are often referred to

as optical clocks and optical rulers. Most of these applications are related to their metrological

performance, which is ultimately limited by their noise-driven fluctuations. For this reason, it

is of high importance to understand the influence of random noise on the microcomb dynamics.

In this chapter, we theoretically investigate a model where Gaussian white noise is added to the

coupled-mode equations governing the microcomb dynamics. This stochastic model allows us to

characterize the noise-induced broadening of the spectral lines, and permits us to determine the

Laser

Kerr comb generator

CombPump

WGMR
PD

Phase noise 

analyzer

Microwave

Figure 3.1: Schematic illustration of microwave generation using a Kerr optical frequency comb.
WGMR: Whispering-gallery mode resonator; PD: Photodiode.

31



phase noise spectra of the microwaves generated via comb photodetection. In this latter case, our

analysis indicates that the low-frequency part of the spectra is dominated by pattern drift while

the high-frequency part is dominated by pattern deformation.

3.2 Stochastic Model for Additive Gaussian White Noise

In order to study the stochastic dynamics of Kerr microcombs, we rewrite Eq. (2.23) with

Gaussian white noise terms, following

Ėl = −κ El + i

[
σ −

+∞∑
k=2

ζk
k!
lk

]
El + δ(l)

√
2κe/TR

√
P

L

+ivgγ
∑
m,n,p

δ(m− n+ p− l) EmE∗
nEp

+
√
2κΛ vl(t) , (3.1)

where the complex-valued additive Gaussian white noise terms vl(t) have the correlation proper-

ties:

⟨vl(t)⟩ = 0 , (3.2)

⟨vl(t)v∗l′(t′)⟩ = δl,l′ δ(t− t′) . (3.3)

The delta symbols with subscripts represent Kronecker functions while the delta symbols fol-

lowed by parentheses represent Dirac delta functions. These noise terms are weighted in each

mode by the real-valued amplitude Λ which is such that Λ2 is in unit of watts. This approach

of coupled-mode equations with additive Gaussian noise was introduced in Ref. [57], and some

statistical properties of the comb had been analytically derived.

32



The output optical field Eout,l is related to the intra-cavity field by the same relation we

derive before

Eout,l =
√

2κeTR
El −

√
P

L
δ(l) , (3.4)

still with |Eout,l|2 in units of watts. A schematic representation of the experimental system is

presented in Fig. 3.1.

We use the following parameters for our numerical simulations: The single-frequency laser

is assumed to be ideal with an infinitely narrow linewidth at wavelength λ0 = 1550 nm. The

refraction index of the resonator at that frequency is ng = 1.43, corresponding to Ω
R
= 2π ×

13.4 GHz. The nonlinear parameter is set to γ = 1.0 W−1km−1. The intrinsic quality factor

is Qi = 109, and the coupling quality factor is Qe = 0.25 × 109. The additive Gaussian noise

amplitude is Λ = 5× 10−6 W1/2. We neglect dispersion beyond the second order.

Simulation of Eq. (3.1) yields a variety of solutions depending on the system parameters

and initial conditions. These solutions for El lead to distinctive intracavity spatiotemporal pat-

terns, as shown in Fig. 3.2. In this chapter, we will focus on four different patterns, namely bright

solitons, dark solitons, roll patterns with single-FSR spacing, and roll patterns with multiple-FSR

spacing. We also assume that the system is in a stationary and noise-free state for t < 0, while

the noise is switched on at t = 0 and the system thereby becomes stochastic.

Numerical simulations are performed by integrating the stochastic coupled-mode equations

in Eq. (3.1) using the fourth-order Runge–Kutta method in MATLAB. We first compute the sta-

tionary solution without noise, and then add Gaussian white noise terms to the simulations of the

microcomb dynamics.

33



(a)

a
0

−a

−a

0

a

0

20

40
|E
(θ
)|2

[a
.u
.]

(b)

a
0

−a

−a

0

a

0

20

40

|E
(θ
)|2

[a
.u
.]

(c)

a
0

−a

−a

0

a

0

20

40

|E
(θ
)|2

[a
.u
.]

(d)

a
0

−a

−a

0

a

0

20

40

|E
(θ
)|2

[a
.u
.]

(e)

a
0

−a

−a

0

a

0

20

40

|E
(θ
)|2

[a
.u
.]

(f)

−30 −15 0 15 30

0

−60

−120

l

|E
o
u
t,
l|2

[d
B
m
]

(g) Comb 1

−200 −100 0 100 200

0

−60

−120

l
|E

o
u
t,
l|2

[d
B
m
]

(h) Comb 2

−200 −100 0 100 200

0

−60

−120

l

|E
o
u
t,
l|2

[d
B
m
]

(i) Comb 3

−30 −15 0 15 30

0

−60

−120

l

|E
o
u
t,
l|2

[d
B
m
]

(j) Comb 4

−30 −15 0 15 30

0

−60

−120

l

|E
o
u
t,
l|2

[d
B
m
]

Figure 3.2: Schematic representation of stochastic Kerr microcombs. The combs are in space
[(a)–(e)] and frequency [(f)–(j)], and are at a given time t (snapshot), after simulation of Eq. (3.1)
with the same parameters as Fig. 2.3. The upper row displays the total intracavity field power
|E(θ)|2, while the lower row displays the corresponding stem plot for the modal intensities |Eout,l|2
in a logarithmic scale (note that they are spectrotemporal snapshots and not Fourier spectra). The
noise in the spatiotemporal patterns has been accentuated for the sake of visual clarity. (a) and (f)
Flat state (no comb). (b) and (g) Bright soliton. (c) and (h) Dark soliton. (d) and (i) Roll pattern
of order L = 1. (e) and (j) Roll pattern of order L = 4. The numerical simulations will involve
the four cases from Comb 1 to Comb 4.

3.3 Optical Fluctuations Under Threshold

In the sub-threshold regime, all the sidemodes l ̸= 0 are fluctuating around zero, while

the amplitude of the central mode l = 0 fluctuates around a non-zero amplitude. The fields can

therefore be rewritten as

El(t) =


E0(0) + δE0(t) for l = 0

δEl(t) for l ̸= 0

, (3.5)

where δEl denotes the fluctuation fields, while the amplitude E0(0) is the solution of the nonlinear

algebraic equation

0 = −κ E0(0) + i

[
σ − ζ2

2
l2
]
E0(0) (3.6)

+ivgγ|E0(0)|2E0(0) +
√
2κe/TR

√
P

L
.

34



If we considerN sidemode pairs around the central mode, the stochastic dynamics Eq. (3.1)

can be rewritten as a set of 2N + 1 differential equations after ignoring the multiplicative fluctu-

ation terms, in the form

δĖl = Rl δEl + Sl δE∗
−l +

√
2κΛ vl(t) (3.7)

for the fluctuations in each mode, with

Rl = −κ+ i

[
σ − ζ2

2
l2
]
+ 2ivgγ |E0(0)|2 , (3.8)

Sl = ivgγ E2
0 (0) , (3.9)

being complex-valued parameters. Here we assume the number N is very large so all of the other

sidemodes can be ignored. Since the fluctuations are pairwise coupled, Eq. (3.7) can be further

rewritten as N + 1 independent sets of 2× 2 noise-driven linear flows, in the form

 δĖl

δĖ∗
−l

 = Jl

 δEl

δE∗
−l

+
√
2κΛ

 vl(t)

v∗−l(t)

 , (3.10)

where

Jl =

 Rl Sl

S∗
l R∗

l

 (3.11)

35



is a 2× 2 Jacobian matrix. These equations can be solved in the Fourier domain as

 δẼl(ω)

δẼ∗
−l(−ω)

 = −
√
2κΛ[Jl − iωI2]

−1

 ṽl(ω)

ṽ∗−l(−ω)

 , (3.12)

where I2 is a 2nd order identity matrix, and the solution is given by

δẼl(ω) =
√
2κΛ

(R∗
l − iω)ṽl(ω)− Slṽ

∗
−l(−ω)

Dl(ω)
, (3.13)

where

Dl(ω) = (Rl − iω)(R∗
l − iω)− SlS∗

l (3.14)

=
(
κ2 + η2l − v2gγ

2|E0(0)|4 − ω2
)
+ 2iκω ,

ηl = ℑ[Rl] = σ − ζ2
2
l2 + 2vgγ |E0(0)|2 . (3.15)

The spectral density of the fluctuations for the output field in the mode l then can be analytically

calculated as

SδE,l(ω) = 2κeTR

〈
δẼl(ω)δẼ∗

l (ω)
〉

(3.16)

= 4κeκTRΛ
2

ω2 + 2ηlω + κ2 + η2l + v2gγ
2|E0(0)|4(

κ2 + η2l − v2gγ
2|E0(0)|4 − ω2

)2
+ 4κ2ω2

.

Using the above equation and Parseval’s theorem, the variance of the fluctuations can be written

36



(a)

−100 −50 0 50 100

0

−50

−100

l

|E
o
u
t,
l|2

[d
B
m
]

(b)

103 106 109

−250

−200

f [Hz]

S
δ
E,
1
8
[d
B
/H

z]

(c)

−100 −50 0 50 100

−136

−139

−142

l

〈|δ
E o

u
t,
l|2
〉[
dB

]

Figure 3.3: Several computation and simulation results of a Kerr microcomb below the threshold.
The parameters are set to: P

L
= 1.65 mW, σ = −κ, and ζ2 = 2π × 2.9 kHz. (a) Snapshot of the

frequency comb pattern under noise. (b) Spectral density of the sidemode l = 18 obtained after
computation of Eq. (3.16) (red line) and simulation of Eq. (3.1) (blue dots). (c) Average square
amplitude of the fluctuations in each mode: The red line shows the computation results from Eq.
(3.17), while the blue dots show the simulation results from Eq. (3.1).

as

〈
|δEout,l|2

〉
=

1

2π

∫ +∞

−∞
SδE,l(ω) dω

=
2κeTR

Λ2(κ2 + η2l )

κ2 − v2gγ
2|E0(0)|4 + η2l

. (3.17)

Figure 3.3 displays the result of our numerical simulations and they are found to be in excellent

agreement with the analytical results. The above formula is valid as long as the sidemodes re-

main noise-driven, i.e., when the real parts of eigenvalues of the Jacobian Jl are strictly negative

following the condition v2gγ
2|E0(0)|4 < κ2 + η2l .

We note that the noise amplitude parameter Λ is defined such that ⟨|El(t)|2⟩ → Λ2 when√
P

L
→ 0, so Λ2 equals the intracavity power of the sidemodes in the hypothetical case where

they are exclusively driven by external noise. Above the threshold, this trivial solution becomes

unstable and Roll patterns may emerge from the noisy background.

37



3.4 Microcomb Above Threshold: Fluctuations from Pattern Drift and Pattern

Deformation

When the power of the external laser pump field is set above the threshold, some sidemode

pairs can be deterministically excited, thereby leading to Kerr comb formation. These combs

correspond to various spatiotemporal patterns inside the cavity, as shown by Fig. 3.2. Here, we

aim to analyze the dynamics of fluctuations generated by the additive Gaussian white noise on

these patterns.

The first step is to note that the unitary symmetry of Eq. (3.1) indicates that the phase of

the stationary side modes is arbitrary because of the invariance under phase rotation

El(0) → El(0) exp(ilϑ) (3.18)

for all l and with arbitrary ϑ. In other words, the system is neutrally stable in the azimuthal

direction. This corresponds to a shift of the entire azimuthal profile by a constant angle ϑ. As a

result of this invariance, the Jacobian of the coupled mode equations has a null eigenvalue.

The effect of this mode is to induce a pattern drift that undergoes a random walk in the

azimuthal direction. In other words, whenever a pattern is stochastically perturbed in the az-

imuthal direction, the perturbation is undamped and can diverge asymptotically – Indeed, further

analysis will show that this is a Wiener process. This phenomenon is referred to as timing jitter

in systems where the output signal is analyzed as an unwrapped timetrace instead of a stationary

pattern [58, 59, 60]. As a result of this drift, the traditional perturbation methods that assume

small deviations are invalid.

38



To circumvent this problem, we rewrite the field in each mode as

El(t) = El(0) exp [ilϑ(t)] + δEl(t) , (3.19)

where ϑ(t) is the global phase deviation of the comb with ϑ(0) = 0 – Note that the complex-

valued amplitude fluctuations are also initialized as δEl(0). The complex-valued perturbations

δEl are referred to as orbital deviations, and they represent the deviation from the orbit of the

stationary solution [61]. These perturbations are associated to eigenvalues with strictly negative

real parts, and for that reason, are damped. As a consequence, the fluctuations δEl are small at

all time and they lead to pattern deformation, i.e., small stochastic fluctuations of the shape of

the spatiotemporal pattern. The effect of pattern drift and deformation will be evaluated via their

corresponding phase noise spectra. A schematic representation of pattern drift and pattern defor-

mation is displayed in Fig. 3.4.

3.4.1 Optical Stochastic Fluctuations of Solitons

In the case of bright and dark solitons, all of the sidemodes adjacent to the pumped mode

are excited. As a consequence, all the related field fluctuations have to be accounted for, and

they will be globally coupled. The variable El(0) can be numerically calculated by solving the

39



𝜋−𝜋 0 𝜋−𝜋 0

Pattern Drift
(Global Phase Deviation)

Pattern Deformation
(Orbital Deviation)

(a) (b)

Figure 3.4: Schematic illustration of pattern drift and pattern deformation. (a) Schematic illustra-
tion of pattern drift for a soliton. The red curve represents the initial unperturbed soliton, while
the blue curves represent the perturbed solitons with pattern drift. (b) Schematic illustration of
pattern deformation for a soliton. The red curve represents the initial unperturbed soliton, while
the green curves represent the perturbed solitons with pattern deformation. Pattern drift corre-
sponds to a translational change of the soliton’s azimuthal position within the cavity. On the other
hand, pattern deformation corresponds to a change of the soliton’s shape. Our analysis shows that
pattern drift is the leading contribution to phase noise at a low frequency offset from the carrier,
while pattern deformation is the leading contribution to phase noise at high frequency offset.

equations

0 = −κ El(0) + i

[
σ − ζ2

2
l2
]
El(0)

+ivgγ
∑
m,n,p

δ(m− n+ p− l) Em(0)E∗
n(0)Ep(0)

+δ(l)
√

2κe/TR

√
P

L
(3.20)

and choosing an appropriate phase.

Substituting Eq. (3.19) into Eq. (3.1), and ignoring the terms for multiplication of fluctua-

tions, we get the set of 2N + 1 differential equations again:

ilEl(0) exp(ilϑ)ϑ̇+ δĖl (3.21)

=
N∑

p=−N

Rl,p δEp +
N∑

p=−N

Sl,p δE∗
p +

√
2κΛ vl(t)

40



for the fluctuations in each mode, where

Rl,p(t) =

[
−κ+ i(σ − ζ2

2
l2)

]
δ(p− l) (3.22)

+2ivgγ
∑
m,n

δ(m− n+ p− l) Em(0)E∗
n(0) exp [i(l − p)ϑ(t)]

Sl,p(t) = ivgγ
∑
m,n

δ(m+ n− p− l) Em(0)En(0) exp [i(l + p)ϑ(t)] (3.23)

are complex-valued parameters determined by the initial stationary solutions and the global phase

deviation.

We can rewrite Eq. (3.21) under the form of a noise-driven flow, in the form

 X

X∗

 ϑ̇+

 δΨ̇

δΨ̇∗

 = J(t)

 δΨ

δΨ∗

+

 V(t)

V∗(t)

 , (3.24)

where the Jacobian J(t) is given by the following square matrix of order (4N + 2)

J(t) =

 R(t) S(t)

S∗(t) R∗(t)

 , (3.25)

while the (2N + 1)-dimensional vectors involved in this flow are explicitly defined as

δΨ(t) =


δE−N(t)

...

δE+N(t)

 , V(t) =
√
2κΛ


v−N(t)

...

v+N(t)

 , (3.26)

41



and

X(t) =



−iNE−N(0) exp [−iNϑ]

−i(N − 1)E−(N−1)(0) exp [−i(N − 1)ϑ]

...

−iE−1(0) exp [−iϑ]

0

iE1(0) exp [iϑ]
...

i(N − 1)EN−1(0) exp [i(N − 1)ϑ]

iNEN(0) exp [iNϑ]



. (3.27)

It should be noted that the Jacobian matrix J(t) has to be determined numerically because its

components depend on the stationary state values of the semi-classical modal fields El(0). This

Jacobian matrix is diagonalizable as

J(t) = W(t)DP(t) , (3.28)

where W(t) is the matrix composed by the eigenvectors of J(t), D is a constant diagonal matrix

composed of the eigenvalues of J(t), and P(t) is the inverse of W(t). It is important to note that

owing to the symmetry properties of J(t), D is a time-independent matrix, while P(t) and W(t)

are time-dependent.

42



Using Eq. (3.18), one can see one of the eigenvectors of the Jacobian J(t) is

w1(t) = C

 X(t)

X∗(t)

 exp(iα0) , (3.29)

where

C =
1√

2
∑N

l=−N l
2|El(0)|2

(3.30)

is the normalization factor that leads to |w1|2 = 1, while α0 is a real-valued constant. This

solution associated to the eigenvalue 0 is generally referred to as the Goldstone mode. Without

loss of generality, we define w1(t) as the first column vector in matrix W(t), 0 as the first diagonal

element in the matrix D and p1(t) as the first row vector in the matrix P(t).

It now appears that each row of the matrix W(t) and each column of the matrix P(t) have

the same temporal dependence. For example, the time dependence of any element in the first row

of W(t) is the same as the time dependence ofW1,1, which is exp [−iNϑ] as shown in Eq. (3.27).

Correspondingly, the time dependence of any element in the first column of P(t) is the same as

the time dependence of P1,1, which is exp [iNϑ]. Since P(t) is the inverse of W(t), the relations

between these vectors are given by:

pm(t) ·wn(t) =


1 for m = n

0 for m ̸= n

. (3.31)

Substituting Eqs. (3.28) and (3.29) into Eq. (3.24), the dynamical flow can be linearly transformed

43



into the new form

exp(−iα0)

C
P(t)w1(t)ϑ̇+ δΦ̇ = D δΦ+ Z(t) , (3.32)

where we have introduced

δΦ = P(t)

 δΨ

δΨ∗

 =


δϕ1(t)

...

δϕ(4N+2)(t)

 (3.33)

and

Z(t) = P(t)

 V(t)

V∗(t)

 =


z1(t)

...

z(4N+2)(t)

 . (3.34)

The correlation properties of these new stochastic variables can be calculated as

Ξm,n = ⟨z∗m(t)zn(t)⟩ = 2κΛ2

4N+2∑
i=1

〈
P∗

m,iPn,i

〉
. (3.35)

Since all columns in the matrix P(t) have the same dependence on time t, we can see that the

quantity P∗
m,iPn,i is invariant, i.e.

P∗
m,i(t)Pn,i(t) = Constant . (3.36)

Equation (3.32) describes a 4N+2-dimensional set of dynamical equations in the diagonal space.

44



The first one of these equations can be written as

exp(−iα0)

C
ϑ̇+ δϕ̇1 = z1(t) , (3.37)

which is a Wiener process, because the first element of the diagonal matrix D is null. Since δϕ1

represents the weight of the orbital deviation on the stationary solution orbit, we have to select

ϑ̇ = C exp(iα0)z1(t) and δϕ̇1 = 0 , (3.38)

and set the value of δϕ1 to 0. This means that although we initially introduced 4N + 3 variables

(2N + 1 orbital deviations El, their 2N + 1 conjugates, and global phase deviation ϑ), the total

number of degrees of freedom for this system is still 4N + 2.

The other dynamical equations embedded in Eq. (3.32) describe a set of (4N + 1) Orn-

stein–Uhlenbeck processes driven by additive Gaussian white noises. The correlation properties

of the linearly transformed modes are given by

Πm,n = ⟨δϕ∗
m(t)δϕn(t)⟩ = 0 (3.39)

for the cases that m = 1 or n = 1, and

Πm,n = ⟨δϕ∗
m(t)δϕn(t)⟩ = − Ξm,n

D∗
m,m +Dn,n

(3.40)

for any other cases, where Di,i is the ith diagonal element of the matrix D.

The orbital deviation field δEl(t) can now be recovered from δΦ via an inverse linear trans-

45



(a) Comb 1

−200 −100 0 100 200
−110

−105

−100

−95

l

⟨|δ
E l
|2 ⟩

[d
B
]

(b) Comb 2

−200 −100 0 100 200
−110

−105

−100

−95

l

⟨|δ
E l
|2 ⟩

[d
B
]

(c) Comb 3

−30 −15 0 15 30
−110

−105

−100

−95

l

⟨|δ
E l
|2 ⟩

[d
B
]

(d) Comb 4

−30 −15 0 15 30
−110

−105

−100

−95

l

⟨|δ
E l
|2 ⟩

[d
B
]

Figure 3.5: Variance of orbital deviation of the intra-cavity field. The figure shows the results of
(a) bright soliton (Comb 1), (b) dark soliton (Comb 2), (c) roll pattern of order L = 1 (Comb
3), and (d) roll pattern of order L = 4 (Comb 4) displayed in Fig. 3.2 as a function of the
sidemode order l, respectively. The thick gray curves correspond to analytical results obtained
from Eqs. (3.41), (3.64), and (3.65), while the thin red curves result from the simulations of
Eq. (3.1). One can note the excellent agreement between analytical and numerical results.

formation, and the variance of this fluctuation is given by

〈
|δEl(t)|2

〉
=

∑
m,n

W∗
l+N+1,mWl+N+1,nΠm,n . (3.41)

The results of this orbital deviation analysis in solitons are presented in Fig. 3.5(a-b) and display

an excellent agreement with numerical simulations.

3.4.2 Optical Stochastic Fluctuations of Rolls

For roll patterns of order L, the main steps of the previous analysis as developed for solitons

remain valid. Accounting for the fact that not all of the sidemodes are excited, the modal fields

of the roll pattern can be written as

El(t) =


El(0) exp(ilϑ(t)) + δEl(t) for l = nL

δEl(t) for l ̸= nL

, (3.42)

46



where n is an integer. After neglecting higher-order stochastic terms, we obtain an equation

analogous to Eq. (3.24), but now the (4N + 2)-dimensional Jacobian

J(t) =

 R(t) S(t)

S∗(t) R∗(t)

 (3.43)

is defined with the complex-valued elements

Rl,p(t) =

[
−κ+ i(σ − ζ2

2
l2)

]
δ(p− l) (3.44)

+2ivgγ
∑
m,n

δ(mL− nL+ p− l) EmL(0)E∗
nL(0)

× exp [i(l − p)ϑ(t)] ,

Sl,p(t) = ivgγ
∑
m,n

δ(mL+ nL− p− l) EmL(0)EnL(0)

× exp [i(l + p)ϑ(t)] , (3.45)

which are determined by the deterministic (ideally noiseless) stationary roll pattern.

Equations (3.44) and (3.45) indicate that the Jacobian becomes block diagonal since only

sidemodes l that are multiples of L are excited, and we can therefore reorder this Jacobian to

make the qth block Jq correspond to the Bloch modes with wavenumber q, namely δEpL±q with p

being an integer [8]. Without loss of generality, we consider M excited modes on each side and

the non-excited adjacent modes, which includes the mode l when it fulfills the condition

−(ML+ L/2) < l ≤ML+ L/2 . (3.46)

If L is an odd number, there are (L + 1)/2 blocks in the Jacobian matrix, with the label q being

47



from 0 to (L − 1)/2. The Jacobian J0 is a matrix of order 2(2M + 1), and the other blocks are

matrices of order 4(2M + 1). If L is an even number, there are L/2 + 1 blocks in the Jacobian

matrix (from 0 to L/2), where J0 and JL/2 are matrices of order 2(2M +1), and the other blocks

are matrices of order 4(2M + 1).

We can now treat the equations associated with the various blocks separately. Equa-

tion (3.21) can therefore be rewritten under the form of several noise-driven flows, in the form

δ(q)

 X0

X0
∗

 ϑ̇+

 δΨ̇q

δΨ̇∗
q

 = Jq

 δΨq

δΨ∗
q

+

 Vq(t)

V∗
q(t)

 .

(3.47)

From these equations, we can obtain vectors of order 2(2M + 1) for the field and external noise

fluctuations either for the case q = 0 or for the case q = L/2 when L is even, in the form

δΨq(t) =


δE−(ML+q)(t)

...

δEML+q(t)

 , Vq(t) =
√
2κΛ


v−(ML+q)(t)

...

vML+q(t)

 . (3.48)

For all the other cases, the 4(2M + 1)-dimensional fluctuation and noise vectors are explicitly

48



expressed as

δΨq(t) =



δE−(ML+q)(t)

...

δEML+q(t)

δE−(ML−q)(t)

...

δEML−q(t)



, Vq(t) =
√
2κΛ



v−(ML+q)(t)

...

vML+q(t)

v−(ML−q)(t)

...

vML−q(t)



. (3.49)

The Jacobian and all the blocks Jq(t) are diagonalizable. Hence, they can be written as

Jq(t) = Wq(t)Dq Pq(t) , (3.50)

where Dq is a constant diagonal matrix, while Wq(t) and Pq(t) are time-dependent.

The real part of the eigenvalues for every block is always negative, except for the one for

the 0th block. We can see that 0 is also one of the eigenvalues of J0, and the corresponding

eigenvector w(0)
1 can be written as the following vector of order 4M + 2:

w
(0)
1 = C

 X0(t)

X∗
0(t)

 exp(iα0) , (3.51)

49



where

X0(t) =



−iMLE−ML(0) exp [−iMLϑ]

−i(M − 1)LE−(M−1)L(0) exp [−i(M − 1)Lϑ]

...

−iLE−L(0) exp [−iLϑ]

0

iLEL(0) exp [iLϑ]
...

i(M − 1)LE(M−1)L(0) exp [i(M − 1)Lϑ]

iMLEN(0) exp [iMLϑ]


(3.52)

is defined in a similarly way as Eq. (3.27), while

C =
1√

2
∑M

m=−M m2L2|EmL(0)|2
(3.53)

is the normalization factor with α0 being a real-valued constant phase. We still assume that w(0)
1

is the first column vector in the matrix W0, where W0 is the inverse of P0. The first diagonal

element in matrix D0 is null, and the first row vector in matrix P0 is p(0)
1 .

We can now process each block in the same way as in the previous subsection. In the

diagonal space, the 0th flow is linearly transformed into

exp(−iα0)

C
P0(t)w

(0)
1 (t)ϑ̇+ δΦ̇0 = D0 δΦ0 + Z0(t) , (3.54)

50



and the other flows (q = 1, 2, ...) are transformed into

δΦ̇q = Dq δΦq + Zq(t) , (3.55)

with

δΦq = Pq(t)

 δΨq

δΨ∗
q

 =


δϕq,1(t)

δϕq,2(t)

...

 , (3.56)

and

Zq(t) = Pq

 Vq(t)

V∗
q(t)

 =


zq,1(t)

zq,2(t)

...

 . (3.57)

In order to make orbital deviations orthogonal to the orbit of the stationary solution, the first

dynamical equation in Eq. (3.54) leads to

ϑ̇ = C exp(iα0)z0,1(t) (3.58)

and

δϕ0,1 = 0 . (3.59)

51



The correlation properties of the linearly transformed stochastic variables in Zq(t) are:

Ξ(q)
m,n =

〈
z∗q,m(t)zq,n(t)

〉
= 2κΛ2

∑
i

〈
P(q)∗

m,i P(q)
n,i

〉
, (3.60)

with

P(q)∗
m,i (t)P(q)

n,i (t) = P(q)∗
m,i (0)P(q)

n,i (0) . (3.61)

As a consequence, the correlation properties of the linearly transformed modes are obtained as

Π(0)
m,n =

〈
δϕ∗

q,m(t)δϕq,n(t)
〉
= 0 (3.62)

for the cases that m = 1 or n = 1, and as

Π(q)
m,n =

〈
δϕ∗

q,m(t)δϕq,n(t)
〉
= − Ξ

(q)
m,n

D(q)∗
m,m +D(q)

n,n

(3.63)

for the other cases, where D(q)
i,i is the ith diagonal element of Dq.

Since the eigenvalue 0 only belongs to the 0th block, the Goldstone mode only influences

the correlation properties of modes dominated by J0, which are the excited sidemodes in the

stationary solution. Although the global phase deviation ϑ(t) also affects all the other block

matrices, it does not affect the value of the correlations.

After recovering δEl(t) from δΦ via an inverse linear transformation, we get the average

52



variance of the optical field fluctuations as

〈
|δEkL+q(t)|2

〉
=

∑
m,n

W(q)∗
k+M+1,mW

(q)
k+M+1,nΠ

(q)
m,n (3.64)

for the (kL+ q)th mode, and

〈
|δEkL−q(t)|2

〉
=

∑
m,n

W(q)∗
k+3M+2,mW

(q)
k+3M+2,nΠ

(q)
m,n (3.65)

for the (kL − q)th mode. Equation (3.65) only works for the case where q is a positive integer

smaller than L/2, while Eq. (3.64) also works for the case with q = 0 or q = L/2. The results of

this orbital deviation analysis for roll patterns are presented in Fig. 3.5(c-d).

3.5 Optical Power Spectra

The applications of Kerr microcombs are essentially related to their metrological perfor-

mance. The spectral purity of the comb and of the microwaves that can be generated with them

is generally evaluated in terms of phase noise spectra. If we assume that the stationary output

signal is

Eout,l(0) = |Eout,l(0)| exp(iψs,l) , (3.66)

53



the effect of the additive Gaussian white noise is that the instantaneous output optical field exiting

the resonator now becomes

Eout,l(t) = [|Eout,l(0)|+ δEout,l(t)] exp[i(ψs,l + δψl(t))] , (3.67)

where δEout,l(t) is the real-valued deviation from the initial stationary amplitude and δψl(t) is

the deviation from the initial stationary phase. The deviation δψl(t) is therefore defined as the

phase noise and its spectrum is referred to as the phase noise spectrum.

We hereafter show that the overall optical phase noise is the sum of two contributions:

The first one originates from slow pattern drift and is dominant at low frequency offset from the

carrier; the second one is induced by fast pattern deformation and is dominant at high frequency

offset from the carrier.

3.5.1 Contribution of Pattern Drift

We first determine the optical phase noise contribution due to pattern drift, also referred to

as global phase deviation. From the previous discussion and calculations, we know this global

phase deviation originates from the Goldstone mode with eigenvalue 0. This mode induces a free

azimuthal motion for the pattern and, as shown in Eq. (3.38), the dynamics of the global phase

deviation of the soliton is governed by

ϑ̇ = C exp(iα0)p1 ·V(t) . (3.68)

54



This equation is a Wiener process, which here corresponds to a Brownian motion for the phase of

the pattern. It is known that in that case the spectrum of the global phase deviation is expressed

as [62]

Sϑ(ω) =
C2Ξ1,1

ω2
. (3.69)

On the other hand, the phase spectrum of the mode l in the comb is expressed as

LE,l(ω) = l2Sϑ(ω) =
l2C2Ξ1,1

ω2
, (3.70)

while the corresponding power spectral density of the output field Eout,l is

SE,l(ω) =
l2C2Ξ1,1|Eout,l(0)|2
(l2C2Ξ1,1)2/4 + ω2

. (3.71)

When the additive noise is very small, the phase noise spectrum is related to the power spectrum

by

LE,l(ω) ≃
SE,l(ω)

|Eout,l(0)|2
(3.72)

except when the offset frequency is too close to the carrier. Meanwhile, since both spectra have

very large values near the carrier frequency, they are expected to be accurate at low offset fre-

quencies.

Since the effect of Goldstone mode is limited to the fluctuations dominated by the 0th block

of the Jacobian for roll patterns, the solutions of the phase spectrum and power spectral density

55



can be calculated as

LE,l(ω) =
l2C2Ξ

(0)
1,1

ω2
(3.73)

and

SE,l(ω) =
l2C2Ξ

(0)
1,1|Eout,l(0)|2

(l2C2Ξ
(0)
1,1)

2/4 + ω2
≃ |Eout,l(0)|2LE,l(ω) , (3.74)

which are exactly the same as the one previously obtained after replacing Ξ1,1 with Ξ
(0)
1,1.

3.5.2 Contribution of Pattern Deformation

We now focus on the optical phase noise contribution due to pattern deformation induced

by the sidemode fluctuations. Although the optical phase noise can be obtained from Eq. (3.21)

numerically, we aim at deriving an analytical formula by only accounting for the interaction

terms between mode l and 0, while neglecting the remaining ones. We therefore assume that the

dynamical equations are

δĖl =
[
−κ+ i

(
σ − ζ2

2
l2 + 2vgγ |E0(0)|2

)]
δEl (3.75)

+
√
2κΛ vl(t) .

This approximation is valid because |E0(0)| ≫ |El(0)| for l ̸= 0, and therefore dominates the

entire high frequency spectrum. This approximation also indirectly assumes that the modal fluc-

tuations are independent of each other, i.e., uncorrelated.

56



We can write the lth nonzero stationary output field as

Eout,l(0) = |Eout,l(0)| exp(iψs,l) (3.76)

and the corresponding fluctuation as

δEout,l(t) = |δEout,l(t)| exp(iψf,l(t)) . (3.77)

As a consequence, the optical phase noise contribution due to pattern deformation at time t can

be approximately written as the projection of δEout,l(t) in the vertical direction of Eout,l(0), in the

form

ψl(t) ≃
|δEout,l(t)|
|Eout,l(0)|

sin(ψf,l − ψs,l) . (3.78)

The Fourier transform of the projection |δEout,l(t)| sin(ψf,l−ψs,l) can now be solved analytically,

and we can obtain the spectrum of ψl, in the form

LE,l =
2κeκTR

Λ2

|Eout,l(0)|2
κ2 + η2l + ω2

(κ2 + η2l − ω2)2 + 4κ2ω2
. (3.79)

The corresponding power spectral density of the field Eout,l may be derived from Eq. (3.75) as

SE,l(ω) =
4κeκTR

Λ2

κ2 + (ω − ηl)2
. (3.80)

One can note that there is a term of |Eout,l(0)|2δ(ω) because of the direct current part of Eout,l(t).

57



Although most of the additive noise drives the orbital deviation, we know that the pattern drift

noise is large at low-frequency offset and this model is only intended to be valid for high-

frequency offset. Therefore, this zero-frequency term is omitted in the above equation. This

approach is valid for both soliton and roll patterns.

3.5.3 Comparison Between Numerical and Analytical Results

In the numerical simulations, the global phase deviation ϑ(t) is tracked in the calculation by

numerical integration of Eq. (3.68), while the phase noise is derived by calculating the difference

between the phase of the output signal and the phase of the initial stationary solution. The phase

noise spectra are simulated with a Hanning window

Hw(t) =
√
8/3 sin2 (πt/T ) , (3.81)

where T is the temporal width of the window [62].

Figure 3.6 displays the excellent agreement between our analytical predictions and the nu-

merical simulations. It can be seen that for small frequency offsets, the numerical simulations

agree well with the pattern drift contribution of phase noise. Conversely, for large frequency off-

sets, the numerical simulations agree instead with the pattern deformation contribution of phase

noise. The pattern drift sidemode spectra in Eq. (3.74) depends on the sidemode amplitude, while

the pattern deformation sidemode spectra in Eq. (3.80) does not. As a consequence, there is a

notable difference when we compare the spectra for different sidemodes of the roll pattern. How-

ever, for solitons, the amplitude of the sidemode decreases slower with the increase of |l|, so that

the change is not so large.

58



(a)

Comb 1

−250

−200

−150
S
E,
1
[d
B
ra
d2
/H

z]

(b)

Comb 2

−250

−200

−150

S
E,
1
[d
B
ra
d2
/H

z]

(c)

Comb 3

−250

−200

−150

S
E,
1
[d
B
ra
d2
/H

z]

(d)

Comb 4

−250

−200

−150

S
E,
4
[d
B
ra
d2
/H

z]

(e)

Comb 1

102 104 106 108

−250

−200

−150

Frequency Offset [Hz]

S
E,
4
[d
B
ra
d2
/H

z]

(f)

Comb 2

102 104 106 108

−250

−200

−150

Frequency Offset [Hz]

S
E,
4
[d
B
ra
d2
/H

z]

(g)

Comb 3

102 104 106 108

−250

−200

−150

Frequency Offset [Hz]

S
E,
4
[d
B
ra
d2
/H

z]

(h)

Comb 4

102 104 106 108

−250

−200

−150

Frequency Offset [Hz]

S
E,
1
6
[d
B
ra
d2
/H

z]

Figure 3.6: Sidemode spectra of the output optical field. The figure corresponds to (a, e) bright
soliton (Comb 1), (b, f) dark soliton (Comb 2), (c, g) roll pattern of order L = 1 (Comb 3) and (d,
h) roll pattern of order L = 4 (Comb 4) displayed in Fig. 3.2. (a-d) Spectra of excited sidemodes
that are closest to the right of the center mode. (e-f) Spectra of excited sidemodes that are fourth
closest to the right of the center mode. The blue lines correspond to the noise contribution of
pattern drift [the computation of Eq. (3.74)]. The green lines correspond to the noise contribution
of pattern deformation [the computation of Eq. (3.80)]. The red lines correspond to the spectra
obtained after the numerical simulation of the stochastic modal equations. These simulations
clearly indicate that pattern drift noise dominates the spectra for low offset frequencies, while
pattern deformation noise dominates the spectra for high offset frequencies.

3.6 Microwave Phase Noise Spectra

Microwave generation using Kerr combs is obtained via a process of photodetection. The

spectral purity of these microwaves is generally evaluated in terms of phase noise spectra [7, 63,

64, 65, 66]. For the sake of simplicity, we here consider an ideal photodetector, so that the noise

properties of the generated microwave can be solely attributed to the comb.

The photodetector extracts the power envelope of the output field and generates a radio-

frequency (RF) signal that is proportional to the incoming optical power (in units of Volts) fol-

lowing

VPD(t) = S|Eout|2 , (3.82)

59



where S is the responsivity of the photodetector (in units of VW−1), while

Eout =
∑
l

Eout,l exp (ilΩR
t) . (3.83)

This signal can be Fourier-expanded in the form

VPD(t) =
1

2
M0 +

+∞∑
n=1

[
1

2
Mn exp(inΩR

t) + c.c.

]
, (3.84)

where the spectral components

Mn = 2S
∑
m

E∗
out,mEout,m+n (3.85)

are the complex-valued envelopes of the microwave harmonics of frequency n × Ω
R

, and c.c.

stands for the complex conjugate [67, 68].

3.6.1 Contribution of Pattern Drift

We can now evaluate the effect of pattern drift on microwave phase noise. The slowly

varying envelope of the microwave spectral component of frequency nΩ
R

in a noiseless situation

can be written as

Mn(0) = 2S
∑
l

Eout,l+n(0)E∗
out,l(0) (3.86)

60



with

Eout,l(0) =
√
2κeTR

El(0)−
√
P

L
δ(l) . (3.87)

Pattern drift induced by global phase deviation modifies the slowly-varying envelope of the nth

microwave component as

Mn = Mn(0) exp [inϑ(t)] . (3.88)

For bright and dark solitons, the parameter n could be any positive integer and the spectrum of

the microwave phase noise is

LM,n(ω) = n2Sϑ(ω) =
n2C2Ξ1,1

ω2
, (3.89)

while the corresponding power spectrum is

SM,n(ω) =
n2C2Ξ1,1|Mn(0)|2
(n2C2Ξ1,1)2/4 + ω2

≃ |Mn(0)|2LM,n(ω) . (3.90)

For roll patterns, n is restricted to be a positive multiple of the order L. Similar to our previous

analysis of the sidemode spectrum, the formulas for the drift-induced microwave spectrum and

the corresponding power spectrum can be derived as

LM,n(ω) = n2Sϑ(ω) =
n2C2Ξ

(0)
1,1

ω2
(3.91)

61



and

SM,n(ω) =
n2C2Ξ

(0)
1,1|Mn(0)|2

(n2C2Ξ
(0)
1,1)

2/4 + ω2
≃ |Mn(0)|2LM,n(ω) , (3.92)

which are the same equations we derived before, except that the parameter Ξ1,1 is replaced by

Ξ
(0)
1,1.

3.6.2 Contribution of Pattern Deformation

With respect to the effect of pattern deformation on microwave phase noise, one can first

note that the effect of modal fluctuations on the slowly-varying envelope of the microwave spec-

tral component of frequency nΩ
R

is expressed as

Mn = 2S
∑
l

Eout,l+nE∗
out,l (3.93)

≃ 2S
∑
l

[
Eout,l+n(0)δE∗

out,l + δEout,l+nE∗
out,l(0)

]
+Mn(0) .

Using the formula for optical phase noise spectra as derived in Eq. (3.79), we obtain the following

formula for the microwave phase noise spectrum:

LM,n(ω) =
8κeκTR

Λ2S2

|Mn|2
× (3.94)∑

l

[
|Eout,l+n(0)|2 + |Eout,l−n(0)|2

] (κ+ η2l + ω2)

(κ2 + η2l − ω2)2 + 4κ2ω2
.

62



while the corresponding power spectrum is given by

SM,n(ω) = 16κeκTR
Λ2S2 × (3.95)∑

l

[ |Eout,l+n(0)|2
κ2 + (ω + ηl)2

+
|Eout,l−n(0)|2
κ2 + (ω − ηl)2

]
.

The term of |Mn(0)|2δ(ω) is omitted in the above equation because this model is intended to

describe high frequency fluctuations. Here again, these equations are valid for both solitons and

roll patterns. For solitons, n could be any positive integer, while for roll patterns, n must be a

positive multiple of L.

3.6.3 Comparison Between Numerical and Analytical Results

The numerical simulations of microwave phase noises are performed at the same time with

the optical phase noise, after taking into account the relation between the microwave and optical

modes in Eq. (3.82).

The solution of the simulation and analytical results for the generated microwaves are

shown in Fig. 3.7. It can be seen that the pattern drift spectra are always very accurate at low

offset frequencies because the global phase deviation dominates in that frequency range. On the

other hand, the numerical simulations converge to pattern deformation spectra at high frequencies

because it is the range corresponding to the fast field fluctuations. The pattern deformation spec-

tra model predicts the potential existence of one or several peaks near the frequency of overall

detuning ηl while the pattern drift spectra do not have any. Similarly to Fig. 3.6, by compar-

ing phase noise spectra for different microwave components in the same comb, we can see the

difference is larger in roll patterns.

63



(a)

Comb 1

−200

−150

−100
L M

,1
[d
B
c/
H
z]

(b)

Comb 2

−200

−150

−100

L M
,1
[d
B
c/
H
z]

(c)

Comb 3

−200

−150

−100

L M
,1
[d
B
c/
H
z]

(d)

Comb 4

−200

−150

−100

L M
,4
[d
B
c/
H
z]

(e)

Comb 1

102 104 106 108

−200

−150

−100

f [Hz]

L M
,4
[d
B
c/
H
z]

(f)

Comb 2

102 104 106 108

−200

−150

−100

f [Hz]

L M
,4
[d
B
c/
H
z]

(g)

Comb 3

102 104 106 108

−200

−150

−100

f [Hz]

L M
,4
[d
B
c/
H
z]

(h)

Comb 4

102 104 106 108

−200

−150

−100

f [Hz]

L M
,1
6
[d
B
c/
H
z]

Figure 3.7: Phase noise spectra of the microwaves. The figure corresponds to the results obtained
via photodetection of (a, e) bright soliton (Comb 1), (b, f) dark soliton (Comb 2), (c, g) roll
pattern of order L = 1 (Comb 3)