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

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 [1] to the thermodynamic limits of computation [2].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 [3].Supplementing optical with electrostatic forces [4], and combining with an optical cavity [5] results in a reconfigurable experimental system with widely tunable reactive and dissipative forces, capable of supporting dynamical effects across multiple physical regimes [6].
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 [7][8][9][10] with the future promise of macroscopic entanglement [11][12][13].
In the classical domain, the harmonic potentials associated with optical tweezers can confine cooled particles forming high Q oscillators with exquisite force sensitivity [14,15].However, optomechanical systems can exhibit far richer behaviour.Optical forces, including interaction forces, are, in general, non-conservative [16][17][18][19], and can be holographically sculpted [20,21].In combination, arbitrarily structured non-conservative and nonlinear 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 [22], coherence resonance [23], stochastic bifurcations [24] and stochastic synchronization, all of which are exploited by nature in the sense apparatus of animals [25,26].
In this article we demonstrate an archetypal nonequilibrium effect: synchronization of a pair of noisy limit cycle oscillators [27,28].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 [29][30][31], and coupled through weak hydrodynamic and optical interactions.Analogous effects can be induced via birefringence [32], or phase difference [33,34].Below a critical, threshold power, we observe biased Brownian motion, featuring correlations which strengthen with increasing optical power.At threshold, a bifurcation occurs [35] and the stable trapping points are replaced with noisy limit cycles that form robust synchronized states with characteristic detuning behaviour [23].
The ability to form coherent, coordinated nonequilibrium states, such as the synchronized states described here, could have applications in sensor arrays [37], suppressing phase noise and natural variations in fundamental frequencies [38,39].These linearly nonconservative oscillators are particular examples of a broader class of non-Hermitian oscillator [40,41], 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 [37,[42][43][44].These phenomena have been successfully realized in the classical domain with the use of micro-robotics [45,46].Optical forces, such as those considered here, offer a route to realize similar effects,  c) spontaneously, in the mesoscopic regime.Moreover, under appropriate conditions, such systems may also act like classical time crystals [47,48].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 [13,49,50].

A. 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 quasimodes (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 c and P 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 c and P b , see Fig. 1(b,c).For example, when P c < P 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 b < P c , it is BR that becomes dominant.As described further below, the difference between the threshold powers, P c − P b , oscillates about zero as the beam separation, d, is increased so that P c < P b for some separations and P b < P 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 theoret-ical principles before applying them to experimental investigations of the sub-threshold (Fig. 1(c,d)) and above threshold regimes (Fig. 1e).

B. 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.
where q = r − r 0 are small displacements with respect to the coordinates of the mechanical equilibrium, r 0 .The stiffness matrix, K, is proportional to the optical power such that K = P k, with k the power normalized stiffness.
Two qualitatively distinct cases emerge: a. Linearly conservative forces: In this case 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 b T /2 for any value of the optical power, P .
b. Linearly non-conservative forces: In this case, K is non-symmetric and its eigenvalues can occur in complex conjugate pairs [51].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, where i indexes the QM, λ i is the associated complex eigenvalue of k and ξ i is the effective drag, directly proportional to the effective viscosity, µ.In the absence of thermal fluctuations, these frequencies, ω i± , 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 ω i± , which describes motional damping.By increasing the power, P , (ω 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 (ω i− ) = 0 is, (3) As the power is increased towards this threshold, the interaction between thermal fluctuations and the nonconservative force causes the instantaneous variance, a 2 i and decay time of the autocorrelation of the QM to increase, where a i is the amplitude of the QM, and (ω i− ) → 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 [35] 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.

C. 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 q = (q 1 , q 2 ) = (x 1 , y 1 , x 2 , y 2 ) and the stiffness matrix for this system have the form, Here, 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 φ , are connected with nonconservative, azimuthal forces deriving from inhomogeneous optical spin [31].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.K (1) ≈ K (1) , while A describes the relatively weak coupling between the two particles.Note that 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 K, and their dependence on beam separation, d, and particle radius, a, is provided in Supplementary Note.
The overall form of 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 (q c and q b respectively) with, where q c/b = (x c/b , y c/b ).This transformation decouples the system stiffness, K, according to, 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, q c or q 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, where ξ 0 = 6πµa is the Stokes drag on a single particle and δ is a complex scalar quantity derived from elements of 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, δ, 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 c − P b > 0 or P c − P 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 c/b , y c/b ), Eq. ( 6)].These results correspond to a separation of d = 8.6 µm, for which P b < P 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 b and x c , for increasing optical power.Written in terms of the components of q c/b the auto-correlation of these QMs, Eq.4, is, Equation (9a) describes the increased amplitude and coherence of the stochastically driven oscillations of x c and x b , while the cross correlation, Eq. 9b, describes the growing tendency of q c or q b to circulate about the origin [52].
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 b , x b (t + τ )x b (τ ) , and the cross correlation of x b with y b , x b (t + τ )y b (τ ) are shown in Fig. 2d, demonstrating the tendency for q 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 π rads [31,52], as illustrated in Fig. 1b.
Figures 2(e,f) gives an analysis of the statistical behaviour of the azimuthal coordinates of the particles, φ 1/2 , see Fig. 1(a), in the form of PDFs at discrete bins of ∆φ k = φ 1 −φ 2 , p(∆φ k ), as the optical power is increased.For d = 8.6µm, BR grows to dominate the motion, while CoM is emphasized in for d = 8.9µ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 = − N k=1 p(∆φ k ) ln p(∆φ k ) [53].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.

D. 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 selfsustained, periodic trajectories or limit cycles [54].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, Ω = 2π/T .The position of the particle on the limit cycle can be associated with a single scalar coordinate, the phase, φ.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 [54].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 [55], in which the interaction is more direct and the correlation is directly proportional to a coupling constant.
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, Φ, has a variance that increases linearly with time [54].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 [56,57] 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 [58,59].
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, φ i with i = 1, 2, which specifies the position of particle i on its limit cycle and takes values in the interval (−π, π), and the accumulated phase, Φ 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, ∆φ = (φ 1 − φ 2 ) (restricted again to the interval (−π, π) and the accumulated phase difference, ∆Φ = (Φ 1 − Φ 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, Φ 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 ∆φ 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 µm and w 0 = 1.06 µm, holding the first constant at w 0 = 1.03 µm.This has the effect of continuously varying the fundamental frequency of the second limit cycle oscillator (see [31]).Figs.4a,b show the accumulated phase difference, ∆Φ, 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π radians, or an integer multiple.As the detuning increases the time between phase slips decreases until synchronization becomes impossible.The effect on p(∆φ) 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 dynam-ical 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 nonconservative 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 [60].We note that hydrodynamic interactions have been inferred experimentally in similar systems [33,55].These considerations imply a fundamental difference between this system, and the paradigmatic, Kuramoto model for synchronization [28], in which the underlying mechanism relies on reactive forces alone.
More generally, the combination of structured nonconservative 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 [40,42], 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.IV.METHODS

A. 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 ±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 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 µ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 µm, see red lines in Fig. 5) to increase the axial trapping stability [61] (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 ∼ 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 µ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 µ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 [62].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.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.

FIG. 1 :
FIG. 1: Overview of the experiment.(a) Basic geometry and coordinate system.Two pairs of counter-propagating circularly polarized laser beams (red arrows) (wavelength λ = 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 µ ≈ 1.15 µPa s [36]).Their synchronized motion emerges fromve the breathing (QM 1 , QM 2 ) and center of mass (QM 3 , QM 4 ) quasi-modes.Cumulative azimuthal angles φ 1 and φ 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 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 QM 2 for beam separation d = 8.6 µm, corresponding to a condition where the threshold power P c for the center of mass mode QM 4 is larger than P b for the breathing mode.(d) At larger powers approaching the threshold, fluctuating limit cycles occasionally begin to be formed.At d = 8.6 µm the breathing quasi-mode QM 2 fulfilling P c − P b > 0 (red) is further developed stable while for d = 8.9 µm the center of mass quasi-mode QM 4 with P c − P b < 0 (green) further grows.(e) Short time trajectories of particles for both cases P c − P b ≶ 0 from (d) demonstrating the rise of anti-phase (left) and in-phase synchronized motions (right).
FIG. 2: Quantitative analyses of the experimental results below threshold.(a) The growth in the auto-correlation of the x b component of the excited breathing mode x b (t + τ )x b (τ ) with increasing optical power for the beam distance d = 8.6 µm giving P c − P 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 c component CoM auto-correlation, x c (t + τ )x c (τ ) , also grows slightly, but the effect is far weaker.(c) The cross-correlation of the x c,b components of CoM and BR, x c (t + τ )x b (τ ) shows weak coupling caused by minor imperfections in the system.(d) Cross-correlations of the BR coordinates are π/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), ∆φ = φ 1 − φ 2 for beam separations of d = 8.6 µm (e) and d = 8.9 µ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).

1
FIG.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, ∆Φ = Φ 1 − Φ 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(Φ 1 ) (blue) with the variance in the accumulated phase difference between the synchronized oscillators i.e.Var(∆Φ).(e) Probability density function PDF of the relative phase ∆φ of the oscillators obtained by a fast CMOS camera and mapped on the interval −π, π).

FIG. 4 :
FIG. 4: Behavior of the relative phase of two interacting limit cycle oscillators having detuned fundamental frequencies ∆f = f 2 − f 1 .(a) Long-time traces of the accumulated phase difference ∆Φ for a series of detuned oscillators.(b) Shorter-time traces of records (a) illustrating phase slips of an integer multiple of 2π 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, ∆φ for larger detuning.(d) Quantitative characterization of weakening of synchronization using the relative Shannon entropy S r .

FIG. 5 :
FIG. 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.