The transition to synchronization of networked systems

We study the synchronization properties of a generic networked dynamical system, and show that, under a suitable approximation, the transition to synchronization can be predicted with the only help of eigenvalues and eigenvectors of the graph Laplacian matrix. The transition comes out to be made of a well defined sequence of events, each of which corresponds to a specific clustered state. The network’s nodes involved in each of the clusters can be identified, and the value of the coupling strength at which the events are taking place can be approximately ascertained. Finally, we present large-scale simulations which show the accuracy of the approximation made, and of our predictions in describing the synchronization transition of both synthetic and real-world large size networks, and we even report that the observed sequence of clusters is preserved in heterogeneous networks made of slightly non-identical systems.

We study the synchronization properties of a generic networked dynamical system, and show that, under a suitable approximation, the transition to synchronization can be predicted with the only help of eigenvalues and eigenvectors of the graph Laplacian matrix.The transition comes out to be made of a well defined sequence of events, each of which corresponds to a specific clustered state.The network's nodes involved in each of the clusters can be identified, and the value of the coupling strength at which the events are taking place can be approximately ascertained.Finally, we present large-scale simulations which show the accuracy of the approximation made, and of our predictions in describing the synchronization transition of both synthetic and real-world large size networks, and we even report that the observed sequence of clusters is preserved in heterogeneous networks made of slightly non-identical systems.
From brain dynamics and neuronal firing, to power grids or financial markets, synchronization of networked units is the collective behavior characterizing the normal functioning of most natural and man made systems [1][2][3][4][5][6][7].As a control parameter (typically the coupling strength in each link of the network) increases, a transition occurs between a fully disordered and gaseous-like phase (where the units evolve in a totally incoherent manner) to an ordered or solid-like phase (in which, instead, all units follow the same trajectory in time).
The transition between such two phases can be discontinuous and irreversible, or smooth, continuous, and reversible.The first case is known as Explosive Synchronization [8], which has been described in various circumstances [9][10][11][12][13], and which refers to an abrupt onset of synchronization following an infinitesimally small change in the control parameter, with hysteresis loops that may be observed as in a thermodynamic first-order phase transition.The second case is the most commonly observed one, and corresponds instead to a second-order phase transition, resulting in intermediate states emerging in between the two phases.Namely, the path to synchrony [14] is here characterized by a sequence of events where structured states emerge made of different functional modules (or clusters), each one evolving in unison.This is known as cluster synchronization (CS) [15][16][17], and a lot of studies pointed out that the structural properties of the graph are responsible for the way nodes clusterize during CS [18][19][20][21].In particular, it was argued that the clusters formed during the transition are to be connected to the symmetry orbits and/or to the equitable partitions of the graph [19,22].Precisely, partitions of the network's nodes into orbits (for the whole, or a subgroup, of the automorphism group of the graph) or, more generally, equitable partitions, have been identified as potential clusters that can synchronize before the setting of complete synchronization in the network [23,24].The fact that a given graph partition can sustain cluster synchronization of its cells is different, however, from the problem faced in our study i.e., that of determining whether it will indeed be realized for some range of the coupling strength parameter and, more importantly, in which order do clusters synchronize during the transition between the two main phases (i.e., as the coupling strength increases from zero).
On the other hand, clusters due to symmetries are found rather ubiquitously, as most of real-world networks have a large number of localized symmetries, that is, a large number of small subgraphs called symmetric motifs [25,26], where independent (formally, supportdisjoint) symmetries are generated, with each symmetric motif made of one or more orbits.For instance, Ref. [26] analyzes several real-world networks and, in all cases, gives evidence of a large number of symmetric motifs: from 149 motifs in the protein-protein interaction network of the yeast (a network of 1,647 nodes) to over 245,000 motifs in the LiveJournal social network (more than 5 million nodes).More evidences of symmetries in real-world networks can be found in Refs.[27][28][29].
In our work we assume that, during the transition, the synchronous solution of each cluster does not differ substantially from that of the entire network and, under such an approximation, we elaborate a practical technique which is able to elucidate the transition to synchronization in a generic network of identical systems.Namely, we introduce a (simple, effective, and limited in computational demand) method which is able to: i) predict the entire sequence of events that are taking place during the transition, ii) identify exactly which graph's node is belonging to each of the emergent clusters, and iii) provide a well approximated calculation of the critical coupling strength value at which each of such clusters is observed to synchronize.We also demonstrate that, under the assumed approximation, the sequence of events is in fact universal, in that it is independent on the specific dynamical system operating in each network's node and depends, instead, only on the graph's structure.Our study, moreover, allows to clarify that the emerging clusters are those groups of nodes which are indistinguishable at the eyes of any other network's vertex.This means that all nodes in a cluster have the same connections (and the same weights) with nodes not belonging to the cluster, and therefore they receive the same dynamical input from the rest of the network.As such, we prove that synchronizable clusters in a network are subsets more general than those defined by the graph's symmetry orbits, and at the same time more specific than those described by equitable partitions (see the Supplementary Information for a detailed description of the exact relationship between the observed clusters and the graph's symmetry orbits and equitable partitions).Finally, we present extensive numerical simulations with both synthetic and real-world networks, which demonstrate the high accuracy of our approximations and predictions, and also report on synchronization features in heterogeneous networks showing that the predicted cluster sequence can be maintained even for networks made of non-identical dynamical units.

The synchronization solution
The starting point is a generic ensemble of N identical dynamical systems interplaying over a network G.The equations of motion are where x i (t) is the m-dimensional vector state describing the dynamics of each node i, f : R m −→ R m describes the local (identical in all units) dynamical flow, d is a real-valued coupling strength, L ij is the (ij) entry of the Laplacian matrix associated to G, and g : R m −→ R m is the output function through which units interact.L is a zero-row matrix, a property which, in turn, guarantees existence and invariance of the network's synchronized solution x s (t) = x 1 (t) = . . .= x N (t).The necessary condition for the stability of such solution can be assessed by means of the Master Stability Function (MSF) approach, a method initially developed for pairwise coupled systems [30], and later extended in many ways to heterogeneous networks [31], and to timevarying [32] and higher-order [33] interactions.When possible, the MSF formalism can be complemented by other global approaches (such as the Lyapunov functionals) which provide sufficient conditions for stability.As the MSF is of rather standard use, all details about the associated mathematics are contained in the Methods section.Notice, moreover, that our Eq.( 1) is in fact written in the simplest as possible form for the easiness of the mathematical passages that will be described below.However, our approach can be extended straightforwardly (see details in Refs.[33,34] and references therein) to the much more general case of a coupling term given by d N j=1 A ij g(x i , x j ), where A ij is the (ij) entry of the graph's adjacency matrix, under only two assumptions: (i) the Laplacian matrix L associated to A is diagonalizable, and the set of eigenvectors form an orthonormal basis of R N , and (ii) the coupling function g(x i , x j ) is synchronization non-invasive i.e., it strictly vanishes at the synchronization solution (g(x s , x s ) = 0).
The full details behind the MSF approach can be found in the Methods section.We here concisely summarize them to pave the way for the cluster synchronization analysis which will follow.In essence, one considers perturbations around the synchronous state, and performs linear stability analysis of Eq. (1).Due to the properties of L, the perturbations can be expanded as linear combinations of its eigenvectors v i (which span the space T tangent to the synchronization manifold M).The coefficients η i of this expansion obey variational equations involving the Jacobians of f and g, which only differ to one another for the value of the corresponding eigenvalue λ i .This fact allows one to parameterize the problem, and to refer to a unique variational equation which is now an explicit function of the parameter ν = dλ i .The maximum Lyapunov exponent Λ of such a parametric equation can then be computed for increasing values of ν.The function Λ(ν) is the Master Stability Function (MSF), and its values assess the expansion (if positive) or contraction (if negative) rates in the directions of eigenvectors v i .One needs all of these rates to be negative in order for the synchronization manifold M to be stable.
As noticed for the first time in Chapter 5 of Ref. [4], and as illustrated in Fig. 1, all possible choices of (chaotic) flows and output functions, are in fact categorizable into only three classes of systems: • Class I systems, for which the MSF does not intercept the horizontal axis, like the yellow curve in Fig. 1.These systems intrinsically defy synchronization, because all directions of T are always (i.e., at all values of d) expanding, no matter which network architecture is used for connecting the nodes.Therefore, neither the synchronized solution x s (t) nor any other cluster-synchronized state will ever be stable.
• Class II systems, for which the MSF has a unique intercept with the horizontal axis at a critical value ν = ν * , like the violet line of Fig. 1.The scenario here is the opposite of that of Class I. Indeed, given any network G, the condition dλ 2 > ν * warrants stability of the synchronized solution.These systems, therefore, are always synchronizable, and the threshold for synchronization is d c ≡ ν * λ2 i.e., is inversely proportional to the second smallest eigenvalue of the Laplacian matrix.
• Class III systems, for which the MSF intercepts instead the horizontal axis at two critical points ν = ν * 1 and ν = ν * 2 , like the brown V-shaped curve of Fig. 1.In order for the synchronization solution to be stable, it is required in this case that the entire spectrum of eigenvalues of L falls (when multiplied by d) in between ν * 1 and ν * 2 .In other words, the two conditions dλ N < ν * 2 and dλ 2 > ν * 1 must be simultaneously verified, and this implies that not all networks succeed to synchronize Class III sys-tems.In fact, the former condition gives a bound d max = ν * 2 λ N for the coupling strength above which instabilities in tangent space start to occur in the direction of v N , the latter provides once again the threshold d c ≡ ν * 1 λ2 for complete synchronization to occur.
In general, one should point out that the region of stability could even be formed by the union of several intervals on the ν axis.However, the aim of present work is to describe the first transition to synchronization (i.e., the process that starts at d = 0 and occurs when progressively increasing the coupling strength) for which the three mentioned classes are indeed the only possible scenarios.When a system will show multiple regions of stability in the ν axis, instead, a series of synchronization and de-synchronization transitions (as many as the stability regions of the system) will occur, and this latter scenario will be reported by us elsewhere.
Finally, it is crucial to remark that all the above results are formally valid only for the whole network's synchronous solution.
The trajectories followed by the nodes' dynamics in each cluster-synchronous state slightly differ, instead, from those which are followed in the global solution, as they rigorously depend on the quotient network (and therefore on topology, node dynamics, and clusterization, see the Methods section for a detailed description of such differences), and several ad-hoc methods have been proposed to assess the stability of synchronization in each specific clustered state [19,22,[35][36][37][38].
However, since we are focusing on the transition to synchronization (i.e., on a regime where the coupling strength is normally very small), in the following we will adopt the approximation that such a difference will only be tight, and will not give rise to substantial changes in the calculation of maximum Lyapunov exponents.As a consequence, we will refer to the MSF calculated for the full network's synchronous solution as a an approximation of the values of the maximum Lyapunov exponents corresponding to the eigenvectors orthogonal to each cluster's synchronous solutions.Notice, furthermore, that at each value of the coupling strength one might have a distinct cluster-synchronous state.

The path to synchronization
With all this in mind, let us now move to describe all salient features characterizing the transition to the synchronization solution (as d increases from 0), and in particular to predict all the intermediate events that are taking place during such a transition.Since now, we anticipate that our results are valid for all systems in Class II, as well as for those in Class III (up to the maximum allowed value of the coupling strength i.e., for d < There are three conceptual steps that need to be made. The first step is that, as d progressively increases, the eigenvalues λ i cross the critical point (ν = ν * in Class II, or ν = ν * 1 in Class III) sequentially.The first condition which will be met will be, indeed, dλ N > ν * in Class II (dλ N > ν * 1 in Class III), while for larger values of d the other eigenvalues will cross the critical point one by one (if they are not degenerate) and in the reverse order of their size.
Therefore, one can use this very same order to progressively unfold the tangent space T .In particular, at any value of d, T can be factorized as T + (d) ⊗ T − (d), where T + (d) [T − (d)] is the subspace generated by the set of eigenvectors {v i } whose corresponding dλ i are below (above) the stability condition, i.e. for which one has 1 ) in Class III.In other words, the subspace T + (d) [T − (d)] contains only expanding (contracting) directions, and therefore the projection on it of the synchronization error δX will exponentially increase (shrink) in size.
The second step consists in taking note that, if one constructs the matrix V having as columns the eigenvectors then the rows of matrix (2) provide an orthonormal basis of R N as well, since V T V = I implies V T = V −1 , and hence also V V T = I.Therefore, one can now examine the eigenvectors component-wise and, for each eigenvector v i , define the following matrix Furthermore, following the same sequence which is progressively unfolding T , one recursively defines the following set of matrices S n It is worth discussing a few properties of such matrices.First of all, as v 1 is aligned with M, all its components are equal, and therefore E λ1 = 0 and S 1 = S 2 .Second, all the E-matrices, and thus all the S-matrices, are symmetric, non negative and have all diagonal entries equal to zero.In fact, the off diagonal ij elements of the matrix S n (n = 1, ..., N ) are nothing but the square of the norm of the vector obtained as the difference between the two vectors defined by rows i and j of matrix (2), truncated to their n last components.As so, the maximum value that any entry (ij) may have in the matrices S n is 2, which corresponds to the case in which such two vectors are orthonormal.In particular, all off-diagonal entries of S 2 = S 1 are equal to 2.
The third conceptual step consists in considering the fact that the Laplacian matrix L uniquely defines G, and as so any clustering property of the network G has to be reflected into a corresponding spectral feature of L [26,27].In this paper, we can prove rigorously that the synchronized clusters emerging during the transition of G can be associated to a localization of a group of the Laplacian eigenvectors on the clustered nodes.Namely, let us first define that a subset S ⊆ {v 2 , v 3 , ..., v N } consisting of k − 1 eigenvectors forms a spectral block localized at nodes {i 1 , ..., i k } if • each eigenvector belonging to S has all entries (except i 1 , i 2 , ..., i k ) equal to 0; • for each other eigenvector v not belonging to S, the entries i 1 , i 2 , ..., i k are all equal i.e., Moreover, all eigenvectors {v 2 , v 3 , ..., v N } are orthonormal to v 1 , and therefore the sum of all their entries must be equal to 0. The main theoretical result underpinning our study is the Theorem stated below.
Theorem.The following two statements are equivalent: 1.All k nodes belonging to a cluster defined by the indices {i 1 , . . ., i k } have the same connections with the same weights with all other nodes not belonging to the cluster, i.e. for any p, q ∈ {i 1 , . . ., i k } and j ̸ ∈ {i 1 , . . ., i k } one has L pj = L qj .
A group of nodes satisfying condition (1) of the theorem is also called an external equitable cell [39].
The reader interested in the mathematical proof of the theorem is referred to our Supplementary Information.We here concentrate, instead, on the main concepts involved.Conceptually, the first statement of the Theorem is tantamount to assert that the nodes belonging to a given cluster are indistinguishable to the eyes of any other node of the network, but puts no constraints on the way such nodes are connected among them within the cluster.Therefore, fulfillment of the statement is realized by (but not limited to) the case of a network's symmetry orbit.In other words, the first statement of the theorem says that the clustered nodes receive an equal input from the rest of the network, and therefore (for the principle that a same input will eventually -i.e. at sufficiently large coupling -imply a same output) they may synchronize independently on the synchronization properties of the rest of the graph.Therefore, the intermediate structured states emerging in the path to synchrony of a network are more general than the graph's symmetry orbits, but more specific than the graph's equitable partitions.
However, the most relevant consequence of the theorem is that the localization of the eigenvectors' components implies that the matrices S n may actually display entries equal to 2 also for n strictly larger than 2! Indeed, the (ij) entry of the matrices S n is just equal to Now, suppose that S is a spectral block localized at i, j and some other nodes.Then, if v k does not belong to S, the term (v kj − v ki ) 2 = 0, and one has therefore that the (ij) entry of S n has contributions only from those eigenvectors v n belonging to S. Now, all the times that a localized spectral block S is contained in the set of those eigenvectors generating T − (d) the corresponding cluster of nodes will emerge as stable synchronization cluster, because the tangent space of the corresponding synchronized solution (where synchrony is limited to those specific nodes) can be fully disentangled from the rest of T and will consist, moreover, of only contractive directions.

Complete description of the transition
By exploiting the outcomes of the three conceptual steps described above, an extremely simple (and computationally low demanding) technique can then be introduced, able to monitor and track localization of eigenvectors along the transition, and therefore to describe the path to synchronization.
The method consists in the following steps: • given a network G, one considers the Laplacian matrix L, and extracts its N eigenvalues λ i (ordered in size) and the corresponding eigenvectors {v i }; • one then calculates the N matrices E λi and S n (i = 1, ..., N ; n = 1, ..., N ); • one inspects the matrices S n in the same order with which the Laplacian's eigenvalues (when multiplied by d) crosses the critical point ν * (i.e., N, N −1, N − 2, ..., 2, 1), and looks for entries which are equal to 2; • when, for the first time in the sequence (say, for index p) an entry in matrix S p is (or multiple entries are) found equal to 2, a prediction is made that an event will occur in the transition: the cluster (or clusters) formed by the nodes with labels equal to those of the found entry (entries) will synchronize at the coupling strength value ν * /λ p .The inspection of matrices S n then continues, focusing only on the entries different from those already found to be 2 at level S p ; • after having inspected all S n matrices, one obtains the complete description of the sequence of events occurring in the transition, with the exact indication of all the values of the critical coupling strengths at which each of such events is occurring.
By events, we here mean either the formation of one (or many) new synchronized cluster(s), or the merging of different clusters into a single synchronized one.
Once again, we here remark that our results and methods are valid under the approximation that all clustersynchronous solutions are well described by the synchronization state that each node in the cluster would display in the complete synchronization scenario i.e., when the entire network synchronizes.In the following, we will demonstrate the accuracy of such an approximation with respect to synthetic and real-world networks.We begin with illustrating the method in a simple case, in order for the reader to have an immediate understanding of the consequences of the various steps that have been discussed so far.For this purpose, we consider the network sketched in panel (a) of Fig. 2, which consists of an all-to-all connected, symmetric, weighted graph of N = 10 nodes (see the Supplementary Information for the adjacency matrix of the network).By construction, the graph is endowed with three symmetric orbits: the first being composed by the pink nodes 1, 2 and 3, the second containing the blue nodes 4, 5 and 6, and the third being made of the four green nodes 7, 8, 9 and 10.The 10 eigenvalues of the Laplacian matrix, when ordered in size, are {0, 1, 1, 1, 4, 4, 4, 6, 6, 6}.
After calculations of the corresponding eigenvectors, the matrices E λi and S n are evaluated.Then, one starts inspecting matrices S n in the reverse order of the size of the corresponding eigenvalues.Panel (b) of Fig. 2 shows S 10 , which corresponds to λ 10 = 6, and it can be seen that there are no entries equal to 2 in such a matrix.Nor entries equal to 2 are found in S 9 (not shown).However, when inspecting S 8 (which corresponds to λ 8 = 6), one immediately identifies [panel (c) of Fig. 2] many entries equal to 2, which clearly define a cluster formed by nodes 7, 8, 9 and 10.A prediction is then made that the first event observed in the transition will be the synchronization of such nodes in a cluster, occurring at 2 Node 1 0 0 0 0.9 0.9 0.9 0.9 0.9 0.9 0.9 2 2 2 0.9 0.9 0.9 0.9 0 0 0 0.9 0.9 0.9 0.9 0.9 0 0 0 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0 2 2 2 2 0.9 0.9 0.9 2 0 2 2 2 0.9 0.9 0.9 2 2 0 2 2 0.9 0.9 0.9 0.9 0.9 2 2 2 2 0 0.9 2 2 2 0.9 0.9 0.9 2 2 2 0.9 0.9 0.Then, one continues inspecting the matrices S n and concentrates only on all the other entries.No further entry is found equal to 2 in S 7 , nor in S 6 .When scrutinizing S 5 [panel (d) of Fig. 2] one sees that other entries becomes equal to 2, and they indicate that the cluster formed by nodes 4, 5 and 6 will merge at d 2 = ν * /λ 5 = ν * /4 = 3 2 d 1 with the already existing group of synchronized nodes, forming this way a larger synchronized cluster.No further features are observed in S 4 and S 3 , while the analysis of S 2 [panel (e) of Fig. 2] reveals that at d 3 = ν * /λ 2 = ν * = 4d 2 also nodes 1, 2 and 3 join the existing cluster, determining the setting of the final synchronized state, where all nodes of the network evolve in unison.
Panel (f) of Fig. 2 schematically summarizes the pre-dicted sequence of events: if a dynamical system is networking with the architecture of panel (a) of Fig. 2, its path to synchrony will first (at d = d 1 ) see the formation of the synchronized cluster {7, 8, 9, 10}, then (at d = d 2 ) nodes {4, 5, 6} will join that cluster, and eventually (at d = d 3 ) all nodes will synchronize.Already at this stage it should be remarked that all predictions made are totally independent on f and g: changing the dynamical system f operating on each node and/or the output function g will result in the same sequence of events.The only difference will be that the estimated values d i at which the i th event will occur will be, under the adopted approximation, rescaled with the corresponding value of ν * .In order to show how factual is our prediction, we mon-  It is worth noticing that we here concentrate on chaotic systems, and do not consider instead periodic dynamics.The reason is that, in the MSF, ν = 0 corresponds to λ 1 = 0 i.e., to the eigenvector v 1 aligned with the syn-chronization manifold.Therefore, Λ(0) is the maximum Lyapunov exponent of the isolated system which is equal to 0 for a periodic dynamics, leading to the vanishing of the critical value ν * as well as of the threshold for complete synchronization.In turn, this would imply that intermediate clusters could not be observed as the entire network would synchronize for any infinitesimal coupling strength.
In order to properly quantify synchronization, one calculates the synchronization error over a given cluster, as where the sum runs over all nodes i forming the cluster, < ... > ∆T stands for a temporal average over a suitable time span ∆T , xcl is the average value of the vector x in the cluster, and N cl is the number of nodes in the cluster.
In addition, the synchronization error is normalized to its maximum value, so as to range from 1 to 0. The results are reported in Fig. 3, where the normalized synchronization errors are reported as a function of d (for both the Lorenz and the Rössler case) for the entire network (black dotted line) and for each of the three orbits in the network (yellow, orange and red lines).It is clearly seen that all predictions made are fully satisfied.Panels (b1-b4) of Fig. 3 report temporal snapshots of the dynamics of each of the 10 network's nodes, and illustrate visually the different collective phases which are observed during the transition to synchronization reported in Fig. 3b.
Looking at Fig. 3, it must be remarked that the transition to synchronization is identical for the Rössler and Lorenz systems, with the only difference being the different values ν * (0.179 for Rössler and 7.322 for Lorenz).In other words, the horizontal axis in Fig. 3b just corresponds to that of Fig. 3a multiplied by the ratio 0.179/7.322 of the two critical values.This depends on the fact that the set of clusters that are synchronizing during the transition, the nodes that compose them, and the order in which clusters synchronize are, under the approximation adopted in our study, completely independent on the specific dynamical system that is ruling the evolution of the network's nodes (as these features only depend instead on the spectrum of the Laplacian matrix, and consequently only on the topology of the graph).This is a strong result, because it implies that in all practical situations (i.e., when there are uncertainties in the model parameters, or even when the knowledge of the dynamics of the nodes is fully unavailable) the entire cluster sequence forming the transition to synchronization can still be predicted.

Synthetic networks of large size
The next step is to test the accuracy of the method in the case of large size graphs.To this purpose, we consider 2 networks that were synthetically generated in Ref. [42], where a specific method of generating graphs endowed with desired clusters was introduced that initially considers ensembles of disconnected nodes and then connects each one of the nodes in each of such ensembles to the same group of nodes of an external core network, this way forming clusters containing nodes of equal degree.While we refer the reader to Ref. [42] for the full details of the generating algorithm, we here limit ourselves to report the main attributes of the considered networks.The first network G 1 is of size N = 1, 000 nodes, and is endowed with two symmetry orbits that generate two distinct clusters of sizes 20 nodes (Cluster 1) and 10 nodes (Cluster 2), respectively.The second network G 2 is made of N = 10, 000 nodes and is endowed with four symmetry orbits generating four clusters of sizes that span more than an order of magnitude (Clusters 1 to 4 contain, respectively, 1,000, 300, 100, and 30 nodes).
Following the expectations which are detailed in the theorems of our Supplementary Information, the calculation of the Laplacian eigenvalues of G 1 allows one to identify a first group of 19 degenerate eigenvalues λ 277,278,...,295 = 4 and a second group of 9 degenerate eigenvalues λ 59,60,...,67 = 2.The analysis of the matrices S n then reveals that the first event in the transition is the synchronization of Cluster 1 at d 1 = ν * /4, followed by the synchronization of Cluster 2 at d 2 = ν * /2 = 2d 1 , this time in a state which is not synchronized with cluster 1, thus determining an overall state where two distinct synchronization clusters coexist.Eventually, the entire network synchronizes at d 3 = ν * /λ 2 = ν * /0.4758.
In the case of G 2 , four blocks of degenerate eigenvalues are indeed found in correspondence with the four clusters.The transition predicted by inspection of the matrices S n is characterized by the sequence of five events.First, the cluster with 30 nodes synchronizes at d 1 = ν * /5.Second, at d 2 = ν * /4, the cluster with 100 nodes synchronizes.At d 3 = ν * /3 (d 4 = ν * /2) also the cluster with 300 nodes (with 1,000 nodes) synchronizes.The four clusters evolve in four different synchronized states.Eventually, at d 5 = ν * /0.6025 the entire network synchronizes.
We then simulated the Rössler system on G 1 and G 2 and reported the results in Fig. 4, which are actually fully confirming the predicted scenarios.

Real-world and heterogeneous networks
We move now to show three applications to real-world networks.
For the first application, we have considered the network of the US power grid.The PowerGrid network consists of 4,941 nodes and 6,594 links, and it is publicly available at https://toreopsahl.com/datasets/#uspowergrid.
It was already the object of several studies in the literature, the first of which was the celebrated 1998 paper on small world networks [43].In the PowerGrid network, nodes are either generators, transformers or substations forming the power grid of the Western States of the United States of America, and therefore they have a specific geographical location.Recently, it was proven that a fraction of 16.7% of its nodes are forming non trivial clusters corresponding to symmetry orbits which are small in size, due to the geographical embedding of the graph [26].
Application of our method detects that the synchronization transition is made of a very well defined sequence of events, which involves the emergence of 381 clusters.The clusters that are being formed are all small in size, because of the constraints made by the geographical embedding.In particular, 310 clusters contain only 2 nodes, 49 clusters are made of 3 nodes, 14 clusters are formed by 4 nodes, 4 clusters have 5 nodes, 2 clusters appear with 6 nodes, 1 cluster has 7 nodes, and 1 cluster is made of 9 nodes.The overall number of network's nodes which get clustered during the transition is 871.A partial list of these clusters (spanning about one order of magnitude in the size of the corresponding eigenvalues) and the various values of coupling strengths at which the different events are predicted is available in Table 1 of the Supplementary Information.
We have then simulated the Rössler system on the PowerGrid network, and monitored the synchronization error on 6 specific clusters (highlighted in red in the list of Table 1 of the Supplementary Information) that our method foresees to emerge during the path to synchrony in a well established sequence and at well specific values of the coupling strength (d 1 = 0.04475, d 2 = 0.0596, d 3 = 0.0895, d 4 = 0.1294, d 5 = 0.179, d 6 = 0.3056) .The values d 1 , ..., d 6 are explicitly calculated in the Supplementary Information, and are marked as filled dots in the horizontal axis of Fig. 5 with the same colors identifying the corresponding clusters.Looking at Fig. 5 (a) one sees that, once again, the observed sequence of events matches the predicted one, with an excellent fit with the values d 1 , ..., d 6 , thus validating ex-post the approximation adopted in our study.
We also tested how robust is the predicted scenario against possible heterogeneities present in the network.To this purpose, one again simulates the Rössler system on the PowerGrid network, but this time one distributes randomly different values of the parameter b to different network's nodes.Precisely, for each node i of the network, one uses a parameter b i in the Rössler equations which is randomly sorted from a uniform distribution in the interval [0.1 − ϵ, 0.1 + ϵ], with the extra parameter ϵ that now quantifies the extent of heterogeneity in the graph.The results are reported in panel (b) of Fig. 5 for ϵ = 0.01, corresponding to 10 % of the value (b = 0.1) which was used for all nodes in the case of identical systems, thus representing a case of a rather large heterogeneity.
It has to be remarked that, when the networked units are not identical, the very same synchronization solution ceases to exist, and therefore it formally makes no sense to speak of the stability of a solution that does not even exist.Indeed, if one selects no matter which ensemble of nodes, the synchronization error never vanishes exactly in the ensemble.Nonetheless, it is still observed that the values of the normalized synchronization errors fluctuate around zero for some sets (clusters) of network's nodes which, therefore, anticipate the setting of the almost completely synchronized state (wherein all nodes evolve almost in unison).In Fig. 5(b) it is clearly seen that, while the synchronization errors approach zero at values that are obviously different from those predicted in the case of identical systems, the sequence at which the different clusters emerge during the transition is still preserved.Similar scenarios characterize the evolution of the network also when the heterogeneity is affecting the other two parameters (a and c) entering into the equations of the Rössler system.
The second (third) application is given with reference to a real-world biological (social) network.Precisely, we consider the Yeast protein-protein interaction network (a dataset made of N = 1, 647 nodes and E = 2, 518 edges [5]) and the ego-Facebook network (a dataset containing N = 2, 888 nodes and E = 2, 981 edges [6]).The results are reported in Fig. 6, and refer to ensemble averages over 100 different numerical simulations of identical Rössler systems (same parameters and initial conditions as in the case reported in Fig. 5) coupled by the structure of connections of the two graphs.For the biological network, it is found that the transition to synchronization includes 188 clusters, and Fig. 6a reports the synchronization error of 4 of them: C 1 and C 2 (which are both 2 nodes clusters), C 3 (a cluster containing 3 nodes), C 4 (a cluster of 7 nodes).For the social network, 18 clusters are found during the transition to synchronization, and also in this case Fig. 6b reports the synchronization error of 4 of them: C 1 and C 2 (which are both made of 2 nodes), C 3 (a cluster containing 5 nodes), C 4 (a subset of 3 nodes -out of a cluster of 280 nodes -which is forming a cluster by itself).In both cases, one easily sees that the numerically observed sequence of events perfectly matches the predicted one, with an excellent fit

DISCUSSION
In conclusion, our work provides an approximated, yet very accurate, description of the transition to synchronization in a network of identical systems.We unveil, indeed, that the path to synchronization is made of a sequence of events, each of which can be identified as either the nucleation of one (or several) cluster(s) of synchronized nodes, or to the merging of multiple synchronized clusters into a single one, or to the growth of an already existing synchronized cluster which enlarges its size.In our study, all nodes in a cluster have the same connections (and the same weights) with nodes not belonging to the cluster, and therefore they receive the same dynamical input from the rest of the network.
By combing methodologies borrowed from stability of nonlinear systems with tools of algebra and symmetries, and under the approximation that the systems' trajectories in all cluster-synchronous states do not substantially differ from those featured at complete synchronization, we have been able to introduce a simple and effective method able to provide the complete prediction of the sequence of such events, to identify which graph's node is belonging to each of the emergent clusters, and to give a forecast of the critical coupling strength values at which such events are taking place.
While the local node dynamics and the coupling function are entirely responsible for the specific synchronization class of the Master Stability Function, the sequence of events that are taking place during the transition to synchronization depends, instead, only on the graph's structure and, more precisely, on the only knowledge of the full spectrum of eigenvalues and eigenvectors of the Laplacian matrix.
Our method does not require an a-priori knowledge of the network's symmetries.We moreover clarify once forever the intimate nature of the clusters that are being formed along the transition path: they are formed by those nodes which are indistinguishable at the eyes of any other network's vertex.This implies that nodes in a synchronized cluster have the same connections (and the same weights) with nodes not belonging to the cluster, and therefore they receive the same dynamical input from the rest of the network.Synchronizable clusters in a network are therefore subsets more general than those defined by the graph's symmetry orbits, and at the same time more specific than those described by equitable partitions.
Finally, our work gave evidence of several extensive numerical simulations with both synthetic and real-world networks, and demonstrates how high is the accuracy of our predictions.Remarkably, the synchronization scenario in heterogeneous networks (i.e., networks made of non-identical units) preserves the predicted cluster sequence along the entire synchronization path.
Our results, therefore, call for a lot of applications of general interest in nonlinear science, ranging from synthesizing networks equipped with desired cluster(s) and modular behavior, until predicting the parallel (clustered) functioning of real-world networks from the analysis of their structure.

METHODS
In what follows we consider a connected weighted undirected graph G with N nodes and uniquely identified by its adjacency matrix A and its Laplacian matrix L. Furthermore, we will call λ 1 = 0 < λ 2 ≤ λ 3 . . .≤ λ N the (ordered in size) real and nonnegative eigenvalues of L, and The MSF for the network's synchronous solution Let us recall here that the equation governing a generic ensemble of N identical dynamical systems interplaying over a network G is Eq.( 1) i.e., ẋi = f (x i ) − d N j=1 L ij g(x j ), where x i (t) is a m-dimensional vector state describing the dynamics of each node i, f : R m −→ R m is the local (identical in all units) dynamical flow describing the evolution of the isolated systems, d is a real-valued coupling strength, L ij is the ij entry of the Laplacian matrix, and g : R m −→ R m is an output function describing the functional way through which units interplay.
L is a zero-row matrix, a property which, in turn, guarantees existence and invariance of the synchronized solution In order to study stability of such a solution, one considers perturbations δx i = x i − x s for i = 1, ..., N , and performs linear stability analysis of Eq. ( 1).The result are the following equations: where Jf (x s ) and Jg(x s ) are, respectively, the m × m Jacobian matrices of the flow and of the output function, both evolving in time following the synchronization solution's trajectory.One can, in fact, consider the global error δX ∈ R N m ≡ (δx 1 , δx 2 , ..., δx N ) T around the synchronous state, and recast Eqs.(4) as where I is the identity matrix, and ⊗ stands for the direct product.Moreover, L is a symmetric, zero row sum, matrix.As so, it is always diagonalizable, and the set of its eigenvectors forms an orthonormal basis of R N .The zero row sum property of L implies furthermore that λ 1 = 0 and that v 1 = 1 √ N (1, 1, ...., 1, 1) T .Therefore, all components of v 1 are equal, and this means that v 1 is aligned, in phase space, to the manifold M defined by the synchronization solution, and that an orthonormal basis for the space T tangent to M is just provided by the set of eigenvectors v 2 , v 3 , ..., v N .For the synchronization solution to be stable, the necessary condition is then that all directions of the tangent space be contractive.
One can now expand the error δX as a linear combination of the eigenvectors {v i } i.e., Then, plugging the expansion in Eq. ( 5) and operating the scalar product of both the right and left part of the equation times the eigenvectors v i , one obtains that the coefficients η i ∈ R m obey the equations ηi = [Jf (x s ) − dλ i Jg(x s )] η i .
Notice that the equations for the coefficient η i are variational, and only differ (at different i ′ s) for the eigenvalue λ i appearing in the evolution kernel.This entitles one to cleverly separate the structural and dynamical contributions, by introducing a parameter ν ≡ dλ, and by studying the m-dimensional parametric variational equation The kernel K(ν), indeed, only depends on f and g (i.e., on the dynamics), and the structure of the network is now encoded within a specific set of ν values (those obtained by multiplying d times the Laplacian's eigenvalues).
The maximum Lyapunov exponent Λ [i.e., the maximum of the m Lyapunov exponents of Eq. ( 6)] can then be computed for each value of ν.The function Λ(ν) is called the Master Stability Function, and only depends on f and g.At each value of d, a given network architecture is just mapped to a set of ν ̸ = 0 values.The corresponding values of Λ(ν) provide the expansion (if positive) or contraction (if negative) rates in the directions of the eigenvectors v 2 , v 3 , ..., v N , and therefore one needs all these values to be negative in order for M to be attractive in all directions of T .
Finally, notice that ν = 0 corresponds to λ 1 = 0 i.e., to the eigenvector v 1 aligned with M. Therefore, Λ(0) is equal to the maximum Lyapunov exponent of the isolated system ẋ = f (x).In turn, this implies that the Master Stability Function starts with a value which is strictly positive (strictly equal to 0) if the networks units are chaotic (periodic).
The 3 different Classes of systems supported by the Master Stability Function are illustrated in Fig. 1, and largely discussed in the Manuscript.

The stability properties of the clusters' synchronous solution
It is important to remark that all the above results are formally valid only for the whole network's synchronous solution.The trajectories followed by the nodes' dynamics in each cluster-synchronous state slightly differ, instead, from those which are followed in the global solution, as they rigorously depends on the quotient network, and therefore on topology, node dynamics, and clusterization.Let us indeed focus on a given cluster C l , and let us call C l the set of indices identifying the nodes that belong to C l .For each node i (i ∈ C l ), Eq. ( 1) becomes where the overall coupling is now split into the sum of an intra-cluster term and of a term accounting for the connections of the cluster to the rest of the network.Eq. ( 7) can be rewritten as where A is the adjacency matrix.This is because all the elements of the Laplacian matrix considered in the second coupling term are just the opposite of the corresponding terms of the adjacency matrix.The second coupling term is, indeed, limited to j / ∈ C l and therefore, by definition, it does not contain the diagonal element of the Laplacian (j = i) which is instead contained in the intra-cluster coupling term.
We now recall that the main theorem of our study asserts that synchronizable clusters are those formed by nodes which are equally connected to, and receive an equal input from, the rest of the network.Therefore, as the last term of Eq. ( 8) accounts for the total input received by node i from all nodes that do not belong to the cluster, our theorem ensures that it is a term which is formally independent on i.The cluster synchronous solution is x i (t) = x j (t) = x C l (t), ∀i ∈ C l and ∀j ∈ C l (j ̸ = i), and obeys the equation This is because the diagonal terms of the Laplacian matrix are equal to the sum of the number of intracluster connections and of the number of extra-cluster connections, the latter ones being now incorporated in the remaining coupling term, which once again does not depend on i.Once again, one can consider perturbations δx i = x i − x C l (∀i ∈ C l ), and perform linear stability analysis of Eq. ( 8).The result is where δX ∈ R N l m is the global error vector, and N l is the number of nodes forming the cluster C l .It is immediately seen that the linearized Eq. ( 10) is formally identical to Eq. ( 5), and therefore the same expansion of the error can be made with the eigenvectors of the Laplacian.The difference, however, is that the evaluation of the maximum Lyapunov exponents requires now to calculate the Jacobians of the functions f and g over the cluster-synchronous solution x C l , which obeys Eq. ( 9).In other words, while the MSF formalism calculates the Maximum Lyapunov exponents using the trajectories experienced by the whole network's synchronous solution (a state in which each node of the network repeats the same dynamical behavior of a single, isolated, system), the trajectories experienced by the cluster synchronous state are perturbed by an extra term and depend therefore explicitly on the entire network's dynamics, and on the specific coupling function.This fact leads to two main consequences: • The very same cluster-synchronous solution is not invariant, as K(t) explicitly depends on d.In particular, at each value of the coupling strength one would have a distinct cluster-synchronous state, and therefore the entire framework of the MSF would make no sense in this case as it would not be possible to rigorously separate dynamics and structure; • The perturbation K(t) may lead the synchronous trajectories to visit areas of the phase space which are instead never visited by a single isolated system, and therefore it may determine slight variations in the calculations of the maximum Lyapunov exponents, and consequently slight variations in the determination of the critical coupling strength value at which the cluster synchronous state becomes stable.
On the other hand, the term K(t) is directly proportional to d, and therefore it has to be expected that such a perturbation will in fact be small across the transition to complete synchronization, where d starts from 0 and only slightly increases to values which are normally smaller than 1.In addition, K(t) consists of the sum of terms which are in general uncorrelated, as there are no constraints on the dynamics of the nodes which do not form part of the synchronous cluster.This latter fact would also contribute to determine smallness of the perturbation.
In our study, we decided therefore to adopt a practical approximation to the problem, by assuming that the perturbation K(t) is always negligible and consequently that the trajectories visited by the clustered synchronous nodes are always equal to those that characterize complete synchronization (i.e., those of a single, isolated, system).This allows one to refer to a unique MSF (the one constructed in the whole network's synchronous state) for determining the stability properties of all cluster synchronized states.

Simulations
Simulations were performed with an adaptative Tsit integration algorithm implemented in Julia.In each trial, the network is simulated for a total period of 1,500 time units, and the synchronization errors are averaged over the last ∆T = 100 time units.
As one is only interested to monitor the vanishing of E cl , with the purpose of saving calculations in all our simulations the synchronization error is computed by only taking into account those variables of the system's state where the coupling is acting.This implies that, when referring to the Rössler (Lorenz) system, E cl has been evaluated taking into account only the difference y i − ȳcl (x i − xcl ).The results are, indeed, identical to those obtained when all state variables of the systems are accounted for in the evaluation of E cl , as its formal definition of Eq. ( 3) would instead require.
By hypothesis, we have that L ji = L ij = a for all i ∈ {1, . . ., k}, and hence It means that every entry of Lu after the k-th entry is 0. The second part (the sum of all the entries of the vector Lu is zero) holds for any vector: the column (or row) sum of L is zero, that is, 1 N L = 0 N , which implies 1 N Lu = 0 N for any vector u, where 1 N (respectively 0 N ) is a N -dimensional vector made of all entries equal to 1 (respectively 0).Consider now L ′ , the k × k principal submatrix of L given by the first k rows and columns.If u ∈ U , then u = (u ′ |0 N −k ) T and we can write, in matrix block form, where Bu ′ = 0 since Lu ∈ U .In particular, Lu = λu if and only if L ′ u ′ = λu ′ .Consider now the graph G ′ induced by the vertices 1, . . ., k.Its Laplacian L(G ′ ) equals L ′ except for the fact that it has diagonal elements Let u ′ 1 , . . ., u ′ k be an orthonormal basis of the Laplacian L(G ′ ) with u ′ 1 a constant vector.In particular, the sum of the entries of each u ′ 2 , . . ., u ′ k must be zero.If we define u i = (u ′ i |0 N −k ) T for all i = 2, . . ., k, then such vectors are in the subspace U and, following the discussion above, u 2 , . . ., u k must be eigenvectors of L.Moreover, the eigenvalues corresponding to u 2 , . . ., u k depend only on the structure of the induced graph G ′ and the external degree d.
Next, we show that every vector v ∈ R N orthogonal to u 2 , . . ., u k must be constant on its first k entries.Indeed, since u ′ 2 , . . ., u ′ k are linearly independent and U has dimension k −1, these vectors are also generate U .In particular, v must be orthogonal to all vectors in U .Recall (see above) that U has a basis {e 2 , . . ., e k }, and note that v • e j = v 1 − v j = 0 implies v 1 = v j for all k = 2, . . ., k.
2 =⇒ 1 Without loss of generality, we assume that the last k − 1 eigenvectors form the spectral block localized at nodes {1, . . .k}.Then, the matrix V having the eigenvectors as columns can be written in block form as where V ′ and V ′′ are k × (k − 1) and (N − k) × (N − k + 1) matrices respectively, C is a k × (N − k + 1) column-constant matrix (C pj = C qj for all p, q, j), and If Σ is the diagonal matrix of the eigenvalues λ 1 , λ 2 , . . ., λ N , we can write the Laplacian as where D 1 (respectively D 2 ) is the diagonal matrix of the eigenvalues λ 1 , . . ., λ N −k+1 (respectively λ N −k+1 , . . ., λ N ), and 0 represents a zero matrix of the appropriate size.This implies that the connections between the nodes {1, . . ., k} and the nodes {k + 1, . . ., N } are described by the top right submatrix Finally, note that if C is a column-constant matrix, so is CM for any matrix M where the product exists.Or, explicitly, if 1 ≤ p, q ≤ k and j > k, then For the final part of the statement, let us remind the (i, j) entry of S n is equal to Let i, j be any two different nodes at which the spectral block S is localized.Then, if v k does not belong to S the term (v kj − v ki ) 2 is equal to 0. Hence, the (i, j) entry of S n changes only when v n belongs to S. So, if λ n is greater than the maximum eigenvalue of S then the (i, j) entry of S n is 0. On the other hand, if λ n is less than the minimal eigenvalue of S the (i, j) entry of S n is 2, as was to be shown.
Actually, to be mathematically correct, the formulation of the second statement of the theorem should be There is an eigenbasis v 1 , . . ., v N and a spectral block S localized at nodes i 1 , . . ., i k , since it could happen that eigenvectors not from U could have the same eigenvalues as eigenvectors in S, so they could be not orthogonal to U .However, such a possibility does not affect the synchronization scenario, so we decided to omit to mention explicitly this extra case for the sake of simplicity.
One important example of equitable cells are symmetry orbits, that is, orbits under the action of the automorphism group of the graph.Note that, in real-world networks, most orbits are either complete or empty subgraphs with all permutations of the vertices realized as network symmetries [1,4].This particular case can be characterized as follows.
Theorem .3.Let G be a connected network with N vertices and Laplacian L and let {i 1 , . . ., i k } ⊆ {1, . . ., N }.Then the following three statements are equivalent: 1.For any pair of vertices indexed in the set {i 1 , . . ., i k } a permutation of this pair preserves the Laplacian of the network.
3. There is a spectral block S localized at nodes i 1 , . . ., i k and all k − 1 eigenvectors of S have the same, degenerate, eigenvalue.
In the proof of this theorem we will use the following proposition: Proposition .4.Let σ be a permutation of nodes {1, . . ., N } with permutation matrix P .Then the following two statements are equivalent 1. Permutation σ preserves L, i.e.P −1 LP = L.
2. For any eigenvector v with eigenvalue λ, the vector P v is also eigenvector of L with the same eigenvalue λ.
Proof of Proposition .4. 1 =⇒ 2 Let us consider any eigenvector v of L with an eigenvalue λ.Recall that a permutation matrix is invertible (in fact, orthogonal).Then 2 =⇒ 1 Let λ 1 , . . ., λ N be the eigenvalues of L with corresponding orthonormal eigenvectors v 1 , . . ., v N .Consider the associated spectral decomposition of L, that is, where Σ is a diagonal matrix of the eigenvalues and V the matrix with columns the corresponding eigenvectors.By hypothesis, {P v 1 , P v 2 , . . ., P v N } are also eigenvectors with the same eigenvalues λ 1 , λ 2 , . . ., λ N , and they are also orthonormal, since P is a permutation, hence orthogonal, matrix.As the matrix having P v 1 , P v 2 , . . ., P v N as columns is precisely P V , we can write and we are done.
Proof of Theorem .coupling parameter is high enough (indeed, when its critical value is reached), while being a SEC relative to clusters C 1 , . . ., C m guarantees cluster synchronization when C 1 to C m are synchronized and the coupling parameter is high enough.(Note that, for simplicity, this discussion only considers Type II systems, see Fig. 1 of the main text.) We can now clarify the relation between SEC and (external) equitable partition.If C 1 , . . ., C r is an equitable partition with r cells, each cell C i is a SEC relative to the other cells.This guarantees synchronization when the critical values of all the cells are reached, but not necessarily before.In particular, it doesn't explain the order at which the cluster synchronization occurs, nor partial synchronization of, nor merger between, subsets of the equitable cells, as detected by the S n algorithm presented in the main text.
In the theoretical and numerical analysis of the present article, all the cluster synchronization examples found correspond to SECs, or in SECs relative to synchronized clusters, confirming our intuition (see Supplementary Figs. 8, 9).We conjecture that cluster synchronization can only occur in those cases.

THE SYNCHRONIZATION CLUSTERS DESCRIBED IN THE MAIN TEXT
This final Section contains the details of the various networks considered in the main text.First, we report here below the adjacency matrix A of the network depicted in Figure 2a of the main text, which is Second, we report the output of the application of the method described in the main text to the PowerGrid network.In the first column of the list, we report the values of λ at which an event of cluster formation occurs along the transition, and the corresponding values of 1/λ which (once multiplied by ν * ) give the critical coupling strength's values at which such event has to be observed.Each event consists of the simultaneous emergence of synchronization clusters, and the nodes forming each one of them are reported in the second column of the list.
The 6 clusters that are highlighted in red are those whose synchronization error is reported in Fig. 5 of the main text, for a Rössler chaotic system with parameters and coupling function such that it belongs to Class II, displaying (a) A toy network with 8 nodes and 11 links, reproduced from [3].(b) Our Sn algorithm correctly predicts cluster synchronization on clusters {7, 8}, {1, 2} and {5, 6}, and {1, 2, 3}, in that order.Note that all clusters are SECs, and all except {1, 2, 3} are also SOs.For SECs, the critical eigenvalue λcrit (shown as λ + d, where λ is the second smallest Laplacian eigenvalue, and d the external degree, of the induced subgraph) of each cluster, shown, determines the order at which cluster synchronization occurs (largest to smallest).The Laplacian eigenvalues of this toy network, λ1 ≤ . . .≤ λ8, are 0, 0.63, 1, 3, 3, 4, 4, 6.37 (up to 2 decimal places) and, as predicted in Theorem .2,λcrit coincides with λn at the point when the corresponding 2-block in the matrix Sn appears, as shown on the table above.Complete synchronization corresponds to the second smallest Laplacian eigenvalue of the whole network, λ2 = 0.62 (2 decimal places).(c) Synchronization error for each of the predicted clusters (color-code in the legend).Simulations are performed with identical Lorenz chaotic oscillators, and with the same parameters and initial conditions used for generating Fig. 3a of the main text.showing small examples of typical real-world symmetric motifs (subgraphs where the network symmetry is locally generated [1,4]), shown inside the dotted lines.Each symmetric motif is composed of one or more vertex orbits (SOs in our terminology), shown by color.Of the 10 SOs in this toy example, 4 are SECs (dark green, brown, orange, and turquoise), corresponding to symmetric motifs with one orbit, and 6 are relative SECs: the red orbit is a SEC relative to the dark blue orbit (and vice-versa), the purple orbit relative to the green orbit (and vice-versa), and the yellow orbit relative to the pale blue orbit (and vice-versa).Notably, the yellow orbit contains two SECs, namely the clusters {30, 31} and {32, 33}, separately.The Laplacian eigenvalues of this toy network, λ1 ≤ . . .≤ λ33, are 0.0, 0.04, 0.11, 0.14, 0.23, 0.27, 0.38, 0.42, 0.59, 1.0, 1.0, 1.0, 1.0, 1.0, 1.37, 1.7, 2.0, 2.09, 2.49, 2.62, 2.72, 3.0, 3.41, 3.73, 3.97, 4.0, 4.0, 4.42, 5.09, 5.25, 5.53, 6.11, and 7.3 (up to 2 decimal places).(b) Our Sn algorithm correctly predicts cluster synchronization on the clusters shown, in the order shown.Note that the four SOs that are SECs become synchronized at their critical value λcrit (shown as λ + d, as in Supplementary Fig. 8), plus {30, 31} and {32, 33}, which are SECs but not SOs.As expected, the critical value λcrit of each SEC corresponds to the nth smallest Laplacian eigenvalue of the whole network, λn, with λ2 = 0.04 (2 decimal places) achieving global synchronization.The three relative SECs also become synchronized after, or at the same time as (depending on the relative values of the critical values) the clusters they are SECs relative to.Indeed, the red cluster becomes synchronized at the same time as the dark blue cluster, and the same for the purple and green clusters.
For the yellow orbit, note it contains two (structurally identical) SECs ({30, 31}, and {32, 33}), which synchronize separately when their critical eigenvalue is reached, but their union, which is not a SEC, rather a SEC with relative to the light blue orbit, only synchronizes when the light blue orbit does.The predicted sequence of events is fully confirmed in our simulations.
The list reported here below is limited to the first 11 predicted events.

Figure 1 .
Figure 1.Classes in the Master Stability Function.The maximum Lyapunov exponent Λ as a function of the parameter ν (see text for definitions), for a chaotic flow (Λ(0) > 0).The curve Λ(ν) is called the Master Stability Function (MSF).Given any pair of f and g, only three classes of systems are possible: Class I systems (yellow line) for which the MSF does not intercept the horizontal axis; Class II systems (violet line), for which the MSF has a unique intercept with the horizontal axis at ν = ν * ; Class III systems (brown line), for which the MSF intercepts the horizontal axis at two critical points ν = ν * 1 and ν = ν * 2 .

Figure 2 .
Figure 2. Predicting the transition to synchronization.(a) An all-to-all connected, symmetric, weighted graph of N = 10 nodes is considered.The graph is endowed with three symmetry orbits: the one composed by the red nodes {1, 2, 3}, the one made of the orange nodes {4, 5, 6}, and the one made of the four yellow nodes {7, 8, 9, 10}.In the sketch, the widths of the links are proportional to the corresponding weights, and the sizes of the nodes are proportional to the corresponding strengths.(b) The entries of S10, which corresponds to λ10 = 6.(c) S8 (associated to λ8 = 6), where the entries equal to 2 clearly define a cluster formed by nodes {7, 8, 9, 10}.The first predicted event in the transition then consists in the synchronization of such nodes at d1 = ν * /λ8 = ν * /6.(d) S5 (related to λ5 = 4), where additional entries become equal to 2, indicating a second foreseen event in which nodes {4, 5, 6} join the existent synchronization cluster at d2 = ν * /λ5 = ν * /4.(e) S2 (corresponding to λ2 = 1) where it is seen that nodes {1, 2, 3} also join the existing synchronized cluster at d3 = ν * /λ2 = ν * in a third predicted event where complete synchronization of the network takes place.(f) The expected events (and their exact sequence) occurring in the path to synchrony of the network's architecture depicted in panel (a).The bar at the bottom of the Figure gives the color code used in panels (b-e) for matrices' entries.

Figure 3 .
Figure 3.The numerical verification of the predicted transition.Panels (a,b): The normalized synchronization errors E cl (see text for definition) as a function of d, for the Lorenz [panel (a)] and the Rössler [panel (b)] case.Data refers to ensemble averages over 500 different numerical simulations of the network sketched in panel (a) of Fig. 2. Cluster 1 (yellow line) is formed by nodes {7, 8, 9, 10}, Cluster 2 (orange line) is formed by nodes {4, 5, 6}, and Cluster 3 (red line) is formed by nodes {1, 2, 3}.The black dotted line refers to the synchronization error of the entire network(EN).In both panels it is seen that the expected sequence of events taking place during the transition is verified.Furthermore, the values d1 = 7.322/6 = 1.220 (d1 = 0.179/6 = 0.0298), d2 = 7.322/4 = 1.8305 (d2 = 0.179/4 = 0.04475) and d3 = 7.322 (d3 = 0.179) are marked in the horizontal axis respectively with a yellow, orange and red filled dot in panel (a) (in panel (b)), indicating how accurate are the predictions and approximations made on the corresponding critical values for the coupling strength.For each interval, the arrow points to the composition of the synchronized cluster that is being observed, once again in perfect harmony with the predictions made.Finally, we have verified that no extra synchronization features emerge during the transition, other than those explicitly foreseen in Fig.2.Panels (b1-b4): Temporal snapshots illustrating the evolution of the y variable of each of the 10 network's nodes (see color code at the bottom of the four panels) during the transition to synchronization reported in panel (b).At d = 0.01 (panel b1) the nodes display a fully uncorrelated dynamics.At d = 0.035 (panel b2) the yellow nodes(7,8,9,10) are clustered and display a synchronous motion, whereas all other nodes feature a uncorrelated dynamics.At d = 0.1 (panel b3) the violet nodes(4,5,6)  have joined the clustered evolution, while nodes (1,2,3) remains unsynchronized.Finally, at d = 0.2 (panel b4), all network's nodes are synchronized.

c 1 c 2 c 3 c 4 EN c 1 c 2 EN d 1 d 2 d 3 d 1 d 2 d 3 d 4 d 5 Figure 4 .
Figure 4. Applications to large size synthetic networks.E cl (see text for definition) vs. d, for the Rössler system (see the differential equations in the text).Data in panel (a) [in panel (b)] refer to ensemble averages over 50 (150) different numerical simulations of the graph G1 (G2) described in the main text.In both panels, the legend sets the color code for the curves corresponding to each of the existing clusters Ci and to the Entire Network (EN).Once again, the two predicted transitions are verified.(a) Cluster 1 synchronizes at d1 = 0.179/4 = 0.04475 (marked with a yellow dot), Cluster 2 synchronizes at d2 = 0.179/2 = 0.0895 (orange dot), and the entire network synchronizes at d3 = 0.179/0.4758= 0.376 (black dot), as predicted.(b) The cluster with 30 nodes synchronizes at d1 = 0.179/5 = 0.0358 (yellow dot), the cluster with 100 nodes at d2 = 0.179/4 = 0.04475 (orange dot), the cluster with 300 nodes at d3 = 0.179/3 = 0.0597 (red dot), the cluster with 1,000 nodes at d4 = 0.179/2 = 0.0895 (violet dot), and the entire network at d5 = 0.179/0.6025= 0.297 (black dot).

6 c 1 c 2 c 3 c 4 c 5 c 6 EN c 1 c 2 c 3 c 4 c 5 c 6 5 Figure 5 .
Figure 5. Applications to the PowerGrid network.E cl (see text for definition) vs. d, for the Rössler system (see the differential equations in the text).Data in panel (a) [in panel (b)] refer to ensemble averages over 850 (200) different numerical simulations of the PowerGrid network.As in Fig.4, the legends of both panels set the color code for the curves corresponding to each of the reported clusters Ci and to the Entire Network (EN).Panel (a) reports the case of identical Rössler systems, and the error of 6 specific clusters is plotted (see Table 1 of the Supplementary Information for the composition of each of the 6 clusters Ci).The observed sequence of events perfectly matches the predicted one, with an excellent fit with the values d1, ..., d6.In panel (b) the effects of heterogeneity in the network are reported.Namely, for each node i of the PowerGrid network, the parameter bi in the Rössler equations is randomly sorted from a uniform distribution in the interval [0.1 − ϵ, 0.1 + ϵ].The curves plotted refer to ϵ = 0.01.

1 Figure 6 .
Figure 6.Application to biological and social networks.Synchronization error for the considered clusters C1, C2, C3 and C4 (see the color-code at the top of the panel) for the Yeast protein-protein interaction network[5]  (panel a), and the ego-Facebook network[6]  (panel b).In both panels, we report also the synchronization error of the entire network (EN, black dotted line).
Figure 8.(a) A toy network with 8 nodes and 11 links, reproduced from[3].(b) Our Sn algorithm correctly predicts cluster synchronization on clusters {7, 8}, {1, 2} and {5, 6}, and {1, 2, 3}, in that order.Note that all clusters are SECs, and all except {1, 2, 3} are also SOs.For SECs, the critical eigenvalue λcrit (shown as λ + d, where λ is the second smallest Laplacian eigenvalue, and d the external degree, of the induced subgraph) of each cluster, shown, determines the order at which cluster synchronization occurs (largest to smallest).The Laplacian eigenvalues of this toy network, λ1 ≤ . . .≤ λ8, are 0, 0.63, 1, 3, 3, 4, 4, 6.37 (up to 2 decimal places) and, as predicted in Theorem .2,λcrit coincides with λn at the point when the corresponding 2-block in the matrix Sn appears, as shown on the table above.Complete synchronization corresponds to the second smallest Laplacian eigenvalue of the whole network, λ2 = 0.62 (2 decimal places).(c) Synchronization error for each of the predicted clusters (color-code in the legend).Simulations are performed with identical Lorenz chaotic oscillators, and with the same parameters and initial conditions used for generating Fig.3aof the main text.

Figure 9 .
Figure 9. (a) A toy network with 33 nodes and 39 links, reproduced from[4], showing small examples of typical real-world symmetric motifs (subgraphs where the network symmetry is locally generated[1, 4]), shown inside the dotted lines.Each symmetric motif is composed of one or more vertex orbits (SOs in our terminology), shown by color.Of the 10 SOs in this toy example, 4 are SECs (dark green, brown, orange, and turquoise), corresponding to symmetric motifs with one orbit, and 6 are relative SECs: the red orbit is a SEC relative to the dark blue orbit (and vice-versa), the purple orbit relative to the green orbit (and vice-versa), and the yellow orbit relative to the pale blue orbit (and vice-versa).Notably, the yellow orbit contains two SECs, namely the clusters {30, 31} and {32, 33}, separately.The Laplacian eigenvalues of this toy network, λ1 ≤ . . .≤ λ33, are 0.0, 0.04, 0.11, 0.14, 0.23, 0.27, 0.38, 0.42, 0.59, 1.0, 1.0, 1.0, 1.0, 1.0, 1.37, 1.7, 2.0, 2.09, 2.49, 2.62, 2.72, 3.0, 3.41, 3.73, 3.97, 4.0, 4.0, 4.42, 5.09, 5.25, 5.53, 6.11, and 7.3 (up to 2 decimal places).(b) Our Sn algorithm correctly predicts cluster synchronization on the clusters shown, in the order shown.Note that the four SOs that are SECs become synchronized at their critical value λcrit (shown as λ + d, as in Supplementary Fig.8), plus {30, 31} and {32, 33}, which are SECs but not SOs.As expected, the critical value λcrit of each SEC corresponds to the nth smallest Laplacian eigenvalue of the whole network, λn, with λ2 = 0.04 (2 decimal places) achieving global synchronization.The three relative SECs also become synchronized after, or at the same time as (depending on the relative values of the critical values) the clusters they are SECs relative to.Indeed, the red cluster becomes synchronized at the same time as the dark blue cluster, and the same for the purple and green clusters.For the yellow orbit, note it contains two (structurally identical) SECs ({30, 31}, and {32, 33}), which synchronize separately when their critical eigenvalue is reached, but their union, which is not a SEC, rather a SEC with relative to the light blue orbit, only synchronizes when the light blue orbit does.The predicted sequence of events is fully confirmed in our simulations.