Dynamics of Many-Body Quantum Synchronisation

We analyse the properties of the synchronisation transition in a many-body system consisting of quantum van der Pol oscillators with all-to-all coupling using a self-consistent mean-field method. We find that the synchronised state, which the system can access for oscillator couplings above a critical value, is characterised not just by a lower phase uncertainty than the corresponding unsynchronised state, but also a higher number uncertainty. Just below the critical coupling the system can evolve to the unsynchronised steady state via a long-lived transient synchronised state. We investigate the way in which this transient state eventually decays and show that the critical scaling of its lifetime is consistent with a simple classical model.

Several different ways of quantifying the synchronisation of quantum oscillators have been proposed [11,13,20,21,24] and the connection between synchronisation and entanglement [11,14,16,20,24,25,30] has been examined in a variety of different systems. Comparisons of synchronisation in quantum and semiclassical oscillator models have revealed significant quantitative and qualitative differences in behaviour [13,[25][26][27]. Detailed proposals have also been made for experiments which could probe synchronisation in the quantum regime using trapped ions [13,21], optomechanical systems [12] or superconducting circuits [31].
Synchronisation in quantum models of many coupled limit-cycle oscillators [12,13,23,32] has so far received rather less attention than few-oscillator systems. However, quantum many-oscillator models form a novel class of many-body system and the synchronisation transition they undergo makes an interesting comparison not just with classical or semiclassical oscillator systems, but also with the rich variety of nonequilibrium transitions which have been studied extensively in other types of many-body quantum systems [33][34][35][36], including e.g., the driven dissipative Bose-Hubbard model [37][38][39].
A particularly simple model system for studying many-body synchronisation consisting of coupled quantum van der Pol oscillators was introduced by Lee and Sadeghpour [13]. In their work Lee and Sadeghpour compared the predictions of quantum and semiclassical versions of the model, and found that the transition which the oscillators undergo between unsynchronised and synchronised states consistently occurs at a lower value of the inter-oscillator coupling strength in the quantum case. Here we examine the properties of the synchronised and unsynchronised states of this system together with its critical dynamics. We find that the phase-ordering which occurs when the oscillators synchronise is accompanied by other changes in the state of the system: a decrease in the average occupation number and an increase in number uncertainty. For coupling strengths just below the transition, the system displays critical slowing down with a long-lived transient synchronised state emerging and then eventually decaying. We look at exactly how the long-lived transient state decays and show that the critical scaling fits a simple classical model. The rest of this paper is organised as follows. We begin by introducing the many van der Pol-oscillator model in section 2 and then go on to describe the synchronisation transition it undergoes in section 3. We compare the quantum properties of the synchronised and unsychronised states of the system in section 4. We explore the dynamics close to the transition in section 5 and then analyse the critical scaling in section 6. Finally, we present our conclusions in section 7.

Model oscillator system
The van der Pol (vdP) oscillator is a simple limit-cycle oscillator and is a popular choice of model to study synchronisation [1,2,40]. In the quantum regime the vdP model is described by a harmonic oscillator which gains individual photons at a rate κ 1 (linear anti-damping) whilst also losing two photons at a time with rate κ 2 (nonlinear damping) [13]. The master equation in the interaction picture takes the form, where a is the oscillator lowering operator. The ratio R=κ 2 /κ 1 controls the size of the limit cycle in the system [13,26,41]. For extremely small values of R the average photon number becomes very large, growing as 1/R, and semiclassical methods [13] provide an accurate description. However, for larger R values the behaviour of quantum and semiclassical versions of the model become markedly different [13] and eventually, in the limit R  ¥ the steady state of the system becomes entirely confined to just the two lowest number states.
We are interested in the behaviour of an ensemble of N identical vdP oscillators with all-to-all coherent couplings of strength ε, for which the Hamiltonian is As is typically the case with quantum many-body problems, the extremely large state space involved precludes an exact numerical treatment and so we follow Lee and Sadeghpour [13] in assuming N is large and adopting an approximate self-consistent mean-field approach [12,13,37,38,39,42] The master equation is solved self-consistently via numerical integration, starting from a chosen state, with the value of a á ñ updated at each time step [12]. In each case we carried out integrations of the master equation using a Runge-Kutta algorithm and generally chose a coherent state as an initial condition. We worked in the number state basis, using a cut-off, nmax, chosen to be large enough not to influence the results.
The self-consistent mean-field approach is very commonly used in studies of non-equilibrium quantum many-body systems [33][34][35][36][37][38][39], often as an approximate description for a lattice of optical or microwave cavities in which individual cavities are coupled to a small number of their nearest neighbours. In such situations meanfield calculations provide a useful starting point though they are not expected to describe the behaviour faithfully in low-dimensional systems [33,43,44]. In this case we have in mind a large number of individual vdP oscillators with all-to-all coupling and so expect that the mean-field approach will work increasingly well as that number is increased. Systems of many nonlinear quantum oscillators with all-to-all couplings are of interest in other contexts and detailed proposals have been made for schemes which could realise such systems [45].
Note that even after we have assumed a large number of coupled vdP oscillators, the value of the damping rate ratio R can still be varied. This parameter controls the size of the limit-cycles of the individual oscillators with R 0  setting the semiclassical limit for an uncoupled oscillator [46].

Synchronisation transition
The long-time state of the many-body vdP system given by equation (4) displays either synchronised or unsynchronised behaviour depending on the coupling strength ε, the ratio of rates of the 2-and 1-photon processes, R, and the initial state of the system. The synchronised state is characterised by a clear phase preference, signalled by a non-zero value of [12,13] a á ñ, which oscillates periodically in time with a magnitude that settles down to a constant value. In contrast, the unsynchronised state has no preferred phase so that a 0 á ñ = , leading to a time-independent state which matches the steady state of the corresponding uncoupled vdP oscillator. The value of a á ñ | |therefore provides a natural order parameter for the system. Initial states always exist which allow the system to reach an unsynchronised state 1 , but the synchronised state can only be accessed for couplings beyond a certain critical value ε c which depends on R. The behaviour of ε c as a function of R was mapped out by Lee and Sadeghpour [13]. The value of ε c tends towards zero for small R, but grows rapidly with increasing R, before apparently diverging at a finite value of R. The behaviour of the vdP oscillator becomes more strongly quantum mechanical as R is increased and this is reflected in the value of the critical coupling [13] which always takes a lower value in the quantum model compared to its semiclassical counterpart with the difference between the two growing rapidly with R. Figure 1 shows examples of the dynamics of a á ñ | |and a á ñ as a function of ε with everything else kept fixed. For ε>ε c , a á ñ rapidly reaches a periodically oscillating state whose amplitude and period depend on ε. For the smallest values of ε the value of a á ñ | |decays exponentially, but for values of ε just below ε c the character of the decay is very different and it instead occurs in two distinct stages: a slow part followed by a very rapid part. The development of a very slow relaxation time in the dynamics is a precursor to the emergence of the synchronised state at ε=ε c and we investigate its properties in detail in sections 5 and 6 below. However, this is not the only interesting feature in the dynamics: for fairly weak couplings there is a regime in which the decay of the order parameter actually speeds up with increasing couplings (see figure 1(a)).
When the coupling exceeds the critical value, ε>ε c , the behaviour of the system in the limit of long times is determined by the initial state of the system. Using an initial coherent state with an eigenvalue, init a , which is varied, we find a transition from unsynchronised to synchronised final states at a critical value, α crit , as shown in figure 2(a). The critical value of α init decreases with increasing ε, mapping out a separatrix between initial conditions that lead to synchronised and unsynchronised states as shown in figures 2(b) and (c). We note that although a dependence of the long-time density operator on initial conditions is not usually expected for open quantum systems, such behaviour does emerge when they are treated approximately, using e.g. self-consistent mean-field methods, as we do here [35,39,43].
Much of the behaviour of the order parameter seen in figure 2 can be captured using a simple classical effective-potential model [35], V(r), such as that sketched in figure 3. In such a model the dynamics of the order parameter is given by and hence as soon as the gradient in the potential becomes zero r will stop changing, giving rise to a fixed point in the dynamics. Although we don't have a way of deriving the form of V(r), its basic properties are clear. It will always have a stable fixed-point solution at r=0 (the unsynchronised state), and above a critical value of ε a second stable c e e >  , signifying a synchronised state. We adopt units of time such that κ 1 =1 throughout. 1 For example, if instead of a coherent state, we chose an initial Fock (i.e. number) state, for which a 0 á ñ = , the coupling term will start off being zero. Without this term, off-diagonal terms in the number state basis, which are all zero initially, remain zero and hence a á ñ will be zero for all later times within the mean-field description, so that the system will always evolve to the unsynchronised state.
fixed-point solution with r>0 should emerge (see figure 3(a)), reflecting the coexistence of unsynchronised and synchronised (r>0) states. Furthermore, a third, unstable fixed point in the form of a peak between the two minima in the potential must emerge at the same time as the r>0 fixed point: the effective potential model undergoes a saddle-node bifurcation [47] at a critical coupling. The unstable fixed point corresponds to the separatrix between synchronised and unsynchronised states seen in figure 2. Figure 3(b) shows the potential close to, but below, the critical value of ε. Here there is only one fixed point, at r=0, but the proximity to the critical point means that the gradient of the potential becomes very shallow, and correspondingly a bottleneck appears [47] in the dynamics which means the time taken for r to decay becomes extremely long. As figure 3(c) illustrates, for small enough ε the long timescale decay is expected to disappear. One example of a potential which incorporates these features is where the functions of R, f and g, always take positive values. This kind of effective potential model predicts the phase structure of the system by design, but it also predicts more subtle features of the system such as the scaling behaviour of the slow decay time that emerges for couplings below the critical value. We look in detail at the scaling behaviour of the relaxation time in the vdP system in section 6 where we compare it to the prediction of the classical bifurcation model.

Synchronised versus unsynchronised states
We now turn to the properties of the synchronised and unsynchronised states which emerge in the limit of long times, beyond the value of the order parameter r a = á ñ | | that we have focussed on so far. The number distribution, P n n n r = á ñ ( ) | | (the diagonal elements of the density operator in the number state basis), shown in figure 4(a) reveals clear differences between the synchronised and unsynchronised states which emerge in the limit of long times. The average occupation number, n á ñ, is reduced in the synchronised state compared to the unsynchronised case, something which is not surprising as a reduction in the oscillation amplitude is a common accompaniment to synchronisation when classical limit-cycle oscillators are coupled together [2]. However, the P(n) distribution is clearly much broader when the oscillators are synchronised and has a much larger vacuum state occupation probability. This suggests an interplay between phase and number fluctuations at the transition which is something we might expect intuitively. Figure 4(a) also illustrates the way in which larger values of R push the P(n) distribution down to lower occupation numbers.
To quantify how the number and phase properties of the quantum state of the vdP system change at the transition we need to adopt suitable measures of the uncertainty in both of these quantities. There are a number of different ways of defining phase uncertainty [48], starting from the phase distribution   (12)]. We know already that the synchronisation order parameter a á ñ | |is zero until the critical value of ε, at which point it takes a finite value (for For ε<ε c the system is always in the unsynchronised state (which is independent of ε); for ε>ε c the system reaches a synchronised state whose properties are weakly dependent on ε.
an appropriate choice of initial conditions) indicating (partial) phase synchronisation. The phase uncertainty naturally drops at the transition point: it is maximal for the unsynchronised state ln 2 f p D = ( ) and must be lower when a phase preference emerges. The conjugate variable, the average occupation number n á ñ, drops significantly at the transition and the uncertainty in n goes up by a small amount. That a rise in number uncertainty should accompany a drop in phase uncertainty is perhaps not surprising, given their conjugate relationship, but the combined uncertainty in this case always remains well above the lower bound for the entropic uncertainties, ln 2p ( ). What is interesting here is that the phase-ordering that occurs at the transition is accompanied by number-disordering even when the latter is not required to ensure that the uncertainty relation is satisfied.

Dynamics of metastable decay
In this section we look at the dynamics of the system just below the critical coupling where the system eventually evolves to a single unsynchronised state, but the time it takes to reach that state can become extremely large. We focus in particular on the question of how exactly the system makes its transition to the final unsychronised state. Then in the next section we will analyse the way in which the lifetime of the transient state scales with the distance from the critical coupling.
For couplings just below ε c , the order parameter exhibits a long period of very slow decay followed by an abrupt drop in its value. The final drop in r is accompanied by rapid changes in the number distribution of the oscillator, which are clearly seen in the behaviour of n mp , the n value corresponding to the peak in the P(n) distribution, as shown in figure 5(a). In contrast to the order parameter, the evolution of n mp is not monotonic. Instead it displays a kind of latching behaviour: as the value of r drops abruptly n mp first dips before rising again to a value that is larger than its initial value. Figure 5(a) also reveals an interesting dependence on the value of R. As expected from what we saw of the long-time steady states of the system (e.g. figure 4), for smaller limit cycles (i.e. larger values of R) where the quantum noise is strong (and semiclassical descriptions provide a less accurate description of the system [13]), the slow relaxation process involves a substantial readjustment of both the number and phase properties of the system's state. Indeed, for smaller limit-cycles the value of n mp goes through very large variations-actually dropping to zero before growing again to a value larger than before. In contrast, for larger limit-cycles (smaller R values) there is only a modest dip in n mp .
The drop and subsequent growth of n mp that occurs at the same time that r undergoes rapid decay isn't sensitive to the overall length of the decay time. Figure 5(b) shows that the variation of n mp during the final rapid decay always takes the same form. Figure 6 provides a detailed illustration of the latching dynamics during the final rapid decay of the transient synchronised state. As the state decays, the number distribution P(n) initially broadens and its peak drops to zero (for large enough R). The distribution then narrows as its peak moves to higher occupation numbers. The Wigner function [49] (which we calculated numerically using QuTiP [51]) provides additional information about the phase during this transition. Throughout the slow decay stage it has a single well-defined peak centred away from the origin which indicates that the system has a preferred phase. During the final rapid decay, the Wigner function peak becomes smeared out around a central dip. Finally, the Wigner function reaches a limitcycle state with no phase preference-matching the steady state of an equivalent uncoupled vdP oscillator [13].

Critical scaling
The critical scaling of the relaxation time in a classical dynamical system just below a saddle-node bifurcation (such as that described by the simple potential model given by (6) and (7)) takes the universal form [35,47] with b=−0.5. To see whether this matches the behaviour of our vdP system we looked in detail at the way in which the relaxation times of the order parameter grows with time as the coupling approaches the critical value. As figure 7 shows, the critical scaling in the relaxation time of the vdP system for different values of R is indeed consistent with the predictions of the classical bifurcation model. Since we had no a priori knowledge of the precise value of the critical coupling strengths, we obtained critical exponents by assuming the relation given by (13). We then carried out linear fits for a range of choices of critical coupling, making no assumption about the value of the exponent b. The values of ε c used in figure 7 are those for which the best fits to the assumed relation were obtained and the corresponding critical exponents were almost exactly −0.5 in each case. Since the relaxation behaviour for weak couplings is very different to that near the critical coupling (the relaxation time actually decreases with increasing coupling), we carried out our fits using only couplings for which the relaxation time was about at least as long as in the weak-coupling limit (the full circles in figure 7).
One very surprising feature of the dynamics is the wide range of couplings which seem to be well-described by (13). Typically, one expects critical scaling to apply in a rather narrow range around the critical point, but figure 7 shows that it can extend over a substantial range of coupling strengths. Furthermore, this range seems to increase with R, i.e. as the amplitude of the underlying limit-cycle gets smaller.

Conclusion and discussion
We have explored the synchronisation dynamics of a model describing a large number of quantum vdP oscillators with all-to-all couplings. Adopting a quantum mean-field description, we carried out numerical integrations of the resulting effective nonlinear quantum master equation. For couplings above a critical value the system evolves towards either a synchronised state or an unsynchronised state, depending on initial conditions. We found that these states differ not just in their phase properties, but also in the properties of their number distributions. In a sense the synchronisation transition can be thought of as involving both phaseordering and number-disordering: the synchronised state has a well-defined phase, but its number distribution is always broader than that of the corresponding unsynchronised state which has a flat phase distribution.
We found that the dynamics of the system is rather rich with a number of interesting features. Just below the critical coupling the system always evolves towards the unsynchronised state, but the time taken to reach it can become extremely large. Looking at the dynamics of the slow relaxation process in detail we find that the system displays an interesting 'latching' behaviour in which the average occupation number drops before rapidly rising again as the system approaches the unsynchronised state. The relaxation time for couplings below the critical value displays a scaling which is consistent with that predicted by a simple classical model describing a system undergoing a saddle-node bifurcation. However, in the regime where the occupation numbers of the system are relatively small (and the behaviour is not well-described by semiclassical models [13]), the scaling seems to apply over an unexpectedly broad range of couplings below the critical value.
Although we have focussed on a specific oscillator model involving vdP oscillators, we expect our results to apply rather generally to quantum limit-cycle oscillators coupled via a simple coherent coupling. For example, very similar results are found [52] for many-body synchronisation in an analogous oscillator system which is instead based on the micromaser [25].
Our findings also suggest some promising directions for future work. It would be interesting to explore the dynamics within the vicinity of phase transitions in other non-equilibrium many-body quantum systems within the mean-field approximation. It will also be worth investigating whether the dynamics is qualitatively similar in models in which the quantum mean-field approximation is relaxed (e.g. cluster mean-field descriptions based on plaquettes with two or more oscillators [53]) or indeed in systems containing several rather than many oscillators (perhaps using models based on spin-systems which naturally have a much smaller state space [30,54]). Finally, the question of why the behaviour of the many-body vdP system can end up being well described by critical scaling for a wide range of parameters is surely worth pursuing.