Synchronization in complex networks characterizes what happens when an ensemble of oscillators in a complex autonomous system become phase-locked. We study the Kuramoto model with a tunable phase-lag parameter $\alpha $ in the coupling term to determine how phase shifts influence the synchronization transition. The simulation results show that the phase frustration parameter leads to desynchronization. We find two global synchronization regions for $\alpha \u2208[0,2\pi )$ when the coupling is sufficiently large and detect a relatively rare network synchronization pattern in the frustration parameter near $\alpha =\pi $. We call this frequency-locking configuration as “repulsive synchronization,” because it is induced by repulsive coupling. Since the repulsive synchronization cannot be described by the usual order parameter $r$, the parameter frequency dispersion is introduced to detect synchronization.

Networks of coupled oscillatory units have received wide concern because it helps understand the synchronization in a variety of systems such as biology, finance, and transportation. Most of the studies focused on attractive couplings between the oscillators, in which all nodes can be synchronized to be either identical or just detuned within a small parameter range. In this paper, we investigate the competing effect of attractive and repulsive interactions in oscillator ensembles. We study the Kuramoto model with a tunable parameter $\alpha $, and an intrinsic frequency of each oscillator is fixed to be proportional to its degree. It is found that there are two types of synchronization, normal and repulsive ones. Since the latter is not described as the Kuramoto global order parameter that is based on the phase-locking, we, therefore, introduce a parameter that measures the frequency dispersion. When $\alpha =\pi $, even though phases are not locked, the frequencies are well entrained. This repulsive synchronization is known as the frequency-locking state, which was only found in regular lattices and grids that repulsive coupling leads to stable anti-ferromagnetic-type states. Here for the first time, the repulsive synchronous state is observed on disordered networks.

## I. INTRODUCTION

Complex dynamical systems usually exhibit macroscopic coherent states when their coupling strength exceeds a threshold. This behavior is ubiquitous in such far-from-equilibrium real and artificial systems as the synchronous rhythmic light pulse in a large population of flashing fireflies,^{1} the optical-forced oscillatory Belousov–Zhabotinsky reaction,^{2} and the synchronization of neuronal activities.^{3} Synchronization in complex networks is an extensively explored topic, in particular the examination of the interaction between network structure and synchronized dynamics, because studying the relationship between dynamic behavior and coupling structure is essential to understand, predict, and control real-world complex systems.^{4–15}

Although much has been learned about the synchronizations of large amount of oscillators with attractive coupling, the effect of repulsive coupling on the synchronization is not so clear. Actually, the systems of repulsive interactions commonly exist in ecological systems^{16} and human society. For instance, in the traffic systems, the vehicles are trying to maintain a safe gap between neighbors, and this speed adaption mechanism is found to form a synchronized flow state.^{17,18}

In order to investigate the effect of repulsive interactions between oscillators, we use the Sakaguchi–Kuramoto model,^{19} which introduced a phase shift parameter in the original Kuramoto model,^{20} and the coupling interactions between oscillators can change from attraction to repulsion when tuning the phase shift. Previous research has shown that phase shifts strongly impact the behavior of many systems—e.g., they can induce frustration in such disordered systems as the Josephson series array^{21} and the XY spin-glass model^{22}—because of the competition between the difference in a pair of oscillators and the phase shift that leads to either a positive or negative weak coupling. Futhermore, it has shown that with the repulsive interactions, the system will turn to partial synchronization from fully synchronized state, and also traveling wave is observed.^{23} Recently, some works have investigated about synchronizing phase frustrated Kuramoto oscillators^{24–26} and found a complex relationship between the phase shifts and the synchronization phase transition. In these works, the phase shift is focused on the range of $[0,\pi /2]$. In addition, recent research has focused on first-order phase transitions in complex systems;^{27–33} however, how frustration affects the first-order synchronization phase transition has received less attention.

In the current work, we focus on the change of the synchronous state when coupling of the oscillators turns from attraction to repulsion when the phase shift changes continuously from 0 to $2\pi $. We find that phase frustration can induce desynchronization, and most interestingly, repulsive coupling can lead to the frequency-locking synchronization with homogeneous phase distribution in the disordered systems, which is called repulsive synchronization. Also, the repulsive synchronization behavior in heterogeneous networks is investigated and found the first-order phase transition in the system with a heterogeneous topological structure. This helps to understand thoroughly the synchronization state in different coupling structures. The investigation of this repulsive synchronous state may give a hint in the regulation of the complex systems with the repulsive coupling effect.

## II. MODEL

In a complex network of $N$ coupled oscillators in which each has a phase $\theta i\u2208[0,2\pi )$ and a natural frequency $\omega i$, then the instantaneous phase velocity is $vi=\theta \u02d9i$. Thus, the dynamical equation of phase oscillators is

This model with a constant phase lag $\alpha $ is well known as Sakaguchi–Kuramoto model.^{19} The coupling strength of two arbitrary oscillators $i$ and $j$ is $\lambda $, when they are connected such that $Aij=1$; otherwise, they are not directly coupled when $Aij=0$. Here, $\alpha $ is a phase frustration parameter that can introduce a constant phase shift between connected oscillators.

Previous research^{34,35} defines the order parameter $r$ that quantifies the degree of synchronization

When the system is fully synchronized, $r\u21921$ and all oscillators are locked in a common average phase. In contrast, when $r\u21920$, the system is incoherent. Thus, the synchronous state can be indicated by the $r$ value, but as an order parameter $r$ cannot be applied to all coherent states, instead, is valid only for global phase-locked synchronization, and cannot give accurate partial synchronization results with certain symmetries. For instance, if there are a variety of synchronized clusters in a synchronous system and the oscillating phases vary between different clusters, then $r$ may be quite small. In this paper, we propose a new order parameter $\sigma \Omega $ to deal with this kind of coherent state that numerically describes the phase transition of synchronization including the cluster synchronization.

It was found that for scale-free networks, there is an explosive phase transition from $r=0$ to $r=1$ with increasing $\lambda $, if the natural frequencies are in positive correlation with node degrees $ki$.^{31} Here, we set $\omega i=ki$, in order to focus on the effect of frustration in the first-order synchronization transition.

## III. RESULTS

The Kuramoto synchronization process is simulated on a Barabási-Albert (BA)^{36} network of $N=1000$ oscillators with an average degree $\u27e8k\u27e9=6$. The value of phase shift $\alpha $ is tuned in the range of $[0,2\pi )$ due to the periodicity of trigonometric functions.

In Fig. 1, the synchronization diagrams of six different $\alpha $ are presented in the six panels, respectively. Here, we focus on $\alpha <\pi $ because the synchronization behavior of $[\pi ,2\pi )$ is similar to $[0,\pi )$, as we will see below. From $\alpha =0$ to $1.0$, $r$ keeps dropping until zero when $\alpha $ exceeds $1.0$, where no synchronous transition is observed from the synchronization diagram. This phenomenon is referred to as desynchronization, because fully synchronization state is destroyed by the phase frustration, and only a cluster of oscillators can come up to the phase-locking state.

In addition, we can also recognize that the critical coupling strengths $\lambda $ of both forward and backward continuations in the synchronization phase transition are raised with increasing $\alpha $. Also, it is noted that the width of hysteresis in the synchronization diagram (i.e., the difference between forward and backward critical couplings) increases when $\alpha $ is growing from 0 to 0.6, but declines after $\alpha >0.6$.

To examine the frequency dynamics as a function of $\alpha $, we calculate the effective frequency $\Omega i$ to be

In contrast with the natural frequency, $\Omega i$ is the observed frequency of an oscillator. With this parameter, we explore how the effective frequencies of the system change with phase frustration $\alpha $ at certain coupling strengths and form the synchronization tree shown in Fig. 2.

Figure 2(a) shows that, when the coupling strength is sufficiently weak, the effective frequency of all oscillators almost remains unchanged and is consistent with its node degree at every $\alpha $. Nevertheless, in Fig. 2(b), two synchronized regions emerge at the vicinity of $\alpha =0$ and $\alpha =\pi $, respectively, as the coupling value is increased to $\lambda >\lambda c$, where $\lambda c$ is the critical point of synchronized transition in the absence of phase shift. Beyond the two synchronized regions, however, when $\alpha $ deviates from $0$ or $\pi $, a number of the oscillators will be separated from the synchronous cluster, leading to first partial synchronization and finally to complete desynchronization. Figure 2(c) shows that when $\lambda $ is further increased, the position of the synchronization area changes slightly, but the width of the area remains approximately the same. Note that there is an inversion in the curve where the effective frequencies $\Omega i$ convert from positive to negative values. This is because some of the coupling terms in Eq. (1) are negative, causing a negative total coupling strength, and thus negative effective frequency. Moreover, effective frequency $\Omega $ of nodes with higher degrees are more negative than that with less degrees, which can be clearly observed from Fig. 2(c), when $\alpha \u2208(\pi ,2\pi )$, The top three lines with highest degrees (blue, red, and purple) goes to the bottom. The reason is that the sum of the negative coupling terms is more negative for higher degree nodes.

There are two synchronous regions in the synchronization tree described above, one near $\alpha =0$ and the other near $\alpha =\pi $. Figure 3 shows a plot of the order parameter $r$ in the steady state $rmax$ as a function of $\alpha $ that verifies this. The line indicates when $\alpha $ is sufficiently greater than 0, at which point the desynchronization is induced and the system partially synchronized, consistent with Fig. 2. Note, however, that although frequency entrainment is realized when $\alpha =\pi $ from the synchronization tree, no phase-locking is present because $rmax$ is approaching 0. Although for both $\alpha =0$ and $\alpha =\pi $ the system is clearly synchronized with entrained frequencies among all oscillators, $rmax$ is different. This indicates that the synchronization states of the two areas of entrained effective frequency have different patterns.

We find that the phase distributions of the oscillating system near $\alpha =0$ and $\alpha =\pi $ differ. Figure 4 reveals that the phase of the oscillators is approximately equal near $\alpha =0$. They form a phase-locked cluster that moves synchronously, and order parameter $r$ reaches its maximum. In contrast, when $\alpha $ is near $\pi $, the distribution of the oscillator phases is homogeneous in $[0,2\pi )$, and the value of order parameter $r$ is low. But the entire whole system is frequency locked, i.e., the frequency of all oscillators is the same. This anomalous phase synchronization, which has been rarely observed in network synchronization, differs from the traditional phase-locked synchronization and thus cannot be described using the order parameter $r$ because it is characterized by a phase distribution instead of a frequency distribution. We thus introduce a new indicator to describe synchronization by defining the dispersion of effective frequencies to be

Here, $\u27e8\Omega \u27e9$ is the average value of the effective frequencies. If the oscillators are phase-locked, the observed frequencies are locked in a common average frequency, the dispersion of effective frequency is thus considerably small.

Figure 5 compares the abilities of normalized $\sigma \Omega $ and $r$ to detect a synchronized phase transition. The normalized $\sigma \Omega $ is obtained by dividing each original $\sigma \Omega $ by the value of $\sigma \Omega $ at $\lambda =0$. Figure 5(a) shows that when $\alpha =0$, the value of $r$ increases from 0 to 1 when $\lambda $ exceeds the critical coupling. Meanwhile, normalized $\sigma \Omega $ decreases from 1 to 0 during the synchronization transition, in accordance with $r$. Figure 5(b) shows that when $\alpha =\pi $, the order parameter $r$ cannot detect a synchronization transition, but normalized $\sigma \Omega $ detects in the system a phase transition to the frequency-locking state that forms a hysteresis loop between the forward and backward continuations.

Thus, phase frustration can induce an anomalous type of global synchronization state which has not been observed in disordered systems. The difference between the two synchronization states is in their coupling patterns. Near $\alpha =0$, the coupling between adjacent oscillators is an attractive interaction that concentrates the phase and increases the $r$ value. In contrast, when $\alpha $ is near $\pi $, the weak coupling turns into a repulsive force between each pair of oscillators. Suppose there are only two oscillators, the phase difference is $\pi $, and as the population grows, they become uniformly distributed and the phase differences among the oscillators converge. When $\alpha \u2208(0,\pi )$, the system is either partially synchronized or desynchronized, as a result of the competition of the attractive and repulsive coupling. Figure 5(c) shows a typical example in which $\alpha =0.6$ and $r<1$ and $\sigma \omega >0$ are stationary, which is a partially synchronous state.

Therefore, in a large system at $\alpha =\pi $, the states of all oscillators enter a regime of “repulsive synchronization”, which is engendered by negative or repulsive coupling. This repulsive synchronization explains the mechanism in some synchronization behaviors. For instance, the synchronized flow state in the three-phase traffic theory^{17,18} occurs when moving vehicles adapt their speed to other vehicles and decelerate to avoid collision. This speed adapation mechanism corresponds to the phase frustration in the Kuramoto model, and this is more intuitive when the traffic pattern forms a circle.^{17} If an ensemble of self-driving cars set to have a strong speed adaption to maintain a safe gap between vehicles move along a circular road, a synchronized flow state emerges and all vehicles are approximately distributed in an homogeneous lane with a common angular velocity or frequency. Thus, the Kuromoto model with phase frustration explains this traffic problem. Other examples of repulsive synchronization include competitive coevolution.^{16,37} More importantly, the synchronized transition due to the negative coupling should be described by the parameter related with the observed frequency, rather than $r$.

Note that the frequency locking state with a uniformly distributed phase among the oscillators is related to the phase frustration $\alpha $ and is prevalent in many complex networks other than the regular lattice as it initially found. Figure 6 presents the synchronization transition diagram in an Erdös-Rényi (ER) random network, which is similar to that in a BA network. This result shows the survival of repulsive synchronization in disordered networks. Also, the frustration effect is reflected in Fig. 7, in which $\alpha \u2212\lambda $ diagram showing the critical forward as well as backward transition points for frequency-locking transition. Similar to the BA network, there are two synchronized regions, one is near $\alpha =0$ and the other is near $\alpha =\pi $. When $\alpha $ is in the neighborhood of 1.2, no synchronization is observed.

Furthermore, in order to find out how the disorder of the networks influences the synchronization, we conduct simulation on Watts-Strogatz (WS) small-world networks generated with different rewiring probabilities $p$. In Fig. 8, the repulsive synchronous state at $\alpha =\pi $ [shown in Figs. 8(f)–8(j)] emerges for all networks, and repulsive synchronization phase transition remains almost unchanged with different $p$’s. It demonstrates that the properties of repulsive synchronization have little relationship with the randomness or disorder of the network. Notice that in Figs. 8(a)–8(b), $r$ is very small. This is because natural frequency $\omega i$ of each node $i$ equals to its degree $ki$, and the degree distribution is more homogeneous for smaller $p$; therefore, the frequency will get locked before the phases of more oscillators become synchronized.

To examine the explosive behaviors associated with repulsive synchronization, we plot the phase transition diagram on a scale-free (SF) network, shown in Fig. 9. These SF networks have a degree distribution $P(k)\u223ck\u2212\gamma $ with $\gamma =2.2,2.6$, and $3.4$. We obtain the three left diagrams at $\alpha =0$ and the right at phase shift $\alpha =\pi $. For the sake of comparison, we measure $r$ and $\sigma \omega $ as order parameters. When $\alpha =0$ we see that when coupling strength $\lambda $ is increased the synchronization from the order parameters is $r\u21921$ and $\sigma \omega \u21920$. In addition, there is a transformation from a first-order synchronization phase transition to a second-order transition when the degree exponent $\gamma $ increases from top to bottom, which is consistent with previous findings.^{31}

When $\alpha =\pi $, we see repulsive synchronization in all SF networks irrespective of $\gamma $. As it increases the coupling strength, $\sigma \omega $ evolves from 1 to 0, revealing the ubiquity of this frequency-locked configuration. In addition, as $\gamma $ increases, the synchronization changes from a first-order to a second-order phase transition and the width of the hysteresis loop decreases. This is associated with the heterogeneity in the degree and natural frequency distribution. When the network is heterogeneous, some of the synchronous clusters of oscillators first emerge and then different clusters become entrained, leading to explosive behaviors. Note that this repulsive synchronous state is not reflected in the order parameter $r$.

The $\alpha \u2212\lambda $ diagram in Fig. 10 gives the critical forward as well as backward transition points for the frequency-locking transition. The two synchronized regions near $\alpha =0$ and $\pi $ are also presented, and with different $\alpha $, there is a transition from first-order to second-order synchronization phase transition.

## IV. CONCLUSION

In this paper, we have investigated the synchronization phase transition in complex systems with repulsive interactions among the oscillators. By tuning the phase shift $\alpha $ in the Sakaguchi–Kuramoto model, we find that a large phase shift $\alpha $ causes desynchronization. Moreover, two distinct synchronization patterns emerge when the phase frustration parameter $\alpha $ is changed. One is caused by attractive coupling when $\alpha $ nears 0, and a synchronized cluster emerges in which both phases and frequencies of the oscillators are entrained. The other frequency locking pattern is caused by repulsive coupling, which has a uniform distribution of oscillating phases. Because this repulsive synchronization, therefore, cannot be characterized by the order parameter $r$, we introduce the dispersion of effective frequencies $\sigma \Omega $ to describe it. However, since $\sigma \Omega $ is defined statistically, the classical analysis^{25} cannot be used, hence a new approach is expected in the future. The repulsive synchronization is found to exist in complex networks with different structures, such as BA, ER, scale-free networks, and small-world networks, so may be prevalent in complex networks. Since the synchronization with repulsive coupling widely existed in ecological systems and human society, our results may be helpful to understand the mechanism in these systems.

## ACKNOWLEDGMENTS

This work is supported by the National Natural Science Foundation of China (NNSFC) (Grant Nos. 71601030 and 61673086) and the Fundamental Research Funds for the Central Universities (Grant No. ZYGX2016J058). The Boston University Center for Polymer Studies is supported by the National Science Foundation (NSF) grants (Nos. PHY-1505000, CMMI-1125290, and CHE-1213217) and by the DTRA grant (No. HDTRA1-14-1-0017).

## REFERENCES

*International Symposium on Mathematical Problems in Theoretical Physics*, Lecture Notes in Physics, edited by H. Araki (Springer, Berlin), Vol. 39, p. 6023.