Topological Control of Synchronization Patterns: Trading Symmetry for Stability

TL;DR

This paper reveals that breaking symmetry in network structures can enhance synchronization stability, challenging traditional views, and provides algorithms for topological control of synchronization patterns.

Contribution

It demonstrates that symmetry breaking can improve synchronization stability in networks, supported by experimental validation and optimization algorithms.

Findings

Symmetry breaking enhances cluster synchronizability.

Experimental validation confirms robustness against noise.

Optimization algorithms improve network structure for desired synchronization.

Abstract

Symmetries are ubiquitous in network systems and have profound impacts on the observable dynamics. At the most fundamental level, many synchronization patterns are induced by underlying network symmetry, and a high degree of symmetry is believed to enhance the stability of identical synchronization. Yet, here we show that the synchronizability of almost any symmetry cluster in a network of identical nodes can be enhanced precisely by breaking its structural symmetry. This counterintuitive effect holds for generic node dynamics and arbitrary network structure and is, moreover, robust against noise and imperfections typical of real systems, which we demonstrate by implementing a state-of-the-art optoelectronic experiment. These results lead to new possibilities for the topological control of synchronization patterns, which we substantiate by presenting an algorithm that optimizes the structure of individual clusters under various constraints.

Figures

Tables

Table 1. Table 1: Connected symmetry clusters of 6 nodes and optimal clusters embedded within them. Some symmetry clusters have more than one embedded optimal network, in which case we show one that can be obtained through a minimal number of link deletions.

Symmetry clusters
Eigenratio	4	2.5	2	1.5	1
Optimal clusters
Eigenratio	1	1	1	1	1

Equations

$$\dot{{\mathbf{x}}_i}=\mathbf{f}({\mathbf{x}}_i)-\sigma \sum _{j=1}^n{L}_{ij}\mathbf{h}({\mathbf{x}}_j),$$

$$\begin{split}\dot{\bm{x}_{i}}&=\bm{f}(\bm{x}_{i})-\sigma\sum_{j=1}^{m}L_{ij}\bm{h}(\bm{x}_{j})+\sigma\sum_{j=m+1}^{n}A_{ij}\bm{h}(\bm{x}_{j})\\ &=\bm{f}(\bm{x}_{i})-\sigma\sum_{j=1}^{m}L_{ij}\bm{h}(\bm{x}_{j})+\sigma\bm{I}\big{(}\{\bm{x}_{j}\}_{j>m}\big{)},\end{split}$$

$$L_{ij}=\{\begin{array}{ll}{\tilde{L}}_{ij},& 1\le i\ne j\le m,\\ {\tilde{L}}_{ij}+\tilde{\mu },& 1\le i=j\le m,\end{array}$$

$$\delta \dot{\mathbf{X}}=\left[\text{mathds}{1}_m\otimes J\mathbf{f}({\mathbf{s}}_I)-\sigma \hat{\mathbf{L}}\otimes J\mathbf{h}({\mathbf{s}}_I)\right]\delta \mathbf{X},$$

$$\dot{\bm{\eta}}_{i}=\big{[}J\bm{f}(\bm{s}_{I})-\sigma\widehat{v}_{i}J\bm{h}(\bm{s}_{I})\big{]}\bm{\eta}_{i},$$

$$x_i(t+1)=\beta I[x_i(t)]-\sigma \sum _{j=1}^n{L}_{ij}I[x_j(t)]\text{mod}2\pi ,$$

$$\Vert {\mathbf{\eta }}_i\Vert \to 0\text{ in lx@crefcreftype refnumeq:3 for all transverse modes.}$$

$$\{\begin{array}{l}\delta \dot{\mathbf{X}}=\left(\text{mathds}{1}_m\otimes J\mathbf{f}({\mathbf{s}}_{I_X})-\sigma {\hat{\mathbf{L}}}_X\otimes J\mathbf{h}({\mathbf{s}}_{I_X})\right)\delta \mathbf{X}+\sigma \mathbf{C}\otimes J\mathbf{h}({\mathbf{s}}_{I_Y})\delta \mathbf{Y},\\ \delta \dot{\mathbf{Y}}=\left(\text{mathds}{1}_{m^{\prime }}\otimes J\mathbf{f}({\mathbf{s}}_{I_Y})-\sigma {\hat{\mathbf{L}}}_Y\otimes J\mathbf{h}({\mathbf{s}}_{I_Y})\right)\delta \mathbf{Y}+\sigma \mathbf{D}\otimes J\mathbf{h}({\mathbf{s}}_{I_X})\delta \mathbf{X}.\end{array}$$

$$\{\begin{array}{l}\delta \dot{\mathbf{X}}=\left(\text{mathds}{1}_m\otimes J\mathbf{f}({\mathbf{s}}_{I_X})-\sigma {\hat{\mathbf{L}}}_X\otimes J\mathbf{h}({\mathbf{s}}_{I_X})\right)\delta \mathbf{X},\\ \delta \dot{\mathbf{Y}}=\left(\text{mathds}{1}_{m^{\prime }}\otimes J\mathbf{f}({\mathbf{s}}_{I_Y})-\sigma {\hat{\mathbf{L}}}_Y\otimes J\mathbf{h}({\mathbf{s}}_{I_Y})\right)\delta \mathbf{Y}.\end{array}$$

$$\Vert \delta \mathbf{X}\Vert \to 0\text{ and }\Vert \delta \mathbf{Y}\Vert \to 0\text{ in lx@crefcreftype refnumeq:s7}$$

$$\Vert \delta \mathbf{X}\Vert \to 0\text{ and }\Vert \delta \mathbf{Y}\Vert \to 0\text{ in lx@crefcreftype refnumeq:s8}.$$

$$\delta \dot{\mathbf{X}}=\left[\text{mathds}{1}_m\otimes J\mathbf{f}({\mathbf{s}}_I)-\sigma \hat{\mathbf{L}}\otimes J\mathbf{h}({\mathbf{s}}_I)\right]\delta \mathbf{X},$$

$$\mathbf{B}={\mathbf{P}}^{-1}\hat{\mathbf{L}}\mathbf{P}=\left(\begin{array}{cccc}\tilde{\mu }& & & \\ & {\mathbf{B}}_1& & \\ & & \ddots & \\ & & & {\mathbf{B}}_q\end{array}\right),{\mathbf{B}}_j=\left(\begin{array}{cccc}{\hat{v}}_{j+1}& & & \\ 1& {\hat{v}}_{j+1}& & \\ & \ddots & \ddots & \\ & & 1& {\hat{v}}_{j+1}\end{array}\right),$$

$${\dot{\mathbf{\eta }}}^{(j)}=\left[\text{mathds}{1}_k\otimes J\mathbf{f}({\mathbf{s}}_I)-\sigma {\mathbf{B}}_j\otimes J\mathbf{h}({\mathbf{s}}_I)\right]{\mathbf{\eta }}^{(j)},j=1,\dots ,q.$$

$$\begin{array}{rl}{\dot{\mathbf{\eta }}}_1^{(j)}=& \left[J\mathbf{f}({\mathbf{s}}_I)-\sigma {\hat{v}}_{j+1}J\mathbf{h}({\mathbf{s}}_I)\right]{\mathbf{\eta }}_1^{(j)},\\ {\dot{\mathbf{\eta }}}_2^{(j)}=& \left[J\mathbf{f}({\mathbf{s}}_I)-\sigma {\hat{v}}_{j+1}J\mathbf{h}({\mathbf{s}}_I)\right]{\mathbf{\eta }}_2^{(j)}-\sigma J\mathbf{h}({\mathbf{s}}_I){\mathbf{\eta }}_1^{(j)},\\ \cdots & \\ {\dot{\mathbf{\eta }}}_k^{(j)}=& \left[J\mathbf{f}({\mathbf{s}}_I)-\sigma {\hat{v}}_{j+1}J\mathbf{h}({\mathbf{s}}_I)\right]{\mathbf{\eta }}_k^{(j)}-\sigma J\mathbf{h}({\mathbf{s}}_I){\mathbf{\eta }}_{k-1}^{(j)}.\end{array}$$

$$x_i(t+1)=r{x}_i(t)-\sigma \sum _{j=1}^n{L}_{ij}{x}_j(t)\text{mod}2\pi ,$$

$$\{\begin{array}{l}{x}_i(t+1)=1-a{x}_i^2(t)+y_i(t)-\sigma {\sum }_{j=1}^n{L}_{ij}{y}_j(t),\\ y_i(t+1)=b{x}_i(t),\end{array}$$

Full text

DOI: 10.1103/PhysRevLett.122.058301

Symmetry and synchronization are interrelated concepts in network systems. Synchronization, being a symmetric state among oscillators, has its existence and stability influenced by the symmetry of the network stewart2003symmetry ; nicosia2013remote ; aguiar2011dynamics . For example, recent research has shown that network symmetry can be systematically explored to identify stable synchronization patterns in complex networks pecora2014cluster . Different work has shown that structural homogeneity (and hence a higher degree of network symmetry) usually enhances synchronization stability donetti2005entangled ; denker2004breaking ; nishikawa2003heterogeneity . Any given network of identical oscillators can always be partitioned into so-called symmetry clusters macarthur2008symmetry , characterized as clusters of oscillators that are identically coupled, both within the cluster and to the rest of the network, making them natural candidates for cluster synchronization pecora2014cluster ; sorrentino2016complete . Cluster synchronization has been investigated in numerous experimental systems, including networks of optoelectronic oscillators pecora2014cluster ; sorrentino2016complete ; williams2013experimental , semiconductor lasers nixon2012controlling ; argyris2016experimental , Boolean systems rosin2013control , neurons vardi2012synchronization , slime molds takamatsu2001spatiotemporal , and chemical oscillators totz2015phase . Many of these experiments explicitly investigated the beneficial impact of network symmetries on cluster formation totz2015phase ; takamatsu2001spatiotemporal ; pecora2014cluster ; sorrentino2016complete ; hart2016experimental . Taken together, previous results support the expectation that oscillators that are indistinguishable on structural grounds are also more likely to exhibit indistinguishable (synchronous) dynamics.

In this Letter, we investigate the relation between symmetry and synchronization in the general context of cluster synchronization (including global synchronization). We show that, in order to induce stable synchronization, one often has to break the underlying structural symmetry. This counterintuitive result holds for the general class of networks of diffusively coupled identical oscillators with a bounded and connected stability region, and it follows rigorously from our demonstration that almost all clusters exhibiting optimal synchronizability are necessarily asymmetric. In particular, the synchronizability of almost any symmetry cluster can be enhanced precisely by breaking the internal structural symmetry of the cluster. These findings add an important new dimension to the recent discovery of parametric asymmetry-induced synchronization PhysRevLett.117.114101 ; zhang2017asymmetry ; zhang2017nonlinearity , a scenario in which the synchronization of identically coupled identical oscillators is enhanced by assigning nonidentical parameters to the oscillators. Here, we show that synchronization of identically coupled identical oscillators is enhanced by changing the connection patterns of the oscillators to be nonidentical. We refer to this effect as structural asymmetry-induced synchronization (AISync). We confirm that this behavior is robust against noise and can be found in real systems by providing the first experimental demonstration of structural AISync using networks of coupled optoelectronic oscillators. In excellent agreement with theory, the experiments show unequivocally that both intertwined and nonintertwined clusters can be optimized by reducing structural symmetry.

We consider a network of $n$ diffusively coupled identical oscillators,

[TABLE]

where $\bm{x}_{i}$ is the state of the $i$th oscillator, $\bm{f}$ is the vector field governing the uncoupled dynamics of each oscillator, $\bm{L}=\{L_{ij}\}$ is the Laplacian matrix describing the structure of an arbitrary unweighed network, $\bm{h}$ is the interaction function, and $\sigma>0$ is the coupling strength. We are interested in the dynamics inside a symmetry cluster. To facilitate presentation, we first assume that the cluster is nonintertwined pecora2014cluster ; cho2017stable ; that is, it can synchronize independent of whether other clusters synchronize or not. The general case of intertwined clusters—in which desynchronization in one cluster can lead to loss of synchrony in another cluster—requires considering the intertwined clusters concurrently, and this important case is addressed after our analysis of nonintertwined clusters.

Numbering the oscillators in that cluster from $1$ to $m$, we obtain the dynamical equation for the cluster:

[TABLE]

where $L_{ij}=\delta_{ij}\mu_{i}-A_{ij}$, $\bm{A}=\{A_{ij}\}$ is the adjacency matrix of the network, $\mu_{i}$ is the indegree of node $i$, and the equation holds for $1\leq i\leq m$. Here, we denote the input term from the rest of the network $\sum_{j=m+1}^{n}A_{ij}\bm{h}(\bm{x}_{j})$ by \bm{I}\big{(}\{\bm{x}_{j}\}_{j>m}\big{)} to emphasize that this term is independent of $i$ and hence equal for all oscillators $1,\dots,m$. This term is zero only when the cluster receives no connection from the rest of the network, such as the important case in which the entire network consists of a single symmetry cluster (i.e., $m=n$).

For $m<n$, if we regard the cluster subnetwork consisting of oscillators $1,\dots,m$ as a separate network (by ignoring its connections with other clusters), then its $m\times m$ Laplacian matrix $\widetilde{\bm{L}}$ is closely related to the corresponding block of the $n\times n$ Laplacian matrix $\bm{L}$ of the full network:

[TABLE]

where $\widetilde{\mu}\geq 0$ is the number of connections each oscillator in the cluster receives from the rest of the network. It is then clear that there are two differences in the dynamical equation when the cluster subnetwork is part of a larger network [i.e., as a symmetry cluster, described by Eq. 2] rather than as an isolated network. First, the Laplacian matrix $\widetilde{\bm{L}}$ in the dynamical equation is replaced by $\widehat{\bm{L}}=\{L_{ij}\}_{1\leq i,j\leq m}=\widetilde{\bm{L}}+\widetilde{\mu}\mathds{1}_{m}$; that is, the diagonal entries are uniformly increased by $\widetilde{\mu}$. Second, each oscillator now receives a common input \sigma\bm{I}\big{(}\{x_{j}\}_{j>m}\big{)} produced by its coupling with other clusters, which generally alters the synchronization trajectory $\bm{s}_{I}\equiv\bm{x}_{1}=\cdots=\bm{x}_{m}$, causing it to be typically different from the ones generated by the uncoupled dynamics $\dot{\bm{s}}=\bm{f}(\bm{s})$. This has to be accounted for when calculating the maximum Lyapunov exponent transverse to the cluster synchronization manifold to determine the stability of the cluster synchronous state.

Despite these differences, a diagonalization procedure similar to the one used in the master stability function approach pecora1998master can still be applied to the variational equation in order to assess the cluster’s synchronization stability. The variational equation describing the evolution of the deviation away from $\bm{s}_{I}$ inside the cluster can be written as

[TABLE]

where $\delta\bm{X}=(\delta\bm{x}_{1}^{\intercal},\cdots,\delta\bm{x}_{m}^{\intercal})^{\intercal}=(\bm{x}_{1}^{\intercal}-\bm{s}_{I}^{\intercal},\cdots,\bm{x}_{m}^{\intercal}-\bm{s}_{I}^{\intercal})^{\intercal}$ and $\otimes$ denotes the Kronecker product. The rest of the network does not enter the equation explicitly, other than through its influence on the coupling matrix $\widehat{\bm{L}}$ and the synchronization trajectory $\bm{s}_{I}$. If $\widehat{\bm{L}}$ is diagonalizable (as for any undirected network), the decoupling of Eq. 4 results in $m$ independent $d$-dimensional equations corresponding to individual perturbation modes:

[TABLE]

where $d$ is the dimension of node dynamics, $J$ is the Jacobian operator, $\bm{\eta}=(\bm{\eta}_{1}^{\intercal},\cdots,\bm{\eta}_{m}^{\intercal})^{\intercal}$ is $\delta\bm{X}$ expressed in the new coordinates that diagonalize $\widehat{\bm{L}}$, and $\widehat{v}_{i}=\widetilde{v}_{i}+\widetilde{\mu}$ are the eigenvalues of $\widehat{\bm{L}}$ in ascending order of their real parts [with $\left\{\widetilde{v}_{i}\right\}=\mbox{eig}(\widetilde{\bm{L}})$]. If $\widehat{\bm{L}}$ is not diagonalizable nishikawa2006maximum , the analysis can be carried out by using the Jordan canonical form of this matrix to replace diagonalization by block diagonalization, as explicitly shown in the Supplemental Material SM . In both cases the cluster synchronous state is stable if $\Lambda(\sigma\widehat{v}_{i})<0$ for $i=2,\dots,m$, where $\Lambda$ is the largest Lyapunov exponent of Eq. 5 and $\widehat{v}_{2},\dots,\widehat{v}_{m}$ represent the transverse modes; the maximum transverse Lyapunov exponent (MTLE) determining the stability of the synchronous state is $\max_{2\leq i\leq m}\Lambda(\sigma\widehat{v}_{i})$. Moreover, for the large class of oscillator networks for which the stability region is bounded and connected barahona2002synchronization ; li2010consensus ; flunkert2010synchronizing ; huang2009generic , as assumed here and verified for all models we consider comment1 , the synchronizability of a cluster can be quantified in terms of the eigenratio $R=\text{Re}(\widetilde{v}_{m})/\text{Re}(\widetilde{v}_{2})$: the smaller this ratio, in general, the larger the range of $\sigma$ over which the cluster synchronous state can be stable. The cluster subnetwork is most synchronizable when $\widetilde{v}_{2}=\cdots=\widetilde{v}_{m}$, which also implies that all eigenvalues are real and in fact integers if the network is unweighted nishikawa2010network , as considered here. It is important to notice that the optimality of the cluster subnetwork is conserved in the sense that if $\widetilde{v}_{2}=\cdots=\widetilde{v}_{m}$ for the isolated cluster, then $\widehat{v}_{2}=\cdots=\widehat{v}_{m}$ will hold for the cluster as part of a larger network. Since the analysis above does not invoke the continuity of the equations anywhere, it holds for discrete-time systems as well. In this case one can simply replace $\delta\dot{\bm{X}}$ and $\delta\bm{X}$ in Eq. 4 by $\delta\bm{X}(t+1)$ and $\delta\bm{X}(t)$, respectively.

Now we can compare symmetry clusters with optimal clusters and show rigorously that almost all optimally synchronizable clusters are asymmetric. Without loss of generality, we consider an unweighted cluster in isolation and assume it has $m$ nodes and $\ell$ directed links internal to the cluster. In a symmetry cluster, because the nodes are structurally identical, the in- and outdegrees of all nodes must be equal. Thus, $\ell$ must be divisible by $m$ if the cluster is symmetric. In an optimal cluster, because $\widetilde{v}_{2}=\cdots=\widetilde{v}_{m}\equiv\widetilde{v}$ and thus tr$(\widetilde{\bm{L}})=(m-1)\widetilde{v}$, it follows that $\widetilde{v}=\ell/(m-1)$. The fact that $\widetilde{v}$ is an integer implies that $\ell$ must be divisible by $m-1$ if the cluster is optimal. Since $\ell\leq m(m-1)$, the two divisibility conditions can be satisfied simultaneously if and only if $\ell=m(m-1)$ (i.e., when the isolated cluster is a complete graph). But there are numerous optimal clusters for $\ell<m(m-1)$ nishikawa2006maximum ; nishikawa2010network . Therefore, for any given number $m$ of nodes, all optimal clusters other than the complete graph are necessarily asymmetric, meaning that (with the exception of the complete graph) the synchronization stability of any symmetry cluster can be improved by breaking its structural symmetry comment . This general conclusion forms the basis of structural AISync and holds, in particular, when an entire network consists of a single symmetry cluster.

When viewed as isolated subnetworks, symmetry clusters are equivalent to the vertex-transitive digraphs in algebraic graph theory, defined as directed graphs in which every pair of nodes is equivalent under some node permutation biggs1993algebraic ; mckay2014practical . Thus, in order to improve the synchronizability of any nonintertwined symmetry cluster from an arbitrary network, we only need to optimize the corresponding vertex-transitive digraph by manipulating its (internal) links. In particular, this can always be done by removing links inside the symmetry cluster nishikawa2006synchronization ; nishikawa2010network , despite the fact that sparser networks are usually harder to synchronize. For concreteness, we focus on clusters that are initially undirected and consider the selective removal of individual directional links. As an example, we show in Table 1 all connected undirected symmetry clusters of 6 nodes and their embedded optimal networks. Apart from the complete graph, which is already optimal to begin with, the synchronizability of the other symmetry clusters as measured by the eigenratio $R$ is significantly improved in all cases.

Because in practice it can be costly or unnecessary to fully optimize a symmetry cluster, it is natural to ask whether its synchronizability can be significantly improved by just modifying a few links. We developed an efficient algorithm for this purpose and summarize the statistical results based on all connected undirected symmetry clusters of sizes between $m=8$ and $17$ in the Supplemental Material SM . On average, only about $14\%$ of the links need to be rewired to reduce $R-1$ by half and thus significantly improve synchronizability of symmetry clusters. This illustrates the potential of structural AISync as a mechanism for the topological control of synchronization stability. Our simulated annealing code to improve cluster synchronizability is available in Ref. SA . This algorithm can also be used to demonstrate structural AISync in global synchronization, as shown in the Supplemental Material SM .

Having established a theoretical foundation for our main finding, we now turn to our experimental results. The experiments are performed using networks of identical optoelectronic oscillators whose nonlinear component is a Mach-Zehnder intensity modulator. The system can be modeled as

[TABLE]

where $t$ is now a discrete time, $\beta$ is the feedback strength, $I(x_{i})=\sin^{2}(x_{i}+\delta)$ is the normalized intensity output of the modulator, $x_{i}$ is the normalized voltage applied to the modulator, and $\delta$ is the operating point (set to $\pi/4$ in our experiments). Each oscillator consists of a clocked optoelectronic feedback loop. Light from a $780$ nm continuous-wave laser passes through the modulator, which provides the nonlinearity. The light intensity is converted into an electrical signal by a photoreceiver and measured by a field-programmable gate array (FPGA) via an analog-to-digital converter (ADC). The FPGA is clocked at 10 kHz, resulting in the discrete-time map dynamics of the oscillators. The FPGA controls a digital-to-analog converter (DAC) that drives the modulator with a voltage $x_{i}(t+1)=\beta I[x_{i}(t)]$, closing the feedback loop. The oscillators are coupled together electronically on the FPGA according to the desired Laplacian matrix. Specifically, the experimental system uses time multiplexing and time delays to realize a network of coupled oscillators from a single time-delayed feedback loop, as described in detail in Ref. hart2017experiments . A schematic illustration of the experimental setup can be found in the Supplemental Material SM .

We first consider the network configuration shown in Fig. 1(a), which is a complex network with five symmetry clusters. The symmetry cluster highlighted in magenta is nonintertwined, and can be optimized by removing the red dashed links. The MTLE calculation in Fig. 1(b) predicts AISync to be common in the parameter space. Fixing $\beta=6$, we performed $8$ runs of the experiment starting from different random initial conditions, and measured the normalized voltages $x_{i}$ for $8196$ iterations at each fixed coupling strength before increasing $\sigma$ by 0.015. The synchronization error is defined as $\Delta=\sqrt{\sum_{1\leq i\leq m}\|x_{i}-\bar{x}\|^{2}/m}$, where $\bar{x}$ is the mean inside the cluster. The data points in Fig. 1(c) correspond to the average synchronization error $\langle\Delta\rangle$, defined as $\Delta$ averaged over the last $5000$ iterations for each $\sigma$ and then further averaged over the $8$ experimental runs. The error bars corresponding to the standard deviation across different runs are smaller than the size of the symbols. One can observe AISync over a wide range of the coupling strength $\sigma$, matching the theoretical prediction. Structural AISync is also common for different oscillator types and network structures and is robust against noise and parameter mismatches, as demonstrated systematically in the Supplemental Material SM .

We now turn to the case of intertwined clusters. Consider two intertwined clusters $X$ and $Y$ subject to transverse perturbations $\delta\bm{X}$ and $\delta\bm{Y}$, respectively. The variational equation for $\delta\bm{X}$ has the same form as Eq. 4 except for an additional cross-coupling term $\sigma\bm{C}\otimes J\bm{h}(\bm{s}_{I_{Y}})\delta\bm{Y}$ added to the right, where $\bm{C}$ is the adjacency matrix describing the intercluster coupling from cluster $Y$ to cluster $X$. The variational equation for $\delta\bm{Y}$ is defined similarly. Now, if $\delta\bm{X}$ ($\delta\bm{Y}$) does not converge to zero according to Eq. 4, then the cross-coupling term must not vanish and $\|\delta\bm{Y}\|$ ($\|\delta\bm{X}\|$) must stay away from zero in order for $\|\delta\bm{X}\|\rightarrow 0$ ($\|\delta\bm{Y}\|\rightarrow 0$) in the full variational equation. Thus, in order to stabilize synchronization in intertwined clusters, the following condition must be satisfied for each cluster:

[TABLE]

In other words, $\|\delta\bm{X}\|$ and $\|\delta\bm{Y}\|$ converging to zero in Eq. 4 is a necessary condition for stable synchronization in $X$ and $Y$. Because optimizing the clusters independently (as if they were nonintertwined) is guaranteed to expand the region satisfying the condition in Eq. 7, such independent optimization is an effective strategy for improving synchronization in intertwined clusters. For more details on this analysis, see Supplemental Material SM .

We demonstrate the strength of our approach on a random network containing two intertwined clusters, which are highlighted in Fig. 2(a). Each cluster is optimized by removing the red dashed links, which breaks the structural symmetry but reduces the eigenratio of the cluster to 1. The orange shade in Fig. 2(b) indicates the region where the condition in Eq. 7 is satisfied by the original clusters. The region satisfying this condition is expanded to include the purple region when the clusters are optimized. Direct simulations allow us to identify a large parameter region exhibiting AISync, which is highlighted in dark shades in Fig. 2(b) and is included mainly in the expanded (purple) region. A small portion of the AISync region also extends into the orange region, which follows from the condition in Eq. 7 being necessary but not sufficient for synchronization in the original clusters. To validate the theory and the numerics, we perform experiments with parameters varied quasistatically along the dashed line in Fig. 2(b). As shown in Fig. 2(c), the symmetry clusters are both incoherent for the entire range of parameters studied. The two optimized clusters exhibit perfectly synchronized dynamics except at the very edge of the AISync region, where the noise in the ADC has a marked impact on the dynamics (nevertheless, they are still much more synchronized than the symmetry clusters). It is interesting to mention that although both optimized clusters are in synchrony themselves, they are not synchronized with each other.

In summary, we established the role of structural asymmetry (or structural heterogeneity) in promoting spontaneous synchronization through both theory and experiments. Our theory confirmed the generality of the phenomenon, while our experiments demonstrated its robustness. Because symmetry clusters arise naturally in complex networks, our findings are applicable to a wide range of coupled dynamical systems. In particular, since identical synchronization in a symmetry cluster is the basic building block of more complex synchronization patterns, our results can be used for the targeted topological control of cluster synchronization in complex networks, which echoes the positive effect of structural asymmetry on input control whalen2015observability .

The authors thank Alex Mercanti, Takashi Nishikawa, Don Schmadel, and Thomas E. Murphy for insightful discussions. This work was supported by ONR Grant No. N000141612481 (J.D.H. and R.R.) and ARO Grant No. W911NF-15-1-0272 (Y.Z. and A.E.M.).

J.D.H. and Y.Z. contributed equally to this work.

S1 Experimental implementation of the coupled optoelectronic oscillators

S2 Optimizing intertwined clusters

In this section, we provide more details on the optimization of intertwined clusters. When two clusters are intertwined, desynchronization in one cluster will in general lead to the loss of synchrony in the other cluster (an example would be two equal-sized rings coupled in one-to-one fashion). This is because the symmetry group acting on the two clusters does not admit a geometric decomposition; that is, symmetry permutations cannot be applied to each cluster independently. As a consequence, a desynchronized cluster sends incoherent signals to nodes in the other cluster, causing its intertwined counterpart to desynchronize as well. The irreducible representation transformation introduced in Ref. pecora2014cluster is a powerful formalism that enables stability analysis on many cluster synchronization patterns. In that framework, the presence of intertwined clusters is reflected in nontrivial transverse blocks (i.e., blocks with dimension greater than 1) in the transformed coupling matrix, whereas nonintertwined clusters only give rise to $1\times 1$ transverse blocks. Unfortunately, the high dimensionality of the transverse blocks makes the effect of topological perturbations on cluster synchronizability opaque, and thus the analysis of the transformed matrix offers little insight into how to optimize the clusters to support desired synchronization patterns.

We developed a new perspective that gives a simple necessary condition for the synchronization in intertwined clusters. This in turn points to an extension of the previous optimization scheme that is no longer limited to nonintertwined clusters.

Consider two intertwined clusters $X$ and $Y$ subject to transverse perturbations $\delta\bm{X}$ and $\delta\bm{Y}$, respectively. Their variational equation reads

[TABLE]

Here, $\bm{C}_{ij}=1$ if the $i$-th oscillator in cluster $X$ receives an input from the $j$-th oscillator in cluster $Y$ and $\bm{C}_{ij}=0$ otherwise. The intercluster coupling matrix $\bm{D}$ is similarly defined with the role of two clusters exchanged ($\bm{D}=\bm{C}^{\intercal}$ if the intercluster coupling is undirected). Without the cross-coupling term, Eq. S1 reduces to the nonintertwined case discussed in the main text

[TABLE]

Because of the intertwined nature of the two clusters, they must be considered concurrently when synchronization is desired in either of them. That is, $\widehat{\bm{L}}_{X}$ and $\widehat{\bm{L}}_{Y}$ should be optimized to ensure that $\delta\bm{X}$ and $\delta\bm{Y}$ both vanish in Eq. S1.

It is difficult to establish a synchronizability measure on two clusters based on Eq. S1, but we can see the following connection between Eqs. S1 and S2:

[TABLE]

$\Downarrow$

[TABLE]

That is, $\|\delta\bm{X}\|$ and $\|\delta\bm{Y}\|$ going to zero in Eq. S2 is a necessary condition for the synchronization in intertwined clusters. For example, if $\|\delta\bm{X}\|$ does not vanish in Eq. S2, then $\|\delta\bm{Y}\|$ must be away from zero in order for $\|\delta\bm{X}\|\rightarrow 0$ in Eq. S1. This connection between Eqs. S1 and S2 implies that we can promote synchronization in the intertwined clusters by optimizing each of the two clusters independently, using the same method originally developed for nonintertwined clusters. In particular, such optimization is guaranteed to expand the region in parameter space satisfying the necessary condition in Eq. S4 (i.e., the condition in Eq. 7 in the main text). Inside this expanded region, one is likely to observe structural AISync, as experimentally demonstrated in the main text. It is worth mentioning that the same argument still holds when more than two clusters are intertwined.

S3 Stability analysis of nondiagonalizable clusters

When dealing with directed networks, one must be aware of the possibility of nondiagonalizable coupling matrices, which can be the case even for symmetric networks godsil1982eigenvalues . Here, we present details of how the analysis in the manuscript also applies to nondiagonalizable networks. To demonstrate that, our key observation is that the treatment of nondiagonalizable networks in Refs. nishikawa2006synchronization ; nishikawa2006maximum can be generalized to the case of a cluster subnetwork in which each oscillator receives a common input from the rest of the network.

We start from the variational equation of the system in the form of Eq. 4 in the main text,

[TABLE]

but this time we lift the assumption that the matrix $\widehat{\bm{L}}$ is diagonalizable. For such systems, in general we can not find $m$ independent eigenvectors for $\widehat{\bm{L}}$. Nevertheless, this matrix can always be transformed into a Jordan canonical form through a similarity transformation defined by an invertible matrix $\bm{P}$, such that

[TABLE]

where $\widehat{v}_{j+1}$ is the eigenvalue of $\widehat{\bm{L}}$ corresponding to the Jordan block $\bm{B}_{j}$, and the matrix entries not shown are zero. The eigenvalues are numbered from $2$ to $q+1$ for consistency with the eigenvalue notation in the main text, and are thus ordered as in the rest of the paper but now without relabeling the (identical) eigenvalues associated with the same Jordan block (which is why $q+1<m$ in the nondiagonalizable case). The special case in which $\widehat{\bm{L}}$ is diagonalizable is also included in this transformation, and it merely corresponds to the case in which all Jordan blocks are $1\times 1$.

Equation S5 can now be decoupled into $q+1$ independent equations accounting for the Jordan blocks. The central difference between the case of an isolated network, as considered in Refs. nishikawa2006synchronization ; nishikawa2006maximum , and the cluster subnetworks considered here is the entry $B_{11}=\widetilde{\mu}$, which is zero for isolated networks. However, this term corresponds to a perturbation mode parallel to the cluster synchronization manifold and hence has no influence on the stability of the synchronization state. (The input connections from the rest of the network to the cluster also impact the synchronization state $\bm{s}_{I}$ and shift the eigenvalues $\widehat{v}_{j}$, but those are not material differences since the same also occurs in the diagonalizable case.) Thus, to analyze the transverse modes, we focus on the $q$ block-decoupled equations associated with the Jordan blocks $\bm{B}_{1},\cdots,\bm{B}_{q}$:

[TABLE]

Assuming that $\bm{B}_{j}$ is $k\times k$, the corresponding equation can be written explicitly for each mode as

[TABLE]

Starting from the first equation in Eq. S8, we notice that $\bm{\eta}^{(j)}_{1}$ does not depend on any other $\bm{\eta}^{(j)}_{i}$ and its equation is exactly the master stability equation [Eq. 5 in the main text]. If Eq. 5 is stable for $\widehat{v}_{j+1}$, then $\bm{\eta}^{(j)}_{1}$ converges to zero exponentially. Turning to the second equation in Eq. S8, we can see that the influence of $\bm{\eta}^{(j)}_{1}$ on $\bm{\eta}^{(j)}_{2}$ vanishes and $\bm{\eta}^{(j)}_{2}$ will also approach zero as $t\rightarrow\infty$ (under the reasonable assumption that $J\bm{h}(\bm{s}_{I})$ is bounded). Applying the same argument iteratively, it follows that the stability of Eq. S8 is entirely determined by the stability of Eq. 5 for the eigenvalue $\widehat{v}_{j+1}$ (with $\bm{\eta}_{i}$ demoted by $\bm{\eta}^{(j)}_{1}$). The same applies for all $j$ and leads to the conclusion that, even if $\widehat{\bm{L}}$ is nondiagonalizable, the condition for the cluster synchronous state to be stable is that $\Lambda(\sigma\widehat{v}_{j+1})<0$ for $j=1,\dots,q$, where $\Lambda$ is the largest Lyapunov exponent of Eq. 5 and $\widehat{v}_{2},\cdots,\widehat{v}_{q+1}$ represent the eigenvalues associated with the transverse modes. Therefore, our analysis of synchronizability presented in the main text (including the use of the eigenratio $R$) applies equally well to nondiagonalizable networks.

S4 Improving synchronizability through minimal link rewiring

In this section, we consider the optimization of symmetry clusters by rewiring a small number of links. One rewiring consists of removing an existing link and adding a different link not yet present in the cluster. Specifically, we developed an algorithm to optimize synchronizability by rewiring intra-cluster connections SA , which preserves the nonintertwined nature of the clusters. This allows us to investigate how many directional links need to be rewired to reduce the eigenratio gap $R-1$ by half. Figure S2 summarizes results for all connected symmetry clusters that are undirected of sizes between $m=8$ and $17$, where the rewiring percentage $p=h/\ell$ is the ratio between the minimal number of link rewiring $h$ that halves $R-1$ and the total number $\ell$ of internal directed links of the cluster. Figure S2(a) shows that on average only about $14\%$ of the links need to be rewired to significantly improve synchronizability of symmetry clusters, and it is largely size independent. Our algorithm also works for link addition and link removal. In the case of link addition, link density needs to increase by about $20\%$ on average to reduce the eigenratio gap to half; for link removal, about $40\%$ of the links need to be removed to achieve the same effect.

Figure S2(b) shows the rewiring percentage $p$ as function of the eigenratio $R$ and link density $D=\frac{\ell}{m(m-1)}$, where each data point represents one symmetry cluster. It is clear that clusters that are small in both $D$ and $R$ require the highest percentage of links to be rewired in order to significantly reduce the eigenratio gap. This confirms the intuition that if a network achieves a small eigenratio with a relatively small number of links, then its organization is efficient and its synchronizability is relatively hard to improve. Conversely, a dense non-optimal network or a network with a relatively large eigenratio is easy to optimize with a small number of link modifications.

S5 Application of the minimal-rewiring algorithm to global synchronization

In this section, we apply the algorithm from the last section to a case in which the full network is symmetric and we seek to optimize global synchronization. In Fig. S3 we study a $16$-node symmetric network and show explicitly through our experiments that it becomes more synchronizable with less symmetry. In the original network [Fig. S3(a)], all nodes play exactly the same structural role. After seven directional link rewiring [marked in Fig. S3(b)], the symmetry of the network is largely broken and almost all nodes are now structurally different: the original 16-node symmetry cluster is reduced to $14$ single-node clusters and only $2$ nodes occupying symmetric positions. The eigenratio, however, reduces from $R=4.62$ to $R=2.80$ and thus improves significantly.

The experimental results are presented in Fig. S3(c), where we show the average synchronization error as a function of the coupling strength for both networks. The experimental data clearly demonstrates that synchronization is only achieved for the network with reduced symmetry. The experimental result is consistent with the MTLE determined from numerical calculations of the variational equation of the model in Eq. 6 [color-coded curves in Fig. S3(c)]. Indeed, for values of $\sigma$ close to the boundary of linear stability, synchronization is not observed in experiments due to noise in the ADC hart2017experiments , but synchronization is consistently observed once the MTLE becomes sufficiently negative.

S6 Prevalence of structural AISync

To further demonstrate that the phenomenon we describe is common across different nodal dynamics and network structure, we present two additional examples. For both examples we consider a random network with five symmetry clusters, as shown in Fig. S4a(a). Within this network, we focus on the highlighted symmetry cluster (magenta nodes), which in isolation corresponds to the second symmetry cluster in Table 1, and we contrast its synchronizability with that of the non-symmetric cluster generated by removing a subset of its links (red links).

We first consider this system when the nodes are equipped with dynamics of a Bernoulli map,

[TABLE]

which, for being piecewise linear and one dimensional, is arguably one of the simplest possible nodal dynamics that one can consider in an oscillator network. Despite its simplicity, this system exhibits a rich stability diagram in the $r\times\sigma$ parameter space, including a wide region in which synchronization is stable for the non-symmetric cluster but unstable for the symmetric one, as shown in Fig. S4a(b). For $r\geq 5$, in particular, synchronization in the symmetric cluster is unstable for any coupling strength $\sigma$. Topological control is particularly valuable in this case as it allows for stability that would be impossible by merely adjusting $\sigma$ in the original cluster.

As an illustration of higher dimensional nonlinear nodal dynamics, we also consider the system in Fig. S4a(a) when equipped with the dynamics of a Hénon map,

[TABLE]

where the variables $x_{i}$ and $y_{i}$ are defined on a torus and limited to $[-2,2]$; the coupling between oscillators are through the $y_{i}$ variables. As shown in Fig. S4a(c), fixing $b=0.3$ and calculating the stability diagram in the $\sigma\times a$ parameter space, once again we identify a wide region in which the non-symmetric cluster exhibits stable synchronization whereas the symmetric one does not.

As illustrated by these and other systems we have studied in detail, in general a significant portion of the parameter space is occupied by a region in which synchronization is not stable for the symmetric cluster but it becomes stable when the structure of the cluster is optimized, which in turn goes in tandem with breaking its symmetry under the given constraints. These examples also further illustrate the excellent agreement between direct simulations and theoretical predictions observed throughout.

S7 Robustness of structural AISync

In order to demonstrate the robustness of structural AISync, we perform direct simulations of three different systems in the presence of Gaussian noise or random oscillator heterogeneity; the results are summarized in Fig. S5a. The three systems include the Bernoulli maps and Hénon maps studied in Sec. S4, as well as the optoelectronic oscillators from the main text.

In Fig. S5a(a), we fix the parameter of the Bernoulli map to be $r=5$ and slowly increase the coupling strength $\sigma$ from 0.3 to 1. For the trajectories in the upper left panel, a random mismatch of magnitude $\xi=10^{-3}$ is introduced to the oscillator parameter $r$; For the trajectories in the middle left panel, the oscillators are subject to Gaussian noise with zero mean and standard deviation equal to $\xi=10^{-3}$ (approximately the noise intensity in the experimental system). Despite the noise and oscillator heterogeneity, the synchronization error $\Delta$ match well with the prediction based on the MTLE calculations shown in the lower left panel. We investigate the dependence of the time-averaged synchronization error $\langle\Delta\rangle$ on the magnitude $\xi$ of noise/mismatch in the right panel, where $\sigma$ is fixed at 0.85 (corresponding to the dashed line on the left).

The same analysis is performed for the Hénon maps in Fig. S5a(b) and for the optoelectronic oscillators in Fig. S5a(c). For the Hénon maps, mismatch is introduced in the parameter $b$, whose homogeneous value is set to $b=0.3$, for coupling strength fixed at $\sigma=0.5$. For the optoelectronic oscillators, mismatch is introduced in the parameter $\beta$, whose homogeneous value is set to $\beta=6$. It can be seen that in all three cases structural AISync is robust against both noise and oscillator heterogeneity.