Symmetry-breaking mechanism for the formation of cluster chimera patterns

Full understanding of the dynamics of complex systems is probably one of the biggest challenges of physics this century. The emergence of collective behavior 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,19,20 the simultaneous flashing of male fireflies during the mating season,²¹ and the synchronization of the alternating current in power grids.²² 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,²³ the synchronization of chaotic oscillators and the phenomenon of cluster synchronization when the set of oscillators is divided into subsets of synchronized ones,¹⁵ or the emergence of exotic states where identically coupled oscillators synchronize in clusters.²⁴ 

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 (Usually by state in the dynamical systems theory, it is referred to as an equilibrium solution. In this sense, the long transient of the chimera behavior makes it effectively a state, although its stability might not be truly asymptotic.), and their sensitivity to the initial conditions that make their experimental verification difficult to achieve.^4,9,25

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, especially, 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.⁷ So far, however, there is no exhaustive answer to the stability of chimera states^5,6 or a robust mechanism that explains how chimeras are generated.¹² 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 sets of nodes that can be stable or unstable.^10,11 Other recent alternatives have been based on symmetry-breaking principles where chimeras emerge from an earlier uniformly synchronized regime.^12,13 Aside from theoretical results, experimental evidence²⁶ proves the existence of cluster synchronization and chimera states for the case of coupled Belousov–Zhabotinsky oscillators.¹⁶ 

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^16–18 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 1970s, 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,25 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.^17,18

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²⁹ or multispecies interactions,³² 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. It is also worth mentioning that network modularity has already been found to be a good setting where chimera states can be observed.⁶ Nevertheless, the chimera states studied in Ref. 6, as is the case of chimera states, in general, suffer from a high dependence on the initial conditions. In fact, it has been shown that only a few of the randomly chosen initial conditions converge to a chimera state.^2,9,25 To complement this finding and, remarkably, in addition to it, we propose a systematic method, based on pattern formation theory, to obtain chimeras robust to initial conditions, an important aspect lacking in the paper above and the related literature.

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.

We start by considering the general setting of a network of coupled dynamical systems through the equations

where each of the $N$ individuals, identified as nodes of the network, is represented by $M$-dimensional vector state variables $xi$⁠.² The local dynamics of the $i$th node is given by the nonlinear function $F(xi)$⁠, $σ$ is the vector of the coupling strength, $Wij$ are the weighted entries of the adjacency matrix, and $G(xj)$ is the coupling function. Notice that in the above equation, $⊙$ represents the element-wise product, and we have generalized the global coupling term $σ$ into a vector, a generalization that will prove useful in the following. We start by considering the same synchronized equilibrium state for all the nodes $x∗(t)$⁠, which is stable in the absence of coupling. Then, by perturbing this state by a small random term $xi(t)=x∗(t)+ξi(t)$⁠, we can obtain at the first perturbative order the variational equations $ξ˙i=DF(x∗(t))ξi(t)+σ⊙∑i=1NWijDG(x∗(t))ξj(t)$⁠. Now, starting from here and assuming that the linearized spatial operator $DG(x∗(t))$ 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. Therefore, 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 diffusive coupling through the Laplacian matrix $L$ whose entries are $Lij=Wij−kiδij$ with $ki=∑j=1NWij$⁠.

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 a 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; therefore, $xi∗(t)=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.¹⁸ Since a homogeneous stationary equilibrium lacks the oscillatory behavior that characterizes synchronization dynamics, we will recover final time-varying patterns by considering a nonlinear function constituted by three variables, known in the biomathematical literature as species¹⁸ $F(xi)=[f(xi),g(yi),h(zi)]$⁠. For each node, we write the three dimensional state vector $xi=[ui,vi,zi]$ with a little abuse of notation. In other words, we have

It is well known that a three species system is a sufficient requirement to have an oscillatory instability and, as a consequence, nonlinear oscillatory patterns.^17,18 This is the key factor in producing the global oscillations of our system, as is shown in Fig. 1 [panels (a$1$⁠)–(a$3$⁠)]. As an alternative to this model, oscillatory patterns can be obtained also for a two 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, emphasize 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^17,28 [see Figs. 1(c₁) and 1(d₁)]. 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 variational equations, we look for solutions of the form $ξi=∑α=1NcαeγαΦi(α)$⁠, where $γα$ and $Φ(α)$ represent the $α$th growth mode and the Laplacian eigenvector, respectively. Therefore, 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(a1), then the entries of vectors of the three species will be scaled by the entries of the eigenvector $Φ(α)$⁠. Also, in a three species model, the imaginary part is different from zero allowing an oscillatory instability to occur.¹⁷ 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 resulting patterns will be a combination, depending on the number of the unstable (Fourier) modes, of time-varying sinusoidal shaped traveling waves.^17,18,28

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 want to note that the concept of nodes being in synchrony with each other is not binary, either fully synchronized or completely asynchronous, and it is quantified by an order parameter (initially introduced by Kuramoto^1,2), which varies from zero to one.]. We start considering the set of eigenvectors $[Φ(0),Φ(1),…,Φ(N)]$ 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 $Φ(0)=a1=[a,a,…,a]T$⁠, where $a$ 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. Therefore, 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, $C$ such eigenvalues $1Λ(0)=⋯=CΛ(0)=0$⁠. This means that our network has $C$ 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=diag[1L,…,nL,…,CL]$ and $Φ(α)=[1Φ(α),…,nΦ(α),…,CΦ(α)]$⁠. Here, we have divided the $N$ nodes into $C$ disconnected clusters with the same number of nodes. Observe that the eigenvector of the $n$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 $nΦ(α)$ is the $α$th eigenvector of the $n$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 $C$ disconnected clusters through a minimal number of links. From a purely perturbative point of view, 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 $C−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 $C−1$ equilibrium eigenvectors? Their entries will be slightly deformed, too; therefore, the $C−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.

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.³ They are widespread in many biological and social settings,³ making potentially either the cluster synchronization or the cluster chimera states’ typical emergent behaviors in a wide class of networks. It is also important to emphasize that modularity is the only criterion requested for having clustered Laplacian eigenvectors, but no other characteristics on the modules’ structure are needed. Throughout this paper, modules are identically generated, in the sense that they are outcomes from the same ensemble where we have fixed only the number of modules, their nodes’ number, and the (average) number of intra- and inter-edges. In Fig. 1, we show how the oscillations in a three 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). Notice that the patterns shown throughout this paper have been obtained integrating the system for $106$ integration steps with a time step of $10−2$⁠. In the temporal evolution of the patterns [panels (c) in Figs. 1–3], the transient will die out soon so that we have shown only the last $500$ or 1000 time steps when the patterns have already reached the equilibrium. 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 such behavior can be found in the eigenvectors of the first $C−1$ non-zero eigenvalues, previously introduced, that we denote here as modular eigenvalues. The first of them that we have selected to be unstable in Fig. 1(a₁) is known as the Fiedler eigenvalue, and it is an indicator of the algebraic connectivity of the network.³ As shown in Sec. II, 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 $N−C$ 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. Therefore, 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 behavior, 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 behaviors will be present to similar levels. In fact, as can be seen from the pattern snapshot taken at some given time, Fig. 2(b), the second module is non-homogeneous in the concentration (of the first species in this case). Consequently, in panel (c), 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 intra- and inter-connections, but other scenarios where different or more incoherent modules are also possible, as can be seen in  Appendix B. Also, both modular and non-modular modes must be unstable because the former is responsible for generating the oscillations in the synchronized clusters and the latter for destroying the coherence in one (or more) modules. 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.

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, we 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 $NW$ topology of $200$ nodes divided into two modules with a number of inter-edges between the modules few on average compared to the intra-edge 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 $100$ nodes each as shown in panels (a$1$⁠)–(a$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 shown throughout this work stand firmly, as confirmed by the clear gap between the modular and non-modular set of eigenvalues, panel (a$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 (a$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 (a$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 a module was larger, i.e., of $100$ nodes. Nevertheless, the results of Fig. 3, in particular, panels (a$2$⁠), (b$2$⁠), and (c$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.

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. Our approach drastically differs from the standard one widely spread in the literature, where synchronized oscillations arise due to coupled oscillators. Alternatively, we bypass such restriction showing that coherent or incoherent clusters of synchronized oscillations are also possible when the system initially sets on a uniform fixed-point state. In fact, due to a global instability, the system departs from a static homogeneous state and the perturbations associated with each node will be shaped according to the eigenvectors corresponding to the unstable modes of the Laplacian matrix. This linear regime behavior 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 organized in groups of entries close to each other. Based on a perturbative approach, we prove that these networks have $C−1$ such eigenvectors, denoted as modular ones, where $C$ is the number of their clusters. On the other hand, the remaining $N−C$ non equilibrium eigenstates, here called non-modular ones, present irregular shapes. At this point, we could select (at least) two unstable eigenvalues, one modular and one non-modular in the master stability function such that, due to their competition during the pattern formation 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. Therefore, modular networks are ideal candidates for possessing clustered eigenvectors of the Laplacian, making our method immediately applicable to this particular setup where chimera states have been numerically observed.^35,36 However, in principle, the proposed technique can extend to other families of networks with spectral clustering features. Let us also notice that the symmetry breaking in our model is mainly due to the model’s parameters and its dynamical properties and not simply due to the modular structure of the network. For instance, in Fig. 1 (left column panels), symmetry breaks even in a ring, although there is a total lack of modularity in this case. If symmetry does not break, i.e., there are no unstable modes in the dispersion relation, then the system remains in a stationary homogenous state, independent of whether the network is modular or the modules are identical or of different sizes.

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.²⁵ The only requirement in this sense is that the systems must be initiated with a random distribution of initial values centered around the stable fixed-point of the aspatial (reactions) system. The latter does not represent a restriction (it varies continuously with the parameters) since it influences neither the spatial nor the temporal properties of the chimera patterns, which depend exclusively on the coupling topology and strengths. Second, patterns that follow such a global instability principle are well known in the literature^17,18 to be stable and not transient, which is also observed in our numerical simulations. These two aspects 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 also be 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.²⁶ 

B.A.S. acknowledges funding from the Irish Research Council under Grant No. GOIPG/2018/3026. The work of J.P.G. and M.A. is partly funded by the Science Foundation Ireland (Grant Nos. 16/IA/4470, 16/RC/3918, and 12/RC/2289 P2) and cofunded under the European Regional Development Fund. A.A. acknowledges financial support from the Spanish MINECO (Grant No. PGC2018-094754-B-C21), the Generalitat de Catalunya (Grant No. 2017SGR-896), the Universitat Rovira i Virgili (Grant No. 2019PFR-URV-B2-41), the Generalitat de Catalunya ICREA Academia, and the James S. McDonnell Foundation (Grant No. 220020325).

The authors have no conflict of interest to declare.

M.A. and B.A.S. contributed equally to this work. All authors discussed the research and approved the manuscript.

The data that support the findings of this study are available within the article.

Starting from a perturbative setting is possible to show that the eigenvector of the $n$th null mode of $C$ disconnected clusters can select only one of the possible solutions, the one that has constant entries of alternating signs for the blocks of each cluster,

In fact, the most intuitive eigenvector of the null mode is the canonical one,

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 two clusters, but the same easily follows for the case of $C$ clusters. The perturbative eigenvalue problem for the canonical eigenvector is now written as

where $E12$⁠, $E21$ 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 $nΦ(α)$ is the $α$th eigenvector of the $n$th disconnected cluster. In fact, if for the two cluster case, we consider as candidate $[nΦ(α),1]T$⁠, where instead of $1$⁠, we could chose any other nonzero vector, from the same reasoning as above, we have

therefore,

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 synchronization and chimera states emerge due to the modularity of directed networks. A major contribution of directed networks is due to the fact that the spectrum of the Laplacian matrix is complex, which causes a Turing–Hopf instability in the linear regime, known as topology-driven instability,²⁹ yielding this way oscillatory patterns. Let us remind here that such patterns are not otherwise possible for a two species reaction–diffusion system. Also, the directedness of the networked structure contributes in the increment of chances for the 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 Refs. 29 and 30. It is important also to remind that, in general, a directed networked structure might also be strongly non-normal as observed in many real-world networks.³⁷ In this case, the pattern formation process follows different mechanisms induced by the non-normality of the Laplacian matrix.³⁸ However, we decided to skip such a scenario that goes beyond the scope of this paper and focus on the topology-driven case only.

We start by considering a two species activator–inhibitor model of a reaction–diffusion system. The concentration for the two species is denoted as $ui$ for the activator and $vi$ for the inhibitor and where the index $i$ refers to one of the $N$ nodes. The corresponding equations take the following form:

where $Du$ and $Dv$ are the diffusion constants and $f(⋅,⋅)$⁠, $g(⋅,⋅)$ are the two nonlinear functions representing the interaction of the species inside each node $i$⁠. The spatial interactions on the other side are represented by the diffusion operator with entries $Lij$⁠, which are a reflection of the topology of the networked support.

To proceed with the stability analysis, we linearize 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

and the diffusion one for both species is

where $O$ is the $N×N$ zero-valued matrix. We consider the basis of the eigenvectors $Φiα$⁠, with $α=1,…,N$⁠, 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, it has been shown in Refs. 29 and 30 that for balanced directed networks, the numbers of incoming edges is equal to the outgoing ones $kiout=kiin$ such that a basis always exists. This will also be the case in the following discussion. In order to solve this linearized system, we expand the perturbations on the basis of the eigenvectors, i.e., $δui=∑α=1NbαΦiα$ and $δvi=∑α=1NcαΦiα$⁠. Following this, it is possible to reduce the $2N×2N$ system to a $2×2$ eigenvalue problem for each value of the node index $α=1,…,N$⁠,

where $Jα$ is the $2×2$ extended Jacobian, i.e., $Jα≡J+DΛα$ where $D=diag(Du,Dv)$⁠, $Λα=diag(Λα,Λα)$⁠. If the real part of $λα$⁠, known as the dispersion relation, has at least a mode $Λα≠0$ for which it takes positive values, then the steady state becomes unstable. Mathematically, this is formalized as

where

In the above expression, we have denoted by $sgn(⋅)$ the sign function and $()Re$⁠, $()Im$ 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 $z=a+bι$⁠, where $ι=−1$ is the imaginary unit, then

The instability sets in when $|(trJα)Re|≤γ$⁠, a condition that can be expressed by the following inequality:²⁹ 

where $(ΛRe,ΛIm)$ span the complex plane where lie the eigenvalues of the Laplacian operator. The explicit form of the polynomials $S1$⁠, $S2$ is described by

The above constants are given by

The above system undergoes a Turing–Hopf instability, induced by the (directed) topology of the interaction network. Such a mechanism responsible for the emergence of patterns in directed networks is known as the topology-driven instability.²⁹ In Fig. 4, we show the cluster synchronization phenomenon in a directed modular network where the patterns 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 destabilize one of the non-modular ones, then chimera states would appear in this case as shown in Fig. 5.