Quickly exploring weakly coupled oscillators Synchronizing oscillators have been useful models for exploring coupling in dynamic systems. However, many macroscopic platforms such as pendula evolve on slow time scales, which can limit the observation of states that emerge after many cycles. Matheny et al. fabricated a ring of eight nanoelectromechanical oscillators resonating at ∼2.2 megahertz with quality factors of ∼4000 that could be rapidly controlled and read out. Analysis of these large datasets revealed exotic synchronization states with complex dynamics and broken symmetries. Theoretical modeling showed that emergent higher-order interactions (such as biharmonic and next-nearest neighbor) stabilized complex dynamics, despite the network having weak nearest-neighbor coupling. Science, this issue p. eaav7932

Structured Abstract INTRODUCTION A paramount contemporary scientific challenge is to understand and control networks. General studies of networks are essential to a variety of disciplines, including materials science, neuroscience, electrical engineering, and microbiology. To date, most studies are observational or “top-down,” relying on phenomenological models of nodal behavior deduced from data extracted from observations of the entire network. On the other hand, oscillator synchronization provides a popular “bottom-up” experimental paradigm for studies of network behavior. Synchronization occurs when a large number of networked oscillators tend to phase-lock and reach global consensus, despite the presence of internal disorder (such as differences in oscillator frequencies). However, the state of global consensus is not the only dynamical state manifested within networks of coupled oscillators. Recent work has discovered long-lived states that spontaneously break the underlying symmetries of the network, even when its constituent nodes are identical. Understanding the mechanisms that underlie these exotic symmetry-breaking dynamics will be of general benefit to network science and engineering. To enable experimental studies with unprecedented control and resolution, we developed an oscillator network based on nonlinear nanoelectromechanical systems (NEMS). RATIONALE Complex oscillator dynamics emerge in a variety of settings. Previous experimental studies on symmetry breaking in oscillator networks have manifested such complex behavior only by implementing complicated coupling mechanisms designed for that purpose. By contrast, real-world networks are often dominated by simple coupling, so these previous experimental results cannot readily be generalized. Here, exotic states are seen to emerge within NEMS oscillator networks, just beyond the weak coupling limit in simple settings. RESULTS Using an array of coupled nonlinear NEMS oscillators, we observed spontaneous symmetry breaking in a simple and general network setting. These NEMS oscillator nodes were made to be nearly identical, and their fastest dynamical time scales were short enough to generate large datasets, permitting observation and statistical analyses of exotic, slowly emerging network phenomena. In addition, NEMS are very stable, so that transient effects within the network were able to relax within experimental time scales. We examined the network dynamics in detail throughout parameter space and showed experimentally that a simple network of oscillators can reproduce the predictions of theoretical models with explicit complex interactions. We fully explained these phenomena by applying a higher-order phase approximation to the full oscillator model. As strong evidence of the symmetry-breaking argument, we delineated the symmetry subgroups associated with each state. CONCLUSION Our results show that simple real-world oscillator networks display complex and exotic system states without the need for complex interactions. Our findings can be applied to observational studies of the behavior of natural or engineered oscillator networks. This work clearly elucidates how a 16-dimensional system, comprising eight magnitude-phase oscillators, collapses into synchronized states that evolve within lower-dimensional subspaces. Real-world networks based on our nonlinear NEMS oscillator platform will enable further insight into the mechanisms by which such behavior emerges in more complex topologies at even larger network scales. States stabilized by emergent interactions. (A) Oscillator network showing the physical connections and the emergent phase interactions. NN, nearest neighbor; NNN, next-nearest neighbor. (B) Stable fixed points for states with no magnitude variations. Colors correspond to colors of interactions from (A) required for stability. (C to E) Experimental heat maps of time-domain data for different combinations of phase differences (using the same data). Exotic states appear as bands in the plots (2-precess, WC-I, 2-TW-I, 2-TW-II).

Abstract Synchronization of oscillators, a phenomenon found in a wide variety of natural and engineered systems, is typically understood through a reduction to a first-order phase model with simplified dynamics. Here, by exploiting the precision and flexibility of nanoelectromechanical systems, we examined the dynamics of a ring of quasi-sinusoidal oscillators at and beyond first order. Beyond first order, we found exotic states of synchronization with highly complex dynamics, including weak chimeras, decoupled states, traveling waves, and inhomogeneous synchronized states. Through theory and experiment, we show that these exotic states rely on complex interactions emerging out of networks with simple linear nearest-neighbor coupling. This work provides insight into the dynamical richness of complex systems with weak nonlinearities and local interactions.

The mutual entrainment of interacting oscillators arises in both natural (1–3) and engineered systems (4–10). The associated phenomena range from simple locked states, in which all oscillators have the same phase, to the general dynamics of inhomogeneous phase configurations associated with pattern formation (11). To explore these phenomena, we used oscillators formed from nanoelectromechanical systems (NEMS), a hybrid of electronic and mechanical degrees of freedom. We exploited the capability of high-speed electronics, yielding individual readout and control, to construct systems with minimal parameter disorder. These innovations allow us to make detailed comparisons with theoretical models.

We examined the dynamics of a ring of eight NEMS limit-cycle oscillators with linear nearest-neighbor coupling. Relative to chemical oscillators (12, 13) or pendulums (7), the time scales over which these dynamics evolve are shorter by at least three orders of magnitude (14); this allows for rapid collection of large datasets that reduce statistical errors and permit emergence of phenomena evolving over very long time scales. This system has allowed us to explore subtle aspects of complex dynamics, given the ability to control the nodes and edges over a large range of parameters (15–17) while simultaneously capturing all degrees of freedom in real time. For example, we could quench the system to examine attractors that are unlikely to be found from random initial conditions.

The foundational model of synchronization is the Kuramoto equation (18, 19), which describes a network of N oscillators with phases ϕ j and frequencies ω j and equations of motion ϕ ˙ j = ω j − ( K / N ) ∑ i = 1 N sin ( ϕ i − ϕ j ) , where K is the strength of the coupling. This model displays dissipative phase dynamics and exhibits a transition to a synchronized state in which a macroscopic fraction of the oscillators evolves at the same frequency. The transition occurs when the inter-oscillator coupling overcomes the disorder of the oscillator frequencies. Kuramoto and Sakaguchi later extended this basic model to include a phase lag θ that induces dispersion, so the sum in the coupling function is over sin(ϕ i – ϕ j + θ) (20).

For weak coupling, as described below, we experimentally found synchronized states with equal oscillator phases or with patterns of fixed phase differences. We present data that can be quantitatively understood from analytical modeling. For stronger coupling, a broad array of stable states emerged showing richer structure and dynamics. We developed experimental and analytical tools to classify and understand these exotic states and demonstrate how they arise in a group-theoretic analysis of the system symmetries (21, 22). Unlike studies of chemical oscillators (23), optoelectronic oscillators (24), or electronic oscillators (25), specialized inter-oscillator or intra-oscillator feedback was not needed to fine-tune network parameters to create these exotic states. By contrast, we show here that quasi-sinusoidal oscillators with linear nearest-neighbor coupling, just beyond the weak coupling limit, are sufficient to manifest exotic states.

Previous studies of the symmetry breaking of oscillator networks focused on explaining exotic dynamics with the use of either phase models with explicit biharmonic phase interactions or complex amplitude models with nonlinear coupling of nodal amplitudes (26–29). Using a model with linear nearest-neighbor coupling of complex amplitudes, we show how an expansion in the coupling parameter yields an approximation for the dynamics of the oscillator phases with emergent biharmonic, next-nearest-neighbor, and triadic phase interactions. This expansion was sufficient to yield almost all of the exotic states we see. The emergent properties are in rough quantitative agreement with experiments and with the full model (which describes complex amplitudes). Our results suggest that exotic symmetry breaking and dynamics are much more prevalent than earlier studies indicate, and should be of direct relevance to biological or socio-technological networks with low connectivity formed from pairwise linear interactions.

Coupled nanoelectromechanical oscillator model Our system can be described by a coupled set of N saturated oscillators (17, 30–32) with equations of motion for the complex amplitudes A j of each oscillator j given by: d A j d T = − A j 2 + A j 2 | A j | + i [ ω j A j + α | A j | 2 A j ] − i β A j + i β 2 ( A j − 1 + A j + 1 ) (1)Equation 1 describes the dynamics on time scales longer than the relaxation time of the resonator, τ slow = Q m /f m , where f m is the resonant frequency of the mechanical cavity and Q m is the quality factor (see supplementary materials). The natural frequency of each oscillator ω j allows for dispersion within our oscillator equations, β is a real number representing the strength of the coupling between nearest neighbors, and α gives the nodal nonlinearity that couples frequency to amplitude. Except where explicitly stated, all of the natural oscillator frequencies are set to be equivalent, making a uniform ring with ω j = 0 in the rotating frame. The variables A j , T and parameters α, β, ω j of Eq. 1 are unitless. The scaled time T is defined in terms of the physical time t via T = 2π × t/τ slow . From this, the relation of α, β to the physical oscillator parameters can be deduced. The oscillator dynamics of Eq. 1 are similar to the dynamics described by the Stuart-Landau equation (33) but with the quadratic nonlinear dissipation in that equation replaced by feedback saturation. Periodic boundary conditions (j = j + N) yield a ring topology. Equation 1 can be divided by A j and separated into real and imaginary parts to give the dynamics of the magnitude |A j | = a j and phase ∠A j = ϕ j , d a j d T = 1 − a j 2 − β 2 [ a j + 1 sin ( ϕ j + 1 − ϕ j ) + a j − 1 sin ( ϕ j − 1 − ϕ j ) ] (2a) d ϕ j d T = ω j + α a j 2 − β + β 2 a j [ a j + 1 cos ( ϕ j + 1 − ϕ j ) + a j − 1 cos ( ϕ j − 1 − ϕ j ) ] (2b)Solutions of this set of equations with a j = 1 are known to be stable (34). For solutions with a j ≠ 1, less is known. For magnitude variations induced by phase dynamics that is slow relative to the magnitude relaxation rate (½ in Eq. 2a), the deviations from a j = 1 can be written to order O ( β ) in terms of the instantaneous phases δa j ≈ –β[sin(ϕ j +1 – ϕ j ) + sin(ϕ j –1 – ϕ j )]. Inserting this into an expansion of Eq. 2b up to first order in the magnitude perturbations δa j , akin to (35), yields d ϕ j d T = ω j + α − β − 2 α β [ sin ( ϕ j + 1 − ϕ j ) + sin ( ϕ j − 1 − ϕ j ) ] ︸ Kuramoto + β 2 [ cos ( ϕ j + 1 − ϕ j ) + cos ( ϕ j − 1 − ϕ j ) ] ︸ Sakaguchi + β 2 4 [ sin ( ϕ j + 2 − ϕ j ) + sin ( ϕ j − 2 − ϕ j ) ] ︸ next-nearest neighbor − β 2 2 { sin [ 2 ( ϕ j + 1 − ϕ j ) ] + sin [ 2 ( ϕ j − 1 − ϕ j ) ] } ︸ biharmonic − β 2 4 sin ( ϕ j + 2 − 2 ϕ j + 1 + ϕ j ) ︸ triadic + − β 2 4 sin ( ϕ j − 2 − 2 ϕ j − 1 + ϕ j ) ︸ triadic − + β 2 2 sin ( ϕ j + 1 − 2 ϕ j + ϕ j − 1 ) ︸ triadic 0 (3)where we label individual coupling terms by their underlying nature for later reference. At O ( β ) this reduces to the Kuramoto-Sakaguchi (K-S) equation. We can rewrite it in a more convenient form (36), d ϕ j d T = ω j + α + K ∑ i = j − 1 , j + 1 { sin ( ϕ i − ϕ j ) + γ [ 1 − cos ( ϕ i − ϕ j ) ] } (4)with K = –2αβ and γ = (4α)−1. Here, γ is a measure of the phase lag θ. For the small coupling regime where Eq. 4 applies, increasing the nonlinear parameter α limits the system to Kuramoto-like dynamics. The dissipative sine coupling in this limit arises from indirect interactions in the phase equation via the nonlinearity α, a parameter derived from the nanoscale physics of the mechanical resonator (30, 37). Our experimental system is unique in that we can independently shift α and β over a wide range spanning two orders of magnitude. Equation 3 does not account for dynamics with rapidly varying magnitudes. If the magnitude dynamics are dominated by a single frequency, for weak coupling we can expand the magnitude perturbations in harmonics. At O ( β 2 ) this introduces reactive, cosine-like next-nearest-neighbor, biharmonic, and triadic terms into Eq. 3. On one hand, this gives quantitative differences at O ( β 2 ) between the phase model, Eq. 3, and the full magnitude-phase model, Eq. 2. On the other hand, if the dissipative sine-like terms dominate stability considerations, qualitative correspondence between Eq. 2 and Eq. 3 for these dynamic situations should be evident. States with static magnitudes should quantitatively agree with Eq. 3. As seen below, we make a comparison between the phase model and experimental data with the inhomogeneous synchronized states. Note that in the full complex amplitude equations, the transformation β → –β, A j → –A j (for all even values of j) gives the same set of equations up to a uniform frequency shift. One consequence is the equivalence between in-phase synchronization for attractive interactions (negative β) and antiphase synchronization for repulsive interactions (positive β). More generally, on the ring topology, the same types of dynamical states exist whether we choose to use an attractive or a repulsive interaction. This symmetry is conserved in the phase approximation, Eq. 3. We study the case of positive β so that antiphase synchronization is expected. We begin by discussing the experimental setup and the system’s symmetries. After this, we examine behavior as the coupling parameter β is increased, starting from the K-S limit where quantitative comparisons can be made. We then divide the discussion of exotic behavior into decoupled states, weak chimeras, and the inhomogeneous synchronized states. To provide additional insight, all of the states presented in this work are shown as animations (movies S1 to S32) using real-time data; see tables S1 to S3 for a guide to the movies.

Experimental setup The experiment is configured to be a controllable network of eight oscillators in which the dynamics of each node are described by a limit cycle in the phase space of mechanical displacement and velocity. Each oscillator was coupled to two neighbors in a chain. We imposed periodic boundary conditions to create the ring topology (Fig. 1A). In an example of the time-domain data (Fig. 1B), we extracted the magnitudes and phase differences of the complex amplitudes A j , with a j = |A j | and Δ j = ϕ j +1 – ϕ j . We show the dynamics of the eight oscillators before, during, and after synchronization into the antiphase state. At fixed α = 0.2, we abruptly changed β from 0.02 to 0.7 within 1 ms to induce synchronization. An animation of the data presented in Fig. 1 is shown in Movie 1. Fig. 1 Experimental setup and representative data. (A) Each oscillator node rotates on a limit cycle (orange on black circle), with coupling edge to nearest neighbors (springs). (B) An example of a synchronization transition. The real-time network dynamics shows the emergence of synchronization as the coupling is initialized. It depicts the oscillator magnitudes a j and the phase differences between nearest neighbors Δ j = (ϕ j +1 – ϕ j ) mod(2π). (C) Phases of the eight oscillators on four 2D tori before and after synchronization. Movie 1. The transition from an unsynchronized state to an antiphase synchronized state after the inter-oscillator coupling is induced (composed of experimental data). If we assume that the magnitudes can be ignored, then the state of the eight oscillators evolves upon an eight-dimensional torus. We show the data from Fig. 1A on four two-dimensional (2D) tori in Fig. 1C. In the top panel, the oscillators were uncoupled and the phase diffused randomly across the tori. In the bottom panel, after coupling was initialized, the system synchronized and its dynamics collapsed to a 1D state space.

Symmetries of the oscillator ring Group theory allows classification of the states of a complex network by considering its symmetry. Ashwin and Swift (21) performed this analysis for rings of weakly coupled oscillators, which applies to our system in most settings, and we follow their notation. Our eight-oscillator ring has the rotational and reflection symmetries associated with an octagon, described by the dihedral group D 8 . Note that Eq. 1 is also invariant under the action of an overall phase shift (represented by the 1-torus T). Combining these symmetries gives the overall symmetry group, D 8 × T. The dynamical states of the system may show reduced symmetries given by the isotropy subgroups of this group. Each isotropy subgroup has an associated subspace that is invariant under its action. The dimension of this invariant subspace gives the number of phase variables required to predict the dynamics (29). Thus, a subgroup with a 1D subspace is described by a single phase, giving a synchronized state with a pattern of locked phases that breaks the symmetry of the ring. A 2D subspace requires two independent instantaneous phases. If these phases remain unlocked, the state will show cluster synchronization, with two clusters of phases locked into a pattern that evolves with different frequencies on a toroidal attractor. Higher-dimensional invariant subspaces typically correspond to more complex dynamics. Table 1 provides the isotropy subgroups of D 8 × T, the dimensions of the invariant subspaces, the group generators, and the associated oscillator patterns. Here, κ is a reflection (with an axis through two nodes), σ is a rotation by one element, and ωp is a phase shift by p × π/4. A “–1” with the generator means that all phases are increased by π in the transformation. The pattern notation is written so that letters a, b, … represent exp(iϕ 1 ), exp(iϕ 2 ), … . As an example, we explain the symmetry represented by the generator (σ3κ, –1) acting on the pattern {a, b, –b, –a, a, b, –b, –a}. First, a reflection about the 2-6 axis gives {–b, b, a, –a, –b, b, a, –a}. Then, there is a rotation by three oscillators counterclockwise, taking the pattern to {–a, –b, b, a, –a, –b, b, a}. By increasing the phase by π, the original pattern is recovered. Table 1 Isotropy subgroups of D 8 × T. View this table:

Splay states The simplest states of identical oscillators coupled on a ring are the states with equal phase differences Δ j = (ϕ j +1 – ϕ j ) mod(2π) along each edge. These phase differences add up to an integral number of phase rotations around the ring given by the winding number k = 1 2 π ∑ j = 1 N Δ j , k = 0 ... N − 1 (5)Note that k = 0 gives the uniform synchronized state and, for even N, the winding number k = N/2 corresponds to the antiphase state. The remaining values of k give splay states (33, 34). Because all of the phase differences are the same for one of these states, we can define Δ(k) ≡ Δ j (k) = 2πk/N for winding number k. From Eqs. 2a and 2b, these states have unit magnitude and the oscillators form a single frequency cluster. For our eight-oscillator ring, the eight unity-magnitude states k = 0, 1, …, 7 comprise the in-phase state (k = 0), the antiphase state (k = 4), and the splay states of rotating phase (k = 1, 2, 3, 5, 6, 7). The in-phase state has all of the spatial symmetries of D 8 . The antiphase state has the same structural symmetry but with twist, D 8 (+, –). The splay states have the symmetry of the twisted cyclic group Z 8 (p), where p ∈ {1,2,3,5,6,7} gives a state with winding number k = {1,2,3,5,6,7}, respectively. Note that the subgroups with {7,6,5} are conjugate to the groups with {1,2,3}. In Movie 2, we show how the twist of the isotropy subgroup relates to a splay state (k = 7). Movie 2. Highlighting the spatiotemporal symmetry of a synchronized splay state, k = 7 (composed of experimental data). For unity magnitudes, Eq. 2b gives a simple prediction for the normalized frequency in the rotating frame Ω(k) = α – β(1 – cos Δ(k)). Note that Eq. 4 exactly captures this result because it is approximated to O ( β ) . Figure 2A shows a plot of the experimental result for the frequency of each state with winding number k against its phase difference (green spheres). We fit this with the equation for Ω(k) and found a difference of 1% in β and an agreement in α better than 0.5%, compared with the values from the initial calibration. The fit is shown as a purple line. Fig. 2 Experimental data on splay states. (A) Frequencies Ω(k) of the unity-magnitude synchronized states k ∈ {0, 1, … 7}, for which the neighboring phase differences are Δ(k) = 2πk/8. At α = 0.074 and β = 1, the data (green spheres) agree well with the fit from the theoretical predictions (purple line). (B) Probabilities of finding the antiphase state from 300 experimental trials for each point with random initial phases, as a function of the phase-lag parameter γ from Eq. 4 (green error bars). We simulated the full set of equations (Eq. 2) with a stochastic white noise term of strength D added to both magnitude and phase. (C) Experimental data on a perturbation of the antiphase state. We perturbed the natural frequency of a single oscillator (oscillator 1) by an amount δω 1 and examined its magnitude a 1 and the deviation of the mean frequency of the ring δω 1 /8, hence δΩ(k=4) – δω 1 /8 is shown. We numerically determined the effect of the perturbation; the plot shows the result from the full model (solid lines) and the Kuramoto-Sakaguchi model (dashed lines). Horizontal error bars are given by digital step resolution of the circuit used for frequency control. Vertical error for the oscillator frequencies is given by the average resolvable frequency over eight oscillators within a ~100-ms time window. For small values of β, we observed the antiphase state k = 4 and the two splay states k = {3,5}. With Δ j = Δ(k), this result is consistent with the linear stability analysis of Eq. 4, which shows that the k = {3,4,5} states are stable for any α, β combination. This analysis also shows that the k = {2,6} states have neutral linear stability (see below). For large β (1.0) and small α (0.074), the in-phase state and k = {1,7} states were also observed (Fig. 2A). Splay states that are not stable within a repulsive K-S model with K < 0—that is, k = {0,1,7}—are stable within simulations of the extended phase model (Eq. 3). From simulations of Eq. 3, it turns out that only the K-S and biharmonic terms were required to stabilize the in-phase state. Likewise, only the K-S and triadic terms in Eq. 3 were required to stabilize k = {1,7}. Thus, different terms in the extended phase model stabilize different states. In the absence of noise, the probability to settle into a particular synchronized state from a random initial set of phases can be calculated from the volumes of the basins of attraction of the different solutions (38). In Fig. 2B, we show the probability of synchronizing into the antiphase state as a function of phase-lag parameter γ in our experiment. The probability of finding this state in the ring with an even number of repulsively coupled oscillators is equivalent to finding the in-phase state in the same size ring with attractively coupled oscillators. Thus, our results can be directly compared to previous results for the Kuramoto limit γ = 0 (38). We ran the experiments at a fixed value of coupling K = −0.2 (Fig. 2B, green error bars) and swept γ from 0.04 to 0.75, starting with 300 different random initial conditions for each value. At γ ≈ 0.75, β was increased to a point where magnitudes varied strongly and states other than the splay states with winding number k = 3, 4, 5 started to stabilize. We plotted the ratio of the number of trials ending in an antiphase state (k = 4) to the total number of trials. The error is given by the uncertainty generated by sampling the probabilities with only a finite number of trials. At smaller values of γ, the probabilities appear to approach a constant value before sharply deviating from this trend for γ < 0.07. Fitting the data for γ > 0.07 to a quadratic function of γ2 (Fig. 2B, red curve) yielded the γ → 0 intercept of p(k=4)(γ = 0) = 0.76 ± 0.02. This result proved consistent with previous numerical simulations on larger rings in the Kuramoto limit by Wiley et al. (38) when extrapolated to our eight-oscillator ring. We associate the behavior for γ < 0.07 with effects of noise. To confirm this, we performed stochastic numerical simulations of Eq. 2 (with the oscillators having identical natural frequencies ω j ) by adding independent stochastic forces ξ j (t) to both the magnitude and phase equations. We used a two-step stochastic Runge-Kutta scheme for additive Gaussian white noise (39), with 〈ξ i (t)ξ j (t′)〉 = Dδ(t – t′)δ ij , and determined the network’s state at the end of these simulations, carrying out more than 1000 trials. The results (Fig. 2B) show that the large-γ behavior was effectively noise-independent, whereas the small-γ behavior was strongly noise-dependent. The simulations without noise for γ = 0 agree with the result extrapolated from Wiley et al. (38). Although we do not expect the noise in our experimental system to be completely white (16), the experimental data for γ < 0.07 were roughly consistent with the results of numerical simulations with a noise strength D ~ 4 × 10–4. Magnitude noise fed into the phase equation via the nonlinear frequency tuning for small values of γ. Nonlinear frequency tuning is characterized by the parameter α = 1/(4γ); it enhances the effective phase noise strength at low frequencies by a factor of (1 + γ­–2). Thus, at low γ, phase noise became much stronger, and noise-induced switching between attracting states was much more likely. We studied the robustness of the synchronized state to dispersion in the natural oscillator frequencies by suddenly perturbing the system after it had settled into the antiphase state. To accomplish this, we shifted the natural frequency of the first oscillator by an amount δω 1 and increased the size of the perturbation until the oscillator broke loose and synchronization was destroyed. For a weak perturbation, linear response gave a frequency shift of the array equal to the perturbation of the single oscillator divided by the number of oscillators in the array, δΩ(4) = δω 1 /8. In general, this result broke down for large perturbations approaching the limit of synchronization. It remains true, however, in the Kuramoto limit even for large perturbations. Figure 2C shows the frequency shift of the array after subtracting off the small-perturbation prediction of δω 1 /8 for three different values of K-S parameter γ. A constant residual frequency shift of δΩ = 0.005 common to all three datasets was also subtracted. This could be accounted for by drift in fast oscillator frequency δf m /f m ~ 10−6 over the experimental time scale of ~100 ms, which was chosen to acquire data with error ~0.001. In addition to the mean frequency, we report the magnitude variable of the perturbed oscillator. Numerical solutions of Eqs. 2 and 4 are plotted as solid and dashed lines, respectively. We compare the experiment to the fixed-point solutions found in a numerical solver for both the full set of magnitude-phase equations (Eq. 2) and the K-S equation. The experimental uncertainty in the x axis originated from the uncertainty in the frequency shift caused by the discreteness of the digital-to-analog converter used to implement control over nodal frequencies. When oscillator 1 desynchronized, the system did not settle into a fixed point, and the solver failed to find a valid solution, so we included no numerical data after desynchronization. In this case, we have manually modified the parameters α and β (reported as K and γ, respectively, for the K-S solution) in the solver to obtain a match to the data of Fig. 2C. Although good agreement between experiment and theory was attained, more disagreement was found in this case than for the frequency data from the uniform ring (Fig. 2C versus Fig. 2A). This difference may be due to stray nonlinear coupling arising from parametric actuation of the mechanical resonators (see supplementary text). Note that even at β = 0.1, the solution to the K-S model shows some deviations from the solution of the full complex amplitude equations. It is noteworthy that we found quantitative agreement with theory in a high-dimensional nonlinear system, even for subtle phenomena such as response near the desynchronization threshold, where linear response breaks down. Also, through careful control of frequencies, we were able to attain a nearly uniform network. This enabled us to test analytically intractable problems, such as quantifying the state probabilities arising from random initial conditions in the presence of noise. Estimating these requires the calculation of the volume of the basins of attraction. These results demonstrate that NEMS oscillators enable general and detailed studies of synchronization.

Inhomogeneous synchronization In addition to the aforementioned frequency-synchronized splay states exhibiting identical phase differences between nearest-neighbor oscillators, we discovered frequency-synchronized states with a spread of phase differences. We accessed these states by taking advantage of our experimental capabilities that permit nonadiabatic movement through the parameter space by quickly changing the values of α, β, and ω j while carrying out the simultaneous readout of the states of all eight oscillators. The inhomogeneous synchronized states we explored in this case arose from the pattern-forming instability of the spatially uniform in-phase state. This instability is a Fourier mode of phase perturbations, which grows and ultimately saturates to form a stationary state with periodic spatial modulation of nearest-neighbor phase differences. We expect that the period of the spatial pattern that develops from the uniform state is associated with the wave vector maximizing the growth rate from linear stability (eq. S5). For the eight-oscillator ring and values of α, β we used, this yielded the inhomogeneous synchronized pattern of phases with a single spatial period. Figure 3A shows an example of the patterns that emerged from the uniform state, where the phase differences are depicted as a function of time. We prepared the system without natural frequency disorder ω j ≈ 0, and then switched on the coupling to β = 1 at t ≈ 0.3 s. At t ≈ 0.53 s, we then suddenly shifted the coupling to β = 0.28. At these parameter values (α = 0.09, β = 0.28), the in-phase state was calculated to be unstable, and the phases subsequently exploded into a patterned, inhomogeneous state. This state was frequency-synchronized because the phase differences were constant in time. However, a large spread in the phase differences between nearest neighbors appeared, larger than could be accounted for by residual natural frequency dispersion in ω j . In Fig. 3B, we show the growth of the spatial pattern as a function of time; the color of the slice corresponds to the time axis. After the quench (green dotted line), the spatial pattern grew into a rough sinusoidal shape. As shown, the patterns exhibited a single spatial period, fulfilling the linear stability prediction. In Movie 3, we show the pattern formation through the quench shown in Fig. 3B. Fig. 3 Inhomogeneous synchronization. (A) Oscillator nearest-neighbor phase differences showing pattern formation. At t ≈ 0.53 s, we suddenly shifted the coupling from a value of β = 1 to β = 0.28. (B) Phase differences of the inhomogeneous state as the system is quenched from the data in (A). The color corresponds to the time axis, with more blue representing later times. (C) Numerical data on the state space of the inhomogeneous synchronization. The plot shows the magnitude of the Kuramoto order parameter Z ip characterizing the steady state found in these simulations. (D) Data showing the overlap of the uniform and inhomogeneous synchronized states as a function of coupling parameter β at α = 0.09. The plot shows the mean (in trial number) of the absolute value of the order parameter with a threshold for each trial given by Δ ¯ j < 0.02 and M n < 0.94 (blue stars) or M n > 0.94 (maroon triangles). A clear separation of the states appeared when this threshold was applied. The experimental results are compared with the numerical simulations of Eq. 2 (black dots) and Eq. 3 (green circles). Movie 3. Formation of the spatial pattern of phases in the inhomogeneous synchronized state through a quench of the coupling parameter β (composed of experimental data). We compared the experimental data to the results of numerical simulations, both for the full magnitude-phase model (Eq. 2) and for the second-order phase model (Eq. 3). Figure 3C displays the results of a deterministic simulation (a simulation without noise) of Eq. 1, starting the system at the fixed point a j = 1 and Δ j = ξ j , where ξ j is a pseudorandom number sampled from a uniform distribution over the interval (–π/20, π/20). We plotted the resulting absolute value M of the order parameter Z ip = M exp(iΦ) = ( 1 / N ) ∑ j = 1 N exp ( i ϕ j ) over the parameter space defined by α ∈ [0.04, 0.18] and β ∈ [0.1, 1]. We observed three regions: a region of high synchrony (“uniform in-phase”), a region of reduced order parameter (“inhomogeneous synchronized”), and a region where neither the uniform nor inhomogeneous synchronized states is stable (“other synchronized states”). The black region represents values of M < 0.5. In addition, this region did not have a near-zero phase difference (calculated on the interval –π, π). The boundary between the uniform and inhomogeneous regions was consistent with the threshold of linear stability α = β sin2(π/8) (34) of the in-phase state. That is, at values to the lower right of this line, the in-phase state was stable and M = 1 (shown in white). To obtain quantitative comparison with theory, we performed the following experiment. First, we set the system to (α, β) = (0.09, 1) to initialize it near the fixed point (a j , Δ j ) = (1, ξ j ) with ξ j < 0.05, and then reinitialized it until the in-phase state was detected. Second, we rapidly quenched the system down to a lower value of β, along the dashed line in Fig. 3C. We applied this procedure, stepping β between 0.3 and 0.9 with an increment of 0.03, and repeating it 20 times for each step. Occasionally, the system was observed to randomly switch into neither the uniform nor the inhomogeneous state. Accordingly, we preselected states with a mean phase difference Δ ¯ j < 0.002 when averaged over both space and time. Each trial was binned into one of two groups by thresholding at M = 0.94, which represents the greatest value for the states found in simulation along the dashed line in Fig. 3C. From each of these two bins, we plotted the average of M over the trials against β in Fig. 3D. We also show the results from numerical simulations of the full magnitude-phase model (Fig. 3C) and the approximate phase model (Eq. 3) as black and green dots, respectively. We found that the phase model agreed with the magnitude-phase model for values of β < 0.3. The in-phase state for both models was found to destabilize at the same value of coupling. The experimental data show qualitative agreement with the numerical results (black dots) of the full model—for example, manifesting a discontinuous transition at β ~ 0.6. In this plot, we could not readily explain the systematic quantitative disagreement between simulations of the full model (Eq. 2) and our experimental data. Synchronized states with a spread of phase differences have been discussed by Pikovsky and Rosenblum (27) in networks of Stuart-Landau oscillators with global nonlinear coupling. By contrast, for our system, such spread occurred without need for explicit nonlinear coupling. Also, no spatial pattern of the phases was observed in their previous work because simulations were performed in a globally connected topology.

Decoupled states The following sections describe states that can form multiple frequency-synchronized clusters. To find the symmetries that are preserved, we had to introduce several measures that identify frequency clustering. This is illustrated in Fig. 4, A and B, which depicts (from top to bottom) several neighboring phase differences; space-time plots of oscillator frequencies that identify the frequency clusters; pairwise mutual information of the phases within the state I(ϕ i ; ϕ j ) (40); and experimental snapshots of the oscillator phases. In these lowermost panels, each oscillator node is colored to identify members of a specific synchronized node cluster, arrow orientation depicts each node’s relative instantaneous phase, symmetry operations corresponding to the group generators are shown as dashed arrows, and the corresponding isotropy subgroup is shown at the lower right. Fig. 4 Decoupled states. (A) Phase differences, frequencies, mutual information matrix, and a snapshot of the 2-precess state (see text). Note that the next-nearest neighbors (e.g., oscillators 1 and 3) are in the same group and are antiphase-locked. (B) Phase differences, frequencies, mutual information matrix, and a snapshot of the 2-TW-I state. (C) Mean frequency difference of the clusters (within the decoupled state) as a function of the difference of the mean natural frequencies of the clusters. The data show that for two different values of coupling, the intercept is near zero. Thus, in a ring without disorder, no frequency difference between clusters exists in the decoupled state. (D) The in-phase synchronized state, the 2-precess state, and the 2-TW-I state are plotted in a reduced phase space of Δ 8 , Δ 1 , a 1 . Decoupling occurs for the special case of a ring containing 4m oscillators (where m is a positive integer) with nearest-neighbor linear interactions (21): In states with A j +1 = –A j –1 , the oscillator A j becomes decoupled from its neighbors (see Eq. 1). In our system, for the k = 2 splay state, Δ j +1 + Δ j = π, which gives A j +1 = –A j –1 . A consequence of the decoupling is that the k = 2 and k = 6 splay states became two points on a continuous line within state space, corresponding to an arbitrary phase difference between the even- and odd-numbered clusters. For a nonideal system with a nonzero difference in the average natural frequency within a cluster, 1 / 4 ∑ j = 1 4 ω 2 j ≠ 1 / 4 ∑ j = 1 4 ω 2 j − 1 , we expect to observe two clusters precessing at different frequencies. This is depicted as what we term the 2-precess state in Fig. 4A, first row, where two neighboring phase differences drift. Note that the decoupling condition, Δ j +1 + Δ j = π, is satisfied at all times. The second row of Fig. 4A depicts the frequencies as a function of time. The frequency clusters show a difference of ~20 Hz (0.04 in normalized units). The third row demonstrates that the mutual information is appreciable between next-nearest neighbors but is negligible between nearest neighbors. Two phase variables were required to describe this state. Here, evidently, symmetry was not completely broken and four-fold rotational symmetry persisted. The symmetry subgroup associated with this state is the twisted cyclic group Z 4 (p = 2). In Movie 4, we directly show the precession of the two groups of oscillators within this decoupled state. Movie 4. The 2-precess state with two sets of decoupled oscillators (composed of experimental data). The magnitudes within the observed 2-precess state manifested small, slow variations. Thus, the phase model (Eq. 3) should provide a good approximation to the state’s dynamics. Indeed, we found that this 2-precess state was reproduced in our simulation of the phase model. In these simulations, we varied the coefficients of the higher-order phase interactions and thereby could identify the nearest-neighbor coupling term as most important for this state’s stability. The biharmonic and triadic terms appeared to have little effect on this state, although they smoothed the trajectories of the phase differences. We show these simulations in fig. S6A. Figure 4C demonstrates that the two clusters of the decoupled state have identical frequencies in the ideal ring. We extracted data from the time-domain records of the natural frequencies ω j and the oscillator frequencies within the state Ω j (as in the lower panel of Fig. 1C before and after t = 0.25 s, respectively). We plotted the difference between the frequencies of the oscillator clusters, 1 / 4 ∑ j = 1 4 Ω 2 j − 1 / 4 ∑ j = 1 4 Ω 2 j − 1 , against the difference of the means of their natural frequencies, 1 / 4 ∑ j = 1 4 ω 2 j − 1 / 4 ∑ j = 1 4 ω 2 j − 1 . At two values of β, the intercept nearly vanished. Thus, when the difference of the mean of the natural frequencies was near zero, the frequency of two clusters became nearly indistinguishable. We note that the type of decoupling we observed experimentally—a network of D 4 m symmetry with linear, pairwise coupling—was previously discovered theoretically by Alexander and Auchmuty (41). The 2-precess state was not the only state that was found near the k = 2 and k = 6 fixed points. We also observed two traveling-wave states that oscillate around each of the k = 2 and k = 6 fixed points (Fig. 4B). These states were not observed in simulations of the original set of equations (Eq. 2). However, we were able to reproduce these states in subsequent simulations of the complex amplitude equation (Eq. 1) in which we included additional nonlinear nearest-neighbor coupling proportional to A j + 1 2 and A j − 1 2 (specifically using eq. S9). These additional terms likely arose from harmonics at 2f m ~ 4.4 MHz, which could serve to parametrically excite nearest neighbors. The traveling waves can manifest a spatial wave vector of π/2 (the 2-TW-I or first traveling-wave state around the k = 2 state shown in Fig. 4B) or π/4 (see fig. S2A). For the 2-TW-I state, the symmetry subgroup is Z 2 . Finally, in Fig. 4D we display the decoupled and traveling-wave states of Fig. 4, A and B, together with the in-phase k = 0 state in a 3D plot. These demonstrate how a single oscillator magnitude varies according to its two neighboring phase differences. The in-phase state is ideally a point within this phase space and hence can be described by a single phase variable. However, the other two states are lines in this cross section of phase space and need at least two phase variables to be described. Note that in the absence of noise, the 2-precess and the 2-TW-I states would not intersect each other. In other words, the volumes in phase space where these states can reside will not overlap. However, with experimental noise, we observed some degree of overlap of the respective regions of phase space; in this situation, the system could indeed transition between these two states.

Discussion We demonstrated that a simple ring of eight self-sustained nanoelectromechanical oscillators with linear, nearest-neighbor coupling exhibits exotic states of synchronization with complex dynamics and broken symmetries. These states include weak chimeras, decoupled states, traveling waves, and inhomogeneous synchronized states. We developed a theoretical formulation starting from the equations of motion for the oscillators that yielded an approximation for the dynamics of the oscillator phases with emergent higher-order interactions including biharmonic, next-nearest neighbor, and triadic interactions. More generally, we showed that emergent higher-order interactions arise in simple networks of weakly nonlinear elements with linear, nearest-neighbor coupling. We found excellent agreement between our model describing the coupled NEMS oscillators (Eq. 1) and the experimental data for splay states, including the frequencies of the various states. More subtle issues were also elucidated, such as the response to perturbing the natural frequency of a single oscillator until desynchronization occurs and the probability of falling into a particular state from random-phase initial conditions. At high values of coupling, we also found approximate quantitative agreement for one of the weak chimera states. Qualitatively, all states (except the decoupled traveling wave) of the system are reproduced in our simulations based on both the full model (Eq. 2) and a reduced phase model (Eq. 3). These results demonstrate the promise of performing detailed quantitative tests of synchronization in large arrays of NEMS, even far from the weak-coupling limit. We categorized the observed states of the system in terms of their symmetries. Of the 18 isotropy subgroups within the oscillator ring’s D 8 × T symmetry group, we found distinct states associated with 11 of them. These include the in-phase, antiphase, and splay states associated with a K-S ring. Also included are complex states that require more than one phase variable for their description. Although higher-dimensional states have been analyzed previously with phase oscillators that interact through exotic couplings, we found similar behaviors in magnitude-phase oscillators with simple linear nearest-neighbor interactions. We found that the emergent second-order phase interactions stabilizing the higher-dimensional states provide intuition about the correspondence between our work and previous work that has modeled phase oscillators using far more complicated coupling schemes. Our discussion on the robustness of the antiphase state (Fig. 2C) demonstrated control over individual nodes. Our calibration procedure for ensuring identical nodes [see supplementary materials and (30)] and edges can be modified to induce explicit symmetry breaking. This itself constitutes a viable strategy for network control (47, 48). Our results illustrate the promise of nanomechanical oscillator arrays for exploring control of complex and large-scale networks. Combined with the full access to detailed real-time dynamics of the state variables we have achieved, feedback control is an exciting possibility with this singular experimental system. Larger arrays of NEMS oscillators (~100) are possible with minimal changes to the system architecture discussed herein. These results reveal the diversity of exotic, self-organized patterns of synchronization that can occur in simple oscillator networks. Practically, they also provide insights for real-world systems where such dynamically self-organized patterns of oscillation can be either disruptive or beneficial. For instance, circulating power flows have been observed in the North American high-voltage electric grid. Driven by oscillating generators, power flows have been observed to unintentionally self-organize into closed loops around geographically constrained regions, subsuming transmission line capacity without delivering useful power to consumers (49, 50). Coherent patterns of local and cluster synchronization can also lead to instabilities in models of the electric grid (22, 51). Likewise, even in 1D models of automobile traffic flow, small disturbances often lead to dynamically self-organized soliton pulses of congestion (52–54). Indeed, spontaneous symmetry breaking in the form of chimeras (43) is under intense investigation in neuronal dynamics (55, 56) and has been implicated as a mechanism for Parkinson’s tremor, epilepsy, and unihemispheric slow-wave sleep (57). Further, new patterns of synchronization from self-organized dynamics exploiting differences between the attributes of agents (58) will unravel the consensus processes across engineering applications, including sensor networks and robot swarms (59). More broadly, our system will help with the design of intrinsically resilient technological networks at multiple scales, including rapidly evolving energy systems (60).

Methods summary The assembled network comprised individual oscillators constructed from high-quality piezoelectric nanoelectromechanical membrane resonators, which were clamped on all sides and driven via electronic feedback. Nodal parameters were controlled via electronic circuitry assembled on printed circuit boards (PCBs) onto which the NEMS resonators themselves were bonded. Details of the individual nonlinear NEMS oscillator nodes were presented by Fon et al. (30). These individual oscillator PCB nodes were inserted into a larger motherboard that formed the coupling edges (connections) between network nodes. Control signals from a Raspberry Pi were routed to the motherboard and multiplexed to individual oscillators by means of digital addressing. Oscillator signals were read out by coaxial connections on the motherboard to achieve radio-frequency bandwidths. Power was fed into the oscillators via the motherboard. All experiments, numerical simulations, and data processing were performed with custom Python and C scripts. Signals were acquired through use of a simultaneous-sampling 8-channel oscilloscope. Oscillator phase was extracted via a Hilbert transform of raw time records from the oscilloscope. Oscilloscope sampling was set to 200 kHz, which is faster than the slow-time oscillator dynamics by about three orders of magnitude. See supplementary materials for further details.

Supplementary Materials www.sciencemag.org/content/363/6431/eaav7932/suppl/DC1 Materials and Methods Supplementary Text Figs. S1 to S6 Tables S1 to S3 Movies S1 to S32 References (61–64)

http://www.sciencemag.org/about/science-licenses-journal-article-reuse This is an article distributed under the terms of the Science Journals Default License.

Acknowledgments: We thank D. Abrams for fruitful comments regarding the traveling-wave states, and J. E. Sader and R. Lifshitz for useful comments regarding the manuscript. We also thank CEA-LETI (Grenoble, France) for providing piezoelectric multilayers enabling this research. We acknowledge critical support and infrastructure provided for this work by the Kavli Nanoscience Institute at Caltech. Funding: This material is based on work supported by, or in part by, the U.S. Army Research Laboratory and the U. S. Army Research Office under MURI award W911NF-13-1-0340 and W911NF-18-1-0028 and Intel Corporation support of CSC as an Intel Parallel Computing Center. Author contributions: M.H.M. and J.L. fabricated the nanomechanical devices; M.H.M., W.F., and M.L.R. designed and constructed the experimental apparatus; M.H.M. wrote the control/measurement software and performed the measurements; M.H.M. analyzed the experimental data with input from J.E., W.F., M.C.C., J.P.C., and M.L.R.; M.H.M. composed and narrated all videos with input from J.P.C.; J.E., M.H.M., A.C., A.S., and M.C.C. performed theoretical modeling with input from R.M.D., J.P.C., M.P., and M.M.; J.E., M.H.M., A.S., and A.C. performed numerical simulations with input from J.P.C., M.R., M.H.d.B., M.L.R., M.M., M.C.C., M.P., and R.M.D.; M.H.M. and M.L.R. conceived the experiment with input from W.F., M.C.C., J.E., A.C., M.M., L.D.-O., R.M.D., and J.P.C.; and M.H.M., M.C.C., W.F., J.P.C., J.E., R.M.D., and M.L.R. prepared the manuscript with input from all authors. All authors discussed the results and their implications equally at all stages. Competing interests: None declared. Data and materials availability: All (other) data needed to evaluate the conclusions in the paper are available in the following online database: 10.5281/zenodo.2543765.