{'fraction_non_alphanumeric': 0.04151935021500239, 'fraction_numerical': 0.011562350692785476, 'mean_word_length': 4.449362145274668, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low dimensional behaviour of generalized Kuramoto model

We study the global bifurcations of frequency weighted Kuramoto model in low-dimension for network of fully connected oscillators. To study the effect of non-zero-centered frequency distribution, we consider two symmetric Lorentzians as an example. We derive the stability diagram of the system and show that the infinite-dimensional problem reduces to a flow in four dimensions. Using the system symmetries, it can be further reduced to two dimensions. Using this analytic framework, we obtain bifurcation boundaries of the system, which is compatible with our numeric simulations. We show that the system has three types of transitions to synchronized state for different parameters of the frequency distribution: (1) a two-step transition, representative of standing waves, (2) a continuous transition, as in the classical Kuramoto model, and (3) a first-order transition with hysteresis. Numerical simulations are also conducted to confirm analytic results.


Introduction
Synchronization is a phenomenon representing the collective behavior in populations of interacting units. Examples are the synchronous flashing of groups of fireflies 1 , oscillation of neuronal communication causing circadian rhythms 2 , chemically reacting cells 3 , Josephson junction circuits 4 , etc. A key contribution in this area was introduced by Kuramoto 5,6 . The Kuramoto model describes the system's transition to collective dynamics with the increase in coupling strength bypassing a threshold. The original Kuramoto model considers a solid interaction between dynamical units. Even though it describes the collective phenomena, it misses the interplay between structure and dynamics. However, the generalized model, inspired by natural systems, considers the correlation between structure and dynamics by considering connections as a function of the intrinsic dynamics of units. Properties of transition from disordered (asynchronous) to ordered (synchronous) state are of particular interest. Synchronization is commonly found to follow a smooth, second-order phase transition, as shown generally by Kuramoto. However, it has been shown that in some cases, there is a first-order, or discontinuous synchronization transition, i.e., an abrupt and irreversible phase transition to the synchronized state, which is called explosive synchronization (ES). The first ES was reported in scale-free network topology with a positive correlation between natural frequencies of the oscillators and their degrees 7 . It has been shown that even if there is a partial correlation between natural frequencies of nodes and their degrees, the transition would be first-order 8 . This type of transition has been studied from different perspectives in [9][10][11][12][13][14][15][16][17][18][19][20][21] . Recently, experimental evidence of explosive synchronization in electrical, chemical, and neuronal systems are reported 11,22,23 . The first-order transition causes a hysteretic behavior in the system. The origin of hysteretic behavior has been investigated in 24,25 . In particular, ES can also happen when there is a positive correlation between the coupling strengths of the oscillators and the absolute value of the difference in their natural frequencies 26 . In 27 they have obtained the exact value of the critical coupling for transition to the synchronized state. To study the transition of the system to the synchronized state, we use the Ott-Antonsen ansatz 28 to lower the dimension of the infinite-dimensional system. In this way, one can investigate the effect of frequency distribution on the evolution of the order parameter and phase transition as a function of the coupling strength.

The frequency weighted Kuramoto model
The Kuramoto model with interactions weighted by the absolute value of frequencies, considers an ensemble of N coupled limit-cycle oscillators. The phase of i'th oscillator is given by θ i , and its evolution is governed by Where {ω i }, the natural frequency of the ith oscillator, is drown from a frequency distribution. The connectivity pattern is given by A i j s that are the elements of adjacency matrix, representing the topology of interacting network, with A i j = 1 for connected nodes and A i j = 0 otherwise. Therefore, k i = ∑ j A i j is the degree of node i. λ stands for the strength of coupling. The complex order parameter, as a measure of synchrony, can be defined as: The significant characteristic of this model is the frequency-weighted coupling, which causes a positive correlation between the coupling strength of oscillators and the absolute value of their natural frequencies, providing a heterogeneous interaction between units. This kind of coupling is characteristic of some natural and synthesis systems, such as power grids and social interactions. For example, introverted people less frequently interact with the rest of society. Hence, the rate of adaptation is lower for them.

Continuum limit
In the limit of the infinite number of oscillators (N → ∞), for each natural frequency ω, there is a continuum of oscillators distributed around the unit circle. We define this distribution with ρ(θ , ω,t), so that ρ(θ , ω,t)dθ is the fraction of oscillators of natural frequency ω, whose phase is between θ and θ + dθ at time t. ρ(θ , ω,t) is a 2π periodic function and follows the normalization condition: and frequency distribution also is supposed to be normalized: Considering equations 3 and 4 we have: Having conservation of oscillators of frequency ω, the evolution of ρ is governed by the continuity equation: The velocity in this continuous representation can be stated as: or alternatively as: with r, complex order parameter, defined with the following non-linear integro-partial-differential equation: Which is the continuous version of Eq.(2). In the incoherent state, which corresponds to the case of small coupling strength, the oscillators are distributed around unit circle ρ 0 (θ , ω,t) = 1 2π , having their own intrinsic dynamic. By finding the point that this solution loses it's stability, one can find the critical coupling in which the system does transition to synchronized state, which is the point that the system looses it's incoherent state. To do so, we add a perturbation (η) to incoherent density of oscillators: where ε 1, and η is real, 2π-periodic, and the normalization condition Eq. (3) requires that 2π 0 η(θ , ω,t)dθ = 0. We restrict ourselves to linear analysis about the fixed point and we expand η in it's Fourier modes η(θ , ω,t) = ∑ n η n (ω,t)e inθ and η n (ω,t) = α(ω,t) n , as considered by Ott-Antonsen 28 . |α(ω,t)| ≤ 1 to avoid divergence of series. This ansatz provides a family of distribution functions that provide information about the synchronized and non-synchronized dynamic naturally. By plugging this ansatz into Eq. (10), the density can be written as: Therefore, the continuity equation (Eq. (6)) becomes: Since different frequency modes evolve independently, one can study them separately. By considering the lower order terms in the sum in both sides of the equation, we obtain: and For frequency distribution, we consider bimodal Lorentzian distribution with four poles at ω = ±ω 0 ± i∆. The distribution in partial fraction form, can be written as: To make sure that g(ω) is concave in the center and therefore can be considered bimodal, it has to fulfill the condition ω 0 /∆ > 1/ √ 3. Otherwise, the tails of the two distributions merge and they form a unimodal distribution. We consider α(ω,t) to be continuous in the complex ω plane. Evaluating r * = ∞ −∞ α(ω,t)g(ω)dω, by deforming the integration from the real ω axis to Im(ω) −→ −∞, the order parameter becomes: r = r 1 +r 2 2 , where r 1,2 (t) = α * (±ω 0 − i∆,t). By plugging these values

3/12
into Eq. (14), the system reduces to a set of two coupled complex ODEs for two sub order parameters: using α * = r 1 , ω = −ω 0 − i∆, we have: and for α * = r 2 , ω = ω 0 − i∆, we get: Which provides information about the dynamics of the two clusters in the system. Since α is the first mode of perturbation when it is non-zero, it shows deviation from the random state.

Linear stability
We analyse the linear stability of the incoherent state by perturbing the system around the fixed points. By adding perturbations ( r 1 → r 1 + dr 1 and r 2 → r 2 + dr 2 ) to Eq. (16), and keeping the linear terms we have: Which can be written in a matrix form: The eigenvalues of this matrix (Λ 1 , Λ 2 ) can be found: For the second square root to be real, it has to fulfill the condition: Which indicates the bifurcation boundary of the system at Shown with solid black lines in FIG. (1). For ∆ ω 0 λ 2 , both eigenvalues are real. Otherwise, they have imaginary parts as well. Since it is symmetric around zero, we just plot the upper half plane. In the next section we will study the dynamic of the flow field in these two regions.

Reduction of the system
We define complex order parameters in polar coordinate, and a phase difference between two sub-order parameters as: Plugging in Eqs. (16) and comparing real parts in both sides of the equation, we get a differential equation for ρ, the radius of the deviation from the random state, and by comparing imaginary parts, we get a differential equation for the evolution of phase (φ ) of each cluster. Using the expressions obtained forφ 2 andφ 1 we can find an expression forψ. So, we have a set of three 4/12 differential equations for the amplitude and phase difference of suborder parameters.
We look for solutions that satisfy the symmetry condition ρ 1 (t) = ρ 2 (t) = ρ(t). Considering this condition, the dimension of the system reduces to two.
By defining: q = ρ 2 thenq = 2ρρ, and multiplying both sides of equation forρ with 2ρ The first equation becomes: We end up with a two dimensional phase plane: We use equations in (26) to plot dynamical flow in the q − ψ plane. In FIG. 1, the colored region shows the parameter values (∆, ω 0 ), which have the same type of dynamic, same for the white region. In this case, zero is a saddle node in the colored regions and, in the white regions, an unstable fixed point. In shadowed regions, both Λ 1 and Λ 2 are real, and zero is a saddle node. While, in the white regions, they have imaginary parts as well, leading to an unstable fixed point at zero. Different colors represent different values of λ . The color green stands for λ = 2.1, the color blue for the case of λ = 2.4, and red for λ = 3. There is a good match between the solid line, obtained from Eq. (21), and the boundary of the shaded region obtained from flow field characteristic using Eq. (26). The little mismatch is due to the higher-order terms, which we have ignored.

5/12
Using the trigonometric relation sin 2 + cos 2 = 1, one obtains a third order polynomial equation for fixed points of q: Using the characteristic equation of the third order polynomial δ = 18abcd − 4b 3 d + b 2 c 2 − 4ac 3 − 27a 2 d 2 , which gives information about the roots of the polynomial, we can find the bifurcation in the system.
which changes sign when: and simplifies to: Meaning that ω 0 /∆ = √ 3λ 2 − 1 is the bifurcation line of the system, and by crossing this line, the dynamic of the system changes.

Order parameter
By simulating the generalized Kuramoto model, Eq. (1), we investigate the effect of the correlation between coupling constant and intrinsic dynamics of each node in the transition to the synchronized state. We consider a fully connected network. To understand the effect of intrinsic frequencies, we take them from a bimodal Lorentzian distribution; by fixing the mean of the distributions (ω 0 = 2), we investigate the impact of the width of the distributions (∆) on transition to synchronized state.
We compute the stationary value of order parameter r, Eq. 2, by progressively changing the coupling strength in an adiabatic way and recording the average value at a time window in the steady state. The forward path is obtained by increasing the coupling strength, using the outcome of the previous step as the initial condition. The same progress is performed for backward continuation by decreasing the amount of coupling strength. In all cases, the size of the network is considered to be N = 1000, and initial phases are taken randomly between 0 and 2π. In FIG. 2, we show the results for different values of ∆. For narrow distributions (∆ < 1), the system has a two-step transition to the synchronized state FIG. 2.(a). In this intermediate step, there are two giant synchronous clusters formed in the system. Inside each cluster, the oscillators are synchronized with each other, but two clusters are not in sync with each other. One of them has a center frequency ω 0 and the other one −ω 0 . So they form standing waves in the system. This region remains stable even with increase in the λ . At the final transition (λ 2), these two synchronous clusters join together and make the whole population synchronized.
By simulating for different values of the ∆, we have observed that the smaller the width of the distribution, ∆, the larger this stable region is. By getting closer to ∆ = 1, the size of this intermediate region shrinks, and close to ∆ = 1, the system shows a continuous transition to the synchronized state, FIG. 2.(b), similar to the classical Kuramoto model. The reason for such a behaviour is that for small ∆, two clusters are far apart from each other and focused around their natural frequencies ω 0 and −ω 0 . The narrower the distribution, the larger coupling is required to combine both clusters. By further increasing the width of the distribution, we observe a first-order transition with hysteresis to the synchronized state for ∆ > 1 (FIG. 2.(c)). The existence of this hysteresis is indicative of memory in the system. Having a first-order transition to the synchronized state and a hysteresis loop is representation of the explosive synchronization. Increasing the width of the distribution ∆, causes the forward transition to happen in larger λ , and causing the system to have a larger bistable region FIG. 2.(d). When there is a hysteretic transition (∆ > 1), the backward transition, shown with green, is always happening at λ * = 2. While the forward transition, in blue, happens in larger λ for wider distribution. The reason is that in the backward transition, the dynamic starts from a synchronized state, so the width of the distribution does not play a role.
In FIG. 3, we plot the critical value of the transition (λ * ) in the parameter space. To do so, we simulate the system for values of ω 0 , and ∆; the same procedure as FIG. 2. For simplicity, we fix ω 0 = 2, and obtain the critical value of the transition (λ * ) for different values of ∆. Connecting these values gives the solid blue line for the forward critical value and the green lines for the backward transition. The light green shows the initial step of the transition, which leads to the intermediate region. The dark green shows the final step in the transition to the synchronized state. The region between them is the region of the standing wave. In between the blue and green lines, there is a region that corresponds to the hysteresis behavior in which the system shows bistable dynamic for ∆ > 1. We use the Eq. (26) again to plot the flow field in the q − ψ plane for a given value of parameters. The parameters 7/12 correspond to the value at the center of the circle for each flow field plot. We observe different types of flow field characteristics in different regions of the parameter space. In all cases, zero is an unstable node. In FIG. 3.(a),(b),(c), there is an unstable non-zero node in the second quarter, which means the synchronized state is unstable. In addition, there exists a stable limit cycle centered at zero. It indicates the fact that there are two giant dynamical clusters around the centers of the distributions. The phase difference between them continues to increase. In FIG. 3.(d), there is a fixed point solution which is representation of the synchronized state. We have observed the same type of dynamic until close to the boundaries of the magenta line. The magenta line obtained from Eq. (31), defines the boundaries of the transition. Close to this line, the zero-centered limit cycle disappears, and the fixed point becomes stable (d). Stable fixed point reflects a stable synchronized state. This region corresponds to the bistable dynamic. In the lower part, shown in FIG. 3.(e), there exists a stable fixed point for whichq = 0,ψ = 0 indicating a stable synchronized state.
We did simulation for bimodal Gaussian distribution as well and observed similar behavior.

Correlation matrix
To shed more light on the dynamic of the system, we plot the frequency distribution in the network for two types of dynamic. We select one example from a two-step transition ∆ = 0.8, and one example from an explosive transition ∆ = 2. for different values of coupling strength (λ ). Given that frequencies are symmetric, in FIG.4(a), the overall frequency distribution cancels around zero. By increasing the coupling strength, the distribution becomes asymmetric at (c), which indicates one cluster becoming stronger than the other one. Finally, at (d) λ = 2.04, one cluster takes over.
In order to study the collective dynamics in more detail, we look at the correlation matrix. The local phase configurations of the system can be determined by the correlation matrix D defined as: where t s is the time needed for reaching to stationary state. The correlation matrix element D i j is a measure of coherency between the pair of oscillators at i and j positions and takes a value in the interval 1 ≤ D ≤ −1. D i j = 1 when there is a full synchrony between nodes i and j (θ i = θ j ), and D i j = −1 when they are in the anti-phase state (θ i = θ j ± π). Here we used 5 × 10 5 time steps for the stationary time t s and 3 × 10 4 time steps for the averaging window ∆t.

Correlation matrix, frequency distribution (forward path)
We  .6 (a), the system has pick around zero. The same corresponding value (λ = 0.02) in ∆ = 2, the distribution is narrower, meaning even less correlation. For narrow distribution, by increase in the coupling constant, small correlation rises close to one. By increasing the coupling strength (λ ), the tail at one becomes stronger as in (d) while the central pick becomes weaker with shifts towards one, corresponding to a shift in the correlation matrix toward red. There is no gradual increase in correlation in the wide distribution, but the system has a first-order transition from one state to other (k),(l).

Correlation matrix, frequency distribution (backrward path)
To see the difference between the dynamic in forward and backward transition, we present the frequency distribution and correlation matrix and it's probability distribution in FIG.7, 8, 9 respectively. Similar to the forward plots (FIG. 4, 6, 5). The first row corresponds to the case of ∆ = 0.8, and the second row ∆ = 2. For ∆ = 0.8, we observe three type of dynamics. For λ ≥ 2.08, (a)-(c), the system is in the high order parameter state, with large population of oscillators being in frequency-locked state which shows itself as an asymmetric distribution in frequency FIG.7 and as a pick at one in FIG. 9. with a small population being out of sync. Increasing the coupling constant, (d)-(e), the probability distribution around zero gets stronger in FIG. 9. Increasing the coupling constant, the system goes to a full unsynchronized state (d), with a symmetric distribution around zero.   In this case, the frequency distribution also has become sharp around zero, implying a symmetric distribution for the oscillators. For ∆ = 2, on the other hand, there is no intermediate state, and the system shows a sudden transition to the unsynchronized   state. At the transition point (h) to (g), there is a sudden change of the distribution from a pick at one to a distribution around zero.

Summary and concluding remarks
In conclusion, a low dimensional description of the generalized Kuramoto model has been studied analytically and checked numerically. We have demonstrated that the Ott-Antonsen ansatz is valid for the generalized Kuramoto model. It gives a reasonable estimate of the separation of the dynamic in the parameter space. This means that except for the bistable region, all the attractors of the infinite-dimensional system lie in the low-dimensional manifold obtained by this ansatz. Our numeric simulations show that the transition to the synchronized state can be a two-step transition, a continuous transition, or an explosive transition with hysteresis, depending on the width of the distribution. Our result is also compatible with analytic results obtained for the exact critical coupling of the transition 27 , in case that there is hysteresis.