Critical exponents in coupled phase-oscillator models on small-world networks
aa r X i v : . [ n li n . AO ] J u l Critical exponents in coupled phase-oscillator models on small-world networks
Ryosuke Yoneda, Kenji Harada, and Yoshiyuki Y. Yamaguchi
Graduate School of Informatics, Kyoto University, 606-8501 Kyoto, Japan
A coupled phase-oscillator model consists of phase-oscillators, each of which has the natural fre-quency obeying a probability distribution and couples with other oscillators through a given periodiccoupling function. This type of models is widely studied since it describes the synchronization tran-sition, which emerges between the non-synchronized state and partially synchronized states, andwhich is characterized by the critical exponents. Among them, we focus on the critical exponentdefined by coupling strength dependence of the order parameter. The synchronization transition isnot limited in the all-to-all interaction, whose number of links is of O ( N ) with N oscillators, andoccurs in small-world networks whose links are of O ( N ). In the all-to-all interaction, values of thecritical exponent depend on the natural frequency distribution and the coupling function, classi-fied into an infinite number of universality classes. A natural question is in small-world networks,whether the dependency remains irrespective of the order of links. To answer this question wenumerically compute the critical exponent on small-world networks by using the finite-size scalingmethod with coupling functions up to the second harmonics and with unimodal and symmetric nat-ural frequency distributions. Our numerical results suggest that, for the continuous transition, theconsidered models share the critical exponent 1 /
2, and that they are collapsed into one universalityclass.
I. INTRODUCTION
Ever since Huygens found that two pendulum clockshanging on a wall swung in the opposite direction fromeach other, many illustrations of synchronization havebeen established in various fields of nature, such as frogchoruses [1], flashing of fireflies [2, 3], metronomes [4],and circadian rhythms [5]. It is natural to try to under-stand synchronization theoretically, and a coupled phase-oscillator model is one of successful models to describesynchronization [6]. This model consists of many coupledoscillators, and the coupling is expressed by a periodiccoupling function, Each oscillator has the so-called natu-ral frequency, randomly drawn from a natural frequencydistribution. When the coupling strength K increases,the oscillators exhibit the synchronization transition fromthe non-synchronized state to (partially) synchronizedstates. The synchronization transition is continuous ordiscontinuous, depending on the natural frequency dis-tribution and the coupling function [7–17].The critical phenomena have been extensively studiedin statistical mechanics. One of their remarkable fea-tures is the existence of universality classes; the systemsin a universality class share the critical exponents de-fined around the critical point K c of a continuous tran-sition. One of the critical exponents is β , defined by r ∼ ( K − K c ) β , where r is the order parameter. Thus, itis natural to ask the universality classes in the coupledphase-oscillator models through values of the critical ex-ponent β .For the all-to-all and uniform coupling, extended re-searches have computed the value of β , and it dependson the coupling function and the natural frequency dis-tribution [7–17]. For simplicity, we focus on couplingfunctions which have two harmonics at most, and reviewvalues of the critical exponent β for the following threecases: (i) the second harmonics is absent, (ii) the second harmonics has the opposite sign with the leading har-monics, and (iii) the second harmonics has the same signwith the leading harmonics. We assume that the naturalfrequency distribution is unimodal and symmetric, andthat the second-leading term of its Maclaurin expansionis of the order 2 n .In the case (i) and (ii), the model shows a continu-ous transition, whereas in the case (iii), a discontinu-ous transition occurs, hence we cannot define the criti-cal exponent β [12]. In the case (i), the model becomesthe Kuramoto model [7], a paradigmatic coupled phase-oscillator model. Several researches have pointed outthat the critical exponent β = 1 / (2 n ) [7–10]. This n dependence is a strong feature of the Kuramoto modeland gives a sharp contrast with the case (ii). In the case(ii), the critical exponent β becomes 1 for n = 1 [11–14],and this value is suggested to be universal irrespective of n ∈ N [12].The universality of the coupled phase-oscillator mod-els depends also on the manner of coupling, that is, thenetwork of the coupling. For example, oscillator mod-els on complex networks are of interest to many naturalphenomena [18]. In particular, a previous research [19]has studied the continuous synchronous transition of thecoupled phase-oscillator model on a small-world network,one of complex networks, and claims that β = 1 / n = 1 in the case (i). However, the research lacks to con-sider other cases. In this paper, we numerically study thesynchronous transitions on small-world networks in thecases (i)-(iii) above. Our results suggest that in the case(i) and (ii), a continuous transition occurs, and β = 1 / n ∈ N .This paper is organized as follows. In Sec. II, we showthe algorithm of constructing the small-world network,then introduce coupled phase-oscillator models on small-world networks. We also introduce a family of the naturalfrequency distributions, whose second-leading term is of FIG. 1. Comparison between the all-to-all network (left) anda small-world network (right) with 20 nodes. In the small-world network, we set k = 3 and the rewiring probability p = 0 . the order 2 n . In Sec. III, we show a way to calculatethe critical exponent β using finite-size scaling, then wenumerically obtain β . Finally, in Sec. IV, we summarizethis paper and note some future works. II. COUPLED PHASE-OSCILLATOR MODELSON SMALL-WORLD NETWORKS
In this section, we define the coupled phase-oscillatormodel on a small-world network, and introduce the or-der parameter to observe synchronization. A small-worldnetwork possesses the property of a small diameter and alarge clustering coefficient despite its sparsity. This net-work can be seen in various fields of the real world, suchas human relationships, World Wide Web, citations ofscientific papers, and so on. In 1998, Watts and Stro-gatz proposed a breakthrough network model showingthe property of a small-world network, which is createdin the following algorithm [20]. We first make a periodic k -nearest neighbor network with N nodes, which resultsin kN edges. Then we rewire each edge with probability p , keeping in mind that we do not allow self-loops or linkduplications. Also, we use only connected small-worldnetworks: if a generated network is disconnected, we dis-card it and generate another one until connected one iscreated. See Fig. 1 for a comparison between the all-to-all network and a small-world network. In this paper, weset k = 3 and p = 0 .
2, following the previous researchon the coupled phase-oscillator model on a small-worldnetwork [19].We consider the coupled phase-oscillator model with N oscillators on a small-world network, which evolves intime with the following equations,d θ i d t = ω i + K k X j ∈ Λ i f a ( θ j − θ i ) ,f a ( θ ) = sin θ + a sin 2 θ, (1)for i = 1 , · · · , N . θ i and ω i are the phase and the naturalfrequency of the i th oscillator respectively, and ω i is ran- − − ω . . . . . . g n ( ω ) n = 1 n = 2 n = 3 n = 10 n = ∞ FIG. 2. Graphs of g n ( ω ) with n = 1 , , , , and ∞ , wherewe set ∆ = 1. We see that the g n ( ω ) converges to g ∞ ( ω ) as n becomes larger. domly drawn from a natural frequency distribution g ( ω ). K > i is the index set,which contains the indexes of oscillators connecting tothe i th oscillator. For each oscillator, the average num-ber of edges is 2 k , and we normalize the coupling termby 2 k . The coupling function is f a ( θ ), and when a = 0,this model is identical to the one proposed in [19].As the natural frequency distribution g ( ω ), we intro-duce a family of distributions, parametrized by a naturalnumber n ∈ N , g n ( ω ) = n Γ(1 / (2 n ))∆ e − ( ω/ ∆) n , (2)where Γ( z ) = R ∞ t z − e − t d t is the Gamma function de-fined on ℜ ( z ) >
0. Here, ∆ > n = 1 givesthe Gaussian distribution. For every n , this distributionis unimodal and symmetric with respect to ω = 0, andits Maclaurin expansion has the following form, g n ( ω ) = g n (0) − C n ω n + · · · , (3)where C n = n/ (Γ(1 / (2 n ))∆ n +1 ) is positive. We remarkthat the generalized Lorentzian distribution introducedin [21] also has the same expansion form. In the limit n → ∞ , g n ( ω ) converges to g ∞ ( ω ) in the L -norm, g ∞ ( ω ) = (cid:26) / (2∆) , ω ∈ ( − ∆ , ∆) , , otherwise . (4)This distribution is a uniform distribution on a compactsupport. See Fig. 2 for the graph of the distribution g n ( ω )and its convergence to g ∞ ( ω ).To visualize the extent of synchronization of oscillators,we introduce the order parameter r N ( K ) defined for eachcoupling strength K , r N ( K ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N N X j =1 e iθ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (5)This order parameter represents the centroid of the oscil-lators moving on the complex unit circle S . When theoscillators are uniformly distributed on S , which corre-sponds to the non-synchronized state, r N ( K ) gets closeto 0. On the other hand, when the oscillators gather at apoint on S , which corresponds to the synchronized state, r N ( K ) equals 1. Calculating r N ( K ) is therefore use-ful for monitoring synchronization of the coupled phase-oscillator models. In the next section, we will look intothe relationship between the order parameter r N ( K ) andthe coupling strength K .We remark a difference between the model consideredin [22, 23] and ours. In the literature, the authors con-sider the coupled phase-oscillator models on the small-world network with the equivalent equations as (1), butthey construct the small-world network with each nodeconnected to k = ⌊ sN ⌋ neighbors for s ∈ (0 , . k is constant regardless of N in our model. Here, ⌊ x ⌋ is the floor function that takes the greatest integer lessthan or equal to x . As a result, their network has O ( N )edges, and is much more dense than our network, whichhas O ( N ) edges. An advantage of their model is that itcan be analyzed through the equation of continuity us-ing the graphon [24], the continuous limit of the matrixrepresenting the couplings, of the small-world network.We note the dependency on p for a synchronizationtransition. For p = 0, the small-world network be-comes the periodic one-dimensional k -nearest-neighbornetwork. In [19], the authors numerically showed thata coupled phase-oscillator model on a k -nearest-neighbornetwork does not show a synchronization transition; r ≈ K . However, if p is greater than zero, a syn-chronization transition occurs because of emergence ofshortcuts. III. NUMERICAL SIMULATIONS
In the large population limit N → ∞ , the coupledphase-oscillator model (1) is expected to show a syn-chronization transition around a critical point K c . For K < K c , the order parameter r ( K ) := lim N →∞ r N ( K )is zero, which corresponds to the non-synchronized state.On the other hand, for K > K c , the model shows a par-tial synchronized state, in which r ( K ) exhibits power lawbehavior close to the critical point in the form of r ( K ) ∼ ( K − K c ) β , (6)where β is one of the critical exponents. The criticalexponents are crucial to describe critical phenomena,and models are classified into universality classes, eachof which shares the same critical exponents. Therefore calculating the critical exponents, including β , is an im-portant problem from theoretical and numerical perspec-tives. In this section, we numerically calculate the criti-cal exponent β of the coupled phase-oscillator models onsmall-world networks.The critical exponent β is defined in the large popula-tion limit N → ∞ , and this cannot be achieved throughthe numerical simulations. To overcome this difficulty,we use the finite-size scaling, that is, we assume thatthere exists a function F , so called the scaling function,and that the order parameter reads r N ( K ) = N − β/ ¯ ν F (( K − K c ) N / ¯ ν ) . (7)Here, the critical exponent ¯ ν is defined through the cor-relation length ξ , which diverges at K = K c with ξ ∝ ( K − K c ) − ¯ ν . (8)The finite-size scaling is widely used for numerical studiesof critical phenomena in second-order phase transitions,including coupled phase-oscillator models [19, 25–30].In this paper, we use a new statistical method [31, 32] toestimate the values of K c , β and ¯ ν in the finite-size scalingrelation. Using a new technique in the field of machinelearning, we can automatically find the best parameterset with a good collapse of data points onto a scalingfunction F .Numerical simulations of the coupled phase-oscillatormodel on the small-world networks are performed by us-ing the fourth-order Runge–Kutta algorithm with thetime step δt = 0 .
1. We first make a Watts–Strogatz’ssmall-world network, and for each K , we numericallyintegrate (1) by setting random initial phases and ran-dom natural frequencies obeying g n ( ω ), and we take thetime average of the order parameter in the time interval t ∈ [300 , N = 1600 to N = 25600. The order parameter r N ( K ) and its errorbar are evaluated by the resamplingtechnique: We randomly choose 200 samples out of 400realizations, and we calculate the mean of the order pa-rameter, which we refer to { r ( i ) N ( K ) } Si =1 . We take thisstep for S = 1000 times, and r N ( K ) is given by the meanof { r ( i ) N ( K ) } Si =1 , and the confidence interval of r N ( K ) isgiven by the statistical deviation of { r ( i ) N ( K ) } Si =1 . Wewill show results of a = 0 , − .
2, and 0 .
5, chosen fromthe neighborhood of a = 0, as the case (i), (ii), and (iii),respectively. The order of the second-leading term in thenatural frequency distribution are taken as n = 1 , , ∞ . See Fig. 3 for the numerical results with n = 1and a = 0 , − . r N ( K ) for several system size N , weevaluate the critical exponents β and ¯ ν and the criticalpoint K c from the scaling relation (7). We determinethese values by finding the best values so that points of r N ( K ) N β/ ¯ ν as a function of ( K − K c ) N / ¯ ν are tightlycollapsed to a scaling function F . Here we again usethe resampling method to obtain the critical exponents . . . . K . . . . . . r N ( K ) (a) N = 1600 N = 3200 N = 6400 N = 12800 N = 25600 1 . . . K (b) N = 1600 N = 3200 N = 6400 N = 12800 N = 25600 FIG. 3. Graphs of order parameter r N ( K ) with its confi-dence interval for the model (1), where we take the couplingfunction f a ( θ ) with (a) a = 0 and (b) a = − .
2. As a nat-ural frequency distribution, we use g ( ω ) with ∆ = 1, and N = 1600 , , , , and 25600 from top to bottom. r N ( K ) and its confidence interval are evaluated by the resam-pling technique. Errorbars are so small that they may not bevisible. . . . β . . . . . . . K c (a) N min = 1600 N min = 3200 N min = 6400 2 . . . ν (b) N min = 1600 N min = 3200 N min = 6400 FIG. 4. Scattering plots of computed parameters K c , β and¯ ν , evaluated by the Bayesian scaling analysis. Here, we use( a, n ) = (0 , N min to 1600 , and their confidence intervals: By introducing the symbol N min to indicate the smallest system size for the finite-size scaling analysis, we use r ( i ) N min ( K ) , r ( i )2 · N min ( K ), and r ( i )4 · N min ( K ) to evaluate β, ¯ ν , and K c by the Bayesian scal-ing analysis [31, 32]. Figure 4 shows β, ν , and K c esti-mated from each resampling data set. Mean values andstandard deviations of β, ¯ ν , and K c give us the best-fitparameters and confidence intervals. We carried out thisprocedure for N min = 1600 , r N ( K ) N β/ ¯ ν versus ( K − K c ) N / ¯ ν are shown in Fig. 5 forchecking the validity of the estimated parameters β, ¯ ν ,and K c .We check the hysteresis of (1) with a = 0 . { θ i } Ni =1 : (i) We startwith the random initial phases { θ i } Ni =1 at K = K start ,and the final phases at t = 500 is used as the initial − −
10 0 10( K − K c ) N / ¯ ν r N ( K ) N β / ¯ ν N = 3200 N = 6400 N = 12800 N = 25600 FIG. 5. Graph of scaled order parameter r N ( K ) N β/ ¯ ν versusscaled coupling constant ( K − K c ) N / ¯ ν for ( a, n ) = (0 , β, ¯ ν and K c , obtained by the Bayesian scal-ing analysis for N min = 6400. The values of β, ¯ ν , and K c are shown in Table I. We see that the scaled data are wellcollapsed to a scaling function F . phases at the successive value of K in the increasing di-rection. We call this process the “forward” process, and r (forward) N ( K ) denotes its order parameter. (ii) Contraryto the “forward” process, we start with the random ini-tial phases { θ i } Ni =1 at K = K end , and the final phases at t = 500 is used as the initial phases at the successive valueof K in the decreasing direction. We call this process the“backward” process, and r (backward) N ( K ) denotes its orderparameter. We have executed the numerical simulationsof (1) for a = 0 , − .
2, and 0 .
5, and n = 1 , ,
3, and ∞ ,and N = 25600 with two different initial phases as shownabove, and confirmed that there exists a hysteresis onlyfor a = 0 . n . See Fig. 6 for an exampleof a hysteresis with ( a, n ) = (0 . , a = 0 . N min = ∞ by plottingthe values of K c , β and ¯ ν as a function of 1 /N min ,and use the least square method to find the values atlim N min →∞ /N min = 0. The extrapolated values arelisted in the line of N min = ∞ in Table I. See Fig. 7for the extrapolation of β with a = 0 and − . β converge to β ≃ , (9)which is independent of a and n , unlike the coupledphase-oscillator on the all-to-all network. We noticethat the critical exponent ¯ ν obtained in [19], claimingthat ¯ ν ≃ a, n ) = (0 , TABLE I. Critical exponents β, ¯ ν and the critical point K c of (1) depending on the coupling function f a ( θ ) = sin θ + a sin 2 θ and the natural frequency distribution g n ( ω ) in (2), for a = 0 and − . n = 1 , ,
3, and ∞ . For each pair of ( a, n ), we use N min = 1600 , , N min = ∞ by using the least square method, and they are listed in the line of N min = ∞ . f a ( θ ) g n ( ω ) N min K c β ¯ νa = 0 n = 1 1600 2 . . . . . . . . . ∞ . ( ) . ( ) . ( ) n = 2 1600 1 . . . . . . . . . ∞ . ( ) . ( ) . ( ) n = 3 1600 1 . . . . . . . . . ∞ . ( ) . ( ) . ( ) n = ∞ . . . . . . . . . ∞ . ( ) . ( ) . ( ) a = − . n = 1 1600 2 . . . . . . . . . ∞ . ( ) . ( ) . ( ) n = 2 1600 2 . . . . . . . . . ∞ . ( ) . ( ) . ( ) n = 3 1600 2 . . . . . . . . . ∞ . ( ) . ( ) . ( ) n = ∞ . . . . . . . . . ∞ . ( ) . ( ) . ( ) ¯ ν ≃ /
2. They first find the best fit β/ ¯ ν and K c atwhich r N ( K ) N β/ ¯ ν crosses at K = K c by varying thesystem size N . Then they use the following formula,log (cid:20) d r N d K ( K c ) (cid:21) = 1 − β ¯ ν log N + const ., (10)which is calculated by taking the derivative of (7) withrespect to K , and they calculate (1 − β ) / ¯ ν by computingthe slope of log[ d r N d K ( K c )] with respect to log N . Howeverthis method has a disadvantage that the estimated crit-ical point K c is needed for calculating (1 − β ) / ¯ ν . Also,this method does not take into account the numericalresults r N ( K ) for K = K c . We believe that our estima-tion of the critical exponents is more reliable because ituses the Bayesian scaling analysis which overcomes thesedisadvantages. IV. CONCLUSION AND DISCUSSION
We calculated the critical exponents β and ¯ ν for thecoupled phase-oscillator models on the small-world net-work in (1) by using the finite-size scaling method. Weset the coupling function as f a ( θ ) = sin θ + a sin 2 θ , andthe natural frequency distribution as g n ( ω ) in (2), andwe studied the ( a, n )-dependency of the critical expo-nent β . Our results suggest that β is 1 / g n ( ω )and coupling function f a ( θ ) with a = 0 and − .
2, whichmeans that all such systems are in the same universal-ity class. This universality shows a sharp contrast withthe Kuramoto model, which has various values of β de-pending on the coupling function and the natural fre-quency distribution. Remarking that the coupled phase-oscillator models on a small-world network with O ( N )edges [22, 23] shares the same critical exponent as theKuramoto model, whose network also has O ( N ) edges,we believe that this contrast comes from the number ofedges of networks. We also find that the model (1) shows . . . . . K . . . . . . r N ( K ) r (forward) N ( K ) r (backward) N ( K ) . . K . FIG. 6. Graphs of r N ( K ) and its errorbar of (1) for ( a, n ) =(0 . ,
1) with two different types of initial phases, where weset the number of oscillators N = 25600. We see that r N ( K )takes a different value depending on the choice of the initialphases around K ∈ (1 . , . r (backward) N ( K ) − r (forward) N ( K ). . . . . . β (a) a = 0 n = 1 n = 2 n = 3 n = ∞ / ∞ / / / /N min . . . . . β (b) a = − . n = 1 n = 2 n = 3 n = ∞ FIG. 7. Graphs of β as a function of 1 /N min for (a) a = 0 and(b) a = − . /N min = 1 / ∞ , which are calculated by the leastsquare method. For each a , the resulting linear regressionlines are drawn with the solid line for n = 1, the dashed linefor n = 2, the dot-dashed line for n = 3, and the dotted linefor n = ∞ . a discontinuous transition for a = 0 .
5. This discontinuityis shared between the two types of networks: networkswith O ( N ) edges and O ( N ) edges.We end this paper commenting on five future works.Firstly, we mainly focused on the dependence of criticalexponents of the natural frequency distribution, and wepicked representative points of a from around a = 0 tocheck universality in the coupled phase-oscillator modelson a small-world network. Studying a global phase dia-gram on the ( K, a )-plane is a subject for future research.Secondly, we note a relationship with the noisy Ku-ramoto model. In the Kuramoto model, the value of β depends on the natural frequency distribution, but byadding some noises, β = 1 / /N extrapolation performed in Fig. 7. These coinci-dence thus suggests introducing the small-world networkplays a similar role with applying noises because of itsinhomogeneous couplings.Thirdly, in addition to the critical exponent β , thereare still other critical exponents to calculate, γ and δ con-cerning with response for instance. Once we obtain threecritical exponents, we can check if the Widom equality γ = β ( δ −
1) holds. For the Kuramoto model, the Widomequality is obtained by analyzing the self-consistent equa-tion with respect to r [10]. However, we cannot derive theself-consistent equation for the coupled phase-oscillatormodel on the small-world network, and it is not obviousif the Widom equality holds.Fourthly, it remains to clarify α -dependence of the crit-ical exponent β in coupled phase-oscillator models with O ( N α )-edge networks (1 < α < α -dependencemay make a bridge to the gap between α = 2 and 1,which has been revealed in this paper.Finally, since we numerically obtained β , it is naturalto check the value theoretically. To investigate criticalphenomena in the model (1), one has to consider thelarge population limit N → ∞ , but to the best of ourknowledge, we do not have a tool for taking the largepopulation limit with network of O ( N ) edges. It is chal-lenging to attack this problem. ACKNOWLEDGMENTS
In this research work we used the supercomputer ofACCMS, Kyoto University. Y.Y.Y. acknowledges thesupport of JSPS KAKENHI Grant No. 16K05472. [1] I. Aihara, T. Mizumoto, T. Otsuka, H. Awano, K. Nagira,H. G. Okuno, and K. Aihara, Spatio-Temporal Dynam-ics in Collective Frog Choruses Examined by Mathemat-ical Modeling and Field Observations,
Sci. Rep. , 3891(2014).[2] H. M. Smith, Synchronous flashing of fireflies, Science , 151 (1935).[3] J. Buck and E. Buck, Mechanism of rhythmic syn-chronous flashing of fireflies, Science , 1319 (1968).[4] J. Pantaleone, Synchronization of metronomes,
Am. J.Phys. , 992 (2002).[5] A. T. Winfree, Biological Rhythms and the Behavior ofPopulations of Coupled Oscillators, J. Theoret. Biol. ,15 (1967).[6] Y. Kuramoto and H. Nakao, On the concept of dynamicalreduction: the case of coupled oscillators, Phil. Trans. R.Soc. A , 20190041 (2019).[7] Y. Kuramoto, Self-entertainment of a population of cou-pled non-linear oscillators,
International Symposium onMathematical Problems in Theoretical Physics, LectureNotes in Physics , 420 (1975).[8] S. H. Strogatz, From Kuramoto to Crawford: exploringthe onset of synchronization in populations of coupledoscillators, Physica D , 1 (2000).[9] H. Chiba, A proof of the Kuramoto conjecture for a bi-furcation structure of the infinite-dimensional Kuramotomodel,
Ergod. Th. & Dynam. Sys. , 762 (2015).[10] H. Daido, Susceptibility of large populations of coupledoscillators, Phys. Rev. E , 012925 (2015).[11] J. D. Crawford, Scaling and Singularities in the Entrain-ment of Globally Coupled Oscillators, Phys. Rev. Lett. , 4341 (1995).[12] H. Chiba and I. Nishikawa, Center manifold reduction forlarge populations of globally coupled phase oscillators, Chaos , 043103 (2011).[13] M. Komarov and A. Pikovsky, Multiplicity of SingularSynchronous States in the Kuramoto Model of CoupledOscillators, Phys. Rev. Lett. , 204101 (2013).[14] M. Komarov and A. Pikovsky, The Kuramoto model ofcoupled oscillators with a bi-harmonic coupling function,
Physica D , 18 (2014).[15] L. Basnarkov and V. Urumov, Phase transitions in theKuramoto model,
Phys. Rev. E , 057201 (2007).[16] D. Paz´o, Thermodynamic limit of the first-order phasetransition in the Kuramoto model, Phys. Rev. E ,046211 (2005).[17] H. Daido, Intrinsic Fluctuations and a Phase Transitionin a Class of Large Populations of Interacting Oscillators, J. Stat. Phys. , 753 (1990).[18] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes,Critical phenomena in complex networks, Rev. Mod.Phys. , 1275 (2008).[19] H. Hong, M. Y. Choi, and B. J. Kim, Synchronization onsmall-world networks, Phys. Rev. E , 026139 (2002).[20] D. J. Watts and S. H. Strogatz, Collective dynamics of‘small-world’ networks, Nature , 440 (1998).[21] B. Pietras, N. Deschle, and A. Daffertshofer, First-orderphase transitions in the Kuramoto model with com-pact bimodal frequency distributions,
Phys. Rev. E ,062219 (2018).[22] H. Chiba, G. S. Medvedev, and M. S. Mizuhara, Bifur-cations in the Kuramoto model on graphs, Chaos ,073109 (2018).[23] G. S. Medvedev, Small-world networks of Kuramoto os-cillators, Physica D , 13 (2014).[24] L. Lovsz, Large networks and graph limits, (AmericanMathematical Society, Providence, Rhode Island, 2012).[25] A. Pelisseto and E. Vicari, Critical phenomena andrenormalization-group theory,
Phys. Rep. , 549(2002).[26] M. Hasenbusch, Finite size scaling study of lattice modelsin the three-dimensional Ising universality class,
Phys.Rev. B , 174433 (2010).[27] A. Pikovsky and M. Rosenblum, Dynamics of globallycoupled oscillators: Progress and perspectives, Chaos ,097616 (2015).[28] H. Hong, H. Chat´e, L-H. Tang, and H. Park, Finite-sizescaling, dynamic fluctuations, and hyperscaling relationin the Kuramoto model, Phys. Rev. E , 022122 (2015).[29] T. Coletta, R. Delabays, and P. Jacquod, Finite-size scal-ing in the Kuramoto model, Phys. Rev. E , 042207(2017).[30] R. Juh´asz, J. Kelling, and G. ´Odor, Critical dynamicsof the Kuramoto model on sparse random networks, J.Stat. Mech. , 053403 (2019).[31] K. Harada, Bayesian inference in the scaling analysis ofcritical phenomena,
Phys. Rev. E , 056704 (2011).[32] K. Harada, Kernel method for corrections to scaling, Phys. Rev. E , 012106 (2015).[33] H. Sakaguchi, Cooperative Phenomena in Coupled Oscil-lator Systems under External Fields, Prog. Theo. Phys.79