Synchronization is a very common phenomenon in different fields, such as physics, chemistry, biology and so on. It is also the basic mechanism to reveal the collective phenomenon in nature and man-made systems. The traditional discussion of these phenomena is mainly based on the coupled limit-cycle system. Theoretically, synchronization has been successfully studied in models of coupled oscillators  . Among these models, Kuramoto model is the simplest in which the dynamics of an oscillator are described by only a phase variable   . The Kuramoto model not only introduces statistical methods, but also links the synchronous behavior in the field of dynamics with the non-equilibrium phase transition in statistical physics.
In recent years, the influence of the frustration factors on the synchronization transition of coupled oscillators has attracted extensive attention  - . The frustration factor is one of the most important parameters in phase transition models such as Heisenberg XY model, Frenkel-Kontorova model and so on. For example, in frustrated large square Josephson junction arrays, the frustration factor plays a crucial role. In the field of engineering technology, the frustration factor can be easily realized by adding a magnetic field.
Inspired by the law of nature, we introduce the frequency mismatch rules into the Sakaguchi-Kuramoto model, in which only the two oscillators with large frequency difference can be connected directly. This idea is based on the following facts. In daily life, two people with different personalities complement each other and make up for the missing part of their personalities, so they often become good friends. It is a natural law that differences complement each other and opposites attract.
The rest of this paper is organized as follows. In Section II, some basic notions on Sakaguchi-Kuramoto model with frequency mismatch rules are introduced. Then in Section III, we give the main results of numerical simulations. The synchronization dynamics of the two-oscillator system are numerically simulated and theoretically analyzed in Section IV. Section V is devoted to the concluding remarks, in which we analyze the limitation of current work and put forward the next work plan.
2. Basic Notions on Sakaguchi-Kuramoto Model
In the current work, our aim is to study the synchronization dynamics of coupled oscillators on complex networks. Therefore, the first task is to build a substrate network, which follows the steps below  . 1) The natural frequency of the ith oscillator is set to . 2) Each edge added randomly must stand the test of inequality (1), otherwise it will be discarded. Adding an edge repeatedly is prohibited. 3) The process of adding edges does not stop until the given L edges are obtained.
where and are the natural frequencies of the nodes locating at both ends of the jth edge, is the threshold value of average frequency gap of the whole system.
After completing the substrate network, each node is embedded with a Sakaguchi-Kuramoto oscillator whose natural frequency has been specified. The dynamical behaviors of N coupled oscillators are governed by the following equations ,
where is the phase of the ith oscillator, is its natural frequency, is a positive coupling strength, N is the number of oscillators in the entire network. The oscillators are coupled with each other according to the elements of the network’s adjacency matrix , with if there is a connection between oscillator i and oscillator j, and otherwise. Here, is a phase frustration factor that can prevent the phases of two connected oscillators from getting too close.
We now fix and for all simulations in the present work, since the results described below do not change qualitatively for larger system size or greater frequency mismatch. We keep as a constant throughout this work. Therefore, the average degree of each oscillator is . Typically, the collective behavior in Equation (2) can be monitored by the phase order parameter  :
where is an index of phase coherence. When all the oscillators have the same phase, reaches one, and for the completely random state, it drops to zero. In the following numerical simulations, Equation (2) is integrated by the fourth-order Runge-Kutta method with time step 0.01. The initial phase of each oscillator is chosen from the interval at random.
3. Main Results
In this part, we mainly focus on the influence of frustration factors on the synchronization transition of coupled oscillators. In view of the periodicity of sine function and its symmetry with respect to , the value of frustration factor is tuned in the interval .
In Figure 1, both synchronization and desynchronization diagrams are plotted for four different values of . From to , the peak value of decreases to zero first and then rises slightly. The main reason is that the frustration factor forces the phase difference between the directly connected oscillators to decay, thus suppressing the emergence of phase synchronization.
Furthermore, it is noted that with the increase of , the area of hysteresis loop shows the same change rule. When is growing from 0 to , the mode of phase synchronization changes from explosive synchronization to continuous phase transition. When the value of reaches , no matter how strong the coupling strength is, it can not wake up the phase synchronization between the oscillators. When exceeds , the synchronous ability of network rises slightly under weak coupling strength, but with the increase of coupling strength, the phase synchronization gradually annihilates. This abnormal
Figure 1. The phase order parameter vs. the coupling strength for different values of . (a) , (b) , (c) , (d) . The forward (backward) continuations shown in the panels refer to the simulations of adiabatic increase (decrease) of coupling strength.
phenomenon naturally gives rise to the question: with the increase of coupling strength, how do the phases of oscillators evolve in the system?
The natural frequencies and instantaneous phases of all oscillators in the system are plotted on a series of unit circular surfaces, in which each black dot represents an oscillator, as shown in Figure 2. The polar radius and the polar angle of the black dot represent the natural frequency and the instantaneous phase of the oscillator, respectively.
We find that the phase distributions of the oscillators show diversity for different frustration factors at the same coupling strength. A typical partial phase-locked state is illustrated in Figure 2(a). The phases of the oscillators are completely random and evenly scattered on the unit circular surface, as shown in Figure 2(b). In Figure 2(d), the phases of oscillators are diametrically opposed, separated by an angle , which is commonly known as a state. Interestingly, all the oscillators spontaneously split into two subgroups bounded by the average frequency of the system. Figure 2(c) is a mixture of Figure 2(b) and Figure 2(d), in which a completely incoherent state and a state coexist. By comparing Figure 1 with Figure 2, one may notice that although the phase order parameters are the same on the macro level, the phase distributions of oscillators in the system are quite different.
Next, we reveal the frequency dynamics of coupled oscillators from a microscopic point of view. The effective frequency of each oscillator in steady state is defined as ,
Figure 2. Phase distribution of four hundred oscillators on the unit circular surfaces for different values of . (a) , (b) , (c) , (d) . A black dot represents an oscillator. The natural frequency and instantaneous phase of an oscillator are described by the polar radius and polar angle of the black dot, respectively. The remaining parameter is set to 0.10.
where T is an average time and is a discarded transient time.
Figure 3(a) shows that with the increase of coupling strength the average effective frequency of the oscillators becomes smaller and smaller and diverges in the negative direction when is less than . The frequency synchronization transition shown in Figure 3(c) is mirror-symmetric with that shown in Figure 3(a). Figure 3(b) shows that when is equal to , the divergence degree of the effective frequency of the oscillators in the system increases with the increase of coupling strength, rather than the expected decrease. What is most striking is that all the oscillators stick to their own natural frequencies until the emergence of an abrupt synchronization, as shown in Figure 3(d). There is no premonition for this kind of synchronization transition.
This phenomenon of frequency synchronization cannot be detected by the phase order parameter. Therefore, it is necessary to define a new index to characterize the degree of frequency synchronization of coupled oscillators. The frequency order parameter is given by,
Figure 3. The panels show the evolution of the effective frequencies of four hundred oscillators with the increase of coupling strength for different values of . (a) , (b) , (c) , (d) . Each data point in the panels is the average of the 2000 time steps after the transient process.
where represents the standard deviation of instantaneous frequency of all oscillators in the system at time t. According to the definition, is between 0 and 1. The higher the value of , the higher the degree of frequency synchronization. Specifically, indicates a totally random frequency distribution, while represents that the instantaneous frequencies of all oscillators converge to the same value.
Figures 4(a)-(c) show a similar frequency evolution law, in which the degree of frequency synchronization decreases with the increase of coupling strength. However, the most surprising result is observed for the case of , in which an abrupt frequency synchronization transition appears, as shown in Figure 4(d). In addition, we find that there is a hysteresis loop between the forward and backward continuations via changing of the coupling strength . The region of hysteresis loop in Figure 4(d) is consistent with that in Figure 1(d). Because of the introduction of the frequency order parameter, the explosive synchronization of frequencies is easily detected.
Figure 5a and Figure 5(b) display the dependence of ( ) on and . By comparing Figure 5(a) with Figure 5(b), one may find that both panels are axisymmetric. Their axes of symmetry are both . The phase locking in Figure 5(a) inevitably induces the frequency synchronization in Figure 5(b).
Figure 4. The frequency order parameter vs. the coupling strength for different values of . (a) , (b) , (c) , (d) . The forward (backward) continuations shown in the panels refer to the simulations of adiabatic increase (decrease) of coupling strength.
Figure 5. (Color online). Color mark plots representing the phase synchronization degree (a) and the frequency synchronization degree (b) in the space of and , respectively. The greater the value of the color scale, the higher the degree of coherence.
On the contrary, the conclusion is not true. Figure 5(b) has one more synchronization area than Figure 5(a), in which the instantaneous phases of the oscillators are locked in reverse and the effective frequencies are the same.
4. Theoretical Analysis
In this section, we concentrate ourselves on the simplest case, i.e., . In this way, the complex multi-body system can be simplified into two-body problem, which is convenient for theoretical analysis. The dynamic equations of the two oscillators are written as
Once the phase difference is inserted in Equation (5) and Equation (6), one obtains
By subtracting Equation (9) from Equation (8), the evolution of phase difference with time is obtained,
In the phase-locked state, the left-hand side of Equation (10) is equal to 0. Thus, the sine of the phase difference is given by
Considering the boundedness of sine function, one can obtain
By adding Equation (1) to Equation (2), the effective frequency is obtained as follow
Equation (11) is taken into Equation (13), and the expression of effective frequency is rewritten as
where the minus sign “−’’ and plus sign “+’’ correspond to two different intervals and , respectively. This theoretical result successfully explains why the divergence directions of the effective frequencies in Figure 3(a) and Figure 3(c) are opposite. Except for the case of , when the coupling strength exceeds the critical point of the system, the curves obtained by theoretical analysis pass through each discrete point of numerical simulations, as shown in Figure 6. The theoretical analysis is in good agreement with the numerical simulation.
We have studied phase locking and frequency entrainment of coupled oscillators in the Sakaguchi-Kuramoto oscillator network with frequency mismatch rules.
Figure 6. The panels show the evolution of the effective frequencies of two oscillators with the increase of coupling strength for different values of . (a) , (b) , (c) , (d) . Each data point in the panels is the average of the 2000 time steps after the transient process.
The results show that the frustration factor plays an important role in synchronization transition of networks. The frustration factor can control the types of synchronization transition. Appropriate adjustment of frustration factors can eliminate the hidden dangers of harmful synchronization. We propose the frequency order parameter, which makes up for the defect that the phase order parameter can not detect the frequency synchronization. In order to verify the correctness of numerical simulations of the multi-oscillator system, we degenerate the multi-body system into a two-oscillator system. We obtain the analytical solution of the two-oscillator system. Theoretical analysis supports the simulation results. Unfortunately, we are incapable of obtaining analytical solutions of multi-oscillator system. This will be the key part of our next study. In the future research, we will use the stochastic noise to replace the fixed frustration factor, which is more suitable for the real world. We believe that the method developed in this paper also provides a useful tool to study the analogous problems such as time-delay effects. Suppose this model is applied to the study of time-delay effects, one only needs to rewrite the frustration factor as the ratio of the distance between oscillators to the speed of information transmission. The effect of variable frustration factors on the synchronization dynamics of complex networks will be an interesting topic.
This work is supported by the National Natural Science Foundation of China under Grant Nos. 61563054 and 11875031, the Natural Science Foundation of Guangxi under Grant No. 2019JJA110069, the open fund of Guangxi Colleges and Universities Key Laboratory of Complex System Optimization and Big Data Processing under Grant No. 2015CSOBDP0101, Initiation Fund of Doctoral Research of Yulin Normal University under Grant No. G20150003.
 Boccaletti, S., Almendral, J.A., Guan, S., et al. (2016) Explosive Transitions in Complex Networks’ Structure and Dynamics: Percolation and Synchronization. Physics Reports, 660, 1-94.
 Acebrón, J.A., Bonilla, L.L., Vicente, C.J.P., et al. (2005) The Kuramoto Model: A Simple Paradigm for Synchronization Phenomena. Reviews of Modern Physics, 77, 137.
 Gao, Y.-C., Fu, C.-J., Cai, S.-M., et al. (2019) Repulsive Synchronization in Complex Networks. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29, Article ID: 053130.
 Brede, M. and Kalloniatis, A.C. (2016) Frustration Tuning and Perfect Phase Synchronization in the Kuramoto-Sakaguchi Model. Physical Review E, 93, Article ID: 062315.
 Khanra, P., Kundu, P., Hens, C., et al. (2018) Explosive Synchronization in Phase-Frustrated Multiplex Networks. Physical Review E, 98, Article ID: 052315.
 Kundu, P. and Pal, P. (2019) Synchronization Transition in Sakaguchi-Kuramoto Model on Complex Networks with Partial Degree-Frequency Correlation. Chaos, 29, Article ID: 013123.
 Leyva, I., Navas, A., Sendiña-Nadal, I., et al. (2013) Explosive Transitions to Synchronization in Networks of Phase Oscillators. Scientific Reports, 3, Article No. 1281.
 Zhu, L.H., Tian, L. and Shi, D.N. (2013) Criterion for the Emergence of Explosive Synchronization Transitions in Networks of Phase Oscillators. Physical Review E, 88, Article ID: 042921.
 Sakaguchi, H. and Kuramoto, Y. (1986) A Soluble Active Rotater Model Showing Phase Transitions via Mutual Entertainment. Progress of Theoretical Physics, 76, 576-581.
 Gomez-Gardenes, J., Gomez, S., Arenas, A., et al. (2011) Explosive Synchronization Transitions in Scale-Free Networks. Physical Review Letters, 106, Article ID: 128701.