Noise-Induced Transitions in Optomechanical Synchronization

We study how quantum and thermal noise affects synchronization of two optomechanical limit-cycle oscillators. Classically, in the absence of noise, optomechanical systems tend to synchronize either in-phase or anti-phase. Taking into account the fundamental quantum noise, we find a regime where fluctuations drive transitions between these classical synchronization states. We investigate how this"mixed"synchronization regime emerges from the noiseless system by studying the classical-to-quantum crossover and we show how the time scales of the transitions vary with the effective noise strength. In addition, we compare the effects of thermal noise to the effects of quantum noise.


I. INTRODUCTION
The field of cavity optomechanics deals with systems where light in an optical cavity couples to mechanical motion [1]. Recent experimental progress allows to move on from the investigation of a single optomechanical system to several coupled optomechanical systems. There are already first experiments that involve a few mechanical and optical modes [2][3][4][5], exploiting them for wavelength conversion, phonon lasing and efficient cooling. Such few-mode optomechanical setups have been the subject of an increasing number of theoretical proposals, on topics such as efficient state transfer [6], two-mode squeezing [7], back-action evading measurements [8], entanglement [9][10][11] or Landau-Zener dynamics [12,13].
Most notably, optomechanical arrays provide a platform to study synchronization of mechanical oscillators. This was initially pointed out in Ref. [31]. Synchronization is a well known phenomenon in many different branches of science [32] and typically arises whenever there are stable limit-cycle oscillations. However, we note that synchronizationlike phenomena have been recently studied also in the context of linear oscillators dissipating into a common bath [33]. Optomechanical systems exhibit a Hopf bifurcation and can be optically driven into mechanical limit-cycle oscillations [34][35][36][37]. These self-oscillations have also been analyzed theoretically in the quantum regime [38][39][40]. The theoretical description of the synchronization dynamics in optomechanical arrays has initially focussed on the classical regime [31,41]. More recent insights into this regime include the pattern formation of the mechanical phase field in larger optomechanical arrays [42]. In the quantum regime, it was found [43] that quantum noise can drive a sharp nonequilibrium transition towards an unsynchronized state in an extended array, even for optomechanical systems with identical frequencies. Further general insights into quantum synchronization were gained in a model system of one van-der-Pol oscillator coupled to an external drive [44] or two coupled van-der-Pol oscillators [45], which can serve as a rough approximation to an actual optomechanical system. A number of more recent works have explored quantum synchronization on the more conceptual level [46], as well as in various other physical systems, such as e.g. trapped atoms and ions [47][48][49][50], qubits [51][52][53], and superconducting devices [54]. The relation of synchronization and correlations in the quantum-to-classical transition was studied in a system of coupled cavities containing a nonlinearity [55]. Notably, it is still challenging to define a good measure for quantum synchronization [49,56,57].
Only recently synchronization of two nanomechanical oscillators in the classical regime was demonstrated experimentally. This has been achieved using optomechanical systems with coupled micro-disks [14] (Fig. 1(b)), as well as in an experiment involving an optical racetrack cavity coupled to two mechanical oscillators [58] (Fig. 1(c)), and also in a setup using nanoelectromechanical systems [59]. In a recent first step towards larger arrays, up to seven optically coupled micro-disks were used to demonstrate the expected phase noise reduction due to synchronization [60]. This 1/N phase noise reduction with the number of coupled systems N is considered to be one of the main prospects of synchronized nanomechanical arrays. Indeed, recognized from the very beginning with the synchronization of pendulum clocks [61], synchronization has the potential to improve time-keeping and frequency stability. Examples where different types of synchronization have been applied or suggested for application are for frequency stabilization of high power lasers by coupling to a more stable, low power laser [62] and for secure communication in connection with chaos [63]. For a more complete overview and also the many applications to biology see e.g. [32,64,65]. Micro-disk oscillators that support optical whispering gallery modes which couple evanescently to each other [14,60], (c) optical racetrack resonator coupled to two nanomechanical oscillators [58], (d) optomechanical crystal structure with two optomechanical cells in the vicinity of each other, allowing for optical and mechanical coupling [20].
In this work we study the effects of quantum and thermal noise on the synchronization of two optomechanical systems. We focus on a bistable synchronization regime that either exists already in the absence of noise, or is induced by it. Bistabilities in quantum systems have been investigated before [69][70][71] and noise-induced bistabilities are known in several other systems in biology and chemistry [66][67][68], as well as in physics [72][73][74]. We find and discuss noise-induced bistable behaviour now in the context of optomechanical (quantum) synchronization.
This manuscript is organized as follows: We begin with a brief review of classical synchronization in the absence of noise, Sec. II. Then, we introduce our model in Sec. III, explain our methods in Sec. IV, and state our main results in Sec. V. We note that both quantum and thermal noise lead to similar effects, but start out with the investigation of quantum noise effects. Therefore, in Sec. VI, we simulate the full quantum behaviour of the system, in contrast to the previous investigation presented in Ref. [43]. We find a regime of "mixed" synchronization with two stable synchronization states, and we explore its classical-to-quantum crossover in Sec. VII. In Sec. VIII, we give an overview of the different synchronization regimes. Finally, in Sec. IX, we discuss the effects of thermal noise as compared to quantum noise. This is important to gauge the potential of observing the quantum noise effects discussed here in future experiments.

II. BRIEF REVIEW: CLASSICAL SYNCHRONIZATION OF OPTOMECHANICAL OSCILLATORS
In this section we briefly review the concepts of classical synchronization of optomechanical systems in the absence of noise. This will set the stage for the discussion of our results on synchronization in the presence of thermal and quantum noise. A widely studied model for synchronization is the Kuramoto model [75,76], which describes a set of coupled phase oscillators. Each phase oscillator has a phase φ i and an intrinsic frequency ω i , and couples to the other oscillators via the phase difference. For two oscillators, the two corresponding phase equations collapse into one equation for the relative phase δφ = φ 2 − φ 1 , The two oscillators are synchronized if their respective phase velocities become equal,φ 1 =φ 2 , i.e. δφ = const. From Eq. (1) one finds the synchronization threshold: The two oscillators are synchronized if the coupling k exceeds the natural frequency difference, |k| > |ω 2 − ω 1 |. Following this condition, oscillators with identical intrinsic frequencies ω i are always synchronized. There is only a single stable value of δφ in the synchronized regime. For k → +∞, this value approaches zero, δφ → 0, whereas δφ → π for k → −∞. Although synchronization appears in a large range of systems which are very different in terms of microscopic parameters, their behaviour can often be captured by effective phase equations of the Kuramoto-type [32,76]. In the context of optomechanics, the mechanical motion of a single optomechanical system near the Hopf bifurcation can be described with a phase and an amplitude equation [36]. Starting from these equations, an effective Kuramototype model for coupled optomechanical systems has been derived [31]. This model describes arrays of arbitrary many optomechanical cells with arbitrary intrinsic frequencies. For two oscillators, the phase equations of this model reduce, again, to a single equation for the phase difference Here, S 1 and S 2 are effective parameters depending on the microscopic parameters of the underlying optomechanical systems [31,43]. Note that the S 1 -term was added to the model only later [43], and accounts for the change of the intrinsic frequencies ω i with the mechanical oscillation amplitude. Recently, the resulting full Hopf-Kuramoto model has been used to study pattern formation in 2D arrays of optomechanical systems [42]. In contrast to the original Kuramoto model Eq. (1), the optomechanical Hopf-Kuramoto model Eq.
(2) includes a term that involves sin(2δφ). Rewriting Eq. (2) in terms of an effective potential, δφ = −U (δφ) [31], this term corresponds to the appearance of a second minimum close to δφ = π which can co-exist with the minimum close to δφ = 0. This allows optomechanical systems to synchronize not only in-phase (0-synchronization), δφ → 0, but also anti-phase (π-synchronization), δφ → π. For the effective potential there are three different possibilities, depending on the parameters S 1 and S 2 : (i) It has a single minimum close to δφ = 0 which leads to 0-synchronization only, (ii) it has a single minimum close to δφ = π which leads to π-synchronization only, or (iii) both minima appear simultaneously in the effective potential, such that the initial conditions determine whether 0-or π-synchronization occurs. These three regimes are schematically shown in Fig. 2(b) and (c) and discussed in Sec. III. In general, for optomechanical arrays an effective potential for the phases φ i does not exist [42]. However, in the case of two oscillators only, the system can be described with one degree of freedom, i.e. the relative phase δφ, and an effective potential for δφ can always be constructed. Below, we make use of the existence of this effective potential to give an intuitive understanding of synchronization even in the presence of noise.

III. MODEL
We now introduce the model that we investigate throughout the rest of this manuscript. We study two optomechanical systems which are mechanically coupled, see Fig. 1(a). Our goal is to analyze the synchronization behaviour in the presence of quantum and thermal noise. Each optomechanical system consists of a driven optical mode (â j ) coupled to a mechanical mode (b j ) via radiation pressure. In a frame rotating at the laser drive frequency ω L , the Hamiltonian of each optomechanical system is Here, ∆ = ω L − ω c denotes the detuning from the cavity resonance ω c , Ω is the resonance frequency of the mechanical mode, g 0 denotes the optomechanical single photon coupling strength, and α L is the laser driving strength. Note that the optical and the mechanical systems experience damping at a rate κ and Γ, respectively, and these are described by adding to the Hamiltonian a system-bath coupling in the usual manner:Ĥ . Driving a single optomechanical system with a blue-detuned laser causes self-oscillations of the mechanical resonator when the optomechanically induced negative damping Γ opt is larger than the intrinsic damping of the oscillator Γ [36]. In Fig. 2(a) we show this threshold of self-oscillations. Throughout this work we only consider parameters such that the single optomechanical systems are above this threshold. Self-oscillations at ∆ < 0 occur due to the static optomechanical shift which leads to an effective blue detuning, i.e. ∆ eff > 0. These limit-cycle oscillations, in the absence of noise, can effectively be described by a fixed amplitude and a phase and are treated as a prerequisite for synchronization throughout this work. We consider two self-oscillating optomechanical systems that are coupled mechanically with strength K, such that the total Hamiltonian of the system (except for the dissipative part) readŝ Experimentally, this coupling between the mechanical oscillators can also be mediated by an optical coupling, cf. Fig. 1 In recent experiments [14,58,60], a single joint optical mode was employed to couple the mechanical oscillators. However, this mode served a dual purpose in creating the limit cycles via a blue-detuned drive and providing the coupling. Using an additional, independently driven optical mode for the coupling would allow to tune the effective (optically induced) coupling independently from the laser drive used to create the limit cycles.
In experiments, the typical mode of operation is to have the two optomechanical oscillators at slightly different intrinsic mechanical frequencies. These start out un-synchronized but can synchronize upon changing some parameter (e.g. the laser drive strength, the detuning, or potentially the coupling). This is schematically shown in Fig. 2 where we indicate the synchronization regimes in the absence of noise. In accordance to the Hopf-Kuramoto model, cf. Sec. II, there are unsynchronized regions and three different synchronization regimes: (i) 0-synchronization, (ii) π-synchronization, and (iii) classical bistable synchronization where the type of synchronization depends on the initial conditions. Note that for different intrinsic mechanical frequencies, δΩ = Ω 2 − Ω 1 = 0, the relative phase δφ is not exactly 0 or π but only close to one of these values and varies within the synchronization regime.
However, in the presence of noise it is already interesting to investigate the behaviour even for identical frequencies. In particular, for large-scale optomechanical arrays of identical oscillators, it has been found that there is a synchronization transition as a function of noise strength [43]. Moreover, in the present article we will focus on noise-induced transitions between various synchronization states. The observation of this physics does not rely on whether there is an actual synchronization transition at lower values of the coupling. Therefore, in most of our analysis, we will focus on identical systems, i.e. we assume all the parameters to be equal in both systems. In Fig. 2(c) we schematically show the synchronization regimes for identical optomechanical systems and for a larger mechanical coupling K than in Fig. 2(b), but still in the absence of noise. It is important to note that both the mechanical detuning δΩ and the coupling K have an influence on this diagram: Not synchronized regions can become synchronized and synchronization regime borders are shifted.
We will comment on the dynamics of two coupled optomechanical oscillators with different frequencies in Sec. VIII.

IV. METHODS
In this work we use Langevin equations and quantum jump trajectories to study the system described by Hamiltonian (4) in the presence of (quantum) noise. Here we want to briefly present both approaches and discuss their respective advantages and problems.
Most of our results are computed with semi-classical Langevin equations. They are obtained by first deriving quantum Langevin equations from Hamiltonian (4) using input-output theory [88]. We then adopt a semi-classical approach by turning the quantum Langevin equations into classical Langevin equations for the complex amplitudes α j and β j , where the noise terms mimic the quantum-mechanical zero-point fluctuations. This can be understood as a variant of the "truncated Wigner approximation". The semi-classical equations are: The corresponding equations for the second optomechanical system can be obtained from Eqs. (5) by exchanging the indices 1 ←→ 2. Here, α jin (t) and β jin (t) represent the optical and mechanical input noise, given by Gaussian stochastic processes. Since complex numbers commute, they obviously cannot correctly fulfill the input-output quantum noise correlators, â † jin (t)â jin (t ) = 0 and â jin (t)â † jin (t ) = δ(t − t ). Instead, α jin and β jin are made to mimic quantum noise by fulfilling α jin (t)α * jin (t ) = α * jin (t)α jin (t ) = δ(t − t )/2 (and likewise for β jin for T = 0). This approach allows to study also large parameter ranges with relatively low computational effort.
Deep in the quantum regime it is initially not clear that these semi-classical Langevin equations describe the correct physical behaviour. In order to verify the qualitative effects observed with Langevin equations, we also present a few results obtained with quantum jump trajectories [77,78], i.e. an "unraveling" of the Lindblad master equation. Applying this method, the fully quantum system is simulated on an appropriately truncated Hilbert space. Notably, it allows to work with wave functions, in contrast to the formalism of the Lindblad master equation which requires the use of density matrices. Hence, simulations of larger Hilbert spaces become feasible. In addition, quantum jump trajectories give access to additional observables, for instance the full counting statistics. For these reasons, quantum jump trajectories have been applied in many different contexts. In the field of quantum synchronization, such methods have been used to study synchronization of qubits [51]. In the field of cavity optomechanics, they have been employed for instance to discuss QND measurements, photon statistics and single photon optomechanics [79][80][81][82]. Another motivation is the recent experimental detection of individual phonons in optomechanical systems [83] -a step towards monitoring full quantum jump trajectories in experiments.
To explain this approach [77,78], let us consider for a moment photon decay in cavity 1. The "unraveling" discussed here corresponds to the physical setup of placing a single photon detector at the output port of the cavity. At each time step δt, the probability of a single photon to leak out of the cavity (at temperature T = 0) is given by p = κδt â † 1â 1 . In the case of a photon loss it is detected at the output port and the wave function is updated |Ψ(t + δt) =â 1 |Ψ(t) . If no photon was lost to the environment, the wave function evolves in time according to a (non-Hermitian) Hamiltonian that is obtained by adding a term −i (κ/2)â † 1â 1 . This additional term accounts for the information gained about the system by not observing a photon [84]. In both cases the state |Ψ(t + δt) has to be normalized before proceeding to the next time step. The treatment of photon loss in cavity 2 and of phonon losses works analogously. This simulation approach naturally accounts for quantum fluctuations.
Although it would be favourable to use quantum jump trajectories throughout the whole study, this approach is only computationally feasible whenever the needed Hilbert space remains sufficiently small. This leads to severe limitations in the choice of parameters and especially gives no access to the full classical-to-quantum crossover that is studied in Sec. VII. In Ref. [38] it has been shown for a single optomechanical system that semi-classical Langevin equations produce good agreement with the full quantum theory. A systematic comparison for two coupled optomechanical systems is not possible, since the number of required Hilbert space dimensions is squared as compared to the single optomechanical system. We comment on this in Sec. VII A (cf. Fig. 4(a) and the discussion below).
Note that quantum jump trajectories serve here as a numerical approach only. In experiments, single photon and phonon detection is not necessary. Instead, the mechanical oscillators could be measured using standard homodyning techniques.

V. MAIN RESULTS
In this section we briefly state our main results, which will be discussed in the following sections. Most notably, it is known even from the classical theory that two coupled optomechanical oscillators can be in either one of two synchronization states (with a phase difference near 0 or near π). We find a regime of "mixed" synchronization, where transitions between 0-and π-synchronization occur (Sec. VI). These transitions are driven by (quantum or thermal) noise and cannot be found in the classical, noiseless situation. The average residence times in the two synchronization states can differ and their ratio varies with the system parameters (Sec. VII). Investigating the classical-to-quantum transition, we find that mixed synchronization can evolve from two different regimes in the classical, noiseless limit: (i) there are already two stable synchronization states but in the absence of noise there are no transitions, (ii) there is only one stable synchronization state and only the presence of noise leads to a second stable solution. Although the first sections are devoted to the investigation of quantum noise effects, we note that we find similar effects for thermal noise acting on the mechanical resonator. However, quantitative differences remain due to the different nature of the noise source (Sec. IX). We find that quantum noise effects should dominate over thermal noise effects if the optomechanical cooperativity is sufficiently large, and a large value of g 0 is not necessarily required.

VI. MULTISTABLE QUANTUM SYNCHRONIZATION
First, we start by analyzing two coupled optomechanical systems, cf. Fig. 1(a), deep in the quantum regime. Quantum jump trajectories are used to initially investigate the full quantum dynamics. In the following, we consider a small-amplitude limit cycle and a large single-photon coupling strength, g 0 /κ = 1. This ensures that quantum fluctuations can potentially have a large impact on the system's dynamics. Furthermore, small photon and phonon numbers are necessary to keep the numerical simulations tractable, since they determine the size of the truncated Hilbert space.
To study synchronization we focus on the relative phase δφ = φ 2 − φ 1 between the optomechanical oscillators. Classically, δφ allows to identify synchronization (δφ = const) and distinguish the different synchronization regimes (0and π-synchronization). To extract the relative phase at each time step of a quantum jump trajectory, we use that the motion of the uncoupled self-oscillators in the absence of noise can effectively be described by b j (t) ≈ |b|e −i(Ωt+φ0j ) , with an initial random phase φ 0j . Thus, the correlator b † 1b 2 ∼ e −iδφ can be interpreted as a measure for the relative phase δφ. In Fig. 3 we show the probability density of the relative phase obtained from the distribution of δφ and a typical sample of the corresponding phase trajectory. We find a 0-synchronized regime, Fig. 3(a), where the relative phase is predominantly close to δφ = 0. Similarly, we find a π-synchronized regime, Fig. 3(b), for different parameters. As expected, noise prevents perfect synchronization and thus the corresponding maximum in the probability density has a finite width. Likewise, the trajectories show fluctuations around either δφ = 0 or δφ = π. These two synchronization regimes have an analog in the classical, noiseless limit, cf. Sec. II. In addition, we find another regime where the probability density of the phase difference has maxima close to both δφ = 0 and δφ = π, see Fig. 3(c). The corresponding trajectory shows that in this regime transitions between the two synchronization states occur. We call this regime "mixed" synchronization. This is in contrast to the classical, noiseless result, where the system does not change its synchronization state with time, not even in the regime of bistable synchronization.
Similar to the classical, noiseless system [31] we can understand these results in terms of an effective potential for the relative phase δφ. This is illustrated in the bottom row of Fig. 3. It offers an intuitive understanding of our results. The noise makes the phase fluctuate around the stable point(s) near δφ = 0 (and)or δφ = π. In the case of mixed synchronization the effective potential has two minima, near δφ = 0 and δφ = π, and quantum noise drives transitions between those two states. The probability to be in either the 0-or π-synchronized state is given by the area of the corresponding maximum in the probability density. It is associated to the depth of the minimum in the effective phase potential. The ratio between the probabilities of the two states can assume arbitrary values, depending on the parameters of the system. Note that the absence of two distinct peaks in the probability density does not necessarily mean that the system explores the region around one synchronization state only. In fact, sometimes the trajectories themselves can reveal short stretches of phase dynamics in the vicinity of the other state. Nevertheless, if the fraction of time spent in that other state is short, a second peak will not be visible.
From this interpretation in terms of an effective potential, it is clear that the classical bistable synchronization regime where two potential minima already exist (but no transitions can occur), turns into a mixed synchronization regime (showing transitions) in the presence of noise. However, our analysis in the following section reveals that mixed synchronization can also appear for parameters where classically there is only a single stable synchronization state.

VII. CLASSICAL-TO-QUANTUM CROSSOVER
At the moment, the single-photon coupling strength g 0 is still comparatively small in almost all experiments. As a consequence, quantum effects have only been observed in the linearized regime where only Gaussian states are produced. In the most promising cases [85,86], g 0 /κ can take values up to 10 −2 and g 0 /Ω up to above 10 −4 . Much larger values have been reported for experiments with cold atoms (up to around g 0 /κ ∼ 1) [87], but these do not operate in the "good cavity limit", i.e. one has κ Ω in those experiments, precluding the observation of single-photon strong coupling effects. Nevertheless, as experiments are approaching the single-photon strong coupling regime, they will gradually see increasingly strong effects of quantum fluctuations even in non-linear dynamics. In this section, it is our aim to explore the crossover between the classical regime (small g 0 /κ) and the quantum regime (large g 0 /κ) with respect to optomechanical quantum synchronization. We will focus our investigations mostly on the mixed synchronization regime which is the most interesting one, as we can have noise-induced transitions between the synchronization states. In this section, we will disregard thermal noise, i.e. we assume temperature T = 0, such that only quantum noise is present. The "classical" regime we are discussing here is therefore the noiseless limit of the classical equations of motion. We will later remark on the effects of thermal noise (Sec. IX).
To explore the classical-to-quantum crossover, we want to effectively vary while making sure to keep all the classical predictions unchanged. Notably, the optomechanical coupling strength g 0 = ∂ω c /∂x depends on . For the simplest case of a Fabry-Pérot cavity of length L with one static and one movable mirror the optomechanical coupling is g 0 = ω c (0)x ZPF /L ∼ √ , where x ZPF = /2mΩ denotes the zero-point fluctuations of the mechanical oscillator of mass m. As discussed in Ref. [38], the "quantum parameter" g 0 /κ ∼ √ can thus be varied to effectively change the quantum noise strength. This implies that all classical ( -independent) parameters (κ/Ω, Γ/Ω, ∆/Ω, K/Ω, g 0 α L /Ω 2 ) are kept fixed while g 0 is modified. To see that the quantum parameter has indeed this anticipated effect, it is very helpful to rescale the amplitudes α j and β j , Eqs. (5), such that they tend to a well-defined finite value in the classical limit g 0 /κ → 0 [38]. This can be achieved by definingα j = g 0 α j andβ j = g 0 β j . While |α 1 | 2 gives the number of photons (a "quantum-mechanical" quantity), the rescaled version |α 1 | 2 = g 2 0 |α 1 | 2 ∼ |α 1 | 2 is proportional to the energy ω c |α 1 | 2 inside the cavity (i.e. a classical quantity). A similar argument applies toβ j . With this rescaling, we find the following equations:α Hereα L = g 0 α L is the rescaled laser-driving amplitude that we keep fixed while varying g 0 . The important observation here is that g 0 has been completely eliminated from the equations and now only appears in the strength of the quantum noise: we now have α jin (t)α * jin (t ) = g 2 0 δ(t − t )/2 (and likewise forβ jin ) which indeed vanishes in the classical limit of g 0 ∼ √ → 0.   Fig. 3(c), a rotating wave approximation has been used for the mechanical coupling.
If we consider mechanical thermal noise at finite temperature T (as discussed in Sec. IX), we have β jin (t)β * jin (t ) = g 2 0 (n th + 1/2)δ(t − t ). Here, n th denotes the thermal occupancy of the bath coupled to the mechanical oscillator. The product g 2 0 n th ∼ n th becomes independent of in the classical limit k B T Ω, where n th ≈ k B T /Ω and k B is Boltzmann's constant. In summary, the rescaled equations (6) nicely show how an increase of g 0 can indeed be viewed solely as an increase of the strength of "quantum noise" in our system.
Note that if g 0 → 0 we have to increase the laser driving strength α L to keepα L = g 0 α L constant which means that the total light energy circulating inside the optical cavity is constant. Due to → 0, this corresponds to an increasing number of photons inside the cavity and thus increases drastically the size of the Hilbert space necessary for the full quantum simulation. Therefore, quantum jump trajectories are not suitable to explore the full quantum-to-classical crossover and we have to apply Langevin equations.

A. Synchronization as a Function of Quantum Noise Strength
In all studies of synchronization, one needs to select suitable quantities that measure the degree of synchronization. In this context, it is important to note that for any finite noisy system (subject to quantum and/or thermal noise), there is no sharp synchronization transition and correspondingly no unambiguous measure that displays nonanalytic behaviour at any parameter value. Before proceeding to our results, we summarize and discuss the synchronization measures adopted here which have to be combined to obtain a full picture: (i) the probability density of δφ, (ii) the average phase factor e −iδφ , and (iii) individual trajectories.
The probability density of the relative phase δφ is either mostly flat (no synchronization) or, as shown in Fig. 3, centered predominantly around 0 or π, or it may have two peaks, according to the synchronization regime. In the following, we aim to compress the information contained in the phase distribution into one quantity and calculate the Its real value, Re[C] = cos δφ , distinguishes the three different synchronization regimes: (i) cos δφ ≈ 1 for 0-synchronization, (ii) cos δφ ≈ −1 for π-synchronization, (iii) intermediate values of cos δφ for mixed synchronization. However, this measure has its limitations: When δφ is more or less evenly distributed (no synchronization) cos δφ ≈ 0, this cannot be distinguished from a mixed synchronization situation where almost equal time is spent in the 0-and π-synchronized states. Furthermore, even in the absence of synchronization (and even in the noiseless case) the phase δφ may spend an increased amount of time around certain values. This leads to a finite value of cos δφ , and similarly would also show up in the phase distribution. A solution to this problem is to simultaneously look at a part of the corresponding trajectory, where synchronization can easily be distinguished from an unsynchronized state. Instead, one could start to use more complicated correlators, e.g. b † 1b † 1b 2b2 ∼ e −i2δφ . This correlator allows to distinguish unsynchronized states from synchronized states, but an additional measure is needed to distinguish 0-from π-synchronization. Note that the imaginary part of the above defined correlator, Im[C] = sin δφ , can be used as well. However, in the special case of identical optomechanical systems Im[C] ≈ 0 due to the symmetry of the system.
In Fig. 4(a) we show how cos δφ varies as a function of the quantum parameter g 0 /κ. We chose parameters that lead to a mixed synchronization regime for larger values of g 0 /κ. In the deep quantum regime (g 0 /κ → 1) the probability P 0 to find the system in the 0-synchronized state is larger than the probability P π to find the phase around π, such that cos δφ > 0. Going towards smaller values of g 0 /κ, the ratio P π /P 0 increases, such that eventually cos δφ < 0. Finally, we should reach the classical (noiseless) limit, when g 0 /κ → 0. It turns out that, for the parameters adopted here, the classical solution always ends up in the π-synchronized state, independent of initial conditions. This implies that there is only one minimum in the effective potential. We conclude that the system has turned from a mixed synchronization regime into a purely π-synchronized regime as the quantum parameter was reduced. This cannot be understood in the simple picture of a noise-independent phase potential. We will discuss this kind of behaviour in more detail later on (Sec. VII C).
Note that the calculations for Fig. 4 have been performed using Langevin equations; although in the deep quantum regime, two data points were also acquired with quantum jump trajectories (red triangles). They are shifted as compared to the Langevin results, but show the same trend. We expect that the difference between the Langevin and quantum jump results decreases for smaller quantum parameters g 0 /κ, as it is the case for a single optomechanical system [38]. Since smaller g 0 /κ require a significantly larger Hilbert space for the quantum jump simulations, we cannot compute this for our coupled system. For large quantum parameter g 0 /κ ∼ 1, qualitative differences between the full quantum model and the Langevin equations have been already observed in Ref. [38] as well. Especially a shift of the detuning ∆ was reported, that could be determined numerically also for our system. Taking this detuning shift into account would improve the agreement of our results, although differences remain. Here, we show the uncorrected outcomes of both approaches.

B. Residence Times in the Mixed Synchronization Regime
The measure cos δφ quantifies the fraction of total time spent in 0-synchronized parts as compared to πsynchronized parts. However, it does not provide any information about the rate of transitions between the two synchronization states. Based on the effective potential picture, one might expect the transition rates to be determined by the barrier height and the noise strength. In particular, for larger effective noise strengths g 0 /κ, we expect more frequent transitions. This behaviour is qualitatively visible in Fig. 4(b) and (c). However, as concluded in the previous section, the potential picture is not sufficient to explain all observations. Thus, we now turn to a quantitative analysis and discuss how the transition rates between 0-and π-synchronized states (i.e. the typical residence timesτ ) change during the classical-to-quantum crossover. We extract the fluctuating residence times from the phase trajectories and obtain their distribution. The results are shown in Fig. 5(a) and (b) for the 0-and π-synchronized states, for two different quantum parameters g 0 /κ. In all cases the probability densities decay exponentially with time, ∼ e −τ /τ , and the average residence timeτ is obtained from a fit to the distribution. The extracted average residence times τ 0 and τ π for the two states are shown in Fig. 5(c) as a function of the quantum parameter. Note that the ratio of residence times equals the ratio of probabilities, τ 0 /τ π = P 0 /P π . Nevertheless, the dependence of the times τ 0 , τ π on g 0 /κ reveals new information.
We have chosen parameters such that at g 0 /κ = 1 both 0-synchronized and π-synchronized parts have almost equal average residence times. This corresponds to cos δφ ≈ 0 and P 0 ≈ P π . Furthermore, in the classical limit g 0 /κ = 0 the system is π-synchronized only. When the classical limit g 0 /κ → 0 is approached, we find that both τ 0 and τ π increase. As expected, τ π increases much faster than τ 0 and eventually diverges for g 0 /κ → 0, as the system gets trapped forever in the π-state. In contrast, τ 0 increases first when decreasing g 0 /κ, but then saturates at a finite level. Such a behaviour is unexpected based on the simple phase potential picture, where a fixed potential would imply diverging residence times for both states in the noiseless limit. The behaviour observed here hinges on the fact that the synchronization regime switches from "mixed" to "π" as one reduces the quantum parameter (i.e. reduces the quantum noise). The observations would change significantly for different parameters, where the system always stays in the mixed regime, for any value g 0 /κ > 0. Then, one expects the simple picture of a fixed phase potential to be approximately correct and both residence times to diverge as the noise is becoming weaker.

C. Noise-induced Synchronization Bistability
In the previous sections we have explained that some basic features of the classical-to-quantum crossover, like the increase of the residence times with decreasing quantum parameter, can be understood as effects of a decreasing quantum noise strength. The decrease of noise strength leads to less frequent transitions across an energy barrier in the effective potential. However, we made also less easily explained observations: (i) the reverse in the order of τ 0 and τ π as g 0 is being reduced, (ii) the saturation of τ 0 at low noise levels, (iii) the disappearance of stable 0-synchronization (for the applied parameters) in the classical limit g 0 /κ = 0. Whereas (i) and (ii) could originate from more complicated potential shapes with a combination of broad and narrow minima, (iii) suggests that the effective potential itself changes when the quantum parameter is varied. In Fig. 6 we show how the distribution of δφ evolves as a function of g 0 /κ. For very small values of g 0 /κ, there is only a single peak close to δφ ≈ π, in accordance with the single stable solution of the classical limit. While increasing the quantum parameter, this peak is first broadened. The increasing quantum noise strength allows the system to explore more of the effective potential around the minimum. A significant accumulation close to δφ ≈ 0 appears only for rather large quantum parameters, signaling the appearance of a second stable solution, i.e. a second minimum in the effective potential at δφ ≈ 0.
In addition to the above described appearance of a second stable solution, there are also parameter regions where already in the classical regime the effective potential has minima at both δφ = 0 and δφ = π (classical bistable synchronization). In this case quantum noise naturally drives transitions between the two synchronization states as soon as it is added to the description of the system. The number of observed transitions then naturally depends on both the noise strength as well as the potential shape.

VIII. OVERVIEW OF SYNCHRONIZATION REGIMES
In the previous sections we have shown examples of the different synchronization regimes in the presence of quantum noise and studied the properties of mixed synchronization in more detail. In the following, we map out the different synchronization regimes as a function of the mechanical coupling K and the quantum parameter g 0 /κ. Furthermore, we also discuss the case of detuned mechanical oscillators. Figure 7(a) shows the synchronization measure cos δφ for resonant oscillators. We have indicated the synchronization regimes (which are not sharply delineated). For the applied parameters, we find 0-synchronization for large mechanical coupling K in both the quantum and classical regime. In the classical, noiseless limits cos δφ → 1, indicating less fluctuations around the synchronization state. In contrast, at smaller mechanical coupling, we find more complicated behaviour: there is mixed synchronization for g 0 /κ ∼ 1, while the classical limit g 0 /κ → 0 selects either 0− or π−synchronization, depending on the mechanical coupling K. The "pixelated" region in Fig. 7(a) indicates that the system is multistable even in the classical limit. There, the residence times have become so large that the system is stuck in a random synchronization state depending on initial conditions and the transient behaviour.
Notably, the closer the system is to a border of synchronization regimes in the classical limit (this can be seen in Fig. 7(a) for small g 0 /κ when K is varied), the smaller the noise strength (g 0 /κ) that is needed to lead to mixed synchronization. As an example, in the "middle" of the classically π-synchronized regime (at about K/Ω 1 ≈ 0.1) similar mixed synchronization as compared to the system close to the regime border (K/Ω 1 0.15) appears only for larger values of the quantum parameter. An exception is of course the classically bistable regime, where mixed synchronization appears naturally as soon as there is noise. An interesting feature appears close to K/Ω ≈ 0.18, where the measure cos δφ shows a sharp dip in the middle of a 0-synchronized region. We suspect a non-linear resonance, since the oscillator trajectory x j (t) is no longer simply sinusoidal and period-doubling is observed. At the same time, the oscillation amplitude increases. For larger values of the mechanical coupling K, the trajectories are simply sinusoidal again, with the same frequency as for coupling strengths below the feature.
Up to now, we have studied the ideal case of identical mechanical oscillators. We now turn to the case where the mechanical oscillators have slightly different resonance frequencies, i.e. δΩ = Ω 2 − Ω 1 = 0. This is a typical situation in experiments, since fabrication inaccuracies lead to deviations between two systems. Figure 7(b) shows how δΩ affects synchronization. A finite δΩ corresponds to a tilted effective potential in the classical, noiseless limit, where a finite threshold for synchronization appears [31]. This classical threshold is indicated in Fig. 7(b) with a dashed line. Similar behaviour is visible for quantum parameters g 0 /κ > 0. However, for large g 0 /κ it is not possible to determine the onset of synchronization, because the measure cos δφ cannot distinguish between no synchronization and mixed synchronization. Instead, we now also have to analyze the trajectories in more detail, which reveal mixed synchronization for sufficiently large mechanical coupling. A significant deviation from the classical threshold cannot be observed at the given resolution. We expect the threshold to slightly increase for larger quantum parameter g 0 /κ. At the same time, however, the threshold is also smeared out due to quantum noise.
Here, we have chosen to show the dependence of the synchronization regimes on the mechanical coupling strength K and the quantum parameter g 0 /κ. Note that other parameters influence the synchronization type as well. In Fig. 2 we have already seen that the laser driving strength α L and the detuning ∆ influence the classical synchronization regimes. Also the mechanical damping Γ can affect the observed synchronization, cf. Fig. 3b and c. However, first of all it already influences the limit cycles of individual optomechanical systems by modifying the threshold to selfsustained oscillations. In addition, Γ has an influence on the mechanical noise strength. Note that, when changing these parameters, care has to be taken to remain on a stable limit cycle for each optomechanical system.

IX. THERMAL NOISE
So far, we assumed zero temperature environments for both the optical and the mechanical mode. In this section we investigate the effects of thermal mechanical noise on synchronization.
For our study we use the Langevin equations (5) with modified mechanical noise correlators to account for the coupling of the mechanical oscillators to a finite temperature bath, β jin (t)β * jin (t ) = β * jin (t)β jin (t ) = (n th + 1/2)δ(t − t ), where n th denotes the thermal occupancy of the bath. Hence, both quantum and thermal noise are included. For optical frequencies in the visible spectrum the effective thermal occupation of the optical bath is very small. Thus, the assumption of an optical bath at zero temperature is valid and the optical input noise terms are not modified.
In Fig. 7(c) we show an overview of the synchronization regimes as a function of thermal noise n th . Here, we chose a comparatively small quantum parameter g 0 /κ = 0.01 in order to observe mainly effects due to thermal noise. For small n th the results are similar to Fig. 7(a) at small g 0 /κ. In both cases the influence of the quantum or thermal noise is still weak. For increasing thermal noise strength n th , we find qualitatively the same behaviour as for increasing quantum noise. However, quantitative differences appear with increasing n th : Even though a mixed synchronization regime can develop for both quantum and thermal noise, the evolution of the relative weight of both synchronization states is different.
In the following, we want to estimate the critical thermal noise strength n * th at which thermal and quantum noise should have a comparable effect. At lower temperatures, quantum noise will dominate. The main source of quantum fluctuations is the laser shot noise (for the parameters explored here). Thus, we estimate the effect of optical quantum fluctuations on the mechanical oscillator, with the (symmetrized) shot-noise spectrum evaluated at the mechanical resonance frequency Ω [1,89], We expect similar effects from quantum and thermal noise if the shot-noise spectrum at the mechanical resonance frequency Ω becomes equal to the thermal force spectrum, S SN F F (Ω) = S th F F . The thermal noise spectrum at temperatures k B T Ω is where T is the temperature of the thermal mechanical bath. Setting S SN F F and S th F F equal, we find (in the resolved sideband regime κ Ω at ∆ ≈ Ω): where we used the optomechanical cooperativity [1] This approach suggests that observing quantum noise phenomena does not necessarily require a large g 0 /κ and very low temperatures. Instead, if the cooperativity is sufficiently large (comparable to values that enable ground-state cooling, C > n th ), quantum noise should dominate the behaviour of the system even in the presence of thermal noise.
However, depending on parameters, we find large deviations from this simple expectation. In these cases, the real shot noise spectrum is no longer well described by the weak-coupling expression of S SN F F (Ω) given above, and the actual noise strength may have a much larger value. Consequently, the transition between behaviour dominated by quantum noise vs. that dominated by thermal noise takes place at much larger values of n * th than those predicted by Eq. (9). In other words, it should be even easier to observe quantum noise in optomechanical synchronization than the naive ansatz would lead one to expect.
Here, we don't observe that quantum noise can be exactly mapped to thermal noise. This is already evident in the rescaled Langevin equations (6). Physically, the thermal force spectrum acting on the mechanical oscillator is flat in frequency, whereas the optical shot-noise spectrum is frequency dependent and is also modified by the dynamics of the system.

X. CONCLUSIONS
We have investigated the effects of quantum and thermal noise on two coupled optomechanical limit-cycle oscillators. One usually expects that noise prevents strict synchronization, i.e. exact phase locking and a sharp transition to synchronization. Here we have shown that fluctuations additionally drive transitions between 0-and π-synchronization, i.e. the two synchronization states that can appear in the absence of noise. We have discussed the residence times of these states and observed a smooth crossover between different synchronization regimes. Finally, we have compared the effects of quantum and thermal noise. We have argued that it should be possible to experimentally reach the regime where quantum noise dominates. This should happen when the optomechanical cooperativity is large enough for ground state cooling.
For further investigations it would be useful to identify a measure that can genuinely distinguish between an unsynchronized regime and 0-, π-and mixed synchronization. Finally, it will be very interesting to extend the insights obtained here to large optomechanical arrays.