Cloning of Chimera States in a Large Short-term Coupled Multiplex Network of Relaxation Oscillators

A new phenomenon of the chimera states cloning in a large two-layer multiplex network with short-term couplings has been discovered and studied. For certain values of strength and time of multiplex interaction, in the initially disordered layer, a state of chimera is formed with the same characteristics (the same average frequency and amplitude distributions in coherent and incoherent parts, as well as an identical phase distribution in coherent part), as in the chimera which was set in the other layer. The mechanism of the chimera states cloning is examined. It is shown that the cloning is not related with synchronization, but arises from the competition of oscillations in pairs of oscillators from different layers.


INTRODUCTION
Study of the formation of chimera states, i.e., peculiar types of hybrid states consisting of oscillators with coherent and incoherent behavior is one of the hot problems of the modern non-linear dynamics. To date, the chimera states have been discovered not only in a variety of theoretical papers [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18], but also in experimental systems of various natures, for example, mechanical [19][20][21][22], optical [23,24], chemical [25][26][27][28][29], and radiotechnical ones [30][31][32][33]. Similar states have also been registered in the neural activity of animal brain networks [34,35]. At present, great attention is paid to studying of interaction of chimera states. The effects of generalized synchronization of chimera states [36], synchronization of chimera states in ensembles with asymmetrical connections [37], synchronization of chimera states in multiplex networks with delays [38], synchronization of chimera states in a two-layer multiplex network with adaptive connections in each layer [39], synchronization of chimera states in modular networks [40,41], interaction of chimera states with fully coherent or fully incoherent states [42], etc. were explored. Note that in all these works the interaction of chimera states led to the formation of new chimera states (in some cases with synchronous coherent parts) that are different from the pre-existing chimera state. Recently we presented an example [43] of a two-layer multiplex network with seven oscillators in each layer where due to short-term interaction, one more chimera is emerged which is identical to the initial one (excluding the phase distribution of the incoherent part). We called this effect the chimera states cloning. In this article, we generalize the results of Dmitrichev et al. [43] for the case of a multiplex network with an arbitrary dimension of the layer and give a theoretical justification for the cloning effect based on the study of the fast-slow dynamics of the model using the methods of Geometric Singular Perturbation Theory (GPST) [44,45].

MODEL OF MULTIPLEX NETWORK
We consider a two-layer multiplex network with the topology illustrated in Figure 1A. Each layer of the network is a ring of locally and linearly coupled relaxational oscillators with phase portrait shown in Figure 1B. The dynamics of the network is described by the following system:  (B) Qulitative phase portrait of a single oscillator. There exist two stable limit cycles (bold black lines) separated by an unstable limit cycle (bold blue line) and an unstable equilibrium state (bold blue dot) located at the origin.
If the oscillators do not interact with each other, i.e., d r = 0; d m (t) ≡ 0 , then the dynamics of each oscillator is described by second-order equation. Two stable limit cycles with "low" and "high" amplitudes exist on the (u, v) phase plane (see Figure 1B). The regions of attraction of these cycles are separated by an unstable limit cycle. An unstable equilibrium state is located at the coordinate origin (u = 0, v = 0). Thus, an oscillator can be either in the regime of low-amplitude oscillations with a dimensionless frequency of 0.0039 or, in the case of initial conditions "outside" the unstable cycle, in the regime of highamplitude oscillations with a dimensionless frequency of 0.0021.
In our previous paper [32] we showed that various chimera states exist in a separate layer of system (1) at the chosen parameters and obtained chimera states at the experimental system consisting of seven bistable self-exciting oscillators with linear couplings. The example of chimera state is presented in Figure 2A by black crosses (distribution of instant phases ϕ, average amplitudes < A >, and average frequencies < ω >). Coherent part of the chimera state is formed by first 99 (j = 1 − 99) oscillators with high-amplitude oscillations, and incoherent part is formed by oscillators with j = 101 − 199 that demonstrate alternately the low-and high-amplitude oscillations. The oscillation frequencies and phases of oscillators were calculated as follows. For definiteness, we consider the jth oscillator in the ith layer. Let {t n j } be a sequence of time instants at which the output voltage of the oscillator increases and intersects the straight line u i j = 0; i.e., Then, the oscillation phase of the jth oscillator at the time t is given by the expression and ω n j = 1 t n j −t n−1 j is the instantaneous oscillation frequency.
This definition of the phase is meaningful only if w n j are constant or quite close to each other. In the former case, the oscillator undergoes regular oscillations and the phase always varies at the same rate. In the latter case, oscillations can be, in particular, irregular and the phase is a piecewise linear function of the time. Furthermore, the phase introduced in Equation (2) describes the dynamics of only an individual oscillator decoupled from the remaining system. For this reason, to describe the dynamics of the system as a whole, it is more convenient to use the parameter ϕ k j = φ j − φ k describing the oscillation phase of the jth oscillator with respect to oscillations of the reference kth oscillator. If ϕ k j is independent of time, this means the phase matching of oscillations of the kth and jth oscillators. In the general case of interaction between oscillators, instantaneous frequencies and amplitudes are not constant. For this reason, here and below, frequencies and amplitudes are calculated with averaging over a quite long time series  and Notice that in addition to the coherent and incoherent parts, the chimera state also contains two isolated oscillators at j = 100 and j = 200. Such solitary states exist due to the bistability of the network oscillators. The explanation of this phenomenon was given in Nekorkin et al. [46]. The bistability leads to formation of amplitude clusters with highand low-amplitude oscillations and so to an amplitude gap (see the amplitude distribution in Figure 2A). The amplitude gap, in turn, causes a frequency gap. Since the magnitudes of gaps are quite large, diffusive coupling between the oscillators leads to emergence of solitary oscillators whose dynamics "smooths out" the dynamics of clusters (amplitude and frequency) having significantly different characteristics.

CLONING OF CHIMERA STATES
Let the chimera state exist in the first layer at the initial instant of time when there is no interaction between the layers (black crosses in Figure 2A). And the initial conditions in the second layer correspond to the low-amplitude oscillations with random initial phases (red circles in Figure 2A). Now if we switch on the interaction between the layers, and then after some time switch the interaction off, then we can obtain a clone of the initial chimera state in the second layer. For certain values of strength and time of multiplex interaction, in the second layer, a chimera state is formed with the same average frequency and amplitude distributions in coherent and incoherent parts, as well as an identical phase distribution in coherent part, as in the chimera which was set in first layer. The example of such cloning for interaction strength d = 0.06 and interaction time T c = 1, 000 is shown in Figure 2. Notice that, by definition, instantaneous phases of the incoherent parts in the chimera state should be random, and integral characteristics, e.g., average frequencies and amplitudes, are of fundamental importance for the incoherent part. For this reason, we believe that the coincidence of the phases of the incoherent parts is not necessary for the cloning of chimera states. Note that the cloning of chimera states occurs with certain initial conditions The role of initial conditions in more detail is discussed in section 4.
Next we show that the cloning effect is structurally stable. To do this, we introduce a reference layer with the index "0." Let us set in the reference layer a chimera state with the same characteristics as the original chimera that we set in the first layer. Then after the interaction of the first and the second layers we compare states formed in those layers with one in the reference layer using the following characteristics: where < ω i,j > and < A i,j > are averaged frequencies and amplitudes of the oscillators with the number j of the ith layer; accordingly, < ω 0,j > and < A 0,j > are those in the reference layer.
The results of such calculations for T c = 1, 000 are presented in Figure 3. They indicate that there is an interval of multiplex coupling strength 0.37 ≤ d ≤ 0.96 where the states in all layers have the same averaged characteristics as the reference chimera state. Thus, the cloning occurs in the large enough interval of strength and so the effect is structurally stable.

THE CLONING MECHANISM
We showed above that cloning of chimera states takes place when strength of coupling between the elements of the same layer is much smaller than that between the elements of different layers (d r << d, see Figure 3). So in the first approximation we can assume that the key role in cloning is played by the dynamics of (multiplex) pairs of elements taken from different layers. Next we consider the dynamics of a pair in more detail. It is described by the following system of equations: where Notice also that to realize the cloning, initial conditions in non-interacting layers must be formed in a special way. In particular, a chimera state is set in the first layer with coherent part formed by the oscillators, demonstrating high-amplitude oscillations, and incoherent part formed by the oscillators demonstrating alternately low-and high-amplitude oscillations.
In the second layer all oscillators are set in the regime of lowamplitude oscillations whose phases are randomly distributed. Moreover, after interaction, the elements of the second layer should switch to the regimes the corresponding elements of the first layer were in initial moment. Thus, we need to consider the evolution of a pair only for two types of initial conditions: (I.C.) 1 An oscillator of the first layer is in the regime of highamplitude oscillations, while that of the second layer is in the regime of low-amplitude oscillations; (I.C.) 2 The oscillators of both layers are in the regime of low-amplitude oscillations.

Dynamics of a Pair of Constantly Coupled Oscillators
First, we study dynamics of a pair for the case when interaction between its oscillators is not limited by time. Then we obtain the following system of equations: Since 0 < ε ≪ 1, the system (4) belongs to the class of fast-slow systems. Such systems are characterized by the presence of two timescales (or speeds), namely, fast and slow ones. In the result, the trajectories of the systems have epochs of a slow and a fast movements. In our system u 1 and u 2 are fast variables, while v 1 and v 2 are slow variables. Next to study the dynamics of the system we use GPST theory. According to the GPST, the partition of phase space R 4 of system (4) into trajectories can be established by studying two subsystems. As ε → 0, the trajectories of system (4) converge during fast epochs to the trajectories of the fast subsystem (or layer equations) where t = ετ . During slow epochs the trajectories of (4) converge to the trajectories of the reduced system (or the slow flow) The goal of GPST is to use the fast and slow subsystems (5) and (6) to understand the dynamics of the full system (4) for 0 < ε ≪ 1.

Dynamics of the Fast Subsystem
From system (5) one can see that v 1 = const and v 2 = const, so they play the role of additional parameters (denote them by v 0 1 and v 0 2 correspondingly). Thus, system (5) can be rewritten in the following gradient form where the trajectories of system (7) [and so system (5)], except for equilibrium states, relax to one of the stable equilibrium states. Moreover, since the system is gradient their trajectories relax to the equilibrium states by the fastest way. The number and type of equilibrium states depend on the parameters and may change due to saddle-node bifurcations. For example, for v 0 1 = v 0 2 = d = 0 there are 49 equilibrium states, among which there are 16 stable and 9 unstable nodes and 24 saddles. The qualitative phase plane of system (7) in this case is shown in Figure 4.

Dynamics of the Reduced Subsystem
The first two algebraic equations in system (6) define a critical manifold S{(u 1 , Let us rewrite (6) in terms of the fast variables u 1 and u 2 to obtain the slow flow on S. For this, we differentiate algebraic Equation (6) with respect to t and combine the result with the equations forv 1 andv 2 : System (8) is a system of linear inhomogeneous algebraic equations for derivatives du 1 dt and du 2 dt . Determinant of its coefficient matrix If ∆ = 0, then system (8) has the only solution Note that in the four-dimensional phase space R 4 of system (4), algebraic equation ∆ = 0 and two algebraic equations of system (6) define set S ∆ = ∪ j S ∆ j . Any its element, S ∆ j , defines a curve, where the fast and slow trajectories are stitched. According to GPST, each equilibrium state of fast subsystem (5) [and so (7)] in the phase space R 4 of system (4) corresponds to a submanifold, whose stability with respect to the trajectories of the fast subsystem coincides with the stability of the corresponding equilibrium state. Since the coordinates of the equilibrium states of system (5) depend on two parameters v 0 1 and v 0 2 , any such equilibrium state corresponds to a two-dimensional submanifold in the phase space of (4). The boundary of the submanifold is given by one of the curves S ∆ j . Hence the curves S ∆ j decompose the critical manifold S into some number of submanifolds of different stability where S a i denotes one of the stable submanifolds, S r k are the unstable submanifolds, and S sad n are the saddle submanifolds of the slow flow. Note that ∆ > 0 on submanifolds S a i and S r k , ∆ < 0 on submanifolds S sad n and ∆ = 0 on submanifolds S ∆ l . The number of these submanifolds depends on the parameter d. For example, Figure 5A for d = 0.06 on the phase plane (u 1 , u 2 ) depicts submanifolds S a i , i = 1, 16, corresponding to stable nodes O a j of the fast subsystem and curves S ∆ l . Since the system (4) on submanifolds S a i has no equilibrium states and limit cycles, any trajectory starting on the submanifolds eventually leave them. For example, Figure 5B shows the behavior of trajectories on S a 8 starting on entering part of its boundary S ∆ 8 (blue line) and going till the exit part of the boundary (red line). Note that, the time the trajectories stay on S a i varies depending on the initial conditions. The dependence for trajectories on S a 8 is shown in Figure 5C. Note that the dependence is a monotonically increasing function asymptotically tending to the value T ≈ 220.

Dynamics of System (4) for Initial Conditions (I.C.) 1
Let us study the dynamics of system (4) for initial conditions (I.C.) 1 . Note, the dynamics of system (4) is formed by the alternating dynamics of fast and slow epochs, which results in a "stitched" trajectory.

Slow epoch of motion
The initial conditions (I.C.) 1 corresponds to one of the stable submanifolds of slow motions S a 5 , S a 8 , S a 9 , S a 12 (see Figure 5A). First let the initial conditions belong to S a 8 . The motions on S a 8 are defined by system (9). Since ∆ > 0 on S a 8 , we can de-singularize the slow flow near S ∆ 8 by rescaling time with the factor ∆. This gives the following system: where dt = ∆dt n . System (6), and hence system (10), have the only equilibrium point at the origin. Therefore, the trajectories starting on S a 8 leave the submanifold at some points of its boundary S ∆ 8 . Let O sn (u 1 = u 0 1 , u 2 = u 0 2 ) ∈ S ∆ 8 be one of such exit points (see Figure 5B). Since the fast and slow trajectories of system (4) are glued together on S ∆ 8 , the point O sn is also the equilibrium state of the fast system (5). From this condition we find By using Equation (11) fast system (5) can be rewritten in the form The eigenvalues of the Jacobian matrix of (12) in the point O sn are given by Because of Equation (13)

Fast epoch of motion
Consider the level curves of the function G(u 1 , u 2 ) = C = const, satisfying the condition C ≤ C sn , where C sn = G(u 0 1 , u 0 2 ). Obviously, the curves of the level {G(u 1 , u 2 ) = C sn } pass through the saddle-node O sn . In Figure 6A, the curve of this level has a maximum value, i.e., the curves of the higher levels are not indicated (white color), and the curves of the level corresponding to the values of the lower levels (C < C sn ) are marked with different colors. Each color corresponds to the same value of C. Solid black lines in Figure 6 show the critical lines At the points of these lines the following conditions are satisfied: Consider the asymptotic behavior of the separatrix of the saddlenode O sn . Taking into account Equation (8), the location of the level curves of the function {G(u 1 , u 2 ) = C sn } and lines (14), we establish that for the parameter value d = 0.0075 ( Figure 6A) the separatrix W u (O sn ) asymptotically tends to the equilibrium state O a 12 . Without changing the coordinates of the point O sn , we increase the value of the parameter d = 0.0115 (Figure 6B). For this value of the parameter d the lines L h min , L v max , L v min merge at one point, a saddle-node bifurcation of equilibrium states occurs, and a new equilibrium state O * appear (see Figure 6B). With the further increase in the parameter, the equilibrium state O * disappears, and following the arrangement of lines (14) and the level curves of the function G, we find that in this case the separatrix W u (O sn ) tends to the node O a 9 (Figure 6C). A further increase in the parameter d leads to the merging of the lines L h min , L h max , L v min (d = 0.020909 = d * ). This corresponds to the merging of the node O a 9 and the saddle O sad (Figure 6C) and emerging of the saddle-node equilibrium state. For d > d * , this equilibrium state disappears, and the separatrix W u (O sn ) tends to the equilibrium state O a 13 ( Figure 6D). It is clear that d * depends on the coordinates of the point O sn .
To describe such possible transitions, we introduced the distance R on the plane (u 1 , u 2 ) of system (9) from the origin to

S a
i . Note that the largest value of R corresponds to the points of the submanifolds S a 1 , S a 13 , S a 4 , and S a 16 . Figure 7A depicts the results of analyzing the behavior of the separatrix W u (O sn ) for different values of u 0 1 belonging to the line of escape (u 0 1 , u 0 2 ) ∈ S ∆ 8 with the submanifold S a 8 (see Figure 5A). For the values of u 0 1 , d from the region marked by red color, the separatrix W u (O sn ) of any of the saddle-nodes in the fast system (5) asymptotically tends to the stable node O a 13 . Since in the phase space R 4 the equilibrium state corresponds to a stable manifold S a 13 , high-amplitude oscillations are established in both elements in the system (4). Note that the trajectories of the submanifold S a 9 of the slow system (10) have a similar behavior. Since the function ∂f ∂u is even, system (10) does not change when converting (u 1 , u 2 ) → (−u 1 , −u 2 ). Thus, transitions from S a 9 to S a 16 exist in R 4 , and high-amplitude oscillations are established in system (4). The initial conditions found by us do not exhaust the entire set of initial conditions under which the oscillation amplitude changes from low to high in the second oscillator. Let us study the dynamics of system (4) for initial conditions (I.C.) 2 . These conditions correspond to one of the stable submanifolds of slow motions S a 6 , S a 7 , S a 10 , and S a 11 . Similar to the case of (I.C.) 1 we have analyzed the behavior of the trajectories leaving the submanifolds. Figure 7B depicts the transitions of trajectories starting from the points at the boundary S ∆ 7 of submanifold S a 7 obtained for different strength of coupling d. One can see that there are no transitions to submanifolds corresponding to high-amplitude oscillations. We have established that such behavior is also typical for other submanifolds, namely S a 6 , S a 10 , and S a 11 .

Dynamics of a Pair of Short-term Coupled Oscillators
So far, we have considered the dynamics of system (4) without any restrictions on the interaction time. However, in the initial model, the layers interact only during the time T c . We numerically investigated the dynamics of system (3). For the initial conditions such as (I.C.) 1 , Figure 8A shows the behavior of 1, 000 pair of oscillators, or in other words, 1, 000 different initial conditions (I.C.) 1 type, interacting during T c = 300 and d = 0.06. For all initial conditions, after some transition process, highamplitude oscillations are established in the interacting pairs. Figure 8B illustrates competition of oscillations in the pairs in the case of initial conditions such as (I.C.) 2 type. Here, interaction of the pairs does not lead to high-amplitude oscillations, and the regime of low-amplitude oscillations persists. The occurrence of high-amplitude oscillations for the initial conditions (I.C.) 1 depends on the values of the parameters T c and d. We examine this dependence numerically and the results are presented in Figure 9. The color gradation on the plane (d, T c ) shows the dependence of the probability of establishing high-amplitude oscillations. The area highlighted in black corresponds to the establishment of high-amplitude oscillations from any initial conditions. The area highlighted in shades of gray corresponds to the establishment of highamplitude oscillations from only some initial conditions. And finally the area marked in white corresponds to those values of the parameters for which high-amplitude oscillations are not established at all. Note that there are threshold values for both parameters d and T c . The existence of threshold value for d has already been discussed above. The presence of threshold value for T c is associated with the motion time T over stable submanifolds (see Figure 6C). Note that the value of the critical value T c is determined by the dynamics of both oscillators and it is not related to the periods of high and low oscillations of the isolated oscillator.

Dynamics of the Multiplex Network
Thus, it has been established that in the case of initial conditions (I.C.) 1 , there are values of the parameters d and T c corresponding to the emergence of high-amplitude oscillations in the pairs of interacting oscillators belonging to the different layers. On the other hand, it has been shown that in the case of initial conditions (I.C.) 2 , when the oscillators belonging to the different layers do not change their initial regimes after the interaction and keep demonstrating low-amplitude oscillations. Now let us consider the dynamics of multiplex network (1) based on the findings of the previous subsections. It can be divided into two main stages.
(a) In the time interval 0 < t < T c , oscillators of different layers interact with each other through inter-layer couplings with strengths, d c , greatly exceeding those of the diffusive intra-layer couplings, d r . Therefore, in this stage the main contribution to the dynamics of the system due to the dynamics of the interacting pairs. We have established that as a result of this dynamics, the pairs of oscillators with high-amplitude oscillations are formed in (1) from the initial conditions (I.C.) 1 . On the other hand, the pairs of oscillators with low-amplitude oscillations are formed in (1) from the initial conditions (I.C.) 2 . This means that the average amplitude distribution in the first layer does not change, while that of the second layer becomes the same as in the first one.
(b) For t > T c , there are no inter-layer couplings, and the oscillators interact only through diffusive intra-layer ones. Under the influence of these couplings, the neighboring oscillators with similar amplitudes become phase-locked with each other at some average frequency and form the coherent part of the chimera state. The neighboring oscillators with different amplitudes do not become phase-locked with other oscillators and form the incoherent part of the chimera state with distinguished bellshaped distributions of average frequencies and amplitudes. Thus, the same chimera state is formed in the second layer as in the first one. Note that a finite interaction time is required to stop the competition of oscillations of pairs of oscillators. Otherwise, new complex states arise in the layers and they differ from the initial chimera.

CONCLUSIONS
In a large two-layer multiplex network with short-term couplings, a new phenomenon of the chimera states cloning, has been discovered and studied. Each layer of the system has a ring topology and consists of relaxation oscillators having two stable limit cycles on their phase planes. The oscillators inside the layers interact through diffusive couplings, while those of different layers interact by means of multiplex couplings. When the chimera state existing in one of the layers interacts for a while with oscillations of the other layer having a random distribution of phases, the same chimera state appears in the latter layer. Note that the time of occurrence of the chimera state in the second layer is less than the minimal partial oscillation period. We have found that the phenomenon is not related with synchronization of oscillations existing in the layers, but instead is determined by the competition of high-and low-amplitude oscillations. Using GPST, we showed that competition of oscillations in each (multiplex) pair of oscillators in the multiplex network is controlled by switching four-dimensional slow-fast dynamics. We have analytically established the initial conditions leading to the trajectories in phase space which start from stable "competitive" submanifolds of slow motions and then transit to stable "winner" submanifolds. The "competitive" submanifold corresponds to the case where oscillations in different layers have different (low and high) oscillation amplitudes. The "winner" submanifold corresponds to the case where oscillations in different layers have high amplitudes. Transitions between stable submanifolds occur along the trajectories of a two-dimensional fast subsystem. The given initial conditions belong to the basin of attraction of both the initial chimera state and the clone. We found that strength, as well as time of multiplex interaction, play a crucial role in the existence of the cloning effect of chimera states. A chimera clone is formed with 100% probability if the strength and time of multiplex interaction exceed certain threshold values.
Below these threshold values a chimera clone occurs with a certain probability. Note that the effect of chimera state cloning does not depend on the choice of boundary conditions, since the dynamics of pairs of oscillators plays a crucial role in its existence. We hope also that the cloning effect is not specific to considered model and exists in other models, since the conditions necessary for it to take place are fairly general.