Structure preserving schemes for the continuum Kuramoto model: phase transitions

The construction of numerical schemes for the Kuramoto model is challenging due to the structural properties of the system which are essential in order to capture the correct physical behavior, like the description of stationary states and phase transitions. Additional difficulties are represented by the high dimensionality of the problem in presence of multiple frequencies. In this paper, we develop numerical methods which are capable to preserve these structural properties of the Kuramoto equation in the presence of diffusion and to solve efficiently the multiple frequencies case. The novel schemes are then used to numerically investigate the phase transitions in the case of identical and non identical oscillators.


Introduction
Synchronization phenomena of large populations of weakly coupled oscillators are very common in natural systems, and it has been extensively studied in various scientific communities such as physics, biology, sociology, etc. [2,4,29,52]. Synchronization arises due to the adjustment of rhythms of self-sustained periodic oscillators weakly connected [2,47], and its rigorous mathematical treatment is pioneered by Winfree [53] and Kuramoto [4]. In [4,53], phase models for large weakly coupled oscillator systems were introduced, and the synchronized behavior of complex biological systems was shown to emerge from the competing mechanisms of intrinsic randomness and sinusoidal couplings. Since then, the Kuramoto model becomes a prototype model for synchronization phenomena and various extensions have been extensively explored in various scientific communities such as applied mathematics, engineering, control theory, physics, neuroscience, biology, and so on [28,47,52].
Given an ensemble of sinusoidally coupled nonlinear oscillators, which can be visualized as active rotors on the unit circle S 1 , let z j = e iϑ j be the position of the j-th rotor. Then, the dynamics of z j is completely determined by that of its phase ϑ j . Let us denote the phase and frequency of the j-th oscillator by ϑ j and θ j , respectively. Then, the phase dynamics of Kuramoto oscillators are governed by the following first-order ODE system [4]: sin(ϑ i − ϑ j ), i = 1, · · · , N, t > 0, (1.1) subject to initial data ϑ i (0) =: ϑ i0 , i = 1, · · · , N, where K is the uniform positive coupling strength, and ω i denotes the natural phase-velocity (frequency) and is assumed to be a random variable extracted from the given distribution g = g(ω) satisfying R g(ω) dω = 1.
We notice that the first term and second term in the right hand side of the equation (1.1) represent the intrinsic randomness and the nonlinear attraction-repulsion coupling, respectively. The system (1.1) has been extensively studied, and it still remains a popular subject in nonlinear dynamics and statistical physics. We refer the reader to [2] and references therein for general survey of the Kuramoto model and its variants. In [4], Kuramoto first observed a continuous phase transition in the continuum Kuramoto model (N → ∞, see (1.2) below) with a symmetric distribution function g(ω) by introducing an order parameter r ∞ which measures the degree of the phase synchronization in the mean-field limit. More precisely, the order parameter is given by In fact, Kuramoto showed that the order parameter r ∞ as a function of the coupling strength K changes from zero (disordered state) to a non-zero value (ordered state) when the coupling strength K exceeds a critical value K c := 2/(π g(0)), i.e., r ∞ (K ) = 0 (incoherent state) for K ∈ [0, K c ), 1 ≥ r ∞ (K ) > 0 (coherent state) for K > K c , and r ∞ (K ) increases with K .
Later, it is also observed that the continuum Kuramoto model can exhibit continuous or discontinuous phase transition by taking into account different types of natural frequency distribution functions [11,12,45]. Here, continuous phase transitions refer to the continuity at K = K c of the order parameter r ∞ (K ). It is known to be continuous for the Gaussian distributed in frequency oscillators [52] and discontinuous for the uniformly distributed in frequency oscillators [45].
A very relevant issue from the application viewpoint is how stable these stationary states and phase transitions are by adding noise to the system [4,50]. The mean-field limit equation associated to (1.1) with standard Gaussian noise of strength 5) subject to the initial data (1.3). Phase transitions have also been found in terms of the order parameter r ∞ as a function of K for D > 0. In fact, it was proven in [50] that for identical oscillators there is a continuous phase transition. This phase transition is also continuous for Gaussian and uniformly coupled noisy oscillators [2] and discontinuous for bimodal distributions [16,25]. Stability of the coherent and incoherent states and the asymptotic dynamics for the Kuramoto model with noise (1.5) were analyzed in [31,32,37,38]. Noise-driven phase transitions in interacting particle systems are a hot topic of research from the numerical and theoretical viewpoints, see [5][6][7][8][9][10]25,33] and the references therein.
In the present work, we focus on the construction of effective numerical schemes for the continuum Kuramoto system with diffusion obtained in the limit of a very large number of oscillators and in presence of noise. The numerical solution of such system is challenging due to the high dimensionality of the problem and the intrinsic structural properties of the system which are essential in order to capture the correct physical behavior. The literature dealing with the study of numerical schemes to the continuum Kuramoto model (1.2) is mainly focused on deterministic spectral and stochastic methods. We refer to [2, Section VI. B] and the references therein for general discussions on the stochastic and deterministic numerical methods used in Kuramoto models, see also [10,35,36,47] for Monte Carlo related methods and [1,3,46] for spectral methods.
Structure preserving numerical schemes for mean-field type kinetic equations have been previously constructed in [2,14,17,18,20,22,44,51]. We refer also to [34] for related approaches for general systems of balance laws and to [27] for an introduction to numerical methods for kinetic equations. Here, following the approach in [18,22,44] we introduce a discretization of the phase variable such that nonnegativity of the solution, physical conservations, asymptotic behavior and free energy dissipation are preserved at a discrete level for identical oscillators. This approach is then coupled with a suitable collocation method for the frequency variable based on orthogonal polynomials with respect to the given frequency distribution. Similar collocation strategies have been used for kinetic equations in the case of linear transport and semiconductors [39,41,49].
The rest of the manuscript is organized as follows. First, in Section 2 we recall some basic properties of the continuum Kuramoto model, in particular concerning some relevant existence theory and some useful estimates underpinning our numerical strategy. We also discuss the asymptotic behavior of the Kuramoto model with noise, stationary states and related free energy estimates in the case of identical oscillators. The discretization of Kuramoto systems is discussed in Section 3. Structure preserving schemes are developed in combination with a collocation method for the frequency variable. The asymptotic behavior of the schemes as well as the other relevant physical properties are discussed. In Section 4 we present several numerical tests, with a particular emphasis on the study of phase transitions. We will also compare our structure preserving methods to Monte Carlo and spectral methods. Future research directions and conclusions are reported in the last Section.

Stability of the Kuramoto model with respect to the frequency distribution
In this subsection, we briefly review some relevant issues related to the well-posedness theory and several useful estimates for the continuum Kuramoto equation (1.2) and its version with noise (1.5). We refer to [19,21,40] for further details. We first provide a priori estimates for the equation (1.2).
Proof. The proof of (i) is straightforward. For the identity (ii), by using the definition of r in (1.4) we rewrite u[ρ] as We next discuss global existence and uniqueness of measure-valued solutions. Let us denote by M(T × R) the set of nonnegative Radon measures on T × R, which can be regarded as nonnegative bounded linear functionals on C 0 (T × R).
Then our notion of measure-valued solutions to (1.2) is defined as follows.
) be a measure-valued solution to (1.2) with an initial measure μ 0 ∈ M(T × R) if and only if μ satisfies the following conditions: • μ is weakly continuous, i.e., μ t , h is a continuous as a function of t, ∀h ∈ C 0 (T × R), We also introduce the definition of the bounded Lipschitz distance. Definition 2.2. Let μ, ν ∈ M(T × R) be two Radon measures. Then the bounded Lipschitz distance d(μ, ν) between μ and ν is given by where the admissible set S of test functions is defined as We now present the results on the global existence and stability of measure-valued solutions to (1.2). We refer the reader to [19,21,40] for details of the proof.
Remark 2.1. One can also characterize the solution μ of the initial-value problem for (1.2) as the push-forward of the initial data μ 0 through the flow map generated by u[μ], i.e., for any h ∈ C 1 Note that the characteristic system (2.2) is well defined since the velocity field u[μ] is globally Lipschitz in (ϑ, ω) and continuous on t even for general measures μ ∈ C w (0, T ; M(T × R)).

Remark 2.2.
As a simple extension of Proposition 2.1, see [21], we also obtain that the solution μ t can be approximated as a sum of Dirac measures of the following form This can be accomplished by choosing a particle approximation of the initial data. One procedure to construct particle approximations is the following. Let us choose the smallest square containing the support of the initial data μ 0 , i.e., supp(μ 0 ) ⊂ R. For a given n, we divide the square R into n 2 subsquares R i , i.e., Let ϑ 0 i , ω i be the center of R i for i = 1, · · · , n 2 . Then we construct the initial approximation μ n 0 as Then we can easily show that d(μ n 0 , μ 0 ) ≤ C /n → 0 as n → ∞.
From now on, we assume that the initial measure is a smooth absolutely continuous with respect to Lebesgue density with connected support. This assumption can be removed by the stability property in Proposition 2.1. We next show stability of the order parameter r with respect to approximated frequency distributions and initial data by again using the stability estimate presented in Proposition 2.1.  Proof. We first estimate the 1-Wasserstein distance between ρ n and μ. Observe that ρ n (t) ∈ S for all t ≥ 0 by the pushforward characterization in Remark 2.1. Moreover, it also follows from Remark 2.1 that with n (s; s, ϑ, ω) = ϑ for all n ≥ 1 and satisfies (2.2). For notational simplicity, we denote n (t; 0, ϑ, ω) = n (t) and (t; 0, ϑ, ω) = (t). Then a straightforward computation yields where C > 0 is independent of n. This together with the boundedness of n and where C > 0 is independent of n. We now show the strong convergence of r n to r. For this, we use the following identities to obtain |r n (t) − r(t)| ≤ T×R cos ϑ ρ n t g n − μ t g (dϑ, dω) + T×R sin ϑ ρ n t g n − μ t g (dϑ, dω) Here I 1 can be estimated as This completes the proof. 2

Remark 2.3.
A similar result to Proposition 2.2 can be obtained if we approximate the initial data μ 0 ∈ M(T × R) and the frequency distribution g n (ω) ∈ L ∞ (R) as follows: μ n

Remark 2.4.
A similar result to Proposition 2.2 can also be proved for the Kuramoto model with noise (1.5). The strategy of the proof is analogous but it uses stochastic processes techniques to write the corresponding stochastic differential equation systems. Moreover, one needs to resort to Wasserstein distances instead of the bounded Lipschitz distance above to make the stability argument of the solutions. We do not include the details of the proof since the technicalities lie outside of the scope of the present work, see [15] for related problems.
Notice that Proposition 2.2 and Remark 2.4 allow us to work with continuous frequency distributions for both Kuramoto models (1.2) and (1.5) by approximation. More precisely, we can approximate Gaussian and uniform frequency distributions by sums of Dirac Deltas at a finite number of frequencies while approximating the initial data by smooth functions if they are not regular enough. This fact will be used in the numerical schemes in Section 3 and the simulations in Section 4.

The Kuramoto model with diffusion: stationary states and free energies
We first derive an explicit compatibility condition for smooth stationary states ρ ∞ of the equation (1.5). We first easily find from (1.5) and (2.1) that where r ∞ and ϕ ∞ are given by We notice that we can set ϕ ∞ = 0 in the stationary state without loss of generality by choosing the right angular reference system. This yields for some function c 0 (w). Solving the above differential equation, we find On the other hand, since ρ ∞ (2π , ω) = ρ ∞ (0, ω), we get and subsequently, this implies Let us discuss further properties in the case of noisy identical Kuramoto oscillators, which are governed by the equation Then we can set ω = 0 without loss of generality in (2.3) to conclude that its stationary state is given by Note that for both identical and non-identical oscillators if r ∞ = 0, i.e., incoherence state, we obtain , for all ω ∈ R, due to the normalization. In this case the equation (1.5) has the structure of a gradient flow. More precisely, if we set with W (ϑ) = cos ϑ , then it is easy to check that Using this observation, we now estimate the following free energy where the dissipation rate D E (t) is given by We next provide the monotonicity of the order parameter for the identical case.
In particular, if the strength of noise D is strong enough such that D ≥ K , then ṙ ≤ 0 for all t ≥ 0.
Proof. By definition of the order parameter r, we geṫ where I 1 vanishes since For the estimate of I 2 , we find Combining the above two estimates concludes the desired result. 2 Before passing to the construction of numerical methods the following remark should be made.

Remark 2.5.
Several works for the continuum Kuramoto model are based on the g-weighted kinetic density f (ϑ, ω, t) := ρ(ϑ, ω, t)g(ω) (see [13] for example). In this way, one can rewrite the Kuramoto model (1.5) as Note that the distribution of the natural frequencies is now given by Even if, in the continuation, we will use the form (1.5) for the construction of our structure preserving numerical methods, they can be easily reformulated to the form (2.4). In Section 4 we will show both ρ(ϑ, ω, t) and f (ϑ, ω, t) in some numerical tests.

Structure preserving methods
The goal now is to propose finite volume numerical schemes preserving the structure of gradient flow to the case of identical oscillators and that are generalizable for oscillators with natural frequencies given by a distribution function g(w).
To start with, from (2.1) and (1.5), it follows that the density ρ satisfies the following continuity equation 3.1. Semi-discrete structure preserving schemes Inspired by [17,18,22,44,51], we construct a discrete numerical scheme in the variable ϑ for the above equation as follows. For i = 1, · · · , N, we first introduce a uniform spatial grid ϑ i ∈ T such that ϑ i+1 − ϑ i = ϑ and ϑ N+k = ϑ k for k ∈ R. Without loss of generality, we set ϑ 1/2 = ϑ N+1/2 = 0 ≡ 2π , and we then define We consider the following approximations for (3.1) As in [22,44] a suitable choice of the weight functions δ i+1/2 yields a method that maintains nonnegativity of the solution (without restrictions on ϑ ) and preserves the steady state of the system with arbitrary order of accuracy. We will refer to the schemes obtained in this way as Chang-Cooper type schemes [22]. First, observe that when the numerical flux (3.4) vanishes we get Similarly, if we consider the exact flux (3.2), by imposing F [ρ] ≡ 0, we have D∂ ϑ ρ = uρ.
Integrating the above equation on the cell Therefore, by equating (3.6) and (3.7) we recover We can state the following. Proof. The latter result follows from the simple inequality exp(x) ≥ 1 + x. 2 On the other hand, the above expression for u i+1/2 is not that useful in practice since it contains r which depends on ρ.
Thus, it would be more technically useful to write where we set ρ s Remark 3.2. The resulting scheme is second order accurate in ϑ for D > 0, and degenerate to simple first order upwinding in the limit case D = 0. In fact, it is immediate to show that as D → 0 we obtain the weights In the lemma below, we show that the numerical scheme conserves the mass.
Proof. It follows from the periodicity of domain that ρ N+1 = ρ 1 and ρ N = ρ 0 . (3.10) Similarly, we get u N+1/2 = u 1/2 and subsequently this implies ξ N+1/2 = ξ 1/2 , δ N+1/2 = δ 1/2 , and ρ N+1/2 =ρ 1/2 . From the above properties, we can easily obtain This, together with the following straightforward computation gives the desired result. 2 We next provide the positivity preservation whose proof can be obtained by using almost same argument as in [44,Proposition 1]. However, for the completeness of this work, we sketch the proof in the proposition below. For this, we introduce the time discretization t n = n t with t > 0 and n ∈ N 0 and consider the following forward Euler method Then the explicit scheme (3.11) preserves nonnegativity, i.e., ρ n+1 Proof. It follows from (3.11) that On the other hand, we easily find We also notice that max 0≤i≤N This, together with the property of convex combination, yields that the nonnegativity is preserved if the time step t satisfies The parabolic stability restriction t = O( ϑ 2 ) which appears in Proposition 3.2 can be avoided if we use a semi-implicit scheme as in [44]. To be more precise, let us consider It follows from (3.12) that Thus, the matrix M n is invertible and we havē Then, by the same arguments as in [44], if the time step t satisfies t < ϑ 2C 0 , (3.14) where C 0 > 0 appeared in Proposition 3.2, then the semi-implicit scheme (3.12) preserves the nonnegativity. It is also clear that the scheme (3.12) conserves the mass. Finally, concerning the computational cost, it is worth to mention that the usual Thomas algorithm for the solution of a tridiagonal system in O (N) operations can be extended to cyclic tridiagonal matrices of the form (3.13) thanks to the Sherman-Morrison formula. We refer to [48], Section 2.7.2, for more details.

Discrete free energy dissipation
Next, we present the discrete free energy estimate for the case of identical oscillators. Set Taking the time derivative to E (t) gives d dt (3.10). This together with (3.3) yields d dt On the other hand, it follows from Remark 3.1 that We then combine (3.15) and (3.16) to find d dt We notice that if ρ i+1/2 is given through the following relation we have the discrete dissipation of free energy estimate, however, this is not the case for On the other hand, if we choose [18,44] ρ E (3.18) then ρ E i+1/2 satisfies the relation (3.17), and subsequently this yields d dt where the discrete dissipation rate D (t) is given by (3.19)

Proposition 3.3. Let us consider in (3.3) the numerical flux function
with D (t) given by (3.19).
Remark 3.4. The weights defined by (3.18) are such that 0 < δ E i+1/2 < 1, moreover, it is easy to verify that the numerical flux function (3.20) vanishes when the corresponding flux (3.2) is equal to zero over the cell [ϑ i , ϑ i+1 ]. On the other hand, nonnegativity of the solution, is satisfied only under more restrictive conditions. In fact, similar to central differences, we have a restriction on the mesh size ϑ which becomes prohibitive for small values of the diffusion function D. We refer to [44] for further details.

Frequency discrete schemes
In order to derive fully discrete schemes we must discuss the discretization of the frequency variable ω ∈ R. Since, the computation of the order parameter and the flux function are obtained as integrals in ω with respect to the probability density g(ω) it is natural to consider Gaussian quadrature formulas based on orthogonal polynomials. This technique, which corresponds to a collocation method, was previously used for linear transport equation and semi-conductors models and is closely related to moment methods [27,39,41,49]. Moreover, from the perspective of Proposition 2.2, we are approximating g by a suitable weighted average of Dirac Deltas at chosen frequencies as specified below.
Here, for reader's convenience, we recall some basic facts concerning orthogonal polynomials and Gaussian quadrature.
Let us consider an orthogonal system of polynomials H n (ω) n∈N , with respect to the probability density function g(ω) where δ mn is the Kronecker delta and c n are normalization constants Of particular interest for applications to the Kuramoto system are the Legendre polynomials in the case of a uniform distribution of frequencies and the Hermite polynomials for normally distributed frequencies (see for example [54]). Now let ω k , k = 1, . . . , M, be the zeros of H M and let l k (·) be the (M − 1)th degree Lagrange polynomial through the nodes ω k , k = 1, . . . , M. Therefore, the integrals ρ s j (t) and ρ c j (t) in (3.9) can be approximated taking (3.21) where the weights g k are given by It is well-known that (3.21) becomes an equality if ρ j (ω, t) is any polynomial of degree less than or equal to 2M − 1 in R. Notice, as mentioned earlier, that for general ρ j (ω, t), these formulas are equivalent to say that we approximate our frequency distribution by The resulting numerical scheme for the Kuramoto model is therefore a collocation method of the form

Numerical examples
In this section, we present several numerical experiments showing phase transitions of the order parameter r ∞ , which is defined by for the continuum Kuramoto equation (1.2) solved with the structure preserving introduced in Section 3. It is worth mentioning that the semi-implicit scheme described in Remark 3.3 requires the time step t = O( ϑ), see (3.14) whereas the parabolic restriction t = O( ϑ 2 ) is needed for the explicit scheme. For that reason, we employ the semi-implicit scheme to investigate the large time behavior of solutions numerically. We refer to this method as Implicit Structure Preserving (ISP) scheme. We only use the explicit scheme presented in Proposition 3.2, hereafter denoted as ESP scheme, for the time dependent solutions of ρ in Figs. 1, 7, and 10 since it provides second order accuracy in both t and ϑ . Of course, the stationary behavior is computed exactly by both methods. For comparison purposes we compute a direct particle simulation of the noisy Kuramoto system where the diffusion part is solved using a standard random Brownian motion [43]. We refer to this method as Particle Monte Carlo (PMC) scheme. For identical oscillators, to emphasize some of the properties of the new schemes, we compute also the solution with the spectral solver developed in [2] which essentially consist in a standard Fourier-Galerkin approximation to the continuum Kuramoto equation (1.2). We refer to this method as Fourier-Galerkin Spectral (FGS) method (see [2, Section VI. B] for more details). In all the numerical results, thanks to their favorable stability properties, we used the Chang-Cooper type fluxes defined in (3.4) and (3.8). Analogous results are obtained using the fluxes in (3.20) and (3.18). In all graphs in Figs. 1, 7, 10 and 13, the black curves are the reference solutions obtained with N = 510 grid points.

Identical Kuramoto oscillators
First, to test the validity of our method, we present numerical simulations for the identical Kuramoto oscillators with diffusion.
It is known that the identical Kuramoto equation (1.5) with g = δ 0 exhibits a phase transition at K = K c := 2D, i.e., the oscillators become uniformly distributed which implies r ∞ ≡ 0 for K ≤ K c , on the other hand, for K > K c the oscillators converge to a phase-locked state which have a positive order parameter r ∞ > 0, see [2,31,32] for more discussion. In order to observe the phase transition, we consider as initial data a symmetric sum of two Gaussians: where c 0 > 0 is fixed by the normalization, i.e., This yields that if we fix the strength of noise D and varies the coupling strength K , then we get the opposite behavior of solutions.
In Fig. 1    We also compare the time evolution of solutions ρ(t) with the two different numerical fluxes, (3.8) and (3.18), presented in Section 3.2 in Fig. 2. As expected the two fluxes are in good agreement and provide essentially the same result. In the same figure, we show the time evolution of free energies with the two different fluxes. To emphasize the differences we used a very rough mesh with N = 11. They share the nonincreasing property of the free energy and converge to the same value as the mesh is refined. Next, to emphasize the structure preserving properties of the new schemes, in Fig. 3, we report the same kind of computations for decreasing values of the diffusion coefficient D together with the results obtained with the Fourier-Galerkin spectral method. Clearly, the FGS method is not capable to compute the solution for small values of D unless the grid size resolves the size of the diffusion coefficient and we can observe the formation of oscillations which produce negative values of the solution. On the contrary the ISP method is robust even for small values of D and provides nonnegative solutions under the stability condition (3.14).
In Table 1 we report the discrete L 1 error at the steady state for the ISP scheme and the FGS method for K = 1 and D = 0.1. Note that 64 modes are required to the spectral method to match the accuracy of the structure preserving scheme at the steady state. Here the reference solution was computed with the spectral method using N = 4096 grid points. Clearly, for smaller values of D a larger number of modes is required by the FGS method to have a comparable accuracy with the ISP scheme.
The robustness of the ISP method for small values of D is also shown in Fig. 4 (left), where the diffusion coefficients range from 0 to 1. The figure shows that there is a phase transition around K c = D = 1/2 from coherent to incoherent states, see [2,32] for detailed discussion on the critical value K c = D = 1/2. We next fix the strength of noise D = 1 and varies the coupling strength K , and investigate the behavior of r ∞ (K ). As expected, the phase transition occurs at K c = 2D = 2, see Fig. 4 (right). In the pictures we took a mesh spacing of 0.01 for D and K , in both cases the structure preserving schemes are capable to capture very well the phase transition. As a comparison, we also report the results obtained with the PMC method on a mesh spacing of 0.1 in D and K , where N p = 5 × 10 4 particles and 1000 averages at the steady state are taken. It is clear, that a far superior accuracy can be obtained with the deterministic approach. Let us mention that, even if a careful comparison of the computational cost of the two approaches was outside of the scopes of the present manuscript (we did not try to optimize the coding), for the given parameters the overall computational cost of the PMC method was considerably larger than the ISP method. Beside the different cost of PMC due to the number of particles versus grid points: O (N p ) for PMC and O (N) for ISP for a single value of K ; we also observed that a smaller time step has to be used at the Monte Carlo level in order to achieve the desired accuracy.
Finally, let us comment that in order to speed up the computations with the deterministic scheme, we proceed by continuation forward or backward decreasing the value of D or K respectively and taking as initial seed the steady state of the previous computation. This numerical strategy is done in all reported cases of phase transitions below. Using as transition parameter K , when continuous phase transitions are expected, as it is the case for identical oscillators, our backward iterations were faster and more accurate to give the right transitions values. However, as we will see in Section 4.2.3, for discontinuous phase transitions, our forward iterations on K by continuation were more accurate. This procedure allows us to continue the bifurcation branches even if their stability basins are becoming very small near the critical order parameters K .
Lastly, we numerically investigate the long time behavior of the order parameter with different values of K in Fig. 5 showing that the convergence to steady state gets slower and slower around the critical value K = 2. This critical slowdown of the convergence near the phase transition critical value is a well-known phenomenon, see [26] for instance. It is worth mentioning that in this case using an implicit method is of paramount importance since the explicit CFL condition is too restrictive for such a long time computation close to the critical value.

Nonidentical Kuramoto oscillators
Next, we numerically study the challenging case of the phase transition of the order parameter for nonidentical oscillators. For the initial data, similarly as before, we set

Uniform distribution
In this test case, as a function g(ω) for the natural frequencies, we consider the following uniform distribution centered about 0 with variance σ g : We discretize the velocity field using the nodes of the Gauss-Legendre quadrature approximation. In Fig. 6, we again take D = 0.5, N = 200, M = 10 to investigate the phase transition of the order parameter r ∞ (K ) numerically. Note that in this case, the critical coupling strength K c (D) can be computed from the formula [50] K c (D) = 2 ⎛ ⎝ R g(Dω) The expression above is derived under the symmetry assumption on the distribution function g around its mean value.   √ K − K c near and above the critical coupling strength, see [50] for instance. In order to confirm that numerically, in the small figure in Fig. 6, we compare between r ∞ (K ) and c √ K − K c , where the constant c is computed by using the built-in polyfit MATLAB command. As a comparison, we also report the results obtained with the PMC method on a mesh spacing of 0.1 in K , where N p = 5 × 10 5 particles and 10000 averages at the steady state are taken. Even in the case of nonidentical oscillators, the deterministic ISP scheme, which has a cost of O (MN) against O (N p ) for a single value of K , was far more efficient than the corresponding PMC. We also investigate the time evolution of the solution in Fig. 7 (left) where we reported at various times the behavior of the average value of the solution ρ(ϑ, t) := R ρ(ϑ, ω, t)g(ω) dω, for N = 51, M = 10, K = 1, and D = 0.5. The evolution of the averaged steady state ρ ∞ with D = k/10, k = 1, · · · , 4 is demonstrated in Fig. 7 (right).
In Fig. 8, we show the steady state for the density ρ ∞ (ϑ, ω) for D = 0.5, N = 100, and M = 30 with different values of K . We choose K = 0.5, which is in the region of subcritical coupling strength and K = 2 in the region of supercritical one.

Gaussian distribution
Next, we choose the function g(ω) for the natural frequencies as the following Gaussian distribution centered in 0 with variance σ g : We compute the definite integral in the velocity field (1.2), see also (3.5), using the Gauss-Hermite quadrature approximation, therefore the M grid points for the natural frequencies corresponds to the nodes of such quadrature formula. In , where the phase transition of the order parameter is expected to happen. The comparison between r ∞ (K ) and c √ K − K c is shown in the zoomed image in Fig. 9. Again, we also report the results obtained with the PMC method on a mesh spacing of 0.1 in K , where 5 × 10 5 particles and 10000 averages at the steady state are taken. Next, in Fig. 10 we plot the time evolution of ρ(t) and the evolution ρ ∞ (D) with different values of D.
Finally, in Fig. 11, we show the steady state ρ ∞ (ϑ, ω) and f ∞ (ϑ, ω), which is the weighted kinetic density introduced in (2.4), for K = 2, D = 0.5, N = 100, and M = 30, where the coupling strength K = 2 lies in the region of supercritical coupling strength. Clearly, if the coupling strength K lies in the subcritical region, we have the constant state ρ ∞ as in the steady state of Fig. 8 (left), so we omit the result.

Bimodal distribution
In the last example, we consider a double Gaussian (bimodal) distribution function g(ω) for the natural frequencies: Note that, in this case, it is observed that a discontinuous phase transition of the order parameter occurs, see [16,25].
Similarly as before, we investigate this discontinuous phase transition r ∞ (K ) with N = 200, M = 10, and D = 0.5 in Fig. 12.
We took different mesh spacing in K , we use a step of 0.005 in the interval [1.5, 2], a step of 0.0005 in the zoomed image in Fig. 12. Again, we also put the PMC method result with a mesh spacing of 0.05 in the interval [1.5, 2] and 0.002 in the zoomed image, where 5 × 10 4 particles and 10000 averages at the steady state are taken. The critical coupling strength K c is given by K c 1.72584. As mentioned before, in this bimodal case, the forward ISP method is capable to describe very well the discontinuity in the phase transition. Note that, the backward ISP and the PMC methods fail to capture the correct jump location in the discontinuous phase transition. This might be due to the propagation of numerical errors leading to a jump outside the stability basin of the homogeneous state that is getting smaller and smaller as one approaches from the left the critical value of K . This is a strong indication of hysteresis phenomena usually characteristic of discontinuous phase transitions since both the homogeneous state and the second bifurcating branch have small basins of attraction coexisting in a range of the order parameter K .
Similarly as before, we also show the time evolution of ρ(t), the evolution ρ ∞ (D) with different values of D in Fig. 13, and the steady states ρ ∞ (ϑ, ω) and f ∞ (ϑ, ω) for K = 2, D = 0.5, N = 100, and M = 30 in Fig. 14. Note that it follows from

(4.3)
It is known that the above system (4.3) with h ∈ (0, 1) exhibits a hysteresis phenomenon. We refer to [23,24] for a description of the synchronization transition as a bifurcation problem. In order to analyze the phase transition as a function of the order parameter, we set the same initial data with (4.1) and take the Gaussian distribution g(ω) defined in Section 4.2.2.
In Fig. 15, we take D = 0.1, N = 200, M = 10, h = 0.5 to investigate the phase transition of the order parameter r ∞ (K ) numerically. We also take a mesh spacing in K of 0.01. As observed in [24], we find a different type of bifurcation with respect to Fig. 9. Fig. 15 shows an explosive jump from the incoherent state to the coherent one when the coupling strength is increased continuously. On the other hand, it also shows that a drop from the coherent state to the incoherent one when the coupling strength is decreased progressively. The critical coupling strengths are different from each other; a larger coupling strength is required to reach the coherent state from the incoherent one. This is a clear indication of the hysteresis phenomena as proven in [24, Figure 1(b)].

Conclusion
In this manuscript we focused our attention to the construction of effective numerical schemes for the solution of the continuum Kuramoto model with diffusion. For such system we introduce a discretization of the phase variable such that nonnegativity of the solution, physical conservations, asymptotic behavior and free energy dissipation are preserved at a discrete level. This approach is then coupled with a suitable collocation method for the frequency variable based on orthogonal polynomials with respect to the given frequency distribution. The method has been tested against the study of phase transitions in terms of the order parameter r ∞ as a function of K for D > 0. The numerical results agree very well with the theoretical basis, showing that the phase transition is continuous in the identical case and for Gaussian and uniformly coupled noisy oscillators [2], and discontinuous for bimodal distributions [16,25]. Hysteresis phenomena are also well-captured in the Kuramoto-Daido type model [23,24]. Future research directions and extensions of the present approach will consider the case of the Kuramoto model with inertia [2,10].