Multiple-parameter bifurcation analysis in a Kuramoto model with time delay and distributed shear

In this paper, time delay effect and distributed shear are considered in the Kuramoto model. On the Ott-Antonsen's manifold, through analyzing the associated characteristic equation of the reduced functional differential equation, the stability boundary of the incoherent state is derived in multiple-parameter space. Moreover, very rich dynamical behavior such as stability switches inducing synchronization switches can occur in this equation. With the loss of stability, Hopf bifurcating coherent states arise, and the criticality of Hopf bifurcations is determined by applying the normal form theory and the center manifold theorem. On one hand, theoretical analysis indicates that the width of shear distribution and time delay can both eliminate the synchronization then lead the Kuramoto model to incoherence. On the other, time delay can induce several coexisting coherent states. Finally, some numerical simulations are given to support the obtained results where several bifurcation diagrams are drawn, and the effect of time delay and shear is discussed.


I. INTRODUCTION
The Kuramoto model was first proposed in [1,2], consisting of a group of phase oscillators.
It is now a mathematical model used most to describe synchronized phenomena [3]. More specifically, it models the behavior of the phase of a large set of weakly coupled, near identical oscillators [4][5][6][7]. Its formulation was motivated by the behavior of systems of chemical and biological oscillators, and it has been found widespread applications such as in neuroscience [8][9][10], oscillating flame dynamics [11,12], and some physical systems, such as Josphson junctions [13].
The nonisochronicity, or shear we say, is an important factor in forming spatial-temporal patterns in oscillatory system, which quantifies the dependence of the frequency on the amplitude of oscillations, and can also induce very rich dynamical behaviors in ensembles of identical or near identical oscillators [14,15]. Montbrio and Pazo [16] have analyzed the Kuramoto model with distributed shear and showed that the onset of collective synchronization is impossible if the width of the shear distribution exceeds a precise threshold. Kuramoto model with shear was also studied in [17,18] and the references there in.
Time delay is another important factor in coupled systems. In the real world, every items always takes some time lag to interact with others. As we all know, time delay is ubiquitous in many fields, such as the control theory, ecology, mechanics, management science, physics, neural, etc [20,21,37]. Recently, there also has been extensive interest in the delay effect on the dynamical behavior of the Kuramoto model [10,23,24,[26][27][28][29][30]. So far, Kuramoto model with shear and time delay effect has not been well studied yet.
Motivated by such considerations, we combine time delay effect and shear in a Kuramoto system with the following forṁ where ω i 's are randomly distributed frequencies, K is the coupling strength, and q j 's are randomly distributed shears which characterize the impact of the distance that oscillators move away from the unit circle on the frequencies of the oscillators themselves.
As stated in [16], the synchronization transition in (1) fails in case of large spread of shear distribution in the absence of time delay. Also as known to all, time delay will also induce synchronization transition in Kuramoto model [27,28]. Thus in this paper we will study the effect of delay together with shear on the synchronization transition from the bifurcation approach. The idea is followed by [19] based on the Ott-Antonsen's manifold reduction [32,33]. Kuramoto model (1) can be reduced into a functional differential equation on a submanifold of the phase space. Trivial equilibrium stands for incoherent state whereas nontrivial periodic oscillations indicate synchronized oscillations of the Kuramoto model.
The transition can be described by Hopf bifurcations which will be discussed in the parameter space consisting of delay and shear.
Through some rigorous bifurcation analysis, we find that time delay can have very important impact on the system dynamics: compared with the previous results in [16][17][18], time delay will further prevent synchronization in some cases and will also induce several coexisting coherent states. The bifurcation behavior in multiple parameters spaces are also discussed. These bifurcation phenomena are investigated both theoretically and numerically in this paper.
The rest of this paper will be organized as follows. In Section 2, we shall use the Ott-Antonsen method to reduce the original Kuramoto model and obtain the corresponding functional differential equation. In Section 3, we shall consider the stability of the trivial equilibrium standing for the incoherence of Kuramoto model and the existence of the local Hopf bifurcation. In Section 4, the stability and direction of periodic solutions bifurcating from Hopf bifurcations are investigated by using the normal form theory and the center manifold theorem due to [31], which gives clearly the location where coherent states appear and their stability. In Section 5, inspired by the method given in [19,24,25], numerical simulations are carried out to support the obtained results.

II. OTT-ANTONSEN'S MANIFOLD REDUCTION
In this section, Ott-Antonsen's method is employed to study the delayed Kuramoto model with shear. The approach is actually a direct extension of the method in [16-18, 32, 33] thus we only give some key steps. To analyze model (1) we adopt its thermodynamic limit N → ∞, then drop the indices and introduce the probability density for the macroscopical phases f (θ, ω, q, t) [34], which represents the ratio of oscillators with phases between θ and θ + dθ, natural frequencies between ω and ω + dω, and shear between q and q + dq. Hence the density f must obey the continuity equation where c.c. stands for complex conjugate of the preceding term, and the complex-valued order parameter is If the phases are uniformly distributed, we know r(t) vanishes to zero. This state is customarily referred to as incoherence or incoherent state. Similarly |r(t)| = 1 stands for the completely synchronized state, and |r(t)| ∈ (0, 1) partially synchronized states or coherent state. Since f (θ, ω, q, t) is a real function and 2π-periodic in the θ variable, we know it admits the Fourier expansion where f l = f * −l , f 0 = 1, the " * " stands for the complex conjugate, and p(ω, q) is the joint probability density function of ω and q. The first order Fourier term f 1 is important because it determines the order parameter Inserting the Fourier series (4) into the continuity equation we obtain the following equation which governs the dynamics of Kuramoto model. Ott and Antonsen [32,33] found that the ansatz is a particular and important solution of the Kuramoto model which can eliminate the ω-variable without loss of generality.
In this paper, we further assume the frequency and the shear are independent and both obey the Lorentzian distribution with mean ω 0 and q 0 , and the width δ and γ, respectively.
This is a delay differential equation, which characterizes the dynamical behavior of (1), as proved in [32,33]. Investigating stability conditions for the incoherent state r(t) = 0 and exploring the existence of Hopf bifurcations will be used to detect the appearence of coherent states.

III. STABILITY AND BIFURCATION ANALYSIS
The characteristic equation of the linearization of (8), around the trivial equilibrium (incoherence), is, Based on the stability theory of the functional differential equations, if all the roots of Eq. (9) have negative real parts, then the zero solution of system (8) is always asymptotically stable, which means that the Kuramoto model will remain in a incoherent state in the sense of small perturbations. If system (8) has an orbitally asymptotically stable periodic solution which is bifurcated from Hopf bifurcation after Eq.(9) has purely imaginary roots, then the Kuramoto model will exhibit partial synchronization (coherent states). This Hopf bifurcation gives the critical value of the synchronization transition of the Kuramoto model.
From the discussion above we know that, under the assumption 2) has a purely imaginary root when τ = τ ± j , respectively. We reorder {τ + j } {τ − j } as {τ j } so that τ 0 < τ 1 < τ 2 < . . . . Clearly, τ 0 = min{τ + 0 , τ − 0 }. In order to state the stability and Hopf bifurcation results of Eq.(9) clearly, we make the following assumptions In fact, when investigating the effect of shear, assumptions (P1) can be cast into a different Thus, for fixed γ, (P1) means that we are using a small q 0 , the mean value of shear, or a small K, the coupling strength, or a large δ, the spread of the frequency distribution, whereas (P2) means the contrary. Furthermore, if (P1) holds, so does (P4). Similarly if (P3) holds true, so does (P2). Thus three cases are considered in the following part: (P1) holds, (P3) holds, and (P2,P4) holds, respectively. In the coming part, we will show that these three cases correspond to that the Kuramoto model exhibits incoherent states, that the Kuramoto model is synchronized when delay exceeds a threshold, and that the model exhibits synchronization windows, i.e., incoherent and coherent states are interlacing.

IV. DIRECTION AND STABILITY OF HOPF BIFURCATION
In Section 3, we have obtained the stability of the incoherent state and a group of conditions which guarantee that the equation undergoes Hopf bifurcation at some critical values of τ . In this section, we shall study the direction and stability of the bifurcating periodic solutions. The method we used is based on the normal form method and the center manifold theory presented in Hassard et al [31]. We shall compute the center manifold construction of system (8) at τ =τ ∈ {τ ± j , j = 0, 1, . . .} while the characteristic equation has a root iβ 0 . The following variable can be calculated as shown in the Appendix: Based on [31] and [38], we have the following results.
Remark 2 Assume the incoherence is stable when τ = 0, and Hopf bifurcation occurs atτ . By using the global Hopf bifurcation results (see [19,22]), we know there will always exist hysteresis loop near the subcritical bifurcations. This will be shown in the simulation part such as Figure 3(a).

V. NUMERICAL EXAMPLES
In this section, with the aid of theoretical results obtained in the previous sections, we will carry out several groups of numerical simulations to illustrate the complicated dynamical behavior of the Kuramoto model with time delay and shear. All the illustrations we will show are about the bifurcation branches of the model, because bifurcation branches are widely used in literatures [19,24,25], which are also a quite unambiguous way to show the change of numbers of steady states or periodic oscillations.

A. Single-parameter bifurcations
We first give some numerical simulations about the reduced model (8)  increases, |r(∞)| first reduces to 0, which means the system becomes incoherent, then the system switches between the coherent state and the incoherent state, i.e., the synchronization switches appear. We find it supports results in the Theorem 1(2)(b). In Figure 2 Similarly, if we fix τ = 2.5 and let K vary, we obtain simulation results in Figure 3. We find that the system (8) exhibits incoherent state when K is low and it exhibits partially coherent state when K exceeds a precise threshold, which supports the obtained result in literature [16]. Moreover, if we further increase k more than one branches of bifurcating solutions appear which may also be originated from a backward bifurcation (shown in Figure   3(a) with red circles). It is well known that the system exhibits globally stable incoherent state when K = 0, thus all branches of bifurcating solutions exists globally for K → ∞.
This means that more than one branches of coherent states coexists (with similar order parameters shown in Figure 3(a) with black circles) but have different periods (shown in Figure 3 (b)).
If we fix τ = 2 and 4.4 respectively, and let γ vary, we obtain simulation results in Figure   4. We find that the system (8) exhibits partially coherent state when γ is low and it exhibits incoherent state when γ exceeds a precise threshold. This indicates that large spread of shear will significantly eliminate the synchronization of the Kuramoto model, which is an extension of the results given in [16,17] in the absence of delay. Larger τ induces more branches of coherent states as shown in these simulations, which may also be subcritical branches (red circles in Figure 4 (c)).
If we choose K = 1, ω 0 = 3, δ = 0.1, γ = 1, τ = 1, and let q 0 vary, we obtain simulation results in Figure 5(a-b). We find that the system (8) exhibits incoherent state when |q 0 | is low and it exhibits partially coherent state when |q 0 | exceeds a precise threshold. This indicates that large value of average shear induces partial synchronization in the Kuramoto model. Using larger τ = 4, we give similar results in Figure 5(c-d), where we find this also brings more branches of coherent states. The two branches of coherent states are both stable when γ is away from the bifurcation points, e.g., there exists one branch marked by red circles in Figure 5(c) which turns to be stable for small γ.

B. Two-parameter bifurcations
In this section, we will show the effect of shear and delay in two-parameter plane.
When ω 0 = 3, δ = 0.1, q 0 = 0.5 and γ = 0.5, we draw the Hopf bifurcation values by curves shown in Figure 6(a). The red curves stands for τ + j and the blue curves stands for τ − j . We can find that system (8) is incoherent when K is less than the critical value and system is in partial synchronization when K crosses the Hopf bifurcation values, which means the coupling strength enhances the synchronization. Also, the synchronization switches are also observed when varying time delay τ .
When K = 1, ω 0 = 3, δ = 0.1 and q 0 = 0.5, we draw the Hopf bifurcation values by curves shown in Figure 6(b). The red curves stands for τ + j and the blue curves stands for τ − j . We can find that system (8) is in partial synchronization when γ is less than the critical Stable periodic orbits are labeled by black circles, unstable (1/3 Floquet exponents with positive real part) periodic orbits are labeled by blue/red circles, respectively.
value and system is incoherent when γ crosses the Hopf bifurcation values, which means the large spread width of the shear can weaken the synchronization. In this figure, we can clearly find that increasing τ leads to more branches of Hopf bifurcations which is also a theoretical explanation of Figure 4(a) and 4(c).
When K = 1, ω 0 = 3, δ = 0.1 and γ = 1, we draw the Hopf bifurcation values by curves shown in Figure 6 (c). The red curves stands for τ + j and the blue curves stands for τ − j . We can find that system (8) is incoherent when |q 0 | is less than the critical value and system is in partial synchronization when |q 0 | crosses the Hopf bifurcation values, which means the absolute value of the mean of shear can strengthen the synchronization. Again, as shown in this figure, we find increasing τ leads to more branches of Hopf bifurcations.

C. Three-parameter bifurcations
When ω 0 = 3, δ = 0.1 and q 0 = 0.5, we draw the Hopf bifurcation values by surfaces shown in Figure 7. In this figure, we combine the previous one-parameter or two-parameter bifurcation results. As discussed in the previous section, we find that the incoherent state loses its stability and coherent states appear, when parameters crosses the surfaces along the direction that γ decreases. We can find that system (8)  We found that increasing time delays would lead the Kuramoto model to synchronization switching. On one hand, time delay may induce synchronization thus lead the system to order. On the other, time delay can eliminate synchronization thus lead the system to disorder (see Theorem 1). The effect of mean value and variance of the shears on the dynamics are also obtained. As expected, decreasing the variance will lead the system to coherent states. Also, increasing the mean value will lead the system to coherent states.
Because of the existence of time delay, there may exist more than one branches of bifurcated coherent states. Some local and global bifurcation diagram are given in multiple-parameter spaces, which indicate large time delay and suitable shears could induces several stable coherent states in Kuramoto model simultaneously.