Phase synchronization in coupled bistable oscillators

We introduce a simple model system to study synchronization theoretically in quantum oscillators that are not just in limit-cycle states, but rather display a more complex bistable dynamics. Our oscillator model is purely dissipative, with a two-photon gain balanced by single- and three-photon loss processes. When the gain rate is low, loss processes dominate and the oscillator has a very low photon occupation number. In contrast, for large gain rates, the oscillator is driven into a limit-cycle state where photon numbers can become large. The bistability emerges between these limiting cases with a region of coexistence of limit-cycle and low-occupation states. Although an individual oscillator has no preferred phase, when two of them are coupled together a relative phase preference is generated which can indicate synchronization of the dynamics. We find that the form and strength of the relative phase preference varies widely depending on the dynamical states of the oscillators. In the limit-cycle regime, the phase distribution is $\pi$-periodic with peaks at $0$ and $\pi$, whilst in the low-occupation regime $\pi$-periodic phase distributions can be produced with peaks at $\pi/2$ and $3\pi/2$. Tuning the coupled system between these two regimes reveals a region where the relative phase distribution has $\pi/2$-periodicity.


INTRODUCTION
The last few years has seen rapid progress in engineering and probing the properties of nonlinear oscillators in the quantum regime [1][2][3][4]. This has stimulated renewed theoretical interest in the properties of such systems, both at the level of individual oscillators and for more complex many-body realizations [5][6][7][8]. Of particular interest are phenomena, such as synchronization, which result from an interplay between nonlinearity and nonequilibrium features arising from a combination of gain and loss processes. In the classical regime, the standard paradigm for synchronization involves a nonlinear oscillator in a limit-cycle state which has a well-defined amplitude, but no preferred phase [9,10]. Such oscillators have a tendency to adjust their rhythm to match either a weak external drive or that of other oscillators to which they may be coupled, typically leading to the emergence of a definite phase (or relative phase) for the oscillations. Although synchronization has been studied for a very long time in classical oscillators [9], the systematic study of this behavior in quantum oscillators outside the regime where semiclassical approximations work well [11] is quite recent [12][13][14][15][16][17].
Significant efforts have also been devoted to proposing ways in which the behavior could be probed in experiment using systems such as trapped ions [16,22] or superconducting circuits [32].
Although studies of synchronization in the quantum regime have employed a wide range of model systems, they have focussed (with occasional exceptions [33]) on systems whose dynamics is essentially just a limit-cycle. Here we instead explore how weak coupling generates synchronization, in the form of particular phase preferences, in a quantum oscillator which has a more complex bistable dynamics. We do this by proposing a minimal model for a quantum oscillator that displays a limit-cycle state as well as a low occupation-number state (in which the oscillator simply fluctuates about the origin) and can be tuned to a bistable regime in which both of these states coexist.
Our model involves only dissipative processes in which photons (quanta) are lost or gained, ensuring that an isolated oscillator never has a preferred phase. The key ingredient of the model responsible for generating bistability is a two-photon gain process. This is balanced by two channels of photon loss in which either a single photon or three photons are annihilated in the oscillator. Different dynamical states of the oscillator (low occupationnumber regime, limit-cycle and bistability) are achieved by tuning the relative sizes of the gain and loss rates. In many ways the model is a logical extension of the muchstudied QvdP oscillator which combines one-photon gain and two-photon loss [16]. However, the QvdP model always displays a limit-cycle, albeit one whose size depends on the ratio of loss and gain rates.
We investigate in detail the phase synchronization that occurs when two of the bistable model oscillators are coupled via a weak photon exchange process. This leads to a rich range of behavior in the relative phase distribution, with a different pattern of phase preferences emerging depending on the underlying dynamical states of the oscillators. When the gain is strong and the oscillators are in limit-cycle states, maximal values of the relative phase distribution form at 0 and π, matching what is usually found for such states(as seen, e.g. in the QvdP model [16]). For oscillators with low occupation numbers and no limit-cycle, however, relative phases of π/2 and 3π/2 are preferred instead, a result which we argue can be understood as a result of the two-photon gain. All four peaks emerge simultaneously for a small, intermediate parameter regime. Furthermore, this π/2-periodic behavior is strongest when the bistability is most pronounced. Interestingly, and in strong contrast to the QvdP oscillator, we find that the strength of the synchronization in the limit-cycle state does not increase with increasing photon numbers.
The outline of the rest of this paper is as follows. In Sec. II, we start by introducing our bistable quantum oscillator model and exploring its steady-state and dynamical properties. Then we investigate the behavior of two coherently coupled oscillators in Sec. III, focussing in particular on the way in which the pattern and strength of features in the relative phase distribution of the system depend on the underlying dynamical properties of each oscillator. We conclude in Sec. IV and the Appendixes provide details about aspects of the calculations employed.

A. 321-Oscillator Model
Our oscillator model involves three dissipative processes, as illustrated in Fig. 1(a). A two-photon gain process with rate κ 2 drives the oscillator to higher photon numbers, whilst a one-photon loss process damps it at a rate κ 1 ; an additional three-photon loss process at rate κ 3 is included to stabilise the system, ensuring that it has a steady state for any strength of the gain. The master equation for a single oscillator in the interaction picture is given by [16,34,35], whereâ is the oscillator lowering operator and we have defined D[Ĉ] (ρ) =ĈρĈ † − 1 2 {Ĉ †Ĉ , ρ}. Our model makes an interesting contrast with the QvdP oscillator [16,17], where one-photon gain is balanced by two-photon loss. The presence of a non-linear gain process in our model leads to important features, such as bistability, not seen for the QvdP.
The steady-state properties are readily found by exploiting the fact that the system is purely dissipative, so that the dynamics of the diagonal and off-diagonal matrix elements of the density operator in the number (Fock) basis are decoupled [36,37]. The master equation can therefore be rewritten as a set of k equationṡ where ρ (k) n = n|ρ|n + k , with |n the n-th number state, and M (k) a matrix. For the diagonal elements, writing out Eq. (1) explicitly leads to the coupled set of equationṡ for the probabilities, P n = n|ρ|n , with G n = [κ 1 n + κ 2 (n + 1) (n + 2) +κ 3 n (n − 1) (n − 2)] , A n+1 =κ 1 (n + 1) , B n−2 =κ 2 n (n − 1) , C n+3 =κ 3 (n + 1) (n + 2) (n + 3) , from which the form of M (0) follows. In the steady state, the off-diagonal terms ρ (k =0) all go to zero and the eigenvector of M (0) with zero eigenvalue gives the P n distribution.

B. Phase diagram
The steady-state of the oscillator can be characterised by the behavior of the P n distribution along with the Wigner distribution [35], W(α r , α i ). Figures 1(b)-(e) show how the state of the system evolves as the ratio κ 1 /κ 2 is changed for a small (fixed) value of κ 3 /κ 2 . When the non-linear gain dominates (κ 1 /κ 2 1) the oscillator is driven to large phonon numbers with an almost Gaussian P n distribution centered at a value n = n nP n 1 [see Fig. 1(b)]. The corresponding Wigner distribution reveals a distribution with a ring of maxima [ Fig. 1(c)], we classify this as a limit-cycle (LC) state, as it has a well-defined average amplitude, but no preferred phase [16]. In the opposite limit of dominant loss (κ 1 /κ 2 1), the oscillator is damped to the lowest photon number states, leading to a sharp peak in the P n distribution at n = 0. In this regime the Wigner distribution displays a single maximum at the origin [ Fig. 1(e)] and we call this a fixed point (FP) state. In between these limits we find bistability (B) where features from both the LC and FP states can be found in the Wigner distribution [see Fig. 1(d)] and two peaks of similar size feature in the P n distribution [38].
The bimodality of the P n distribution is captured by a sharp peak in the second moment [39] µ (2) = n 2 − n 2 , as shown in Fig. 1(b). The distributions with the highest values of µ (2) are found to be those with the most pronounced bistability, i.e. with two peaks of comparable area that are separated by a significant gap.
We use a standard mean-field (or semiclassical) approach to understand the origin of the oscillator's dynamical states (the details are described in Appendix A). Two physically relevant mean-field solutions for the average photon number are found: a zero photon solution n 0 = 0 (corresponding to the fixed point), and a nonzero solution n + = κ2 3κ3 1 + 1 + 3κ3 κ 2 2 (4κ 2 − κ 1 ) (corresponding to the limit-cycle). Linear stability analysis reveals that n 0 is stable for κ 1 > 4κ 2 , whilst n + is both non-negative and stable for 3κ 1 κ 3 < 12κ 3 κ 2 +κ 2 2 . Hence the mean-field approach predicts a region of bistability associated with the coexistence of these stable solutions in the parameter regime These predictions are compared with the behavior of the P n distribution in Fig. 1(b). The n + solution matches the average photon number n very well deep inside the limit-cycle regime, where photon numbers are large. However, the mean-field calculation doesn't describe the extent of the bistable region very accurately, which is not surprising given the strong quantum non-linearity in this regime.
The behavior of the oscillator as a function of both κ 1 /κ 2 and κ 3 /κ 2 is summarised in Fig. 2. Classifying the oscillator state as either FP, B, or LC, based on the corresponding Wigner function leads to the 'phase diagram' shown in Fig. 2(a) whilst the behavior of the average occupation number is shown in Fig. 2(b). A sharp transition occurs between the FP and LC states via a bistable region which is narrow and very well defined (e.g. via a clear peak in the second moment) for κ 3 /κ 2 1 as κ 1 /κ 2 is reduced. Average photon numbers tend to decrease as κ 3 /κ 2 is increased, and although we can still use the Wigner distribution to distingusih the different states, the differences between them become much less distinct.

C. Dynamical Properties
We now turn to the dynamical properties of the system to understand whether the bistable states we have identified are also metastable in the sense that they display slow switching between the two states. To investigate this, we start by calculating how the eigenvalues with largest real parts (i.e. least negative, but non-zero) of the system behave [40][41][42].
The largest eigenvalue of the matrix M (k) in Eq. (2) can be used to obtain the slowest timescale associated with the dynamics of ρ (k) , Fig. 3(a). The k = 0 case, τ 0 , describes the relaxation of the diagonal elements and it becomes very large for a range of κ 1 /κ 2 . The other timesscales τ 1,2 , describe the relaxation of phase preferences in the system. Although they never become as large as the peak values of τ 0 and display no obvious signature of the bistability, they do change significantly as the system evolves from FP to LC states, becoming orders of magnitude larger in the latter case.
We can use the emergence of a single very slow timescale as a criterion to identify metastability in the system [42]. We calculate the ratio of the differences between the three largest eigenvalues as shown in Fig. 3(b), and take (somewhat arbitrarily) a separation of λ 1 and λ 0 an order of magnitude larger than λ 2 and λ 1 , i.e. M ≥ 10, as a threshold for metastability. The parameter range identified via this criterion is, unsurprisingly, a small subset of that labelled on the basis of the steady-state as bistable in Fig. 2(a). Nevertheless, the peaks in M and the second moment plotted in Fig. 2(a) match up well. We can get more insight into the dynamical properties of the system by looking at quantum jump trajectories [43], obtained by unravelling the master equation.
The system is evolved in time with the non-Hermitian Hamiltonian For each short time step, quantum jumps of three different kinds (one-photon loss, two-photon gain, or three-photon loss process) can occur with a probability that depends on the state of the system (e.g. the probability of two photon gain occurring over the interval δt is given by The frequency of the different jump processes is illustrated in Fig. 4 for sample trajectories obtained for parameters corresponding to the three different states of the system (LC, FP and B). Within the LC state the oscillator has a large number of photons and consequently displays a high level of activity (i.e. frequent jumps). In contrast, all jump processes are strongly suppressed (and the three-photon loss especially so) in the FP state because of the very low occupation numbers in this regime. Within the region which is both bistable and metastable (based on the behavior of the eigenvalues), the oscillator switches between periods in which it exhibits high and low levels of activity; the switching continues indefinitely, never settling into one state or the other. This intermittency in the dynamics of the trajectories is consistent with our interpretation of this regime as metastable [42].

III. SYNCHRONIZATION OF COUPLED OSCILLATORS
We now explore how phase ordering and synchronization occurs when two of these oscillators are coupled together weakly. For simplicity, we consider two identical oscillators and assume a coherent (photon-exchange) interaction of the form [16,17,22] whereâ j is a lowering operator for oscillator j and J is the strength of the coupling. The master equation of the coupled system is given bẏ where the dissipation terms follow from Eq. (1): j ](ρ). Classical limit-cycle oscillators typically synchronize when they are coupled together weakly, developing a preference for one or more relative phase values [9,44]. In the quantum regime, the phase states |ϕ = n e iϕ |n /2π, can be used to construct a relative phase probability dis- 2 /J 2 ) as a function of κ1/κ2 (calculated with perturbation theory), with J/κ2 = 10 −2 and κ3 = κ2 × 10 −1 . The range of κ1/κ2 shown spans the three states (FP, B, and LC). Very weak single-photon loss (κ1 κ2) leads to coupled limit-cycles with peaks at 0 and π. As the single-photon loss rate is increased, this π-periodic pattern vanishes before then reappearing with peaks at π/2 and 3π/2. Eventually, very strong single-photon loss (κ1 κ2) again suppresses the pattern. (b) The dominant Fourier coefficients F2 (red) and F4 (blue) [defined in Eq. (9)] (here units are chosen such that κ2 = 1). F2 accounts for the π-periodic component of the relative phase distribution and its sign determines the position of the peaks; F4 is the next largest, though it is much less important than F2, except for the region in which F2 passes through zero. The corresponding π 2 -periodic P (φ) distribution for the case where F2 is zero is shown in the inset; the peak-to-peak height is 10 −4 in units of (J/κ2) 4 .
tribution [19,20,22,[45][46][47] n + k, m|ρ|n, m + k , where φ = ϕ 1 − ϕ 2 is the relative phase of the two oscillators. If the two oscillators are uncoupled, their phases are independent and the relative phase distribution is uniform, P (φ) = 1 2π . In the following we use the relative phase distribution to explore the impact of coupling on the behavior of our oscillator model in each of the dynamical states it displays.
Since we are interested in the effect of a very weak coupling between the oscillators, we use a perturbation method [18,19,28] to calculate the relative phase distribution, details of which are described in Appendix B. Figure 5(a) shows how the relative phase distribution evolves as a function of κ 1 /κ 2 and the system passes from FP to LC via the bistable region. Interestingly, the system displays π-periodic phase distributions with peaks which reach a similar size in both the LC and FP regimes. However, the location of the peaks are different in the two cases and in between (the bistable region) the phase distribution appears to flatten. A more detailed picture of the phase behavior is obtained by looking at the Fourier coefficients of the relative phase distribution, F k = Re[ ∞ n,m=0 n + k, m|ρ|n, m + k ]. For our coupled oscillator system, these coefficients are non-zero for even k and typically get smaller very rapidly with increasing k. The two most important coefficients F 2 and F 4 are shown in Fig. 5(b) (calculated to second and fourth order in J respectively). The πperiodicity of the distribution apparent in Fig. 5(a) stems from the fact that the magnitude of F 2 is almost everywhere much larger than that of F 4 . However, as the system transitions from LC to FP (via the bistablity) F 2 changes sign to produce the shift in the locations of the peaks; as F 2 passes through zero F 4 dominates, giving rise to an unusual π/2 periodic distribution.
The different π-periodic patterns that arise in the phase distribution can be understood using simple arguments that exploit the specific characteristics of the system in the LC and FP states. Well-within the LC regime photon occupation numbers are large and semiclassical approaches work well as Fig. 1(b) illustrates for the single oscillator case. A straightforward calculation described in Appendix C recovers the preference for relative phases of 0 and π which is generic for coherently coupled limit-cycle oscillators [16,19].
The pattern of peaks in the relative phase distribution at π/2 and 3π/2 seen for the FP can be understood by focussing on the limit where κ 1 /κ 2 1. In this regime photon occupation numbers are very small, the value of κ 3 becomes irrelevant, and we can simplify the perturbation theory calculation by assuming only the three lowest Fock states of the oscillators have non-negligible occupations (see Appendix B for details). In this limit we find where f (κ 1 , κ 2 ) is a (positive-valued) function of the two rates, whilst P 0 , P 1 and P 2 are the occupation probabilities of the three lowest Fock states for the corresponding uncoupled oscillator. Clearly the behavior is always π-periodic, but the positions of the maxima depend on the occupation probabilities which in turn depend on the details of how the system is driven and damped in this regime. In our oscillator model, the two-photon driving gives a boost to P 2 , generating a steady-state with a photon distribution where P 2 P 0 /P 2 1 > 1 which leads to the peaks at π/2 and 3π/2 seen in Fig. 5. Thus the fact that our oscillator model displays a π-periodic relative phase distribution with peaks at π/2 and 3π/2 can be seen as being due to the particular form of the lowest order nonlinearity in the system, two-photon gain, which shapes the steady-state number distribution.
Applying a very similar analysis to oscillators with different gain and loss processes casts light on the different ways relative phase preferences can develop (within the same low-occupation number limit). In particular, for an oscillator coupled to a thermal bath [with onephoton gain and loss related by a ratio of rates n/(n+1), with thermal occupation number n] and the QvdP oscillator [16] (one-photon gain and two-photon loss) one obtains expressions like that in Eq. (10). In each case f takes a form that is different in detail, but the sign of P 2 1 − P 0 P 2 still determines the location of the peaks. For the QvdP oscillator in the strongly damped limit P 2 P 0 < P 2 1 , hence one finds the same phase behavior as in the large photon number regime [16] (peaks at 0 and π). In contrast, for a thermal oscillator P 2 P 0 = P 2 1 and so no phase preference is expected.
Finally, we look at how the phase behavior of the system varies over the full phase diagram. Figure 6 shows how the components F 2 and F 4 behave for a range of κ 3 /κ 2 as well as κ 1 /κ 2 . The transition between negative and positive values of F 2 is sharpest, and the corresponding peak in F 4 strongest, in the region where κ 3 /κ 2 1, which is where the system also displays metastability. Figure 6 also reveals a somewhat surprising feature of the phase behavior. The strongest phase preferences are not associated with the limit-cycle regime. Furthermore, deep within the limit-cycle regime the magnitude of F 2 starts to saturate, even as the average photon number continues to grow, i.e. as κ 3 /κ 2 is further reduced in the bottom left quadrant of Fig. 6(a). This is in contrast to other quantum limit-cycle oscillators, such as the QvdP oscillator, for which phase synchronization effects are enhanced by increasing the photon number [16]. We can gain some insight into why this occurs from the phase dynamics of an uncoupled oscillator. At least in the semiclassical regime of large photon numbers, we expect to see stronger synchronization effects emerge when coupling is introduced in systems where the phase diffusion is weaker [9,19].
Using an approximate analytic approach (see Appendix D for details), we find that the phase diffusion rate is simply proportional to κ 2 in the regime where gain dominates over losses and photon numbers are large. This contrasts with the behavior of the QvdP oscillator where the phase diffusion is ∝ 1/ n , which is why for this system synchronization effects get stronger as photon numbers are increased.

IV. CONCLUSIONS
We have introduced a simple oscillator model with a two-photon gain process balanced by one-and threephoton losses that can be used to engineer a bistable oscillator state. The bistability occurs when the gain/loss rates are tuned between regimes where the oscillator displays a limit-cycle state characterised by a large amplitude (but no preferred phase) and a fixed-point state where occupation numbers are low. Quantum trajectory simulations show clear evidence of metastability in the bistable region, signalled by intermittency in the frequency of the quantum jumps.
When two such oscillators are coupled together weakly their relative phase distribution displays a rich pattern of behavior. Whilst the system has the usual predominantly π-periodic distribution with peaks at 0 and π within the limit-cycle regime, in the limit of low occupation numbers the peaks appear at π/2 and 3π/2. In between these two regimes, where the bistability arises, the distribution can instead be π/2-periodic. We identify the presence of twophoton gain as the key factor giving rise to this unusual behavior.
Our goal in this work has been to explore the properties of the relative phase distribution and its connection to the underlying oscillator dynamics in the simplest possible model system displaying limit-cycles and bistability. We plan to address the question of how the complex behavior we uncovered could be realised, and detected, in future work. It would also be interesting to investigate how the dynamics is affected by strong couplings [48] and the patterns of phase behavior that arise in systems where more than two oscillators are coupled together [14,16,49]. Â B . We choose to carry out this approximation after first normal ordering the operators.
Making the substitution a = re iϕ where r and ϕ are classical amplitude and phase variables, allows us to rewrite Eq. (A1) aṡ Evaluating the real and imaginary parts separately, leads toφ = 0 anḋ with steady state solutions n 0 = 0 and These solutions correspond to the average photon number and so are subject to the constraints of being both positive and real. The stability of the solutions can be determined, e.g. through the properties of the relevant Jacobian.
The negative branch, n − , is never both physical and stable, the zero photon solution, n 0 , is stable for κ1 κ2 > 4, and the positive branch, n + , is stable whenever it is physical, i.e. for κ1 κ2 < 4 + κ2 3κ3 . This mean-field calculation results in two stable solutions for the average photon number and therefore a predicted bistability in the photon number for the parameter regime Appendix B: Perturbation theory

General Method
Perturbation theory provides a convenient way of calculating the way in which the relative phase distribution behaves for weak coupling. The steady-state of the uncoupled (J = 0), two-oscillator system [Eq. (7)] only has diagonal terms. Treating the coupling as a perturbation [20,28] allows us to calculate the terms in the first off-diagonal as a function of the uncoupled oscillator terms. Each subsequent off-diagonal can, in turn, be calculated from previous ones [19].
In the steady-state, this reduces to sets of simultaneous equations with the coupling term, ∆ n,m coupling together terms with different p-values. The zeroth-order terms are the diagonal (p = 0) elements, the uncoupled probabilities ρ (0) n,m = P n P m , which necessarily sum to unity. The first-order terms are obtained by substituting the zeroth-order terms into the expression for ∆ (p) n,m , leading to non-zero contributions for p = 1 and the process is continued to higher order in J recursively.
The first-order terms obey the relation ρ  and hence sum to zero [19], which means that they make no contribution to the relative phase distribution [Eq. (9)] since it depends on sums of the off-diagonal elements. The sum of the p = 2 terms, however, is real and finite and so does contribute resulting in a π-periodic relative phase distribution. Continuing to higher orders, we find that all of the odd-p terms sum to zero, and so only the even p sums contribute to the relative phase distribution. In particular, the p = 4 terms, lead to a π/2-periodic contribution which can dominate the phase distribution when the π-periodic terms vanish.

Low Occupation-Number Regime
This calculation can be simplified and solved analytically in the limit of very low photon numbers. We proceed by assuming only the lowest three photon states are appreciably occupied, i.e. P n>2 = 0, and hence truncate the state space to include only |0 , |1 , and |2 . Due to the size of the Hilbert space, only a single term contributes to the relative phase distribution, 0,0 . In the steady-state, Eq. (B1) with p = 2 leads to using the relation ρ This results in the relative phase distribution This is a π-periodic distribution and the position of the peaks is determined by the steady state of the uncoupled oscillators. The two-photon driving in our model ensures P 2 P 0 > P 2 1 , which leads to peaks at π/2 and 3π/2.
In this Appendix we return to the case of a single oscillator and obtain an estimate for the phase diffusion rate in the semiclassical limit where photon numbers are large. In the semiclassical regime at least, the strength of phase synchronization in coupled oscillators is determined by competition between the coupling and the rate of phase diffusion in the individual oscillators, with slower phase diffusion leading to stronger phase preferences [9,19] The phase distribution for a single oscillator takes the form [19,45,50] P (ϕ) = 1 2π ∞ n,m=0 n|ρ|m e i(m−n)ϕ , with Φ (k) = ∞ n=0 ρ (k) n , where ρ (k) n = n|ρ|n + k . Although the behavior is in general quite complex, we can obtain a simple approximate description in the semiclassical limit where the density matrix is tightly peaked around a large average photon occupation number [19,36].
Using (1) [and the notation introduced in (B1)], we can obtain the equation of motion for Φ (k) Φ (k) = κ 2 n −G (k) n + A (k) n + B (k) n + C (k) n ρ (k) n . (D3) In the semiclassical limit, i.e. the strong gain regime where γ = κ 1 /κ 2 1 and Γ = κ 3 /κ 2 1, the photon number saturates to a large value n 2κ 2 /(3κ 3 ) [see Eq. (A4)]. We proceed by assuming we can replace n by n , and expand the square-roots appearing in A (k) n , B (k) n , etc, treating 1/ n together with γ and Γ as small quantities [19,36,51]. This leads to the simplified equa-tionΦ Hence to leading order the relaxation timescale for the k-th component, Φ (k) , is simply proportional to 1/κ 2 , τ LC k κ 2 4/(5k 2 ). Notice that the semiclassical approximation here assumes all the off-diagonal elements (ρ (k) n ) decay at the same rate, hence the decay time for k = 1 is also an approximate linewidth for the oscillator [36]. The slowest timescales τ k associated with the matrices M (k) in Eq.
Finally, using the definition Eq. D2 and Eq. D4, we see that the phase distribution obeys a diffusion equation [19,51] This is very different to what is found obtained in a similar calculation for the QvdP oscillator (or indeed the laser [36]) where the diffusion constant is ∝ 1/ n . Hence for such systems phase diffusion gets weaker (and the linewidth narrower), so that synchronization effects get stronger, as the photon number increases [16].