Asymmetry-Induced Synchronization in Oscillator Networks

TL;DR

This paper introduces a general method for creating asymmetry-induced synchronization in oscillator networks, showing that breaking symmetry by making oscillators nonidentical often stabilizes synchronization in multilayer networks.

Contribution

It presents the first general scheme for AISync systems and demonstrates that AISync is prevalent in multilayer physical networks.

Findings

AISync is common in multilayer networks.

Breaking symmetry stabilizes synchronization.

AISync facilitates cluster synchronization.

Abstract

A scenario has recently been reported in which in order to stabilize complete synchronization of an oscillator network---a symmetric state---the symmetry of the system itself has to be broken by making the oscillators nonidentical. But how often does such behavior---which we term asymmetry-induced synchronization (AISync)---occur in oscillator networks? Here we present the first general scheme for constructing AISync systems and demonstrate that this behavior is the norm rather than the exception in a wide class of physical systems that can be seen as multilayer networks. Since a symmetric network in complete synchrony is the basic building block of cluster synchronization in more general networks, AISync should be common also in facilitating cluster synchronization by breaking the symmetry of the cluster subnetworks.

Figures

Tables

Table 1. Table 1: Number of isomorphically distinct AISync -favoring networks, listed for N = 3 , 4 , 5 𝑁 3 4 5 N=3,4,5 nodes and L = 2 𝐿 2 L=2 layers (with a = b = 1 𝑎 𝑏 1 a=b=1 to enable counting). The numbers are given for both binary (B) and quaternary (Q) choices of internal sublink configurations, as well as for different AISync strength (as measured by r 𝑟 r defined in the text). The network diagrams encode all possible symmetric networks of a given size.

	$N=3$	$N=4$	$N=5$
symmetric networks
Q (optimal)	9	14	21
Q ($r>0.2$)	11	81	254
Q ($r>0.05$)	29	318	2154
B ($r>0.2$)	11	101	204
B ($r>0.05$)	31	400	2406

Table 2. Table 2: Diagrams of symmetric networks with N = 6 𝑁 6 N=6 , 7 7 7 , and 8 8 8 nodes.

$N=6$
$N=7$
$N=8$	and 11 others

Table 3. Table S1: The 9 9 9 optimal heterogeneous systems with N = 3 𝑁 3 N=3 and L = 2 𝐿 2 L=2 .

symmetric network	optimal heterogeneous systems

Table 4. Table S2: The 14 14 14 optimal heterogeneous systems with N = 4 𝑁 4 N=4 and L = 2 𝐿 2 L=2 .

symmetric network	optimal heterogeneous systems

Equations

$${\dot{\mathbf{X}}}_i={\mathbf{F}}_i({\mathbf{X}}_i)+\sum _{\alpha =1}^K\sum _{\begin{array}{c}{i}^{\prime }=1\\ i^{\prime }\ne i\end{array}}^N{A}_{i{i}^{\prime }}^{(\alpha )}{\mathbf{H}}^{(\alpha )}({\mathbf{X}}_i,{\mathbf{X}}_{i^{\prime }}),$$

$$\dot{\bm{x}}^{(i)}_{\ell}=\bm{f}(\bm{x}^{(i)}_{\ell})+\sum_{i^{\prime}=1}^{N}\sum_{\ell^{\prime}=1}^{L}\widetilde{A}^{(ii^{\prime})}_{\ell\ell^{\prime}}\bigl{[}\bm{h}(\bm{x}^{(i^{\prime})}_{\ell^{\prime}})-\bm{h}(\bm{x}^{(i)}_{\ell})\bigr{]},$$

$${\dot{\mathbf{x}}}_i=D\mathbf{f}{\mathbf{x}}_i-\sum _j{\tilde{L}}_{ij}D\mathbf{h}{\mathbf{x}}_j,$$

$$D\mathbf{f}=\left(\begin{array}{cccc}-2& 2& -1& 2\\ -1& 1& 0& 0\\ 0& 0& -3& 4\\ 0& 0& -1& 1\end{array}\right)D\mathbf{h}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 1& 0& -1\\ -1& 0& 1& 0\end{array}\right).$$

$$x(x+3)-y^2-(2x+3{)}^2{y}^2>0,$$

$$x_{i}^{t+1}=\biggl{[}f(x_{i}^{t})-\sum_{j}\widetilde{L}_{ij}f(x_{j}^{t})+\delta\biggr{]}\mod 2\pi,$$

$${\mathbf{X}}_i:=\left(\begin{array}{c}{\mathbf{x}}_1^{(i)}\\ ⋮\\ {\mathbf{x}}_L^{(i)}\end{array}\right),{\mathbf{F}}_i:=\left(\begin{array}{c}{\mathbf{F}}_1^{(i)}\\ ⋮\\ {\mathbf{F}}_L^{(i)}\end{array}\right),{\mathbf{H}}^{(\alpha )}:=\left(\begin{array}{c}{\mathbf{H}}_1^{(\alpha )}\\ ⋮\\ {\mathbf{H}}_L^{(\alpha )}\end{array}\right),$$

$${\mathbf{F}}_{\ell }^{(i)}({\mathbf{X}}_i):=\mathbf{f}({\mathbf{x}}_{\ell }^{(i)})+\sum _{{\ell }^{\prime }=1}^L{\tilde{A}}_{\ell {\ell }^{\prime }}^{(ii)}[\mathbf{h}({\mathbf{x}}_{{\ell }^{\prime }}^{(i)})-\mathbf{h}({\mathbf{x}}_{\ell }^{(i)})],$$

$${\mathbf{H}}_{\ell }^{(\alpha )}({\mathbf{X}}_i,{\mathbf{X}}_{i^{\prime }}):=\sum _{{\ell }^{\prime }=1}^L{B}_{\ell {\ell }^{\prime }}^{(\alpha )}[\mathbf{h}({\mathbf{x}}_{{\ell }^{\prime }}^{(i^{\prime })})-\mathbf{h}({\mathbf{x}}_{\ell }^{(i)})],$$

$${\mathbf{H}}^{(3)}({\mathbf{X}}_1,{\mathbf{X}}_4)=\left(\begin{array}{c}\mathbf{h}({\mathbf{s}}_2)-\mathbf{h}({\mathbf{s}}_1)\\ 0\end{array}\right)\not\equiv 0.$$

$${\dot{\mathbf{x}}}_j=\mathbf{f}({\mathbf{x}}_j)+\sum _{j^{\prime }=1}^n{\tilde{A}}_{j{j}^{\prime }}[\mathbf{h}({\mathbf{x}}_{j^{\prime }})-\mathbf{h}({\mathbf{x}}_j)],$$

$$\dot{\mathbf{\xi }}=[D\mathbf{f}(\mathbf{s})-\lambda D\mathbf{h}(\mathbf{s})]\mathbf{\xi },$$

$$e^2:=\frac{1}{N}\sum _{i=1}^N\sum _{\ell =1}^L||{\mathbf{x}}_{\ell }^{(i)}-{\overline{\mathbf{x}}}_{\ell }|{|}^2,{\overline{\mathbf{x}}}_{\ell }:=\frac{1}{N}\sum _{i=1}^N{\mathbf{x}}_{\ell }^{(i)}.$$

$$\begin{array}{r}\mathbf{f}(\mathbf{x})=\left(\begin{array}{c}\gamma (x_2-x_1)\\ x_1(\rho -x_3)-x_2\\ x_1{x}_2-\beta x_3\end{array}\right),\\ \mathbf{h}(\mathbf{x})=\left(\begin{array}{c}{x}_2\\ 0\\ 0\end{array}\right),\mathbf{x}=\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right)\end{array}$$

$$\left(\begin{array}{c}{\dot{\xi }}_1\\ {\dot{\xi }}_2\\ {\dot{\xi }}_3\end{array}\right)=\left(\begin{array}{ccc}-\gamma & \gamma -\lambda & 0\\ \rho -s_3& -1& -s_1\\ s_2& s_1& -\beta \end{array}\right)\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\\ {\xi }_3\end{array}\right),$$

$$\begin{array}{rl}{\dot{s}}_1& =\gamma (s_2-s_1),\\ {\dot{s}}_2& =s_1(\rho -s_3)-s_2,\\ {\dot{s}}_3& =s_1{s}_2-\beta s_3.\end{array}$$

$${\tilde{A}}_{j{j}^{\prime }}^c:=(1-{\tilde{A}}_{j{j}^{\prime }})(1-{\delta }_{j{j}^{\prime }}).$$

$${\sigma }^c=\frac{\tilde{m}\sigma }{n(n-1)-\tilde{m}},$$

$$\frac{{\sigma }_1}{{\sigma }_2}=\frac{{\sigma }_1^c}{{\sigma }_2^c}$$

Full text

I Introduction

A common assumption in the field of network dynamics is that homogeneity in the local dynamics Restrepo et al. (2004); Sun et al. (2009) and interaction network Nishikawa et al. (2003); Denker et al. (2004); Donetti et al. (2005)—or in the combination of both Motter et al. (2005); Zhou and Kurths (2006)—can facilitate complete synchronization. It has been recently shown, however, that structural heterogeneity in networks of identical oscillators Nishikawa and Motter (2010) or oscillator heterogeneity in structurally symmetric networks Nishikawa and Motter (2016) can stabilize otherwise unstable synchronous states, thus effectively breaking the symmetry of a system to stabilize a symmetric state. These scenarios, which we refer to as asymmetry-induced synchronization (AISync), can be interpreted as the converse of symmetry breaking, and hence as a converse of chimera states Kuramoto and Battogtokh (2002); Abrams and Strogatz (2004). Perhaps the most striking and the strongest form of AISync is the one in which oscillators coupled in a symmetric network (i.e., each oscillator plays exactly the same structural role) can converge to identical dynamics only when they themselves are nonidentical; this has been demonstrated, however, exclusively for rotationally symmetric networks and one type of periodic oscillators Nishikawa and Motter (2016). Whether such AISync behavior can be shown to be common among systems with other symmetric network structures and oscillator dynamics, including experimentally testable ones, has been an open question.

In this article we introduce and analyze a broad class of AISync systems that can have general symmetric network structure with multiple link types and general oscillator dynamics (which can be chaotic, periodic, continuous-time, discrete-time, etc.). This in particular includes physical systems previously used in network synchronization experiments, thus providing a recipe for future empirical studies. For this class, we demonstrate that AISync is indeed common and provide a full characterization of those networks that support AISync behavior, showing that the fraction of such networks is significant over a range of network sizes and link densities.

II Definition of AISync

To formulate a precise definition of AISync, we consider networks of $N$ (not necessarily identical) oscillators coupled through $K$ different types of interactions. The network dynamics is described by

[TABLE]

where $\bm{X}_{i}=\bm{X}_{i}(t)$ is the $M$-dimensional state vector of node $i$, the function $\bm{F}_{i}$ governs the intrinsic dynamics of node $i$, the adjacency matrix $A^{(\alpha)}=(A^{(\alpha)}_{ii^{\prime}})$ represents the topology of interactions through links of type $\alpha$, and $\bm{H}^{(\alpha)}$ is the coupling function associated with the link type $\alpha$. A completely synchronous state of the network is defined by $\bm{X}_{1}(t)=\bm{X}_{2}(t)=\cdots=\bm{X}_{N}(t)$.

To isolate the effect of breaking the homogeneity of oscillators, we consider adjacency matrices $A^{(\alpha)}$ that together represent a symmetric network, defined as a network in which every node can be mapped to any other node by some permutation of nodes without changing any $A^{(\alpha)}$. Thus, the set of links of any given type must couple every node identically (see Fig. 1(a) for an example). The rationale for using symmetric network structures here is to ensure that any stabilization of complete synchronization by oscillator heterogeneity is due to the reduced system symmetry (as required for AISync) and not due to having network heterogeneity and oscillator heterogeneity compensating each other, which may not break the system symmetry.

When restricted to undirected networks with a single link type, our definition of symmetric networks yields the class of vertex-transitive graphs from graph theory Biggs (1993). This rich class encompasses Cayley graphs (defined as a network of relations between elements of a finite group; Appendix A) and circulant graphs (defined as a network whose nodes can be arranged in a ring so that the network is invariant under rotations), which have previously been used to study chimera states Panaggio and Abrams (2015). Enumerating all vertex-transitive graphs of a given size $N$ becomes challenging as $N$ grows and has so far been completed only for $N<32$ oei . The symmetric networks we consider here generalizes vertex-transitive graphs to the even richer class of networks that can be directed and include multiple link types.

Given a symmetric network structure, the system in Eq. (1) exhibits AISync if it satisfies the following two conditions: (C1) there are no asymptotically stable synchronous states for any homogeneous system (i.e., with $\bm{F}_{1}=\cdots=\bm{F}_{N}$), and (C2) there is an asymptotically heterogeneous system (i.e., with $\bm{F}_{i}\neq\bm{F}_{i^{\prime}}$, for some $i\neq i^{\prime}$) for which a stable synchronous state exists. A challenge in establishing AISync is that the form of Eq. (1) does not guarantee the existence of a completely synchronous state. Another challenge concerns the stability analysis of such a state, since Eq. (1) is beyond the framework normally used in the master stability function (MSF) approach and its generalizations currently available Sun et al. (2009); Dahms et al. (2012); Pecora et al. (2014); del Genio et al. (2016): oscillators can be nonidentical (different $\bm{F}_{i}$), and the network can host $K>1$ types of directed interactions.

III Multilayer systems considered

To overcome these challenges, below we propose a multilayer construction that defines a large, general subclass of systems within the class given by Eq. (1). We show that any system in this subclass is guaranteed to have a synchronous state, and the stability of that state can be analyzed by applying the MSF framework to the flattened representation of the system. The MSF approach decouples the oscillator dynamics from the network structure, which enables us to draw conclusions about AISync for general oscillator dynamics.

In our multilayer system, each node is composed of $L$ identical subnodes, belonging to $L$ different layers and connected by a set of internal sublinks. The pattern of these internal sublinks is thus part of the node’s properties and determines the heterogeneity across nodes. For a pair of connected nodes, the type of the connecting link is determined by the pattern of external sublinks between the subnodes of these two nodes. This construction yields a multilayer network Gao et al. (2012); Gomez et al. (2013); De Domenico et al. (2013); Boccaletti et al. (2014); De Domenico et al. (2016) of subnodes and sublinks with $L$ layers; see Fig. 1(b) for an $L=2$ example. Note that in general there is more than one possible multilayer network for a given symmetric network. Networks with such layered structure have been used extensively as realistic models of various natural and man-made systems. The class of systems just defined is broader than most classes of systems used in previous studies of synchronization on multilayer networks Zhang et al. (2015); Gambuzza et al. (2015), since the links between two different layers are not constrained to be one-to-one. The underlying hierarchical organization, in which each node decomposes into interacting subnodes, is shared by many physical systems, such as the multi-processor nodes in modern supercomputers.

Coupling the dynamics of subnodes diffusively in this network, the multilayer system can be described at the subnode level as

[TABLE]

where $\bm{x}^{(i)}_{\ell}=\bm{x}^{(i)}_{\ell}(t)$ is the $m$-dimensional state vector for subnode $\ell$ (i.e., in layer $\ell$) of node $i$, the function $\bm{f}$ determines the dynamics of every isolated subnode, and $\bm{h}$ is the coupling function common to all sublinks. Here, for all links of a given type between different nodes, the corresponding coupling matrix $\widetilde{A}^{(ii^{\prime})}:=(\widetilde{A}^{(ii^{\prime})}_{\ell\ell^{\prime}})$, $i\neq i^{\prime}$, is the same and encodes the subnode connection pattern for that link type. In contrast, the subnode connection pattern within each node $i$ is encoded in the matrix $F^{(i)}:=(\widetilde{A}^{(ii)}_{\ell\ell^{\prime}})$. Since the subnode-to-subnode interactions are diffusive, the synchronous state given by $\bm{x}^{(i)}_{\ell}(t)=\bm{s}(t)$, $\forall i,\ell$ with $\dot{\bm{s}}=\bm{f}\left(\bm{s}\right)$ is guaranteed to exist. Note that the diffusive coupling among subnodes do not necessarily imply that the node-to-node interactions are diffusive, as intralayer synchronization of the form $\bm{x}^{(i)}_{\ell}=\bm{s}_{\ell}$ among subnodes is also valid as a state of complete synchronization among all nodes. The interactions among nodes do not vanish in this case due to the existence of external sublink connections among different layers.

To summarize, Eq. (2) describes a general class of multilayer models of symmetric networks that admit a state corresponding to complete synchronization, $\bm{X}_{i}(t)=\bm{S}(t)$, $\forall i$, when written in the form of Eq. (1) (see Appendix B for details).

IV Establishing AISync

To facilitate the stability analysis required to establish AISync, we flatten the multilayer network representation into a single layer (see Fig. 1(c) for an example). We use $\widetilde{A}=(\widetilde{A}_{jj^{\prime}})$ to denote the adjacency matrix that encodes the structure of the resulting monolayer network (see Fig. 1(d) for an example). This matrix has a block structure in which the matrices $F^{(i)}$ appearing on the diagonal blocks characterize node properties, while $\widetilde{A}^{(ii^{\prime})}$ appearing on the off-diagonal blocks reflect the link types. Since subnodes and sublinks are identical, we can directly apply the MSF analysis Pecora and Carroll (1998) to the monolayer network and obtain the stability function $\psi(\lambda)$ (see Appendix C for details). The maximum transverse Lyapunov exponent (MTLE) of the synchronous state is then computed as $\Psi:=\max_{2\leq j\leq n}\psi(\lambda_{j})$, where $\lambda_{1},\ldots,\lambda_{n}$ are the eigenvalues of the corresponding Laplacian matrix $\widetilde{L}:=(\widetilde{L}_{jj^{\prime}})$, defined as $\widetilde{L}_{jj^{\prime}}:=\delta_{jj^{\prime}}\sum_{k=1}^{n}\widetilde{A}_{jk}-\widetilde{A}_{jj^{\prime}}$, where $\delta_{jj^{\prime}}$ is the Kronecker delta function. Here, $\lambda_{1}$ is the identically zero eigenvalue, which is excluded in the definition of $\Psi$ for corresponding to a mode of perturbation that does not affect synchronization. Thus, the synchronous state is asymptotically stable if $\Psi<0$ and unstable if $\Psi>0$.

To establish AISync for our multilayer system, we first verify that all homogeneous systems have $\Psi>0$ (i.e., synchronous state $\bm{x}^{(i)}_{\ell}=\bm{s}$, $\forall i,\ell$, is unstable), and check numerically that all other synchronous states $\bm{x}^{(i)}_{\ell}=\bm{s}_{\ell}$, $\forall i,\ell$, are also unstable. This establishes condition (C1). We then find a heterogeneous system with $\Psi<0$, which establishes condition (C2). This procedure is detailed in Appendix D.

In the case of linear $\bm{f}$ and $\bm{h}$, which is widely used to study consensus dynamics and encompasses a variety of nontrivial stability regions Li et al. (2010), the problem of verifying AISync is fully solvable. To see this, we first note that in this case the stability function $\psi(\lambda)$ determines the (common) stability of all completely synchronous states of the form $\bm{x}^{(i)}_{\ell}=\bm{s}_{\ell}$, $\forall i,\ell$, where the subnode states $\bm{s}_{\ell}$ can in general be different for different $\ell$. Next, for a given (homogeneous or heterogeneous) system, we sort its Laplacian eigenvalues into two groups: $\lambda_{1},\ldots,\lambda_{j^{*}}$, corresponding only to those perturbations parallel to the synchronization manifold, and $\lambda_{j^{*}+1},\ldots,\lambda_{n}$, corresponding to perturbations that are transverse to the manifold and thus destroy synchronization. The stability (of all completely synchronous states) is then determined by $\Psi^{\prime}:=\max_{j^{*}<j\leq n}\psi(\lambda_{j})$, noting that both $j^{*}$ and $\lambda_{j}$ generally depend on the network structure. This leads to the following solution for the AISync conditions: $\Psi^{\prime}\geq 0$ for all homogeneous systems and $\Psi^{\prime}<0$ for some heterogeneous system (where we include $\Psi^{\prime}=0$ in the first condition because $\Psi^{\prime}=0$ for linear system would exclude asymptotically stable synchronization).

V Examples of AISync

V.1 Consensus dynamics

Here we establish AISync for the system with the symmetric network structure shown in Fig. 2, in which the subnodes follow the consensus dynamics used in Ref. Li et al. (2010):

[TABLE]

where

[TABLE]

This leads to the stability region $\psi(\lambda)<0$ shown in Fig. 2(c), defined by

[TABLE]

where $x$ and $y$ denote the real and imaginary parts of $\lambda$, respectively.

For $L=2$ there are only two possible homogeneous systems, associated with the two possible directions of the internal sublink in each node. The homogeneous system in Fig. 2(a) has Laplacian eigenvalues $\lambda_{1}=0$, $\lambda_{2}=2$, $\lambda_{3,4}\approx 0.5\pm 0.866i$, and $\lambda_{5,6}\approx 1.5\pm 0.866i$, where $\lambda_{1}$ and $\lambda_{2}$ correspond to the perturbations parallel to the synchronization manifold and $\lambda_{3},\ldots,\lambda_{6}$ correspond to those in the transverse directions (i.e., $j^{*}=2$). Since $\psi(\lambda_{j})>0$ for $j=3,4,5,6$ [i.e., all these $\lambda_{j}$’s fall outside the stability region defined by Eq. (5), as indicated by the red squares in Fig. 2(c)], we have $\Psi^{\prime}=\max_{2<j\leq 6}\psi(\lambda_{j})>0$. The other homogeneous system is not synchronizable since all the single-prime subnodes have no incoming sublink. In contrast, for the heterogeneous system in Fig. 2(b), the Laplacian eigenvalues are $\lambda_{1}=0$, $\lambda_{j}=1$ for $1<j\leq 5$, and $\lambda_{6}=2$ (i.e., $j^{*}=1$ in this case). As shown by the blue dots in Fig. 2(c), we have $\Psi^{\prime}=\max_{1<j\leq 6}\psi(\lambda_{j})<0$ for this heterogeneous system. We thus see that $\Psi^{\prime}\geq 0$ (i.e., the synchronous state is not asymptotically stable) for both homogeneous systems and $\Psi^{\prime}<0$ (i.e., the synchronous state is asymptotically stable) for a heterogeneous system, establishing AISync: the agents can reach consensus only when some of them are different from the others.

V.2 Coupled Lorenz oscillators

An example of nonlinear system exhibiting AISync is shown in Fig. 3. The network structure is symmetric and has $N=3$ nodes and $K=2$ types of links representing sublink patterns in the clockwise and counterclockwise directions in Fig. 3(a). Each node consists of $L=2$ subnodes, each of which is a chaotic Lorenz oscillator. The two subnodes are connected by a sublink, the direction of which determines the node type. This gives rise to two node types, and there are four possible distinct combinations of node types for the network—two homogeneous and two heterogeneous. The system has two parameters, $a$ and $b$, representing the coupling strength of internal and external sublinks, respectively. We seek to determine for which values of $a$ and $b$ the system exhibits AISync.

In Fig. 3(b), we show $\Psi_{=}$ (red) and $\Psi_{\neq}$ (blue) as functions of $a$ and $b$, where $\Psi_{=}$ ($\Psi_{\neq}$) are defined to be the smaller value of $\Psi$ between the two possible homogeneous (heterogeneous) systems. In the region shaded purple (where $\Psi_{\neq}>0$ and $\Psi_{=}<0$), the synchronous state $\bm{x}^{(i)}_{\ell}(t)=\bm{s}(t)$, $\forall i,\ell$ is stable for at least one of the heterogeneous systems, but unstable for both homogeneous systems. We further verify in this region that the other possible forms of synchronous states, $\bm{x}^{(i)}_{\ell}(t)=\bm{s}_{\ell}(t)$, $\forall i,\ell$, are unstable for both homogeneous systems (through extensive numerical simulation—see Appendix E for details). This establishes conditions (C1) and (C2), thus confirming that the system exhibits AISync in the purple region.

The AISync behavior of the system for a specific combination of $a$ and $b$ is illustrated by the sample trajectory in Fig. 3(c), which diverges from synchrony while the nodes are kept homogeneous, but re-synchronizes spontaneously after the nodes are made heterogeneous.

V.3 Coupled electro-optic systems

We now present an experimentally testable AISync system based on the discrete-time model of the electro-optic system implemented in Refs. Hagerstrom et al. (2012); Pecora et al. (2014) and given by

[TABLE]

where $f(x)=\beta I(x)$ determines the isolated subnode dynamics and also serves as the coupling function. Here, $I(x)=(1-\cos x)/2$ is the normalized optical intensity, $\beta=1.7\pi$ is the self-feedback strength, $\delta=0.2$ is the offset introduced to suppress the trivial solution $x_{i}=0$, and $\widetilde{L}_{ij}$ is the weighted graph Laplacian [weights controlled by parameters $a$ and $b$, as shown in Fig. 4(a–c)].

Figure 4 shows an example of AISync using these electro-optic maps as subnodes. The internal connections are chosen from the quaternary set (no sublink, one directed sublink in either direction, and directed sublinks in both directions). When the same choice is made for all internal connections, this leads to four different homogeneous systems, but two of them have $\lambda_{2}=0$ (not synchronizable), leaving only two homogeneous systems to consider [Figs. 4(b) and (c)]. For comparison, we take the heterogeneous system in Fig. 4(a), which forms a directed chain network in its monolayer representation. Each of the three systems [Figs. 4(a–c)] has a companion plot showing under what parameters the nodes are synchronizable [Figs. 4(e–g)]. In the latter panels, each pixel is generated from $24$ independent simulations run from random initial conditions. The pixels are then color-coded according to how many times a fully synchronized state was reached after 2500 iterations (“sync-only”: 24 times; “desync-only”: 0 times; “multistability”: all other cases). Here we consider a trajectory to be fully synchronized if the synchronization error $e$ defined in Eq. (13) and averaged over the last $100$ iterations falls below $10^{-3}$. It is worth noting that, in this example, when a homogeneous system is synchronizable the synchronous state is always in the form of cluster synchronization among subnodes (those indexed with prime and double prime form two separate synchronized clusters), since complete synchronization among subnodes is always unstable for both homogeneous systems.

Figure 4(d) shows the stability function $\psi(\lambda)$ for the electro-optic subnode dynamics and coupling function, which has a bounded stable region. The lines are quite dense inside the stable region, meaning that the stability landscape is steep there and the function reaches very deep negative values. This is confirmed in Fig. 4(h), where the AISync regions are shaded purple, with the MTLE of the synchronous state for the heterogeneous system shown as red contour lines.

V.4 Other examples

The three AISync systems considered in this section include both linear and nonlinear coupling functions, as well as discrete- and continuous-time dynamics. While they provide AISync examples with directed network structures and bounded stability regions, we also demonstrate AISync for systems with unbounded stability region (Supplemental Material Sec. S1), and for systems with undirected network structure (Supplemental Material Sec. S2).

VI Propensity for AISync

But how often does a network structure support AISync? To systematically address this question, we use the spread $\sigma$ of the eigenvalues of the Laplacian matrix $\widetilde{L}$ for the monolayer network representation, which is a measure of synchronizability Nishikawa and Motter (2010) defined by $\sigma^{2}:=\sum_{j=2}^{n}|\lambda_{j}-\overline{\lambda}|^{2}/[d^{2}(n-1)]$, where $d:=\sum_{j=1}^{n}\widetilde{L}_{jj}/n$ and $\overline{\lambda}:=\sum_{j=2}^{n}\lambda_{j}/(n-1)$. A smaller $\sigma$ indicates higher synchronizability. Given an external sublink structure corresponding to a symmetric network, we compare the minimum spread $\sigma_{=}$ among all systems with homogeneous $F^{(i)}$ to the corresponding minimum $\sigma_{\neq}$ among all systems with strictly heterogeneous $F^{(i)}$. We call the structure AISync-favoring if $\sigma_{\neq}<\sigma_{=}$, which indicates that heterogeneous $F^{(i)}$ can make the system easier to synchronize than any homogeneous $F^{(i)}$. As a measure of how strongly the structure supports AISync, we define $r:=1-\sigma_{\neq}/\sigma_{=}\leq 1$, where $r>0$ indicates an AISync-favoring structure, and $r=1$ implies $\sigma_{\neq}=0$ (i.e., there is a heterogeneous system with optimal synchronizability). For example, the structure in Fig. 3 has $\sigma_{=}\approx 0.56$ and $\sigma_{\neq}\approx 0.33$, and $r\approx 0.41$.

Using this AISync strength $r$, we first enumerate all networks of a given size supporting AISync (Table I). For each $N$, we generate one or more diagrams representing all $N$-node symmetric networks, which are shown in the first row of Table I for $N=3,4,5$. In these diagrams, each color indicates a set of links that, in any given symmetric network, must all exist together and be of the same type or not exist at all (noting that links from different sets can be of the same type). For example, there are three distinct symmetric networks for $N=3$: a directed ring (cyan or black links), an undirected ring (cyan and black links of the same type), and the superposition of two directed rings in opposite directions (cyan and black links of different types; as in Fig. 3). For a given symmetric network derived from these diagrams, we choose the external sublink pattern for each link type from all possible ways of connecting a subnode pair to another. For the internal sublink patterns, we use either the binary or quaternary choices, where each node has one directed sublink (in either direction) in the binary case, while all four possibilities are allowed in the quaternary case. The rest of Table I lists the total numbers of isomorphically distinct external sublink structures with $r>0.05$, $r>0.2$, and (optimal) $r=1$; see Supplementary Tables S1 and S2 for all optimal networks with $N=3$ and $4$, respectively.

Table II extends the first row in Table I, showing the symmetric network diagrams for $N=6,7$, and 8. In each row, the leftmost diagram is the full representation as in Table I, which is decomposed into multiple components (the partial diagrams in the same row) to make them more clearly visible. The partial diagrams with the same background color indicate identical components appearing in multiple rows. Thus, for $N=6$, we have four different diagrams (rows), each with a different combination of components. There is only one diagram for $N=7$, while we show one representative diagram out of twelve in the case of $N=8$.

Figure 5 shows the statistics of AISync-favoring networks. For numerical feasibility, we focus on those systems whose network structure is a directed circulant graph with multiple link types (which covers all symmetric networks if $N$ is a prime number). Sampling uniformly within this class (Appendix F), we observe that significant fraction of external sublink structures are AISync-favoring over a range of external sublink densities [Fig. 5(a)] and network sizes [Fig. 5(b)]. We also observe that sparse and dense structures favor AISync more often than medium-density ones, despite the expectation that the effect of internal sublink heterogeneity would be smaller with higher external sublink density. This phenomenon is further explored in Appendix G by establishing the approximate left-right symmetry in Fig. 5(a).

VII Discussion

Given a symmetric network of identical oscillators, it is instructive to compare our results above in which the symmetry is broken by making the oscillators nonidentical with the alternative scenario in which the symmetry is broken by making the network structure asymmetric. For directed unweighted networks of diffusively-coupled identical oscillators, it can be shown that: 1) with the exception of the complete graphs, all topologies that optimize synchronizability (i.e., those with $\sigma=0$) are asymmetric; 2) any network topology that can be spanned from a node (i.e., $\min_{i\geq 2}\text{Re}(\lambda_{i})>0$) embeds optimally synchronizable subnetworks generated by deleting a subset of links Nishikawa and Motter (2006, 2010). For example, a synchronous state that is not stable for a directed ring network may become stable for a directed chain formed by removing a link. More generally, introducing structural heterogeneity (breaking the symmetry of the network) can stabilize otherwise unstable homogeneous (symmetric) states.

Finally, we note that the defining characteristic of AISync considered here—that preserving the symmetry of a stable state requires breaking the symmetry of the system—can bear analogs in oscillator networks whose structure is not necessarily symmetric. Such a network can always be partitioned into symmetric subnetwork clusters (structurally equivalent subsets of nodes) that are candidates for cluster synchronization Golubitsky and Stewart (2016); Pecora et al. (2014); Nicosia et al. (2013). Synchronization of one of these clusters plays the role of complete synchronization in a symmetric network, which opens the possibility of exploiting AISync to tune cluster synchronization patterns through oscillator heterogeneity in arbitrary complex networks. We hope that our findings, and future theoretical and experimental studies they will stimulate, will significantly advance understanding of the interplay between symmetry and network dynamics.

Acknowledgements.

The authors thank Thomas Wytock and Young Sul Cho for insightful discussions. This work was supported by ARO Grant No. W911NF-15-1-0272.

Appendix A Definition of Cayley graphs

Given a generating set $S$ of a finite group $G$, the Cayley graph associated with $S$ and $G$ is defined as the network in which a node represents an element of $G$ and a directed link from one node $g\in G$ to another $g^{\prime}\in G$ represents the composition of some element $s\in S$ with $g$ that gives $g^{\prime}$ (i.e., $gs=g^{\prime}$). While such a network is generally directed, it will be undirected if the inverse of every element of $S$ belongs to $S$. Choosing $S$ to be a generating set guarantees that the resulting network is (strongly) connected. A generalization to multiple link types can be obtained if we assign different elements of $S$ to different link types.

Appendix B Details on multilayer models

Since Eq. (2) defines a subclass of systems governed by Eq. (1), it can always be written in the form of Eq. (1) for a given network structure specified by $A^{(\alpha)}$. This can be seen by stacking the $m$-dimensional vectors in Eq. (2) and defining appropriate functions as follows:

[TABLE]

where $B^{(\alpha)}_{\ell\ell^{\prime}}$ is defined to be the value of $\widetilde{A}^{(ii^{\prime})}_{\ell\ell^{\prime}}$ when node $i^{\prime}$ is connected to node $i$ by a link of type $\alpha$. Note that these node-to-node interactions are not necessarily diffusive, since we can have $\bm{H}^{(\alpha)}(\bm{X}_{i},\bm{X}_{i^{\prime}})\neq\bm{0}$ even for $\bm{X}_{i}=\bm{X}_{i^{\prime}}$, if $\bm{x}^{(i)}_{\ell}\neq\bm{x}^{(i^{\prime})}_{\ell^{\prime}}$ for some $\ell\neq\ell^{\prime}$ [which in particular means that the coupling term cannot be written in the form ${\bm{H}}^{(\alpha)}({\bm{X}}_{i},{\bm{X}}_{i^{\prime}})=\widetilde{\bm{H}}^{(\alpha)}({\bm{X}}_{i^{\prime}})-\widetilde{\bm{H}}^{(\alpha)}({\bm{X}}_{i})$]. For example, even when nodes $1$ and $4$ are synchronized in the network of Fig. 1, i.e., ${\bm{X}}_{1}={\bm{X}}_{4}=(\bm{s}_{1}(t),\bm{s}_{2}(t))^{T}$, the coupling term corresponding to the link of type $\alpha=3$ between them is in general not identically zero:

[TABLE]

However, since we assume identical dynamics for subnodes and diffusive coupling between subnodes, a synchronous state of Eq. (2) given by $\bm{x}^{(i)}_{\ell}(t)=\bm{s}(t)$, $\forall i,\ell$ with $\dot{\bm{s}}=\bm{f}\left(\bm{s}\right)$ is guaranteed to exist even if $F^{(i)}$’s are heterogeneous. This corresponds to a global synchronous state of Eq. (1) defined by $\bm{X}_{i}=\bm{S}:=(\bm{s},\ldots,\bm{s})^{T}$, $\forall i$, which can be verified by noting that $\bm{H}^{(\alpha)}(\bm{S},\bm{S})=\bm{0}$ and $\bm{F}^{(i)}_{\ell}(\bm{S}):=\bm{f}(\bm{s})$, $\forall i,\ell$.

Appendix C Details on MSF analysis

Equation (2) can be rewritten as a monolayer network by defining a single index for all the $n:=LN$ subnodes, in which node $i$ has subnodes $j=k_{i1},\ldots,k_{iL}$ with $k_{i\ell}:=(i-1)L+\ell$. This leads to the standard form for a (monolayer) diffusively coupled network of oscillators:

[TABLE]

where $\bm{x}_{j}=\bm{x}^{(i)}_{\ell}$ and $\widetilde{A}_{jj^{\prime}}:=\widetilde{A}^{(ii^{\prime})}_{\ell\ell^{\prime}}$ for $j=k_{i\ell}$ and $j^{\prime}=k_{i^{\prime}\ell^{\prime}}$. In the monolayer adjacency matrix $\widetilde{A}=(\widetilde{A}_{jj^{\prime}})$, the matrix $B^{(\alpha)}=(B^{(\alpha)}_{\ell\ell^{\prime}})$ appears as multiple off-diagonal blocks of size $L$, and the arrangement of those blocks within $\widetilde{A}$ matches with the structure of the corresponding adjacency matrix $A^{(\alpha)}$, reflecting the topology of node-to-node interactions through links of type $\alpha$ [see Fig. 1(d) for an example]. This equation allows application of the MSF analysis Pecora and Carroll (1998) because subnodes and sublinks (and the associated coupling functions) are identical. The stability function $\psi(\lambda)$ is defined as the maximum Lyapunov exponent of the reduced variational equation,

[TABLE]

where $\bm{\xi}$ is an $m$-dimensional perturbation vector, $D\bm{f}(\bm{s})$ and $D\bm{h}(\bm{s})$ are the Jacobian of $\bm{f}$ and $\bm{h}$, respectively, at the synchronous state, $\bm{x}_{j}=\bm{s}(t)$, $\forall j$, and $\lambda$ is a complex-valued parameter.

Appendix D Verifying the AISync conditions

Here we describe our scheme for verifying AISync conditions (C1) and (C2) given a symmetric network structure (adjacency matrices $A^{(\alpha)}$), external sublink configurations (matrices $B^{(\alpha)}$), a set $\mathcal{F}$ of possible internal sublink configurations (from which matrices $F^{(i)}$ are chosen), isolated subnode dynamics $\bm{f}$, and sublink coupling function $\bm{h}$. We first obtain the stability function $\psi(\lambda)$:

Compute a trajectory $\bm{s}(t)$ of an isolated subnode by integrating $\dot{\bm{s}}=\bm{f}\left(\bm{s}\right)$, which determines the synchronous state, $\bm{x}^{(i)}_{\ell}(t)=\bm{s}(t)$, $\forall i,\ell$. 2. 2.

Integrate Eq. (12) and calculate its maximum Lyapunov exponent (MLE), which defines $\psi(\lambda)$ for a range of $\lambda$ in the complex plane.

Note that $\psi(\lambda)$ depends only on $\bm{f}$, $\bm{h}$, and $\bm{s}(t)$. For a given symmetric network structure and external sublink configurations, we can compute the stability $\Psi$ of the synchronous state for any combination of $F^{(i)}\in\mathcal{F}$ by calculating and substituting the Laplacian eigenvalues $\lambda_{j}$ into the formula $\Psi=\max_{2\leq j\leq n}\psi(\lambda_{j})$. To establish the AISync property, we verify the following conditions:

• •

(C1)′: For each matrix $F\in\mathcal{F}$, set $F^{(1)}=\cdots=F^{(N)}=F$ (leading to a homogeneous system) and verify $\Psi>0$.

• •

(C2)′: Identify a combination of (heterogeneous) $F^{(i)}\in\mathcal{F}$ for which $\Psi<0$ (e.g., by checking exhaustively or by using a numerical optimization algorithm to minimize $\Psi$ over $F^{(i)}$).

The verification of condition (C1)′ provides strong support for (C1), since the only other possibility for a stable synchronization of all nodes is a state of the form $\bm{x}^{(i)}_{\ell}=\bm{s}_{\ell}(t)$, $\forall i,\ell$, with at least one $\bm{s}_{\ell}(t)$ different from the others (which we find does not exist in many cases, such as the examples in Fig. 3 and in Supplemental Material Sec. S1). To provide additional support for (C1), we directly simulate Eq. (2) from a set of initial conditions and verify that the synchronization error $e$ does not approach zero whenever $F^{(1)}=\cdots=F^{(N)}$, where $e$ is defined as the standard deviation of the node state vectors, or equivalently,

[TABLE]

Here $||\cdot||$ denotes the $2$-norm in the state space of the subnode dynamics, and $e=0$ is achieved if and only if the system is in a synchronous state of the form $\bm{x}^{(i)}_{\ell}=\bm{s}_{\ell}(t)$. To complete our procedure, we verify condition (C2)′, which rigorously establishes (C2).

Appendix E Details on example in Fig. 2

In the example system from Fig. 3, the coupling matrices for the two link types are $B^{(1)}=(\begin{smallmatrix}b&b\\ 0&0\end{smallmatrix})$ and $B^{(2)}=(\begin{smallmatrix}0&0\\ 0&b\end{smallmatrix})$, where the constant $b$ represents the coupling strength common to all external sublinks. The coupling matrix $F^{(i)}$ for internal sublinks is chosen from the binary set $\mathcal{F}=\{(\begin{smallmatrix}0&a\\ 0&0\end{smallmatrix}),(\begin{smallmatrix}0&0\\ a&0\end{smallmatrix})\}$, corresponding to the two possible sublink directions [and thus to two types of nodes indicated by green and cyan color, respectively, in Fig. 3(a)], where the constant $a$ represents the coupling strength common to all internal sublinks. The Lorenz dynamics of the subnodes and the coupling represented by sublinks are given by

[TABLE]

with the standard parameters, $\gamma=10$, $\rho=28$, and $\beta=8/3$.

The stability function $\psi(\lambda)$ is determined by Eq. (12), which for this system reads

[TABLE]

where $\bm{\xi}:=(\xi_{1},\xi_{2},\xi_{3})^{T}$ is the variation of the state vector $\bm{x}$ and the synchronous state $\bm{s}:=(s_{1},s_{2},s_{3})^{T}$ satisfies the equation for a single isolated Lorenz oscillator:

[TABLE]

For a given $\lambda$ in the complex plane, we compute $\psi(\lambda)$ by numerically integrating Eqs. (15) and (16) for $2\times 10^{4}$ time units and estimating the MLE Wolf et al. (1985) associated with the variable $\bm{\xi}$. Figure 6 shows the resulting estimate of $\psi(\lambda)$, which has a bounded stability region $\{\lambda\in\mathbb{C}\,|\,\psi(\lambda)<0\}$.

For a given combination of $a$ and $b$, we obtain $\Psi_{=}$ and $\Psi_{\neq}$, which are shown in Fig. 3(b). Note that for this example there are only two distinct homogeneous systems and two distinct heterogeneous systems. One of these heterogeneous systems is shown in Fig. 3(a). We also note that $\Psi_{=}>0$ and $\Psi_{\neq}<0$ are equivalent to the conditions (C1)′ and (C2)′ in Appendix D, respectively. For each combination of $a$ and $b$ satisfying both conditions [on a grid covering Fig. 3(b) with a resolution of $0.2$], we additionally run $24$ direct simulations of Eq. (2) for $200$ time units. The initial condition $\bm{x}^{(i)}_{\ell}(0)$ for each subnode is chosen randomly and independently from the uniform distribution in the region $[0,10]\times[0,10]\times[0,10]$ of its state space. The results confirm that the synchronization error $e$ defined in Eq. (13) and averaged over the last $100$ time units does not fall below $10^{-3}$ in all $24$ runs for both homogeneous systems, providing solid evidence that the system satisfies the AISync condition (C1) for these combinations of $a$ and $b$. Since $\Psi_{\neq}<0$ implies (C2)′ and thus (C2), this confirms AISync in the region shaded purple in Fig. 3(b).

The initial condition for the sample trajectory in Fig. 3(c) is chosen randomly within a distance of $10^{-3}$ from the synchronous state. The trajectory is then computed by integrating the system with all nodes green for $t\leq 25$, instantaneously switching the direction of the sublink between subnodes $2^{\prime}$ and $2^{\prime\prime}$, and then continuing to integrate for $25\leq t\leq 50$.

Appendix F Sampling protocol used in Fig. 4

We randomly sample systems whose network structure $A^{(\alpha)}$ is a circulant graph (with directed links of possibly multiple types) of given size $N$ and external sublink in-degree $D$ (i.e., the total number of sublinks received by the subnodes of a given node). Each of the $D$ sublinks coming into node $1$ is chosen randomly; it connects a random subnode chosen uniformly from the other $N-1$ nodes to a random subnode chosen uniformly from node $1$. The incoming sublinks into nodes $2$ to $N$ are then chosen to precisely match those coming into node $1$, which ensures that the network structure is a circulant graph. This simultaneously specifies $A^{(\alpha)}$ and $B^{(\alpha)}$ defining the system. To determine $\sigma_{\neq}$, $\sigma_{=}$, and $r$ for this system, we calculate the eigenspread $\sigma$ of the monolayer network representation for all the possible internal sublink configurations $F^{(i)}$, chosen here from the binary set $\{(\begin{smallmatrix}0&1\\ 0&0\end{smallmatrix}),(\begin{smallmatrix}0&0\\ 1&0\end{smallmatrix})\}$. For each combination of $N$ and $D$, we generate a sample of $4{,}000$ such systems to compute the fraction of AISync-favoring networks.

Appendix G Approximate symmetry in Fig. 4(a)

The approximate symmetry with respect to the vertical line at density $0.5$ observed in Fig. 4(a) can be explained using the notion of network complement. The complement of a given (unweighted) network with adjacency matrix $\widetilde{A}=(\widetilde{A}_{jj^{\prime}})$ is defined as the network having the adjacency matrix $\widetilde{A}^{c}=(\widetilde{A}_{jj^{\prime}}^{c})$ given by

[TABLE]

The external sublink density of a network and its complement add up to one, placing them symmetrically about the vertical line at density $0.5$ in Fig. 4(a). When the nontrivial Laplacian eigenvalues of the network and its complement, which we denote $\lambda_{2},\ldots,\lambda_{n}$ and $\lambda_{2}^{c},\ldots,\lambda_{n}^{c}$, respectively, are both indexed in the order of increasing real part, they are related by $\lambda_{j}+\lambda_{n+2-j}^{c}=n$ Nishikawa and Motter (2010). This implies that, if $\sigma$ is the eigenvalue spread for a monolayer network with given internal sublink configurations $F^{(i)}$, then the spread for its complement is given by

[TABLE]

where $\widetilde{m}:=\sum_{j}\sum_{j^{\prime}\neq j}\widetilde{A}_{jj^{\prime}}$ is the number of directed links in the network $\widetilde{A}$. Now consider two systems with $n$ subnodes and adjacency matrices $\widetilde{A}_{1}$ and $\widetilde{A}_{2}$, whose $\sigma$ values are $\sigma_{1}$ and $\sigma_{2}$, respectively. If we denote the $\sigma$ values of the complement of these systems by $\sigma_{1}^{c}$ and $\sigma_{2}^{c}$, respectively, we have

[TABLE]

when $\widetilde{A}_{1}$ and $\widetilde{A}_{2}$ have the same number of directed links, i.e., $\widetilde{m}_{1}=\widetilde{m}_{2}$. It follows from Eq. (18) that if $\widetilde{A}_{1}$ is the best homogeneous system and $\widetilde{A}_{2}$ the best heterogeneous system for a given external connection pattern with density $x$, then their complements are the best homogeneous and heterogeneous system for an external connection pattern with density $1-x$. Thus, the value of AISync strength $r$ is the same at density $x$ and $1-x$. The symmetry, however, is not perfect between sparse and dense parts of the plot, since we exclude the cases in which the network is not synchronizable (i.e., we require $\min_{j\geq 2}\text{Re}(\lambda_{j})>0$), the effect of which is not symmetric between sparse and dense cases.

Appendix S1 Example of AISync with unbounded stability region

Figure S1 shows another example of AISync using systems in which subnode dynamics is the same Lorenz oscillator as in Fig. 2, but with a different coupling function [leading to a different stability function $\psi(\lambda)$] and a different symmetric network. We use a three-node symmetric network with a single link type [with the corresponding coupling matrix $B^{(1)}=(\begin{smallmatrix}0&b\\ 0&0\end{smallmatrix})$], as can be seen in Fig. S1(a), where the monolayer representation of an example heterogenous system is shown. The coupling matrix $F^{(i)}$ for internal sublinks are chosen from the quaternary set $\{(\begin{smallmatrix}0&0\\ 0&0\end{smallmatrix}),(\begin{smallmatrix}0&a\\ 0&0\end{smallmatrix}),(\begin{smallmatrix}0&0\\ a&0\end{smallmatrix}),(\begin{smallmatrix}0&a\\ a&0\end{smallmatrix})\}$, so there are four distinct homogeneous systems. Two of them have $\lambda_{2}=0$ (and thus are not synchronizable), and one of the remaining ones is always more stable (i.e., smaller $\Psi$) than the other in the range of $a$ and $b$ considered in Fig. S1.

The more stable homogeneous system is the one with the dashed sublink in Fig. S1(a). Among the heterogeneous systems, the one with the smallest value of $\sigma$ is the one without the dashed sublink in Fig. S1(a), which is optimal (i.e., $\sigma=0$). We thus have just one homogeneous system and one heterogeneous system to compare. For the coupling function, we use $\bm{h}(\bm{x})=(x_{1},0,0)^{T}$, which leads to the stability function $\psi(\lambda)$ shown in Fig. S1(b). Unlike the example in Fig. 2, the stable region for this example (the region to the right of the red curve) is unbounded. Using this $\psi(\lambda)$, we calculate the MTLEs for the homogeneous and heterogeneous systems, which are shown in Fig. S1(c) as functions of parameters $a$ and $b$. We confirm that the global stability also follows the same trend by integrating Eq. (2) directly for $200$ time units from a random initial condition and then computing the synchronization error $e$ defined in Eq. (D1). The result is shown in Fig. S1(d), where we see that $e$ behaves similarly as the MTLE. Figure S1(e) shows the synchronization error $e$ for a typical trajectory of the system with $a=b=9$ [marked by black cross in Fig. S1(c)], which switches from homogeneous to heterogeneous at $t=150$ with the removal the sublink inside the pink node in Fig. S1(a).

Appendix S2 Example of AISync in undirected networks

The simplest example of undirected network exhibiting AISync is the network of two nodes, coupled bidirectionally. Consider a system with this network structure in which node $1$ has subnodes $1^{\prime}$ and $1^{\prime\prime}$, and node $2$ has subnodes $2^{\prime}$ and $2^{\prime\prime}$, and the network of sublinks (all with weight equal to one) forms a directed chain connecting the four subnodes as $2^{\prime\prime}\to 1^{\prime}\to 1^{\prime\prime}\to 2^{\prime}$. This network, when described at the node level, is indeed undirected because the pattern of sublink connections from node $1$ to $2$ is identical to the pattern of connections from node $2$ to $1$. The nodes are heterogeneous because subnode $1^{\prime}$ is connected to $1^{\prime\prime}$, while subnode $2^{\prime}$ is not connected to $2^{\prime\prime}$. Since the directed chain is an optimal network with $\lambda_{2}=\lambda_{3}=\lambda_{4}=1$, this heterogeneous system is more synchronizable than any combination of internal sublink configurations that leads to a homogeneous system. Thus, this undirected network exhibits AISync.

Appendix S3 Supplementary Tables

For each of the symmetric networks that can be derived from the network diagrams in Table I for $N=3$ and $4$, we identify all possible two-layer optimal heterogeneous systems with that symmetric network structure. Tables S1 and S2 show all these systems in the monolayer representation for $N=3$ and $N=4$, respectively