Network experiment demonstrates converse symmetry breaking
NNetwork experiment demonstratesconverse symmetry breaking
Ferenc Molnar, Takashi Nishikawa, , Adilson E. Motter , Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL 60208, USA
Symmetry breaking—the phenomenon in which the symmetry of a system isnot inherited by its stable states—underlies pattern formation, superconduc-tivity, and numerous other effects. Recent theoretical work has establishedthe possibility of converse symmetry breaking (CSB), a phenomenon in whichthe stable states are symmetric only when the system itself is not. This in-cludes scenarios in which interacting entities are required to be non-identicalin order to exhibit identical behavior, such as in reaching consensus. Here wepresent an experimental demonstration of this phenomenon. Using a networkof alternating-current electromechanical oscillators, we show that their abilityto achieve identical frequency synchronization is enhanced when the oscilla-tors are tuned to be suitably non-identical and that CSB persists for a range ofnoise levels. These results have implications for the optimization and controlof network dynamics in a broad class of systems whose function benefits fromharnessing uniform behavior.
Nature Physics , 351–356 (2020) The published version is available online at: http://dx.doi.org/10.1038/s41567-019-0742-y a r X i v : . [ n li n . AO ] S e p ynchronization —perhaps the most widely studied phenomenon in network dynamics —has been observed in many contexts, including both natural systems (e.g., circadian clockcells , ecological populations
8, 9 , human menstrual cycles , and crowds of pedestrians ) andengineered systems (e.g., Boolean logic gates , semiconductor lasers , electrochemical andnanomechanical oscillators , and power generators ). Such observations are significantbecause they show that approximately homogeneous dynamics can emerge in heterogeneouspopulations. Yet, until recently, the prevailing view had been that homogeneity in the dynamicsis facilitated by increased homogeneity in the population. This view has now changed withthe theoretical discovery that, in numerous systems, heterogeneity can be required for stableidentical synchronization—even when the entities are identically coupled to the population.The underlying phenomenon, which we term converse symmetry breaking , can be elegantlydescribed using the notion of symmetry—a fundamental property that can characterize a systemand has deep implications for its dynamics . In contrast to the well-known phenomenon ofsymmetry breaking, in which symmetry in the system implies broken symmetry in the stablestates, converse symmetry breaking represents a scenario in which symmetry in the stable statesimplies broken symmetry in the system. The interaction networks of many real systems areinvariant under node permutations and hence possess symmetries . Symmetry breaking innetworks include important examples of chimera states , in which a broken-symmetry statewith coexisting groups of synchronized and non-synchronized nodes is observed even thoughthe system is symmetric. Converse symmetry breaking, on the other hand, has been predictedfor oscillator networks in which the phenomenon can be mediated, for example, by amplitudedynamics , couplings internal to the oscillators , and interaction delays . However, unlikesymmetry breaking, evidence for converse symmetry breaking has thus far remained theoretical.In this Article, we present the first experimental demonstration of converse symmetry break-ing, in which we account for noise and other realistic features. Our experimental system isdesigned to allow for frequency synchronization and consists of alternating current (AC) elec-tromechanical oscillators, which are identically coupled in order to isolate the effect of oscillatorheterogeneity from that of coupling heterogeneity. We show that, within the precision of exper-imental measurements, the optimal stability of frequency synchronization can be enhanced bymaking the values of a tunable parameter of the oscillators—the damping coefficient—suitablydifferent from each other. Our results indicate that we can harness converse symmetry breakingin optimizing network dynamics. In potential applications to networks in which synchroniza-tion is desirable, this would translate to controlling oscillator heterogeneity to enable enhanced2 peedControllers12VDCSupplyComputerInterfaceControl andData CollectionElectricLoads BrushedDC MotorsMechanicalBrakesSynchronousGenerators ∞
12 3 b cd e a ReflectiveInfrared Sensors
Fig. 1:
Experiment involving a network of coupled electromechanical oscillators. a , Main compo-nents of the experimental setup, including three AC generators, three DC motors driving them, and thecomputerized data acquisition system. b , Diagram of the AC electrical circuit connecting the three gen-erators, running at Hz. c , Network representation of the circuit, where the nodes represent generatorsand the links represent the electrical interactions between them. The parameters characterizing the nodesand links are normalized by suitable references (i.e., given in per-unit quantities). d , Predicted stability ofthe frequency-synchronous splay states as a function of the generator parameters β i . Inside the coloredsurface is the region of stability given by λ max < λ max th for noise corresponding to λ max th = − . . e ,Cross section of the stability landscape at the plane shown in d . Color-coded is the value of λ max relativeto λ max th . The optimal uniform assignment ( (cid:101) β ) and the globally optimal non-uniform assignment ( β g ) aremarked by red and yellow crosses, respectively. Also marked by crosses are the projections of β A and β B , the nearby assignments that we realize experimentally. stability and performance.Figure 1a illustrates the main components of our experiment, in which three permanent-magnet generators are mechanically driven by DC motors with adjustable speed and a separate12V DC power supply. The generators are chosen to have identical parameters (e.g., internaldamping coefficient, internal impedance, and terminal voltage at various speeds) within manu-facturing precision; see Methods for details. To allow for heterogeneous configurations of thegenerators, their shafts are equipped with mechanical brakes that can be used to adjust fric-tion. The generators’ output is connected to a set of electric loads (inductors and capacitors)forming the circuit depicted in Fig. 1b. The parameters of the circuit components are cho-sen for the system to be symmetric with respect to rotational permutations of the generators( → → → ). The pattern of coupling among the generators can thus be represented as arotationally symmetric network of three nodes (generators) connected by three identical links,3ach corresponding to a circuit equivalent to the so-called π model . Each generator is fixedon its own platform to minimize mechanical coupling with the other generators via vibrations.We focus on frequency-synchronous states of the system in which the voltage frequenciesof the generators are all equal to a constant frequency ω s . On short timescales, of the order of sor less, the dynamics of the generators when the system is close to such a state can be describedby a coupled oscillator model
44, 45 . When written for an arbitrary number n of generators, themodel equation reads: ¨ δ i + β i ˙ δ i = a i − (cid:88) k (cid:54) = i c ik sin ( δ i − δ k − γ ik ) + εξ i ( t ) , (1)where δ i is the internal electrical angle for generator i (a state variable related to the rotor shaftangle by a factor determined by the number of poles in the generator), relative to a referenceframe rotating at the synchronous frequency ω s ; the constant β i is an effective damping param-eter (capturing both mechanical and electrical damping, normalized by the generator’s inertia); a i is a parameter representing the net power accelerating the generator’s rotor (i.e., the me-chanical power provided by the DC motor driving the generator, minus the power consumed bythe network components and the power lost to damping); c ik and γ ik are the coupling strengthand phase shift characterizing the electrical interactions between the generators; ξ i ( t ) is a ran-dom function representing dynamical noise; and ε is the noise amplitude (see SupplementaryInformation, Sec. S1 for a derivation of the deterministic part, Methods for details on model-ing dynamical noise, and Supplementary Fig. 1 for validation). In the co-rotating frame, thefrequency-synchronous states correspond to the fixed-point solutions of Eq. (1) with ε = 0 ,characterized by ω i ≡ ˙ δ i = 0 , ∀ i . Note that uniform angle shifts of one such solution representthe same state, as all angle differences remain unchanged. The deterministic part of Eq. (1) hasthe same form as that used to model the dynamics of generators in power grids
20, 21 , which haverecently been studied extensively in the network dynamics community , but here the equa-tion is used to model coupled electromechanical oscillators that are tunable and not constrainedto operating states of power grids.Linearizing Eq. (1) for ε = 0 around a fixed point with δ i = δ ∗ i , we obtain ˙ x = Jx , where x = (cid:0) ∆ δ ∆ ω (cid:1) and J = (cid:0) − P − B (cid:1) . Here, ∆ δ and ∆ ω are the n -dimensional vectors of angle andfrequency deviations, δ i − δ ∗ i and ω i − ω i = ˙ δ i , respectively. The n × n matrix P = ( P ik ) is given by P ik = (cid:40) c ik cos( δ ∗ i − δ ∗ k − γ ik ) , i (cid:54) = k, − (cid:80) k (cid:48) (cid:54) = i P ik (cid:48) , i = k, (2)4 is the n × n diagonal matrix with β i as its diagonal elements, and and denote the nulland identity matrices of size n , respectively. The stability of the corresponding frequency-synchronous state of Eq. (1) with ε = 0 is thus determined by the eigenvalues λ i of the Ja-cobian matrix J , excluding the identically zero eigenvalue (which we denote by λ ) presentonly because of the zero row-sum property of P . Specifically, if the maximal real part of theseeigenvalues is negative, i.e., the Lyapunov exponent λ max ≡ max i ≥ Re( λ i ) is negative, then thestate is asymptotically stable, and smaller λ max implies stronger stability. (The zero eigenvalue λ is excluded because it is associated with perturbations that uniformly shift phases, whichlead to equivalent fixed points corresponding to the same frequency-synchronous state.) Theproblem of maximizing this stability with respect to system parameters can thus be formulatedas the optimization of λ max . This is similar to the problem of optimizing the largest real partof the eigenvalues of a matrix, known as the spectral abscissa (relevant when the matrix has noidentically null eigenvalues) , and the problem of optimizing the second largest eigenvalueof a Laplacian matrix, known as the algebraic connectivity
54, 55 . Some previous studies haveconsidered problems concerning the optimization of damping parameters
52, 53 , but they do notfocus on relations to system symmetry and typically exclude the class of non-positive definitenon-symmetric coupling matrices relevant to the experiment considered here.In the presence of dynamical noise (i.e., when ε > ), we can show that, for general classesof discrete- and continuous-time noise models, there is a (negative) threshold λ max value forthe stability of the frequency-synchronous state. We denote this stability threshold by λ max th ; thefrequency synchronization is stable below this threshold and unstable above it (see Methods fordetails).We now use our stability analysis to derive a condition for observing converse symmetrybreaking in this system. Note that λ max is a function of β = ( β , . . . , β n ) , and that, even thoughEq. (1) does not depend on β i when restricted to the synchronous states ˙ δ i = 0 , the corre-sponding variational equation ˙ x = Jx does. This implies that making β i heterogeneous, whichbreaks the symmetry of the system, generally leads to symmetry breaking in the eigenmodesaround the synchronous state (Fig. 2). We see that, if two values of β correspond to distinctvalues of λ max , then there is a range of noise intensities (and thus of λ max th ) for which the stateis stable for one β value and unstable for the other (see Methods for details). Thus, a conditionfor exhibiting converse symmetry breaking is that there exists β ∗ representing a non-uniform β i assignment for which λ max ( β ∗ ) < min { , λ max ( (cid:101) β ) } , (3)5 -15 -10 -5 time time n o r m a li z e d ( d e g ) time n o r m a li z e d ( d e g ) a bc slopeslope r e l a t i v e p e r t u r b a t i o n s i z e , Fig. 2:
Oscillator heterogeneity breaks the symmetry of the dominant eigenmodes.
The predictedslowest decaying eigenmodes of deviations from the splay state are visualized for the system corre-sponding to Fig. 1b–e. a , Magnitude (cid:107) ∆ δ (cid:107) of the eigenmodes, which decay at an overall exponentialrate λ max ( β ) , for β = (cid:101) β (red) and for β = β g (blue). For β = (cid:101) β , the decay shown is for a combi-nation of two oscillatory eigenmodes associated with complex conjugate eigenvalues corresponding to λ max ( (cid:101) β ) , whereas for β = β g , the decay is for a combination of three eigenmodes corresponding to λ max ( β g ) , of which two oscillatory modes are associated with complex conjugate eigenvalues and onenon-oscillatory mode is associated with a real eigenvalue. b , Dynamics of individual oscillators (orange,green, and purple for oscillator i = 1 , , , respectively) given by the eigenmodes for β = (cid:101) β , after nor-malization that removes the exponential decay shown in a . (Left) Dynamics with respect to time. (Right)Amplitude of eigen-oscillations indicated by the bars drawn along the unit circle (normalized such that (cid:107) ∆ δ (cid:107) = 0 . ). The gray dots in the background represent the splay state, in which δ i for different i are120 ◦ apart from each other. c , Same as in b , but for β = β g . In this case, the bars for perturbationamplitudes are shifted by offsets corresponding to the non-oscillatory eigenmode. We observe that thedominant eigenmodes are rotationally symmetric for β = (cid:101) β , but this symmetry is broken for β = β g .For an animated version of the plots in b and c , see https://youtu.be/R_BOIWXYtSk where we define (cid:101) β ≡ ( (cid:101) β, . . . , (cid:101) β ) to represent the (uniform) β i assignment that minimizes λ max under the constraint that β = · · · = β n .In the design of our experiment, the electrical parameters of the AC circuit (indicated inFig. 1b) were chosen to ensure that the circuit is rotationally symmetric and the frequency-synchronous states that inherit that symmetry have λ max < at frequency ω s = 100 Hz (seeMethods for details). The specific states we focus on are known as splay states , which for aring of n phase oscillators are defined as states in which phase differences between consecutive6scillators are all equal to π/n . With the choice of parameters used in the experiment, Eq. (1)describes a reduced -node network system with a i = a , c ik = c , and γ ik = γ (values given inFig. 1c), which has -fold rotational symmetry (as the original AC circuit). This system has twosplay states (with identical λ max ), corresponding to two fixed points of Eq. (1) for which thephase angles of the adjacent generators are exactly degrees apart (i.e., δ − δ = ± ◦ and δ − δ = ∓ ◦ ). Each of these is a rotationally symmetric state of the system, since applyingany rotational permutation of the generators would result in the same state (equivalent up to auniform shift of all δ i ). In these states, the effective coupling pattern expressed in the matrix P takes a form that favors converse symmetry breaking (see Methods for details).The stability landscape—the function λ max ( β ) = λ max ( β , β , β ) —for the frequency-syn-chronous splay states of our experimental system exhibits the same rotational symmetry as thesystem itself, as illustrated in Fig. 1d assuming ideal system components (see Methods for de-tails). Due to this symmetry, the landscape has global minima (all with the same λ max value)at three different locations in the β -space, related by a -degree rotation of the β i -axes. Oneof these minima, which we denote by β g , is located approximately at (3 . , . , . andcorresponds to a globally optimal (non-uniform) assignment of the generator parameters with λ max ( β g ) ≈ − . . In contrast, the optimal uniform assignment (i.e., the optimum underthe constraint β = β = β ) corresponds to the point (cid:101) β = ( (cid:101) β, (cid:101) β, (cid:101) β ) with (cid:101) β ≈ . and λ max ( (cid:101) β ) ≈ − . . These two points are marked on the cross section of the stability landscapeshown in Fig. 1e. Since condition (3) is satisfied for β ∗ = β g , the system is predicted to ex-hibit converse symmetry breaking: the rotationally symmetric, frequency-synchronous splaystates are necessarily unstable for any uniform β i assignment (corresponding to a rotationallysymmetric system), whereas they become stable for the non-uniform β i assignment in β g (cor-responding to a non-rotationally symmetric system) for a range of noise intensities.Specific uniform and non-uniform β i assignments were implemented in our experiment byadjusting the frictional brakes on the generator shafts (see Supplementary Information, Sec. S1on how friction relates to β i ). Due to the physical limitations and finite measurement accuracy,the precision of a β i value that we could realize was ± . (see Supplementary Information,Sec. S2.1). In addition, changes in the parameters of generator–motor units over time due to ex-ternal effects (e.g., heating) as well as inherent heterogeneity of the system components (despitebeing manufactured to be identical) can distort the landscape itself and thus shift the locationsof β g and (cid:101) β . Even though care was taken to minimize deviations from the designed values (seeMethods), it would be experimentally challenging to realize these points exactly. Therefore,7e considered two points on the landscape that we were able to realize and confirm with mea-surements: β A ≡ (3 . ± . , . ± . , . ± . ) ≈ (cid:101) β and β B ≡ (3 . ± . , . ± . , . ± . ) ≈ β g (projections of both points are marked in Fig. 1e). Our theoretical predictions for the stabilityat these points are λ max ( β A ) ≈ − . and λ max ( β B ) ≈ − . , which confirms the property λ max ( β B ) < λ max ( β A ) that is predicted to enable converse symmetry breaking in the experi-ment.To provide experimental evidence for converse symmetry breaking, we performed multipleexperimental runs while measuring the terminal voltage of each generator with a voltage sensorcircuit and concurrently processing the measurements (see Methods for details). These mea-surements yielded multiple time series of terminal voltage phasors (angles and magnitudes) andthe corresponding electrical angles δ i , both recorded at , samples per second. The totalrecorded time span was ≈ . h each for the β A and β B configurations. For our analysis, wefocused on segments of the time series in which the splay states were observed (allowing fora maximum of ± degree deviation in δ i ; see Methods for the full details on this criterion).Since system parameters may fluctuate and shift gradually during an experimental run, we es-timated a (slightly different) steady state for each segment using the averages of the measuredphasor angles δ ∗ i over all data points in the segment. Using these steady-state values, as well asmeasured system parameters, we estimated λ max for each segment of each β configuration (seeMethods for the procedure to calculate λ max ).The λ max estimated from experimental data (shown in Fig. 3a) are distributed around thecorresponding theoretical predictions for the splay states. The mean values of λ max for the twoconfigurations are significantly different, which we confirmed with the paired t -test; the null-hypothesis that the difference of the means is zero is rejected, because the p -value is smallerthan machine precision and is certainly smaller than the significance level of . . The differ-ence between the two means ( ≈ . ) is substantially larger than a typical variation of λ max along the trajectory in a segment (Supplementary Fig. 2a), indicating that our conclusion is notsensitive to uncertainty in the measurement of steady-state δ ∗ i values.In addition to establishing the statistically significant difference in the averaged λ max be-tween the two configurations, we investigated the extent to which changing the configurationfrom β A to β B enhances the stability of the same experimentally realized states. This was doneby recomputing λ max for β = β A and β = β B for each time-series segment, which is justifiedbecause the steady state of Eq. (1) does not depend on the choice of β . Thus, the recomputed λ max determines the stability we would observe if, starting with one β configuration, the brakes8ig. 3: Experimental confirmation of converse symmetry breaking. a , Distribution of the Lyapunovexponents λ max obtained from individual time-series segments for the uniform ( β A , orange histogram)and non-uniform ( β B , blue histogram) configurations. For each segment, we used the measured splaystate and associated parameters, and calculated λ max as an average over that segment. The downwardarrows indicate the predicted values of λ max for the theoretically calculated splay states (in which theangles are exactly 120 degrees apart). b , Distribution of synchronization stability improvement achievedby changing the generator parameters from β A to β B for the same states, as measured by the difference λ max ( β A ) − λ max ( β B ) . c , Inferred splay states and density plot of time-series trajectories. The horizontaland vertical axes are the phase angles δ and δ of the second and third generators, respectively, relativeto the first generator. The trajectory density was estimated for each pixel using all measured time series.The color scale is normalized to the highest measured density. The dots in the enlarged sections markthe steady splay states determined for the time-series segments we considered (275 and 190 segments forthe β A and β B configurations, respectively). on the generators were adjusted to realize the other configuration without changing the state.Finding λ max ( β B ) < λ max ( β A ) < for a given steady state implies that condition (3) for con-verse symmetry breaking is satisfied for β ∗ = β B and (cid:101) β ≈ β A . As shown in Fig. 3b, thiswas indeed observed to be the case for almost all the identified time-series segments (whosecorresponding splay states inferred from data are shown in Fig. 3c). We also verified that theobserved stability improvement is robust against uncertainties in generator parameters (Sup-plementary Fig. 2b). Thus, our measurements and analysis provide experimental evidence ofconverse symmetry breaking—a rotationally symmetric synchronous state of our system be-comes stable when the rotational symmetry of the system is broken by the heterogeneity of thegenerator parameters.Our demonstration of converse symmetry breaking can be interpreted from various differ-ent angles. On the one hand, it establishes a scenario in which stable symmetric states requiresystem asymmetry. On the other hand, it shows that the assumption that increasing uniformityacross individual entities would facilitate uniform behavior is generally false, even when they9re indistinguishable in terms of their interactions. In network systems, this equates to statingthat there are situations in which the nodes need to be suitably non-identical in order for themto stably converge to identical dynamical states even if all nodes occupy structurally equivalentpositions in the network. This leads to the counter-intuitive conclusion that node heterogene-ity across a network can help—rather than inhibit—convergence to a homogeneous dynamicalstate, as required for numerous processes including synchronization, consensus, and herding.But what enables converse symmetry breaking? Symmetry breaking itself can be imme-diately appreciated by considering a one-dimensional Mexican hat potential, where the statewith reflection symmetry (the system’s symmetry) is unstable and the two stable states areasymmetric. This simple example illustrates the fundamental tenet that a symmetric model candescribe asymmetric observations. Conversely, the phenomenon demonstrated here shows thatsymmetric observations may require an asymmetric model. While potentially less immediate tovisualize, which might be the reason it was not demonstrated earlier, converse symmetry break-ing can be interpreted as arising from the following trade-off: system asymmetry can reduce thelikelihood of having a symmetric solution, but it can also increase the likelihood of having onesuch solution stable. This realization opens the door for new control approaches to manipulatesystem parameters and optimize the stability of symmetric states in networks whose functionbenefits from the symmetry of these states. While we purposely designed our experiment withidentically coupled oscillators to isolate the phenomenon, the opportunities for optimizationand control are even more evident when the identical coupling constraint is lifted, since thestabilizing effect of breaking the system symmetry with node heterogeneity is expected to becommon in such cases. We thus suggest that the results presented here will naturally extend toreal systems with tunable node parameters, such as networks of logic gates, neuronal systems,coupled lasers, and networks of mechanical, electrical, and chemical oscillators. Methods
Experimental design.
The AC generators used in the experiment are of . V and . W be-cause they are easily available and provide the desirable low-voltage, low-power output. Whilelarger generators could synchronize with less noise as they are machined with higher precisionrelative to their sizes, they require power electronics, and leave less room for errors that couldcause physical damage. To ensure that the generators are as identical as possible, we obtained10ight generators from the same manufacturer, and selected three of them with the closest ter-minal voltage for a given rotational speed. The procedure we used to measure the generatorparameters is detailed in Supplementary Information, Sec. S2.We chose a specific AC frequency, ω s = 100 Hz, which equals / ≈ . rotations persecond for the generator shaft (our generators have three pole pairs). The generators provide . V terminal (r.m.s.) voltage at this speed (without load), and we verified that the voltage wasproportional to speed, up to ≈ rotations per second. However, we found that the terminalvoltage dropped to ≈ . V when the generator was connected to the network with a typicalconfiguration we considered. We thus used . V as the base (a reference value) for normalizingvoltages in all per-unit calculations. The impedance base ( .
48 Ω ) was determined in the processof making a choice of electrical component parameters, as described below, and the power base( . W) was determined accordingly from the chosen voltage and impedance bases.For any given combination of capacitance, inductance, and resistance that make up a rota-tionally symmetric circuit of the same form as in Fig. 1b (and our choice of ω s = 100 Hz), thetheoretical prediction of λ max for the splay states for any β can be computed using the proce-dure described in the section ‘Calculation of λ max ’ below. The specific combination shown inFig. 1b was identified through multiple iterations of system design in which we made refine-ments in our choice of electrical components, construction of the experimental modules, andmeasurement apparatus. The refinement of component parameters was performed using theirnormalized, per-unit values (for example, the capacitance was chosen such that the shunt sus-ceptance of the π model representing each link is ≈ . in per unit). Ultimately, we chose aparameter combination in which the predicted stability of the splay states is sufficiently strongfor (cid:101) β and improves significantly when changing to β g , while ensuring that system componentsoperate within their rated capacities and the β change is experimentally realizable by frictionadjustment. The actual inductors ( . H and . , according to the manufacturer’s specifi-cation) and capacitors ( µ F) were selected from those that were readily available to closelymatch the values in Fig. 1b.
Modeling dynamical noise.
We analyze the stability of Eq. (1) in the presence of the noiseterm εξ i ( t ) , described either by a discrete-time model or a continuous-time model. In thediscrete-time model, εξ i ( t ) is modeled as random impulse perturbations (i.e., as a sum of Diracdelta functions with random magnitudes located at random times). We can account for anydistributions of perturbation magnitudes and times that are bounded in the sense that there is11 maximum magnitude M for the impulse perturbations and a minimum time interval τ be-tween consecutive perturbations (which corresponds to a maximum rate at which the systemis perturbed). When the system is close to a splay state with maximum Lyapunov exponent λ max < , the dynamics approximately follows the linearized equation (as described in themain text) and is characterized by an exponential convergence to the splay state at a rate of λ max between any two consecutive perturbations. Thus, if the deviation of the system statefrom the splay state immediately after the k th perturbation at time t k is ∆ x ( t k ) , the deviationafter the next perturbation is given by ∆ x ( t k +1 ) = e λ max ∆ t k · ∆ x ( t k ) + ξ k +1 , where ∆ t k isthe (random) time interval between the k th and ( k + 1) th perturbations, and ξ k is the (random)displacement resulting from the k th impulse perturbation. Since ∆ t k ≥ τ and (cid:107) ξ k (cid:107)≤ M , wehave (cid:107) ∆ x ( t k +1 ) (cid:107) ≤ e τλ max (cid:107) ∆ x ( t k ) (cid:107) + M . By recursively applying this inequality, we obtain (cid:107) ∆ x ( t k ) (cid:107) ≤ ( e τλ max ) k (cid:107) ∆ x ( t ) (cid:107) + M · − ( e τλ max ) k +1 − e τλ max → M − e τλ max as k → ∞ , (4)which indicates that the deviation from the splay state is bounded by M/ (1 − e τλ max ) in thelimit of large k (and thus large t ). Using (cid:107) ∆ x ( t k ) (cid:107) ≤ ∆ sync , ∀ k , where ∆ sync is a constant, wesee that the system would stay synchronized in the splay state if M − e τλ max < ∆ sync . (5)It follows that the splay state is stable in the presence of noise if the noise magnitude M issufficiently small to satisfy this condition. On the other hand, if the asymptotic bound M/ (1 − e τλ max ) is larger than ∆ sync , there is a non-zero probability that the deviation (cid:107) ∆ x ( t k ) (cid:107) exceeds ∆ sync for sufficiently large k , implying that the splay state is unstable. Therefore, given M , τ ,and ∆ sync , there is a (negative) threshold λ max value λ max th ≡ τ ln(1 − M/ ∆ sync ) (6)corresponding to the stability transition for the noisy system: the splay state is stable in thepresence of noise if λ max < λ max th and unstable if λ max > λ max th . Note that each value of λ max th represents a range of combinations of M , τ , and ∆ sync . In general, for any pair of configurations β I and β II with λ max ( β I ) < λ max ( β II ) < , there is a range of combinations of M , τ , and ∆ sync for which λ max ( β I ) < λ max th < λ max ( β II ) , i.e., the splay state is stable for the β I configuration,while it is unstable for the β II configuration. 12n the continuous-time model, the term εξ i ( t ) represents the effect of Brownian noise onthe dynamics of the i th oscillator. Since the effect of the noise decays most slowly along theeigenmode associated with the eigenvalue λ max < , we focus on the projection of the full six-dimensional dynamics onto that eigenmode. We model this projected dynamics by the followingstochastic differential equation: dx ( t ) = λ max x ( t ) dt + σdξ ( t ) , (7)where x ( t ) represents the (one-dimensional) deviation from the frequency-synchronous state, ξ ( t ) is the standard Brownian noise (i.e., ξ ( t ) − ξ (0) for fixed t follows the Gaussian distributionwith zero mean and variance t ), and σ > is a constant that can be used to tune the level ofnoise intensity felt by the system. The stochastic process given by Eq. (7) is called the Ornstein–Uhlenbeck (OU) process
57, 58 . Given the constants σ and ∆ sync , we define Pr ( λ max ) to be theprobability that | x ( t ) | ≤ ∆ sync for t → ∞ (i.e., the probability that the system is synchronized inthe large t limit). From the known property of the OU process that x ( t ) for a fixed t follows theGaussian distribution with mean e tλ max x (0) and variance σ (1 − e tλ max ) / ( − λ max ) , we computePr ( λ max ) = erf (∆ sync √− λ max /σ ) , where erf denotes the error function. Note that, as λ max < is increased (and thus √− λ max is decreased), the probability Pr ( λ max ) decreases. Thus, for agiven constant p th (chosen close to one), we define the stability threshold λ max th as the value of λ max for which Pr ( λ max ) = p th . It then follows that λ max th = − σ [ erf − ( p th )] / ∆ sync . Note thatthe threshold λ max th = λ max th ( σ ) depends on the noise intensity level σ . For σ = 0 (noiselesscase), we have λ max th (0) = 0 , recovering the stability threshold for deterministic systems. As σ increases, the threshold λ max th ( σ ) monotonically decreases, indicating that the decay of deviationneeds to be faster to maintain the same probability Pr ( λ max ) for higher noise levels. Therefore,if λ max ( β I ) < λ max ( β II ) < , there is a range of σ for which λ max ( β I ) < λ max th ( σ ) < λ max ( β II ) ,i.e., the system stays near the frequency-synchronous state with high probability as t → ∞ forthe β I configuration, while it does not for the β II configuration. Matrix P for splay states. Assuming ideal system components, the electric circuit we use(Fig. 1b) is rotationally symmetric, and thus the complex node-to-node admittances Y ik areidentical for all i (cid:54) = k . If, in addition, the internal impedances z int ,i of the three generators areidentical, then the effective admittances Y ik are also identical, i.e., Y ik = y = | y | exp( jα ) forall i (cid:54) = k , where j is the imaginary unit. Assuming further that the inertia constants H i are allidentical, the steady splay states satisfy c ik = c for all i (cid:54) = k , and the matrix P in Eq. (2) takes13he form − b − b (cid:48) b b (cid:48) b (cid:48) − b − b (cid:48) bb b (cid:48) − b − b (cid:48) , (8)where b = c (cid:0) − cos γ ∓ √ sin γ (cid:1) , b (cid:48) = c (cid:0) − cos γ ± √ sin γ (cid:1) , and γ = α − π/ . Thetopology of the underlying network corresponds to one favoring converse symmetry breakingin the analysis of coupled oscillators presented in Ref. 41. Calculation of λ max . For both the theoretically predicted stability landscape in Fig. 1 andthe estimates from the experimental measurements in Fig. 3 and Supplementary Fig. 2, wecompute λ max using the following procedure. We first apply the Kirchoff law to compute theactive and reactive power injections at the terminal of each generator ( P i and Q i , respectively)from the network’s inductances and capacitances (at a given frequency) as well as the terminalvoltage magnitudes and angles ( | V i | and θ i , respectively). We then apply the Kirchoff law againto the internal impedance of each generator (assuming the classical model) to compute | E i | and δ ∗ i from P i , Q i , | V i | , and θ i . From these and the values of H i , the parameters c ik and γ ik can be computed. We can then obtain λ max for any given β by computing the eigenvalues of J = (cid:0) − P − B (cid:1) with the matrix P defined in Eq. (2). For the theoretical prediction, we use ω s = 100 Hz and the terminal voltage magnitude | V i | = 1 . V. We set the angles θ i to haveexactly ◦ differences, which ensures splay-state values for δ i as we assume the generatorparameters to be identical: r int ,i = 1 . and x int ,i = 3 .
34 Ω for all i . We compute H i from J i = 5 × − kg · m , which makes H i all identical. For the experimental estimation of λ max from a given time-series segment, we use the average of the measured values over the segmentfor ω s , | E i | , and δ ∗ i , along with the values of J i , r int ,i , and x int ,i measured for each generator (seeSupplementary Information, Secs. S2.2 and S2.3 for our measurement procedure). Conducting the experiment and taking measurements.
Arduino Uno and Arduino Megamicrocontrollers were used to provide sufficient computational capacity and I/O channels tofacilitate measurement and control of the experiment. The Arduino Uno has six analog inputs,of which three were used to measure the terminal voltages of the three generators, and theremaining three were used to measure the currents passing through the generator terminals. Thesame controller also has several digital outputs, four of which were used to control relays thatswitch the voltage inputs between measurement signals and test signals. The test signals wereused to calibrate the voltage levels. We used three analog inputs on the Arduino Mega to take14haft rotational speed measurement from infrared sensors that detect markings on the shafts.This capability was used in measuring β i for each generator (see Supplementary Information,Sec. S2.1 for details). In addition, three digital outputs of this controller were used to driverelays that switch the generators on and off in the AC circuit.Each experimental run was performed as follows. First, while the generator terminal re-lays were open (i.e., disconnected from the circuit), we turned on the DC motors driving thegenerators, and we adjusted their speed manually to / ≈ . rotations per second, yield-ing terminal voltage signals with frequencies approximately at the desired Hz (noting thatthe frequency is three times higher than the shaft’s rotational frequency because our generatorshave three pole pairs). Once this was achieved, the generator terminal relays were closed, es-tablishing the coupling between the generators. This can result in a transient voltage dynamicswith large fluctuations, with the speed of each generator dropping due to the impedance of theconnected network. The speed of each DC motor was adjusted once again to ensure that eachgenerator was running at the desired speed. When the last generator approached the desiredspeed, the synchronous state formed spontaneously. Real-time voltage phasor readouts weredisplayed on the computer interface to confirm that the system had achieved a desired splaystate. At this point, we started recording the voltage phasor readings.External perturbations from vibrations and other sources continuously disturb the synchro-nous state, and the resulting transients prevented us from observing this state for more than afew seconds at a time. Moreover, as a result of components heating up during the experiment,the amount of friction between the shaft and its housing tends to change at different rates fordifferent generators, leading to unbalanced net changes in the mechanical power input of thegenerators that distort the splay state. We mitigated this problem by re-adjusting the speed ofthe DC motors. These adjustments are implemented manually, aided by the real-time display ofthe phase angle differences of the generators. To prevent damage due to rising temperature ofthe components, we were limited to a maximum of 10 min of continuous run when no brakeswere applied (the β A configuration), and a maximum of 5 min when brakes were applied (the β B configuration).The terminal voltage was recorded at , samples per second for each generator usingthe analog-to-digital converters of the Arduino Uno, along with the microsecond-accurate timestamps from the microcontroller’s internal timer. The data were streamed to the computer inter-face for post-processing in which the raw readings were converted to time-dependent frequencyand phasor angle and magnitude with the original resolution. This conversion was done us-15ng custom software for phasor measurement, which implements the following calculations.For each data point, we took the last three oscillatory periods of the voltage signal and calcu-lated the least-squares fit of a sinusoidal function with its amplitude, offset, phase shift, andfrequency as the fitting parameters. The resulting amplitude, phase shift, and frequency wereused as estimates of the terminal voltage peak magnitude, phasor angle, and frequency, respec-tively, for that data point. To eliminate measurement noise from the terminal voltage phasorand frequency time series, we applied a rd -order Savitzky-Golay filter , with a window sizeof . s. The generators’ internal voltage phasors E i = | E i | exp( jδ i ) were then obtained us-ing the Kirchoff laws, along with the calculated terminal voltage phasors, the instantaneoussynchronous frequency (computed as the average of the instantaneous AC frequencies of thegenerator terminals), and each component’s capacitance and inductance. Identifying time-series segments of splay states.
Because some deviation from the exactsplay state is unavoidable, the steady-state phase angles δ ∗ i are not necessarily separated byexactly degrees. We thus identified in each time series the set of maximal segments thatsatisfy the following criteria: 1) for each data point in the segment, one generator is ahead byan angle ∆ + and another is behind by ∆ − relative to the other generator (taken to be generator1 here), while satisfying | ∆ ± − ◦ | < ◦ ; and 2) the length of the segment is at least . s (approximately ten cycles of voltage). In total, we have obtained 275 segments for the β A configuration and 190 segments for the β B configuration. Supplementary Fig. 3 showsthe length distribution of these segments, which are within the timescale for which Eq. (1) isvalid. Within these segments, we verified that the generator frequencies were synchronized:the maximum (instantaneous) frequency difference among the three generators was < Hz,and the standard deviation of each generator’s frequency was < . Hz in % of the
465 =275 + 190 identified segments. We also verified that all parameters in the deterministic part ofEq. (1) were constant within experimental noise (see Supplementary Information, Sec. S1 andSupplementary Fig. 1 for details). This further confirmed the validity of Eq. (1) for describingthe δ i dynamics in each time-series segment. Data Availability.
All data that support results in this article are available from the corre-sponding author upon reasonable request. 16 ode Availability.
The custom code used for the analysis of the data from the experiment isavailable from the corresponding author upon reasonable request.
References Pikovsky, A., Rosenblum, M. & Kurths, J.
Synchronization: A Universal Concept in Nonlin-ear Sciences (Cambridge University Press, Cambridge, 2001). Arenas, A., D´ıaz-Guilera, A., Kurths, J., Moreno, Y. & Zhou, C. Synchronization in complexnetworks.
Phys. Rep. , 93–153 (2008). D¨orfler, F. & Bullo, F. Synchronization in complex networks of phase oscillators: A survey.
Automatica , 1539–1564 (2014). Pecora, L. M. & Carroll, T. L. Synchronization of chaotic systems.
Chaos , 097611 (2015). Yamaguchi, S. et al. Synchronization of cellular clocks in the suprachiasmatic nucleus.
Sci-ence , 1408–1412 (2003). Yamaguchi, Y. et al. Mice genetically deficient in vasopressin V1a and V1b Receptors areresistant to jet lag.
Science , 85–90 (2013). Lu, Z. et al. Resynchronization of circadian oscillators and the east-west asymmetry of jet-lag.
Chaos , 094811 (2016). Ranta, E., Kaitala, V., Lindstr¨om, J. & Linden, H. Synchrony in population dynamics.
Proc.R. Soc. Lond. B , 113–118 (1995). Schwartz, M. K., Mills, L. S., McKelvey, K. S., Ruggiero, L. F. & Allendorf, F. W. DNAreveals high dispersal synchronizing the population dynamics of Canada lynx.
Nature ,520–522 (2002). McClintock, M. K. Menstrual synchrony and suppression.
Nature , 244–245 (1971). Strogatz, S. H., Abrams, D. M., McRobie, A., Eckhardt, B. & Ott, E. Theoretical mechanics:crowd synchrony on the Millennium Bridge.
Nature , 43–44 (2005). Rosin, D. P., Rontani, D., Gauthier, D. J. & Sch¨oll, E. Control of synchronization patterns inneural-like Boolean networks.
Phys. Rev. Lett. , 104102 (2013).17 Fischer, I. et al. Zero-lag long-range synchronization via dynamical relaying.
Phys. Rev. Lett. , 123902 (2006). Zamora-Munt, J., Masoller, C., Garc´ıa-Ojalvo, J. & Roy, R. Crowd synchrony and quorumsensing in delay-coupled lasers.
Phys. Rev. Lett. , 264101 (2010). Argyris, A., Bourmpos, M. & Syvridis, D. Experimental synchrony of semiconductor lasersin coupled networks.
Opt. express , 5600–5614 (2016). Kiss, I. Z., Zhai, Y. & Hudson, J. L. Emerging coherence in a population of chemical oscilla-tors.
Science , 1676–1678 (2002). Kiss, I. Z., Rusin, C. G., Kori, H. & Hudson, J. L. Engineering complex dynamical structures:sequential patterns and desynchronization.
Science , 1886–1889 (2007). Fon, W. et al. Complex dynamical networks constructed with fully controllable nonlinearnanomechanical oscillators.
Nano Lett. , 5977–5983 (2017). Hill, D. J. & Chen, G. Power systems as dynamic networks. In
Proceedings of the 2006 IEEEInternational Symposium on Circuits and Systems , 722–725 (2006). Motter, A. E., Myers, S. A., Anghel, M. & Nishikawa, T. Spontaneous synchrony in power-grid networks.
Nat. Phys. , 191–197 (2013). D¨orfler, F., Chertkov, M. & Bullo, F. Synchronization in complex oscillator networks andsmart grids.
Proc. Natl. Acad. Sci. USA , 2005–2010 (2013). Nishikawa, T. & Motter, A. E. Symmetric states requiring system asymmetry.
Phys. Rev. Lett. , 114101 (2016). Okuda, K. & Kuramoto, Y. Mutual entrainment between populations of coupled oscillators.
Prog. Theor. Phys. , 1159–1176 (1991). Golubitsky, M., Stewart, I., & Schaeffer, D. G.
Singularities and groups in bifurcation theory ,Vol. 2 (Springer-Verlag, New York, 1988). Golubitsky, M., & Stewart, I. Symmetry and pattern formation in coupled cell networks. In
Pattern Formation in Continuous and Coupled Systems , 65–82 (Springer-Verlag, New York,1999). 18 Nicosia, V., Valencia, M., Chavez, M., D´ıaz-Guilera, A. & Latora, V. Remote synchronizationreveals network symmetries and functional modules.
Phys. Rev. Lett. , 174102 (2013). Pecora, L. M., Sorrentino, F., Hagerstrom, A. M., Murphy, T. E. & Roy, R. Cluster synchro-nization and isolated desynchronization in complex networks with symmetries.
Nat. Com-mun. , 4079 (2014). Whalen, A. J., Brennan, S. N., Sauer, T. D. & Schiff, S. J. Observability and controllabilityof nonlinear networks: the role of symmetry.
Phys. Rev. X , 011005 (2015). Sorrentino, F., Pecora, L. M., Hagerstrom, A. M., Murphy, T. E. & Roy, R. Complete char-acterization of the stability of cluster synchronization in complex dynamical networks.
Sci.Adv. , e1501737 (2016). Zhang, L., Motter, A. E. & Nishikawa, T. Incoherence-mediated remote synchronization.
Phys. Rev. Lett. , 174102 (2017). Cho, Y. S., Nishikawa, T. & Motter, A. E. Stable chimeras and independently synchronizableclusters.
Phys. Rev. Lett. , 084101 (2017). Barrett, W., Francis, A. & Webb, B. Equitable decompositions of graphs with symmetries.
Linear Algebra Appl. , 409–434 (2017). MacArthur, B. D., S´anchez-Garc´ıa, R. J. & Anderson, J. W. Symmetry in complex networks.
Discrete Appl. Math. , 3525–3531 (2008). Kuramoto, Y. & Battogtokh, D. Coexistence of coherence and incoherence in nonlocallycoupled phase oscillators.
Nonlinear Phenom. Complex Syst. , 380–385 (2002). Abrams, D. M. & Strogatz, S. H. Chimera states for coupled oscillators.
Phys. Rev. Lett. ,174102 (2004). Panaggio, M. J. & Abrams, D. M. Chimera states: coexistence of coherence and incoherencein networks of coupled oscillators.
Nonlinearity , R67 (2015). Hagerstrom, A. M. et al. Experimental observation of chimeras in coupled-map lattices.
Nat.Phys. , 658–661 (2012). 19 Tinsley, M. R., Nkomo, S. & Showalter, K. Chimera and phase-cluster states in populationsof coupled chemical oscillators.
Nat. Phys. , 662–665 (2012). Martens, E. A., Thutupalli, S., Fourri`ere, A. & Hallatschek, O. Chimera states in mechanicaloscillator networks.
Proc. Natl. Acad. Sci. USA , 10563–10567 (2013). Hart, J. D., Bansal, K., Murphy, T. E. & Roy, R. Experimental observation of chimera andcluster states in a minimal globally coupled network.
Chaos , 094801 (2016). Zhang, Y., Nishikawa, T. & Motter, A. E. Asymmetry-induced synchronization in oscillatornetworks.
Phys. Rev. E , 062215 (2017). Zhang, Y. & Motter, A. E. Identical synchronization of nonidentical oscillators: when onlybirds of different feathers flock together.
Nonlinearity , R1 (2017). Grainger, J. & Stevenson, W.
Power System Analysis . (McGraw-Hill Co., Singapore, 1994). Anderson, P. M. & Fouad, A. A.
Power System Control and Stability (IEEE Press, Piscataway,2003). Nishikawa, T. & Motter, A. E. Comparative analysis of existing models for power-grid syn-chronization.
New J. Phys. , 015012 (2015). Susuki, Y., Mezi´c, I. & Hikihara, T. Coherent swing instability of power grids.
J. NonlinearSci. , 403–439 (2011). Lozano, S., Buzna, L. & D´ıaz-Guilera, A. Role of network topology in the synchronizationof power systems,
Eur. Phys. J. B , 231–238 (2012). Menck, P. J., Heitzig, J., Kurths, J. & Schellnhuber, H. J. How dead ends undermine powergrid stability.
Nat. Commun. , 3969 (2014). Auer, S., Kleis, K., Schultz, P., Kurths, J. & Hellmann, F. The impact of model detail onpower grid resilience measures.
Eur. Phys. J. Special Topics , 609–625 (2016). Sch¨afer, B., Beck, C., Aihara, K., Witthaut, D. & Timme, M. Non-Gaussian power gridfrequency fluctuations characterized by L´evy-stable laws and superstatistics.
Nat. Energy ,119–126 (2018). 20 Burke, J. V., Lewis, A. S. & Overton, M. L. A robust gradient sampling algorithm for nons-mooth, nonconvex optimization.
SIAM J. Optimiz. , 751–779 (2005). Freitas, P. & Lancaster, P. On the optimal value of the spectral abscissa for a system of linearoscillators.
SIAM J. Matrix Anal. A. , 195–208 (1999). Kirillov, O. N. & Overton, M. L. Robust stability at the swallowtail singularity.
Frontiers inPhysics , 1–9 (2013). Boyd, S. Convex optimization of graph Laplacian eigenvalues. In
Proc. ICM , 1311–1319(2006). De Abreu, N. M. M. Old and new results on algebraic connectivity of graphs.
Linear AlgebraAppl. , 53–73 (2007). Zou, W. & Zhan, M. Splay states in a ring of coupled oscillators: from local to global cou-pling.
SIAM J. Appl. Dyn. Syst. , 1324–1340 (2009). Uhlenbeck, G. E. & Ornstein, L. S. On the theory of the Brownian motion.
Phys. Rev. ,823–841 (1930). Gillespie, D. T.
Markov Processes: An Introduction for Physical Scientists (Academic Press,Inc., San Diego, 1991). Doob, J. L. The Brownian movement and stochastic equations.
Ann. Math. , Vol. 43, 351–369(1942). Orfanidis, S. J.
Introduction to Signal Processing (Prentice Hall, Englewood Cliffs, 1996).21 cknowledgements
The authors thank John B. Ketterson for insightful discussions about this research. This researchwas funded by ARO Grant No. W911NF-15-1-0272 and by Northwestern University’s FiniteEarth Initiative (supported by Leslie and Mac McQuown).
Author contributions
F.M., T.N., and A.E.M. designed the research and contributed to the modeling. F.M. performedthe experiments and simulations. F.M., T.N., and A.E.M. analyzed the results and wrote thepaper. All authors approved the final manuscript.
Competing interests
The authors declare no competing interests.
Materials & Correspondence
Correspondence and material requests should be addressed to T.N. or A.E.M.(E-mails: [email protected] or [email protected]).22 upplementary Information
Network experiment demonstrates converse symmetry breaking
S1 Derivation of coupled oscillator network model
Our derivation of the deterministic part of the model in Eq. (1) of the main text is based onthe so-called classical model of a generator. For completeness, we first reproduce a derivationof the classical model starting from Newton’s second law written in terms of the total torqueaccelerating the rotor of generator i : J i ¨ φ m ,i = − d i ˙ φ m ,i + T m ,i − T e ,i , (S1)where φ m ,i is the mechanical (rotor shaft) angle of generator i (relative to a stationary axis), J i isthe total moment of inertia of the rotor, T m ,i is the mechanical torque provided to the rotor, T e ,i isthe electrical torque (load), and d i is the damping-torque coefficient accounting for windage andfriction. Changing to a frame of reference rotating at the synchronous (mechanical) frequency ω sm of the rotor through (cid:101) φ m ,i ≡ φ m ,i − ω sm t , Eq. (S1) becomes J i ¨ (cid:101) φ m ,i = − d i ( ˙ (cid:101) φ m ,i + ω sm ) + T m ,i − T e ,i . (S2)Multiplying this by the angular velocity ω m ,i ≡ ˙ φ m ,i , we can write ω m ,i J i ¨ (cid:101) φ m ,i = − ω m ,i d i ( ˙ (cid:101) φ m ,i + ω sm ) + ω m ,i T m ,i − ω m ,i T e ,i = − ω m ,i d i ( ˙ (cid:101) φ m ,i + ω sm ) + P m ,i − P e ,i , (S3)where P m ,i is the mechanical power supplied to the rotor and P e ,i is the electrical power drawnfrom the rotor. When the system is close to a frequency-synchronous state, we have ω m ,i ≈ ω sm ,and we can write ω sm J i ¨ (cid:101) φ m ,i = − ω sm d i ˙ (cid:101) φ m ,i + P m ,i − ω d i − P e ,i . (S4)Thus, in a synchronous steady state in which ˙ (cid:101) φ m ,i = ¨ (cid:101) φ m ,i = 0 , the mechanical power input mustbalance the electrical power output plus all the losses due to damping and friction. DividingEq. (S4) by a power base ( P base ) allows us to express power in per-unit quantities: ω sm J i P base ¨ (cid:101) φ m ,i = − ω sm d i P base ˙ (cid:101) φ m ,i + P (pu)m ,i − ω d i P base − P (pu)e ,i , (S5)23hich leads to H i ω sm ¨ (cid:101) φ m ,i = − D i ω sm ˙ (cid:101) φ m ,i + (cid:101) P (pu)m ,i − P (pu)e ,i , (S6)where H i ≡ J i ω /P base is the inertia constant (which equals the kinetic energy of the rotorat the synchronous speed), D i ≡ d i ω /P base is the damping coefficient, and (cid:101) P (pu)m ,i ≡ P (pu)m ,i − ω d i /P base is the net power input.We assume that each generator i can be represented by two nodes: an internal node, withvoltage phasor E i = | E i | exp( jδ i ) , whose magnitude | E i | is assumed to be constant; and theterminal node, whose voltage is V i = | V i | exp( jθ i ) and matches the generator’s terminal voltage.Note that j denotes the imaginary unit. The internal and terminal nodes are connected throughan (internal) impedance, which we measured in our experiment (see Sec. S2.3 below). If thegenerator has p pole pairs, the mechanical and electrical angles are related by δ i = p (cid:101) φ m ,i (andthus ˙ δ i = p ˙ (cid:101) φ m ,i , ¨ δ i = p ¨ (cid:101) φ m ,i , and the synchronous voltage frequency ω s = pω sm ), leading to asecond form of the equation of motion: H i ω s ¨ δ i = − D i ω s ˙ δ i + (cid:101) P (pu)m ,i − P (pu)e ,i . (S7)This is equivalent to the classical model. The voltage at the terminal node is related to theterminal voltages of the other generators by the Kirchoff laws as P i = n (cid:88) k =1 | V i V k Y ik | sin( θ i − θ k − α ik + π/ , (S8)where P i is the real power injection from generator i into the network, and Y ik = | Y ik | exp( jα ik ) are the complex entries of the admittance matrix Y = ( Y ik ) representing the AC electric net-work. Through a procedure known as Kron reduction applied to the network (see, for example,Refs. 44 and 45 of the main text), we can also express the real power injection P (pu)e ,i at theinternal node of generator i in a form similar to Eq. (S8), but in terms of the internal voltageangles δ i : P (pu)e ,i = n (cid:88) k =1 | E i E k Y ik | sin( δ i − δ k − α ik + π/ , (S9)where Y = ( Y ik ) is the effective admittance matrix with Y ik = | Y ik | exp( jα ik ) , representing the24ffective coupling between generators. Substituting this into Eq. (S7) and defining β i ≡ D i H i = d i J i ,a i ≡ ω s (cid:2) (cid:101) P (pu)m ,i − | E i Y ii | cos( α ii ) (cid:3) H i ,c ik ≡ ω s | E i E k Y ik | H i ,γ ik ≡ α ik − π/ , (S10)we obtain the oscillator network model in the form of Eq. (1) of the main text (with ε = 0 ).Note that β i is constant because d i and J i are both physical constants characterizing thecorresponding generator. Note also that, for a synchronous state with a given (constant) fre-quency ω s , the parameter H i = J i ω /P base = J i ω / ( p P base ) and the admittances Y ik = | Y ik | exp( jγ ik ) are constant. Thus, if the parameters (cid:101) P (pu)m ,i and E i are constant, then all the pa-rameters defined in Eq. (S10) would be constant. In our experiment, we validated the constancyof ω s , (cid:101) P (pu)m ,i , and E i directly from the measurements for each time-series segment (Supplemen-tary Fig. 1). The linearization of Eq. (1) involving the matrix P in Eq. (2) assumes that allparameters in Eq. (S10) are constant, which is valid when ω s , (cid:101) P (pu)m ,i , and E i are constant. S2 Measurement of generator parameters
We experimentally measured all generator parameters required for our analysis, i.e., the parame-ter β i , the moment of inertia J i , and the internal impedance z int ,i for each generator ( i = 1 , , ).To simplify the notations, the generator index i in the subscript are omitted in the followingsubsections. S2.1 Effective damping parameter β The adjustable parameter β for each generator combines the moment of inertia and all damp-ing effects for the given generator, β = d/J . When both the mechanical and the electricaltorques are turned off, Eq. (S1) becomes ¨ φ m = − dJ ˙ φ m = − β ˙ φ m , which can be expressed as ˙ ω m = − βω m . Thus, the rotor decelerates at an exponential rate of β . We can directly measurethis rate by fitting an exponential decay curve to the experimentally measured rotor speed, afterthe generator is disconnected from the circuit and the DC motor driving it is turned off. We thus25quipped the rotating shaft of each generator with a reflective infrared phototransistor that de-tects markings on the shaft, providing direct measurements of the rotor frequency. The functionwe used to fit the frequency measurement ω m ( t ) is the following: y ( t ) = (cid:40) A, t < t ,A exp( − βt ) , t ≤ t < t , (S11)where the steady-state speed A , the turn-off time t , and the rate of deceleration β are fittingparameters. The turn-off time t needs to be fitted because the switching is initiated manuallyrather than by the microcontroller. The measured values of ω m ( t ) were taken up to t = t , atwhich ω m ( t ) ≈ Hz; below this threshold the exponential decay is not a good approximation.A typical fitting scenario is shown in Supplementary Fig. 4a, and a single set of measure-ments of β in the β A and β B configurations are shown in panels b–d of the same figure. Thestandard deviation of these measurements was ≈ . , indicating the precision of the β valuesthat we could configure and confirm by experimental measurements. At the beginning of eachexperimental run, we adjusted the breaks if necessary to ensure that the average of β (estimatedfrom at least ten measurements) was within ± . of the corresponding designed value for thegiven β configuration. S2.2 Moment of inertia J The measurement of J is based on the measurement of β , which was performed by estimatingthe exponential rate of deceleration as the DC motor driving the generator was turned off, whileno load was attached to the generator. If we repeat this experiment but with a load connectedto the terminal of the generator, then the generator would decelerate slightly faster. Using aresistor with R = 2 Ω as the load (rated for a maximum of W dissipation) would result in atime-dependent electric torque, T e ( t ) = P e ( t ) ω m ( t ) = V ( t ) Rω m ( t ) , (S12)where V ( t ) is the instantaneous r.m.s. voltage measured across the generator terminals. Substi-tuting into Eq. (S1) and using ω m = ˙ φ m and T m = 0 , we obtain ˙ ω m = − βω m − V ( t ) J Rω m ( t ) , (S13)26hose solution is given by ω m ( t ) = exp ( − βt ) (cid:18) ω (0) − J (cid:90) t V ( t (cid:48) ) R exp(2 βt (cid:48) ) dt (cid:48) (cid:19) . (S14)With this and the measured time series of V ( t ) , we can calculate ω m ( t ) numerically for any t .The parameter J can then be estimated by minimizing the sum of squared differences betweenthese calculated values of ω m ( t ) and the corresponding direct measurements for all t (taken upto a time at which ω m ( t ) ≈ Hz).The necessary voltage values V ( t ) for Eq. (S14) were provided by direct measurement ofthe terminal voltages, which were taken simultaneously with the rotor speed measurements(using a reflective infrared phototransistor, as was done for the β measurements; see Sec. S2.1above). Since the voltage measurements were taken at a much higher rate than the rotor speedmeasurements ( , samples per second for the voltage vs. samples per second for therotor speed at ω m = 40 × π rad/s, which decreases as the rotor slows down), we used thesampling rate of the voltages for the time discretization of the integral in Eq. (S14). For the β value in Eq. (S14), we substituted the estimate from Sec. S2.1 above, which is an averageover many measurements (see, for example, Supplementary Fig. 4). The estimate of J for eachgenerator was then obtained by repeating this procedure several times (using the same average β value) and taking the average. A typical fitting curve, as well as the estimated values of J areprovided in Supplementary Fig. 5. S2.3 Internal impedance z int The measurement of the generator’s internal impedance z int = r int + jx int requires the open-circuit voltage ( V oc ) and the short-circuit current ( I sc ) to be measured at the same rotor speed.We can take these readings directly using the voltage and current sensors. Our experimentalsetup has a set of ACS712 current sensors that can measure the current passing through eachgenerator’s terminal. In large industrial synchronous generators the rotating magnetic field isinduced by an adjustable field current, thus making the V oc and I sc values dependent on thiscurrent as well. Our permanent-magnet generators, however, have constant magnetic field, andthus V oc and I sc are functions of the rotor speed only. Both V oc and I sc generally increase linearlywith the speed until the air gap flux starts to saturate the iron parts of the generator. This linearregime was observed up to rotations per second for our generators. We performed all ourexperiments and measurements in this regime.27he resistive component r int of the internal impedance in the classical model is negligiblefor a typical synchronous generator, but we found that it was non-negligible for our generators.We measured both the resistive and reactive components using the following procedure. Wefirst obtained the magnitude of the internal impedance using the relation | z int | = V oc I sc , with V oc and I sc measured at rotations per second. We then measured the real part r int directlyby connecting a multimeter to the generator terminals in resistance measurement mode. Theremaining imaginary part of the internal impedance was calculated using x int = (cid:112) | z int | − r .The measured values of V oc , I sc , and r int , as well as the calculated x int , are as follows: • Generator 1 : V oc = 1 . V, I sc = 0 . A, r int = 1 . , x int = 3 .
34 Ω . • Generator 2 : V oc = 1 . V, I sc = 0 . A, r int = 1 . , x int = 3 .
36 Ω . • Generator 3 : V oc = 1 . V, I sc = 0 . A, r int = 1 . , x int = 3 .
20 Ω .The mean and the standard deviation among the three generators are V oc = 1 . ± . V, I sc = 0 . ± . A, r int = 1 . ± . , and x int = 3 . ± .
087 Ω . Note that the measuredinternal impedances of our generators would be called the synchronous impedance in the powersystems literature, since I sc is measured at a steady speed. This quantity is usually used tomodel the steady-state synchronous dynamics after a short period of transients following ashort circuit, while the so-called transient impedance is used to model the transients. For ourgenerators, however, we found no detectable transient, implying that the transient impedance isessentially equal to the synchronous impedance. Therefore, we used the z int measured by theabove procedure in our analysis of short-term dynamics.28 upplementary Figures Supplementary Fig. 1:
Validation of model assumptions. a – c , Histograms of normalized fluctuationsfor the parameters ω s ( a ), (cid:101) P (pu)m ,i ( b ), and E i ( c ) (defined in Sec. S1). We quantify experimental fluctua-tions of a given parameter q of a given generator in a given time-series segment by σ q / (cid:104) q (cid:105) , where σ q and (cid:104) q (cid:105) are the standard deviation and average, respectively, of the measured instantaneous values of q acrossthe segment. The histograms, taken over all generators and over all time-series segments for each param-eter, indicate that the (relative) magnitude of fluctuations in these parameters is approximately constantin each time-series segment, validating the corresponding assumptions underlying Eq. (1) of the maintext. Robustness of measured stability improvement. a , Distribution of σ λ max , thestandard deviation of λ max as it varies along the trajectory in a given segment, where λ max is computedusing the instantaneous δ ∗ i values at each data point in the segment. b , Histogram of predicted stabilityimprovement λ max ( β A ) − λ max ( β B ) when the dynamical parameters of the generators are randomlyperturbed (using , realizations). The specific parameters perturbed were V oc , I sc , and J for eachgenerator (defined in Secs. S1, S2.2, and S2.3). The perturbations were drawn from the normal distri-bution having the mean and the standard deviation of the measured values of a given parameter of eachgenerator (see Sec. S2.3 and Supplementary Fig. 5). The vast majority of the realizations lead to a mea-surable stability difference between the two β configurations, showing that converse symmetry breakingis observable and robust even under realistic uncertainties in the state and parameters of the system. Segment length (s)
Supplementary Fig. 3:
Length distribution of time-series segments.
The histogram shows the dis-tribution of the lengths of the time-series segments used in the analysis shown in Fig. 3 (for both the β A and β B configurations). bdc Supplementary Fig. 4:
Measurement of parameter β . a , Typical fitting scenario, where the steady-state frequency, the turn-off time (shifted to t = 0 ), and the rate of deceleration β are fitted to thedata from rotor frequency measurement. b , c , Measurements of the β value for generators 1 and 2,respectively. d , Measurements of the β values for generator 3 in the uniform ( β A , red dots) and non-uniform ( β B , blue squares) configurations. In b – d , the mean and standard deviation of the data pointsare indicated by the solid and dashed lines, respectively. bc d decay w/o load Supplementary Fig. 5:
Measurement of parameter J . a , Decay of the measured rotor frequency froma typical run with a resistive load of , along with the corresponding fit of the steady-state frequency,turn-off time (shown as t = 0 here), and parameter J . We also show as a reference the predicted decaywithout the load. b – d , Measurements of J for generators 1, 2, and 3, respectively. In b – d , the mean andstandard deviation of the data points are indicated by the solid and dashed lines, respectively., the mean andstandard deviation of the data points are indicated by the solid and dashed lines, respectively.