## Abstract

Synchronization is an important and prevalent phenomenon in natural and engineered systems. In many dynamical networks, the coupling is balanced or adjusted to admit global synchronization, a condition called Laplacian coupling. Many networks exhibit incomplete synchronization, where two or more clusters of synchronization persist, and computational group theory has recently proved to be valuable in discovering these cluster states based on the topology of the network. In the important case of Laplacian coupling, additional synchronization patterns can exist that would not be predicted from the group theory analysis alone. Understanding how and when clusters form, merge, and persist is essential for understanding collective dynamics, synchronization, and failure mechanisms of complex networks such as electric power grids, distributed control networks, and autonomous swarming vehicles. We describe a method to find and analyze all of the possible cluster synchronization patterns in a Laplacian-coupled network, by applying methods of computational group theory to dynamically equivalent networks. We present a general technique to evaluate the stability of each of the dynamically valid cluster synchronization patterns. Our results are validated in an optoelectronic experiment on a five-node network that confirms the synchronization patterns predicted by the theory.

- synchronization
- complex networks
- stability
- computational group theory

## INTRODUCTION

Synchronization of oscillators in large networks has been an interesting problem for many years. It is a phenomenon that shows up in many natural and man-made systems (*1*–*3*). Global synchronization is a particular type of synchronization in which the oscillators all follow the same trajectory in state space. A network where this is desirable would be generators in a power grid, as well as some control networks and swarming autonomous vehicles. Although global synchronization has a well-developed theory (*4*–*6*), a more recently studied, more complex phenomenon we will call cluster synchronization (CS) has attracted considerable attention (*7*–*16*). In CS, the network evolves into subsets of oscillators in which members of the same cluster are synchronized to the same trajectory, but members of different clusters are not. Such synchronized clusters may show up in swarms of animals, where the network is the simple visual link to one’s neighbors, or swarms of unmanned autonomous vehicles that are connected by a local communication network. Clusters may also show up in power grids, where they would be a sign of a problem, that is, loss of global synchronization.

Given the increase in man-made networks and the growing use of network theory to describe natural systems (for example, food webs, neuronal and genetic networks), it is important to develop a basic approach to determine what cluster structures are possible in a given network. Here, we show what methods can be employed to do this using the concept of oscillators coupled through connection to other (not necessarily all) oscillators in the network. We extend these methods to the currently unsolved problem of finding clusters in networks of oscillators that have a self-coupling to balance incoming signals from other oscillators, often called Laplacian coupling (more on this below). This allows synchronization clusters to be found that elude other methods of finding cluster patterns. We show how to use and extend symmetry methods to find all possible clusters in such networks of Laplacian-coupled oscillators. First, we show how network symmetries can lead to CS. Then, we show how we can go beyond this to analyze CS resulting from Laplacian coupling.

In the Symmetries and Clusters in Networks section, we review the concepts of symmetries and clusters in networks of coupled oscillators. In the Analyzing CS Patterns section, we discuss methods to uncover all of the possible CS patterns in a given network. Our main results are contained in the Stability Analysis and Experimental Validation section, where we present a stability analysis that applies to any CS pattern. We also present an elecro-optic experiment that confirms the patterns of synchronization predicted by the theory before presenting our conclusions.

## SYMMETRIES AND CLUSTERS IN NETWORKS

Figure 1A shows a four-oscillator or four-node network, where the oscillators are identical, as are the couplings between them, which are bidirectional, meaning that signals flow in both directions to the oscillators by the same amount and influence on the connected oscillators. This network has a total of six symmetries. We show two. In Fig. 1B, we show the result of a reflection [shown in Fig. 1A] interchanging nodes 1 and 2. The structure of the network remains the same. In Fig. 1C, we show a rotation of the network by 120°, which also leaves the network indistinguishable from the original in Fig. 1A. These symmetries are manifest in the symmetries of a matrix that describes the network, the adjacency matrix. This matrix, denoted by *A*, is set so that *A*_{ij} = 1 if *i* and *j* nodes are connected and *A*_{ij} = 0 otherwise. For the network in Fig. 1A, the adjacency matrix *A* is(1)

The adjacency matrix plays a crucial role in modeling the dynamics of many networks of identical nodes because it provides the coupling between the oscillators at each node. We will write the dynamics of the networks as follows(2)*i* = 1,…,*N*, where **x**_{i} is the vector of dynamical variables for the *i*th oscillator, is the time derivative of the *i*th node’s variables, *N* is the number of oscillators (nodes), **F** is the vector field of each node (governing each oscillator’s isolated dynamics), and **H** is a coupling function for each link in the network of each oscillator to another. Several papers (*8*, *11*, *12*, *15*, *17*–*20*) have used Eq. 2 to model the dynamics of a network.

Network symmetry applied to the adjacency matrix leaves it unchanged. We recall that symmetries of an object (a network in this case) form a mathematical group 𝒢. Any element in this group, for example *g*, is represented in the space of network nodes as a permutation matrix, for example *R*_{g}. The invariance of *A* under the action of the symmetry immediately implies *R*_{g}*A* = *AR*_{g}; *A* commutes with all group actions, and the equations of motion for nodes that are mapped into each other are the same. For example, in the network in Fig. 1A, nodes 1, 2, and 3 have the same equations of motion, so if they are started from the same initial conditions, they will remain synchronized indefinitely. Node 4 cannot be permuted into any of the other nodes, and it will not synchronize with the others; hence, we color it differently in Fig. 1B and Fig. 1C. This is the intimate relationship between symmetry and dynamics in networks. We see that it immediately separates the network into two clusters: 1,2,3 and 4.

Many of the recent studies of CS are on particular networks that are known or engineered to exhibit cluster synchrony. However, as we have shown above, group theory provides a general approach to finding clusters in arbitrary networks. Steps toward more general approaches using group theory to analyze CS began with the work of Golubitsky and Stewart (*21*) and Golubitsky *et al*. (*22*), where network symmetries are known and can be shown to support CS. Recently, computational methods have been used to study simple symmetric networks, and an approach has been developed that relates the symmetries with the emergence of the CS states (*23*). Finally, we showed that such approaches can be applied to more complex networks with hundreds of oscillators using computational group theory (*24*).

An alternative description of the network dynamics is the following in the case of Laplacian coupling(3)where the coupling from oscillator *j* to oscillator *i* is given by the difference between the output functions **H**(**x**_{j}) and **H**(**x**_{i}). Several papers (*2*, *5*, *7*, *10*, *20*, *23*, *25*) have used Eq. 3 to model the dynamics of a network. Equation 3 can be rewritten as follows(4)which has the same structure as Eq. 2, but now the adjacency matrix has been replaced by the Laplacian matrix, *L* = {*L*_{ij}}, *L*_{ij} = *A*_{ij} *−* δ_{ij} Σ_{j} *A*_{ij} , where δ_{ij} is the Kronecker delta. By construction, then, the sums of the rows of the matrix *L* are equal to zero, that is, the inputs to the *i*th node are balanced by the diagonal self-coupling.

Golubitsky and colleagues (*26*) have shown in their work that, for networks that have balanced coupling (all nodes receive the same cumulative input weights, accounting for adjacent nodes and self-coupling, an example of which is Eq. 4), CS can emerge in many patterns that are not directly the result of symmetries. An example of this is global synchronization, which is not a result of symmetries in the network.

In general, the patterns of cluster synchrony that can be observed in a network are not unique (*27*, *28*); hence, an important problem is that of determining the parameter ranges for stability and multistability for the observed patterns. Although it is known that the stability of the global synchronization state for an arbitrary network can be characterized by using the master stability function formalism (*5*), a corresponding analysis that applies to the CS patterns that may emerge in a network is not available. Here, we address this problem by providing necessary and sufficient conditions for stability of each CS pattern under very general assumptions. Our analysis applies to systems for which the functions describing the individual dynamics (possibly chaotic) and the interactions between the systems are arbitrary, whether the network is described by an adjacency matrix (Eq. 2) or by a Laplacian matrix (Eq. 4), both of which are used to model network interactions.

In what follows, we show how to find all of the CS patterns that may emerge in a given network topology, described either by an adjacency matrix or by a Laplacian matrix and how to evaluate the stability of each allowed pattern. We demonstrate all the above phenomena in an optoelectronic experiment on a five-node network that displays all of the possible CS patterns predicted by the theory.

## ANALYZING CS PATTERNS

Here, we attempt to address the following problem: Given a network structure (either in terms of an adjacency matrix in Eq. 2 or of a Laplacian matrix in Eq. 4), can we find all of the CS patterns that are allowed? For simplicity, we will proceed under the assumption that the network dynamics are described by Eq. 4, but all of our results include the simpler case that the dynamics are described by Eq. 2. Indeed, as we will see, the case of the Laplacian matrix is, in general, more complex to deal with than that of the adjacency matrix. So, we consider the most difficult case.

First, we note that a symmetry of the adjacency matrix is also a symmetry of the corresponding Laplacian matrix and vice versa (*24*). Suppose 𝒢 is a group of permutations of the nodes of the network that leaves the coupling matrix *L* invariant. Then, for each *g* ∈ 𝒢, we have a permutation matrix *R*_{g} that operates on the set of all node vectors **x** = (**x**_{1},…,**x**_{N})^{T} . Because *R*_{g}*L* = *LR*_{g}, this means *d*(*R*_{g}**x**_{i})/*dt* = **F**(*R*_{g}**x**_{i}) + Σ_{j} *L*_{ij}**H**(*R*_{g}**x**_{j} ), that is, the symmetry operation leaves the equations of motion unchanged. Hence, the subset of nodes permuted among each other by the group will remain synchronized if started in a synchronized state. We will refer to these subsets of nodes as clusters. The synchronized states for each cluster are flow-invariant.

An approach to the construction of all allowed CS patterns in a network has been proposed by Kamei and Cock (*29*). Although the method presented in their work (*29*) is general, because it applies to any network topology, it is computationally expensive. Here, we argue that, for the case of symmetric networks, a faster approach can often be followed that takes advantage of computational group theory, which is quite efficient. At each step, we show the results of this method applied to a particular case.

We first decompose the symmetry group 𝒢 into subgroups ℋ_{i}; *i* = 1,…,ν, each of which acts only on some subset of clusters (often only one) but not on any of the others (*30*, *31*). For this reason, we will refer to these subgroups as cluster groups, and the original group is a direct product of all cluster groups 𝒢 = ℋ_{1}×…×ℋ_{ν}. We further decompose each cluster group ℋ_{i} into all of its possible subgroups, which will give us a natural set of symmetry-breaking paths. These subgroups provide the full range of possible symmetry clusters, from the original full symmetry clusters to subclusters to the trivial case where each node is in its own cluster, that is, no symmetries. Whereas for the case of the adjacency matrix, this allows one to find all of the possible CS patterns, in the case of the Laplacian matrix, these patterns are certainly valid, but there may be other valid patterns that are not predicted by the computational group theory analysis. Extra steps are thus required to find these additional patterns.

We examine a particular case, a five-node network (Fig. 2), with adjacency matrix(5)and Laplacian matrix(6)Figure 2 shows all the allowed patterns, where nodes belonging to the same cluster are colored the same.

Cluster patterns A1 to A5 in Fig. 2 can be found by performing an analysis of the group 𝒢 and all of its subgroups. These are all the patterns that may emerge for the adjacency matrix (Eq. 5) and the only patterns that may emerge by symmetry for the Laplacian matrix (Eq. 6). The orbits of the original symmetry group are {1,2,3,4} and {5} by itself and are associated with two cluster groups: ℋ_{1}, which permutes the first, and ℋ_{2}, which is only the identity for the second single-node cluster. Other patterns (A2 to A5) are possible based on subgroups of the original group. These are the results of symmetry-breaking bifurcations.

Now, we create new potential clusters by first choosing a set of cluster groups (ℋ_{i}) and/or their subgroups (one for each cluster group). Together, these determine a subgroup 𝒢′ of the original group 𝒢. Second, we combine or merge some of these clusters as candidates for new synchronized clusters (it is our choice which to try). Because we are merging clusters or subclusters from different original symmetry clusters, the resulting CS patterns will not be the result of symmetries but may be dynamically valid when the coupling is Laplacian. Third, we can set those dynamical variables **x**_{i} equal to others in the merged cluster and see whether their equations of motion are the same (guaranteeing flow invariance). However, examining equations of motion for large networks by eye can become prohibitive.

There is a more direct way of doing the third step above, using the power and efficiency of group theory and computational discrete algebra tools (*32*, *33*). Examining the coupling term in Eq. 4, when clusters are synchronized, the diagonal feedback term [**H**(**x**_{i})] of a node will cancel the coupling terms of nodes from the same merged cluster. Hence, in the synchronized state, the dynamics behave as though nodes from the same merged cluster are not connected. We therefore form a dynamically equivalent coupling matrix, which is the original Laplacian matrix with off-diagonal components between nodes from merged clusters set to 0 and the diagonal values then set to the negative of the new row sums. We then perform the cluster group decompositions and subgroup constructions on the new, dynamically equivalent Laplacian (note that if the original Laplacian matrix is symmetric, so are the dynamically equivalent matrices obtained through this construction). If some set of subgroups of this new coupling matrix has symmetries yielding clusters that are our merged clusters, then their dynamics are flow-invariant in the synchronized state and cluster merging is possible. In this case, we call the new clusters Laplacian clusters. All the above can be automated in software. Hence, in Laplacian networks, even when a CS pattern is not the direct result of a symmetry of the original coupling matrix, it is the result of a symmetry of a dynamically equivalent coupling matrix. This is particularly relevant in terms of computational complexity. Although the problem of finding all of the symmetry operations of a given matrix has not been proven to be polynomial, efficient discrete algebra routines have been developed that make these computations possible, even for very large networks [for example, the Internet at the autonomous system level, for which *N* = 22,332; see the work by MacArthur *et al*. (*31*)]. This is in contrast with the method proposed by Kamei and Cock (*29*), which is not based on evaluation of the symmetries, and has been shown to become inefficient even for networks of moderate size (for example, *N* = 15).

To conclude, when the network is described by an adjacency matrix, a full characterization of all the dynamically valid patterns can be obtained by taking advantage of available computational discrete algebra tools (*32*, *33*). These output the cluster groups ℋ_{1},ℋ_{2},…,ℋ_{ν} and a decomposition of each cluster group ℋ_{i} into all of its possible subgroups, which provides a natural set of symmetry-breaking paths. As discussed above, available computational group theory routines perform these tasks very efficiently when compared to other possible methods [for example, (*29*)]. For the case of the adjacency matrix, our approach is always better than the state of the art. When the coupling is in Laplacian form, all of the (symmetry-related) CS patterns of the associated adjacency matrix are maintained. Moreover, it is possible for some of the cluster groups ℋ_{1},ℋ_{2},…,ℋ_{ν} and some of their subgroups to merge to form new dynamically valid clusters. To find all of the dynamically valid mergings, it can be helpful to use the clusters and subclusters provided by the symmetry analysis as the building blocks of our algorithm. If the number of these different clusters and subclusters is equal to μ, the number of tests that need to be performed (for pairs, triplets, and so on) is upper-bounded by *B*_{μ}, the μ-th Bell number (*34*); hence, it grows combinatorially with μ. For all the networks that do not display many symmetries, that is, for which μ << *N*, our approach based on computational group theory will be much faster than the one presented by Kamei and Cock (*29*). Hence, it is generally a good idea to run a preliminary analysis of the symmetries of the network to assess whether μ << *N* and choose the most convenient approach based on this outcome. It should be noted that there is no need to test mergings between subgroups of the same group, which reduces the complexity of the calculations. In general, finding all of the dynamically valid patterns for Laplacian networks may be substantially harder than for networks described by a symmetric adjacency matrix. Similar limitations were observed in the method by Kamei and Cock (*29*).

## STABILITY ANALYSIS AND EXPERIMENTAL VALIDATION

Stability of CS has been investigated for phase oscillators (*23*) and Stuart-Landau oscillators (*35*) and for lattices of coupled systems (*27*, *28*). However, a general approach to analyzing stability and multistability of CS patterns in arbitrary networks has not been developed. In (*24*), we studied a particular CS pattern corresponding to a minimum number of clusters (that is, maximal symmetry) for the case that the connectivity of the network is in the form of an adjacency matrix.

We now develop variational equations for the merged-cluster system so we can calculate the stability of each one of the allowed CS patterns. Although a number of papers (*23*, *26*–*28*, *35*–*37*) have dealt with CS in networks, only Nicosia *et al*. (*23*) and Poel *et al*. (*35*) have considered the problem of stability for particular dynamics of the individual systems. Belykh *et al*. (*27*), Belykh *et al*. (*28*), and Belykh and Hasler (*37*) have emphasized that for arbitrary systems, this is a difficult problem. To analyze stability, we start with the subgroup 𝒢′ of the original group that generated the clusters that we want to merge. It is formed by a direct product of the subgroups that we choose to use in our merged system. Using *C*_{m} to represent each cluster of nodes, *m* = 1,…,*K*, where *K* = number of clusters in 𝒢′, we have the variational equation of Eq. 2
(7)where the *Nn*-dimensional vector **x**(*t*) = [**x**_{1}(*t*)^{T},**x**_{2}(*t*)^{T},…,**x**_{N}(*t*)^{T} ]^{T} , *L* = the Laplacian-coupling matrix, and *E*^{(m)} is an *N*-dimensional diagonal indicator matrix for each cluster such that is equal to 1 if *i* ∈ *C*_{m} and is equal to 0; otherwise, *i* = 1,…,*N*. Note that we must use the original Laplacian matrix and not the dynamically equivalent one, which is used only for detecting synchronization flow invariance.

As we showed in (*24*), we can first block-diagonalize the coupling matrix *L* using the irreducible representations (IRRs) of 𝒢′, which yields the transformed coupling matrix *L*′ for A3 of Fig. 2(8)

The upper left block represents the variations on the synchronization manifold. It is 3 × 3 because there are three different trajectories in A3 (three clusters). The lower right block represents the variations in the transverse manifold. These are the perturbations that take the system out of synchrony, so it is these whose stability we want to calculate. Suppose we merge the center node (*5*) with nodes 1 and 3 to form the first merged state shown in L3. Geometrically, what must happen is that the dimension of the synchronization manifold must decrease by 1 (from 3 to 2) and the transverse manifold must increase by 1 (from 2 to 3).

To obtain the new coordinates on the synchronization manifold, we note that basis vectors of a cluster on the synchronization manifold have a 1 in the position of each node that belongs to the cluster. For example, in A3, the cluster (1,3) will have a (unit) vector of the form (1,0,1,0,0) and (5) will have the (unit) vector of the form (0,0,0,0,1). The merged cluster (1,3,5) will have the synchronization manifold vector, which is their sum, (1,0,1,0,1). The transformation of this new synchronization vector to the IRR coordinates of *L*′ provides the new synchronization direction, and its orthogonal complement provides the new transverse direction. We use these two new vectors to transform the 2 × 2 subblock associated with the (1,3) and (5) clusters in the 3 × 3 synchronization block to reduce the 3 × 3 synchronization manifold and increase the transverse manifold. This results in the final variational equation for the L3 casewhere we have linearized about the new synchronized merged cluster states {*s*_{1},…,*s*_{m}}, **η** is the vector of variations of all nodes transformed to the merged coordinates as above, *D***F** and *D***H** are the Jacobians of the nodes’ vector field and coupling function, respectively, *J*^{(m)} are the transformed *E*^{(m)}, and(9)In Eq. 9, the new synchronization block (in the upper left-hand corner) represents the new clusters (1,3,5) and (2,4), and the new transverse direction is associated with the new diagonal value −5. Obviously, this can be generalized to more complex cluster mergings. An example of cluster merging is shown geometrically in Fig. 3. We note that this method can be used to analyze the stability of any dynamically valid CS pattern for which knowledge of the block-diagonalized matrix *L*″ is available, and it applies whether the connectivity is given by an adjacency matrix or by a Laplacian matrix. The generalization of the above procedure simply uses the synchronization vectors for all the new clusters as the rows of a matrix for which all other components are zero. An application of a singular value decomposition then gives a basis for the synchronization block (the original synchronization vectors) and a basis for the orthogonal complement that represents all the transverse directions. In some cases, the latter can be simplified by a block diagonalization when all members of the block are in the same merged cluster. This makes it possible to automate this block diagonalization to evaluate stability of all the CS patterns that can emerge in a given network topology.

We show symmetry breaking and the existence of Laplacian clusters using the experimental system described in detail in (*24*, *38*) and also in Methods and Materials herein. The dynamics of the system is modeled by a map according to(10)where σ is the coupling strength, which we will vary. Equation 10 is a map version of Eq. 2 and can show fixed-point, periodic, or chaotic dynamics depending on the values of the parameters. Here, δ = 0.525, β = 1.45π (which guarantees chaotic behavior), and σ was decreased from π to 0. The dynamics of Eq. 10 describe the experiment, which uses a spatial light modulator (SLM) and laser system to display the behavior of the nodes of the five-node system shown in Fig. 2.

Figure 4 shows the experimental synchronization error for each synchronization pattern shown in Fig. 2 as a function of σ. The phases *x*_{i} were not reinitialized when σ was updated. The synchronization error for each pattern is computed as <|*x*_{i}(*t*) – (t)|> , where the symbol <•> indicates an average both in time and over the nodes *i* in a cluster and the clusters in a pattern. is the average state for the nodes in the cluster to which node *i* belongs. In the lower panel of Fig. 4, we plot the results of our stability analysis for each one of the CS patterns, where a colored dot labels the values of σ for which the corresponding pattern is stable. In particular, a CS pattern was indicated to be stable when (i) all the numerically computed maximum Lyapunov exponents corresponding to the transverse blocks were found to be negative and (ii) the synchronous pattern was asymptotically valid, that is, the CS pattern was observed after integrating its equations for a long time. The equations that were used to run these stability calculations are shown in Table 1.

In general, when two or more clusters merge into one, there are two independent effects on stability, as can be seen from the structure of the block-diagonalized matrix *L*″: the first one is that the dimension of the synchronization block decreases, which determines the motion in the synchronization manifold; the second one is that new transverse blocks appear, with the other (preexisting) transverse blocks remaining the same. As a consequence, because each one of these transverse blocks needs to be individually tested for stability, we expect that as the number of transverse blocks increases, the range of stability will decrease. This is confirmed by our experimental results plotted in Fig. 4, showing that the σ range of stability becomes smaller for CS patterns that are characterized by higher symmetry. Exceptions to this rule are possible because the motion in the synchronization manifold (on which the transverse Lyapunov exponents depend) may also affect stability in ways that cannot be predicted by the analysis of the transverse blocks only.

Figure 5 shows the dynamics on the synchronization manifold for the symmetry pattern A3 and the two merged patterns L3 and L1. In each transition A3 **→** L3 **→** L1, the dimension of the synchronization manifold decreases by 1 (from 3 to 2 to 1) and the transverse manifold increases by 1 (from 2 to 3 to 4). Figure 6 shows three snapshots of the experimental dynamics for each one of the patterns A3, L3, and L1. The parameter values were (i) β = 2.2π, σ = 0.9π; (ii) β = 0.5874π, σ = 0.486π; and (iii) β = 2.2π, σ = 0.5π. Videos (movies S1 and S2) of the synchronization manifolds and patterns of Fig. 6 are available in the Supplementary Materials.

## CONCLUSION

We have studied the emergence of CS in networks of coupled oscillators. First, we show how to obtain all of the possible dynamically valid CS patterns for an arbitrary network topology. The methods we illustrated above are general and will apply to nodes with any type of dynamics (ordinary differential equations, maps, etc.). This method also extends to directed networks, weighted networks, and labeled nodes (to represent different dynamics on each node), all of which are handled by a software package (*33*). An important point is that these techniques work for any subgroups of the cluster groups, which include the trivial subgroups (*24*) (only the identity in each). Hence, we can approach merging as a “bottom up” process and analyze any arbitrary merging of nodes into a cluster to determine whether the dynamics allow a synchronized state and whether they are stable in some parameter range. This means that all possible clusters can be analyzed using our approach. We note that combining clusters as a “top down” approach would provide clusters that most likely are not from symmetries, whereas bottom up would be a process that would include clusters easily obtained from symmetries, although it would be useful for cases where one has particular clusters in mind to analyze.

Our main result is a technique to evaluate stability of all the dynamically valid CS patterns for both networks for which the connectivity is given by an adjacency matrix and by a Laplacian matrix. We predict that the range of stability typically becomes smaller for CS patterns that are characterized by higher symmetry, which is confirmed in our experimental system.

## MATERIALS AND METHODS

### Experimental design

The system uses an SLM and a camera in a feedback configuration. The camera has a focal plane array (FPA) of 320 × 256 pixels and an area of 8 × 6.4 mm^{2}. The SLM has a resolution of 512 × 512 pixels and an active area of 7.68 × 7.68 mm^{2}. A light-emitting diode with a wavelength of 1550 nm is used to illuminate the modulator. The light passes through a polarizing beam splitter and a quarter-wave plate (QWP), so that circularly polarized light is incident on the SLM. The SLM imparts a programmable spatially varying phase shift *x* between the two polarization components of the reflected light. The reflected light passes through the QWP and polarizer and is imaged onto a 256 × 256–pixel square region of the camera’s FPA. The relationship between the phase shift *x* applied by the SLM and the normalized intensity *I* recorded by the camera is *I*(*x*) = (1 – cos *x*)/2.

Each oscillator corresponds to a square patch of 16 × 16 pixels on the SLM, which is imaged onto an 8 × 8–pixel region of the camera’s FPA. Using a computer, the phase shift of the *i*th region, **x**_{i}, is iteratively updated on the basis of the intensity measured by the camera according to Eq. 10, where σ is the coupling strength, which we will vary from π to 0.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/4/e1501737/DC1

movie S1. Video of light intensity dynamics of node areas for the five-node network analyzed in the text, displayed in Fig. 2, and whose equations are shown in Table 1 for different values of parameters β and σ as shown in Fig. 6.

movie S2. Video rotating the view of state space trajectories of the synchronized node clusters for the five-node network analyzed in the text, displayed in Fig. 2, and whose equations are shown in Table 1 for different values of parameters β and σ as shown in Fig. 5.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We acknowledge help with computational group theory from D. Joyner of the U.S. Naval Academy and information and computer code from B. D. MacArthur and R. J. Sánchez-García, both of the University of Southampton.

**Funding:**This work was funded by the Office of Naval Research at the Naval Research Laboratory. We also gratefully acknowledge the Office of Naval Research for support of the experimental research at the University of Maryland, College Park, through grant no. N000141410443. F.S. was supported by NSF grants CMMI-1400193 and CRISP-1541148.

**Author contributions:**F.S. conceived cluster merging, analyzed stability, and performed simulations. L.M.P. performed group theory analysis. A.M.H. performed experiments and simulations. T.E.M. and R.R. contributed to experimental design and guidance. All authors contributed to writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2016, The Authors