Twisted states in nonlocally coupled phase oscillators with frequency distribution consisting of two Lorentzian distributions with the same mean frequency and different widths

In globally coupled phase oscillators, the distribution of natural frequency has strong effects on both synchronization transition and synchronous dynamics. In this work, we study a ring of nonlocally coupled phase oscillators with the frequency distribution made up of two Lorentzians with the same center frequency but with different half widths. Using the Ott-Antonsen ansatz, we derive a reduced model in the continuum limit. Based on the reduced model, we analyze the stability of the incoherent state and find the existence of multiple stability islands for the incoherent state depending on the parameters. Furthermore, we numerically simulate the reduced model and find a large number of twisted states resulting from the instabilities of the incoherent state with respect to different spatial modes. For some winding numbers, the stability region of the corresponding twisted state consists of two disjoint parameter regions, one for the intermediate coupling strength and the other for the strong coupling strength.


Introduction
It is well known that natural frequency distribution g(ω) plays a critical role in displaying rich synchronous dynamics in globally coupled phase oscillators. For unimodal frequency distribution, the partial synchronous state steps in through a continuous transition when the coupling strength is above a critical coupling strength [1]. Further increasing the coupling strength, the partial synchronous state may transit to a global synchronization where all oscillators oscillate at a same frequency. For a bimodal frequency distribution, increasing coupling strength always first leads to a standing wave state and then to a traveling wave state [2]. When the peak distance in the bimodal frequency distribution is narrow, the discontinuous transitions between different dynamical states are possible. When the frequency distribution becomes more complicated, for example a trimodal one, the synchronous dynamics in globally coupled phase oscillators may display collective chaos through a cascade of period-doubling bifurcations [3]. In the above works, incoherent state is always unstable when the coupling strength is above a threshold and further increasing coupling strength always enhances synchronization among oscillators. Recently, some authors studied the synchronous dynamics when the frequency distribution is a superposition of two unimodal frequency distributions with the same mean frequency in a globally coupled oscillator system [4]. They found a non-universal synchronization transition in which the incoherent state may be restored for the coupling strength above the threshold.
The impacts of the frequency distribution on the synchronous dynamics possibly exist in other types of coupling topologies, for example in nonlocally coupled phase oscillators. A system of nonlocally coupled phase oscillators has become an important platform for investigating collective dynamics among oscillators. One interesting dynamics is the chimera state in which synchronous clusters and desynchronous clusters coexist [5][6][7][8][9][10][11][12]. Twisted state is another interesting type of collective dynamics and is a simpler one [13][14][15][16][17]. Twisted state was first identified in a ring of coupled identical phase oscillators and is characterized by the linear dependence of the phase of oscillator on oscillator's spatial position. A twisted state with winding number q refers to a state in which the phase of oscillator increases linearly along the ring and varies 2πq through the ring. Wiley et al. found coexistence of a large number of twisted states in nonlocally coupled identical phase oscillators and the probability finding a twisted state with winding number q follows a Gaussian distribution with respect to q [13]. Lately, the investigation on twisted states was generalized to non-identical phase oscillators with unimodal frequency distribution and the concept of partially coherent twisted states was proposed [15]. Beyond a critical coupling strength, the incoherent state becomes unstable and partially coherent twisted states step in. In partially coherent twisted states, the existence of asynchronous oscillators breaks the linear dependence of oscillator's phase on position. Nevertheless, by introducing a spatially-dependent complex order parameter, it can be found that the argument of the order parameter varies 2qπ through the ring in a partially coherent twisted state. Recently, we studied nonlocally coupled phase oscillators with bimodal frequency distribution and found twisted standing waves besides twisted traveling waves [17]. In these works on twisted states, incoherent state remains unstable once twisted states appear and the coherence among oscillators, characterized by the amplitude of spatially-dependent order parameter, always increases with the coupling strength. Now, it is interesting to ask whether the superposition of two unimodal frequency distributions may dispute these observations and how it influences the synchronous dynamics in nonlocally coupled phase oscillators.
In this work, we investigate the effects of the superposition of two unimodal frequency distributions on the incoherent state and twisted states in a ring of nonlocally coupled phase oscillators. We find that, under proper parameters, increasing the coupling strength may encounter multiple stability islands of the incoherent state. We also find that, for some winding numbers, the stability region of the corresponding twisted state consists of two disjoint parameter regions, one for the intermediate coupling strength and the other for the strong coupling strength, and all of these twisted states are related to the instabilities of the incoherent state with respect to different spatial modes. Moreover, we discover the existence of complicated spatiotemporal patterns which have not been observed in previous works mentioned above.
The paper is organized as follows. Firstly, we will describe the model and derive a reduced model in the continuum limit by using the Ott-Antonsen (OA) ansatz. Then, we theoretically investigate the stability of the incoherent state based on the reduced model and we present numerical results on different types of twisted states. The connection between different twisted states and the instabilities of the incoherent state with respect to different spatial modes is proposed. Finally, we will present a brief summary.

Materials and methods
We consider N phase oscillators which sit evenly on a ring with the length 2π and nonlocally interact with each other with the coupling strength K. The model is described as The subscript j denotes the unit index, which has to be taken modulo N. M denotes the coupling range and, in the following, we use the coupling radius σ = M/N as a controlling parameter and 0 < σ < 0.5. The phase lag α is in the range of (0, π/2). ω j is the natural frequency of oscillator j located at x j = 2πj/N, which follows the frequency distribution g(ω). In this work, we let g(ω) be a superposition of two Lorentzian distributions with p in the range of [0, 1]. γ 1,2 are the half widths of the two Lorentzians. When p = 0 or p = 1, the frequency distribution reduces to a normal Lorentzian one. Unless specified, γ 1 = 1, γ 2 = 0.01, and p = 0.9.
Partially coherent twisted states can be characterized by the spatially-dependent complex order parameter Z j ¼ R j e iY j which is defined as When R j becomes nonzero, partially coherent states step in. On the other hand, the dependence of Θ j on oscillator's position distinguishes different twisted states.
It has been shown in Refs. [15,17] that partially coherent twisted states and their stabilities can be theoretically investigated in the continuum limit based on the OA ansatz [18,19]. In the limit N ! 1, we introduce the probability density function f such that f(ω, x, θ, t)dθdω is the fraction of oscillators with phases between θ and θ + dθ and natural frequencies between ω and ω + dω at time t and position x. The spatially-dependent complex order parameter Z(x, t) is reformulated as with The evolution of the density f(ω, x, θ, t) obeys the continuity equation Using the OA ansatz [18,19], we have f ðo; x; y; tÞ ¼ gðoÞ ½aðo; x; tÞ� n e iny þ c:c: where c.c. is the complex conjugate of the previous term. By substituting Eq (7) into Eq (6) and following the procedures in Refs. [15,17], we have , measuring the coherence in the subpopulations of oscillators with natural frequency following Gðx À yÞ ½pu 1 ðy; tÞ þ ð1 À pÞu 2 ðy; tÞ�dy.

Results and discussion
The dynamics in Eq (1) can be reproduced by Eq (8). Therefore, we study the stability of the incoherent state in the model (1) by linearizing Eq (8) around it and the synchronous dynamics in the model (1) by performing numerical simulations on Eq (8). Traditionally, the coherence in Eqs (1) and (8) is described by the global complex order parameter defined by However, R cannot provide correct information on twisted states with nonzero q. Considering that, in Eq (8), u 1,2 (x, t) is driven by spatiallydependent complex order parameter Z(x, t), we measure the coherence in Eq (8) by another global quantity h|Z|i, defined as hjZji ¼ R 2p 0 jZðx; tÞjdx=2p. h|Z|i 6 ¼ 0 indicates a synchronous state while h|Z|i = 0 the incoherent state.

Incoherent state
We start with the stability of the incoherent state. The incoherent state is characterized by u 1,2 (x, t) = 0, and hence by Z(x, t) = 0. Linearizing Eq (8) around the incoherent state, we have the evolution of perturbation δu 1,2 with dZ ¼ R 2p 0 Gðx À yÞ½pdu 1 ðy; tÞ þ ð1 À pÞdu 2 ðy; tÞ�dy. Supposing that the perturbation takes the form δu 1,2 (x, t) = v 1,2 e iqx e λ q t with the wave number q an arbitrary integer (that is, the perturbation on the spatial mode e iqx ), we find the spectrum equation whereĜðqÞ ¼ sin ð2psqÞ=2psq is the Fourier transform of the coupling function G. For any q, there are two solutions to Eq (10). The one with the larger Re(λ q ), termed as the growth rate of the perturbation with the wave number q, is relevant for the stability of the incoherent state. For given parameter sets, the incoherent state is stable if and only if the growth rates Re(λ q ) for all possible q are non-positive. For convenience, we first treat q in Eq (10) as a real number. Fig 1 shows the growth rate Re(λ q ) against real number q for several coupling strengths and coupling radius with other parameters fixed. With the increase of |q|, the growth rate Re(λ) displays a trend of damping oscillation whose period is determined by the joint parameter σq.
Depending on the coupling strength, the maximum of the growth rate may occur at either q m = 0 or a nonzero K-dependent q m . That is, the perturbation taking the form of the spatial modes with the wave numbers close to q m are the most dangerous to the stability of the incoherent state. The right insets in Fig 1 identifying q m against K show nonzero q m for intermediate coupling strength K and zero q m otherwise. It is interesting to note that the intermediate coupling strength occupies the range of (1, 2.5) and is insensitive to the coupling radius σ.
In the range of intermediate coupling strength K, q m increases with the coupling strength. As we have mentioned, the wave number q is treated as a real number in Eq (10) for the convenience of solving the equation. However, a periodic boundary condition with respect to x in the model (1) requires q to be integers, therefore, the right insets in Fig 1 show that the most dangerous spatial modes could be q m = 0, 1 for σ = 0.4, q m = 0, 1, 2 for σ = 0.2, while q m � 20 for σ = 0.02. In other words, decreasing the coupling radius leads that more spatial modes should be considered for the stability of the incoherent state. Besides, we present Re(λ 0 ), Re(λ 1 ) and Re(λ 2 ) against K by the left insets in Fig 1, from which we can determine the stability of the corresponding partially coherent twisted states. Then we present the stability diagram of the incoherent state in the plane of the phase lag α and the coupling strength K. We take σ = 0.4 as an example where the spatial modes q = 0 and q = 1 should be considered. The critical curves with the growth rates Re(λ 0 ) = 0 and Re(λ 1 ) = 0 are plotted in Fig 2(a) and the intersection of the parameter regimes on the right side of critical curves gives the stability regime of the incoherent state. To be addressed, the critical curve of Re(λ 0 ), independent of the coupling radius σ, actually gives rise to the boundary of the stability regime for the incoherent state in globally coupled phase oscillators where σ = 0.5. The Sshaped critical curve Re(λ 0 ) = 0 suggests that, in a certain range of α, there is another stability island for the incoherent state in the globally coupled oscillators besides the coupling strength K close to zero [4]. The resurgence of the incoherent state exists when σ = 0.4 (the nonlocally coupled oscillators), which is exemplified by the inset where h|Z|i against K at α = 0.9 suggests one stability island at intermediate coupling strength. To be noted, with the assistance of the spatial mode q m = 1, there may be two stability islands for the incoherent state with the increase of α at σ = 0.4. With the decrease of σ, more spatial modes endangering the incoherent state have to be considered. Consequently, more stability islands for the incoherent state appear and the range of α supporting the non-universal synchronization transitions narrows [see Fig 2(b) and 2(c)]. In the limit σ ! 0, the stability island of the incoherent state is lost and the forward critical coupling strength for synchronization transition undergoes a jump at a certain α, which is hinted in Fig 2(d) where σ = 0.02. It is worth mentioning that the stability

Partially coherent twisted state
When the incoherent state becomes unstable, the synchronous states pop up. We investigate the transition scenario from the incoherence to synchronous behaviors with the coupling strength K at a typical phase lag, for example α = 0.8, with other parameters the same as those in Fig 2. For each coupling strength, we consider several realizations, each with initial conditions taking the form of u 1,2 (x, t)*e iqx with a specific integer q and different ones with different q. That is, we use different initial conditions for each fixed K. Fig 3 presents h|Z|i against the coupling strength K for different coupling radius σ at α = 0.8 and, for better illustrations on the behavior of h|Z|i, we zoom in two regions: K 2 (0, 2.5) and K 2 (2.4, 2.9). As shown in  To identify these different synchronous states and understand the dependence of h|Z|i on K, we take σ = 0.2 as an example. Following the line (shown by the dashed line α = 0.8 in magenta) in Fig 2(b) when K increases from zero, we encounter several critical coupling strengths, from K A to K G , at which one of Re(λ q ) becomes zero and signals the alternation between the stability and the instability of the incoherent state in response to the perturbation along the spatial modes e iqx . When K 2 (K A , K D ), Re(λ 0 ) becomes positive and the incoherent state is unstable with respect to the homogeneous perturbation. As a result, in Fig 3(b1), h|Z|i rises from zero at K A and falls to zero at K D . Fig 4(a) where K = 0.5 in the range of K A to K D shows a partially coherent twisted state with q = 0 characterized by ϕ 1,2 (x, t) = ϕ 1,2 (t), a homogeneous time-periodic argument, and |u 1,2 (x, t)| = |u 1,2 |, an amplitude independent of time and space. In this twisted state, |u 2 | � |u 1 | with |u 1 | close to zero and two global quantities, R and h|Z|i, are almost the same. When K 2 (K B , K E ), Re(λ 1 ) becomes positive, which gives rise to a twisted state with q = 1. In this range, h|Z|i becomes nonzero in response to the perturbation taking the form of e ix . Fig 4(b) shows a twisted state with q = 1 at K = 1 where ϕ 1,2 (x, t) varies 2π along the ring and |u 1,2 (x, t)| = |u 1,2 |. In this q = 1 twisted state, R ' 0 due to the fact that the contributions from the locations x and x + L/2 are cancelled out since The twisted state with q = 0 and |u 1 (x, t)| far away from 0 at K = 3. (f) The twisted state with q = 1 and |u 1 (x, t)| far away from 0 at K = 6. The first column shows the evolutions of the global quantities R (in green) and h|Z|i (in blue). The second column shows the evolutions of ϕ 1 (x = 0.1, t) (in black) and ϕ 2 (x = 0.1, t) (in red) and the third column shows |u 1 (x = 0.1, t)| (in black) and |u 2 (x = 0.1, t)| (in red). The fourth column shows the snapshots of ϕ 1 (x, t) (in black) and ϕ 2 (x, t) (in red) and the fifth column the snapshots of |u 1 (x, t)| (in black) and |u 2 (x, t)| (in red). γ 1 = 1, γ 2 = 0.01, α = 0.8, and p = 0.9. e iϕ 1,2 (x, t) + e iϕ 1,2 (x + L/2, t) = 0. Furthermore, Re(λ 2 ) becomes positive from negative at K C , which may be confirmed in Fig 3(b1), and induces a twisted state with q = 2 presented in Fig 4(c) where ϕ 1,2 change 4π along the ring. Re(λ 1 ) becomes positive again when the coupling strength is increased to K = K G . However, the new q = 1 twisted state [see Fig 4(f)] pops up till K ' 5.4. It is because that the new q = 1 twisted state may exist since K = K G but it becomes stable only till a larger K, which can be proved by the stability analysis of the twisted state following Ref. [15]. To be stressed, Re(λ 2 )>0 is kept even when K = 5. However, the twisted state with q = 2 may undergo a secondary bifurcation and leads to a modulated twisted state. In the modulated twisted state, as shown in Fig 4(d) where K = 2.55, |u 1,2 | become space-dependent with ϕ 2 varying 4π along the ring and ϕ 1 undergoing a small amplitude oscillation with spatial period L/2. More interestingly, |u 1,2 (x, t)| become periodic in time in this modulated twisted state. Actually, the bifurcation to the modulated twisted state with q = 2 occurs at K ' 2.52 and Fig 3(b2) suggests a continuous transition between the time-independent and the modulated twisted states with q = 2. When increasing K to a strong coupling strength K F , Re(λ 0 ) becomes positive again, which leads to another twisted state with q = 0. As shown in Fig 4(e) with K = 3, this new twisted state with q = 0 is characterized by ϕ 1,2 (x, t) = ϕ 1,2 (t) and |u 1,2 (x, t)| = |u 1,2 |. The new state is quite different from the q = 0 state in Fig 4(a) in the sense that, here, both |u 1 | and |u 2 | are far away from 0. Due to the frequency distribution Eq (2), we may partition phase oscillators into two groups, one with the Lorentzian distribution with half width γ 1 and the other with the Lorentzian distribution with half width γ 2 . Consider that u 1,2 (x, t) = a(−iγ 1,2 , x, t) play the order parameters measuring the coherence in groups of oscillators with the half widths γ 1,2 in the frequency distributions, the difference between these two twisted states can be apprehended as following. In the twisted state with q = 0 in Fig 4(a), |u 1 (x)'0| suggests that only oscillators in γ 2 group get partially synchronized while those in the γ 1 group are desynchronized. In contrast, in the twisted state with q = 0 in Fig 4(e), |u 1,2 | far away from zero indicate that oscillators in both the γ 2 group and the γ 1 group get partially synchronized. The similar cause results in the two different twisted states for q = 1. Noting that |u 1,2 | are far away from zero in Fig 4(f), we know that oscillators in both the γ 1 and the γ 2 groups are in partially synchronous state, which is different from the twisted state with q = 1 in Fig 4(b) where oscillators in the γ 1 group are asynchronous.
Accordingly, the transition scenarios in Fig 3(a), 3(c) and 3(d) can be analyzed in the same way. Briefly, the synchronous states presented in Figs 3 and 4 originate from the instabilities of the incoherent state to the perturbations in different spatial modes. Based on Fig 2, the multistability in Eqs (1) and (8) may be observed when there exist several spatial modes with positive Re(λ q ) at the same parameter set, and it could be understood that decreasing σ leads to more and more synchronous states. In addition, it has to be addressed that there are two types of twisted states with the same winding number q, one with u 1 ' 0 in weak and intermediate coupling strength and the other with u 1 far away from 0 in strong coupling strength. In the former one, only oscillators in the γ 2 group get partially synchronized while, in the latter one, oscillators from both γ 1 and γ 2 groups get partially synchronized.

Complicated spatiotemporal patterns
We have presented several twisted states in correspondence to the response of the incoherent states to perturbation on different spatial modes. However, there are some parameter regions in which there are several different spatial modes unstable to the incoherent state. If we consider arbitrary initial conditions, the competition among different spatial modes may lead to complicated spatiotemporal patterns. Here, we present some of them. Fig 5(a) shows a periodic spatiotemporal pattern. The amplitudes |u 1,2 (x, t)| display characteristics of standing wave state with spatial periodicity L while ϕ 1,2 (x, t) perform periodic oscillation in the presence of some phase slips. Fig 5(b) shows an interesting spatialtemporal patterns. In the respect of | u 1,2 (x, t)|, the pattern looks like a localized structure traveling along the ring. At the same time, the space is divided into two moving regions in the respect of ϕ 1,2 (x, t) and, in different regions, ϕ 1,2 (x, t) increases at different velocities. The model also displays spatialtemporal chaos. Fig 5(c) shows a weak spatialtemporal chaos in which |u 1,2 (x, t)| display a little regularity while Fig 5(d) shows a more turbulent state for both |u 1,2 (x, t)| and ϕ 1,2 (x, t). According to the classification in a recent work [20], the states in Fig 5(c) and 5(d) belong to the amplitude turbulence since |u 1,2 (x, t)| = 0 may be encountered at any spatial location.

Conclusion
In summary, we studied a ring of nonlocally coupled phase oscillators in which the frequency distribution is made up of two Lorentzians with the same center frequency but with different half widths. Using OA ansatz, we derived a reduced model in the limit of infinite number of oscillators. Based on the reduced model, we studied analytically the stability of the incoherent 52. The first column shows the evolution of |u 1 (x, t)|, the second column shows the evolution of |u 2 (x, t)|, the third column shows the evolution of ϕ 1 (x, t), and the fourth column shows the evolution of ϕ 2 (x, t). γ 1 = 1, γ 2 = 0.01, α = 0.8, and p = 0.9. https://doi.org/10.1371/journal.pone.0213471.g005 Twisted states in nonlocally coupled phase oscillators state and found that the most unstable spatial mode to the incoherent state depends on the coupling strength, the homogeneous perturbation with the winding number q m = 0 for the weak and the strong coupling strength and inhomogeneous perturbation with nonzero winding number q m 6 ¼ 0 for intermediate coupling strength. The critical curves of the different spatial modes to the incoherent state may be S-shaped ones, which may give rise to multiple stability islands of the incoherent state. By numerically simulating the reduced equations, we found a large number of twisted states which result from the response of the incoherent state to the perturbations in the form of different spatial modes. Especially, there are two types of twisted states, one for the intermediate coupling strength and the other for the strong coupling strength. By partitioning oscillators into two groups each of which obeys a unimodal frequency distribution, we found that, for the twisted state for the intermediate coupling strength, the group of oscillators with fat frequency distribution remains in desynchronization and, for the states for strong coupling strength, both groups are in partial synchronization.