# Synchronization of spin-driven limit cycle oscillators optically levitated
in vacuum

Oto Brzobohatý<EMAIL_ADDRESS>Martin Duchaň  Petr Jákl  Jan Ježek  Pavel
Zemánek  Stephen H. Simpson<EMAIL_ADDRESS>The Czech Academy of Sciences,
Institute of Scientific Instruments, Královopolská 147, 612 64 Brno, Czech
Republic

###### Abstract

We explore, experimentally and theoretically, the emergence of coherent
coupled oscillations and synchronization between a pair of non-Hermitian,
stochastic, opto-mechanical oscillators, levitated in vacuum. Each oscillator
consists of a polystyrene microsphere trapped in a circularly polarized,
counter-propagating Gaussian laser beam. Non-conservative, azimuthal forces,
deriving from inhomogeneous optical spin, push the micro-particles out of
thermodynamic equilibrium. For modest optical powers each particle shows a
tendency towards orbital circulation. Initially, their stochastic motion is
weakly correlated. As the power is increased, the tendency towards orbital
circulation strengthens and the motion of the particles becomes highly
correlated. Eventually, centripetal forces overcome optical gradient forces
and the oscillators undergo a collective Hopf bifurcation. For laser powers
exceeding this threshold, a pair of limit cycles appear, which synchronize due
to weak optical and hydrodynamic interactions. In principle, arrays of such
Non-Hermitian elements can be arranged, paving the way for opto-mechanical
topological materials or, possibly, classical time crystals. In addition, the
preparation of synchronized states in levitated optomechanics could lead to
new and robust sensors or alternative routes to the entanglement of
macroscopic objects.

##  I Introduction

Over the preceding decade, levitational optomechanics has emerged as a
versatile platform for addressing crucial questions in the physical sciences,
ranging from the macroscopic limits of quantum mechanics Millen and Stickler
(2020) to the thermodynamic limits of computation Konopik  _et al._ (2020). It
makes use of optical forces, which are generated when light scatters from
small particles. These forces can confine or suspend isolated particles in
vacuum, or induce structured interactions, known as optical binding forces,
amongst collections of them Dholakia and Zemánek (2010). Supplementing optical
with electrostatic forces Tebbenjohanns  _et al._ (2019), and combining with
an optical cavity Chang  _et al._ (2010) results in a reconfigurable
experimental system with widely tunable reactive and dissipative forces,
capable of supporting dynamical effects across multiple physical regimes
Gonzalez-Ballestero  _et al._ (2021, 2021).

Most significantly, these techniques have recently enabled motional cooling of
nanoparticles towards and into their quantum mechanical ground state, in
single and multiple degrees of freedom Delić  _et al._ (2020); Piotrowski  _et
al._ (2022); Arita  _et al._ (2022); Liu  _et al._ (2022) with the future
promise of macroscopic entanglement Chauhan  _et al._ (2020); Rudolph  _et
al._ (2020, 2022).

In the classical domain, the harmonic potentials associated with optical
tweezers can confine cooled particles forming high Q oscillators with
exquisite force sensitivity Monteiro  _et al._ (2020); Hempston  _et al._
(2017). However, optomechanical systems can exhibit far richer behaviour.
Optical forces, including interaction forces, are, in general, non-
conservative Sukhov  _et al._ (2015); Roichman  _et al._ (2008); Simpson and
Hanna (2010); Sukhov and Dogariu (2017), and can be holographically sculpted
Woerdemann  _et al._ (2013); Zupancic  _et al._ (2016). In combination,
arbitrarily structured non-conservative and non-linear forces, thermal
fluctuations and dissipation provide the necessary ingredients for numerous
stochastic and dynamic phenomena which are not only of intrinsic interest, but
could also have novel applications in sensing and metrology. Examples include
autonomous stochastic resonance Gang  _et al._ (1993), coherence resonance
Pikovsky and Kurths (1997), stochastic bifurcations Tanabe and Pakdaman (2001)
and stochastic synchronization, all of which are exploited by nature in the
sense apparatus of animals Moss  _et al._ (2004); Laing and Lord (2009).

In this article we demonstrate an archetypal non-equilibrium effect:
synchronization of a pair of noisy limit cycle oscillators Acebrón  _et al._
(2005); Gupta  _et al._ (2018). The oscillators comprise polystyrene
microspheres trapped in circularly polarized, counter-propagating optical
beams, in vacuum. Each oscillator is linearly non-conservative, due to
azimuthal components of momentum associated with inhomogeneous optical spin
Antognozzi  _et al._ (2016); Bliokh  _et al._ (2014); Svak  _et al._ (2018),
and coupled through weak hydrodynamic and optical interactions. Analogous
effects can be induced via birefringence Arita  _et al._ (2020), or phase
difference Rieser  _et al._ (2022); Simpson and Hanna (2006). Below a
critical, threshold power, we observe biased Brownian motion, featuring
correlations which strengthen with increasing optical power. At threshold, a
bifurcation occurs Simpson  _et al._ (2021) and the stable trapping points are
replaced with noisy limit cycles that form robust synchronized states with
characteristic detuning behaviour Pikovsky and Kurths (1997).

Figure 1:  Overview of the experiment. (a) Basic geometry and coordinate
system. Two pairs of counter-propagating circularly polarized laser beams (red
arrows) (wavelength $\lambda=1064$ nm and beam waist $900$ nm) form two
Gaussian standing waves. Two polystyrene particles (nominal radius $a=425$ nm)
are localized in the standing waves in the axial direction ($z$-axis), while
in the lateral direction, they tend to orbit due to the azimuthal spin force
at the ambient pressure 17 mbar (equivalent to an effective viscosity
$\mu\approx 1.15\,\mu$Pa s Beresnev  _et al._ (1990)). Their synchronized
motion emerges fromve the breathing ($\mathrm{QM_{1},\,QM_{2}}$) and center of
mass ($\mathrm{QM_{3},\,QM_{4}}$) quasi-modes. Cumulative azimuthal angles
$\phi_{1}$ and $\phi_{2}$ are used to quantify the level of their
synchronization. (b) Schematic describing the stochastic quasi-modes in the
linear (sub-threshold) regime. Increasing the laser power leads to slightly
larger radii of the lateral particles’ trajectories and suppression of two
quasi-modes $\mathrm{QM_{1},\,QM_{3}}$. (c) Experimentally observed biased
stochastic motion, showing a tendency towards orbital rotation illustrated
using spatial probability densities (PDFs) in the particle displacement basis
and in the quasi-mode basis. The picture depicts excitation of the breathing
mode $\mathrm{QM_{2}}$ for beam separation $d=8.6\,\mu$m, corresponding to a
condition where the threshold power $P_{\mathrm{c}}$ for the center of mass
mode $\mathrm{QM_{4}}$ is larger than $P_{\mathrm{b}}$ for the breathing mode.
(d) At larger powers approaching the threshold, fluctuating limit cycles
occasionally begin to be formed. At $d=8.6\,\mu$m the breathing quasi-mode
$\mathrm{QM_{2}}$ fulfilling $P_{\mathrm{c}}-P_{\mathrm{b}}>0$ (red) is
further developed stable while for $d=8.9\,\mu$m the center of mass quasi-mode
$\mathrm{QM_{4}}$ with $P_{\mathrm{c}}-P_{\mathrm{b}}<0$ (green) further
grows. (e) Short time trajectories of particles for both cases
$P_{\mathrm{c}}-P_{\mathrm{b}}\lessgtr 0$ from (d) demonstrating the rise of
anti-phase (left) and in-phase synchronized motions (right).

The ability to form coherent, coordinated non-equilibrium states, such as the
synchronized states described here, could have applications in sensor arrays
Afek  _et al._ (2022), suppressing phase noise and natural variations in
fundamental frequencies Matheny  _et al._ (2014); Zhang  _et al._ (2015).
These linearly non-conservative oscillators are particular examples of a
broader class of non-Hermitian oscillator Xin  _et al._ (2020); Bergholtz  _et
al._ (2021), characterised by broken time reversal symmetry and a capacity to
exchange energy with the environment. Arrays of such non-Hermitian units can
form topological phases and exhibit exponential sensitivity through the well
know skin effect McDonald and Clerk (2020); Okuma  _et al._ (2020); Budich and
Bergholtz (2020); Afek  _et al._ (2022). These phenomena have been
successfully realized in the classical domain with the use of micro-robotics
Brandenbourger  _et al._ (2019); Fruchart  _et al._ (2021). Optical forces,
such as those considered here, offer a route to realize similar effects,
spontaneously, in the mesoscopic regime. Moreover, under appropriate
conditions, such systems may also act like classical time crystals Yao  _et
al._ (2020); Shapere and Wilczek (2012). Developing appropriate cooling
techniques could take these effects towards the quantum regime and provide
experimental access to mesoscopic quantum dynamic phenomena such as quantum
synchronization or entanglement Roulet and Bruder (2018); Rudolph  _et al._
(2022); Witthaut  _et al._ (2017).

##  II Results

###  II.1 Qualitative Experimental Observations

The geometry of the experiment is graphically represented in Fig. (1)a. The
optical field consists of a parallel pair of counter-propagating Gaussian
beams each of which has the same power, and forms a standing wave normal to
its axis. Circular polarization gives rise to azimuthal components of optical
spin momentum which swirl around the axis of each beam. One particle is
confined in each counter-propagating beam, where it is subject to non-
conservative, azimuthal spin-forces which drive its motion out of
thermodynamic equilibrium. In addition, light scattering between the particles
induces optical binding forces which, in addition to dissipative hydrodynamic
interactions, couple their stochastic motion. The relative strength of these
coupling interactions varies with the separation, $d$ between the beams,
allowing us to tune the form of behaviour manifested in the experiment.

We observe a range of quintessentially non-equilibrium effects ranging from
biased stochastic motion to the formation of synchronized limit cycle
oscillations, Fig 1(b-f). The experimentally observed stochastic motion can be
described in terms of two interconnected pairs of quasi-modes (QMs), whose
properties are described in detail below and in the Supplementary Information.
These pairs of QMs are referred to as the Centre of Mass (CoM) and breathing
(BR) QMs. In combination, they describe in-phase (CoM) and anti-phase (BR)
stochastic orbital rotation, in clockwise and counter-clockwise directions,
Fig. 1b. Each pair of QMs has a threshold optical power, $P_{\mathrm{c}}$ and
$P_{\mathrm{b}}$ for CoM and BR, respectively. Thermal fluctuations combine
with non-conservative forces and excite the QMs to different degrees: the
closer the optical power is to the threshold power of a QM, the more strongly
it is excited and the greater is its mean squared amplitude. In our
experiments, we continuously increase the optical power and make observations
of the stochastic motion produced. The observed behaviour depends, therefore,
on the relative magnitudes of the threshold powers, $P_{\mathrm{c}}$ and
$P_{\mathrm{b}}$, see Fig. 1(b,c). For example, when
$P_{\mathrm{c}}<P_{\mathrm{b}}$, the threshold power of the CoM mode is
approached first as the power increases. This causes the CoM mode to grow most
rapidly until it dominates the observed motion. When
$P_{\mathrm{b}}<P_{\mathrm{c}}$, it is BR that becomes dominant. As described
further below, the difference between the threshold powers,
$P_{\mathrm{c}}-P_{\mathrm{b}}$, oscillates about zero as the beam separation,
$d$, is increased so that $P_{\mathrm{c}}<P_{\mathrm{b}}$ for some separations
and $P_{\mathrm{b}}<P_{\mathrm{c}}$ for others, Fig. 1(d). We are therefore
able to tune the observed behaviour by adjusting $d$. Increasing the power
above one of the threshold powers results in a sudden bifurcation.
Subsequently, each particle executes a limit cycle oscillation. Weak
interaction forces cause these self-sustained oscillations to synchronize,
Fig. 1(f).

In the following sections we first outline some theoretical principles before
applying them to experimental investigations of the sub-threshold (Fig.
1(c,d)) and above threshold regimes (Fig. 1e).

###  II.2 Theoretical Considerations: Generalized Hooke’s Law, Linear
Stability and Limit Cycle Formation in Stochastic Optomechanics

In Supplementary Note we provide a detailed analysis of the general stability
properties and stochastic motion of multi-particle, levitated optomechanical
systems in the linear regime. Below we summarize the results used throughout
the rest of this article.

For many optomechanical systems, including the one studied here, it is
possible to identify a configuration in which the system is at mechanical
equilibrium i.e. a configuration in which the external optical forces vanish
and are locally restoring. For small displacements, the optical force can be
linearly approximated by a generalized Hooke’s law i.e.

$\displaystyle\mathbf{f}\approx-\mathbf{K}\cdot\mathbf{q}\equiv-P\mathbf{k}\cdot\mathbf{q},$
(1a) $\displaystyle K_{ij}=-\frac{\partial f_{i}}{\partial
x_{j}}\Bigg{|}_{\mathbf{q}=0},$ (1b)

where $\mathbf{q}=\mathbf{r}-\mathbf{r}_{0}$ are small displacements with
respect to the coordinates of the mechanical equilibrium, $\mathbf{r}_{0}$.
The stiffness matrix, $\mathbf{K}$, is proportional to the optical power such
that $\mathbf{K}=P\mathbf{k}$, with $\mathbf{k}$ the power normalized
stiffness.

Figure 2:  Quantitative analyses of the experimental results below threshold.
(a) The growth in the auto-correlation of the $x_{\mathrm{b}}$ component of
the excited breathing mode $\langle
x_{\mathrm{b}}(t+\tau)x_{\mathrm{b}}(\tau)\rangle$ with increasing optical
power for the beam distance $d=8.6\,\mu$m giving
$P_{\mathrm{c}}-P_{\mathrm{b}}>0$. Increasing the power increases both the
mean frequency of the oscillation and the time constant governing the loss of
coherence of the oscillation. (b) Due to minor imperfections, the
$x_{\mathrm{c}}$ component CoM auto-correlation, $\langle
x_{\mathrm{c}}(t+\tau)x_{\mathrm{c}}(\tau)\rangle$, also grows slightly, but
the effect is far weaker. (c) The cross-correlation of the $x_{\mathrm{c,b}}$
components of CoM and BR, $\langle
x_{\mathrm{c}}(t+\tau)x_{\mathrm{b}}(\tau)\rangle$ shows weak coupling caused
by minor imperfections in the system. (d) Cross-correlations of the BR
coordinates are $\pi/2$ phase shifted, indicating a tendency toward circular
motion. (e,f) PDFs of the difference between the azimuthal coordinates of both
particles, see Fig. 1(a), $\Delta\phi=\phi_{1}-\phi_{2}$ for beam separations
of $d=8.6\,\mu$m (e) and $d=8.9\,\mu$m (f) showing the emergence of phase
locking approximating to BR and CoM modes, respectively, shown in Fig. 1(e).
(g) Relative Shannon entropy $S_{r}$ as a function of power for the PDFs (e)
and (f).

Two qualitatively distinct cases emerge:

#### Linearly conservative forces:

In this case $\mathbf{K}$ is symmetric with real eigenvalues. The motion of
the system can be described in terms of a discrete, orthogonal set of normal
modes, each satisfying the equipartition theorem, having energy
$k_{\mathrm{b}}T/2$ for any value of the optical power, $P$.

#### Linearly non-conservative forces:

In this case, $\mathbf{K}$ is non-symmetric and its eigenvalues can occur in
complex conjugate pairs Plemmons (1988). Pairs of such eigenvalues are
associated with quasi-modes (QMs) which are not orthogonal, and do not satisfy
equipartition. Each QM has characteristic frequencies that can be approximated
as,

$\omega_{i\pm}\approx\pm\sqrt{\frac{P}{m}}\Re(\lambda_{\mathrm{i}}^{1/2})+i\Big{(}\pm\sqrt{\frac{P}{m}}\Im(\lambda_{\mathrm{i}}^{1/2})+\frac{\xi_{i}}{2m}\Big{)},$
(2)

where $i$ indexes the QM, $\lambda_{\mathrm{i}}$ is the associated complex
eigenvalue of $\mathbf{k}$ and $\xi_{i}$ is the effective drag, directly
proportional to the effective viscosity, $\mu$. In the absence of thermal
fluctuations, these frequencies, $\omega_{i\pm}$, relate to damped
oscillations in which the coupled coordinates spiral into the fixed point
(Supplementary Note). The rate at which they do so, depends on the imaginary
part of $\omega_{i\pm}$, which describes motional damping. By increasing the
power, $P$, $\Im(\omega_{i-})$ can be decreased towards zero before the
changing sign. As it does so, motional damping turns into exponential growth
transforming the inward spiral to an outward spiral, destabilizing the trap.
The condition for $\Im(\omega_{i-})=0$ is,

$\frac{P}{\xi_{i}^{2}}=\frac{\Re(\lambda_{\mathrm{i}})}{m\Im(\lambda_{\mathrm{i}})^{2}}.$
(3)

As the power is increased towards this threshold, the interaction between
thermal fluctuations and the non-conservative force causes the instantaneous
variance, $\langle a_{i}^{2}\rangle$ and decay time of the autocorrelation of
the QM to increase,

$\langle
a_{i}(t+\tau)a_{i}(\tau)\rangle\propto\frac{k_{\mathrm{b}}T}{\Im(\omega_{i-})}e^{i\Re(\omega_{i-})\tau}e^{-\Im(\omega_{i-})\tau},$
(4)

where $a_{i}$ is the amplitude of the QM, and $\Im(\omega_{i-})\rightarrow 0$
(Supplementary Note ). Eventually, the amplitude of the dominant motion
exceeds the range over which the forces are approximately linear. Given
suitable curvature in the force field, stable limit cycles (i.e. isolated,
closed paths in phase space describing self sustained oscillations), or
orbits, can form Simpson  _et al._ (2021) and, ultimately, the fixed point
(i.e. the mechanical equilibrium) of the system is destabilized. Experimental
observations of such a transition are shown graphically in Fig. 1(c,d). We
next apply these principles to our pair of spin-driven oscillators, Fig. 1a.

###  II.3 Sub-threshold behaviour, optical binding between non-conservative
oscillators

First we consider the sub-threshold behaviour, for which the motion remains
within the linear range of the force field, where the generalized Hooke’s law,
Eq. (1) applies, see Fig. 1(c). A detailed account is provided in
Supplementary Note The main theoretical results are summarized below and
compared with experimental demonstrations. We confine attention to the
$xy-$plane, the $z$ motion corresponding to an uncoupled normal mode. In this
plane the displacement coordinates are
$\mathbf{q}=(\mathbf{q}_{1},\mathbf{q}_{2})=(x_{1},y_{1},x_{2},y_{2})$ and the
stiffness matrix for this system have the form,

$\mathbf{K}=\begin{bmatrix}\mathbf{K}^{\prime(1)}&\mathbf{A}\\\
\mathbf{A}&\mathbf{K}^{\prime(1)}\end{bmatrix},\;\;\;\mathbf{K}^{(1)}=\begin{bmatrix}K_{r}&K_{\phi}\\\
-K_{\phi}&K_{r}\end{bmatrix}.$ (5)

Here, $\mathbf{K}^{(1)}$ is the stiffness of a single oscillator, comprising a
single sphere in a counter-propagating, circularly polarized trap. The
diagonal elements, $K_{r}$, quantify the stiffness of the purely attractive
gradient forces and the off-diagonal terms, $K_{\phi}$, are connected with
non-conservative, azimuthal forces deriving from inhomogeneous optical spin
Svak  _et al._ (2018). $\mathbf{K}$ is the stiffness for the pair of
oscillators: the stiffness of each constituent oscillator is slightly modified
by the proximity of its neighbour i.e.
$\mathbf{K}^{\prime}{(1)}\approx\mathbf{K}^{(1)}$, while $\mathbf{A}$
describes the relatively weak coupling between the two particles. Note that
$\mathbf{A}$ is, itself, non-symmetric indicating that the interaction is
intrinsically non-conservative. Its elements oscillate with the separation
between the beams, as is common with conventional binding interactions. A
parametric study of the elements of $\mathbf{K}$, and their dependence on beam
separation, $d$, and particle radius, $a$, is provided in Supplementary Note.

The overall form of $\mathbf{K}$ derives from the inversion symmetry of the
system, and allows separation into two independent oscillators by transforming
to the centre of mass (CoM) and breathing (BR) coordinates
($\mathbf{q}_{\mathrm{c}}$ and $\mathbf{q}_{\mathrm{b}}$ respectively) with,

$\displaystyle\mathbf{q}_{\mathrm{c}}=(\mathbf{q}_{1}+\mathbf{q}_{2})/\sqrt{2},$
(6a)
$\displaystyle\mathbf{q}_{\mathrm{b}}=(\mathbf{q}_{1}-\mathbf{q}_{2})/\sqrt{2},$
(6b)

where $\mathbf{q}_{\mathrm{c/b}}=(x_{\mathrm{c/b}},y_{\mathrm{c/b}})$. This
transformation decouples the system stiffness, $\mathbf{K}$, according to,

$\mathbf{K}=\begin{bmatrix}\mathbf{K}^{\prime(1)}&\mathbf{A}\\\
\mathbf{A}&\mathbf{K}^{\prime(1)}\end{bmatrix}\rightarrow\begin{bmatrix}(\mathbf{K}^{\prime(1)}+\mathbf{A})&0\\\
0&(\mathbf{K}^{\prime(1)}-\mathbf{A})\end{bmatrix}.$ (7)

We refer to these two separate oscillators as CoM and BR oscillators. Each has
two quasi-modes (QMs), with complex conjugate eigenvalues, which, together,
describe stochastic orbital rotation of the coordinates,
$\mathbf{q}_{\mathrm{c}}$ or $\mathbf{q}_{\mathrm{b}}$, about the origin.
These stochastic motions correspond, respectively, to in-phase (CoM) and anti-
phase (BR) circulation of the individual particles, in clockwise or counter-
clockwise directions, see Fig. 1b.

Treating the optical and the hydrodynamic interactions as perturbations, Eq.
(3) gives the difference between the threshold powers for these oscillators
as,

$P_{\mathrm{c}}-P_{\mathrm{b}}\approx-\frac{\xi_{0}^{2}}{m}\Big{(}\frac{4K_{r}\Im(\delta)}{K_{\phi}^{3}}+\frac{9}{2}\frac{a}{d}\frac{K_{r}}{K_{\phi}^{2}}\Big{)},$
(8)

where $\xi_{0}=6\pi\mu a$ is the Stokes drag on a single particle and $\delta$
is a complex scalar quantity derived from elements of $\mathbf{A}$, that
describes the optical interaction. The first term in Eq. (8) is due to optical
coupling and the second is caused by differences in the effective drag for the
CoM and BR QMs (see Supplementary Note). The optical coupling parameter,
$\delta$, oscillates with beam separation, $d$, while the hydrodynamic
interaction decays monotonically with $d$, serving to systematically reduce
the threshold power of CoM relative to BR. As the beam separation, $d$, is
increased the CoM and BR oscillators alternately have the lowest threshold
power, satisfying $P_{\mathrm{c}}-P_{\mathrm{b}}>0$ or
$P_{\mathrm{c}}-P_{\mathrm{b}}<0$ (Supplementary Note).

At low power all four of these QMs have approximately equal energy and the
preference in the sense of stochastic orbital rotation is negligible. As the
power is increased, the stochastic motion is increasingly biased towards
circulation in the sense dictated by the azimuthal spin forces, although the
energy in the CoM and BR QMs remains comparable. For further increases in
power, the energy in the QMs with the lowest threshold power begins to grow
until it becomes dominant. At this point, the observed motion consists of
stochastic rotations of the microspheres around their respective beam axes
which are either in phase (CoM dominant), or anti-phase (BR dominant), but
always in the direction dictated by the azimuthal spin force.

We investigate these phenomena experimentally in Fig. 1c. Fig. 1c shows two
dimensional spatial probability distribution functions (PDFs) for the
particles. On the top two rows, the PDFs are given in displacement coordinates
[i.e. $(x_{1/2},y_{1/2})$] and, on the lower rows, in the QM coordinates
[($x_{\mathrm{c/b}},y_{\mathrm{c/b}}$), Eq. (6)]. These results correspond to
a separation of $d=8.6\,\mu$m, for which $P_{\mathrm{b}}<P_{\mathrm{c}}$, so
that BR grows to dominate, as is clear from the PDFs in the QM basis.

Figure 2 describes the stochastic motion in more detail. Figs. 2a-c show time
dependent correlation functions of $x_{\mathrm{b}}$ and $x_{\mathrm{c}}$, for
increasing optical power. Written in terms of the components of
$\mathbf{q}_{c/b}$ the auto-correlation of these QMs, Eq.4, is,

$\displaystyle\langle x_{\mathrm{c/b}}(t+\tau)x_{\mathrm{c/b}}(\tau)\rangle$
$\displaystyle\propto\frac{k_{\mathrm{b}}T}{\Im(\omega_{i-})}\cos(\Re(\omega_{i-})\tau)e^{-\Im(\omega_{i-}\tau)}$
(9a) $\displaystyle\langle
x_{\mathrm{c/b}}(t+\tau)y_{\mathrm{c/b}}(\tau)\rangle$
$\displaystyle\propto\frac{k_{\mathrm{b}}T}{\Im(\omega_{i-})}\sin(\Re(\omega_{i-})\tau)e^{-\Im(\omega_{i-}\tau)}.$
(9b)

Equation (9a) describes the increased amplitude and coherence of the
stochastically driven oscillations of $x_{\mathrm{c}}$ and $x_{\mathrm{b}}$,
while the cross correlation, Eq. 9b, describes the growing tendency of
$\mathbf{q}_{\mathrm{c}}$ or $\mathbf{q}_{\mathrm{b}}$ to circulate about the
origin Jones  _et al._ (2009).

Figure 3:  Experimental behaviour of non-conservative oscillators, with the
same fundamental frequencies, above threshold where the limit cycles are
formed and synchronized. (a) Trajectories of both particles in $x-y$ plane
corresponding to 10 limit cycle periods. Positions of each particle were
detected by an independent quadrant photodiode detector (QPD). (b) The
correlated motion of both particles along $x$ axis is visible to the naked eye
in records from a fast CMOS camera plotted against time. Blue and red curves
link the centers of particle images. (c) Accumulated phase difference,
$\Delta\Phi=\Phi_{1}-\Phi_{2}$, plotted against time, shows discrete phase
jumps between phase-locked states. (d) Experimental demonstration of the phase
diffusion properties comparing the linear time dependence of the variance of
the accumulated phase of a single oscillator i.e. Var($\Phi_{1}$) (blue) with
the variance in the accumulated phase difference between the synchronized
oscillators i.e. Var($\Delta\Phi$). (e) Probability density function PDF of
the relative phase $\Delta\phi$ of the oscillators obtained by a fast CMOS
camera and mapped on the interval $\langle-\pi,\pi)$.

In Fig. 2a the increase in amplitude and coherence of the BR QM is shown and,
in Fig. 2b, the relative stagnation of CoM, see Eq. (9a). The coupling between
the CoM and BR oscillators is relatively weak, Fig. 2c, but indicates a slight
departure from the ideal symmetry, assumed in Eq. (7), for which CoM and BR
motions would be completely independent. The time dependent autocorrelation of
$x_{\mathrm{b}}$, $\langle x_{\mathrm{b}}(t+\tau)x_{\mathrm{b}}(\tau)\rangle$,
and the cross correlation of $x_{\mathrm{b}}$ with $y_{\mathrm{b}}$, $\langle
x_{\mathrm{b}}(t+\tau)y_{\mathrm{b}}(\tau)\rangle$ are shown in Fig. 2d,
demonstrating the tendency for $\mathbf{q}_{\mathrm{b}}$ to rotate about the
origin, Eqns (9a,9b). This motion corresponds to stochastic rotation of the
individual particles about their beam axes, with a relative phase shift of
$\pi$ rads Svak  _et al._ (2018); Jones  _et al._ (2009), as illustrated in
Fig. 1b.

Figures 2(e,f) gives an analysis of the statistical behaviour of the azimuthal
coordinates of the particles, $\phi_{1/2}$, see Fig. 1(a), in the form of PDFs
at discrete bins of $\Delta\phi_{k}=\phi_{1}-\phi_{2}$, $p(\Delta\phi_{k})$,
as the optical power is increased. For $d=8.6\mu$m, BR grows to dominate the
motion, while CoM is emphasized in for $d=8.9\mu$m. In Fig. 2(g), we plot the
relative Shannon entropy, $S_{r}=1-S/S_{max}$, where $S_{max}=\ln N$, $N$ is
the number of bins in PDF, and $S=-\sum_{k=1}^{N}p(\Delta\phi_{k})\ln
p(\Delta\phi_{k})$ Tass  _et al._ (1998). In this context, $S_{r}$ measures
synchronization strength, taking values between zero and one, where a value of
one indicates perfect synchronization. For the sub-threshold regime, these
values of $S_{r}$ suggest a form of stochastic synchronization, arising prior
to limit cycle formation, as the instability is approached.

###  II.4 Above threshold behaviour, synchronization and phase locking of
limit cycle oscillators

As a dominant QM grows in amplitude, the particles begin to stray further from
the beam axis, where the forces are non-linear, allowing for the formation of
self-sustained, periodic trajectories or limit cycles Pikovsky  _et al._
(2002). Eventually the fixed point destabilizes and each particle forms its
own limit cycle, resembling a circular orbit. These limit cycles exist
independently of one another, and execute a complete cycle in a well defined
time period, $T$, with fundamental frequency, $\Omega=2\pi/T$. The position of
the particle on the limit cycle can be associated with a single scalar
coordinate, the phase, $\phi$. In our system, two limit cycles are formed,
consisting of approximately circular orbits, Fig. 3a.

In general, collections of weakly interacting limit cycles have a tendency to
synchronize. That is, their slightly differing fundamental frequencies are
drawn together so that the ensemble oscillates collectively with a single,
unique frequency Pikovsky  _et al._ (2002). This process is the consequence of
small phase adjustments which accumulate over the course of many time periods.
In this respect, the mechanisms underpinning synchronization differ
fundamentally from those that generate other forms of highly correlated motion
as found, for instance, in conservative optomechanics Svak  _et al._ (2021),
in which the interaction is more direct and the correlation is directly
proportional to a coupling constant.

Figure 4:  Behavior of the relative phase of two interacting limit cycle
oscillators having detuned fundamental frequencies $\Delta f=f_{2}-f_{1}$. (a)
Long-time traces of the accumulated phase difference $\Delta\Phi$ for a series
of detuned oscillators. (b) Shorter-time traces of records (a) illustrating
phase slips of an integer multiple of $2\pi$ radians. As the detuning
increases the time between phase slips decreases until synchronization becomes
impossible. (c) Broadening of the PDF of the absolute phase difference,
$\Delta\phi$ for larger detuning. (d) Quantitative characterization of
weakening of synchronization using the relative Shannon entropy $S_{r}$.

In the mesoscopic regime, synchronization is always accompanied by significant
levels of thermal noise. Since limit cycles are neutrally stable, the phase
diffuses as it advances. In particular, the total change in phase over a time
interval, $\Phi$, has a variance that increases linearly with time Pikovsky
_et al._ (2002). Synchronization of stochastic systems therefore requires that
the interaction forces are strong enough to overcome phase diffusion.
Nevertheless, fluctuations will still give rise to phase slips, in which the
relative phase of oscillators changes abruptly between phase locked states.
This can be likened to Kramers hopping in a potential Rondin  _et al._ (2017);
Mel’nikov (1991) although, in this case, the transitions take place between
non-equilibrium states and a closer analogy is with stochastic motion in a
tilted periodic potential Jung and Risken (1984); Vollmer and Risken (1983).

The above threshold synchronization of our twin oscillators is explored
experimentally in figures (3) and (4). Accompanying simulations and
theoretical comments are provided in Supplementary Note. We take the azimuthal
coordinate of the particle displacement as the phase of the oscillator. This
is less rigorous than the definition obtained from phase reduction but is far
simpler to describe and sufficiently accurate to capture the required
phenomena (see Supplementary Note. In the following discussion we distinguish
between the absolute phase of an oscillator, $\phi_{i}$ with $i=1,2$, which
specifies the position of particle $i$ on its limit cycle and takes values in
the interval $(-\pi,\pi)$, and the accumulated phase, $\Phi_{i}$, which
specifies the total phase difference covered in a given period of time.
Associated with these quantities are the absolute phase difference between the
oscillators, $\Delta\phi=(\phi_{1}-\phi_{2})$ (restricted again to the
interval $(-\pi,\pi)$ and the accumulated phase difference,
$\Delta\Phi=(\Phi_{1}-\Phi_{2})$. Figure 3 describes an established
synchronized state for oscillators with approximately equal fundamental
frequencies. Experimentally measured trajectories are depicted in Fig. 3(a,
b). The accumulated phase difference shows, in Fig. 3(c), a typical noise-
induced phase slip. Figure 3(d) confirms the established linear relationship
between the time dependence of the variance of the accumulated phase,
$\Phi_{1}$, for a single oscillator (blue line). In contrast, the variance in
the difference of the accumulated phases of a pair of synchronized oscillators
saturates quickly demonstrating that the interactions are strong enough to
suppress phase diffusion in the synchronizing pair. The steady-state PDF of
$\Delta\phi$ appears in Fig. (3)(e), showing a sharply peaked phase difference
with relative Shannon entropy $S_{r}=0.39$.

The detuning behaviour obtained when the limit cycles of the oscillators have
differing fundamental frequencies is described in Figure 4. Figures 4a-c show
the effect of varying the second beam waist radius between $w_{0}=1.03\,\mu$m
and $w_{0}=1.06\,\mu$m, holding the first constant at $w_{0}=1.03\,\mu$m. This
has the effect of continuously varying the fundamental frequency of the second
limit cycle oscillator (see Svak  _et al._ (2018)). Figs. 4a,b show the
accumulated phase difference, $\Delta\Phi$, for a series of detuned
oscillators over different time intervals. Over long times, Fig. 4a shows a
steady increase in the accumulated relative phase of the oscillators as the
faster oscillator pulls ahead of the slower. Fig. 4b shows the detailed motion
over shorter time scales. As described previously, the oscillators synchronize
perfectly for short intervals before fluctuations induce a phase slip of
$2\pi$ radians, or an integer multiple. As the detuning increases the time
between phase slips decreases until synchronization becomes impossible. The
effect on $p(\Delta\phi)$ is to lower and broaden the main peak, shifting it
to slightly greater phase differences on average, Fig. 4c. Over the course of
this variation the Shannon entropy changes from 0.35 to 0.17, Fig. 4d.

These experimental results are supported by dynamical simulations (see
Supplementary Note, which shed light on the synchronization mechanism. In
particular, optical interactions alone cannot account for the observed
effects. Dissipative, hydrodynamic interactions act cooperatively to generate
and stabilize synchronized states.

##  III Discussion

In this article we have described the emergence of archetypal, non-equilibrium
behaviour in a non-conservative optomechanical system consisting of a pair of
non-Hermitian oscillators driven by optical spin momentum.

We have shown that stochastic motion becomes progressively more biased,
coherent and deterministic as the optical power is increased. Particular forms
of motion, described by quasi-modes, begin to dominate. Further increases in
power result in collective a Hopf bifurcation and the formation of limit cycle
oscillations which interact and synchronize. This general behaviour is
representative of a far wider class of systems than the particular example
dealt with here.

In addition, our results suggest that hydrodynamic interactions play a role in
the formation of coordinated motion in both the linear and non-linear regimes.
The dependence of hydrodynamic coupling, and therefore dissipation rate, on
the configuration of the system appears to influence the formation of these
non-equilibrium steady states. This effect could be analogous to the minimal
dissipation principle of Onsager Onsager (1931). We note that hydrodynamic
interactions have been inferred experimentally in similar systems Svak  _et
al._ (2021); Rieser  _et al._ (2022). These considerations imply a fundamental
difference between this system, and the paradigmatic, Kuramoto model for
synchronization Gupta  _et al._ (2018), in which the underlying mechanism
relies on reactive forces alone.

More generally, the combination of structured non-conservative forces and
coupled dissipation open up numerous new themes for continuing research in
levitational optomechanics. These avenues range from the development of novel
forms of mesoscale topological matter McDonald and Clerk (2020); Xin  _et al._
(2020), to the experimental exploration of emerging and controversial issues
in the stochastic thermodynamics, of non-equilibrium states, such as the
synchronized states described here. Application of the cooling protocols,
previously applied to conservative systems, could even push these effects
towards the quantum regime allowing experiment to probe the quantum-classical
interface for dynamic phenomena such as limit cycle formation or
synchronization.

Figure 5: Experimental set-up of two pairs of counter-propagating beams
forming standing wave optical traps with circularly polarized light. Particles
are trapped in a small vacuum chamber, dashed square in the inset, placed
between the focusing aspherical lenses AS1,2. Positions of particles in $x-z$
plane are magnified by a telescope and observed by CAM1. Positions of each
particle in $x-y$ plane are independently but synchronously recorded by
quadrant photo-detectors QPD1,2.

##  IV Methods

###  IV.1 Experimental details

In order to optically confine the particles inside a vacuum chamber and
characterize their optical binding, we used a source of infrared laser light
operating at the vacuum wavelength of 1064 nm with low intensity noise
(Coherent Mephisto). We used Thorlabs achromatic doublets with antireflection
coating ACN254-XXX-C (L1 – L6), dielectric mirrors PF10-03 (M1 – M3) and
aspheric lenses C240TME-C with antireflection coating (AS1).

A collimated Gaussian beam from an infrared laser was expanded by a telescope
formed by lenses L1 ($f_{1}=150$ mm) and L2 ($f_{2}=300$ mm) and projected on
a spatial light modulator (SLM) (Hamamatsu LCOS X10468-07). The phase mask
encoded at the SLM diffracted the beam into the $\pm 1$ diffraction orders
that were used to generate the two counter-propagating trapping beams; the
zeroth and higher orders were blocked by a stop placed in the focal plane of
lens L3 ($f_{3}=400$ mm).

The two transmitted $1^{\mathrm{st}}$ – order beams were reflected from prisms
P1 and collimated by lenses L4 ($f_{4}=200$ mm). These lenses formed
telescopes with the lens L3, projecting the SLM plane on the mirrors M2. The
SLM plane was then imaged onto the back focal planes of aspheric lenses AS1
($f=8$ mm, maximal NA = 0.5) by telescopes consisting of lenses L5
($f_{5}=100$ mm) and L6 ($f_{6}=150$ mm).

Two pairs of horizontal counter-propagating laser beams generated by splitting
a single incident beam with a spatial light modulator (SLM) were focused
inside the vacuum chamber by two aspheric lenses with NA=0.5, leading to the
beam waist radii $w_{0}$ adjustable in the range $1-3\,\mu$m. The focal planes
of the four beams created in the trapping region were slightly displaced from
each other along the beam propagation direction $z$ (by approximately 5
$\mu$m, see red lines in Fig. 5) to increase the axial trapping stability
Tatarkova  _et al._ (2002) (see Fig. 5).

Widths of the focused trapping beams in the sample chamber could be controlled
by adjusting the area of the diffraction grating imposed upon the SLM.

Polystyrene particles (Polysciences, mean diameter 850 nm) were dispersed in
isopropyl alcohol and after $\sim 20$ min sonication of the suspension,
droplets containing the particles were sprayed into the trapping region in the
vacuum chamber employing an ultrasonic nebulizer (Beurer IH 50).

We employed two quadrant photo diodes to record the motion of the particles.
Trajectories were recorded for durations of 2 s with 1 MHz sampling frequency.
At the same time we used a fast CMOST camera (I-speed 5 series from IX Camera,
exposure time was set to 1 $\mu$s and the frame rate was 300 kHz) to record
the motion of the particles in $x-z$ plane.

To enable position tracking of the optically trapped and bound particles, the
sample was illuminated by an independent laser beam (Coherent Prometheus,
vacuum wavelength 532 nm) propagating along the $y$-direction perpendicular to
the imaging $xz$-plane. Large beam waist radius $w_{0}=40\,\mu$m and low power
(approximately 5 mW at the sample) of the green illuminating beam ensured its
negligible contribution to the net optical force acting on the particles.
Typically, we recorded at least 100 000 frames from the studied optically
bound structures to obtain sufficiently long trajectories for the analysis of
their motional dynamics.

The off-line tracking of the particle position from the high-speed video
recordings was based on the determination of symmetries in the particle images
Leite  _et al._ (2018). Briefly, since a spherical particle produces an
azimuthally invariant image, we used the shift property of the Fourier
transform and looked for the best horizontal and vertical symmetries in the
particle image, which provided us with the information about the in-plane $x$
and $z$ coordinates.

##  V Acknowledgement

The Czech Science Foundation (GF21-19245K); Akademie věd České republiky
(Praemium Academiae); Ministerstvo Školství mládeže a tělovýchovy
($\mathrm{CZ.02.1.01/0.0/0.0/16\\_026/0008460}$).

##  VI Author Contributions

SHS, OB, and PZ designed and developed the study from the theoretical and
experimental aspects, OB, MD, PJ, and JJ upgraded the experimental setup and
performed the measurements, SHS, OB, and PZ analyzed the experimental data and
compared them to the theoretical results. SHS, OB, and PZ contributed to the
text of the manuscript.

## References

  * Millen and Stickler (2020) J. Millen and B. A. Stickler, Contemp. Phys. 61, 155 (2020).
  * Konopik  _et al._ (2020) M. Konopik, A. Friedenberger, N. Kiesel, and E. Lutz, Europhys. Lett. 131, 60004 (2020).
  * Dholakia and Zemánek (2010) K. Dholakia and P. Zemánek, Rev. Mod. Phys. 82, 1767 (2010).
  * Tebbenjohanns  _et al._ (2019) F. Tebbenjohanns, M. Frimmer, A. Militaru, V. Jain, and L. Novotny, Phys. Rev. Lett. 122, 223601 (2019).
  * Chang  _et al._ (2010) D. E. Chang, C. Regal, S. Papp, D. Wilson, J. Ye, O. Painter, H. J. Kimble, and P. Zoller, Proc. Natl. Acad. Sci. U.S.A. 107, 1005 (2010).
  * Gonzalez-Ballestero  _et al._ (2021) C. Gonzalez-Ballestero, M. Aspelmeyer, L. Novotny, R. Quidant, and O. Romero-Isart, Science 374, eabg3027 (2021).
  * Delić  _et al._ (2020) U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Science 367, 892 (2020).
  * Piotrowski  _et al._ (2022) J. Piotrowski, D. Windey, J. Vijayan, C. Gonzalez-Ballestero, A. d. l. R. Sommer, N. Meyer, R. Quidant, O. Romero-Isart, R. Reimann, and L. Novotny, arXiv preprint arXiv:2209.15326 (2022).
  * Arita  _et al._ (2022) Y. Arita, G. D. Bruce, E. M. Wright, S. H. Simpson, P. Zemánek, and K. Dholakia, Optica 9, 1000 (2022).
  * Liu  _et al._ (2022) J.-Y. Liu, W. Liu, D. Xu, J.-C. Shi, H. Xu, Q. Gong, and Y.-F. Xiao, Phys. Rev. A 105, 053518 (2022).
  * Chauhan  _et al._ (2020) A. K. Chauhan, O. Černotík, and R. Filip, New J. Phys. 22, 123021 (2020).
  * Rudolph  _et al._ (2020) H. Rudolph, K. Hornberger, and B. A. Stickler, Phys. Rev. A 101, 011804 (2020).
  * Rudolph  _et al._ (2022) H. Rudolph, U. Delić, M. Aspelmeyer, K. Hornberger, and B. A. Stickler, Phys. Rev. Lett. 129, 193602 (2022).
  * Monteiro  _et al._ (2020) F. Monteiro, W. Li, G. Afek, C.-l. Li, M. Mossman, and D. C. Moore, Phys. Rev. A 101, 053835 (2020).
  * Hempston  _et al._ (2017) D. Hempston, J. Vovrosh, M. Toroš, G. Winstone, M. Rashid, and H. Ulbricht, Appl. Phys. Lett. 111, 133111 (2017).
  * Sukhov  _et al._ (2015) S. Sukhov, A. Shalin, D. Haefner, and A. Dogariu, Opt. Express 23, 247 (2015).
  * Roichman  _et al._ (2008) Y. Roichman, B. Sun, A. Stolarski, and D. G. Grier, Phys. Rev. Lett. 101, 128301 (2008).
  * Simpson and Hanna (2010) S. H. Simpson and S. Hanna, Phys. Rev. E 82, 031141 (2010).
  * Sukhov and Dogariu (2017) S. Sukhov and A. Dogariu, Rep. Prog. Phys. 80, 112001 (2017).
  * Woerdemann  _et al._ (2013) M. Woerdemann, C. Alpmann, M. Esseling, and C. Denz, Laser Photonics Rev. 7, 839 (2013).
  * Zupancic  _et al._ (2016) P. Zupancic, P. M. Preiss, R. Ma, A. Lukin, M. E. Tai, M. Rispoli, R. Islam, and M. Greiner, Opt. Express 24, 13881 (2016).
  * Gang  _et al._ (1993) H. Gang, T. Ditzinger, C.-Z. Ning, and H. Haken, Phys. Rev. Lett. 71, 807 (1993).
  * Pikovsky and Kurths (1997) A. S. Pikovsky and J. Kurths, Phys. Rev. Lett. 78, 775 (1997).
  * Tanabe and Pakdaman (2001) S. Tanabe and K. Pakdaman, Phys. Rev. E 63, 031911 (2001).
  * Moss  _et al._ (2004) F. Moss, L. M. Ward, and W. G. Sannita, Clin. Neurophysiol. 115, 267 (2004).
  * Laing and Lord (2009) C. Laing and G. J. Lord,  _Stochastic methods in neuroscience_ (OUP Oxford, 2009).
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
  * Gupta  _et al._ (2018) S. Gupta, A. Campa, and S. Ruffo,  _Statistical physics of synchronization_ , Vol. 48 (Springer, 2018).
  * Antognozzi  _et al._ (2016) M. Antognozzi, C. Bermingham, R. Harniman, S. Simpson, J. Senior, R. Hayward, H. Hoerber, M. Dennis, A. Bekshaev, K. Bliokh, _et al._ , Nat. Phys. 12, 731 (2016).
  * Bliokh  _et al._ (2014) K. Y. Bliokh, A. Y. Bekshaev, and F. Nori, Nat. Commun. 5, 1 (2014).
  * Svak  _et al._ (2018) V. Svak, O. Brzobohatỳ, M. Šiler, P. Jákl, J. Kaňka, P. Zemánek, and S. Simpson, Nat. Commun. 9, 1 (2018).
  * Arita  _et al._ (2020) Y. Arita, S. H. Simpson, P. Zemánek, and K. Dholakia, Sci. Adv. 6, eaaz9858 (2020).
  * Rieser  _et al._ (2022) J. Rieser, M. A. Ciampini, H. Rudolph, N. Kiesel, K. Hornberger, B. A. Stickler, M. Aspelmeyer, and U. Delić, Science 377, 987 (2022).
  * Simpson and Hanna (2006) S. H. Simpson and S. Hanna, J. Opt. Soc. Am. A 23, 1419 (2006).
  * Simpson  _et al._ (2021) S. H. Simpson, Y. Arita, K. Dholakia, and P. Zemánek, Phys. Rev. A 104, 043518 (2021).
  * Beresnev  _et al._ (1990) S. Beresnev, V. Chernyak, and G. Fomyagin, J. Fluid Mech. 219, 405 (1990).
  * Afek  _et al._ (2022) G. Afek, D. Carney, and D. C. Moore, Phys. Rev. Lett. 128, 101301 (2022).
  * Matheny  _et al._ (2014) M. H. Matheny, M. Grau, L. G. Villanueva, R. B. Karabalin, M. Cross, and M. L. Roukes, Phys. Rev. Lett. 112, 014101 (2014).
  * Zhang  _et al._ (2015) M. Zhang, S. Shah, J. Cardenas, and M. Lipson, Phys. Rev. Lett. 115, 163902 (2015).
  * Xin  _et al._ (2020) L. Xin, Y. Siyuan, L. Harry, L. Minghui, and C. Yanfeng, Curr. Opin. Solid State Mater Sci. 24, 100853 (2020).
  * Bergholtz  _et al._ (2021) E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Rev. Mod. Phys. 93, 015005 (2021).
  * McDonald and Clerk (2020) A. McDonald and A. A. Clerk, Nat. Commun. 11, 5382 (2020).
  * Okuma  _et al._ (2020) N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Phys. Rev. Lett. 124, 086801 (2020).
  * Budich and Bergholtz (2020) J. C. Budich and E. J. Bergholtz, Phys. Rev. Lett. 125, 180403 (2020).
  * Brandenbourger  _et al._ (2019) M. Brandenbourger, X. Locsin, E. Lerner, and C. Coulais, Nat. Commun. 10, 4608 (2019).
  * Fruchart  _et al._ (2021) M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Nature 592, 363 (2021).
  * Yao  _et al._ (2020) N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, Nat. Phys. 16, 438 (2020).
  * Shapere and Wilczek (2012) A. Shapere and F. Wilczek, Phys. Rev. Lett. 109, 160402 (2012).
  * Roulet and Bruder (2018) A. Roulet and C. Bruder, Phys. Rev. Lett. 121, 063601 (2018).
  * Witthaut  _et al._ (2017) D. Witthaut, S. Wimberger, R. Burioni, and M. Timme, Nat. Commun. 8, 1 (2017).
  * Plemmons (1988) R. J. Plemmons, “Matrix analysis (roger a. horn and charles r. johnson),” (1988).
  * Jones  _et al._ (2009) P. Jones, F. Palmisano, F. Bonaccorso, P. Gucciardi, G. Calogero, A. Ferrari, and O. Marago, ACS Nano 3, 3077 (2009).
  * Tass  _et al._ (1998) P. Tass, M. Rosenblum, J. Weule, J. Kurths, A. Pikovsky, J. Volkmann, A. Schnitzler, and H. Freund, Phys. Rev. Lett. 81, 3291 (1998).
  * Pikovsky  _et al._ (2002) A. Pikovsky, M. Rosenblum, and J. Kurths, “Synchronization: a universal concept in nonlinear science,” (2002).
  * Svak  _et al._ (2021) V. Svak, J. Flajšmanová, L. Chvátal, M. Šiler, A. Jonáš, J. Ježek, S. H. Simpson, P. Zemánek, and O. Brzobohatỳ, Optica 8, 220 (2021).
  * Rondin  _et al._ (2017) L. Rondin, J. Gieseler, F. Ricci, R. Quidant, C. Dellago, and L. Novotny, Nat. Nanotechnol. 12, 1130 (2017).
  * Mel’nikov (1991) V. I. Mel’nikov, Phys. Rep. 209, 1 (1991).
  * Jung and Risken (1984) P. Jung and H. Risken, Z. Phys. B Con. Mat. 54, 357 (1984).
  * Vollmer and Risken (1983) H. Vollmer and H. Risken, Z. Phys. B. Cond. Mat. 52, 259 (1983).
  * Onsager (1931) L. Onsager, Phys. Rev 38, 2265 (1931).
  * Tatarkova  _et al._ (2002) S. A. Tatarkova, A. E. Carruthers, and K. Dholakia, Phys. Rev. Lett. 89, 283901 (2002).
  * Leite  _et al._ (2018) I. T. Leite, S. Turtaev, X. Jiang, M. Šiler, A. Cuschieri, P. S. J. Russell, and T. Čižmár, Nat. Photon. 12, 33 (2018).


(x1) edge (x2)
(x1) edge (x3)
(x3) edge (x4)
(x3) edge (x5)
(x5) edge (x6)
;   % This semicolon ends the \path command.
\end{tikzpicture}\]

The largest locally minimal strong defensive alliance of $T$ is $\{x_1,x_2,x_5,x_6\}$.

For a leaf node $x=x_2,x_4,x_6$, we have  $A_x(0)=0$, $A_x(\hat{1}_b)=1$;
$A_x(\hat{1}_g)=A_x(1_{mg})=A_x(1_{ob})=A_x(1_{og})=-\infty.$ For non-leaf node
$x_5$, we have $d=1$, and hence $C_{\lceil\frac{d}{2}\rceil-1}=C_0=\emptyset$.
$$A_{x_5}(\hat{1}_b) =1+\sum\limits_{x\in \emptyset}{\max\Big\{ A_x(1_{og}), A_x (1_{ob})\Big\} } +
{A_{x_6}(0)}=1+0+0=1; $$

$A_{x_5}(1_{ob})=-\infty$ as   $x_5$ does not have at least $\lceil\frac{d+1}{2}\rceil+1=2$ overprotected (non-leaf) children and $A_{x_5}(\hat{1}_g)=1$ as $x_5$ has only one child ($d=1$);

\begin{equation*}
\begin{split}
A_{x_5}(1_{mg}) & = 1+ \max\{A_{x_6}(\hat{1}_g), A_{x_6}(\hat{1}_b)\}\\
& + \sum\limits_{x\in \emptyset}{\max\{A_{x}(\hat{1}_g), A_{x}(\hat{1}_b), A_{x}(1_{og}), A_x (1_{ob})\} }
\sum\limits_{x\in \emptyset} {A_x(0)}\\
\end{split}
\end{equation*}

$ A_{x_5}(0)= \max\Big\{A_{x_6}(0), A_{x_6}(1_{og}), A_{x_6} (1_{mg})\Big\}=\max\{0, -\infty, -\infty\}=0$
and  $A_{x_5}(1_{og})=-\infty$ as $x_5$ does not have at least $\lceil\frac{d+1}{2}\rceil+1=2$ non-leaf children.
For non-leaf node $x_3$, we have $d=2$;
$A_{x_3}(0)= \max\Big\{A_{x_4}(0), A_{x_4}(1_{og}), A_{x_4} (1_{mg})\Big\}
+\max\Big\{A_{x_5}(0), A_{x_5}(1_{og}), A_{x_5} (1_{mg})\Big\}=0+2=2$ and the associated locally minimal defensive
alliance is $\{x_5,x_6\}$;

\begin{equation*}
\begin{split}
A_{x_3}(\hat{1}_b)& =1+\sum\limits_{x\in \mathcal{C}_1}{\max\Big\{ A_x(1_{og}), A_x (1_{ob})\Big\} } +
\sum\limits_{x\in \mathcal{C}\setminus \mathcal{C}_{1}} {A_x(0)}\\
&= 1+\max\Big\{ A_{x_4}(1_{og}), A_{x_4}(1_{ob})\Big\}+ A_{x_5}(0)\\
\end{split}
\end{equation*}

\begin{equation*}
\begin{split}
A_v(\hat{1}_g) & =1+ \max\{A_{x_5}(\hat{1}_g), A_{x_5}(\hat{1}_b)\}\\
& + \sum\limits_{x\in \emptyset}{\max\{A_{x}(\hat{1}_g), A_{x}(\hat{1}_b), A_{x}(1_{og}), A_x (1_{ob})\} }
+  {A_{x_4}(0)}\\
\end{split}
\end{equation*}

and the associated locally minimal defensive alliance is $\{x_3,x_5\}$;

\begin{equation*}
\begin{split}
A_{x_3}(1_{mg}) & = 1+ \max\{A_{x_5}(\hat{1}_g), A_{x_5}(\hat{1}_b)\}\\
& + \sum\limits_{x\in \mathcal{C}_{1}}{\max\{A_{x}(\hat{1}_g), A_{x}(\hat{1}_b), A_{x}(1_{og}), A_x (1_{ob})\} }
\sum\limits_{x\in \emptyset} {A_x(0)}\\
&=1+\max\{A_{x_5}(\hat{1}_g), A_{x_5}(\hat{1}_b)\}+ \max\{A_{x_4}(\hat{1}_g), A_{x_4}(\hat{1}_b), A_{x_4}(1_{og}), A_{x_4} (1_{ob})\}+0\\
\end{split}
\end{equation*}

and the associated locally minimal defensive alliance is $\{x_3,x_4,x_5\}$;
and $A_{x_3}(1_{og})=-\infty$.
Finally, for the root
node $x_1$, we have $d=2$ and  $A_{x_1}(0)=3$ and the associated locally minimal defensive alliance is $\{x_3,x_4,x_5\}$;

\begin{equation*}
\begin{split}
A_{x_1}(1_{mg}) & = 1+ \max\{A_{x_2}(\hat{1}_g), A_{x_2}(\hat{1}_b)\}+
\end{split}
\end{equation*}

and the associated locally minimal defensive alliance is $\{x_1,x_2,x_5,x_6\}$; and $A_{x_1}(1_{sg})=-\infty$ as $x_1$ does not have
at least $\lceil\frac{d}{2}\rceil+1=2$ non-leaf children.
Therefore, the size of the largest locally minimal strong defensive alliance is
$\max\{A_{x_1}(0),A_{x_1}(1_{mg}),A_{x_1}(1_{og}) \}=4$; and the
largest locally minimal strong defensive alliance is $\{x_1,x_2,x_5,x_6\}$.

Clearly, we can decide if there exists a colorful locally minimal connected defensive alliance of size exactly $k$ with respect to coloring $\chi$ in $\mathcal{O}(n+m)$ time. It is clear that if the graph $G$ does not contain locally minimal connected defensive alliance of size $k$ then we cannot find any colorful locally minimal connected defensive alliance of size $k$. Next, we see that if $G$ is a yes instance then there exists a colorful locally minimal connected defensive alliance with probability $\frac{1}{2^{|S|+|N(S)|}}$. We observe that this is exactly the probability that all the vertices in $S$ are colored green and all the vertices in $N(S)$ are colored red. Since $S$ is a (Locally minimal connected) defensive alliance of size at most $K$, we know that $|N(S)|\leq k^{2}$. This implies that $\frac{1}{2^{|S|+|N(S)|}}\leq  \frac{1}{2^{(k^{2}+k)}}$. Now, we know that if we repeat the coloring $2^{(k^{2}+k)}$ time then the success probability is at least $1-\frac{1}{e}$.