Chimera states formed via a two-level synchronization mechanism

Chimera states, which consist of coexisting synchronous and asynchronous domains in networks of coupled oscillators, are in the focus of attention for over a decade. Although chimera morphology and properties have been investigated in a number of models, the mechanism responsible for their formation is still not well understood. To shed light in the chimera producing mechanism, in the present study we introduce an oscillatory model with variable frequency governed by a 3rd order equation. In this model single oscillators are constructed as bistable and depending on the initial conditions their frequency may result in one of the two stable fixed points, $\omega _l$ and $\omega _h $ (two-level synchronization). Numerical simulations demonstrate that these oscillators organize in domains with alternating frequencies, when they are nonlocally coupled in networks. In each domain the oscillators synchronize, sequential domains follow different modes of synchronization and the border elements between two consecutive domains form the asynchronous domains. We investigate the influence of the frequency coupling constant and of the coupling range on the chimera morphology and we show that the chimera multiplicity decreases as the coupling range increases. The frequency spectrum is calculated in the coherent and incoherent domains of this model. In the coherent domains single frequencies ($\omega _l$ or $\omega _h$) are observed, while in the incoherent domains both $\omega _l$ and $\omega _h$ as well as their superpositions appear. This mechanism of creating domains of alternating frequencies offers a reasonable generic scenario for chimera state formation.


I. INTRODUCTION
Systems of nonlocally coupled oscillatory elements often split in domains where the elements oscillate coherently and other domains where the oscillators are incoherent. These counterintuitive, hybrid states are known in the literature as "chimera states" and they occur even if the oscillators are identical and identically linked [1][2][3][4]. Chimera states were first reported in 2002 by Kuramoto and Battogtokh [5,6], while the term "chimera state" was proposed two years later by Abrams and Strogatz [7]. In the original articles the Kuramoto phase oscillator was used and, later on, chimera states were confirmed in diverse oscillatory flows such as the FitzHugh-Nagumo, the Van der Pol, the Stuart-Landau, the Hindmarsh-Rose and Integrate-and-Fire models [8][9][10][11][12][13][14][15][16]. After the pioneer studies in oscillatory flows, chimeras are now frequently observed in coupled chaotic oscillators [17,18] and in coupled discreet maps [19] under a variety of coupling schemes and parameter values.
Despite the extensive numerical evidence of chimera states, their confirmation in many dynamical systems and the experimental observations, the mechanism behind the formation of chimeras remains elusive. Earlier studies have focused on analytical approaches to the Kuramoto model [1,4,5,7] while most recent approaches embrace the idea of the presence of bistable elements in the system [13,17,35,36]. Along these lines, in the present study we propose a toy model to explore further the idea of chimera states produced as a result of the presence of bistable elements. To pursue this idea, an oscillatory toy model of minimal complexity is constructed: it consists of a circularly limiting orbit of constant radius, it has constant mean phase velocity around the orbit, while the dynamical approach to this orbit is of purely exponential type. The uncoupled toy model has an explicit analytical solution [37]. When many toy-oscillators are coupled in a network, the numerically integrated system does not lead to chimera states. We show numerically that it is possible to achieve chimera states by allowing the toy-oscillators to choose between two different frequency levels. This type of chimera states presents, by construction, two levels of synchronization while the dynamics together with the coupling cause the formation of alternating domains with different frequency which are mediated by the incoherent domains. These incoherent domains serve as transition regions and are characterized by a gradient in the values of the mean phase velocities. The proposed mechanism is generic and can be the cause behind many of the systems exhibiting known chimera states.
The organization of the work is as follows: In sec. II we introduce the uncoupled two-frequency nonlinear toyoscillator and we discuss its steady state properties. In sec. III we couple the toy-oscillators in a ring network and we discuss synchronization measures. In sec. IV we examine the emerging chimera morphologies and their multiplicity as a function of the coupling range. In sec. V we vary the frequency coupling which influences the formation of frequency domains and we record the corresponding variations in the chimera states. In sec.VI we analyze the frequency spectra of the oscillators belonging to the coherent and incoherent domains and we show that the incoherent oscillations are the result of the superposition of the two bordering frequencies. In the Concluding section we recapitulate our main results and discuss open problems.

II. THE TWO-FREQUENCY OSCILLATOR
The motivation for introducing the present model is the need of a nonlinear oscillator with explicitly controllable frequency. In this respect we have introduce in Ref. [37] a model oscillator whose trajectory tends exponentially to a closed limiting orbit. We will make use of this model which, throughout its trajectory, keeps a constant, explicit frequency ω, externally controllable as desired. In the following we use the terms "frequency", or "mean phase velocity" or "angular frequency", interchangeably, to refer to the values of ω, although the angular frequency is related to the frequency f by a factor of 2π, ω = 2πf . The model, before frequency modulation, has the following form [37]: This system presents exponential relaxation to a circle with radius R, with relaxation exponent a. The solution of Eq. 1 can be explicitly written: In Eq. 2 the system starts from position (x 0 , y 0 ) : where A determines the initial condition inside a circle of radius R. As time increases the term Ae −at decreases exponentially to 0, giving rise to a purely circular orbit. It is possible to use a similar construction for the case on a limiting orbit with non-equal axes, R 1 and R 2 (an ellipse) [37]. Without loss of generality we will study here the case R 1 = R 2 = R, to keep the computations as simple as possible and with a minimum number of parameters.
To modulate the frequency ω we introduce a third equation treating ω as variable. The resulting toy-oscillator then takes the form: dx dt = −ax + aR In Eq. 3c c is a constant, and ω l , ω c , ω h (standing for low-ω, intermediate-ω and high-ω) represent the three fixed points of Eq. 3c. Depending on the value of c we can have: 1. one attracting fixed point (ω c ), and two repulsive ones (ω l , ω h ) if c > 0 or 2. two attracting fixed points (ω l , ω h ), and one repulsive one (ω c ) if c < 0.
We are interesting in the second case of bistability (2), where depending on the initial frequencies the system can fall in either one of the two attracting fixed points. An example is given in Fig. 1, where starting at t = 0 from ω(0) = 2 we end-up having frequency ω l = 1 (black, solid line), while starting from ω(0) = 4 we end-up having frequency ω h = 5 (red, dashed line). In the next sections, when many oscillators starting from random initial conditions (x(t = 0), y(t = 0), ω(t = 0)) will be coupled in the network, a number of them will tend to the ω l fixed point while the rest will end up on the ω h fixed point depending on their initial condition, see sec. III.

III. TWO-FREQUENCY OSCILLATORS COUPLED IN A RING NETWORK
In this section we couple the two-frequency oscillators in a ring network arrangement. In a network containing N elements, we use the simplest nonlocal coupling scheme where each oscillator is coupled linearly to 2S neighbors: S nearest neighbors on its left and S on its right. Moreover, each x i −variable is only coupled to x j -variables, j = i − S, · · · , i + S, with common coupling constant σ. Similarly, each y i −variable is only coupled to y j -variables, j = i−S, · · · , i+S, with coupling constant σ, while each ω i −variable is only coupled to ω j -variables, j = i−S, · · · , i+S, with coupling constant σ ω . Without loss of generality, we set a common coupling constant, σ, to the x i and y i variables while a different one, σ ω , governs the frequency exchanges. The coupled system of equations takes the form: where all indices are taken mod N . All oscillators start from random initial conditions in the (x, y, ω)-variables.
Assuming that the ω l and ω h are the attracting fixed points, the scenario which is now expected to lead to the formation of chimera states has the following logic: 1. As the system integrates, the frequency ω i of each oscillator will fall on one or the other attracting fixed points ω l or ω h , depending on their initial ω i (t = 0).
2. Because of the coupling σ ω nearby oscillators will be organized in domains having common frequency, either ω l or ω h , while the frequencies will alternate in sequential domains.
3. Within each frequency domain the oscillators will have common frequency and due to their coupling, σ, will also synchronize in phase (x-and y−variables).

Elements in the borders between two frequency domains will be influenced both by their left and right neighbors
(have different frequencies ω l and ω h ) and will oscillate asynchronously, creating the asynchronous domains. 5. Such a chimera state arises, with alternating synchronous domains involving two different frequencies (wells), bordered by the asynchronous domains.
As an example, we present in Fig. 2 the chimera state for the working parameter set with specific parameter values: S = 40, σ = 0.5 and σ ω = 3.0. In panel a) we present the x−variable profile, in panel b) the ω-variable profile and in (c) the spacetime plot of the x−variable and in panel d) the space time plot of the ω profile. A chimera having 4 coherent and 4 incoherent domains is formed. Two of the coherent domains have frequencies ω l and the other two have ω h , while the incoherent domains form the borders between the coherent ones. This figure will be used as an exemplary case in sec. VI for the comparative spectral analysis of the nodes belonging to the coherent and incoherent regions.
The scenario realized above might be at the basis of the chimera states observed in other systems. It is possible that the combined effects of the nonlinear terms and the coupling in these systems, may induce bistability in their frequencies. If this is the case, then the two-frequency scenarios apply and thus the chimera states are created.
Without loss of generality, in the following sections we use the working parameter set: c = −1 (to ensure of the existence of one repulsive fixed point, ω c , surrounded by two attractive ones ω l and ω h ), ω l = 1, ω c = 3, ω h = 5, R = 1, S = 40, a = 1.0, σ = 0.5 and σ ω = 3.0. Using these parameter values, in section IV we vary the coupling range S and in section V the frequency coupling constant to study the chimera variations under changes of these parameters.

IV. VARIATIONS WITH THE COUPLING RANGE
In this section we keep all parameters fixed to the working set and we monitor the chimera properties with variation on the coupling range S. . It is clear that small values of the coupling range S give rise to a large number of coherent (and incoherent) regions, while as we increase the size of S the number of coherent (incoherent) regions decrease. For S > 80 full synchronization is achieved for the working parameter set. This is not unexpected. Provided that the system size remains constant, for large coupling ranges the exchanges between elements cover larger distances, causing communications to larger and larger regions with full synchronization as the ultimate state, see Fig. 3m,n and o.
The decrease in the number of coherent (incoherent) domains with increasing S has also been observed in other systems, such as in FitzHugh-Nagumo [8], Leaky Integrate-and-Fire [38], the Van der Pol [39] and other oscillator networks [2]. An intuitive understanding of this effect is that the coupling range defines the region where oscillators interact and thus can synchronize. Therefor, the larger the coupling range, the larger the synchronized regions and consequently fewer synchronous and asynchronous regions can be accommodated by a system of finite, constant length.
In Fig. 4 we present quantitative results on the ratios of elements that belong to the lower frequency domain ω l (solid, black line), the higher frequency domain ω h (dashed, red line) and the total number of synchronous elements (dashed-dotted, blue line) as a function of neighborhood size 2S. Let us denote by r = 2S/N , the ratio of coupled elements 2S over the total number of elements N in the network. We may notice three regions where the behavior is distinct: a) Small sizes of coupling ranges r < 0.03 (or 2S < 30): Here the number of elements that synchronize close to ω l (ω h ) increases (decreases). b) Intermediate sizes of coupling ranges 0.04 < r < 0.15 (or 40 < 2S < 150): Here the number of elements that synchronize close to ω l decreases and so do the elements that synchronize close to ω h . Around 2S = 160, r = 0.16 all elements that have the highest synchronization frequency disappear and the system is left with a synchronous domain at ω = ω l = 1, while the rest of the elements belong to the asynchronous regime. c) Large coupling ranges r > 0.16 or (2S > 160): Here the number of elements that synchronize at ω l gradually increases and reaches full synchronization beyond r > 0.18, 2S > 180. This scenario is fairly generic and for large values of the coupling constant one of the two frequency domains dominates at the expense of the other.

V. VARIATIONS WITH THE COUPLING CONSTANTS
It is interesting to study the variations of this models with the coupling constants. As in most cases of synchronization in the form of chimera states the coupling constant σ which governs the amplitude synchronization does not affect the chimera multiplicity but only the size of coherent and incoherent domains. We shortly discuss this case using an example in the Appendix A. As the exemplary case shows, the size of the incoherent regions decreases as σ increases, leading to full synchronization for large values of the coupling strength σ .
After the brief discussion on the amplitude coupling strength, we now turn to the more interesting case of the variations in the frequency coupling range σ ω . Figure 5 provides a first account on the formation of the two frequency chimera state as we turn on the coupling on the frequency variables. In the absence or for small values of the σ ω all oscillators fall fast in their attracting fixed points (ω l or ω h ) and they perform oscillations with these frequencies.
Neighboring oscillators are not affected and the system remains asynchronous in time. Such an example is presented in the top row of Fig. 5, with σ ω = 1.5. From the spacetime plot of the mean-phase velocity, Fig. 5b, it is evident that all oscillators acquire a constant in time ω i . As frequency coupling increases in the second row of Fig. 5 to σ ω = 1.5, a part of the system develops synchronization in the highest frequency, ω h = 5, while the rest of the oscillators to the left and right of the synchronous regions have mixed frequencies, see Fig. 5e. This is because the asynchronous regions keep their local frequencies which have been shaped as the oscillators were attracted by fixed points ω l or ω h . The corresponding x-profile, Fig. 5d, demonstrates coherent motion in the synchronous region and incoherent outside of it, while the x-spacetime plot, Fig. 5f, indicates that even within the incoherent domain small regions of random sizes are formed with almost coherent temporal behavior. This is not discernible in the x-profile but is visible in the spacetime plot.
By increasing further the frequency coupling to σ ω = 2, see Figs. 5h,i,g, the previously asynchronous domain now synchronizes to the higher frequency, ω l , thus leading to a chimera state consisting of two coherent domains, one in the high frequency, ω h = 5 and one in the low frequency ω l = 1. Two incoherent domains are develop which serve as the borders between the two coherent ones. The two incoherent domains also serve for continuity purposes to bridge the gap between the domains where ω h and ω l dominate. As the frequency coupling increases further to σ ω = 2.5, see Figs. 5j,k,l, the high frequency domain splits giving rise to four, coherent domains bordered by four incoherent ones. By further increasing σ ω = 5.5 the exchange of ω variables become dominant and as time increases the neighboring domains compete and the lower frequency domains expand in expense of the higher ones, see Figs. 5j,k,l. All runs in Fig. 5 start from the same random initial state.

VI. FREQUENCY SPECTRA OF COHERENT AND INCOHERENT ELEMENTS
To investigate the transition from the coherent to incoherent regions and to check which precise frequencies are present in each region we investigate the Fourier spectra of different oscillators belonging to the core of the coherent and others occupying borderline regions.
We analyze here the results in Fig. 2 and plot the Fourier spectra of nodes i = 440, 650 and 740. The first one, i = 440, is centrally located in the coherent domain which oscillates with ω = 1, the second, i = 650, belongs to the coherent domain with ω = 5, while the third one, i = 740, occupies a position in between the two, in the incoherent region which serves as a transition region between the two domains. The results are shown in Fig. 6.
The left panel of Fig. 6 depicts the Fourier amplitude of node i = 440. Only one peak is clearly seen at frequency values f = ω/2π = 0.167402, or ω ∼ 1, as expected for the oscillators which have fallen in the basin of attraction of the lowest frequency, ω l = 1. The middle panel, Fig. 6b, depicts the Fourier amplitude of node i = 650. Also here, only one peak is developed at frequency values f = ω/2π = 0.789951, or ω ∼ 5, as expected for the oscillators which have fallen in the basin of attraction of the lowest frequency, ω h = 5. The right panel, Fig. 6c, depicts the Fourier amplitude of node i = 740. Here four peaks appear at frequency values i) f = ω/2π = 0.167402 (ω ∼ 1), ii) f = ω/2π = 0.79 (ω ∼ 5), iii) f = ω/2π = 0.642 (ω ∼ 4.03) and iv) f = ω/2π = 0.496 (ω ∼ 3.1). In these regions the oscillators are developing mixed behavior and present mixed oscillatory characteristics, drawing both from the low frequency dynamics (Case i) and the high frequency dynamics (Case ii). In addition two more peaks with high amplitude are developed which can be considered as linear combinations of the ω h and ω l values. E.g., the peak with the highest amplitude in the right panel of Fig. 6 which corresponds to ω = 4.03 is approximately equal to (ω l +ω h )/2. Within the incoherent region, the closer the element is to the high (low) frequency domain, the higher the amplitude of the corresponding peak is (images not shown).
These findings corroborate the intuitive argument that the two different types of coherent domains are formed due to the bistability of the oscillator frequency caused by the addition of the linear coupling terms to the nonlinear dynamics. In this view, the incoherent regions serve for the purpose of continuity when passing from the lower to the higher frequency domains (and the opposite). That is the reason why they have a gradient in frequencies, giving rise to the asynchronous incoherent domains in the x-variable profile (and y−variable profile, not shown).
At this point one would argue that a number of chimera states do not present two-level synchronization, but they demonstrate an arc-shaped mean phase velocity profile. Given the present results, it is not possible to conclude whether the arc-shape ω-profile of the chimera states in the Kuramoto [5] or the FitzHugh Nagumo [8] models emerges as an incompletely formed higher (or lower) frequency domain, or if some other phenomenon is responsible for this profile.
As an additional indicator of the different frequencies dominating in the coherent regions, we plot in Fig. 7 the Fourier amplitude of the peak which corresponds to the low mean phase velocity, ω l ∼ 1 (black line) and to the high one, ω h ∼ 5 (red line). Comparing Figs. 7 and 2b we note that the Fourier amplitude of the peak at frequencies ω ∼ 1 dominates in the coherent regions of low ω in Fig. 2b and vanishes gradually in domains of high frequency (see Fig. 7, black line). Similarly, the red line which depicts the Fourier amplitude extracted from the the peak at frequencies ω ∼ 5 reaches maximum values in the coherent domains of high ω in Fig. 2b and Fig. 6c, whose spectra are the linear combinations of the ω h and ω l values. These additional spectral lines, clearly appearing in 6c, are omitted in Fig. 7 for clarity reasons.
The results in this section indicate that the calculation of the Fourier spectra is a laborious but reliable method to identify qualitatively and qualitatively the different coherent and incoherent domains and to assert the presence (or absence) of a chimera state.

VII. CONCLUSIONS
In the present study we first introduced a model oscillator whose trajectory tends exponentially to a closed limiting orbit with well defined frequency and radius. As numerical results indicate, no chimera states arise when such oscillators are coupled, but the system tends to either full synchronization or to full disorder. When a third equation is added which allows the oscillators to choose between two different frequencies, a higher and a lower ones, then domains are formed where the high and low frequencies dominate. The domains of different frequencies are mediated by incoherent domains where the oscillator frequencies gradually increase or decrease to bridge the frequency gap between adjacent regions and to maintain continuity in the system. This scenario can serve as a general mechanism for formation of chimera states. Even in cases where this mechanism is not explicit, the introduction of coupling together with the nonlinearity in the dynamics may induce bistability in the oscillator frequency creating, in this way, indirectly, two-level synchronization and corresponding chimera states.
The present model can be easily extended to form multi-leveled chimeras, by allowing the oscillators to occupy many different frequency levels.
Finally, in Ref. [37] except for the case of the exponential relaxation to the limiting orbit, power law relaxation has been introduced. Power laws take much longer (infinite time) to reach the final oscillatory trajectory. It would be interesting to investigate whether this power law relaxation dynamics leads also to the formation of chimera states and, further-on, if frequency multistability can lead to multi-leveled chimera states.
As in most cases of local synchronization in the form of chimera states, the coupling constant σ which governs the amplitude synchronization does not affect the chimera multiplicity but only the size of coherent and incoherent domains. As an exemplary case we present here simulation results using the following parameter set: c = −1, ω l = 1, ω c = 3, ω h = 5, R = 1, a = 1.0, S = 40 and σ ω = 3.0, with variable size of σ = 0.2, 0.6, 1.0 and 1.5. In Fig. A1 the spacetime plots of the variable x are presented in the four cases. For small σ values the incoherent regions extend to a large number of oscillator, while the size of the incoherent domains decreases for larger σ-values.