# Synchronization of silicon thermal free-carrier oscillators

Gustavo de O. Luiz  Instituto de Física Gleb Wataghin, Universidade Estadual
de Campinas, Rua Sergio Buarque de Holanda, 777, Campinas-SP, Brasil  nanoFAB
Centre, University of Alberta, Edmonton, Alberta, T6G 2V4, Canada  Caique C.
Rodrigues  Thiago P. M. Alegre  Instituto de Física Gleb Wataghin,
Universidade Estadual de Campinas, Rua Sergio Buarque de Holanda, 777,
Campinas-SP, Brasil  Gustavo S. Wiederhecker<EMAIL_ADDRESS>Instituto de
Física Gleb Wataghin, Universidade Estadual de Campinas, Rua Sergio Buarque de
Holanda, 777, Campinas-SP, Brasil

###### Abstract

Recent exploration of collective phenomena in oscillator arrays has
highlighted its potential for accessing a range of physical phenomena, from
fundamental quantum many-body dynamics [1] to the solution of practical
optimization problems using photonic Ising machines [2, 3]. Spontaneous
oscillations often arise in these oscillator arrays as an imbalance between
gain and loss. Due to coupling between array individuals, the spontaneous
oscillation is constrained and lead to interesting collective behavior, such
as synchronized oscillations in optomechanical oscillator arrays [4],
ferromagnetic-like coupling in delay-coupled optical parametric oscillators
[5] and binary phase states in coupled laser arrays [6]. A key aspect of
arrays is not only the coupling between its individuals but also their
compliance towards neighbor stimuli. One self-sustaining photonic oscillator
that can be readily implemented in a scalable foundry-based technology is
based on the interaction of free-carriers, temperature and optical field of a
resonant silicon photonic microcavity [7, 8, 9, 10, 11]. Here we demonstrate
that these silicon thermal-free-carrier oscillators are extremely compliant to
external excitation and can be synchronized up to their 16${}^{\text{th}}$
harmonic using a weak seed. Exploring this unprecedented compliance to
external stimuli, we also demonstrate robust synchronization between two
thermal free carrier oscillators.

##  I Introduction

Coupled oscillator phenomena plays an important role in physical systems, for
instance, coupled oscillator arrays may lay out the synchronous behavior of
complex biological systems [12], or even emulate the behavior of coupled spins
in solids that quickly solves complex optimization problems at low power
consumption [13]. The dramatic evolution of micro and nanofabrication
techniques has enabled the demonstration of these coupled oscillators using
distinct physical systems such as lasers, optical parametric oscillators [2],
and optomechanical oscillators [14, 15]. In these examples, the oscillation
emerges as a spontaneous process resulting from the imbalance between gain and
loss around a resonant optical mode. When building coupled arrays of these
oscillators, a major challenge is to ensure that their dissimilar resonant
frequencies match each other – a constraint that often requires control loops
[16] or precise individual tuning. This is a consequence of their rather stiff
oscillation frequencies that hardly deviate from the bare frequencies set by
geometrical boundary conditions [14, 17].

More recently, a simple yet powerful oscillator has been explored in silicon,
based on the interaction between the optical field, free charge carriers (FC),
and temperature [18, 19, 8]. The resulting oscillator frequency is not only
very dependent on the intra-cavity energy, but is also very nonlinear, easily
generating tens of harmonics [8]. Such compliance of the oscillation frequency
to intra-cavity energy has been recently explored in the synchronization of
this oscillator to an optomechanical oscillator [20].

In this work, we experimentally demonstrate that silicon free-carrier thermal
oscillators exhibit exceptional compliance to external excitations and can be
injection-locked up to their 16${}^{\text{th}}$ harmonic using a weak
modulation of the continuous wave pump laser. Leveraging this remarkable
compliance, we also demonstrated the synchronization of a pair of coupled
oscillators through their evanescent optical field. Furthermore, the
compatibility of these oscillators with foundry-based CMOS technology should
facilitate their scalability to large oscillator arrays.

##  II Thermal free carrier oscillators

Although silicon is transparent for light with wavelengths longer than
approximately 1.1 µm [21], free-carriers (FC) are generated through two-photon
absorption (TPA) [22, 23] at large optical intensities, leading to changes in
the real and imaginary parts of the refractive index [24]. In an optical
cavity, an increase in FC density (N) causes a blue shift of the resonance
frequency ($\delta\omega$) due to a reduction in the refractive index ($\delta
n$),as well as a nonlinear optical loss. Additionally, the temperature
increase ($\theta$) resulting from linear and nonlinear optical absorption
causes an increase in the refractive index due to the thermo-optic effect
[25], thereby red-shifting the optical resonance.

Figure 1: a) Schematic of the suspended silicon resonator, the colormap
represent the optical field that generates heat and free charge carriers
(represented by the tiny blue circles), higher field intensity in the red
shading; b) Pictorial model of thermal-free carrier dynamics and closed loop
diagram of the oscillating system. $P_{in}$: Input laser power, $P_{out}$:
Output laser power; c) Simulated time trace of the transmitted signal; blue
(red) parts indicate when the optical resonance is suffering blue (red) shifts
in the cycle; d) Simulated dynamical FC density versus temperature. The closed
cycle indicates self-sustained oscillations, the arrows indicate the time
evolution.

As indicated in fig. 1(a), the optical field is concentrated on the disk’s
edge, where it is absorbed and generates FC and heat. The generated FC have
typical life-times on the order of nanoseconds [26, 27], most likely due to
surface facilitated recombination [28]. Heat, conversely, is dissipated mainly
through the supporting silicon oxide pedestal, which has poor heat conduction
characteristics, leading to microsecond scale thermal lifetimes. These two
effects cause opposite refractive index changes, which results in resonance
frequency shifts. This frequency shift alters the amount of energy stored
inside the cavity, closes the feed-back loop represented in fig. 1(b) and
amplifies small thermal fluctuations ($k_{b}T$, where $k_{b}$ is the Boltzmann
constant and $T$ ambient temperature), which causes the self-pulsing of the
stored optical energy encoded into the cavity transmitted light.

The model for such interaction in a single resonator is well known [29, 8, 11]
and reproduces all qualitative features of the transmitted signal time-traces.
An example of the simulated time-dependent transmission is shown in fig. 1(c)
(see supplementary information for details). The two distinct time-scales, due
to temperature and FC relaxations, are highlighted as distinct colors in the
traces shown in fig. 1(c,d); there is a fast FC dominated region (blue),
followed by a slow thermal regime (red). Such interplay can be more clearly
seen by plotting the resulting FC density versus temperature variation, which
defines a closed path in fig. 1(d), indicating that self-sustained oscillation
regime was reached.

The closed trajectory in fig. 1(d) also allows us to understand the different
aspects of this oscillation dynamics. For a pump laser ($P_{in}$) blue-detuned
relative to the cavity optical resonance and above the oscillation power
threshold, there is a rapid increase in the FC density that blue-shifts the
resonance. This blue-shift is large enough to leave the laser red-detuned
relative to the instantaneous optical resonance frequency, resulting in the
narrow dip observed in fig. 1(c). At this point in the cycle the cavity
temperature slowly increases, red-shifting the optical resonance frequency,
until the laser is again blue-detuned, forming the second broader dip in the
transmission trace, and the cycle is restarted.

The optical properties of the two microdisk devices explored in this work are
characterized in fig. 2(b) and fig. 2(c). For the single resonator device (2
µm radius, 220 nm thick silicon), as shown in fig. 2(b), the optical resonance
is centered at 1471.4 nm, with a loaded linewidth of $\kappa/2\pi=2.5$ GHz
($Q=81\times 10^{3}$). Due to coupling of counter-propagating modes caused by
surface roughness [30], two resonance peaks are observed split by
$\beta/2\pi=13.6$ GHz. Similarly, the coupled-resonator device (2 µm radii,
220 nm thick silicon, 300 nm coupling gap) has resonances centered at 1472.1
nm and shows 4 distinct peaks. The pair of peaks at each detuning sign are
also due to coupling between counter-propagating modes ($\beta/2\pi=10$ GHz),
while the splitting between the left and right-hand peak pairs is due to the
evanescent coupling between the two cavities ($J/2\pi\approx 28$ GHz). The
fitted linewidth for the coupled cavities is found to be $\kappa/2\pi=2.9$ GHz
($Q=70\times 10^{3}$).

Figure 2: a) Simplified schematics of the experimental setup. TL: tunable
laser; EOM: electro-optical amplitude modulator; VOA: variable optical
attenuator; PC: polarization control; ESA: electrical spectrum analyzer;
Scope: oscilloscope; b) Measured optical resonance of a single cavity. Center
of resonances is at 1471.43 nm, with $\kappa/2\pi=2.5$ GHz linewidth and
$\beta/2\pi=13.6~{}GHz$ splitting due to surface roughness; c) Measured
optical resonance of the coupled resonator system. Center of the graph is at
1472.1 nm with linewidth $\kappa/2\pi=2.9$ GHz, intra-coupling strenght
$\beta/2\pi=10~{}GHz$ and inter-coupling $J/2\pi\approx 28$ GHz. Inset shows
false-colored SEM top view of each device, scale bars are 2 µm long.

##  III Injection locking of a single oscillator

We investigate the compliance of these oscillators by injection locking them
to a weak external signal (at a frequency $f_{\text{inj}}$) that was imposed
on the amplitude of the pump laser. To achieve this modulation, an electro-
optical amplitude modulator (EOM) is inserted in the light path before the
resonators. A simplified schematic of the experimental setup is shown in fig.
2(a), consisting of a pump laser (TL, for tunable laser), a variable optical
attenuator (VOA), a polarization controller (PC) and a photo-detector. The
pump laser is fiber coupled, and a tapered fiber is used to evanescently
couple the light into the cavity [31]. The transmitted light is detected, and
the slow (DC) component of the photocurrent monitors the average cavity
transmission while the fast (AC) components are sent to an oscilloscope and to
an electrical spectrum analyzer (ESA). Additionally, we used a calibrated
acetylene gas cell and a Mach-Zehnder interferometer for calibrating the laser
frequency, which are not shown in the figure.

åç

Figure 3:  Injection locking of a single resonator. a) Single-cavity
transmission radio-frequency spectrogram as the continuous-wave laser
wavelength is swept around an optical resonance, as measured by an ESA; b) DC
optical transmission measured simultaneously with the spectrogram in part (a).
Blue-shaded region indicates thermal free-carrier self-sustaining
oscillations; c) Single-cavity transmission radio-frequency spectrogram,
similar to (a), but with an amplitude-modulated pump serving as an injection
signal at $f_{\text{inj}}=80$ MHz; d) Raw spectogram obtained at fixed pump
wavelength (1473 nm) while the injection frequency ($f_{\text{inj}}$) is
swept. The red-dashed triangle region highlights the high-harmonic
synchronization region, where the injection frequency is larger than the
oscillator frequency, the vertical black-dashed line indicates the fundamental
frequency of the oscillator. e) The same spectrogram as in (d), but with the
frequency axis rescaled by the injection signal frequency. In this rescaled
diagram, the red-dashed triangle region of (d) is mapped into the rectangular
region, the vertical black-dashed line indicates the fundamental oscillator
tone; thin vertical-dashed lines are guides-to-the-eye for the harmonic
relations ($p/q$); f) Phase noise measurement at the fundamental oscillator
tone while it is free-running and injection-locked at the first 5 harmonics,
the corresponding regions where phase-noise is measured is highlighted by the
colored rectangles in the lower right of (d); g) First harmonic spectrum when
system is free-running and locked at the first 5 harmonics. The colors match
part (f).

A simple laser wavelength sweep across the optical resonance of the single
resonator device, without the amplitude modulation, already reveals the
oscillation dynamics described in fig. 1. With an optical input power of only
310 µW, as the laser reaches a threshold detuning (highlighted in blue in the
transmission plot of fig. 3(b)), the radio-frequency spectrogram captured by
the ESA reveals multiple harmonics of the transmitted signal (x-axis is the
ESA calibrated frequency). Moreover, a large optical spring effect is observed
as the frequency of the oscillator varies with the laser wavelength (y-axis in
fig. 1(a,c)), which is a hint that the pulses period and duty-cycle (ratio
between duration and period) strongly depend on the optical energy stored in
the resonator [8].

Injection locking can be readily observed by driving the EOM with a weak
$f_{\text{inj}}=80$ MHz sinusoidal electrical drive (drive power of -10 dBm,
modulation depth of 8%). In the presence of the modulation signal, the laser
sweep through the resonance results in a dramatically distinct RF spectrogram,
as shown in fig. 3(c). Although a wide spring effect is still observed, as in
fig. 3(a), the resonator frequency evolves in discrete steps, each
corresponding to non-trivial fractional sub-harmonics of the RF tone, forming
a sequence of the so-called devil staircases [32]. The first step of the devil
staircase is observed while the oscillator first harmonic is at 5 MHz, which
implies that it is being entrained by a modulation at its 16${}^{\text{th}}$
harmonic.

To explore the injection locking behavior in greater detail, we presented the
RF spectrogram in fig. 3(d), with the x-axis representing the ESA calibrated
frequency and the y-axis representing the linearly swept modulation injection
frequency from zero to 160 MHz, while the laser is fixed at a certain
wavelength. The spectrogram in fig. 3(d) reveals multiple synchronization
regions, identifiable by the sloping harmonic lines to the right, rather than
vertical lines. Such sloped harmonics indicate the fractional ratios where the
oscillator is entrained. To simplify the identification of the fractional
ratios, we normalized the horizontal frequency axis by dividing it by the
injection frequency (vertical axis in fig. 3(d)). The resulting frequency-
normalized map is presented in fig. 3(e), with the same y-axis as in fig.
3(d), which clearly shows the precise fractional synchronization orders
($p/q$, where $p$ is the observed harmonic and $q$ is the externally driven
harmonic). In this representation, the extent of the locking ranges is also
noticeable as vertical lines, corresponding to several MHz even at very high
harmonics. Such robust fractional synchronization orders are considered rare
[12], but not in such thermal free-carrier oscillators. Such broad high-order
synchronization windows have been observed in optomechanical systems, but only
up to the 4${}^{\text{th}}$ harmonic [15]. The greater compliance observed in
these free-carrier thermal oscillators can be attributed to their resemblance
to relaxation oscillators. The reason for this is the rapid relaxation time-
scales exhibited by thermal-free carriers, which produce an oscillator output
similar to a pulse. This results in a broader frequency response when compared
to optomechanics, where the response tends to be more sinusoidal due to the
resonant filtering of the mechanical modes. These nonlinear features of the
thermal-free carrier limit cycle have been extensively discussed by Abrams et
al. [8].

A useful characteristic of high-harmonic synchronization is the possibility to
effectively suppress the oscillator phase-noise (PN) characteristics [33, 34].
The suppression of phase-noise, however, relies on the compliance of the
oscillator to the injected signal. The wide locking ranges observed in figures
3(d,e) suggests that our oscillators are strong candidates to observe such
phase-noise suppression. Indeed, as we measure the phase-noise performance
when the oscillator is entrained at different harmonics, we observe a
consistent suppression of phase-noise at the oscillator’s fundamental tone.
Fig. 3(f) shows the measured PN spectra of the fundamental tone when the
system is free-running (no modulation) and when it is injection-locked up to
the fifth harmonic ($p/q=1/(1,2,...,5)$), demonstrating a clear PN
suppression. The spectral characteristics of the fundamental tone, which
correspond to the PN measurements in fig. 3(f), are depicted in fig. 3(g). The
wide locking range observed for the $p/q=1/5$ excitation in fig. 3(e) suggests
that this excitation should have the best phase-noise suppression, in line
with usual phase-noise modeling [35] and previous studies [35, 15]. The
overall suppression of PN for higher harmonic injection is also expected, as
shown in fig. 3(f).

Such wide locking ranges are advantageous when exploring synchronization
between two or more oscillators, even if they differ slightly due to
fabrication fluctuations. Moreover, the fact that locking can happen at
different harmonics also permits that two oscillators synchronize at different
frequencies, an unusual characteristic that cannot be found in sinusoidal
oscillators [12]. Furthermore, high-harmonic synchronization could be explored
for all-optical radio-frequency division, a fundamental towards frequency
synthesis [15].

##  IV Synchronization of coupled oscillators

Figure 4:  Synchronization of two oscillators. a) Experimental DC optical
transmission of self-oscillating coupled resonators; b) Experimental
wavelength varying spectrogram of coupled resonators; c) Numerical DC optical
transmission of self-oscillating coupled resonators; d) Numerical wavelength
varying spectrogram of coupled resonators. In (a, c), colored regions mark
different regimes of operation. e-g) Experimental time traces at different
modes of operations; h-j) Numerical time traces at different modes of
operation. In (e-j), colors match the regions in (a, c).

The injection locking results reveal a remarkable property of thermal free
carrier oscillators: their extreme compliance to external stimuli. In this
section, this property is explored to demonstrate synchronization between two
oscillators coupled by their evanescent optical field, (see inset in fig.
2(c)). Without employing any external modulation, the pump laser is swept, in
steps of 50 pm, from shorter to longer wavelength around a set of selected
optical resonances (shown in fig. 2(c)). The pump laser power was increased to
approximately 2 mW to ensure that the oscillation threshold is reached, a
consequence of the larger optical mode volume and weaker coupling of each
optical mode to the bus waveguide (tapered optical fiber).

Instead of observing a steady evolution of the oscillator frequency, as in
fig. 3(a), the coupled cavity dynamics exhibits a more complex behavior. Such
behavior, shown in fig. 4(a,b) is very well reproduced by our numerical model
(fig. 4(c,d), see supplementary information). The first striking difference we
observe in the double cavity system spectrogram is the appearance of
qualitatively distinct regions (highlighted in blue,red and green in fig.
4(a,c)). The temporal behavior of the transmitted signal (AC component)
measured by the oscilloscope is shown in fig. 4(e-g); along with their
simulation counterparts in fig. 4(h-j). When the system is just starting to
oscillate (blue region in fig. 4(a,c)) a single low frequency pulse appears,
similar to the single cavity oscillation condition. In this case there is a
single path in both measured and simulated phase diagrams (amplitude versus
amplitude’s time derivative), as shown in fig. 4(e,h)).

As the pump reaches the red region in fig. 4(a,c), a second distinct peak
appears in the time trace (fig. 4(f,i)) and the spectrogram in fig. 4(b) shows
very clear harmonics. Also, the phase diagram presents two distinct paths.
This is the first feature we observe in the time trace that is exclusive of
the double cavity system. Further red-shifting the pump wavelength, in between
the red and green regions in fig.fig. 4(a), the two oscillators lose
synchronization and produce erratic pulses, which manifests as an elevated
noise level in the spectrogram. Finally, once the green region is reached, the
two oscillators reestablish synchronization for a wide range of laser detuning
($\approx 3.8$ nm).

Although the experimentally observed features suggest that the coupled
oscillator dynamics and synchronization are present, it is difficult to grasp
the dynamics of temperature, free-carrier density, and light from the
transmission traces alone. To better understand the dynamics in each of the
indicated regions, we now turn our attention to the numerical simulations of
the coupled system.

###  IV.1 The intricate dynamics of synchronized oscillators

The good qualitative agreement with the theoretical model highlighted in fig.
4 suggests that individual dynamics may be unveiled by exploring the numerical
model. In particular, we extract the instantaneous optical frequency shift
($\Delta\omega_{1,2}(t)$) of each cavity based on their coupling to the time-
dependent free-carrier density ($\text{N}_{1,2}(t)$) and temperatures
variation ($\theta_{1,2}(t)$),

$\displaystyle\Delta\omega_{1}(t)$
$\displaystyle=g_{\text{th}}\theta_{1}(t)+g_{\text{N}}N_{1}(t)$ (1)
$\displaystyle\Delta\omega_{2}(t)$
$\displaystyle=g_{\text{th}}\theta_{2}(t)+g_{\text{N}}N_{2}(t),$ (2)

where $g_{\text{th}}$ and $g_{\text{N}}$ are the thermal and carrier coupling
parameters; $\theta_{1,2}(t)=T_{1,2}(t)-T_{0}$ (see supplementary for
definition of all terms). As detailed in the supplementary material, cavity 1
is the one coupled to the waveguide, while cavity 2 only receives light
through the coupling to its companion.

Figure 5:  Temperature variations and FC density of a self-pulsing double
optical cavity system. a-c) AC temperature variation of each resonator; d-f)
AC free-carrier density variation of each resonator. The AC part of each
contribution is obtained by subtracting the time-average of each
temperature/carrier dynamical variable; g) Individual cavity resonance shifts
around equilibrium ($\overline{\Delta\omega}$); h) Difference between the two
cavities’ resonances. In (g, h) the frequencies have all been normalized by
the inter-cavity coupling rate. In all figures, the colors correspond to the
traces show in fig. 4(h-j)

We investigate these parameters to show that the temporal characteristics
revealed in fig. 4(h-j) are caused by carrier and temperature dynamics that
are not always synchronized in the trivial ($p/q=1/1$) situation. Moreover,
the frequency oscillations ($\Delta\omega_{1,2}(t)$) are comparable to the
evanescent coupling rate $J$ and affect the coupling between the optical
modes. The individual AC components of temperature
($\Delta\theta_{1,2}(t)=\Delta\theta_{1,2}(t)-\overline{\theta}_{1,2}$) and FC
contributions ($\Delta N_{1,2}(t)=N_{1,2}(t)-\overline{N}_{1,2}$, where the
bar indicates average over time) to the resonator frequencies are shown in
fig. 5(a-c) and fig. 5(d-f), with colors corresponding to the transmission
traces highlighted in fig. 4(h-j). Fig. 5(a,d) reveals that the single cycle
behavior observed in fig. 4(e,h) occurs when the dynamics is concentrated only
in cavity 2, as both $\Delta\theta_{1}$ and $\Delta N_{1}$ are significantly
smaller in resonator 1. As a result, the time trace is similar to the single
oscillator case.

In the red region shown in fig. 5(b,e), both oscillators display similar
magnitudes. However, we observed that cavity 2 consistently exhibits peaks in
the carrier density $\Delta N_{2}$, while cavity 1 appears to skip every other
pulse. This indicates synchronization between the two cavities, but not at the
trivial condition of $p/q=1/1$. Rather, cavity 2 oscillates at twice the
frequency of cavity 1 ($p/q=1/2$). We also observed similar behavior for
$\Delta\theta_{2}$, albeit with a smaller contrast between the two cavities.
The RF spectrograms in this region (fig. 4(b)) reflect this temporal behavior
by showing a denser pattern of harmonics. This characteristic is further
highlighted in fig. 5(g), where we examine the behavior of the time-dependent
frequency shifts (eqs. 1,2) in the gray-shaded region of fig. 5(b,e). During
the first pulse, both cavities experience similar resonance frequency shifts,
while during the second pulse, only cavity 2 undergoes a significant shift.

Another important aspect is that the magnitude of the frequency shift is
comparable with the evanescent coupling rate $J$, showing that the very
coupling between the cavities is affected by their dynamical behavior. Indeed,
we show in fig. 5(h) the difference between the cavities’ resonances, showing
that their relative frequency changes as much as
$\Delta\omega_{1}-\Delta\omega_{2}\approx 0.4J$. Such large changes in the
bare resonator frequencies significantly affect the energy stored in each
cavity and, thus, the oscillation dynamics.

Finally, in the green region (fig. 5(c, f)), where the synchronization is very
well defined and the pulsing behavior is regular in fig. 4(g, j), both
temperature (fig. 5(c)) and free-carrier density (fig. 5(f)) oscillate with
the same period. Nevertheless, even though the temperature behavior of both
resonators only shows a small offset, the free-carrier density presents quite
a different behavior. The slightly delayed peak in cavity 1 causes a quick
blue-shift in its optical resonance frequency, while cavity 2 is already
finalizing its cycle and red-shifting. This is what causes the extra knot we
observe in the phase diagrams in fig. 4(j).

##  V Conclusion

In summary, we demonstrated that carrier-thermal oscillations are extremely
compliant to external stimuli, and as such can be synchronized even by weak
signals at a wide range of its harmonics. The low phase-noise performance and
the multi-octave spanning synchronization range could be explored for all-
optical RF frequency division. Furthermore, we explored the extreme compliance
of such oscillator and demonstrated the intricate synchronization dynamics of
a pair of such resonators. The numerical model of coupled resonators allowed a
detailed understanding of the intricate synchronization dynamics. These
synchronized oscillators, albeit nominally identical exhibit synchronization
at distinct harmonics, opening a path for scalable array synchronization of
dissimilar silicon oscillators.

#### Funding

This work was supported by São Paulo Research Foundation (FAPESP) through
grants 19/14377-5, 18/15577-5, 18/15580-6, 18/25339-4, Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES), and CNPq.

#### Acknowledgments

The authors thank A. Von Zuben for technical support. We also acknowledge
CCSNano-UNICAMP for providing part of the microfabrication infrastructure.

#### Disclosures

The authors declare no conflicts of interest.

#### Data Availability

All simulation files and data are available upon request to the authors.

## References

  * Heinrich  _et al._ [2011] G. Heinrich, M. Ludwig, J. Qian, B. Kubala, and F. Marquardt, Collective dynamics in optomechanical arrays, Phys. Rev. Lett. 107, 043603 (2011).
  * Marandi  _et al._ [2014] A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Yamamoto, Network of time-multiplexed optical parametric oscillators as a coherent ising machine, Nature Photonics 8, 937 (2014).
  * Wang  _et al._ [2013] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, Coherent ising machine based on degenerate optical parametric oscillators, Phys. Rev. A 88, 063853 (2013).
  * Zhang  _et al._ [2015] M. Zhang, S. Shah, J. Cardenas, and M. Lipson, Synchronization and phase noise reduction in micromechanical oscillator arrays coupled through light, Phys. Rev. Lett. 115, 163902 (2015).
  * Okawachi  _et al._ [2015] Y. Okawachi, M. Yu, K. Luke, D. O. Carvalho, S. Ramelow, A. Farsi, M. Lipson, and A. L. Gaeta, Dual-pumped degenerate kerr oscillator in a silicon nitride microresonator, Opt. Lett. 40, 5267 (2015).
  * Utsunomiya  _et al._ [2015] S. Utsunomiya, N. Namekata, K. Takata, D. Akamatsu, S. Inoue, and Y. Yamamoto, Binary phase oscillation of two mutually coupled semiconductor lasers, Opt. Express 23, 6029 (2015).
  * Barclay  _et al._ [2005] P. E. Barclay, K. Srinivasan, and O. Painter, Nonlinear response of silicon photonic crystal microresonators excited via an integrated waveguide and fiber taper, Opt. Express 13, 801 (2005).
  * Abrams  _et al._ [2014] D. M. Abrams, A. Slawik, and K. Srinivasan, Nonlinear oscillations and bifurcations in silicon photonic microresonators, Phys. Rev. Lett. 112, 123901 (2014).
  * Navarro-Urrios  _et al._ [2017] D. Navarro-Urrios, N. E. Capuj, M. F. Colombano, P. D. García, M. Sledzinska, F. Alzina, A. Griol, A. Martínez, and C. M. Sotomayor-Torres, Nonlinear dynamics and chaos in an optomechanical beam, Nature Communications 8, 14965 (2017).
  * Arregui  _et al._ [2021] G. Arregui, M. F. Colombano, J. Maire, A. Pitanti, N. E. Capuj, A. Griol, A. Martínez, C. M. Sotomayor-Torres, and D. Navarro-Urrios, Injection locking in an optomechanical coherent phonon source, Nanophotonics 10, 1319 (2021).
  * Borghi  _et al._ [2021] M. Borghi, D. Bazzanella, M. Mancinelli, and L. Pavesi, On the modeling of thermal and free carrier nonlinearities in silicon-on-insulator microring resonators, Opt. Express 29, 4363 (2021).
  * Pikovsky  _et al._ [2001] A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A Universal Concept in Nonlinear Sciences_, Cambridge nonlinear science series (Cambridge University Press, 2001).
  * Hamerly and Mabuchi [2015] R. Hamerly and H. Mabuchi, Optical devices based on limit cycles and amplification in semiconductor optical cavities, Phys. Rev. Appl. 4, 024016 (2015).
  * Zhang  _et al._ [2012] M. Zhang, G. S. Wiederhecker, S. Manipatruni, A. Barnard, P. McEuen, and M. Lipson, Synchronization of micromechanical oscillators using light, Phys. Rev. Lett. 109, 233906 (2012).
  * Rodrigues  _et al._ [2021] C. C. Rodrigues, C. M. Kersul, A. G. Primo, M. Lipson, T. P. M. Alegre, and G. S. Wiederhecker, Optomechanical synchronization across multi-octave frequency spans, Nature Communications 12, 5625 (2021).
  * Matheny  _et al._ [2019] M. H. Matheny, J. Emenheiser, W. Fon, A. Chapman, A. Salova, M. Rohden, J. Li, M. H. de Badyn, M. Pósfai, L. Duenas-Osorio, M. Mesbahi, J. P. Crutchfield, M. C. Cross, R. M. D’Souza, and M. L. Roukes, Exotic states in a simple network of nanoelectromechanical oscillators, Science 363, eaav7932 (2019), https://www.science.org/doi/pdf/10.1126/science.aav7932 .
  * Bagheri  _et al._ [2013] M. Bagheri, M. Poot, L. Fan, F. Marquardt, and H. X. Tang, Photonic cavity synchronization of nanomechanical oscillators, Phys. Rev. Lett. 111, 213902 (2013).
  * Johnson  _et al._ [2006] T. J. Johnson, M. Borselli, and O. Painter, Self-induced optical modulation of the transmission through a high-q silicon microdisk resonator, Opt. Express 14, 817 (2006).
  * Soltani  _et al._ [2012] M. Soltani, S. Yegnanarayanan, Q. Li, A. A. Eftekhar, and A. Adibi, Self-sustained gigahertz electronic oscillations in ultrahigh-$q$ photonic microresonators, Phys. Rev. A 85, 053819 (2012).
  * Navarro-Urrios  _et al._ [2015] D. Navarro-Urrios, N. E. Capuj, J. Gomis-Bresco, F. Alzina, A. Pitanti, A. Griol, A. Martínez, and C. M. Sotomayor Torres, A self-stabilized coherent phonon source driven by optical forces, Scientific Reports 5, 15733 (2015).
  * Green and Keevers [1995] M. A. Green and M. J. Keevers, Optical properties of intrinsic silicon at 300 k, Progress in Photovoltaics: Research and Applications 3, 189 (1995), https://onlinelibrary.wiley.com/doi/pdf/10.1002/pip.4670030303 .
  * Dinu  _et al._ [2003] M. Dinu, F. Quochi, and H. Garcia, Third-order nonlinearities in silicon at telecom wavelengths, Applied Physics Letters 82, 2954 (2003), https://doi.org/10.1063/1.1571665 .
  * Tsang  _et al._ [2002] H. K. Tsang, C. S. Wong, T. K. Liang, I. E. Day, S. W. Roberts, A. Harpin, J. Drake, and M. Asghari, Optical dispersion, two-photon absorption and self-phase modulation in silicon waveguides at 1.5 $\mu$m wavelength, Applied Physics Letters 80, 416 (2002), https://doi.org/10.1063/1.1435801 .
  * Nedeljkovic  _et al._ [2011] M. Nedeljkovic, R. Soref, and G. Z. Mashanovich, Free-carrier electrorefraction and electroabsorption modulation predictions for silicon over the 1–14- $\mu\hbox{m}$ infrared wavelength range, IEEE Photonics Journal 3, 1171 (2011).
  * Lautenschlager  _et al._ [1987] P. Lautenschlager, M. Garriga, L. Vina, and M. Cardona, Temperature dependence of the dielectric function and interband critical points in silicon, Phys. Rev. B 36, 4821 (1987).
  * Turner-Foster  _et al._ [2010] A. C. Turner-Foster, M. A. Foster, J. S. Levy, C. B. Poitras, R. Salem, A. L. Gaeta, and M. Lipson, Ultrashort free-carrier lifetime in low-loss silicon nanowaveguides, Opt. Express 18, 3582 (2010).
  * Aldaya  _et al._ [2016] I. Aldaya, A. Gil-Molina, H. L. Fragnito, and P. Dainese, Simple method to simultaneous characterize two-photon absorption coefficient and carrier lifetime in silicon waveguides, in  _Latin America Optics and Photonics Conference_ (Optica Publishing Group, 2016) p. LTu5D.4.
  * Dimitropoulos  _et al._ [2005] D. Dimitropoulos, R. Jhaveri, R. Claps, J. C. S. Woo, and B. Jalali, Lifetime of photogenerated carriers in silicon-on-insulator rib waveguides, Applied Physics Letters 86, 071115 (2005), https://aip.scitation.org/doi/pdf/10.1063/1.1866635 .
  * Van Vaerenbergh  _et al._ [2012] T. Van Vaerenbergh, M. Fiers, J. Dambre, and P. Bienstman, Simplified description of self-pulsation and excitability by thermal and free-carrier effects in semiconductor microcavities, Phys. Rev. A 86, 063808 (2012).
  * Borselli  _et al._ [2007] M. Borselli, T. J. Johnson, and O. Painter, Accurate measurement of scattering and absorption loss in microphotonic devices, Opt. Lett. 32, 2954 (2007).
  * Knight  _et al._ [1997] J. C. J. C. Knight, G. Cheung, F. Jacques, and T. A. T. A. Birks, Phase-matched excitation of whispering-gallery-mode resonances by a fiber taper, Optics Letters 22, 1129 (1997).
  * Wang  _et al._ [2016] H. Wang, Y. Dhayalan, and E. Buks, Devil’s staircase in an optomechanical cavity, Phys. Rev. E 93, 023007 (2016).
  * Leijssen  _et al._ [2017] R. Leijssen, G. R. La Gala, L. Freisem, J. T. Muhonen, and E. Verhagen, Nonlinear cavity optomechanics with nanomechanical thermal fluctuations, Nature Communications 8, ncomms16024 (2017).
  * Kennedy  _et al._ [2011] M. P. Kennedy, H. Mo, and X. Dong, Experimental characterization of arnold tongues in injection-locked cmos lc frequency dividers with tail and direct injection, in  _2011 20th European Conference on Circuit Theory and Design (ECCTD)_ (2011) pp. 484–487.
  * Plessas [2011] F. Plessas, A study of superharmonic injection locking in multiband frequency dividers, International Journal of Circuit Theory and Applications 39, 397 (2011), https://onlinelibrary.wiley.com/doi/pdf/10.1002/cta.644 .
  * Slack [1964] G. A. Slack, Thermal conductivity of pure and impure silicon, silicon carbide, and diamond, Journal of Applied Physics 35, 3460 (1964), https://doi.org/10.1063/1.1713251 .
  * Glassbrenner and Slack [1964] C. J. Glassbrenner and G. A. Slack, Thermal conductivity of silicon and germanium from 3°k to the melting point, Phys. Rev. 134, A1058 (1964).
  * Asheghi  _et al._ [1998] M. Asheghi, M. N. Touzelbaev, K. E. Goodson, Y. K. Leung, and S. S. Wong, Temperature-Dependent Thermal Conductivity of Single-Crystal Silicon Layers in SOI Substrates, Journal of Heat Transfer 120, 30 (1998).

##  Appendix A Supplementary Information

###  A.1 Theoretical model

For a single cavity, we model this phenomenon with three coupled equations,
one for the optical field (eq. 3), one for the density of charge carriers (eq.
4) and one for the temperature variation (eq. 5) [29, 8]:

$\displaystyle\frac{da(t)}{dt}$
$\displaystyle=i\left[\Delta+\text{g}_{\theta}\theta(t)+\text{g}_{N}N(t)\right]a(t)-\frac{\kappa+\alpha_{TPA}\left|~{}a(t)~{}\right|^{2}+\alpha_{N}\left|N(t)\right|}{2}a(t)+\sqrt{\kappa_{e}P_{in}}$
(3) $\displaystyle\frac{dN(t)}{dt}$
$\displaystyle=-\gamma_{FC}N(t)+\beta_{FC}\left|~{}a(t)~{}\right|^{4}$ (4)
$\displaystyle\frac{d\theta(t)}{dt}$
$\displaystyle=-\gamma_{th}\theta(t)+\beta_{th}\left(\kappa_{lin}+\sigma_{Si}v_{g}N(t)+\alpha_{TPA}\left|~{}a(t)~{}\right|^{2}\right)\left|~{}a(t)~{}\right|^{2}$
(5)

The definition of the parameters of these equations are summarized on table
LABEL:tab:5.1.

Table A.1: Parameters of the self-pulsing equations. Parameter |  Description | Parameter |  Description
---|---|---|---
$\left|~{}a~{}\right|^{2}$ |  Intra-cavity stored energy | $\gamma_{FC}$ |  Free-carriers decay rate
$N$ |  Density of carriers | $\gamma_{th}$ |  Thermal decay rate
$\theta=T-T_{0}$ |  Temperature variation with respect to equilibrium | $\text{g}_{\theta}$ |  Coupling coefficient of optical resonance to temperature variation
$\Delta=\omega_{l}-\omega_{0}$ |  Pump laser frequency detuning with respect to the optical resonance of the cavity | $\text{g}_{N}$ |  Coupling coefficient of optical resonance to the density of carriers
$\kappa$ |  Total linear decay rate of the optical cavity (includes linear absorption, scattering and coupling to the waveguide) | $\alpha_{TPA}$ |  Two-photon absorption loss parameter
$\kappa_{e}$ |  Coupling rate between optical cavity and waveguide | $\alpha_{N}$ |  Free-carrier optical loss parameter
$\kappa_{lin}$ |  Optical decay rate due to linear absorption | $\beta_{FC}$ |  Free-carrier generation parameter
$v_{g}$ |  Optical mode group velocity | $\beta_{th}$ |  Thermal source parameter
$\sigma_{Si}$ |  Silicon free-carrier absorption cross-section | $P_{in}$ |  Pump laser power

The coupling coefficients of the optical resonance to temperature and carrier
density are given by

$\text{g}_{\theta}=\frac{\omega_{0}}{n_{g}}\frac{dn}{dT}$ (6)

and

$\text{g}_{N}=\frac{\omega_{0}}{n_{g}}\frac{dn}{dN},$ (7)

where $n$ is the refractive index, $n_{g}$ is the group index of refraction,
$\omega_{0}$ is the optical resonance angular frequency in the absence of non-
linear effects (cold-cavity) and $T$ is the temperature. Note that the
opposite effect of $\theta$ and $N$ is implicit in this equations through
$dn/dT$, which is positive, and $dn/dN$, which is negative.

The TPA and FC absorption parameters are defined as

$\alpha_{TPA}=\frac{\Gamma_{TPA}\beta_{Si}c^{2}}{V_{TPA}n_{g}^{2}}$ (8)

and

$\alpha_{N}=\frac{\sigma_{Si}c}{n_{g}},$ (9)

where $\beta_{Si}$ is the TPA constant for silicon and $c$ is the speed of
light.

The FC and thermal source parameters are defined as

$\beta_{FC}=\frac{\Gamma_{FC}\beta_{Si}c^{2}}{2\hbar\omega_{0}n_{g}^{2}V_{FC}^{2}}$
(10)

and

$\beta_{th}=\frac{\Gamma_{disk}}{\rho c_{p}V_{disk}},$ (11)

where $\hbar$ is the reduced Planck constant and $c_{p}$ is the material heat
capacity. The parameters $\Gamma_{TPA}$ and $\Gamma_{FC}$ are overlap factors
for TPA and FC absorption, respectively, while $V_{TPA}$, $V_{FC}$ and
$V_{disk}$ are effective volumes available for TPA and FC absorption and the
total disk volume, respectively. These parameters values are determined from
finite element method (FEM) numerical simulations performed using COMSOL® and
are defined as [7]

$\displaystyle\Gamma_{TPA}$
$\displaystyle=\frac{\int_{Si}n^{4}\left(\vec{r}\right)E^{4}\left(\vec{r}\right)d\vec{r}}{\int
n^{4}\left(\vec{r}\right)E^{4}\left(\vec{r}\right)d\vec{r}}\text{,}$ (12)
$\displaystyle\Gamma_{FC}$
$\displaystyle=\frac{\int_{Si}n^{6}\left(\vec{r}\right)E^{6}\left(\vec{r}\right)d\vec{r}}{\int
n^{6}\left(\vec{r}\right)E^{6}\left(\vec{r}\right)d\vec{r}}\text{,}$ (13)
$\displaystyle V_{TPA}$ $\displaystyle=\frac{\left(\int
n^{2}\left(\vec{r}\right)E^{2}\left(\vec{r}\right)d\vec{r}\right)^{2}}{\int
n^{4}\left(\vec{r}\right)E^{4}\left(\vec{r}\right)d\vec{r}}$ (14)

and

$V_{FC}^{2}=\frac{\left(\int
n^{2}\left(\vec{r}\right)E^{2}\left(\vec{r}\right)d\vec{r}\right)^{3}}{\int
n^{6}\left(\vec{r}\right)E^{6}\left(\vec{r}\right)d\vec{r}},$ (15)

where $E(\vec{r})$ is the spacial distributions of the electric field
amplitude.

There are 3 distinct time scales present in equations 3 to 5. The fastest
time-scale is the optical life-time ($\tau_{opt}=1/\kappa$), which is in the
order of tens of picoseconds for cavities with optical quality factors on the
order of 80k. Then comes the FC life-time ($\tau_{FC}=1/\gamma_{FC}$), in the
order of a few nanoseconds. And the slowest time-scale is the thermal life-
time ($\tau_{th}=1/\gamma_{th}$), in the microsecond scale. One can then say
that the optical field inside the cavity responds instantaneously to any FC or
thermal variation. This means that the stored optical energy
($\left|~{}a~{}\right|^{2}$) is always in steady-state at any time, following
any variations of $N$ and $\theta$ adiabatically. We can then write $da/dt=0$
and solve the first equation for $a(t)=a(N(t),\theta(t))$.

To do so it is necessary first to take care of a non-linear absorption term
that depends explicitly on $a(t)$ in equation 3, namely the optical loss due
to TPA. In fact, by solving the steady-state problem for the carriers density
($dN(t)/dt=0$) it is possible to evaluate the contribution of each loss term
in equation 3 at any given static amount of energy inside the cavity (fig.
A.1). For a typical input power of 1 mW to 2 mW, the stored energy inside of
our cavities is on the order of 10 fJ. For this amount of energy both non-
linear absorption coefficients are still smaller than the linear one
($\kappa$), but the terms due to TPA is much smaller than the other two,
unless for very little amounts of stored energy when it is comparable to the
FC absorption term.

Figure A.1: Optical losses due to FC absorption and TPA as a function of the
intra-cavity energy.

Nevertheless, because the TPA term is always much smaller than the total
optical decay rate, this term can be ignored and the steady-state solution of
equation 3 can be simplified, resulting in a single and simple solution for
the intra-cavity field amplitude, given by

$a(t)=\frac{2\sqrt{\kappa_{e}P_{in}}}{2i\left(\Delta++\text{g}_{\theta}\theta(t)+\text{g}_{N}N(t)\right)-\kappa+\alpha_{N}\left|N(t)\right|}.$
(16)

Inserting equation 16 into equations 4 and 5 we obtain two equations of real
arguments for the temperature variation and FC density, given by

$\displaystyle\frac{dN(t)}{dt}$
$\displaystyle=-\gamma_{FC}N(t)+\beta_{FC}U(t)^{2}$ (17)
$\displaystyle\frac{d\theta(t)}{dt}$
$\displaystyle=-\gamma_{th}\theta(t)+\beta_{th}\left(\kappa_{lin}+\sigma_{Si}v_{g}N(t)+\alpha_{TPA}U(t)\right)U(t),$
(18)

where $U(t)=\left|~{}a(t)~{}\right|^{2}$, for simplicity. This approximation
will be of special value when solving the problem of 2 cavities, as it will
lead to a smaller memory consumption and shorter time to solve the problems.
It is important to note that, although the TPA term was removed from the
optical equation (eq. 3), this term is maintained in equation 5, minimizing
the amount of approximations to solve the problem.

The material parameters used to properly solve the problem are obtained from
the literature [7, 36, 37, 38, 30, 27, 18]. The thermal life-time was
estimated by considering the conduction of heat from the top silicon layer
through a truncated cone oxide pedestal, with top diameter given by the
undercut of the disk obtained from optical images and bottom diameter was made
equal to that of the disk. This produced a typical time-scale of 7 µs for the
single cavity and 3.5 µs for the double cavity, with the difference given by
the different undercut size in each sample. The FC life-time for both systems
was assumed to be 4 times the value presented in reference [27], which results
in a timescale of 5.8 ns. We chose reference [27] because their measurements
were performed in waveguides fabricated by the same foundry as the samples
used in this work. The scale factor of the FC lifetime was chosen so that the
time dependent solutions had similar frequencies to the experimental data. The
linear absorption rate ($\kappa_{lin}$) was set to half of the intrinsic
linear decay rate ($\kappa_{i}$), based on the results presented in reference
[30].

The material and other parameters used in the model are listed in table A.2.
Here $t_{\text{Si}}$ is the silicon layer thickness (220 nm) and $r_{disk}$ is
the disk radius (2 µm). The two values given for $\tau_{th}$ are for the
double and single cavity systems, respectively.

Table A.2: Table with parameter values used in the numerical simulations. Parameter | Value
---|---
$\rho$ | 2330 kg/m3
$c_{p}$ | 712 J/(kg K)
$n_{Si}$ | 3.485
$n_{g}$ | 0.99 $n_{Si}$
$\lambda_{0}$ | 1471.63 nm
$\kappa$ |  $2\pi\times 2.55$ GHz
$\kappa_{e}$ |  $2\pi\times 1.17$ GHz
$\kappa_{i}$ | $\kappa-\kappa_{e}$
$\kappa_{lin}$ | $\kappa_{i}/2$
$\tau_{FC}$ |  $5.8$ ns
$\gamma_{FC}$ | $1/\tau_{FC}$
$\tau_{th}$ | (7.0, 3.5) µs
$\gamma_{th}$ | $1/\tau_{th}$
$T_{0}$ | 300 K
$dn/d\theta$ |  $1.86\times 10^{-4}$ K-1
$dn/dN$ |  $-1.73\times 10^{-27}$ m3
$\sigma_{Si}$ | $10^{-21}\text{m}^{2}$
$\beta_{Si}$ |  $8.410^{-12}$ m/W
$\Gamma_{TPA}$ | 1
$V_{TPA}$ | $1.24\times 10^{-18}\text{m}^{3}$
$\Gamma_{FC}$ | 1
$V_{FC}$ | $1.15\times 10^{-18}\text{m}^{3}$
$\Gamma_{disk}$ | 0.97
$V_{disk}$ | $\pi r_{disk}^{2}t_{Si}=2.76\times 10^{-18}\text{m}^{3}$

###  A.2 Numerical solution details

Because the optical resonance shifts when the optical energy inside the cavity
increases, it is necessary to solve equations 17 and 18 with a time dependent
detuning ($\Delta=\Delta(t)$), as the example in figure A.2. In this function
$\Delta(t)$ varies from the blue side of the cold cavity ($\Delta/\kappa>0$),
until a given stop point on the red side of the cold resonance
($\Delta/\kappa<0$). If this ramp is done too slowly, the solutions tend to
show only static temperature and FC density. That is because the model doesn’t
take noise into account, hence there is nothing to perturb the solution out of
its equilibrium.

On the other hand, if the ramp is done fast enough such that oscillations
appear in the solutions, it typically causes artifacts on the solutions due to
the variations of the stored optical energy being on the same time-scale of
the phenomena of interest. To overcome this problem, the system of equations
is solved multiple times, changing the stop detuning in each solution and
storing only the data related to the times when $\Delta$ is static (red line
in fig. A.2). In this way the ramp of detuning can be done very quickly
without affecting the final solution.

Figure A.2: Time dependent detuning function. The black region simulates the
laser frequency sweep, taking the optical resonance non-linear shift into
account. The red region is where time dependent data is stored and the Fourier
transform evaluated. The horizontal axis is time normalized by the thermal
life-time. The vertical axis is the laser detuning with respect to the cold
cavity, normalized by the total cold-cavity linewidth.

An important observation is that the $\Delta$ values presented in figure A.2
are relative to the cold cavity, i.e., without considering the non-linear
resonance shift. As will be clear later in this section, the effective
detuning is always positive with respect to the shifted resonance position.
Also, we note that these numerical solutions presented here were obtained
before proper characterization of the FC and thermal lifetimes of the samples,
what invariably leads to deviations in quantitative results between model and
experiment. Nevertheless, good agreement in the general behavior of the
systems is achieved for both single and double cavity cases.

###  A.3 Single cavity

To reproduce the single cavity experiments the optical modes are modeled with
a splitting due to coupling of counter propagating modes. The expected optical
transmission spectrum for a low power pump is shown in figure A.3, with an
optical decay rate $\kappa=2\pi\times$2.5 GHz and a coupling rate between
modes $\text{g}_{cp}=2\pi\times$13 GHz. The cavity-waveguide coupling rate
($\kappa_{e}$) is 1.1 GHz and the central wavelength was set to approximately
1471 nm. The list of values for all parameters used in this case are listed in
table A.2.

Figure A.3: Numerical optical transmission spectrum of a single cavity with
coupling of counter-propagating modes. The coupling rate between modes is 13
GHz and the total optical decay rate is 2.5 GHz.

The equations are solved in Mathematica® and, for a given final detuning
value, result in time-traces like the ones shown in figure A.4, which show the
traces over one complete cycle of oscillation. By substituting the $\theta(t)$
and $N(t)$ solutions into equation 16 and calculating the normalized time
dependent transmission, the curve presented in figure A.4a is obtained
(analogous to Figure 1(b,c) in the main text). Additionally, Figure A.4b shows
the frequency shift of the resonances relative to the equilibrium position,
which is given by

$\Delta\omega=\text{g}_{\theta}\theta(t)+\text{g}_{N}N(t)-\overline{\Delta\omega},$
(19)

where
$\overline{\Delta\omega}=\text{g}_{\theta}\overline{\theta}+\text{g}_{N}\overline{N}$
is the shift caused by the mean variation of both temperature and FC density,
with the overbar indicating time average. Figures A.4c and A.4d show the time
behavior of FC and temperature variation over a single period. The colors
indicate the regions where a relative blue-shift (blue) or red-shift (red)
occurs, respectively.

Figure A.4: Numerical time-traces of a self-pulsing single optical cavity over
one period. a) Optical transmission normalized by the input power $P_{in}$. b)
Total resonance shift with respect to the mean resonance shift normalized to
$\text{g}_{cp}$. c) Free-carriers density. d) Temperature variation with
respect to the cold-cavity equilibrium value. The colors indicate the regions
where there is relative red and blue-shift of the resonance. Figure A.5:
Numerical time-dependent optical transmission spectrogram. The vertical line
indicates the position of the pump laser frequency.

Observe that the resonance path reproduces exactly the trace presented in
figure A.4b, and now it is possible to see the two moments when the bluer
resonance passes by the laser position, causing the minima in the pulse
pattern. Also, observe that because the two resonances are of the same cavity,
experiencing then the same shifts due to the FC and temperature variations,
both resonances move exactly in phase, never changing their relative distance.

Finally, taking the Fourier transform of the optical transmission time-trace
(fig. A.4a) for different static detuning positions produces a spectrogram,
with the time axis substituted by a pump wavelength axis (fig. A.6). The pump
wavelength dependent transmission (fig. A.6a) is obtained by taking the
average of the oscillating transmission signal for each detuning. The pulses
spectra present various harmonics, spanning up to the GHz scale, which is
expected for the square-like pulses produced in this process. Note that the
detuning for which oscillations start is very well marked by both an abrupt
increase in the average transmitted signal and the appearance of the peaks in
the spectrum.

Figure A.6: Numerical wavelength dependent spectrogram of a self-pulsing
single optical cavity. a) Normalized transmission of a detuning sweeping pump
laser, with the typical triangular shape due to non-linear bi-stability [7].
The blue region is where the self-pulsing occurs. b) Detuning dependent
spectrogram of the transmitted self-pulsing light. The vertical scale follows
the same scale in (a).

###  A.4 Coupled cavities

To obtain the equations for coupled cavities or for the case when a single
cavity presents coupled modes, one must turn to the coupled mode theory of
optical modes. This method uses the same initial equation, but now one for
each mode coupled to each other, which in the rotating frame are given by

$\displaystyle\frac{da_{1}(t)}{dt}$ $\displaystyle=-\mathfrak{i}\Delta_{1}\
a_{1}(t)-\frac{\kappa_{1}}{2}\ a_{1}(t)+\sqrt{\kappa_{e}}\
s_{in}s+\mathfrak{i}\frac{\text{g}_{12}}{2}a_{2}(t)\hskip
14.22636pt\text{and}$ (20) $\displaystyle\frac{da_{2}(t)}{dt}$
$\displaystyle=-\mathfrak{i}\Delta_{2}\ a_{2}(t)-\frac{\kappa_{2}}{2}\
a_{2}(t)+\mathfrak{i}\frac{\text{g}_{12}}{2}a_{1}(t),$ (21)

where $\text{g}_{12}$ is the coupling rate between modes 1 and 2, which can be
either the modes of two cavities or two modes of the same cavity. Rigorously,
this coupling rate can be estimated from the overlap of the modes and the
energy and momentum conservation laws. Here we wrote the more general problem,
with the two modes having different resonant frequencies, resulting in
different detunings of the laser, and different intrinsic decay rates. Note
that for a typical coupled cavity system $\kappa_{2}$ doesn’t include a
contribution from $\kappa_{e}$, as the second cavity only receives light
through the first, and only the latter is coupled to the waveguide. This is
not necessarily true for the case of coupling between two modes of the same
cavity, e.g., when there is coupling of counter-propagating modes. Also, if
there is coupling between different modes of the same cavity, because their
coupling to the waveguide can be different, this must be considered and the
equations modified. Nonetheless, we do not add such complications here because
there is little contribution for the intended discussion.

As shown in figure A.7a, in this case one of the cavities is coupled directly
to the waveguide (which will be called $\text{C}_{1}$), while the other
($\text{C}_{2}$) only receives light that coupled into $\text{C}_{1}$. Figure
A.7b shows the optical transmission spectrum for this system, as it would
appear in a low input power measurement. The coupling rate between cavities is
$\text{g}_{12}=2\pi\times$28 GHz, and the resonances are centered around 1471
nm with total linewidths $\kappa=2\pi\times$2.5 GHz.

Figure A.7: Schematics and transmission spectrum of a coupled cavity system
model. a) Schematics of a coupled cavity system. The coupling between cavity
modes is given by $\text{g}_{12}$. b) Optical transmission spectrum of a pair
of coupled resonances. The coupling rate between modes is 28 GHz and the
loaded optical decay rate is 2.5 GHz. The difference in extinction ratio is
due to a mismatch in resonance frequency of 1 GHz between the two cavities.

Note that there is a slight asymmetry in the peaks extinction ratio; that is
because a slight mismatch between the resonances of $\text{C}_{1}$ and
$\text{C}_{2}$ is taken into account. This mismatch was created by making the
$\text{C}_{1}$ with resonance $\omega_{0}$ and $\text{C}_{2}$ with resonance
$\omega_{0}-\delta_{0}$, where $\delta_{0}$ is a positive constant, making the
resonance of $\text{C}_{2}$ redder than that of $\text{C}_{1}$. In this case
the difference is $\delta_{0}\approx 2\pi\times$1 GHz. Besides this resonance
frequency mismatch, the cavities are considered to be identical in terms of
optical losses and material properties, as well as FC and thermal life-time.

Also, in the single cavity case it was demonstrated that the splitting due to
coupling of counter-propagating modes doesn’t add sensible information to the
problem. At most it would cause split dips in the transmission time-trace, and
only if the shift caused by temperature and FC is large enough for both peaks
to pass by the pump laser frequency. Hence we didn’t take this feature into
account for the double cavity problem, even though the real device presents
two kinds of splitting, a larger one for the coupling between cavities and two
smaller ones for the coupling of counter-propagating modes in each cavity.

It is important to note that the two peaks in the spectrum of figure A.7b are
not those of the individual resonances of $\text{C}_{1}$ and $\text{C}_{2}$.
Instead, these are what are commonly named super-modes of the coupled system,
which are comprised of symmetric and anti-symmetric combinations of individual
cavity modes, and their distance is related to the coupling rate
$\text{g}_{12}$. Throughout this chapter we will refer to these resonances by
coupled resonances, in this manner distinguishing them from the individual
cavities’ resonances.

In the double cavity case, the importance of reducing the number of equations
for each cavity becomes apparent. The whole problem involves three equations
per cavity, and the equation for the optical field amplitude ($a(t)$) is
complex, which in practice increases the number of equations to four per
cavity, because any numerical solver will make $a(t)=A(t)+\mathfrak{i}B(t)$
before starting to solve the system of equations. Removing the equation for
$a(t)$ by the adiabatic approximation (the field responds instantly to the
temperature and FC variations) reduces the problem to a set of two real
equations per cavity, making the problem much simpler and quicker to solve.
The system of equations for two cavities is given by

$\displaystyle\frac{dN_{1}(t)}{dt}$
$\displaystyle=-\gamma_{FC}N_{1}(t)+\beta_{FC}U_{1}(t)^{2}$ (22)
$\displaystyle\frac{d\theta_{1}(t)}{dt}$
$\displaystyle=-\gamma_{th}\theta_{1}(t)+\beta_{th}\left(\kappa_{lin}+\sigma_{Si}v_{g}N_{1}(t)+\alpha_{TPA}U_{1}(t)^{2}\right)U_{1}(t)^{2}$
(23) $\displaystyle\frac{dN_{2}(t)}{dt}$
$\displaystyle=-\gamma_{FC}N_{2}(t)+\beta_{FC}U_{2}(t)^{2}$ (24)
$\displaystyle\frac{d\theta_{2}(t)}{dt}$
$\displaystyle=-\gamma_{th}\theta_{2}(t)+\beta_{th}\left(\kappa_{lin}+\sigma_{Si}v_{g}N_{2}(t)+\alpha_{TPA}U_{2}(t)^{2}\right)U_{2}(t)^{2},$
(25)

where $N_{i}$, $\theta_{i}$ and $U_{i}$ are the FC density, temperature
variation and intra-cavity energy of cavities $\text{C}_{1}$ and
$\text{C}_{2}$. Notice that there is no direct coupling between $N$ and
$\theta$ of the two cavities, instead it is mediated exclusively by the
optical field of the cavities, i.e. implicitly through $U_{i}$, which are
given by

$\displaystyle U_{1}(t)$
$\displaystyle=\frac{4\left(4\Delta_{2}^{2}+\kappa_{2}^{2}\right)\kappa_{e}P_{in}}{\left(4\Delta_{1}^{2}+\kappa_{1}^{2}\right)\left(4\Delta_{2}^{2}+\kappa_{2}^{2}\right)+\text{g}_{12}^{2}\left(2\kappa_{1}\kappa_{2}-8\Delta_{1}\Delta_{2}\right)+\text{g}_{12}^{4}}$
(26) $\displaystyle U_{2}(t)$
$\displaystyle=\frac{4\text{g}_{12}^{2}\kappa_{e}P_{in}}{\left(4\Delta_{1}^{2}+\kappa_{1}^{2}\right)\left(4\Delta_{2}^{2}+\kappa_{2}^{2}\right)+\text{g}_{12}^{2}\left(2\kappa_{1}\kappa_{2}-8\Delta_{1}\Delta_{2}\right)+\text{g}_{12}^{4}},$
(27)

where $\Delta_{1}=\Delta_{0}+\text{g}_{\theta}\theta_{1}+\text{g}_{N}N_{1}$
and
$\Delta_{2}=\Delta_{0}+\text{g}_{\theta}\theta_{2}+\text{g}_{N}N_{2}+\delta_{0}$
are the effective detunings due to the variation of $\theta$ and $N$ in each
cavity, with $\Delta_{0}$ being the detuning with respect to the unperturbed
resonance of $\text{C}_{1}$; and $\kappa_{1,2}=\kappa+\alpha_{2}N_{1,2}$ are
the decay rates of each cavity, contemplating both linear and non-linear
contributions.