Finite-size scaling in globally coupled phase oscillators with a general coupling scheme

We investigate a critical exponent related to synchronization transition in globally coupled nonidentical phase oscillators. The critical exponents of susceptibility, correlation time, and correlation size are significant quantities to characterize fluctuations in coupled oscillator systems of large but finite size and understand a universal property of synchronization. These exponents have been identified for the sinusoidal coupling but not fully studied for other coupling schemes. Herein, for a general coupling function including a negative second harmonic term in addition to the sinusoidal term, we numerically estimate the critical exponent of the correlation size, denoted by $\nu_+$, in a synchronized regime of the system by employing a non-conventional statistical quantity. First, we confirm that the estimated value of $\nu_+$ is approximately 5/2 for the sinusoidal coupling case, which is consistent with the well-known theoretical result. Second, we show that the value of $\nu_+$ increases with an increase in the strength of the second harmonic term. Our result implies that the critical exponent characterizing synchronization transition largely depends on the coupling function.


I. INTRODUCTION
Populations of coupled rhythmic elements can exhibit synchronization and collective behavior via mutual interactions [1]. Such phenomena are observed in a variety of systems such as chemical reactions, engineering circuits, and biological populations [1]. To elucidate the general properties of such phenomena, the phase description of systems has been widely used [1,2]. In particular, there have been many studies on globally coupled phase oscillators [2,3], defined as follows:θ where θ j is the phase of the jth oscillator, ω j is the natural frequency of the jth oscillator, K > 0 is the coupling strength, h is the coupling function, and N is the number of oscillators. When h(x) = sin(x), this model is referred to as the Kuramoto model [2]. One of the main issues in this model has been the scaling property of the order parameter defined as follows [2]: In the thermodynamic limit (N → ∞), the phase oscillator model (1) exhibits a synchronization transition when the coupling strength K surpasses a critical value K c . This transition can be characterized by a change of the order parameter from zero to a non-zero value. We assume that the stationary state (R(t) = 0) in the incoherent regime supercritically bifurcates at the critical coupling strength K = K c , above which the oscillators are synchronized. The behavior of the order parameter is exemplified for a finite-size system in Fig. 1 [1 -6]. The scaling law of the order parameter with respect to a change in the coupling strength K around the synchronization transition point has been well studied in relation to the second order phase transition [2,4,5,[7][8][9]. The scaling property has been fully understood in the case where N is infinite, not only for the sinusoidal coupling function but also for general coupling functions [5,8,9]. However, it is less clear in finite-size systems because the property of fluctuations depends on the system size. The critical exponents of the order parameter, correlation time, susceptibility (corresponding to the product of the variance of the order parameter and N ), and correlation size (representing the number of oscillators which are almost synchronized but not completely) are significant quantities that can characterize the behavior near phase transitions in physical systems. The values of the critical exponents of these statistical quantities can be used to categorize physical systems into minor classes, because the values are independent of the details of the system [10]. Further, in equilibrium systems, the critical exponents of these statistical quantities completely determine those of all the other statistical quantities [10]. Therefore, the evaluation of the critical exponents in the phase oscillator model (1) enables the clarification of differences between equilibrium and non-equilibrium systems. In particular, fluctuations in the systems of large but finite size can be characterized by the exponents of susceptibility, correlation time, and correlation size. Although the critical exponents in the phase oscillator model (1) with finite large N have been obtained for the sinusoidal coupling function [11][12][13][14][15][16], it is known that the critical exponent of the order parameter depends on the coupling scheme [8,9]. This fact motivated us to examine if the critical exponents of other statistical quantities also depend on the coupling function or not.
In the present paper, we employ a non-conventional statistical quantity to evaluate the critical exponent of correlation size, ν + , in the synchronized regime of the phase oscillator model (1) with finite large N . This is because it is difficult to compute the value of ν + using the critical exponent of the order parameter [16]. The statistical quantity that we use is denoted by D, which is given by the diffusion coefficient of the temporal integration of the order parameter, multiplied by system size N [17]. Using the statistical quantity D, we perform the finite-size scaling analysis [10]. First, we confirm that the estimated value of ν + is approximately 5/2 for the sinusoidal coupling h(x) = sin(x), which is consistent with the well-known theoretical result [13,18]. Second, we consider a general coupling function including a negative second harmonic term in addition to the sinusoidal term, i.e. h(x) = sin(x)− q sin(2x) with q > 0. We show that the value of ν + increases with an increase in the strength q of the second harmonic term. Our result means that the critical exponent characterizing synchronization transition largely depends on the coupling function.
Although fluctuations of the order parameter in the phase oscillator model (1) have been studied for the past two decades [11][12][13][14][15][16][17], those for a general coupling function have not been fully understood. In particular, for any general coupling scheme other than the sinusoidal one, the critical exponents of statistical quantities have not been reported except for that of the non-conventional statistical quantity D [17]. Note that the values of critical exponents are independent of the details of systems and they characterize universal structures of the systems [10]. In fact, in the thermodynamic limit N → ∞, it was previously shown that as far as the coupling function includes a second harmonic term, the critical exponent of the order parameter takes the same value [8,9]. In contrast, our finding means that the value of ν + crucially depends on the strength of the second harmonic term. Our result highlights a new relation between the couping schemes and universal structures in the globally coupled phase oscillators (1). We consider the diffusion coefficient of the time integral of R(t), which characterizes long-term fluctuations in R(t) [17]. The variance of t 0 R(s)ds is given by:  (1), where h(x) = sin(x) and K = 1.635 > Kc = 1.59 · · · . The horizontal axis is normalized by D e, where · e is the ensemble average over 300 different initial conditions. We numerically obtained 300 values of D for each N and estimated the probability density function of D by means of kernel density estimation. The standard deviation of D is sufficiently small for each N .

II. STATISTICAL QUANTITY FOR THE FINITE-SIZE SCALING ANALYSIS
where R t represents the time average, defined as follows: The following diffusion law then holds: The value of D characterizes how t 0 R(s)ds differs from its mean value R t t on the ensemble average. The value of D is estimated by fitting σ 2 (t) with a line of slope 2Dt except for a transient period as shown in Fig. 2. Different initial conditions may yield different values of D because of multistability [19], as seen in Fig. 2. However, the dependence of D on initial conditions is sufficiently small as shown in Fig. 3. This result implies that the property on fluctuations is almost similar for the multiple coexisting solutions. All of them simultaneously transit from a non-synchronous state to a synchronous one as the coupling strength K exceeds the critical value K = K c for a large system size N [1][2][3][4][5][6]20].

III. THE FINITE-SIZE SCALING ANALYSIS
This study addresses the coupling function in the form h(x) = sin(x) − q sin(2x) with 0 ≤ q ≤ 0.1. We examine this coupling function because it is considered to be a "generic" coupling function [8,9]; provided that the coupling function includes the second harmonic term in addition to the sinusoidal one, the critical exponent of the order parameter does not depend on other higher-order terms [8,9].
The scaling hypothesis [10] states that any quantity A, which shows critical divergence at K = K c , is scaled for K > K c as follows: where ν + and r + represent the critical exponents of the correlation size and A in the synchronized regime, respectively. The function Ψ is approximated as Ψ(x) ∼ ux −r+ + x −s with s > r + and a constant u for large x [10] so that the orders with respect to system size N in both sides of Eq. (6) are consistent. It was shown that for a general coupling function h(x) in the phase oscillator model (1), the critical exponent of D in the synchronized regime is equal to 0 [17]. Therefore, we replace A by D and then substitute r + = 0 into Eq. (6). Furthermore, we assume u = 0 in the above approximation form of Ψ(x) because Ψ(x) → 0 as x → ∞. As a result, we obtain where g(ω) is the Gaussian distribution with mean zero and variance one. To avoid a situation in which the finite-size effect yields an erroneous value for ν + , we limit the range of K to (K c <)K − < K < K + . We select the lower value K − such that D decreases in a power law fashion with system size N for K > K − . Furthermore, we choose the upper value K + such that the power law decay of D is kept within the range K − < K < K + in order to avoid a large variation in the estimated values of ν + . In addition, we use a bootstrap method. First, we calculate the value of D by using 100 different initial conditions for each pair of (N, K). We randomly choose M overlapping samples of D values from the 100 simulation results and average them. Next, we estimate the value ofν   Figure 4(a) shows the value of D averaged over 100 overlapping samples for each pair of (N, K). The same data points shown in Fig. 4(b) are plotted against (K − K c )N 1/2.49 in a log-log plot. The data are fitted well by the straight line. The estimated exponentν (100) + ∼ 2.49 is consistent with the analytical results ν + = 5/2 [13,18]. Therefore, our method can estimate the exact value of ν + well.
Next, we investigate the model (1) with a more general coupling function with q = 0.1. Figure 5(a) shows the averaged values of D over 100 overlapping samples. The estimated mean value ofν (100) + is given asν (100) + ∼ 2.95, as seen in Fig. 5(b). This means that the value of ν + in this model is larger than that in the Kuramoto model.
Note that M = 100 samples are enough to estimate the accurate value of ν + because the standard deviation of the estimated values ofν term. The mean ofν (100) + is likely to increase with an increase in the value of q as shown in Fig. 7(a). In addition, we obtain the distributions ofν (1) + as shown in Fig. 7(b). The standard deviation ofν (1) + is not so large and the mean value ofν (1) + is also likely to increase with an increase in the value of q. These results imply that the critical exponent ν + depends on the second harmonic term of the coupling function.

IV. DISCUSSION
We explain why our method yields a good estimation for the critical exponent ν + . In the limit N → ∞, D → ∞ as K → K c − 0, whereas D = 0 for K > K c [17]. However, for a finite N , D takes its maximal value in the coherent regime as shown in Fig. 8. It means that D virtually reflects the feature of the incoherent solution of the infinite-size system even for K > K c if K ≃ K c and N is finite. Such a finite-size effect is well known in the literature [10]. As a result, to well estimate the value of ν + by the finite-size scaling analysis, we must avoid the coherent regime that is very close to the transition point K = K c . Our numerical simulations have selected the region in which the finite-size effect is weak enough so that D ∼ O(1/N a ). Note that it is difficult to determine whether the finite-size effect is sufficiently weak by using other statistical quantities.
We remark on another numerical study about the critical exponent of correlation size [16]. Let us denote the critical exponent of the order parameter R by β and that of correlation size in the incoherent regime by ν − . It is known that R ∼ N −β/ν for K = K c if ν + = ν − (=: ν). The value of ν can be computed as ν = 5/4 for the Kuramoto model with deterministically chosen natural frequencies [16]. However, there is little evidence for ν + = ν − . In fact, as we have already mentioned, the value of the critical exponent of D differs depending on whether the system behavior is coherent or incoherent [17].
Our numerical simulations have shown that when the coupling function possesses a negative second harmonic term with strength q in addition to the sinusoidal one, ν + > 5/2 and ν + increases with q. This result is consistent with the following property. Near the synchronization transition point K = K c , R is scaled as R ∼ (1 + 1/q)(K − K c ) + O((K − K c ) 2 ) [5]. This implies that the more the value of q increases, the more slowly the fluctuations of this system decay with the coupling strength K.
Our method has the following potential application. The method in this paper enables us to obtain an approximate value of ν + for a general coupling function. The value of ν + is the same for finite-size scaling analysis of other statistical quantities [10]. Suppose that we analytically obtain the value of ν + or the critical exponents of other statistical quantities such as susceptibility and correlation time for a general coupling function; our methods make it possible to confirm the validity of the analytical results numerically by the finite-size scaling analysis [10] because we can compute an approximate value of ν + .