Symmetry-breaking mechanism for the formation of cluster chimera patterns

The emergence of order in collective dynamics is a fascinating phenomenon that characterizes many natural systems consisting of coupled entities. Synchronization is such an example where individuals, usually represented by either linear or nonlinear oscillators, can spontaneously act coherently with each other when the interactions' configuration fulfills certain conditions. However, synchronization is not always perfect, and the coexistence of coherent and incoherent oscillators, broadly known in the literature as chimera states, is also possible. Although several attempts have been made to explain how chimera states are created, their emergence, stability, and robustness remain a long-debated question. We propose an approach that aims to establish a robust mechanism through which chimeras originate. We first introduce a stability-breaking method where clusters of synchronized oscillators can emerge. Similarly, one or more clusters of oscillators may remain incoherent within yielding a particular class of patterns that we here name cluster chimera states.

Synchronization of the dynamics of the individual entities that constitute complex systems has been observed and systematically studied for over four decades 1-3 . An exotic behavior that has triggered scientists' curiosity for a long time is the emergence of chimera states, longlasting dynamical patterns of coexistence of synchronous and asynchronous clusters of nodes in networked systems [4][5][6][7] . Although many features are now understood, further questions remain unanswered 8,9 . In particular, the stability of such states for finite networks of oscillators and their high sensitivity to the initial conditions have not yet been tackled exhaustively. Several recipes have been suggested recently that aim to reproduce chimera states by decomposing the network into stable and unstable clusters [10][11][12][13][14] . Here we show that the formation of clustered synchronization 15,16 can follow a global symmetry-breaking mechanism, which is also responsible for selecting the nature of the final patterns [17][18][19] . In this way, we demonstrate that it is possible to reconcile the pattern formation process with cluster synchronization, paving the way to a robust mechanism for explaining the formation of chimera states in networks of coupled dynamical systems.

I. INTRODUCTION
The full understanding of the dynamics of complex systems is probably one of the biggest challenge of physics this century. The emergence of collective behaviour and the coexistence of order and disorder are fundamental key players in this process. A very well-known phenomenon of collective behavior is that of synchronization, where a set of coupled oscillators spontaneously synchronize their phases. This phenomenon has been widely studied in recent years and has been applied to explain the coordination of the beats of the heart myocytes and the firing of brain neurons 1,20,21 , the simultaneous flashing of male fireflies during the mating season 22 and the synchronization of the alternating current in power grids 23 . Theoretical studies have identified particular types of synchronization, among which it is worth mentioning the remote synchronization of two or more oscillators through asynchronous ones 24 , the synchronization of chaotic oscillators and the phenomenon of cluster synchronization when the set of oscillators is divided into subsets of synchronized ones 16 or the emergence of exotic states where identically coupled oscillators synchronise in clusters 25 .
The latter class also includes an essential and fascinating phenomenon of coexistence of order and disorder known as chimera states. Chimeras, named from Greek mythology, are particular states where clusters of coherent and incoherent oscillators coexist. Since their discovery, chimera states have aroused scientists' curiosity about the mechanisms that give birth to them, about the effective lifespan of these states ? , and their sensitivity to the initial conditions that make their experimental verification difficult to achieve 4,9 .
Outstanding theoretical results have been obtained in the study of chimeras, and promising methods have been proposed to reveal the mystery surrounding these exotic states, such as whether they will exist, their stability, and specially, their origin. For example, it has been found that chimeras are stable in infinite-size networks and that they can be transient in finite ones 7 . So far, however, there is not an exhaustive answer to the stability of chimera states 5,6 , or a robust mechanism that explains how chimeras are generated 13 . An important advance in the area came from the group-theoretical char-acterization of network symmetries that corresponds to the identification of partitions of nodes that can be stable or unstable [10][11][12][13] . Using this approach the authors are able to design the emergence of chimera states. The most recent approach to understanding the coexistence of both coherent and incoherent oscillators has been that of group-theoretical characterization of network symmetries that, corresponds to the identification of partitions of set of nodes that can be stable or unstable [10][11][12] . Other recent alternatives have been based on symmetry-breaking principles where chimeras emerge from an earlier uniformly synchronized regime 13,14 . Aside from theoretical results, experimental evidence 26 proves the existence of cluster synchronization and chimera states for the case of coupled Belousov-Zhabotinsky oscillators 17 .
In this paper, we aim to analyze the chimera states, proposing an alternative path to the aforementioned cluster partition. More explicitly, we will use a pattern formation approach based on a global symmetry-breaking mechanism 17-19 that yields cluster synchronization patterns as a result. Based on this principle, we will successively generate coexisting clusters of coherent and incoherent oscillations that are robust to initial conditions. We name these states cluster chimera patterns. The great advantage of this method is that it has been extensively studied since the early 1970's, and many groundbreaking rigorous analytical results have been obtained regarding the prediction of the shape and stability of the final nonlinear patterns. In fact, the pattern formation process is robust to the choice of the initial conditions, at odds with classical chimeras that are highly sensitive to the initial setting of variables 4,5,9 . Furthermore, it has been proven that patterns obtained following an instability principle are not transient, but a stability region of parameters (also known as "stability balloon") where they are asymptotically stable exists 18,19 .
A major difference in our approach is the model setting through which oscillatory behavior emerges. Traditionally, the synchronization phenomenon has been studied in systems of coupled oscillators. The main reason is that global coherence is regarded as the emergent behavior in these systems. Thus it is expected that the individual nodes have to oscillate independently before fully coherent or chimera states eventually emerge. At variance with this framework, here we show that cluster synchronization and chimera patterns can also arise when the individual nodes do not manifest any oscillatory behavior. As we will next show in detail, our method is fully based on a global instability mechanism where the oscillations are also generated due to either inter-node coupling 29 or multispecies interactions 32 , further facilitating the formation of synchronized or chimera clusters. In fact, a key aspect of our analysis is the construction of networks whose structural properties yield the segregation of the entries of the eigenvectors of the coupling operator into clusters. We demonstrate that sufficiently modular networks, in particular, possess this property.
The aforementioned characterizing features of the pattern formation process are crucial because they will pave the way and facilitate obtaining chimera states. Thus, instead of tackling the problem frontally to understand how to decompose the chimera in groups of locally stable and unstable nodes, we propose a recipe for rigorously proving the emergence of simultaneously coherent and incoherent oscillatory nodes.

II. SYMMETRY-BREAKING MECHANISM
We start by considering the general setting of a network of coupled dynamical systems through the equations where each of the individuals, identified as nodes of the network, are represented by -dimensional vector state variables x 2 . The local dynamics of the -th node is given by the nonlinear function F (x ), σ is the vector of the coupling strength, are the weighted entries of the adjacency matrix and G x is the coupling function. Notice that in the above equation, represents the element-wise product, and we have generalised the global coupling term σ into a vector, a generalisation that will prove useful in the following. We start by considering the same synchronised equilibrium state for all the nodes x * ( ), which is stable in the absence of coupling. Then by perturbing this state by a small random term x ( ) = x * ( ) + ξ ( ) we can obtain at the first perturbative order the variational equations ξ = DF x * ( ) ξ ( ) + σ =1 DG x * ( ) ξ ( ). Now starting from here and assuming that the linearized spatial operator DG x * ( ) is diagonalizable, we can proceed further to obtain the Master Stability Function (MSF) 2,27 , which relates the stability of the uniformly synchronized state x * to the topology of the connection network. So, assuming that there is a linearly independent basis of eigenfunctions(-vectors) of the spatial operator, we can obtain the diagonalized decoupled equations where by Λ ( ) we have denoted the -th eigenvalue of the spatial operator. To concretize our analysis and without loss of generality, throughout this text, we will consider a diffusive coupling through the Laplacian matrix L whose entries are L = − with = =1 . The formalism described so far allows the local analysis of the linear stability of the originally uniform state. Depending on the nature of the latter, eq. (2) might be non-autonomous when the uniform state is a limit cycle or autonomous when it is a fixed point. The traditional approach in the study of synchronization phenomenon is based on the former case when the system is already in a synchronized limit cycle manifold. In this case, the aim is to keep as stable as possible such a state since an eventual instability would cause the desynchronization of the whole system. As already anticipated, we will adopt here an alternative approach that to the best of our knowledge has not been used so far in the vast literature related to the synchronization dynamics. The first major difference is that we start from an unstable uniform equilibrium of a fixed point, so x * ( ) = x * for all the nodes. Starting from a uniform stationary state will allow the use of the spectral analysis instead of the classical Lyapunov exponent used for non-autonomous systems 2,27 . This makes possible the necessary analytical treatment for finding the instability conditions and the description of the orbits' behavior in the initial linear regime. In the latter setting the MSF formalism is known as the dispersion relation 19 . Since a homogeneous stationary equilibrium lacks the oscillatory behaviour that characterises synchronization dynamics, we will recover final time-varying patterns by considering a nonlinear function constituted by 3 variables, known in the biomathematical literature as species 19 ] with a little abuse of notation. In other words, we have It is well known that a 3 species system is a sufficient requirement to have an oscillatory instability and, as a consequence, nonlinear oscillatory patterns 18,19 . This is the key factor in producing the global oscillations of our system, as is shown in Fig. 1 panels 1 ) − 3 ). As an alternative to this model, oscillatory patterns can be obtained also for a 2 species model in a directed network where the oscillations are as a consequence of the complex spectrum of the Laplacian matrix 29,30 (see Appendix B for a detailed description). When neither of the two aspects above contributing to oscillations occurs, then the instability is stationary, and modular Turing patterns can emerge 33,34 . We should, however, emphasise that in general, in a complex network, the oscillations appear to be asynchronous following the shape of a traveling front in the limit case of a regular lattice 18,28 (see Fig. 1 1 ) and 1 )). To understand the rationale behind such behavior, we should look for the shape of the eigenvectors, which describe the linear evolution of the orbits in a short time and also determine the general shape of the final nonlinear pattern. In fact in order to solve the linear (decoupled) variational equations (2) we look for solutions of the form ζ = =1 Φ ( ) where and Φ ( ) represent the -th growth mode and Laplacian eigenvector respectively. So, in the initial regime, once there is at least one value of for which the real part of , the dispersion relation defined earlier, is positive, as is the case for Fig.   1 1 ), then the entries of vectors of the 3 species will be scaled by the entries of the eigenvector Φ ( ) . Also, in a 3 species model, the imaginary part is different from zero allowing an oscillatory instability to occur 18 . In regular multi-dimensional lattices, such an eigenvector represents a discrete term of multi-dimensional Fourier series. This is why in such discrete domains, the final patterns will be a combination, depending on the number of the unstable (Fourier) modes, of time-varying sinusoidal shaped traveling waves 18,19,28 .

A. Network construction for clustered eigenvectors
Having briefly introduced the mechanism responsible for shaping the stable nonlinear pattern, we next describe how to construct the coupling network to obtain clusters of coherent and incoherent nodes ? . We start considering the set of eigenvectors Φ (0) , Φ (1) , . . . , Φ ( ) of the Laplacian of a connected and undirected network sorted in decreasing order of the respective modes (eigenvalues). From the algebraic connectivity, we know that the network Laplacian will have a single zero eigenvalue Λ (0) = 0 whose eigenvector is uniformly distributed, for simplicity here assumed to be the scaled unit vector is the scaling constant. Let us remark here that from eq. (2) the null mode Λ (0) = 0 represents the a-spatial case meaning that the nodes should be stable once uncoupled. So we should look for the eigenvectors corresponding to the nonnull modes to destabilize our system globally. To obtain eigenvectors with clusters of entries that are close to each other, in total analogy with the equilibrium state, we proceed by construction. We start by considering networks with more than one null mode, for example, such eigenvalues 1 Λ (0) = · · · = Λ (0) = 0. This means that our network has disconnected components, here referring to as disconnected clusters. The Laplacian matrix and its eigenvectors in this case are respectively a block matrix and block vectors L = [ 1 L, . . . , L, . . . , L] and Φ ( ) = 1 Φ ( ) , . . . , Φ ( ) , . . . , Φ ( ) . Here, we have divided the nodes into disconnected clusters with the same number of nodes. Observe that the eigenvector of the -th null mode can include different solutions, but here we chose purposely the one that has constant entries for the blocks of each cluster Instead for the eigenvectors of the non null modes the only possible solution is where Φ ( ) is the -th eigenvector of the -th disconnected cluster. We invite the interested reader to refer to Appendix A for a rigorous derivation of this assertion. For the next step, we weakly connect the    this means that the former spectrum will be weakly deformed. The first indication of this deformation is that due to the connectedness, the graph Laplacian will have only one zero eigenvalue. All the remaining − 1 formerly zero eigenvalues will be different from zero, but very small in value. Also, in terms of equilibrium states, now we have the only one eigenvector with uniformly distributed entries. So the question that arises naturally is what happened to the remaining former − 1 equilibrium eigenvectors? Their entries will be slightly deformed, too, so the − 1 eigenvectors will have almost uniform entries per block. This essential spectral property of weakly disconnected subgraphs will be the basis for developing our method for the emergence of synchronized clusters and cluster chimera states.

III. EMERGENCE OF CLUSTER SYNCHRONIZED AND CHIMERA PATTERNS
As a result of the above reconstruction procedure, we have obtained particular network structures that contain potential candidates of unstable modes with eigenfunctions that have entries within clusters very close to each other. Such networks indeed exist in the literature, and they are known as modular networks 3 . They are widespread in many biological and social settings 3 , making potentially either the cluster synchronization or the cluster chimera states typical emergent behaviors in a wide class of networks. In Fig. 1 we show how the oscillations in a 3 species reaction-diffusion system change from initially asynchronous (in the left column panels) to synchronous per cluster when the modularity increases (central and right columns panels). The pattern is organized in a traveling wave on the ring network where the modularity is absent but where the network is organized in distinct modules, the species segregate their phases following the clusters they belong to. The explanation for that can be found in the eigenvec-tors of the first − 1 non-zero eigenvalues, previously introduced, that we denote here as modular eigenvalues. The first of them that we have also selected to be unstable in Fig. 1 1 ) is known as the Fiedler eigenvalue, and it is an indicator of the algebraic connectivity of the network 3 . As shown in the previous section, the modular eigenvectors have the unique property of being segregated per module, a feature that, when activated by the corresponding unstable modular eigenvalue, is reflected in the final nonlinear pattern. Besides the modular eigenvalues, we also have the non-modular ones, the remaining − eigenvalues. Following the perturbative analysis developed so far, it is not possible that non-modular eigenvectors have segregated nonnull entries per module. This means that if a non-modular eigenvalue would be the only unstable mode then the pattern of such setting would be irregular for a given module and almost not expressed for the rest of them. The linear stability principle, upon which we base our method, allows the extension of pattern selection beyond the single unstable mode case. In fact, if two or more modes are unstable and near the bifurcation threshold, they will compete with each other in the linear regime until the nonlinearities stabilize the pattern shaping. So one should expect that the patterning outcome will reflect to some extent the shape of the eigenvectors of all the unstable modes involved. To show such behaviour, in Fig. 2 we have taken into consideration only two unstable eigenvalues, a modular and a non-modular one. The choice has been such that in the dispersion relation representation, the two modes have comparable real parts. In this way the two different behaviours will be present to similar levels. In fact, as can be seen from the pattern snapshot taken at some given time, Fig. 2 ), the second module is non-homogeneous in the concentration (of the first species in this case). Consequently, in panel ) we can observe that this oscillatory pattern will constitute a new state that we name as a cluster chimera. Let us note that the choice for the second module to be incoherent is simply due to the intraand inter-connections, but other scenarios where different or more incoherent modules are also possible, as can be seen in Appendix B. Although the results presented in both previous figures have been validated for the case of modules of the same size, in the following we show that our method is robust even for the case of networks with modules of considerably different in size.

IV. ROBUSTNESS OF THE SYMMETRY-BREAKING METHOD FOR MODULES OF DIFFERENT SIZES
The method that we proposed throughout this paper has been successful in the generation and prediction of the cluster synchronization, and the cluster chimera states. In this section will discuss the robustness of the symmetry-breaking mechanism when the modules' size is considerably different. We start by considering a simple network model of topology of 200 nodes divided into two modules with a number of inter-edges between the modules few on average compared to the intra-edges ones. Keeping a small number of inter-edges compared to the number of intra-edges is a necessary requisite for the validity of the perturbative approach used to justify our method. Initially, the two modules have 50 nodes each of as shown in panels 1 ) − 4 ) of Fig. 3, and we chose the parameters in such a way to obtain a (weak) cluster chimera where module I is weakly incoherent due to a mix of unstable modular and non-modular eigenvectors. Clearly, in this setting, the results are shown throughout this work stand firmly, as confirmed by the clear gap between the modular and non-modular set of eigenvalues, panel 1 ). Then we tune the network topology by keeping the same total number of links but changing the size of the modules in the 150 vs. 50 ratio. We notice that the first two non-zero eigenvalues' position slightly shifts to the right, panel 2 ), so that the modular eigenvalue increases the spectral gap from the origin yielding a cluster synchronized pattern. When we raise the ratio to 190 vs. 10 nodes, respectively, per module, we can observe a neat shift of the modular eigenvalue and a decrease of the gap with the non-modular one, panel 3 ).The resulting pattern is made by two synchronized clusters where one of them is slightly incoherent. Such behavior is understandable from the spectral perturbative perspective. In fact, when we have smaller modules, i.e., of 10 nodes as in our case, the action of the unitary weighted interconnection links (that remain invariant in number in all the three examples) over the small module can be considered less perturbative compared to the previous scenarios when such module was larger, i.e., of 100 nodes. Nevertheless, the results of Fig. 3, in particular panels 2 ), 2 ), and 2 ) show that the shape of the eigenvectors remains still sufficiently segregated (proportionally to the amount of modularity of the network) to yield an explicit cluster synchronization or chimera states making it an obvious indicator of the robustness of our method.

V. DISCUSSION AND CONCLUSIONS
To conclude, in this paper we have studied the problem of the emergence of a new class of chimera states, which we call cluster chimeras, following a cluster synchronization mechanism based on the modular structure of the underlying network. To understand the mechanism that generates these states we have developed a theoretical framework based on a symmetry-breaking principle. Due to a global instability the system departs from a static homogeneous state and the perturbations associated to each node will be shaped according to the eigenvectors corresponding to the unstable modes of the Laplacian matrix. This linear regime behaviour is saturated by the nonlinear terms in the final stage, inheriting the eigenvector characteristics in the spatially extended pattern. As a first step, we show that is possible to reconstruct networked structures that have particular eigenvectors with certain features: they are organised in groups of entries close to each other. Based in a perturbative approach, we prove that these networks have − 1 such eigenvectors, denoted as modular ones, where is the number of their clusters. On the other hand the remaining − non equilibrium eigenstates, here called non-modular ones, present irregular shapes. At this point we could select (at least) two unstables eigenvalues, one modular and one non-modular in the Master Stability Function such that due to their competition during the pattern forma-tion process the final result yields a cluster chimera state.
On the other hand, if the unstable eigenvalues are exclusively modular ones, then the resulting pattern is made of synchronized clusters. The pattern formation mechanism presents significant advantages; first and in contrast to classical chimeras, the generation process is robust to the choice of initial conditions. Secondly, patterns that follow such a global instability principle are known to be stable and not transient, which is also observed in our numerical simulations. These two aspects are are the main open questions regarding the process through which chimera states emerge. As part of the novelty, in this work we have particularly studied the case when the coupled dynamical systems are fixed points instead of the usual (nonlinear) oscillators. The main reason for that was because in this way it was possible to make use of the solid classical results already existing in the theory of pattern formation. Nevertheless we believe that the results shown here can be also extended to the case of coupled oscillators. We are confident that our results will shed light on a better understanding of the generation mechanisms and stability of chimera states and also fill the theoretical gap of a mathematical explanation for recent experimental observations in this domain 26 . However, this choice is not possible since it does not respect the perturbative approach once the clusters are weakly connected. To prove that we consider the simple case of 2 clusters, but the same easily follows for the case of clusters. The perturbative eigenvalues problem for the canonical eigenvector is now written as where E 12 , E 21 are the intra-cluster weak coupling, 1 , 2 the perturbation of the uncoupled eigenvectors and the pertubation of the null mode. From this relation we obtain the following equalities .
Clearly since the last relation is not possible this choice for the eigenvector is not permitted. On the contrary, for the eigenvectors of the non null modes the only possible solution is where Φ ( ) is the -th eigenvector of the -th disconnected cluster. In fact, if for the 2 cluster case we consider as candidate Φ ( ) , 1 , where instead of 1 we could chose any other nonzero vector, from the same reasoning as above we have This is equivalent to which again cannot be true. This concludes that the only possible choices for the Laplacian eigenvectors of a network with disconnected clusters are defined as above. In this section we will consider the case when the cluster synchronisation and chimera states emerge due to the modularity of directed networks. A major contribute of directed networks is due to the fact that the spectrum of the Laplacian matrix is complex which cause a Turing-Hopf instability in the linear regime, known as topology-driven instability 29 , yielding this way oscillatory patterns. Let us remind here that such patterns are not otherwise possible for a 2 species reaction-diffusion system. Also the directedness of the networked structure contributes in the increment of chances for emergence of patterns due to a higher dispersion relation compared to the symmetric network or the continuous domain. For further details the interested reader should refer to 29,30 . It is important also to remind that, in general, a directed networked structure might be also strongly nonnormal as observed in many real-world networks 35 . In this case the pattern formation process follows different mechanisms induced by the non-normality of the Laplacian matrix 36 . However, we decided to skip such scenario that goes beyond the scope of this paper and focus on the topology-driven case only.
We start by considering a 2 species activator-inhibitor model of a reaction-diffusion system. The concentration for the two species are denoted , for the activator and , for the inhibitor and where the index refers to one of the nodes. The corresponding equations take the following form: where and are the diffusion constants and (·, ·), (·, ·) are the two nonlinear functions representing the interaction of the species inside each node . The spatial interactions on the other side, are represented by the diffusion operator with entries L which are a reflection of the topology of the networked support.
To proceed with the stability analysis, we linearise the reaction-diffusion system (B1) around the uniform steady state (u * , v * ): where here u, v are the vectors of the perturbations. The Jacobian matrix for the reaction terms, evaluated at the equilibrium point is given by J =  I  I  I  I  , and the diffusion one for both species is D = where O is the × zero-valued matrix. We consider the basis of the eigenvectors Φ , with = 1, . . . , , of the Laplacian operator L to the corresponding eigenvalues Λ . Due to the asymmetry of the adjacency matrix A, the spectrum Λ is, in principle, complex. In particular is has been shown in 29,30 that for balanced directed networks where the numbers of incoming edges is equal to the outgoing ones = , such a basis always exists. This will be also the case in the following discussion. In order to solve this linearised system, we expand the perturbations on the basis of the eigenvectors, i.e. = =1 Φ and = =1 Φ . Following this, it is possible to reduces the 2 × 2 system to a 2 × 2 eigenvalue problem, for each value of the node index = 1, . . . , : If the real part of , known as the dispersion relation, has at least a mode Λ ≠0 for which takes positive values, then the steady state becomes unstable. Mathematically this is formalised as: In the above expression we have denoted by sgn(·) the sign function and () , () indicate, the real and imaginary parts of the respective function. To proceed further with the calculations, we make use of the definition of a square root of a complex number. By taking = + where = √ −1 is the imaginary unit, then The instability sets in when | (trJ ) | ≤ , a condition that can be expressed by the following inequality 29 : where (Λ , Λ ) span the complex plane where lie the eigenvalues of the Laplacian operator. The explicit form of the polynomials 1 , 2 is described by: 1 ( ) = 14 4 + 13 3 + 12 2 + 11 + 10 2 ( ) = 22 2 + 21 + 20 .
The above system undergoes a Turing-Hopf instability, induced by the (directed) topology of the interaction network. Such mechanism responsible for the emergence of pattern in directed networks, is known as the topology-driven instability 29 . In Fig. 4 we show the cluster synchronisation phenomenon in a directed modular network where the pattern are triggered by the directedness of the underlying network and the oscillations are as a consequence of a complex spectrum of the Laplacian operator. If apart the unstable modular mode, we could chose also to destabilise one of the non-modular ones then chimera states would appear in this case as shown in Fig. 5.  a) The dispersion relation shows that a modular and a non-modular eigenvalue are simultaneously unstable. In the inset are given the real and the imaginary part of the dispersion relation. b) The chimera states evolution is shown as a function of time where the last module is incoherent compared to the remaining four ones.