Critical Switching in Globally Attractive Chimeras

Yuanzhao Zhang, Zachary G. Nicolaou, Joseph D. Hart, 3 Rajarshi Roy, 3, 4 and Adilson E. Motter 5 Department of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208, USA Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, Maryland 20742, USA Department of Physics, University of Maryland, College Park, Maryland 20742, USA Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742, USA Northwestern Institute on Complex Systems, Northwestern University, Evanston, Illinois 60208, USA

Early work on chimera states focused mainly on networks of phase oscillators in the limit of large system size [34], where dimension reduction is often possible by employing the Ott-Antonsen ansatz [35][36][37]. For finite-size systems, some chimera states have been shown to be long transients [38], while others have been shown to be stable [39,40] using the Watanabe-Strogatz ansatz [41,42]. Recent research has placed increased emphasis on chimeras in finite-size networks of chaotic oscillators [43][44][45][46][47], which are important given the prevalence of chaos in physical systems [48]. In that context, it has been shown that the stability of chimera states can be studied rigorously using cluster synchronization techniques [46,47].
Even for permanently stable chimeras, an important question is how carefully one has to prepare the initial conditions in order to observe them. Early examples of chimera states required specially prepared initial conditions [11,14,49], while more recent examples include chimera states that emerge from a wide range of initial conditions [17,27,[50][51][52][53]. In the presence of global feed-back control, some chimeras have even been observed to attract almost all initial conditions [54,55]. However, whether globally attractive chimeras can emerge in the absence of control is still an open problem.
Because of the symmetry-broken nature of chimera states, another important question concerns the coexistence of multiple chimeras [49] and the possibility of transitions between them [56]. When multiple chimeras coexist, adding fluctuation or mismatch terms may induce switching events between them. This phenomenon has been studied under the name of "alternating chimeras" [20,57,58]. In previous studies, finite transition barriers must be overcome for transitions between otherwise persistent chimeras to occur. Accordingly, the transition rates are expected to scale exponentially with noise intensity.
Here we report on switching chimeras, which are chimera states that both exhibit power-law dependence of the switching frequency on noise intensity and attract almost all initial conditions in the absence of control. A switching chimera is comprised of two symmetric metastable states-referred to as subchimeras-between which the switching occurs. The power-law switching dynamics is a signature of critical behavior and stems from a vanishing quasi-potential barrier between the two metastable states. It follows that the switching persists indefinitely for any nonzero noise intensity. Strikingly, when the noise intensity is strictly zero, the symmetric subchimeras are linearly stable. Thus, the deterministic dynamics settle into one of the two subchimeras and, as in the original studies of chimeras, the state symmetry is broken. For any nonzero noise intensity, however, the long-term dynamical symmetry is restored due to the persistent switching between the two subchimeras. This dependence on noise intensity shares similarities with singular limits [59], in that the asymptotic network dynamics are qualitatively different for zero and small noise. Our analytical and numerical results are further validated by an experimental demonstration of switching chimeras in networks of optoelectronic oscillators. We suggest that switching chimeras can find applications in the study of intermittently alternating dynamics in biological systems and the development of approaches to measure small experimental noise.
The paper is organized as follows. In Sec. II, we introduce a representative system exhibiting switching chimeras. The power-law dependence between the average switching period and noise intensity is presented in Sec. III A. This critical switching behavior is then established and explained from various angles in the subsequent subsections. In Sec. III B, we show that it arises robustly in a first-exit model derived from an extension of the Freidlin-Wentzell theory. In Sec. III C, we further elucidate the mechanism underlying the switching dynamics by describing the dominant transition paths and the role of invariant saddles. In Sec. III D, we relate the scaling in the switching dynamics with the existence of transition paths of arbitrarily small action and compare it to critical phenomena in phase transitions. In Sec. III E, we establish a connection between power-law switching and intermingled basins of attraction. Experiments confirming switching chimeras and their power-law scaling in a network of optoelectronic oscillators are presented in Sec. IV. In Section V, we discuss connections between switching chimeras and other phenomena in physical and biological systems. Finally, we present our concluding remarks in Sec. VI.

II. COMPUTATIONAL OBSERVATION OF SWITCHING CHIMERAS
We consider 2n-node networks formed by two rings of n nodes, with nearest neighbor coupling of strength σ in each ring. The two rings are all-to-all coupled by weaker links of strength cσ for some 0 < c < 1. In this way, all the nodes are identically coupled, as shown by the network diagram in Fig. 1(a). We assume the oscillators are diffusively coupled, so the network can be represented through a Laplacian matrix in the dynamical equation. Adding to each node an uncorrelated Gaussian noise term of zero mean and tunable standard deviation ξ (which we refer to as the noise intensity), and writing down the coupling explicitly, the resulting stochastic dynamical equation for the first ring reads: intrinsic dynamics intracluster coupling Gaussian noise where N (1) i is Gaussian noise with unit standard deviation, and the superscripts indicate which ring the variables are associated with. The dynamical equation for the second ring can be expressed similarly. (We note that it is not essential for the dynamics to be discrete; an example of switching chimeras in systems with continuoustime dynamics is presented in Supplemental Material [60], Sec. S5. ) We first assume that the dynamics of each node is governed by a logistic map f (x) = x(1 − x). For concreteness, we also set n = 6 and c = 0.2 unless mentioned otherwise. Using a generalization of the master stability function formalism developed in Ref. [8], we can calculate the maximum transverse Lyapunov exponent associated with chimera states efficiently (Appendix A). In particular, we find parameters under which i) the two clusters cannot be simultaneously in stable synchronous states (i.e., any solution satisfying x (1) ii) one of the clusters can be in a stable synchronous state if the other cluster is not.
Inside the region where both conditions are satisfied, coherence is induced by incoherence, meaning that synchronization in one cluster is stabilized by desynchronization in the other cluster. Figure 1(b) shows that the system in Fig. 1(a) has a large parameter region (purple) in which these two conditions are satisfied. In that region, chimera states are linearly stable and do not coexist with stable globally synchronized states.
To confirm that the desynchronized ring is indeed in an incoherent state, we run direct simulations [61] from random initial conditions for 10 8 iterations under noise of intensity ξ = 10 −10 . Figure 1(c) shows representative trajectories and associated synchronization errors for σ = 1.7 and r = 3.05. The synchronization error in the j-th cluster is defined as The system exhibits not only chimera dynamics but also persistent transitions in which the coherent and incoherent rings switch roles: as one ring loses synchrony and becomes incoherent, the other ring synchronizes. Moreover, as we show below, the switching observed here is critical-the transition rate depends on noise intensity as a power law and switching can be triggered by arbitrarily small noise. This power-law dependence distinguishes switching chimeras from previously reported "alternating chimeras," in which either the transitions are forced by large fluctuation terms [20,[56][57][58] or rely on heteroclinic dynamics [62][63][64]. In the first case, there are finite barriers separating the different states, while in the second case each state is inherently unstable and switching occurs in the absence of noise.
The persistence of switching chimeras under many transition cycles suggests it is globally attractive. To verify that this is indeed the case, we evolve the system for 10 4 iterations starting from 10 6 different random initial conditions for σ = 1.7 and r = 2.9, 2.95, 3.0, 3.05 [dots in Fig. 1(b)]. In all tests, the oscillators are swiftly attracted to the chimera state and no other attractors are observed. Further evidence of this global attractiveness is presented in Supplemental Material [60], Sec. S1, where we also demonstrate the prevalence of switching chimeras across different cluster sizes, inter-cluster coupling strengths, and intra-cluster coupling range.

III. POWER-LAW SWITCHING
A. Extreme sensitivity to noise Next, we present numerical results characterizing the effect of noise intensity ξ on the average switching period T . Figure 2 shows that, as one approaches the boundary of the chimera region [from the green dot to the orange dot in the bottom right of Fig. 1(b)], T decreases and switching becomes more frequent. For each fixed value of r, the average switching period increases as the noise intensity decreases, with scaling that follows a power law. It is remarkable that even noise of intensity as small as ξ = 10 −15 (the resolution limit of computers using double-precision floating-point format) can induce frequent switching.
This switching between the coherent and incoherent clusters does not contradict the fact that synchronization in one cluster is linearly stable if the other cluster is incoherent. This is because linear stability analysis assumes the perturbations to be infinitesimally small. Finite-size perturbations, no matter how small, can still grow large enough along the unstable portions of a chaotic attractor to disrupt synchrony in the coherent ring and induce switching.
The power-law scaling of the average switching period, and consequently the extreme noise sensitivity of chimera states, makes the switching behavior observed here "anomalous" in the sense that it appears to contradict the Freidlin-Wentzell theory [65]. According to that theory, for a stochastic system with deterministic dynamics F and noise term of intensity ξ, the rate of transition from one metastable [66] state A to another state B scales as exp(−S A→B /ξ 2 ) and the first exit time scales as exp(S A→B /ξ 2 ) [67]. Here, S A→B is the infimum of the Freidlin-Wentzell action among all paths X connecting state A to state B: The infimum of the action measures how much one has to work against the deterministic part of the dynamics to induce a transition from A to B. This quantity is also known in the literature as a quasi-potential barrier [68], and is analogous to a potential barrier for transitions in gradient systems.

B. First-exit problem in log-error space
Although the power-law scaling observed for switching chimeras and the exponential scaling predicted by the Freidlin-Wentzell theory seem incompatible at first glance, we can establish a connection between them. We first note that the synchronization error inside the coherent ring usually fluctuates close to an error floor determined by the noise intensity, but a switching can be triggered by rare events that drive the error all the way to an error ceiling determined by the synchronization error of the incoherent ring [for an example, see the middle panel of Fig. 1(c)]. Moreover, since the variational equation acts multiplicatively on the synchronization error (see Appendix A), the error naturally evolves on a log scale as long as the linearization around the synchronization manifold is still valid.
Motivated by these observations, we focus on an attribute , defined as the logarithm of the synchronization error inside the coherent ring: As a first approximation, the dynamics of can be modeled as a biased one-dimensional random walk confined within two boundaries, corresponding to the error floor and the error ceiling. At each step, has probability p of moving up a fixed distance d 1 and probability 1 − p of moving down a distance d 2 . The random walker starts from the error floor, and it never goes below that boundary. Every time reaches the error ceiling, we consider that a switching event has occurred and reset to the lower boundary. An illustration of this process can be found in Fig. 3(a).
To derive a relation between the average switching period T and the inter-boundary distance D in the random walk model, we note that when pd 1 < (1 − p)d 2 and D d 1,2 this is a first-exit problem. Thus, according to the Freidlin-Wentzell theory, where λ is some constant determined by p, d 1 and d 2 . Now recall that D is determined by the distance between the error floor and error ceiling. The error floor is given by ln(ξ) and, without loss of generality, we set the error ceiling to be 1. Thus, D = ln(1) − ln(ξ) = ln(ξ −1 ), and Eq. (6) leads to This reproduces the power-law relation between the average switching period T and the noise intensity ξ observed in Fig. 2.
We now turn to a more quantitative analysis to support the idea that the switching events in the original system can be inferred from the one-dimensional attribute . Starting with the system in Eq. (1), we compute the growth rate of the synchronization error in the coherent at each iteration. The distribution of this quantity, which we call the local Lyapunov exponent, is shown in Fig. 3(b) for σ = 1.7 and r = 3.05. Of all the local Lyapunov exponents sampled, 35% were negative, with a mean of −0.46; the remaining 65% of the exponents were positive, with a mean of 0.19. Because e is a one-dimensional variable, the Lyapunov exponent that determines its asymptotic stability is given by averaging over the local Lyapunov exponents from t = 0 to t = ∞. Since −0.46 × 0.35 + 0.19 × 0.65 < 0, although 65% of the chaotic attractor is repelling, the chimera state is actually linearly stable. From the above information, we can set p = 0.65, d 1 = 0.19, and d 2 = 0.46 in our ran- FIG. 4. Dominant transition pathway between the two symmetric subchimeras, which consists of the intermediary stages T1 to T4. Only T2 requires activation from noise, which can be arbitrarily small but not strictly zero; all other transitions follow directly from the deterministic dynamics of Eq. (1). In particular, T3 and T4 follow the stable and unstable manifolds of the invariant saddle, respectively.
dom walk model, and calculate the relation between the average switching period T and the noise intensity ξ. The brown circles in Fig. 3(c) indicate how T scales with ξ for this random walk model; they follow a welldefined power law, as expected from Eq. (7). But it is also clear that a random walk is not a very accurate picture for the dynamics of , since the predicted average switching periods are much larger than the ones obtained from simulating Eq. (1) (orange circles). This discrepancy is partially due to the crude approximation we made when fixing the step sizes of the random walk to be constants. If we choose the step size as well as the direction of the random walk according to the distribution in Fig. 3(b), we observe the scaling indicated by cyan circles in Fig. 3(c), which is closer to the true scaling. However, the predicted exponent of −0. 26 is not yet close to the true value of −0.12, which indicates that something is still missing.
The approach we just took is equivalent to shuffling the time series of the local Lyapunov exponents and using the shuffled sequence to generate the random walk. This preserves the full distribution information, but ignores temporal correlations. Because the stable and unstable portions of a chaotic attractor are usually not well mixed, the actual evolution of is a non-Markovian process, and we expect the temporal information to be relevant. This effect will tend to correlate the upward movements of , which in turn makes it more likely for to reach the error ceiling and shortens the average switching period for small noise. When the temporal information is incorporated into the model (by using the original sequence of local Lyapunov exponents rather than randomly sampling them), we arrive at a more realistic model for the switching dynamics, which takes the form of a deterministic walker. The prediction of this refined model (purple circles) is in excellent agreement with the true scaling (orange circles).
It is important to note that the power-law scaling is preserved even after we allow variable step sizes and strong correlation between steps in our model. This sug-gests that Eq. (7) is robust and that power-law switching is expected for a general class of systems. Transitions in such systems can be modeled as a first-exit problem in which the distance to the exit increases linearly with the logarithm of the inverse of noise intensity.

C. Transition pathways
We can gain a deeper understanding of the switching dynamics by investigating the transition paths connecting the two symmetric subchimeras. One natural question concerns whether there is a single pathway or multiple pathways for the switching? If multiple pathways exist, do they intersect at key intermediate states? For the system in Fig. 1(a), with n = 6, it turns out that there is only one dominant pathway when noise is small. We illustrate the key transitions (T1 to T4) and intermediate states of this pathway in Fig. 4. We later analyze an explicit realization of this pathway in Fig. 5, which provides strong numerical support for the following transition sequence: (T1) Starting from one of the subchimeras, the incoherent ring occasionally visits near-synchronized states (temporary clustering in Fig. 4). (T2) The temporary clustering in the incoherent ring strongly correlates with the instability windows in the coherent ring. This is not surprising since states with both rings synchronized are unstable. Within those short windows, small noise or perturbations applied to the coherent ring will be amplified and lead to a short-wavelength bifurcation. That is, the coherent ring partially desynchronizes and splits into two alternating groups with different dynamics (oscillators in the same group remain synchronized). Reaching this "activated state" is the only stage in which noise is needed, even though it can be arbitrarily small. (T3, T4) The state between T3 and T4 lives in an invariant subspace induced by the rotational symmetry in each ring. In fact, the state is an invariant saddle and serves as the key intermediate state connecting the two subchimeras. During T3, the system moves along the stable manifold of the invariant saddle and the 6 oscillators in the upper ring converge to a synchronized state. During T4, the system moves away from the saddle following its unstable manifold, where the partially desynchronized state in the lower ring evolves into an incoherent state. This concludes the entire sequence of transitions from one subchimera to the other.
The short-wavelength perturbation where the i-th component of this vector is to be interpreted as a perturbation to the i-th node in the ring, is the dominant instability in the coherent ring according to our linear stability analysis, and is the one being excited by noise during the transition T2. To further support this claim, we ran direct simulations of Eq. (1), but with ∆ sw filtered out from the noise applied to each ring. This time, for noise intensity ξ ≤ 10 −9 , the average switching period T becomes independent of ξ and always equals the average switching period induced by round-off errors, as shown in Appendix B. This confirms that the overwhelming majority of the switching events must be initiated through a short-wavelength bifurcation in the coherent ring when noise is small [69].
To better visualize the subchimeras and the invariant saddles, we project them onto the mean state of each ring: i /n. Figure 5(a) shows the projection of the two symmetric subchimeras colored in blue and orange, respectively. We can see the fine structure of the subchimeras under this projection, which is indicative of their fractal nature. In Fig. 5(b), we show the projection of the two invariant saddles (red and green).
We now try to explicitly find a least-action path connecting the two subchimeras, which can be challenging even for transitions between fixed points or periodic orbits [67,68]. In our case, the high dimensionality and the chaotic nature of the subchimeras make the optimization of the transition path extremely difficult when using traditional methods. Fortunately, the mechanism presented in Fig. 4 points to an efficient way of finding paths of arbitrarily low action connecting the two subchimeras. We simply wait for the incoherent ring to visit a nearsynchronized state, then introduce a one-time perturbation in the form of ∆ sw (δ) to excite the short-wavelength bifurcation in the coherent ring. If a transition is successfully triggered, the action of the transition path is simply Using this strategy, we can easily find a transition path with action as small as 10 −28 (i.e., δ around 10 −14 ), which is shown in Fig. 5(c) and (d) for different projections.
The coordinate e 1 (e 2 ) in Fig. 5(c) is defined as the sum of the synchronization error among the odd oscillators and the synchronization error among the even oscillators in the first (second) ring. For this projection, the two subchimeras are found in the upper left and lower right corners, while the key invariant saddle connecting the two subchimeras is projected onto the lower left corner (e 1 = e 2 = 0). It is informative to view the projected transition path in Fig. 5(c) in light of the pathway shown in Fig. 4: the first two transitions (T1 and T2) correspond to the upper left corner while the other two transitions (T3 and T4) loop around the lower left corner as they follow the stable and unstable manifolds of the invariant saddle closely. Conversely, the projected path provides strong numerical support for the pathway illustrated in Fig. 4. However, the evidence is not yet conclusive, as states with both rings synchronized also project onto the lower left corner for the coordinates in Fig. 5(c). Could the two subchimeras be connected by an unstable synchronized state instead of the invariant saddles in Fig. 4? The projection to the synchronization errors e 1 and e 2 in Fig. 5(d) excludes this possibility, since the path goes through the upper right corner (both rings desynchronized) rather than the lower left corner (each ring synchronized). Multiple transition paths with action ranging from 10 −30 to 10 −10 were tested, and they were all qualitatively identical to each other under both pro-jections. This evidence further supports the existence of a dominant transition pathway for the observed switching between subchimeras.

D. Connections with critical phenomena
The fact that switching can be induced by arbitrarily small noise but not in the absence of noise implies that i) no matter how small the action of a transition path, we can always find another path with even smaller action, and ii) there is no zero-action path of finite length connecting the two subchimeras. Thus, a least-action path does not exist in our system. Instead, given an arbitrarily small upper bound on the available action, there are always finite-length transition paths that meet that constraint. It follows that the infimum of the action over all transition paths (i.e., the quasi-potential barrier S separating the two subchimeras) vanishes. In Fig. 6(a) we show that this is indeed the case by applying a single perturbation ∆ sw (δ) to the coherent ring, with δ ranging from 10 −5 to 10 −15 . The distribution of times a transition path is found through this procedure shows that the landscape is highly nontrivial for paths of small action: transition barriers of all heights exist and the height distribution follows a power law. This claim is further supported by Fig. 6(b), where we show the action for 1000 different transition paths, each obtained by applying ∆ sw (δ) at a different time t (the same initial condition is used for all simulations). One can see that the landscape varies wildly and the associated action spans many decades. As we include more transition paths, deeper and deeper valleys can be found, bringing the smallest action ever closer to zero.
The power-law distribution of barrier heights in turn gives rise to the power-law scaling of the average switching periods shown in Fig. 2. This is because the only transition paths that matter are the ones with action comparable to the noise intensity. Although there are many more higher-action paths, the probability of crossing those barriers is exponentially smaller. The argument is further supported by the scaling exponents in Figs. 2 and 6, which differ only by a negative sign.
There are intriguing parallels between what we find here and critical phenomena in second-order phase transitions [71,72]. For instance, in site percolation models, the correlation (which quantifies the likelihood of two sites being connected) decays exponentially with distance when the occupation probability is p < p c , but the decay changes to a power law at the critical point p = p c . Here, the average switching period scales exponentially with the inverse square of noise intensity, ξ −2 , when the quasi-potential barrier has S > 0, but it is replaced by a power law when S = 0. There are finite barriers of all heights between the two subchimeras when S = 0; similarly, in percolation, there are finite clusters of all sizes at the critical point p = p c . The power laws uncovered here, however, are more robust than those from the percolation  Fig. 2. Paths with arbitrarily small action exist but small-action paths become increasingly more difficult to find as the available action is decreased, resulting in power-law relationships between the probability p and the perturbation size δ. Notice that the scaling exponents here match those in Fig. 2. (b) Minimum action ( 1 2 δ 2 ) needed to induce a transition by applying ∆sw(δ) at a given time t, for ξ = 0, σ = 1.7, and r = 2.95. This highly structured profile can be regarded as a visualization of the transition-barrier landscape for switching chimeras.
theory. The latter only happens at the critical point and requires fine-tuning, whereas here the power-law switching persists for a wide range of parameters. In this sense, the analogy is perhaps closer with self-organized criticality [73][74][75], in which scale-invariance emerges in the absence of fine-tuning.

E. Intermingled basins
By now we have explained the "anomalous" power-law switching behavior from a first-exit model in log-error space (Sec. III B) as well as by characterizing the action landscape of transition paths (Sec. III D). In those characterizations one can catch glimpses of chaos lurking in the background, but its exact role is still unclear. In this section, we establish a direct connection between power-law switching and riddled basins [76][77][78][79][80][81][82][83], which is only possible for chaotic attractors [84], thus bringing the fundamental importance of the chaotic dynamics to the forefront.
Chaos has long been known to produce power laws There is a symmetry between the two basins with respect to reflections across the diagonal, which originates from the reflection symmetry of the network. The areas marked for magnification are intentionally over-sized to facilitate visualization. The choice of state-space section and system parameters are specified in the text.
by generating fractal structures in state space [70]. For example, in the presence of fractal basin boundaries, a small uncertainty ε in the initial conditions translates to an uncertainty of Aε α percent on the final states, where prefactor A is a constant and α is the uncertainty exponent given by the difference between the state-space dimension and the box-counting dimension of the basin boundary [85]. In the case of riddled basins, the entire basin is its own (fractal) boundary and α = 0. This means that, for any ε, the ε-neighborhood of an arbitrary point in a riddled basin will always include points that are in the basin of some other attractor [70]. In Fig. 7, we show a two-dimensional section of the twelve-dimensional state space to visually illustrate that the attraction basin of each subchimera is riddled. Because the two basins are mutually riddled, they are referred to as intermingled basins. In this figure, the initial conditions for x (1) 6 and x 6 /2, where 1 ≤ i ≤ 5. We then simulated Eq. (1) for σ = 1.7 and r = 2.95 in the absence of noise, and recorded the subchimera attractor each trajectory is attracted to. (There is nothing special about the choice of the parameters or the section of the state space, since other choices lead to similar results.) One can observe intricate fractal-like structures in all parts of the twodimensional section, for all resolutions considered (up to pixels of size 10 −10 × 10 −10 ). There is also a symmetry between the two basins. If an initial condition is in the basin of one subchimera, then its mirror image reflected along the diagonal line must be in the basin of the other subchimera (i.e., if (x a) is orange). This is the result of the reflection symmetry between the two rings in Fig. 1(a).
Because the basins are intermingled, the basin of one subchimera has points arbitrarily close to the other subchimera attractor, and vice versa, which gives rise to arbitrarily small transition barriers in Fig. 6. This means that the subchimeras are attractors in the sense of Milnor [86] (i.e., attracts initial conditions of nonzero measure), but not in the sense of attracting an open neighborhood of initial conditions containing the attractor.
Apart from the Freidlin-Wentzell action, the perturbation magnitude δ in Fig. 6 can also be interpreted as a distance from the closest subchimera attractor. The probability p then measures the fraction of state space that converges to the opposite subchimera, when at distance δ from the subchimera attractor. As the initial conditions are taken further away from one subchimera, it becomes more likely for the system to land in the basin of the other subchimera. Conversely, as δ → 0, the probability of escaping to the opposite subchimera approaches zero algebraically. This is visualized using a transverse section of the intermingled basins that directly connects the two subchimera attractors, as shown in Appendix C.
Although arbitrarily small perturbations can drive the system out of a subchimera attractor, both subchimeras are transversally stable according to linear stability analysis. While seemingly incompatible, these two conditions can coexist when an attractor is transversally stable for the natural measure but unstable for some other invariant ergodic measure. In fact, transversal stability for the natural measure and instability for at least one other invariant ergodic measure are necessary conditions for riddled basins to occur [84]. This mathematical statement is, in its core, similar to the intuitive explanation given in Sec. III A on why a system can be driven away from a linearly stable state by arbitrarily small perturba tions.

IV. EXPERIMENTAL OBSERVATION OF SWITCHING CHIMERAS
Thus far, we have focused on the theoretical analysis of networks of logistic maps, which reveals remarkable features of a new chimera state, including intermingled basins and switching triggered by arbitrarily small noise. To demonstrate that the theoretical results can be observed under realistic conditions and for different oscillator dynamics, we performed experiments on networks of coupled optoelectronic oscillators. As we show next, our experiments confirm the existence of switching chimeras in physical systems.
The experimental setup is schematically shown in Fig. 8(a). A single optoelectronic oscillator draws non-linearity from a Mach-Zehnder modulator, which takes voltage x as an input and outputs light of intensity sin 2 (x + φ). The operation point φ is fixed at π/4 throughout the experiments. Time-multiplexing and delays are used to realize multiple oscillators from a single time-delayed feedback loop, which reduces apparatus costs and allows for the realization of a large number of truly identical oscillators. The oscillators are coupled together by a digital filter implemented electronically on a field-programmable gate array (FPGA) according to a predetermined Laplacian matrix L = {L ij }. In this case, L describes the two-cluster network shown in Fig. 1(a). Further details of the optoelectronic system can be found in Refs. [87,88].
The main source of intrinsic noise comes from the measurement of light intensity, including the noise introduced by the analog-to-digital converter (ADC) due to its finite resolution. To best model the experimental system, we introduce independent Gaussian noise to the oscillators at each iteration: I x . The dynamical equation describing the optoelectronic oscillator network can then be written as where the noise term is implicitly included in I. In our experiments, we again set c = 0.2 and n = 6.
We first sweep the parameter space of feedback strength β and coupling strength σ using direct simulations of Eq. (9). As shown in Fig. 8(b), switching chimeras are predicted to occupy a significant portion of this space. Inside the switching chimera region (purple), the red and green dots denote the parameters to be systematically investigated in the experiments.  Fig. 8(b)]. (b) Portion of the experimentally measured time series used to generate (a). These measurements are performed at the base-noise level of the system, which is estimated to be 0.0019.
The dynamics exhibited by the experimental system is in many ways qualitatively similar to that of coupled logistic maps. In particular, a clear pattern of irregular switching between two subchimeras is observed for suitable parameters, as shown in Fig. 9(b) for a representative time series. To characterize the experimental dynamics quantitatively, we first test whether the powerlaw relationship between the average switching time T and noise intensity ξ holds in the experimental data. An important step in the data analysis is to estimate the level of the intrinsic experimental noise, which we do by simulating Eq. (9) under different ξ to extract T for a range of noise intensities. The simulation results are then compared with the T observed in the experiments. For both parameter sets (β = 1.3, σ = 1.05 and β = 1.3, σ = 1.1), the simulations with noise intensity 0.0019 agree best with the experiments. We thus choose Gaussian noise of intensity ξ 1 to approximate the base-noise level intrinsic to the experimental system. It is worth noting that this technique can in principle be extended to estimate the level of intrinsic noise in other oscillators, even when the noise is extremely small-an outstanding problem for which, to the best of our knowledge, no general approach currently exists.
To implement variable noise in the experiments, we introduce an additional Gaussian noise term of tunable intensity ξ 2 via the FPGA. Assuming that the intrinsic and external noise terms are independent, the experimental system is effectively subject to a Gaussian noise of intensity ξ = ξ 2 1 + ξ 2 2 . Figure 8(c) summarizes the experimentally measured T for different ξ from the lower bound 0.0019 all the way to 0.02. Each data point is averaged over at least 20000 experimentally observed switching events. It can be seen that the power-law relationship holds under realistic noise levels and is robust against the imperfections typical of an experimental system. In addition, we have also performed systematic simulations to further confirm that the power-law scaling persists in the presence of small oscillator heterogeneity (Appendix D). Figure 9(a) shows the distribution of the switching periods T extracted from 45000 switching events, for data collected from multiple experimental runs with β = 1.3, σ = 1.05, and ξ 2 = 0. The distribution of periods is clearly exponential. This is a consequence of the fact that, although the evolution of the synchronization errors e 1 and e 2 are non-Markovian (Sec. III B), the switching events themselves are described by a Poisson process. In particular, the experimental data shows that the waiting period until the next switching event is independent of the previous switching events. For such a memoryless process with a constant transition rate, the time between switching events is guaranteed to be exponentially distributed [89].
Our experimental results are further visualized using an animated spatiotemporal representation of time-series data presented in Fig. 9(b) (Supplemental Material [60], Sec. S2 and associated animation). As in the case of coupled logistic maps, the underlying state-space structure giving rise to this dynamics is the intermingled nature of the attraction basins. Indeed, direct simulations of Eq. (9) confirm that the basins of the two symmetric subchimeras are intermingled (Supplemental Material [60], Sec. S3).

V. CONNECTIONS WITH BIOLOGICAL AND OTHER PHYSICAL SYSTEMS
A switching chimera can be seen as a chimera state whose symmetry is not broken when considering the longterm dynamics-asymptotically, one cannot distinguish between the behavior of the two clusters. With this observation in mind, we can establish an intriguing parallel between the switching chimera and the symmetrybreaking phenomenon of dipole inversion [90]. Many small molecules, such as ammonia, have more than one (symmetry-broken) ground state with non-vanishing dipole moments. However, due to quantum tunneling, an ammonia molecule switches rapidly between its two ground states, canceling out the opposite dipole moments and restoring the broken symmetry. The same can be stated for switching chimeras, since each of the two symmetric subchimeras has a broken parity symmetry but the switching between them restores that symmetry. For larger and heavier molecules, such as sugars or phosphorus trifluoride, dipole inversion is no longer likely to be excited by quantum tunneling or even thermal fluctuations, and thus the symmetry is spontaneously broken and nonvanishing dipole moments persist. We observe that the tendency for transitioning between subchimeras also decreases in larger systems, with the average switching period growing exponentially as the number of nodes is increased (Supplemental Material [60], Sec. S4).
It is instructive to notice that an exponential dependence of the average switching period on system size is also observed for the magnetized states in the Ising model for any nonzero temperature below the critical point [91,92]. However, because there is a finite energy barrier to overcome for transitions between the magnetized states, the dependence of the average switching period on inverse temperature (the analogue of the inverse square of noise intensity in our systems) is not power law but instead exponential.
Switching between symmetry-broken states are not limited to physical systems. In particular, switching chimeras may have implications for aperiodic lateral switching in biological systems, of which interhemispheric switching in songbirds during vocal production is an example [93]. Other examples of lateral switching include alternating eye movement in chameleons and fish [94], switching in neural activity inside the two sinuses of leech hearts [95], and unihemispheric sleep in dolphins, birds, and other animals [96,97]. A common aspect of these various processes is that they involve alternations in the activity between two approximately symmetrical lateral sides. Despite previous progress [98], the underlying mechanism of lateral switching remains elusive. This is especially the case for aperiodic lateral switching, since such cases cannot be easily modeled by hypothesizing the existence of a central pattern generator or propagating wave dynamics, as in previous alternating chimeras [56,62,63]. In the case of the songbird zebra finches, for instance, the inter-hemispheric switching between song-control areas of the brain is highly irregular, characterized by switching intervals ranging from 4 ms to 150 ms [93]. Switching chimeras offer a simple mechanism by which a wide range of switching intervals can emerge naturally, and thus suggest the possibility that aperiodic lateral switching could be generated spontaneously (as opposed to, for example, being forced by neurotransmitter release [99]).

VI. CONCLUDING REMARKS
The theoretical, computational, and experimental results presented here offer a comprehensive characterization of a novel class of chimera states that are globally attractive and exhibit power-law switching dynamics. We extended the Freidlin-Wentzell theory to derive the observed power-law scaling and we demonstrated that there is no finite quasi-potential barrier separating the two symmetric subchimeras. This unexpected scaling behavior, which should be contrasted with the exponential scaling observed for typical noise-induced transitions [100,101], was confirmed under realistic conditions in our experiments using networks of optoelectronic oscillators. We also established a connection between switching chimeras and intermingled basins, which provides insight into both phenomena. In particular, the latter explains why switching between subchimeras occurs for arbitrarily small noise despite each subchimera being linearly stable.
We expect switching chimeras to be a common phenomenon in multilayer networks with symmetry. These networks are generalizations of the two-layer networks considered in Ref. [14]. In particular, switching between symmetric subchimeras is expected to be possible for networks formed by any number of identically-coupled identical layers, where the layers themselves can have an arbitrary structure. Thus, while we focused on networks with two subchimeras, our analysis extends naturally to other states and to a larger number of switching configurations. From the dynamical perspective, we point to the following conditions for the emergence of power-law switching behavior in general: (i) there are two or more attractors and they are embedded in manifolds of dimension lower than that of the state space; (ii) each attractor is chaotic and has transversally unstable periodic orbits embedded within. If the transitions are not restricted to chimera states, the requirement on the network structure can be further relaxed, as these conditions are often satisfied even by single-layer oscillator networks.
Finally, we note that the observed high noise sensitivity of the switching dynamics has far-reaching implications. It can be exploited, for instance, to detect small intrinsic noise in oscillator systems-e.g., by using timemultiplexing to create a network of such systems that exhibits power-law switching. It also offers a potential explanation for irregular switching noticed in biological systems, which suggests that the dynamical behavior described here may be observed in naturally evolved processes.

VII. ACKNOWLEDGMENTS
oscillators, where x i is the state of the i-th oscillator, f is the mapping function governing the uncoupled dynamics of each oscillator, L = {L ij } is the Laplacian matrix describing the structure of an undirected network with two nonintertwined identical clusters, h is the interaction function, and σ > 0 is the coupling strength.
Let L be the n × n Laplacian matrix that encodes the intra-cluster connection inside the coherent cluster, µ be the total strength of input connections each oscillator in the coherent cluster receives from the incoherent cluster, and x 1 = x 2 = · · · = x n = s be the synchronization manifold for the n oscillators in the coherent cluster. The variational equation describing the evolution of the deviation away from s can be written as where L = L + µ1 n , δX = (δx 1 , · · · , δx n ) = (x 1 − s, · · · , x n − s) , and ⊗ denotes the Kronecker product. Although the incoherent cluster does not enter the equation explicitly, it influences the matrix L and the synchronization trajectory s[t] through the inter-cluster coupling. We note that the input from the incoherent cluster faithfully accounts for the state of those oscillators and is time-dependent in general.
Equation (A2) can be decoupled into n independent equations by diagonalizing L: where η = (η 1 , · · · , η m ) is δX expressed in the new coordinates that diagonalize L, and v i = v i + µ are the eigenvalues of L in ascending order. Synchronization in the coherent cluster is stable if and only if Λ(σ v i ) < 0 for i = 2, . . . , n, where is the Lyapunov exponent of Eq. (A3) and v 2 , · · · , v n represent the modes transverse to the synchronization manifold of the coherent cluster. The maximum transverse Lyapunov exponent (MTLE) determining the synchronization stability is max 2≤i≤n Λ(σ v i ). A chimera state is stable for ξ = 0 if the MTLE for synchronization in the coherent cluster is negative under the influence of the incoherent cluster.

Appendix B: Dominant switching route
We now provide more evidence that short-wavelength bifurcation is the dominant mechanism to initiate switching between the two symmetric subchimeras. Again, we  Fig. 1(a) for σ = 1.7 and the noise is Gaussian (but with the short-wavelength component filtered out). The flatness of the fitting lines below ξ = 10 −9 confirms that short-wavelength bifurcation is the dominant route for chimera switching.
simulate Eq. (1) to extract the average switching period T for various levels of noise intensity ξ, but this time the short-wavelength component ∆ sw is filtered out from the noise applied to each ring. If a short-wavelength bifurcation is indeed the dominant route for switching, then one would expect the average switching period to become independent of the noise intensity after filtration. This is exactly the case shown in Fig. 10, where the slope becomes completely flat for each r when the noise intensity goes below 10 −9 (compare with Fig. 2). Due to the presence of round-off errors in our simulations, whose short-wavelength component cannot be filtered, switching can still be observed in the flat region at a rate induced by the round-off errors (noise intensity around 10 −16 ). When the noise intensity goes above 10 −9 , new switching pathways besides the short-wavelength bifurcation start to become available, as demonstrated by the resulting decrease of the average switching period.
Appendix C: Transversal section of intermingled basins Figure 11 shows the intermingled basins for a twodimensional section of the state space for the logistic map system in Eq. (1). This section is defined by where δ max is taken to be 0.2. For δ = 0, the first ring is synchronized and the second ring is incoherent (orange subchimera), while for δ = δ max , the second ring is synchronized and the first ring is incoherent (blue subchimera). Thus, this section of the state space directly connects the two symmetric subchimeras. As one approaches the orange (blue) subchimera, the points become predominantly orange (blue), but no matter how close δ is to zero (δ max ), speckles of blue (orange) dots can always be found. FIG. 11. Transversal section of the intermingled basins that directly connects the two symmetric subchimeras. This corresponds to a different state-space section of the same system considered in Fig. 7.
Appendix D: Robustness against oscillator heterogeneity In Fig. 12, we quantify the effect of oscillator heterogeneity on the switching dynamics, explicitly demonstrating the robustness of the switching chimeras. We start from a system of identical oscillators (the system in Fig. 1 for r = 3 and σ = 1.7) and introduce independent random perturbations to the parameter r of each oscillator, drawn from a Gaussian distribution of zero mean and standard deviation ∆. For ξ ≥ ∆, the average switching periods in the homogeneous and heterogeneous systems become indistinguishable, with both following a well-defined power-law distribution on noise intensity. For ξ < ∆, the effect of heterogeneity dominates the effect of noise; as a result, the average switching period (dashed lines) branch out of the original power-law relation (solid line) and approaches a constant determined by ∆. These results are largely independent of the particular realization of oscillator heterogeneity.