A stochastic approximation for the finite-size Kuramoto-Sakaguchi model

We perform a stochastic model reduction of the Kuramoto-Sakaguchi model for finitely many coupled phase oscillators with phase frustration. Whereas in the thermodynamic limit coupled oscillators exhibit stationary states and a constant order parameter, finite-size networks exhibit persistent temporal fluctuations of the order parameter. These fluctuations are caused by the interaction of the synchronized oscillators with the non-entrained oscillators. We present numerical results suggesting that the collective effect of the non-entrained oscillators on the synchronized cluster can be approximated by a Gaussian process. This allows for an effective closed evolution equation for the synchronized oscillators driven by a Gaussian process which we approximate by a two-dimensional Ornstein-Uhlenbeck process. Our reduction reproduces the stochastic fluctuations of the order parameter and leads to a simple stochastic differential equation for the order parameter.


I. INTRODUCTION
Ever since Huygen's observation of two pendulum clocks, mounted on the same wall a short distance apart, ending up swinging in anti-phase, the phenomenon of collective behaviour and synchronisation of weakly coupled oscillators has fascinated scientists.Synchronization has been observed in a diverse range of natural and engineered systems [1][2][3] , including in pace-maker cells of circadian rhythms 4 , networks of neurons 5 , in chemical oscillators 6,7 and in power grid systems 8 .
The celebrated Kuramoto model of sinusoidally coupled phase oscillators has served as a rich model to study synchronisation 1,2,[9][10][11][12][13][14] .The Kuramoto model was extended to the Kuramoto-Sakaguchi model 15 to include the effect of time-delayed or phase frustrated coupling which was observed in numerous real-world contexts 16 , including in arrays of Josephson junctions [17][18][19] , in power grids 20 and in seismology 21,22 .Apart from the usual transition from an incoherent state at low coupling strength to synchronisation upon increasing the coupling strength, the Kuramoto-Sakaguchi model exhibits a plethora of dynamical behaviours including bi-stability of incoherence and partial synchronisation, transition from coherence to incoherence with increasing coupling strength 23,24 , chaotic dynamics 25 as well as chimera states [26][27][28] .
The possible high dimensionality of networks of oscillators inhibits an understanding of the underlying dynamic mechanisms which give rise to this rich behaviour.Scientists have therefore looked at model reductions of the Kuramoto model and of the Kuramoto-Sakaguchi model.Most methods are restricted to the thermodynamic limit of infinitely many oscillators.In this limit Kuramoto and Sakaguchi established a mean-field theory which determines the order parameter and the non-zero rotation frequency of the synchronised cluster via a self-consistency relationship 15 .Similarly, the celebrated Ott-Antonson ansatz 29 can be employed to obtain a deterministic evolution equation for the order parameter 23,24 .Real a) wyue8667@uni.sydney.edu.aub) georg.gottwald@sydney.edu.auworld networks, however, are of finite size, and in finite-size systems the onset of synchronisation occurs typically not at the critical coupling strength predicted by the thermodynamic limit.While in the thermodynamic limit, the order parameter asymptotically in time approaches a constant value, in finitesize networks the order parameter exhibits persistent temporal fluctuations.Furthermore, finite-size networks exhibit a singularity in the variance of the order parameter at the onset of the transition to synchronisation with a well-defined finite size scaling [30][31][32] .To overcome the restriction of the thermodynamic limit and to tackle the case of finite-size networks, a collective coordinate approach [33][34][35] was applied to study the Kuramoto-Sakaguchi model for finitely many oscillators and derive an evolution equation for the order parameter 36 .It was established that the order parameter and the dynamics of the synchronised cluster is markedly influenced by the dynamics of the non-entrained rogue oscillators.Describing the effect of the rogue oscillators by their average, a deterministic reduced evolution equation for the entrained oscillators was derived.This allowed for the estimation of the onset of synchronisation and the determination of the average behaviour of the order parameter.
Whereas our previous work captured the averaged effect of the rogue oscillators on the synchronised oscillators, in this work we set out to quantitatively describe the fluctuations around this average behaviour and capture the effective stochastic dynamics of the synchronised phase-oscillators and the order parameter.We show how the thermodynamic limit is approached for increasing number of oscillators.We first present numerical results suggesting that the fluctuations constitute a Gaussian process.In a second step we approximate the Gaussian process by an easily computable Markovian Ornstein-Uhlenbeck process.This allows for a stochastic model reduction of the deterministic Kuramoto-Sakaguchi equation.In particular, we will propose a closed set of equations for the entrained synchronized oscillators where the effect of the non-entrained rogue oscillators is modelled by coloured noise, the variance of which decreases as the number of oscillators increases.We show numerically that the statistics of the order parameter and of the mean phase is well recovered by the reduced stochastic equation for the entrained oscillators.The stochastic approximation depends only, as we show numerically, on the number of the non-entrained rogue oscillators N r .For N r → ∞ the variance of the noise goes to zero and the deterministic thermodynamic limit is approached.On the other extreme, the validity of the stochastic approximation requires a sufficiently large number of rogue oscillators for the validity of the central limit theorem.This implies, as we will show, that to observe an effective noise requires different network sizes for different values of the coupling strength: for large values of the coupling strength with a high degree of synchronization there are proportionally fewer rogue oscillators than for smaller values of the coupling strength, requiring a larger total number of oscillators N to allow for sufficiently large numbers N r of rogue oscillators.We further use the stochastic model reduction to formulate a stochastic differential equation (SDE) for the order parameter.SDEs driven by Brownian motion for the order parameter were recently postulated and determined using a data-driven approach 37 .We show here that the SDEs are in fact driven by coloured noise.
There are two dynamic mechanisms capable of generating effective stochastic behaviour in deterministic systems: multiscale dynamics and weak coupling 38 .Both mechanisms generate stochasticity via some (functional) central limit theorems.In multiscale dynamics, slow variables experience on a long diffusive time scale the integrated effect of fast variables.If the fast dynamics is sufficiently chaotic, the integrated effect may lead to Brownian motion [39][40][41][42][43] .In weakly coupled systems in which resolved variables are weakly coupled to a large number of unresolved degrees of freedom, the sum over the many unresolved variables with uncorrelated initial conditions leads to a Gaussian process 38,44,45 .We establish here that the dynamic mechanism responsible for the effective stochastic dynamics of the collective behaviour of large but finite Kuramoto-Sakaguchi models is via the route of weak coupling.
The paper is organised as follows.In Section II we introduce the Kuramoto-Sakaguchi model.Section II A reviews the mean-field theory for the Kuramoto-Sakaguchi model.Section III presents numerical simulations of the Kuramoto-Sakaguchi model illustrating its stochastic order parameter fluctuations.We develop our stochastic model reduction in Section IV, in which we formulate an SDE for the entrained oscillators.We further use the stochastic model equations to formulate an SDE for the order parameter.In Section V we present numerical results of the reduced stochastic model and its ability to capture observed statistical properties of the full Kuramoto-Sakaguchi model.We conclude in Section VI with a discussion.

II. THE KURAMOTO-SAKAGUCHI MODEL
The Kuramoto-Sakaguchi model 1,9,15 θi is a classic paradigmatic model describing the dynamics of N sinusoidally globally coupled oscillators with phases θ i under global phase frustration λ with coupling strength K.Each oscillator is equipped with an intrinsic frequency ω i which is drawn from a specified distribution g(ω).
The Kuramoto-Sakaguchi model (1) displays a transition to synchronisation as the coupling strength K increases.For low values of K, the system is in an incoherent state in which each oscillator evolves approximately with their own intrinsic frequency.As the coupling strength K increases, some of the oscillators become synchronised, oscillating with a common frequency and with their phases staying close to one another.When K increases further, more and more oscillators become synchronised, until eventually global synchronisation occurs.The non-zero phase frustration λ ̸ = 0 induces a collective rotation of the synchronised cluster with a non-zero frequency Ω in the rest frame, as opposed to the Kuramoto model which supports stationary synchronised clusters.We remark that for λ = 0 and uniform intrinsic frequency distributions a firstorder phase transition occurs 46 .
The collective behaviour can be described by the mean-field variables r and ψ with The degree of synchronisation is quantified by the order parameter r with r = lim Perfect phase synchronisation with θ 1 = θ 2 = • • • = θ N implies r = 1 and r ≳ 0 indicates that phases are spread out with r ∼ 1/ √ N indicating incoherence.

A. Classical mean-field theory
Sakaguchi and Kuramoto 15 developed a mean-field theory for the mean frequency Ω of the synchronzied cluster and the order parameter r which we present here for completeness.We follow here their exposition to obtain self-consistency relations for the mean-field variables r and Ω as well as expressions for the stationary density function.The stationary density will be used later to determine the average effect of the rogue oscillators onto the synchronized oscillators.The assumption in the Kuramoto-Sakaguchi model (1) that each oscillator is connected to all other oscillators is crucial for the following derivation.Moving into the frame of reference rotating with the cluster mean frequency Ω = Ω(K) and setting the mean-field phase variable ψ = 0, the Kuramoto-Sakaguchi model (1) is expressed as with frequency Each oscillator θ i only couples to the other oscillators via the mean-field in the form of r and Ω.
From (3) one can readily identify the oscillators which form the synchronised cluster and the rogue oscillators which are not entrained: The former ones have frequencies |ω i − Ω| ≤ Kr for which (3) has stationary solutions with The synchronized oscillators rotate with the collective mean frequency Ω.In contrast, the non-entrained rogue oscillators have frequencies |ω i − Ω| > Kr and satisfy the Adler equation (3).
In the thermodynamic limit N → ∞, the phases can be described by a probability density function ρ(θ ,t; ω) satisfying the continuity equation For stationary states with ∂ ρ ∂t = 0 there are two kinds of stationary density solutions, depending on the intrinsic frequency ω.Entrained oscillators with |ω i − Ω| ≤ Kr which form the synchronised cluster are captured by the stationary probability density function The non-entrained rogue oscillators with |ω i − Ω| > Kr are captured by the stationary probability density function with normalization constant C(ω).In the thermodynamic limit the order parameter for the mean-field Kuramoto-Sakaguchi equation ( 3) is then given by separating the integration over the frequencies into the entrained and non-entrained ranges, with their respective stationary densities (6) and (7).We obtain Separating the real and imaginary parts of (8), we arrive at These are self-consistency relations for the mean-field variables r and Ω and can be numerically solved for any given intrinsic frequency distribution g(ω), coupling strength K and phase frustration parameter λ .In the following we denote the solution of the order parameter of the self-consistency relation by r ∞ .

III. FINITE-SIZE EFFECTS IN THE KURAMOTO-SAKAGUCHI EQUATION
We now present results of numerical simulations of the Kuramoto-Sakaguchi model (1), illustrating finite-size effects.We employ a 4th order Runge-Kutta(RK4) method with time step dt = 0.05.Initial conditions are chosen randomly from the interval [0, 2π]; to eliminate transient behaviour we discard a transient period of t 0 = 5 × 10 4 time units.Statistics such as means and variances are computed from time series of length t max = 5 × 10 4 .In the following we fix the phase frustration parameter to λ = π/4 and consider a Gaussian distribution of the intrinsic frequencies with mean 0 and variance 1.To avoid sampling effects which may lead to small cluster nucleation 47,48 we consider here equiprobable intervals between the N intrinsic frequencies 49 .This allows us to compare the effect of different system sizes N where for each system size N the intrinsic frequencies are uniquely fixed, and eliminates variations due to particular draws of ω i for a given system size N.For a detailed discussion on the effects of random sampling see, for example, 31,47 .Unless specified otherwise we set K = 3 which corresponds to a partially synchronized state of the system where both synchronized and non-synchronized oscillators occur.
Figure 1a shows a snapshot of the phases for N = 160.The phases are labelled according to their intrinsic frequencies such that the oscillator with the most negative intrinsic frequency is assigned the label i = 1 and the one with the largest frequency the label i = N.One can see clearly the synchronised cluster with frequencies close to the mean frequency Ω, and the non-entrained rogue oscillators with more extreme intrinsic frequencies.Note that due to the non-zero phase frustration λ the synchronised cluster is not centred around the oscillators closest to the mean of their respective intrinsic frequencies as is the case for the standard Kuramoto model with λ = 0. To determine the synchronised cluster and its comple-ment, the non-entrained rogue oscillators, we define the effective frequency of each oscillator, ωi = θi (t) t , where the angular brackets denote a temporal average.Those oscillators with approximately the same effective frequencies ωi are considered to form the synchronised cluster C of size N c , and the remaining oscillators are considered to be within the set of non-synchronised rogue oscillators R with size N r .Figure 1b shows the effective frequencies.We estimate the number of rogue oscillators for K = 3 and N = 160 oscillators to be N r = 43.We identify oscillators with indices i ≤ 116 as the synchronized oscillators with N c = 116 and oscillators with indices i > 116 as the rogue oscillators.
The effect of the non-collective behaviour of the rogue oscillators induces a seemingly stochastic behaviour of the order parameter r shown in Figure 2. We show the temporal evolution of the order parameter from a random initial condition at t = 0; after an initial transient the order parameter exhibits persistent fluctuations around a mean value with a constant variance of 4.16 × 10 −4 for N = 160 and 4.42 × 10 −6 for N = 10, 000.The size of the fluctuations is N dependent and we expect that for N → ∞ the variance of the fluctuations approaches zero as demonstrated in Figure 4 further down.Figure 3 shows that the overall mean phase ψ(t) as defined in (2) as well as the mean phase of the synchronized oscillators only, as defined further down in (15), exhibit diffusive behaviour.We checked that their dynamics is (approximately) Brownian with a linear time dependency of the mean-square displacement (not shown).Finite-size induced diffusivity of the phase has been previously reported for a stochastic Kuramoto model with λ = 0 and additive noise [50][51][52][53][54][55][56] .Here the effect is entirely a deterministic finite-size effect as there is no driving noise in the Kuramoto-Sakaguchi model (1).

A. Dependency of fluctuations on the system size N
To investigate the effect of finite system size, we simulate the system at increasing system size N = {40, 80, 160, 320, 640, 1280, 2560} with other settings kept the same.Figure 4a shows the distribution of r(t) for varying system size.As N increases, the distribution of r(t) becomes, as expected, narrower with decreasing variance and its mean value approaches a fixed value r∞ = 0.766, as shown in Figures 4b-c.We find that r − r ∞ ∼ N −0.859 .The variance of r(t) decreases towards 0 as N increases with Var[r] ∼ N −0.976 .The scaling of the variance with approximately with 1/N suggests an underlying Central Limit Theorem describing the finite-size fluctuations of the order parameter around the thermodynamic mean as a Gaussian process.
Figure 5 shows that the number of rogue oscillators N r very closely exhibits a linear dependency on the total number of oscillators N. The numerically estimated slope of 0.2767 is very well approximated by the thermodynamic limit lim where r and Ω are the solutions of the mean-field consistency relations (9- (10).For the parameters used in Figure 5 9)- (10).All parameters are the same as in Fig. 1.
fixed number of oscillators N = 160.Figure 6a shows the order parameter r with the well known second order phase transition for Gaussian intrinsic frequencies from incoherence with r ∼ 1/ √ N to partial synchronisation at K c = 1.93, after which the synchronised cluster continues to grow in size until all oscillators become synchronised at K g = 8.34. Figure 6b shows the variance of the order parameter as a function of the coupling strength.The almost constant variance for small coupling strength around Var(r) = 1.35×10 −3 corresponds to the almost random distribution of the N uncoupled oscillators.The variance exhibits a singularity at K = K c and approaches zero for K ≥ K g .The singular behaviour of the variance is a well-studied phenomenon and known under the name of anomalous enhancement of fluctuations 30,31,[57][58][59][60] .Figure 6c shows the cardinalities of the synchronised cluster C and the set of rogue oscillators R, labelled N c and N r , respectively, with N c +N r = N.For K < K c all oscillators are rogue whereas for K > K g all oscillators are synchronised.The synchronised FIG.3: Time series of the overall mean phase ψ(t) defined in (2) and of the mean phase ψ c of the synchronized cluster defined in (15) obtained from a simulation of the Kuramoto-Sakaguchi model (1) with N = 160 oscillators, exhibiting diffusive behaviour.All parameters are the same as in Fig. 1.
cluster steadily grows for increasing K c < K < K g .

IV. STOCHASTIC MODEL REDUCTION
The numerical results presented in Section III suggest that the non-entrained rogue oscillators exert a stochastic forcing on the synchronised cluster.To develop a stochastic approximation of the Kuramoto-Sakaguchi model we hence aim to establish a closed evolution equation for the phases of the synchronised oscillators in which the driving force exerted by the rogue oscillators is parametrised by a stochastic process.The challenge is how to describe this stochastic process.We begin by separating the coupling terms which only involve the synchronised oscillators and those which contain the rogue oscillators, and write the Kuramoto-Sakaguchi model ( 1) for

FIG. 5:
The number of rogue oscillators N r as a function of the total number of oscillators N. The continuous line shows a best fit with slope 0.2767 in good agreement with the thermodynamic limit relationship (11).All other parameters are the same as in Fig. 1.
The sum over the rogue oscillators can be written as where we introduced the mean phase of the synchronised cluster ψ c and where we define where we defined the mean-field variables pertaining to the synchronized oscillators in C r c e iψ c = 1 and those pertaining to the rogue oscillators in R r r e iψ r = 1 We remark that the particular choice of writing the trigonometric functions in ( 13)-( 14) was motivated by the form of the invariant density of the rogue oscillators (7) and allows for a convenient averaging (cf.( 18)-( 20)).The overall influence of the rogue oscillators on each synchronised oscillator θ i is now captured by S(t) and C(t), while the remaining terms are expressed entirely in terms of the synchronised oscillators, and we arrive at where we defined for compactness We first establish the average effect of the rogue oscillators and calculate the means ⟨S⟩ and ⟨C⟩ of the rogue oscillator drivers S(t) and C(t).We follow here the averaging procedure proposed in our previous work 36 .
In the co-rotating frame of the synchronized cluster, the synchronized oscillators are stationary whereas the rogue oscillators rotate (cf. Figure 1).Envoking ergodicity we can equate the temporal averages ⟨S⟩ and ⟨C⟩ by averages over the stationary density function of the rogue oscillators (7).Averaging equation (13) for S(t) over the invariant densities of the rogue oscillators (7) becomes where and ψ denotes the overall phase of all oscillators in the frame of reference rotating with the mean frequency Ω.Similarly, averaging (14) yields since the overall phase ψ (see (2)) and the phase of the synchronized cluster ψ c are close for sufficiently large coupling strength.We remark that for λ = 0 we have ⟨S⟩ t = ⟨C⟩ t = 0.In practice we estimate ⟨S⟩ t and ⟨C⟩ t numerically from a long time trajectory, which yields ⟨S⟩ t = 0.148 and ⟨C⟩ t = −0.0236,see Section III for the parameters used in the numerical simulations.These numerically estimated values are well approximated by the analytical expressions (18) and ( 20) which yield ⟨S⟩ t = 0.144 and ⟨C⟩ t = −0.0233.
To evaluate the analytical expressions ( 18) and ( 20) the constant values of r and Ω for the fixed finite system size N are required.These would only be available from numerical simulations.To remain purely analytical we use instead the values corresponding to the thermodynamic limit with r = 0.766 and Ω = −1.706,which are calculated via the selfconsistency relations ( 9) and ( 10) of the mean-field theory.
We remark that instead of using the finite size expressions (18) and ( 20) we could consider their thermodynamic limit by averaging over the frequency distribution g(ω).This yields ⟨S⟩ t = 0.149 and ⟨C⟩ t = −0.0240which is also close to the numerically observed values.
Figure 7 shows the empirical histogram of S and C obtained from a single long simulation of the Kuramoto-Sakaguchi model (1) for the same parameter setting as in Section III with N = 160, K = 3, λ = π/4 and g(ω) ∼ N (0, 1).Whereas the mean values are well approximated by the averaging procedure 36 described above, S(t) and C(t) experience significant fluctuations.It is clearly seen that the distribution of these fluctuations is Gaussian.In Figure 8 we show that the variance of S and of C scales approximately as 1/N.This suggests that the scaled mean-subtracted variables are Gaussian processes that are defined entirely in terms of their mean and covariance functions for a and b being either ξ or ζ .Indeed, Figure 9 shows that correlations decay in time with increasing system size N, and that the covariance functions of the scaled variables ξ t and ζ t converge in the limit of N → ∞.Note that the covariance functions of ξ t and ζ t do not scale with N.
The oscillations of the covariance functions for fixed finite N are persistent features which do not decay for increasing length of the time series.To allow for a stochastic description of the fluctuations ξ and ζ we require integrable covariance functions which decay in time.The dominance of the nondecaying large time correlations for small values of N illustrates that the effective stochastic dynamics implied by (21)  and the implied central limit theorem requires a sufficiently large number of oscillators N. We remark that technically one should consider the scaling with respect to the number of rogue oscillator N r rather than with the total number of oscillators N.However, since N r scales linearly with N (cf Figure 5 and equation ( 11)) this notational simplification is justified.

A. Approximation of S(t) and C(t) by a 2-dimensional Ornstein-Uhlenbeck process
To determine an explicit stochastic process that can be used to generate realisations of this Gaussian process, we approximate the Gaussian process Z KS = (ξ t , ζ t ) T by an easily computable two-dimensional Ornstein-Uhlenbeck process dZ = (−ΓZ + ϒZ) dt + Σ dB t (23)   with two-dimensional Brownian motion B t = (B 1 , B 2 ) T .It is well-known that stationary Gauss-Markov processes can be expressed as solutions of Ornstein-Uhlenbeck processes 61,62 .
We consider for simplicity a diagonal matrix Γ = γ I and skewsymmetric rotation matrix ϒ with entries υ 12 = −υ 21 = υ and υ 11 = υ 22 = 0 and a symmetric diffusion matrix Σ with entries σ i j for i, j = 1, 2. The associated covariance function is given by where the Z 0 are random variables drawn from the stationary density of the OU process and the angular brackets ⟨•⟩ OU denote the average over that density.We provide explicit expressions of the covariance function R (OU) (τ) in the Appendix.The parameters of the Ornstein-Uhlenbeck process ( 23) are determined such that its mean is zero and its covariance function (24) best matches the observed covariance function (22) of Z KS = (ξ t , ζ t ) associated with the full deterministic Kuramoto-Sakaguchi model.We perform a nonlinear least-square optimisation minimising the objective function where the sum goes over all entries of the covariance matrix.We choose β = 1, 000 to enforce that the variances of S and C are reproduced.For the parameter values used in Section III with K = 3, λ = π/4, g(ω) ∼ N (0, 1) and N = 160, the fitting yields γ = 0.727 ± 0.0331, υ = 1.739 ± 0.0331 and σ 11 = 0.374 ± 0.0069, σ 12 = σ 21 = 0.0090 ± 0.0011, σ 22 = 0.271 ± 0.0084 with 95% confidence intervals.For different parameter values of the Kuramoto-Sakaguchi model (1), the parameter values of the OU process will be different.Figure 10 shows the four entries of the covariance matrix of Z KS = (ξ t , ζ t ) together with the best fit of the approximating OU process.The fit is reasonable for τ ≲ 2. We remark that one cannot expect that the approximate stochastic reduction can model the dynamics accurately on all time scales.For example, the Kuramoto-Sakaguchi model is deterministic and hence the derivatives of the covariance functions R ξ ξ (τ) and R ζ ζ (τ) are zero at τ = 0 63 (cf. Figure 9b).This feature cannot be reproduced by the stochastic Ornstein-Uhlenbeck process which supports covariance functions that are non-differentiable at τ = 0; this suggest a lower integration bound t min ̸ = 0. We choose here t min = 0.5 and t max = 2.5.Furthermore, the mismatch between the covariance functions obtained form the full deterministic Kuramoto-Sakaguchi model (1) and those corresponding to the best-fit Ornstein-Uhlenbeck process shown in Figure 10 shows that the effective stochastic dynamics of the fluctuations is only approximately a Gaussian process.We have checked that the accuracy of the fit does not significantly improve when increasing the number of oscillators to N = 2, 560.The discrepancy between the observed covariance function and the best-fit Ornstein-Uhlenbeck process reveals that the underlying Gaussian process is more complex and possibly non-Markovian.However, despite the fact that the OU process is only able to coarsely approximate the covariance function of the actual Gaussian process, we will see in Section V that the fitted OU-process serves as an easily computable surrogate process for the interaction term of the rogue oscillators, which is able to reliably reproduce the main statistical features of the system.
B. Reduced stochastic model for the entrained oscillators and associated evolution of the order parameter Expressing the interaction terms S(t) and C(t) in the Kuramoto-Sakaguchi model ( 12) as an OU process with means ⟨S⟩ t and ⟨C⟩ t , we obtain the following closed system of evolution equations for the entrained oscillators θ i with i ∈ C where ξ t and ζ t are the components of a 2-dimensional OU process governed by We can now establish expressions for the order parameters associated with the stochastic reduced system ( 26)- (27).The order parameter r c and the mean phase ψ c can be directly determined from simulations of ( 26)- (27).The full order parameter r can be expressed as Using the definitions for S and C (13)-( 14) which imply we can express (28) in terms of the complex variable Z t = C(t) + i S(t), where we treat Z t now as a complex valued random variable, and obtain We now perform further approximations to formulate a stochastic evolution equation for the order parameter.We separate the mean part and the fluctuations of Z t according to The random variable Z ′ t = ζ t + iξ t is generated from the OU process (27) via the definition (21).Note that this equation for r is only an approximation as we require r(t) ≤ 1 for all times t.
The order parameters r c and r r are both stochastic processes fluctuating around their respective means.However, the variance of the perturbing Ornstein-Uhlenbeck process Z ′ t / √ N, which represents the rogue oscillators, is much larger than that associated with the synchronised cluster.Indeed, we observe numerically a two orders of magnitude larger variance of the rogue oscillators than that of the synchronized oscillators with Var(r c ) ≈ 1.41 × 10 −5 and Var(r r ) ≈ 6.5 × 10 −3 .This suggests that in equation ( 29) for the order parameter r(t) the stochasticity of r c can be neglected to first order against the dominating noise process Z ′ t / √ N, and we can approximate r c ≈ rc = const.We obtain with complex Assuming that the stochastic perturbations are small we can approximate further and expand up to O(1 √ N) to obtain which formally is equivalent to the following evolution equation for the order parameter with The stochastic process Z ′ is not necessarily stationary and we allow for non-equilibrium initial conditions z ′ 0 which are not drawn from the stationary distribution of the Ornstein-Uhlenbeck process, corresponding to random initial conditions for the phases of the rogue oscillators.Note that the order parameter is driven by coloured noise in (33).This is in contrast to Snyder et al 37 who postulated Brownian motion as the driving noise.

C. Dynamic mechanism for the generation of effective stochasticity in the deterministic Kuramoto-Sakaguchi model
The effective stochastic dynamics of S(t) and C(t) are generated by the weak chaoticity of the non-entrained rogue oscillators 64,65 .We show in Figure 11 typical trajectories of rogue oscillators.The dynamics of the weakly chaotic rogue oscillators is characterized by a nearly periodic motion (indeed under the assumption of constant mean-field variables r and ψ the dynamics of each rogue oscillator is approximately governed by the Adler equation ( 3)).Note that the rogue oscillators closest to the synchronized cluster evolve slowly, almost aligned with the synchronized cluster, for long periods FIG.7: Empirical histogram of S(t) and C(t) obtained from a single trajectory of the Kuramoto-Sakaguchi model (1) with N = 160.The continuous curve (online red) shows the corresponding Gaussian best fit.We denote with a filled circle (online orange) the thermodynamic value of ⟨S⟩ t and ⟨C⟩ t as calculated from ( 18) and ( 20), respectively.
interrupted by fast slips.We observed that the times between consecutive phase slips are not constant but exhibit significant variance (not shown).In particular, the variance relative to the mean is largest for the slow rogue oscillators closest to the synchronized cluster and decays monotonically for faster rogue oscillators with higher intrinsic frequencies ω i .The interaction terms ( 13) and ( 14) hence constitute sums of nearly periodic functions of time.It is well known that a sum of many trigonometric functions of linear uncoupled oscillators with uncorrelated random initial conditions approximates a Gaussian processes 38,44,45 .This suggests that the mechanism for the deterministic generation of diffusive behaviour is less a matter of the time-scale separation between faster chaotic rogue oscillators and slower synchronized oscillators (as would be covered by the theory of homogenization), but instead it is given by weak coupling of the synchronized oscillators to (sufficiently) many uncorrelated rogue oscillators.We remark that the randomness stems here entirely from the choice of the initial conditions of the oscillators.For different initial conditions, different realisations of the noise are generated.The weakly chaotic nature, however, allows one to use highly correlated initial conditions, e.g.θ i = const for all i; initially correlated phases will decorrelate after a sufficiently long transient period.This is different to the classical trigonometric approximation which relies on the presence of uncorrelated random initial conditions.

V. COMPARISON OF THE STOCHASTIC REDUCED MODEL WITH THE FULL KURAMOTO-SAKAGUCHI MODEL
We now show that the fluctuations of the order parameter as observed in Figure 2 can be captured by the closed stochastic evolution equation ( 26)-( 27) for the entrained oscillators θ i with i ∈ C .Here we present results again for the parameter setting used in Section III with λ = π/4, g(ω) ∼ N (0, 1) and N = 160.We begin with a coupling strength K = 3 as used in Section III.
The averages ⟨S⟩ t and ⟨C⟩ t in ( 26) can be computed either from long time simulations or from our analytical mean-field theory results (18) and (20).The agreement is up to an er-FIG.8: Scaling of the variance of S and C with varying system size N as calculated from a single trajectory of the Kuramoto-Sakaguchi model (1).The continuous line is a linear best fit indicating a scaling with N −0.977 and N −0.975 for S and C, respectively.
ror of less than 1% (cf. Figure 7).We use in the following the numerically obtained values.The parameters of the Ornstein-Uhlenbeck process (27) were determined by the fitting of the covariance functions according to (25) which we recall here as γ = 0.727, υ = 1.73 and σ 11 = 0.374, σ 12 = σ 21 = 0.0090, σ 22 = 0.271.We remark that in the thermodynamic limit N → ∞ the effect of the stochastic fluctuations ξ t and ζ t in ( 26) decays and we recover the deterministic evolution equation derived in 36 .We simulate the system ( 26)-( 27) using an Euler-Maruyama integrator with a time step of dt = 0.01. Figure 12 shows the empirical histograms of the order parameter r c pertaining to the synchronized oscillators as well as the total order parameter r for all oscillators when simulated using the full Kuramoto-Sakaguchi model (1) and when estimated by the reduced stochastic system (26)-( 27) using ( 15) and ( 29) for the associated order parameters; we remark that results using the approximation (30) for the associated order parameter lead to results indistinguishable by eye.We further show the histogram of the order parameter r c which only takes into account the entrained synchronized oscillators.It is seen that our stochastic reduction captures the observed fluctuations for both r c (t) and r(t) very well.We remark that the histograms are not distinguishable by eye if the order parameter is calculated by an actual time trajectory of r c (t) or by its constant mean rc (cf.( 29) and ( 30)).
Long-time solutions of the reduced stochastic equation for the order parameter (33) generate empirical histograms undistinguishable by eye from those presented in Figure 12.In Figure 13 we show the time evolution of the order parameter r c (t) computed from simulations of the full Kuramoto-Sakaguchi model (cf. Figure 2) and of our reduced stochastic system (26)- (27).It is seen that the qualitative smooth character of the fluctuations observed in the deterministic Kuramoto-Sakaguchi model ( 1) is reproduced by our stochastic system which is driven by coloured OU noise.
To further probe the ability of the reduced stochastic model ( 26)- (27) to capture the collective dynamics of the entrained synchronized oscillators we show results for the fluctuations of the entrained oscillators around the mean phase of the synchronized cluster, θ i − ψ c for i ∈ C .These fluctuations are Gaussian for all entrained oscillators as shown in Figure 14.It is seen that our SDE very well describes entrained oscillators, such as those with indices i = 1 and i = 50, but the degree of approximation is less good for those entrained oscillators at the edge of the cluster with index i = 115.In Figure 15 we show the mean and the variance of θ i − ψ c for all i ∈ C estimated from a long simulation of the full Kuramoto-Sakaguchi model (1) and of the reduced stochastic model ( 26)- (27).It is seen that the mean is very well recovered by the reduced stochastic model.The variance is very well captured for oscillators with an index i that is sufficiently small to ensure that their intrinsic frequency is not too close to that of the closest rogue oscillator.For the simulations we show in Figure 15 the index of the first rogue oscillator is i = 117 for the full Kuramoto-Sakaguchi model.The entrained oscillator with index i = 116, i.e the oscillator on the edge of the synchronized cluster, is fully entrained in the Kuramoto-Sakaguchi model.However, for the approximate stochastic model ( 26)-( 27) the oscillator with index i = 116 experiences rare fast random slips between long periods of entrainment (of the order of 1, 000 time units on average).The rare and fast slips do not affect the mean which is still close to the mean corresponding to the Kuramoto-Sakaguchi model, but the variance is much higher with a value of 0.012 (not shown in Figure 15 where we only show oscillators with index i ≤ 115).
We have so far reported only on a single value K = 3 of the coupling strength.It is pertinent to mention that our approach is not restricted to any value of the coupling strength but only to the number of rogue oscillators present.For each value of the coupling strength, however, the corresponding parameters of the approximating Ornstein-Uhlenbeck process will differ and will need to be determined.For K = 3 we found that the number of rogue oscillators N r = 43 out of a total of N = 160 oscillators was sufficient for the fluctuations to be described as a Gaussian process.The number of rogue oscillators decreases with increasing coupling strength, see Figure 6.For an increased coupling strength with K = 7 there are only N r = 2 rogue oscillators for N = 160 which is insufficient for the Gaussian process approximation.However, increasing the total number of oscillators to N = 2, 560 increases the number of rogue oscillators to N r = 46, and we observe again that the statistics of the order parameter is well recovered by the Gaussian process approximation as shown in Figure 16.Conversely, for coupling strengths below the critical coupling strength K c one can use the Gaussian approximation for any of the N oscillators with the remaining N r = N − 1 oscillators contributing an additive noise for the single oscillator with index i = i ⋆ and ψ c = θ i ⋆ .Figure 16 shows that indeed the order parameter r is well approximated by the Gaussian approximation for the sub-critical case K = 1.25.
The numerical results reported above all used the equiprobable draws of the intrinsic frequencies described in Section III.For random draws the number of rogue oscillators N r fluctuates across realisations, and the statistical features of S and C will vary.We checked that with random draws of the intrinsic frequencies the distributions of S and C are still near-Gaussian and S and C can be effectively approximated by a Gaussian process for most realisations.For random intrinsic frequency realisations which give rise to the formation of smaller interacting clusters, the effective dynamics would need to involve multiple clusters and becomes much more involved.The fluctuations decrease with increasing number of oscillators N approaching the results of the equiprobable draws presented here.

VI. DISCUSSION
We considered the effect of non-entrained rogue oscillators on the collective behaviour of the entrained synchronized oscillators in the Kuramoto-Sakaguchi model.We established an average effect, in the spirit of a law of large numbers, and then considered fluctuations around the mean effect FIG. 13: Temporal evolution of the order parameter r c (t).
Shown are results obtained from the full Kuramoto-Sakaguchi model (1) (online blue), the reduced stochastic model ( 26)- (27) driven by an OU process (online red).All other parameters are the same as in Fig. 1.  26)- (27).All other parameters are the same as in Fig. 1.
akin to the central limit theorem.We presented results from numerical simulations that suggest that the fluctuations can be approximated by a Gaussian process with a stationary density and fluctuations which decay as 1/ √ N. Hence for fluctuations to have a significant effect one needs the number of rogue oscillators to be sufficiently large to allow for a central limit theorem and sufficiently small to have a nonnegligible variance.Gaussian processes are determined entirely by their mean and their covariance function.We used a two-dimensional Ornstein-Uhlenbeck process as a surrogate stochastic Gaussian process by fitting the covariance function.This enabled us to replace the interaction term involving the rogue oscillators by this Ornstein Uhlenbeck  26)- (27).All other parameters are the same as in Fig. 1.
process, and to formulate a closed equation for the entrained synchronized oscillators only, which are driven by this complex Ornstein-Uhlenbeck process.The reduced system of stochastic differential equations showed remarkable capability to reproduce the statistical behaviour of the entrained oscillators such as the probability density function of the order parameter.It is pertinent to mention that the order parameter is driven by coloured Ornstein-Uhlenbeck noise and not, as previously claimed on heuristic grounds, by Brownian motion.The OU noise was fitted by directly using information of the unresolved dynamics (i.e. the statistical behaviour of S and C) rather than by blind fitting of the resolved synchronized behaviour, as is usually done when imposing coarse grained model.
Our numerical experiments were restricted to normally distributed intrinsic frequencies.There are in principle no restrictions to the nature of the frequency distribution g(ω).For unimodal frequency distributions and also for uniform frequency distributions with λ ̸ = 0 there is a range  26) is calculated using (29).
in coupling strengths for which the dynamics is characterized by the interaction of a single partially synchronized cluster with a collection of non-entrained rogue oscillators as discussed here.It would be interesting to see how finite-size induced noise can effect the interaction of several partially synchronized clusters.This would be possible for multimodal intrinsic frequency distributions.
In future work it would be interesting investigate how far the model reduction is able to characterize the fluctuations of the order parameter the phase transition as seen in Figure 6b.Since the parameters the best-fit OU process typically differ for different values of the coupling strength K this may prove to be numerically costly.
On a theoretical level, it would be interesting to show the convergence to an effective stochastic equation for the synchronized oscillators more rigorously and/or find explicit expressions for the parameters of the reduced SDE.This would require to extend the method of trigonometric approximations to an approximation of Gaussian processes by a sum of solutions of the Adler equation (3).This could then allow for a computationally feasible investigation of the fluctuations of the order parameter at the phase transition, mentioned above.
The stochastic reduced equation can now be further reduced, for example via the method of collective coordinates [33][34][35][36]65 . In articular, one can employ the methods developed for the reduction of a stochastic Kuramoto model 56 to the proposed stochastic system here for the Kuramoto-Sakaguchi model.This has the potential to capture finite size effects which cannot be captured via standard mean-field theories, in particular, the diffusivity of the mean phase ψ and ψ c .
where the first expectation is taken with respect to Brownian motion paths and initial conditions and the second expectation with respect to initial conditions, the covariance function is obtained as R(τ) = lim t→∞ R(t,t + τ) = lim t→∞ t 0 e L(t−u) ΣΣ T e L T (t+τ−u) du.

FIG. 1 :
FIG. 1: Snapshot of the phases θ i (a) and effective frequencies ωi (b) for the Kuramoto-Sakaguchi model (1) with N = 160 oscillators and phase frustration λ = π 4 at coupling strength K = 3 with a Gaussian intrinsic frequency distribution with zero mean and unit variance.Circles (online blue) denote oscillators entrained in collective synchronised dynamics.Crosses (online red) denote the non-entrained rogue oscillators.We identify oscillators with indices i ≤ 116 as the synchronized oscillators with N c = 116 and oscillators with indices i > 116 as the rogue oscillators with N r = 43.

1 FIG. 2 :
Figure6shows the well-known transition from incoherence to synchronisation for increasing coupling strength K for a

FIG. 4 :
FIG. 4: (a): Empirical distribution of the order parameter r(t) obtained from a single simulation of the Kuramoto-Sakaguchi model (1) at different system sizes N. (b): Mean value of r(t) for different system sizes N. The dashed line indicates the value at the corresponding thermodynamic limit, r∞ .(c): Variance of r(t) for different system sizes N. The red line represents the best-fit suggesting a scaling law of N −0.976 .All other parameters are the same as in Fig. 1.

FIG. 6 :
FIG. 6: Transition from incoherence to synchronisation of the Kuramoto-Sakaguchi model (1) with increasing coupling strength K with N = 160 for a Gaussian intrinsic frequency distribution g(ω) with mean zero and variance 1. (a): Order parameter r. (b): Variance of the order parameter, quantifying the size of fluctuations.(c): Number of entrained oscillators N c , which are part of the collective synchronised cluster (circles, online blue) and number of non-entrained rogue oscillators N r (crosses, online red).Shown are relative numbers.The vertical lines mark the onset of partial synchronisation at K c = 1.93 and the onset of global synchronisation at K g = 8.34.

FIG. 9 :FIG. 10 :
FIG. 9: Entry R ξ ξ (τ) of the covariance function of Z KS = (ξ (t), ζ (t)) obtained from a single trajectory of the Kuramoto-Sakaguchi model (1) for different system sizes N. We show the covariance function for a large range τ ≤ 40 as well as a close up for τ ≤ 1.5.

FIG. 12 :
FIG. 12: Comparison of the empirical histograms for r c (left) and r (right) obtained from a single trajectory of the Kuramoto-Sakaguchi model (1) for N = 160 and from the dynamics of the reduced stochastic model (26) driven by an OU process.All other parameters are the same as in Fig. 1.

FIG. 14 :
FIG. 14: Empirical density of θ i − ψ c for the entrained oscillators with indices i = 1, i = 50 and i = 115.Shown are results for the full Kuramoto-Sakaguchi model (1) for N = 160 and for the reduced stochastic model (26)-(27).All other parameters are the same as in Fig.1.

FIG. 15 :
FIG. 15: Mean (left) and variance (right)) of the fluctuations of the entrained synchronized oscillators around their mean phase, θ i − ψ c for all i ∈ C .Shown are results for the full Kuramoto-Sakaguchi model (1) for N = 160 and for the reduced stochastic model (26)-(27).All other parameters are the same as in Fig. 1.