Synchronization in a semiclassical Kuramoto model

Synchronization is a ubiquitous phenomenon occurring in social, biological, and technological systems when the internal rhythms of their constituents are adapted to be in unison as a result of their coupling. This natural tendency towards dynamical consensus has spurred a large body of theoretical and experimental research in recent decades. The Kuramoto model constitutes the most studied and paradigmatic framework in which to study synchronization. In particular, it shows how synchronization appears as a phase transition from a dynamically disordered state at some critical value for the coupling strength between the interacting units. The critical properties of the synchronization transition of this model have been widely studied and many variants of its formulations have been considered to address different physical realizations. However, the Kuramoto model has been studied only within the domain of classical dynamics, thus neglecting its applications for the study of quantum synchronization phenomena. Based on a system-bath approach and within the Feynman path-integral formalism, we derive equations for the Kuramoto model by taking into account the first quantum fluctuations. We also analyze its critical properties, the main result being the derivation of the value for the synchronization onset. This critical coupling increases its value as quantumness increases, as a consequence of the possibility of tunneling that quantum fluctuations provide.


I. INTRODUCTION
Synchronization is perhaps the most cross-disciplinary concept of emergence of collective behavior [1] as it is manifested across many branches of natural and social sciences. Ensembles of neurons, fireflies or humans are prone to synchronize their internal rhythms when they become coupled enough, producing a macroscopic dynamically coherent state. In all these seemingly unrelated situations, no matter the precise nature of the coupled units, interaction drives system's components to behave homogeneously. Thus, the study about the microscopic rules that drive ensembles towards synchrony has a long and fruitful history since the seminal observations made by Christiaan Huygens [2][3][4].
The mathematical formulation of the first models showing synchronization phenomena dates back to the 70's when, after some preliminary works by Peskin and Winfree [11], Kuramoto [5] formalized his celebrated model. The Kuramoto model incorporates the minimum dynamical ingredients aimed at capturing a variety of physical phenomena related with the onset of synchronization. In particular, the Kuramoto model links physical concepts such as self-organization, emergence, order in time and phase transitions, thus revealing as the most paradigmatic framework to study synchronization [1,7,11].
Despite the large body of literature devoted to the Kuramoto model and its variants, its study has always been restricted to the classical domain. At first sight, given the usual nature (scale) of the systems in which synchronization is typically observed, it seems superflous thinking of a quantum theory for the Kuramoto model. However, there is not doubt about the fundamental importance of studying quantum fluctuations within the emergence of synchronized states [9][10][11][12]. Moreover, the Kuramoto model has been implemented on circuits and micro and nanomechanical structures [13,14], systems which have already met the quantum domain [15,16]. At the quantum level, synchronization, understood as the emergence of a coherent behaviour from an incoherent situation in the absence of external fields, is reminiscent of the phenomena such as condensation of Bose-Einsten and has been observed in interacting condensates of quasiparticles [17,18]. Thus, moved by its fundamental and applied importance, in this work we provide the quantum version of the Kuramoto model in an attempt for understanding the influence that quantumness has on the emergence of synchronized states.
Our work in this paper consists, as stated by Caldeira and Leggett [19], on finding consistent equations that in the classical limit matches the Kuramoto model. Our derivation relies on the quantization of open systems in the framework of Feynman's path-integral formalism. Being the classical limit well defined, i.e. the Kuramoto model, and once we have derived its quantum counterpart we analyze its critical properties by deriving the critical point from which synchronization shows up and deter-mine how quantum fluctuations affect this synchronization transition. As a semiclassical approach the quantum Kuramoto model derived takes into account the first quantum fluctuations, and its results represent the first step towards the study of synchronization in quantum extended systems.

II. THE CLASSICAL KURAMOTO MODEL.
The original Kuramoto model [5] considers a collection of N phase-oscillators, i.e., it assumes that the characteristic time scale of their amplitudes is much faster than that for the phases. Thus, the dynamical state of the i-th unit is described by an angular variable θ i ∈ (0, 2π] whose time evolution is given by: The above equation thus describes a set of weakly coupled phase-oscillators whose internal (natural) frequencies {ω i } are, in principle, different as they are assigned following a frequency distribution g(ω) that it is assumed to be uni-modal and even around the mean frequency Ω of the population, g(Ω + ω) = g(Ω − ω).
In the uncoupled limit (K = 0) each element i describes limit-cycle oscillations with characteristic frequency ω i . Kuramoto showed that, by increasing the coupling K the system experiences a transition towards complete synchronization, i.e., a dynamical state in which θ i (t) = θ j (t) ∀i, j and ∀t. This transition shows up when the coupling strength exceeds a critical value whose exact value is: To monitor the transition towards synchronization, Kuramoto introduce a complex order parameter: The modulus of the above order parameter, r(t) ∈ [0, 1], measures the coherence of the collective motion, reaching the value r = 1 when the system is fully synchronized, while r = 0 for the incoherent solution. On the other hand, the value of Ψ(t) accounts for the average phase of the collective dynamics of the system. In Figure 2 we have illustrated the synchronization in the Kuramoto model. The panels in the top show, for different values of the coupling K, how the oscillators concentrate as K increases. Below we have shown the usual synchronization diagram r(K) for which the exact value of r for each K is the result of a time average of r(t) over a large enough time window. In this diagram we can observe that K c = 1 as a result of using the distribution g(ω) shown in the right.
Let us note that the all-to-all coupling considered originally by Kuramoto can be trivially generalized to any connectivity structure by introducing the coupling matrix A ij inside the sum in (1) so that each term j accounting for the interaction between oscillator i and j is assigned a different weight. The latter allows for the study of the synchronization properties of a variety of real-world systems for which interactions between constituents are better described as a complex network [20]. The formalism developed in this work is fully general and valid for any form of K ij thus making possible the extension of the large number of studies about the Kuramoto model in any topology [21] to the quantum domain. However, the numerical part of our work will deal with the all-to-all coupling for the sake of comparison with the original Kuramoto work.

III. QUANTIZATION OF THE KURAMOTO
MODEL.
The most important problem when facing the quantization of the Kuramoto model is its non-Hamiltonian character since, as introduced above, equation (1) assumes the steady-state for the dynamical state of the amplitude of the oscillators. Thus, a question arises, how do we introduce quantum fluctuations in the Kuramoto model? One possible choice is to resort to the original microscopic dynamics of amplitude and phases and then identify the underlying Hamiltonian dynamics. However, many different dynamical setups can have the Kuramoto model as their corresponding limiting case of fast amplitude dynamics. Thus, in order to keep the flavor of generality of the Kuramoto model, it is desirable not to resort to any specific situation (Hamiltonian) and introduce quantum fluctuations directly.
A similar problem was faced by Caldeira and Leggett in the eighties [19] when they studied the influence of dissipation in quantum tunneling. In their case, the corresponding classical dynamics dates back to the studies on activation theory by Kramers [22]. Classically, a particle in a potential experience an energy barrier to surmount, that is typically acquired from thermal fluctuations. On the other hand, a quantum particle finds in tunnelling an alternative way to bypass an energy barrier. Caldeira and Leggett were thus interested in quantifying the catalytic effect of tunnelling in (effectively) lowering the energy barriers. However, as in the Kuramoto model, Kramers activation theory is based in Langevin equations, i.e. stochastic equations that are not directly obtained from any Lagrangian. Furthermore, most of reaction rate equations were phenomenological. Therefore, they searched for a consistent way for introducing quantum fluctuations regardless of the microscopic origin of the effective classical evolution. As a byproduct their work opened the field of quantum Brownian motion in the most general way.
We take here the same route followed by Caldeira and ). The color of each oscillator represents its natural frequency. From left to right we observe how oscillators start to concentrate as the coupling K increases. In the panels below we show the synchronization diagram, i.e., the Kuramoto order parameter r as a function of K. It is clear that Kc = 1 as obtained by using the distribution g(ω) shown in the right panel.
Leggett to introduce quantum fluctuations in the Kuramoto model. In order to accomodate our dynamical system (1) to the framework provided in [19] we start by writting its corresponding Langevin equation: with As usual, ξ i is a Markovian stochastic fluctuating force with ξ i (t) = 0 and ξ i (t)ξ j (t ′ ) = 2δ ij Dδ(t − t ′ ). In the limit D → 0, equation (4) reduces to the Kuramoto model in equation (1). Equation (4) is nothing but a Langevin equation in the overdamped limit. It is first rather than second order in time as the inertia term is neglected. Consequently, the Kuramoto model can be viewed as a set of phases evolving in the overdamped limit. The absence of fluctuations in the limit D → 0 means that the system of phases is at zero temperature, D ∼ T . Such identification with a Langevin equation has been already used for generalizations of the original Kuramoto model taking into account noise and/or inertial effects [1]. In particular, in [10] it is shown that the critical value K c reads: which, in the limit D → 0, recovers the Kuramoto critical coupling (2).
The key point of deriving the Langevin equation (4) corresponding to the Kuramoto model is that it can be obtained from a fully Hamiltonian framework by coupling the system, in our case the coupled phases θ i , to a macroscopic bath or reservoir [22]. In this way, both the damping and fluctuations are seen to be caused by the coupling of the system of phases to the bath. The Hamiltonian description properly casted in the systembath approach reads as follows: where the bath is an infinite collection of harmonic oscillators with frequencies {ω α } (note that greek subindexes will denote the oscillators in the bath). In the case we are dealing with the total Hamiltonian reads: where {(θ i , π i )} and {(Q i,α , P i,α )} denote the system and bath canonical coordinates, respectively, while λ α stands for the coupling constant between bath and system coordinates.
Under well defined conditions, the equations of motion for the system coordinates derived from the Hamiltonian (8) lead to the the afore-derived overdamped Langevin equation (4). In particular, one needs to assume: (i) thermalized initial conditions for the bath: (ii) the frequency spectrum of the bath oscillators is flat (this assumption leads to the widely used Ohmic dissipation), and finally (iii) the changes in time of the velocity (acceleration) induced by the energy potentials are far slower than the energy loss induced by the coupling between the system and the bath (this is the situation when the system and the bath are strongly coupled) so that we could neglect the inertial term.
Once we have a Hamiltonian description for the Kuramoto equation (1), we are ready to perform its quantization. First, we associate the phases and their associated momenta together with the positions and momenta for the bath by providing them with the canonical commutation rules. The hardest work is to find an effective quantum evolution depending only on phase operators, i.e. the so-called quantum Langevin equation. It turns out that such operators equation is a non-local differential equation, which makes it extremely difficult to manipulate in general. However, the quantum version of equation (4) in the overdamped limit is a c-number local differential equation [3,5,23,[26][27][28]. The resulting quantum evolution reads as follows: where we have used the compact notation V ′...′ i,...,k ≡ ∂ θi,...,θ k V while Λ is a parameter tuning the degree of quantumness and ξ i is an stochastic force with the same statistics as in (4). Let us note that the quantum model is perturbative in Λ, which has a clear physical meaning. It stands for the quantum corrections to the thermal averages for the variance of the phases θ i . In the Supplementary Material [29] we provide the details of the derivation of the model. Once we derived the quantum version of the Kuramoto equation, it is natural to unveil the effects that quantum fluctuations induce in the transition to synchronization. As introduced previously, to study the synchronization transition one resorts to the order parameter r [introduced in equation (F33)] that reveals the synchronized state of the system. We solve both the classical Kuramoto model (Λ = 0) and the quantum one (Λ > 0) for the classical (Λ = 0) and the quantum (Λ = 0.1) Kuramoto models. In both cases the thermal noise is chosen such that D = 1. The number of oscillators is N = 10 3 and the distribution of natural frequencies is Lorentzian: g(ω; ω0, α) = π −1 α/[(ω − ω0) 2 + α 2 ] centered in ω0 = 0 and α = 0.5. It is clear that the synchronization onset is delayed as soon as quantumness enters into play. In panels (b) and (c) we show the probability P (θ) of finding an oscillator at a given phase θ as a function of K. Note that for each value of K, the phases has been equally shifted so that the mean phase is located at θ = π. A thick grey line indicates the critical values Kc and K q c for classical and quantum dynamics, respectively.

FIG. 3.
System of two coupled Kuramoto oscillators. The top-left part of the figure shows the analogy between the system of two coupled oscillators and an overdamped particle in a washboard potential. The two possible regimes are shown: synchronized state (the particle is at restφ = 0 at a local minimum) and unsynchronized phase (the particle drifts across the potential). Below we illustrate the possibility that tunneling provides to anticipate the drifting state. On the right part we show the result of the computation of the velocityφ as a function of ∆ω/K for the classical (red line) and quantum (blue line) systems.
numerically, extracting from the dynamics the stationary value of r. Through this work, the numerical calculations are performed with N = 10 3 oscillators and the distribution of natural frequencies is Lorentzian: with α = 0.5 and centered around ω 0 = 0. Figure 2.a shows the typical synchronization diagram, namely, the value of r as a function of the coupling strength K. The comparison of the quantum (for Λ = 0.1) and classical curves r(K) evinces that quantum fluctuations delay the onset of synchronization, i.e., the critical point K c is seen to move to larger values with Λ. We have also considered the evolution for the distribution of the phases as a function of K to monitor the microscopic fingerprint of the synchronization transition. These evolutions for the classical and quantum Kuramoto models are shown in figures 2.b and 2.c, respectively.
To explain the delay in the synchronization onset introduced by quantum fluctuations we resort to the simplest situation: two coupled Kuramoto oscillators. In this case the Kuramoto model (4) consists on just two coupled equations for the evolution of θ 1 and θ 2 . By taking the difference of those two equations and introducing as a new variable the phase difference, ϕ := θ 1 − θ 2 , we obtain for its evolution the following equation: The latter equation describes the evolution of an overdamped particle in a washboard potential (see Figure 3). With this image in mind, we map the synchronous movement of the two oscillators (defined as a state in which the frequencies of the oscillators are locked:θ 1 =θ 2 ) with the resting state of the overdamped particle inside a local minimum of the potential energy (φ = 0). On the other hand, when the two oscillators are not synchronized the particle drifts across the potential (φ = 0). Both situations are shown in Fig. 3. The quantum version for the diffusion of an overdamped particle in a periodic potential has been previously studied in Ref. [27]. The main result is that the scape rate of the particle, and thus its unlocking mechanism, is enhanced through quantum fluctuations. This effect can be seen as a consequence of the enhancement of the transition probability for energies below the height of the barrier which is nothing but the well-known tunnel effect [30]. In Fig. 3 we show, for both the classical and quantum (Λ = 0.1) systems of two coupled Kuramoto oscillators, the value ofφ = 0 as a functio of the ratio between the difference of the natural frequencies of the two oscillators |∆ω| and the coupling K. It is clear that, as stated above, quantum tunnelling facilitates the drift or, equivalently, delays the transition to the synchronous state.

V. ANALYTICAL EXPRESSION FOR THE SYNCHRONIZATION ONSET
Coming back to the original model of N interacting oscillators, we now make an analytical estimation of the value for critical coupling at which the synchronization transition occurs. The procedure is a generalization of the one presented in Ref. [10] and takes advantage of the mean field description of the Kuramoto model. The derivation ( presented in the Supplementary Material [29]) yields a rather simple equation for the critical coupling: being K c the classical critical value shown in Eq. (6).
The above result tells that quantum fluctuations act by effectively decreasing the coupling strength with the degree of quantumness Λ. Coming back to the physical image of a particle in a washboard potential, we can consider the effect of the quantum correction by considering the first and third terms in the right hand side of equation (F33). In this way, the quantum corrections can be casted in the form of an effective potential: that in the particular case of the washboard potential reads: The above equation makes clear that tunnelling is formally reflected by an effective barrier reduction that yields the observed shift to higher values for the critical coupling. Our analytical estimation for K q c (Λ) is plotted in figure  2 (vertical arrows) confirming its validity. To corroborate further the correctness of equation (15), we explore the synchronization transition for different values of Λ in figure 4.a. As expected the onset of synchronization shifts to higher values as the degree of quantumness increases. Again, the predicted value for K q c is plotted (vertical arrows) corroborating the validity of equation (15).
To complete our study, we show in 4.b the dependence of the synchronization diagram with the thermal fluctuations, D, both for the classical and quantum (Λ = 0.1) cases. In all the curves explored the coupling K is rescaled by the corresponding critical coupling K c in the classical regime. In this way we show both for the classical and quantum cases, the robustness of the critical value (15) against temperature changes.

VI. DISCUSSION
The search for quantum corrections to classical phenomena has been pervasive in physics. Some examples related to our work are the generalization to the quantum domain of chaos [32], dissipation [33], random walks [34], etc. Each of these examples finds its own difficulties when incorporating quantum fluctuations and unveiling their role. Some of these obstacles are the quantum linearity versus the typical non-linearity of classical systems and the quantization of non-Hamiltonian system or phenomenological equations. Overcoming these obstacles provides with a consistent quantum description that opens the quantum door to a variety of classical problems and their associated physical phenomena.
Among the most studied phenomena in (classical) complex systems is synchronization. This emergent phenomena is as intriguing as beautiful, since it covers from the description of the sympathy of clocks to the neuronal functioning in our brain, thus overcoming the disparately diversity in the spatial and time scales associated to the bunch of systems in which synchronization is observed. However, the concept of synchronization was usually associated to the classical domain as the typical examples of clocks, fireflies or humans are too macroscopic to think about the need of introducing quantum fluctuations in the description of the associated dynamical models.
Recently, some experimental works have shown that synchronization can be observed in the lab within Josepshon Junction arrays [35], nanomechanical [13] or optomechanical systems [36]. All of these systems share one prominent property: they behave quantum mechanically at sufficiently low temperatures. Therefore, adapting the concept of synchronization among coupled entities within the quantum theory is, apart from an interesting theoretical issue, a must imposed by the rapid experimental advances.
A first step consists in taking the most widely used framework for studying synchronization phenomena, the Kuramoto model, and adapting it to the quatum domain. Being a paradigmatic theoretical setup, the quantization of the Kuramoto model opens the door to the theoretical study of quantum synchronization in the widest possible manner. To this end, and to overcome the non-Hamiltonian character of the Kuramoto equations, we have mapped the model to an overdamped Langevin equation which has a Hamiltonian description by embedding the system in a bath of oscillators. In this way, the quantization of the Kuramoto model is straightforward and it includes its classical counterpart as a limiting case: the quantum version incorporates quantum fluctuations for the phases while the strength of these quantum corrections are encoded in a single parameter.
The route chosen here must be understood as complementary to the study of particular models of coupled quantum systems. The reason is twofold. First, we aim to be as general as possible. The essence of an emergent phenomena is its ability of describe very different situations with different microscopic dynamics. This is the goal of the Kuramoto model, as it explains the synchronization without resorting to the specific dynamics. Second, a force brute study of many body quantum entities is a very difficult task that usually implies the reduction of the system to a few coupled systems. However, the observation of a true synchronization transition demands hundreds or thousands of interacting dynamical systems.
Being general, the results obtained allow to make general statements about the impact that quantumness has on the synchronization of coupled dynamical units. The most important one is that quantum fluctuations delay the appearance of a synchronized state. The explanation of this effect relies on the fact that in the quantum domain the phases not only have a different natural frequency but also the fluctuations around the classical trajectories are different depending on those internal rhythms. To illustrate this interpretation we recall the simple case of two coupled Kuramoto oscillators. In this case quantum fluctuations are nothing but thermal assisted tunneling favoring the phase unlocking. Therefore, the coupling needed to synchronize the two oscillators is higher in the quantum limit.
Finally, we want to point out that in a recent publication the question about synchronization in quantum evolutions was also discussed [9]. Under rather general conditions they find bounds for the degree of synchronization based on the Heisenberg uncertainty principle: the phases, derived as averages of non-conmuting operators, cannot take values infinitely close. Instead, in our case, focused on the quantum version of the Kuramoto model, we have discussed, not the maximum degree of synchronization but the critical onset for the appearance of partially synchronized states. In this case quantumness also limits the emergence of a synchronous state. Therefore, pretty much like in what happens in quantum chaos, synchronization seems to be a quasi-classical phenomena [11].

Appendix A: The system-Bath Approach
Through the work we make use of the system-bath Hamiltonian. Let us introduce the main notation and quantities of interest.
In what follows we are interested in thermal environments. In other words, we are imaging a system at some temperature T . On general grounds, the physics can be described by the system-bath Hamiltonian: here we have chosen latin indexes to label the number of particles of the system, while greek ones stand for the bath infinite degrees of freedom. Notice that each particle is coupled to one bath. A key quantity in this formalism is the spectral density: which condenses the relevant information of the bath modes.

Appendix B: The Kuramoto Model as a Langevin Equation
It is shown in many places that by choosing the so-called Ohmic spectral density: the classical description of the dynamics of a system coupled to a thermal reservoir can be formulated by means of the Langevin (L) equation This effective stochastic formulation is driven by the Markovian zero-mean-fluctuating force ξ i (t), which is characterized by By defining the characteristic scales for time, position and energy: together with the dimensionless quantities: the Langevin equation (B2) can be rewritten as: The stochastic noise is now dimensionless with statistics: In the overdamped limit γτ 0 ≫ 1, which means that γ −1 is the fastest time scale, the dimensionless Langevin equation the above can be approximated by:θ Finally, by defining the Kuramoto potential [1], it turns out that the corresponding overdamped Langevin dynamics is nothing than the Kuramoto model with driving noise. At zero temperature, D → 0 the equation (B9) reduces to the original Kuramoto model.

Appendix C: The Fokker-Planck Equation
In Hamiltonian dynamics, the expectation value of a given observable A can be computed by means of a normalized probability density P , dxdẋP (x,ẋ, t) = 1, as The dynamics of P are given by the Liouville equation,

In the presence of a thermal bath the dynamics is not longer Hamiltonian and the Liouville equation needs to be complemented. This type of complementation is what takes the Newtonian equation of motion into the Langevin equation [Cf. equation (B2)]. For the case of the Liouville equation this refinement of the dynamics is incorporated in
the so-called Fokker-Planck (FP) equation [2], which in the overdamped regime leads to the Smoluchowski equation (S). In the notation of phases appearing in K -model, the S equation can be written as For the case on non-interacting particles, V (x) = 0, the S equation leads to the Diffusion equation. For future convenience, we re-write equation (C3) as with where The FP equation and the way we have rewritten it in equation (C5) will be crucial in making a connection with the quantum case. Moreover, it may be useful for analytical investigations. In particular, if one is interested in stationary solutions: ∂ t P = 0, this implies that L cl P = 0, which in some cases can analytically evaluated [2].

Appendix D: Quantum Master Equations
In the next sections we detail the technical steps for obtaining the quantum generalization to the classical Langevinlike equation (B9). In the same way as for the classical Kuramoto model, we assume to be in the overdamped limit. It has been shown that the overdamped regime corresponds to a kind of semi-classics [3]. This means that the equilibrium density matrix is essentially diagonal in the position basis: . Following Ref. [3] our task consists on i) finding the reduced equilibrium density matrix in the system-bath model, ̺ β and ii) propose a Fokker-Planck equation like equation (C4) such that the stationary state is that ̺ β .

The reduced density matrix
We are interested in the dynamics of the phases. Therefore we must integrate out the bath degrees of freedom. In a more quantum language, we are interested in observables acting on the phases' Hilbert space. Therefore, the object of interest will be: Here W is the total density matrix and Tr bath means partial trace over the Hilbert space for the bath. With this definition, the averages over any observable acting on the phases (the system in this case) is: i.e. the averages can be done by considering the reduced density operator ̺ defined in equation (D1).

Appendix E: Equilibrium Density Matrix: Path Integral Formalism
As we anticipated, we need first to compute the equilibrium density matrix. In particular we are interested in the reduced density matrix (at equilibrium): where W β is the total equilibrium density operator, W β ∼ e −β(Hsys+H bath +Hint) . The equilibrium reduced density matrix can be expressed as [4] ̺ β (x, with the effective action which contains the kernel being ν n the Matsubara frequencies, and the Laplace transform of the damping kernel is given by: (E6)

Overdamped Equilibrium
Based on previous works [3,5] for the single particle case, we compute the equilibrium distribution in the overdamped limit. As we already mention in this supplementary material, the overdamped dynamics refers to a regime in the parameter space where the damping is sufficiently strong to suppress the non-diagonal elements, coherences, of the reduced density matrix, i.e., a regime where x 1 , ..., These semiclassical diagonal contributions can be computed perturbatively on the quantum fluctuations, to be defined below. The versatility and power of the of the path-integral approach in the semiclassical regime is what justifies its use here. In fact, the zeroth order calculation corresponds to compute the minimal path action for the diagonal contributions.

a. Minimal path
Let us denote the minimal action (ma) path as x ma i ≡x i . Besides, since we are interested in the diagonal contributions in the imaginary-time path integral in equation (E2), this means for us to take the trajectories with i.e. periodic trajectories with frequencies ν n . The minimal action path satisfies the generalized Lagrange equations [6] The periodic condition in equation (E7) suggests to Fourier expandx i (τ ), such that where the Fourier components satisfy and the inhomogenous term comes from the jumps and cups singularities arising from fact that the Fourier series expansion forx i (τ ) periodically continues the path outside the interval 0 ≤ τ ≤ β [6]. Note that terms like a i =x i ( β) −x i (0) are, in general, expected; however, since we are interested in the diagonal contributions, they do not contribute to the present case. At this point, we first notice that by making n = 0 for b i we obtain Besides, the components x n,i with n = 0, are suppressed by dissipation. Hence where Λ measures the quantumness: being C = 0.577... the Euler-Mascheroni constant. Note that in the limit → 0, Λ → 0, as it must be. Thus, recovering the classical result. The contribution of the minimal action can be further simplified by considering that together with (E7) and replacing equation (E8) in the second term at the r.h.s of equation (E17), such that By using the relation (E13) and by noticing thatx i ∼ = x 0,i [x n,i are suppressed, see equation (E14)], we have that

. Fluctuations around the minimal action path
We study now the fluctuation around the minimal path subjected to the boundary conditions: Consequently, the correction to the path integral reads, where we have used an economical notation,à la Dirac, for the quadratic form y|L|y = L ij y i y j , being L = {{L ij }} defined as where I is the identity matrix and the second-derivative-potential-matrix V ′′ = {{V ′′ ij }} is defined as, We proceed as above and Fourier expand the fluctuations around the minimal path [Cf. equation (E21)], where µ|y n = i µ i y n,i . Therefore, with, A n = Iλ n + V ′′ and λ n = ν 2 n + |ν n |γ.
This a Gaussian integral that can be performed by resorting twice to the formula Up to first order in 1/γ, we get [Cf. equation (E28]: To be consistent, we also need to compute the determinants at first order in 1/γ [7] det A −1 Based on all the consideration above, in the next appendix we explicitly present the thermal equilibrium state with first order corrections in the fluctuations along the semiclassical minimal path results and derive the associated Smoluchowski equation.

Appendix F: The Quantum Master Equation for the Kuramoto Model: q-K
We proceed here as Ankerhold et al. in Refs. [3,5]: (i) We compute the equilibrium density matrix ̺ β in the overdamped limit. (ii) We assume that in the quantum regime the qme will have the form equation (C4) and equation (C5). (iii) Finally, we find the the coefficients in equation (C4) and equation (C5) by forcing the equilibrium density matrix to be stationary, i.e., LP β (q) = 0-note that this a natural condition. This is a consistent equilibriumbased derivation for the QME. Since we are going to be interested in the long time dynamics (stationary), by imposing the correct equilibrium in the dynamics we consider ourselves in a safe place.

Thermal Equilibrium Density Matrix ̺ β
Based on the result obtained in section E, the equilibrium density matrix in the overdamped limit reads: where we have introduced the notation P β (q) since in the overdamped limit only the diagonal elements ̺ β (q, q) matter.

One-Particle Master Equation
Let us consider first the one-particle model. This will be a useful analogy later. Besides, since the one particle case has been studied in Refs. [3,5], we can use it for the purposes of comparison.
In the classical case, the Fokker-Planck equation can be expressed as with, where Γ denotes mγ.
To derive the corresponding quantum Fokker-Planck equation, let us start by considering the quantum equilibrium for the single particle case, which after equation (F1) reads: where Z is the partition function and Up to leading order in Λ, We choose for D 2 : where F = 1−βΛV ′′ is called the quantum correction factor [3,5] and is the same factor that we found in equation (F6). Imposing L̺ β = 0 together with this election for D 2 , we end up with a D 1 coefficient given by It is easy to calculate the term ∂ x S: So the term D 1 is: This yields the QME for the single case in the overdamped limit:

N -Particles Master Equation
Let us obtain the equations in the multivariable case. With N particles the problem is quite more difficult. The first and second derivatives of the potential are not scalars, but a vector V = {V i } and a matrix V = {{V ij }}, whose components are given by where i, j run from 1 to N , the number of particles. Again, we have to find the Fokker-Planck coefficients which leads to a stationary solution of the Fokker-Planck equation. In the multivariable case, the FP equation is the same that equation (F2), but instead of equation (F3) L reads now: whereas the stationary solution P β in equation (F1) can be rewritten as, where Tr(V ′′ ) denotes trace of the matrix V ′′ . The stationary solution has the general form 1 Z F (x)e S . We can choose freely what part of P β is included in F and what part in e S . The form of the D 1,i and D 2,i terms will be determined by our election of F .
There are three immediate choices: F = 1 is the simplest one, but it does not lead to the quantum harmonic oscillator in the one-particle limit. The second choice is an F , like the one-particle, e −βΛTr(V ′′ ) . This has one inconvenient: in the uncoupled case, we will have, for each particle, an equation that depends on the factor F ; moreover, with that choice, F depends on all the oscillators in the system. This problem can be solved with the third choice: The factor F can be seen as a product of factors F i , which are uncoupled from each other. This third option is the one we follow below.
Before giving a particular functional form to F , we will find the general form that the master equation will have depending on our choice. Let us do something similar to what we did in the single-oscillator case. We consider a quantum white noise term and then impose the stationary solution in order to find the appropriate D 1 .
Since we want to simplify the term D 2 P β in the equation LP β , we write D 2,i = D γ 2 Fi . By Imposing the stationary solution P β together with our election for D 2,i , we find that We finally obtain: Now we have the general form of the terms D 1,i and D 2,i . If we now apply our particular choice for F i , we will obtain the master equation of the system. If P β has the form 1 Z F (x)e S , we easily identify the quantum correction factor F : This F can be separated in N factors F i as follows: After substituting these factors in the expression for D 2,i , we obtain: By doing so for D 1,i and by calculating the term ∂ xi S: as well as the the partial derivate F ′ j,i : we finally get that The expressions (F21) and (F24) allows us for explicitly writing the Fokker-Planck equation that describes a system of many quantum particles in the Smoluchowski regime: This expression constitutes the main core of our quantization of the Kuramoto model.

The Langevin Equation
Once we have derived the master equation, we can easily find the associated Langevin equation, in the form following the guidelines explained in Ref. [9]. Here ξ k is Gaussian δ-correlated white noise with zero mean and variance 2D [Cf. Eqs. (B3) an (B4)]. Let us write the master equation in terms of the Langevin coefficients So that the Langevin coefficients reads: The full Langevin equation reads Let us discuss the third and the fourth terms in detail. If we calculate the third partial partial derivate of the Kuramoto potential (note that this happens only to this potential as a particular case), we realize that third and the fourth terms are: Note that when i = j, the sin function vanishes, so the fourth term is exactly one half of the third term. Thus, the equation is simpler than it looks like. We can rewrite the equation in a Euler scheme, for convenience, just multiplying the equation by mγ. This reads:

Classical and One-Particle Limits
The quantum Langevin equation (F33) must reduces the classical Langevin in the limit when Λ → 0 and to the one-particle one when N = 1.
It is easy to note that the equation leads to the classic one: when Λ → 0, the quantum correction factor F tends to 1. Most authors take also m = 1 and γ = 1, so Γ = mγ = 1 just for simplify. So the master equation (F25) trivially leads to: which is the classical master equation for N particles.
Recall now the quantum master equation (F25) and consider N = 1: As we already did in the one-particle case [see equation (F8)], we approximate the exponential functions up to leading order in Λ, so that for F we have The D 1 term can now be read as: whereas the term A i in the Langevin equation reads now Taking a common factor −V ′ and simplifying, this is just:

Nondimensionalization of the Langevin Equation and the Quantum Kuramoto model
We proceed here as in the classical Langevin equation (B2) where some dimensionless quantities (B6) were introduced for finding the Kuramoto model. We start by writing the quantum langevin equation in the overdamped limit Using now the definitions (B5) and (B6) together with the scaling for the quantum parameter: we can arrive to:
1. Periodicity and self-consistency of the master equation Proceeding as in [10] for the classical case, we use a mean field approach for computing the critical coupling strength, K q c . In doing so we recall that the order parameter r is given by: In a mean-field approximation, we have to rewrite our equations in terms the parameter r. Since we have expressed our equations in terms of V and its derivatives, we only need to write them in terms of r V = −ωψ − Kr cos ψ, where ω = ω i − ω 0 and ψ = φ i − ω 0 t − φ 0 . Note that if ψ → ψ + 2π, then all the derivatives of V do not change; however, V →V − 2π. Hence, because of the presence of V , our stationary solution in equation (F16) does not satisfy the periodic boundary condition P (ψ; ω) = P (ψ + 2π; ω). Therefore, we have to find a periodic stationary solution.
Hence, we have chosen the second term in such way that it yields a factor e −2πβω that cancels the contribution from the original stationary solution.
Proving that P (ψ; ω) in equation (G3), in deed, satisfy the periodic boundary conditions, let us define f (φ) = e βV eff and consider Next we introduce the change of variable φ = τ + 2π, so that the last two terms read now As a consequence, a factor e −2πβω has appeared from the V term in the exponentials. Let us introduce a second change of variables. In particular, for the first term consider τ → φ whereas for the second consider τ → φ − 2π (such that we undo the first change for the second term). This leads to e −2πβω Hence, this solution does satisfy the stationary boundary condition.
If we then apply the master equation and consider that (∂ θ P 0 ) = −βV ′ eff P 0 , then we have that the third and sixth terms cancel each other. In the other terms, after taking ξ as a common factor, we end up with the original solution, which we know is 0. Thus, explicitly shows that the periodic solution is stationary.

Critical value
Once we have this solution, we can proceed as Sakaguchi [10] in order to find the critical coupling strength K q c . The order parameter r can be expressed in terms of ψ as: Replacing (G5) above we have a self-consistent equation for r. In the right hand of (G9), the imaginary part is always zero, because g(ω) is symmetric about ω = 0. The real part can be expanded in powers of Kr/D as: Noticing again that g(ω) is an even function the linear term πω/D do not contribute to the integral. Finally, the critical coupling strength, as a function of the temperature, is obtained from (G11), and it is: As K increases, a non-trivial solution branches off the trivial solution r = 0 at K = K c . This solution reduces to the classical one [10,11] when Λ = 0 at the classical critical coupling strength K c c . In fact, from above a simple relation between the classical and the quantum critical values can be obtained Then, we have a simple relation between the quantum and classic cases.