Designing temporal networks that synchronize under resource constraints
DDesigning temporal networks that synchronize under resource constraints
Yuanzhao Zhang and Steven H. Strogatz Center for Applied Mathematics, Cornell University, Ithaca, New York 14853, USA
Being fundamentally a non-equilibrium process, synchronization comes with unavoidable energycosts and has to be maintained under the constraint of limited resources. Such resource constraintsare often reflected as a finite coupling budget available in a network to facilitate interaction andcommunication. Here, we show that introducing temporal variation in the network structure canlead to efficient synchronization even when stable synchrony is impossible in any static networkunder the given budget, thereby demonstrating a fundamental advantage of temporal networks.The temporal networks generated by our open-loop design are versatile in the sense of promotingsynchronization for systems with vastly different dynamics, including periodic and chaotic dynamicsin both discrete-time and continuous-time models. Furthermore, we link the dynamic stabilizationeffect of the changing topology to the curvature of the master stability function, which provides ana-lytical insights into synchronization on temporal networks in general. In particular, our results shedlight on the effect of network switching rate and explain why certain temporal networks synchronizeonly for intermediate switching rate.
I. INTRODUCTION
Synchronization is critical to the function of many in-terconnected systems [1], from physical [2] to technologi-cal [3] and biological [4]. Many such systems need to syn-chronize under the constraint of limited resources. Forinstance, energy dissipation is required to couple molec-ular biochemical oscillators through oscillator-oscillatorexchange reactions, which are responsible for synchro-nization in systems such as the cyanobacterial circadianclock [5]. For multiagent networks with distributed con-trol protocols, including robotic swarms, the synchro-nization performance is limited by the available budgetof control energy [6].Similarly, for networks of coupled oscillators, one im-portant resource is the total coupling budget [7], whichdetermines how strongly the oscillators can influence eachother. For a typical oscillator network, a minimum cou-pling strength σ c is needed to overcome transversal insta-bility and maintain synchronization. The network struc-tures that achieve synchronization with the minimumcoupling strength are optimal , and they are characterizedby a complete degenerate spectrum [8]—all eigenvaluesof the Laplacian matrix are identical, except the trivialzero eigenvalue associated with perturbations along thesynchronization trajectory. Below σ c , there is no networkstructure that can maintain synchrony without violatingthe resource constraint.The results above, however, are derived assuming thenetwork to be static. That is, the network connectionsdo not change over time. Previous studies have shownthat temporal networks [9–15] can synchronize betterthan two of their static counterparts—namely, those ob-tained either by freezing the network at given time in-stants [16–19] or by averaging the network structure overtime [20–22]. But it remains unclear whether there aretemporal networks that can outperform all possible staticnetworks. In particular, can temporal variations synchro-nize systems beyond the fundamental limit set by the optimal static networks? This question is especially in- teresting given that past studies have often focused onthe fast-switching limit, for which the network structurechanges much faster than the node dynamics. These fast-switching networks are equivalent to their static, time-averaged counterparts in terms of synchronization stabil-ity [17, 23, 24]. Thus, no temporal networks can outper-form optimal static networks in the fast-switching limit.In this Article, we show that the full potential of tem-poral networks lies beyond the fast-switching limit, amessage echoed by several recent studies [21, 25, 26].Importantly, by allowing a network to vary in time ata suitable rate, synchronization can be maintained evenwhen the coupling strength is below σ c for all time t . Wealso develop a general theory to characterize the synchro-nizability of commutative temporal networks. The useof commutative graphs in synchronization was pioneeredin Refs. [18, 20] and subsequently adopted in numerousstudies [22, 25, 27] for its potential of generating analyti-cal insights. A new insight provided by our theory is thatthe effectiveness of introducing time-varying coupling de-pends critically on the curvature of the master stabilityfunction [28] at its first zero. Moreover, we demonstrateanalytically that this condition for improved synchroniz-ability in temporal networks is universally satisfied bycoupled one-dimensional maps. II. NETWORKS OF COUPLED OSCILLATORS
We start by considering systems described by the fol-lowing dynamical equations:˙ x i = F ( x i ) − σ n (cid:88) j =1 L ij ( t ) H ( x j ) , i = 1 , . . . , n, (1)where L = ( L ij ) is the normalized Laplacian matrix rep-resenting a diffusively coupled network. Here, L ij = δ ij (cid:80) k A ik − A ij , with δ ij being the Kronecker delta and A ij encoding the edge weight from node j to node i . Anoverall normalization factor is chosen so that the sum of a r X i v : . [ n li n . AO ] M a r all entries in A , (cid:80) ≤ i,j ≤ n A ij , equals n −
1. As a con-sequence, n − (cid:80) ni =1 L ii ( t ) = n − (cid:80) ni =2 λ i ( t ) = 1, wherethe sum over the eigenvalues λ i ( t ) starts from i = 2 be-cause the trivial eigenvalue λ associated with the eigen-vector v = (1 , , . . . , (cid:124) is always 0. As a result of thenormalization, the amount of resources (per node) usedto maintain synchronization can be quantified solely bythe coupling strength σ for networks of different sizesand densities. The d -dimensional vector x i describes thestate of node i , F is the vector field dictating the in-trinsic node dynamics, and H is the coupling functionmediating interactions between different nodes.To determine the stability of the synchronization state x ( t ) = x ( t ) = · · · = x n ( t ) = s ( t ), we study the varia-tional equation˙ δ = [ n ⊗ J F ( s ) − σ L ( t ) ⊗ J H ( s )] δ . (2)Here, δ = ( x − s , . . . , x n − s ) (cid:124) is the perturbation vector, n is the n × n identity matrix, ⊗ represents the Kro-necker product, and J is the Jacobian operator. Whenthe Laplacian matrices L ( t ) and L ( t (cid:48) ) commute for any t and t (cid:48) , following the master stability function formal-ism [18, 28], we can find an orthogonal matrix Q suchthat Q (cid:124) L ( t ) Q is diagonal for all time t , thus decouplingEq. (2) into n independent d -dimensional equations˙ η i = [ J F ( s ) − σλ i ( t ) J H ( s )] η i , i = 1 , . . . , n. (3)Here, { η i } is linked to the original coordinates throughthe relation ( η , . . . , η n ) (cid:124) = Q (cid:124) δ . Each decoupled equa-tion describes the evolution of an independent perturba-tion mode η i . In order for synchronization to be stable,all perturbation modes transverse to the synchronizationmanifold (namely, the modes η to η n ) must asymptoti-cally decay to zero. Since the decoupled variational equa-tions are all of the same form and only differ in λ i ( t ), it isinformative to study the maximum Lyapunov exponentof the equation ˙ ξ = [ J F ( s ) − αJ H ( s )] ξ (4)as a function of α . We refer to this function as the masterstability function and denote it as Λ( α ).As we will show throughout the rest of the paper, ifΛ (cid:48)(cid:48) ( α ) < α ) first becomes negative at α = σ c (Fig. 1), then it is guaranteed that there exist temporalnetworks that outperform optimal static networks. Intu-itively, this is because introducing temporal variation inthe network structure allows all nonzero λ i ( t ) to spenda significant amount of time above 1, the optimal valueachievable by static networks. (For static networks, be-cause (cid:80) ni =2 λ i = n −
1, there must exist 0 < λ i < λ i = 1for all i ≥ (cid:48)(cid:48) ( α ) < λ i ( t ) > λ i ( t ) < FIG. 1. Example master stability function for which temporalnetworks can synchronize stably below the critical couplingstrength σ c . III. TEMPORAL NETWORKS OUTPERFORMOPTIMAL STATIC NETWORKS
In order to illustrate a simple scheme for designingtemporal networks that synchronize for coupling strengthbelow the critical value σ c , we construct a class of Lapla-cian matrices that have the following spectrum (Fig. 2a): λ i ( t ) = i = 1 , n − m A sin( ωt ) i = 2 , . . . , m + 1 , − n − n − m − A sin( ωt ) i = m + 2 , . . . , n. (5)The nonzero eigenvalues split into two groups with atime-varying gap between them, while their sum remainsequal to n − t . Intuitively, some of theperturbation modes borrow resources from the others toremain stable and then return the favor at a later time.As a result, this kind of dynamic stabilization achievesglobal synchronization with very limited resources.One can design networks with a given spectrum byspecifying a set of orthonormal eigenvectors { v i } [18].For our purpose, any choice of { v i } containing v =(1 , , . . . , (cid:124) is valid, which gives rise to a whole rangeof synchronization-boosting temporal networks. Here,for concreteness, we adopt the eigenbasis proposed inRef. [29]: v i = ( 1 (cid:112) i ( i − , · · · , (cid:112) i ( i − (cid:124) (cid:123)(cid:122) (cid:125) i − , − i − (cid:112) i ( i − , , · · · , (cid:124) (cid:123)(cid:122) (cid:125) n − i copies ) , (6)where i ≥
2. Combining Eqs. (5) and (6) using the for-mula L ( t ) = (cid:80) ni =2 λ i ( t ) v i v (cid:124) i gives rise to a temporalnetwork described by the following weighted adjacencymatrix (Fig. 2b): A ij ( t ) = (cid:40) λ n ( t ) n + λ ( t ) − λ n ( t ) m +1 i, j ≤ m + 1 , i (cid:54) = j, λ n ( t ) n i or j > m + 1 , i (cid:54) = j. (7)Substituting Eq. (5) into Eq. (7) shows that edges con-necting the first m + 1 nodes have a time-dependent FIG. 2.
Designing temporal networks that synchronize better than optimal static networks. a
Evolution of thenonzero Laplacian eigenvalues described in Eq. (5), which are split into two degenerate groups. b Temporal network constructedfrom the Laplacian eigenvalues in a . The weight of each edge is represented by its thickness. In addition, edges whose weightis larger than n are colored orange, whereas those with weight less than n are colored cyan. For this network diagram, we set n = 11 and m = 5, and the corresponding weighted adjacency matrix is given by Eq. (9). Visually, we can see that differentparts of the network are being strengthened in an alternating fashion. weight of n + n ( n − − m ( m +1)( n − nm ( m +1)( n − m − A sin( ωt ), while theweight of the other edges evolve according to n − ( n − n ( n − m − A sin( ωt ). The choice of the time-varying termsin( ωt ) is not essential; the sine function can be replacedby any other periodic function p ( t ) with period T thatsatisfies (cid:82) T p ( t ) d t = 0.When assuming n odd and m = n − , we get a particu-larly simple class of temporal networks whose transverseperturbation modes all have the same stability (analo-gous to the defining property of optimal static networks): λ i ( t ) = i = 1 , A sin( ωt ) i = 2 , . . . , n +12 , − A sin( ωt ) i = n +32 , . . . , n, (8) A ij ( t ) = (cid:40) − n +1 ) A sin( ωt ) n i, j ≤ n +12 , i (cid:54) = j, − A sin( ωt ) n i or j > n +12 , i (cid:54) = j. (9) IV. CRITICAL ROLE OF THE SWITCHINGRATE
To demonstrate the effectiveness of our design, weequip the temporal networks described by Eq. (9) withconcrete node dynamics and probe their synchronizabil-ity in depth. Here, we choose Stuart-Landau oscillatorsas our first example, since they represent the canonicaldynamics of systems in the vicinity of a Hopf bifurcation[30]. The oscillators evolve according to the followingdynamical equation:˙ Z j = Z j − (1+i c ) | Z j | Z j − σ n (cid:88) k =1 L jk ( t )(1+i c ) Z k , (10)where Z j = x j + i y j = r j e i θ j ∈ C represents the stateof the j th oscillator. Equation (10) is the discrete-spacecounterpart of the Ginzburg-Landau equation [31] andadmits a limit-cycle synchronous state Z j ( t ) = e − i c t ∀ j . By writing the perturbations in polar coordinates, wefind that the Jacobian terms in Eq. (4) become J F = (cid:0) − − c (cid:1) and J H = (cid:0) − c c (cid:1) , both of which are con-stant matrices. Thus, according to Eq. (4), the masterstability function can be obtained by solving a charac-teristic polynomial equation and has the following form[27]: Λ( α ) = − α − (cid:113) − c c α − c α . (11)Figure 3a shows Λ( α ) for c = − . c = 4, whichclearly has Λ (cid:48)(cid:48) ( α ) < α ≈ η i = B i ( t ) η i , (12)where B i ( t ) = (cid:16) − − σλ i ( t ) c σλ i ( t ) − c − c σλ i ( t ) − σλ i ( t ) (cid:17) is periodic withperiod T = πω (henceforth we drop the subscript i toease the notation). According to Floquet theory [32], thesolution to Eq. (12) must be of the form e µt P ( t ), where P ( t ) has period T . The Floquet exponents µ and µ canbe extracted by finding the principal fundamental matrix,and their real parts are the corresponding Lyapunov ex-ponents [33]. Figure 3b shows the maximum Lyapunovexponent Γ = max { Re( µ ) , Re( µ ) } as a function of ω fordifferent values of the temporal activity A . (It is clearfrom Eq. (8) that all transverse perturbation modes havethe same Γ. Thus, Γ is also the maximum transverse Lyapunov exponent and determines the synchronizationstability.) We set the coupling strength to slightly below σ c at σ = 2 . A is increased, Γ becomesnegative for an increasingly wide range of switching rate ω , signaling that the temporal variation in the networkstructure is successfully stabilizing synchronization underthe given coupling budget.Since the only difference between Eqs. (3) and (4) isthe periodic λ ( t ) vs. the fixed α , it is natural to expect FIG. 3.
Temporal networks enable synchronizationamong Stuart-Landau oscillators. a
Master stabilityfunction for Stuart-Landau oscillators. Parameters are setto c = − . c = 4. b Maximum Lyapunov exponentΓ as a function of the switching rate ω for different valuesof the temporal activity A (solid lines). The dashed linesindicate the slow-switching limit predicted by the averagedmaster stability function Λ. the stability of the temporal network to be related to themaster stability function averaged over a suitable rangeof α . Specifically, one might reasonably associate Γ withthe averaged master stability functionΛ = (cid:90) λ max λ min W ( λ )Λ( σλ ) d λ, (13)where W ( λ ) is the probability distribution of λ (it fol-lows that (cid:82) λ max λ min W ( λ ) d λ = 1). However, it is clear thatΛ cannot be used to predict Γ in general. One immedi-ate observation is that Λ does not depend on the rate inwhich λ ( t ) is changing (it only depends on the distribu-tion of λ ), whereas the curves representing Γ in Fig. 3bclearly depend on the switching rate ω . Indeed, in orderto go from Γ to Λ, we are required to shuffle B ( t ) tempo-rally in Eq. (12). This operation is forbidden when thematrices { B ( t ) | t ∈ R } do not commute (or, equivalently,when { B ( t ) | t ∈ R } cannot be simultaneously diagonal-ized). To see why, we can look at the formal solution toEq. (12) expressed in terms of the matrix exponential: η ( t ) = η (0)e Ω ( t ) , (14)where Ω ( t ) is given by the Magnus expansion [34]: Ω ( t ) = (cid:90) t B ( τ ) d τ + 12 (cid:90) t d τ (cid:90) τ d τ (cid:48) [ B ( τ ) , B ( τ (cid:48) )]+ higher-order terms involving nested matrix commutators . (15)Here, [ B ( τ ) , B ( τ (cid:48) )] = B ( τ ) B ( τ (cid:48) ) − B ( τ (cid:48) ) B ( τ ) is thematrix commutator. Equation (15) makes it clear that { B ( τ ) | < τ < t } can be shuffled without affecting Ω ( t )if and only if [ B ( τ ) , B ( τ (cid:48) )] = 0 for all τ (cid:48) < τ < t , inwhich case everything on the right-hand side except thefirst term vanishes.However, Λ is still extremely informative on whethera given temporal network can synchronize or not. Inparticular, for ω → { B ( t ) | t ∈ R } commute, η can be decomposed into components thatgrow along the eigendirections of B ( t ), whose growth rateis dictated by the corresponding eigenvalues. Eventually,the component along the direction with the largest eigen-value becomes dominant. However, when { B ( t ) | t ∈ R } do not commute, the growth along the eigendirections areoften “interrupted”, since the eigenvectors of B ( t ) are nolonger fixed and will rotate over time. To keep track ofthe growth of the dominant component, we must project η onto the new dominant eigendirection upon rotation.These frequent projections can significantly influence theasymptotic growth rate (this is also why the maximumLyapunov exponent is usually not the mean of the max-imum local Lyapunov exponents). At the slow-switchinglimit, η can grow along an eigendirection uninterruptedfor long enough that the effect of the projections becomesnegligible. In this case, Γ is determined by the averagegrowth rate of η in the dominant direction of each B ( t ),which is exactly Λ.It is worth noting that the equivalence between Γ andΛ at the slow-switching limit is not specific to Stuart-Landau oscillators and can be expected for generic os-cillator models. As a result, Λ (cid:48)(cid:48) ( α ) < λ around 1 to be small and Tay-lor expand Λ( α ) around α . Then the averaged masterstability function for coupling strength σ = σ c isΛ = (cid:90) (cid:15) − (cid:15) W ( λ )Λ( σ c λ ) d λ = (cid:90) (cid:15) − (cid:15) W ( λ ) (cid:104) Λ( σ c ) + Λ (cid:48) ( σ c )( λ − (cid:48)(cid:48) ( σ c )( λ − (cid:105) d λ + O ( (cid:15) )= Λ( σ c ) (cid:124) (cid:123)(cid:122) (cid:125) =0 +Λ (cid:48) ( σ c ) (cid:90) (cid:15) − (cid:15) W ( λ )( λ −
1) d λ (cid:124) (cid:123)(cid:122) (cid:125) =0 + 12 Λ (cid:48)(cid:48) ( σ c ) (cid:90) (cid:15) − (cid:15) W ( λ )( λ − d λ (cid:124) (cid:123)(cid:122) (cid:125) > + O ( (cid:15) ) . (16)Thus, if Λ (cid:48)(cid:48) ( α ) = Λ (cid:48)(cid:48) ( σ c ) <
0, then Λ < σ = σ c and stability is guaranteed to be improved at the slow-switching limit, where Γ = Λ. This improvement is ex-pected to extend into intermediate switching rate due tothe continuity of Γ as a function of ω .At the other limit, for ω → ∞ (i.e., fast-switching net-works), Γ clearly does not match with Λ. In particular,Γ does not depend on A . For the system in Fig. 3b, Γapproaches Λ( σ ) as ω → ∞ , which is the value expectedfor an optimal static network at coupling strength σ (inthis case the time-averaged network is a complete graph FIG. 4.
Temporal networks synchronize only for in-termediate switching rate.
Evolution of the oscillatorstates x i and the synchronization error ∆ for a : ω = 0, b : ω = 1, and c : ω = 100. The oscillator parameters are thesame as in Fig. 3 and the underlying temporal network isillustrated in Fig. 2b. with uniform edge weights). The mapping from a tem-poral network to its time-averaged counterpart at thefast-switching limit is intuitive and well established inthe literature [17, 23, 24].The results above provide new insights into the intrigu-ing phenomenon that certain temporal networks onlysynchronize for intermediate switching rate [21, 25, 26]:When switching is too fast, the temporal network reducesto its static counterpart and one cannot take full advan-tage of the temporal variation in the connections; whenswitching is too slow, although the asymptotic stabilitymight be maximized, the system would have desynchro-nized long before the network experiences any meaningfulchange. Thus, the sweet spot often emerges at an inter-mediate switching rate.In Fig. 4, we show typical trajectories of n = 11Stuart-Landau oscillators on the temporal networks de-scribed by Eq. (9), with the temporal activity set to A = 0 .
15. Systems in all three panels are initiatedclose to the synchronous state, and their only differ-ence lies in the switching rate ω , which allows us tocompare networks with static, moderate-switching, andfast-switching topologies. By monitoring the synchro-nization error ∆( t ), defined as the standard deviationamong Z j ( t ), we see that only the system with an inter-mediate switching rate ( ω = 1, panel b) can maintainstable synchrony. Interestingly, ∆( t ) in that system goesdown non-monotonically and is bounded from above byperiodic envelopes. The width of each envelope is 2 π ,which coincides with the period of the changing networktopology. V. UNIVERSAL STABILIZATION OFLOW-DIMENSIONAL MAPS
The framework developed so far can be readily trans-ferred from differential equations to discrete maps, fromcontinuous variation in network topology to discreteswitching, and from periodic oscillator dynamics tochaotic ones. The discrete-time analog of Eq. (1) canbe written as x i [ t + 1] = β F ( x i [ t ]) − σ n (cid:88) j =1 L ij [ t ] H ( x j [ t ]) . (17)To demonstrate the advantage of temporal networks inthese extended settings, we focus on the following classof coupled one-dimensional discrete maps: x i [ t + 1] = βF ( x i [ t ]) − σ n (cid:88) j =1 L ij [ t ] F ( x j [ t ]) , (18)where F : R → R is the mapping function. As we showbelow, this setup allows us to develop an elegant theorythat offers new insights.Similar to the continuous-time case, the synchroniza-tion stability is determined by the decoupled variationalequations η i [ t + 1] = (cid:104) ( β − σλ i [ t ]) F (cid:48) ( s [ t ]) (cid:105) η i [ t ] . (19)For fixed λ , the Lyapunov exponent of Eq. (19)is given by ln | β − σλ | + Γ s , where Γ s =lim T →∞ T (cid:80) T t =1 ln | F (cid:48) ( s [ t ]) | is a finite constant. Thus,the master stability function has the universal form(illustrated in Fig. 1)Λ( α ) = ln | α − β | + Γ s . (20)Taking the second derivative with respect to α , we seethat Λ (cid:48)(cid:48) ( α ) = − α − β ) < . (21)Thus, synchronization in any system described byEq. (18) can benefit from the temporal networks designedin this paper. In particular, this holds for any mappingfunction F , which encompasses important dynamical sys-tems such as logistic maps, circle maps, and Bernoullimaps.For concreteness, we set F ( x ) = sin ( x + π/
4) and β =2 . s = − . A ij [ t ] = (cid:40) − (cid:98) t/T (cid:99) (6 − n +1 ) An i, j ≤ n +12 , i (cid:54) = j, − ( − (cid:98) t/T (cid:99) An i or j > n +12 , i (cid:54) = j, (22)where (cid:98)·(cid:99) is the floor function. Basically, the networkswitches between two configurations every T iterations,with each configuration being the extremal in the contin-uous scheme described by Eq. (9). Consequently, everynonzero eigenvalue of the temporal network alternatesbetween 1 + 2 A and 1 − A with period T .Again, the averaged master stability function Λ ac-curately predicts the stability of the temporal networkat the slow-switching limit. More interestingly, for sys-tems described by Eq. (18), the connection is muchstronger: Λ determines the stability of the temporal net-work for all switching periods T . To see why, we notethat the synchronization stability is determined by thelimit product (cid:81) ∞ t =1 ( β − σλ [ t ]) F (cid:48) ( s [ t ]). Normally, theseare matrix products and cannot be reordered. However,since 1 × T →∞ T ln (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T (cid:89) t =1 ( β − σλ [ t ]) F (cid:48) ( s [ t ]) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:90) λ max λ min W ( λ ) ln | β − σλ | d λ + lim T →∞ T T (cid:88) t =1 ln | F (cid:48) ( s [ t ]) | = (cid:90) λ max λ min W ( λ ) (ln | β − σλ | + Γ s ) d λ = (cid:90) λ max λ min W ( λ )Λ( σλ ) d λ = Λ . (23)This independence of Γ on T might seem contradictoryto the fact that, at the fast-switching limit, temporalnetworks can be reduced to their static counterparts. Butnotice that there is usually no fast switching in discrete-time systems—even if the network topology changes atevery iteration, it is still evolving at the same timescaleas the node dynamics. Moreover, unlike in continuous-time systems [16, 17, 23, 24], the discrete nature of thedynamics precludes the use of the averaging techniques[16, 24] essential for connecting fast-switching networkswith their time-averaged counterparts. Thus, one cannotmap a temporal network to its time-averaged counterpartin discrete-time systems even when the network topologychanges much more rapidly than the node dynamics.In Fig. 5, we show the maximum transverse Lyapunovexponent Γ of the synchronization state in the optoelec-tronic system for σ = 1, which is slightly below σ c .The dashed line corresponds to the theoretical predic-tion of Γ based on the averaged master stability func-tion Λ = (ln | A − β | + ln | − A − β | ) + Γ s . Asexpected, the static network ( A = 0), despite being opti-mal, is unstable. As the temporal activity A is increased, FIG. 5.
Temporal networks promote synchroniza-tion universally in discrete-time systems with low-dimensional node dynamics.
The maximum transverseLyapunov exponent Γ decreases as the temporal activity A isincreased. The dashed line represents the theoretical predic-tion based on Λ, whereas the solid lines (shifted vertically forvisibility) are calculated directly from Eq. (19) for differentswitching periods T . Without the shift, all curves completelyoverlap (inset), which confirms our prediction that the stabi-lization provided by temporal networks does not depend onthe switching rate for coupled one-dimensional maps. Λ deceases and synchronization is eventually stabilized.On the other hand, the solid lines represent Γ obtainednumerically by evolving Eq. (19) for different switchingperiods T . These lines are shifted vertically by differentamounts in Fig. 5, purely as an aid to the eye. The un-shifted versions are shown in the inset. Notice that allthe lines collapse onto a single curve, demonstrating theexcellent agreement between theory and simulations.An interesting question is what happens when we in-troduce random fluctuations to the network structure ateach time step t , which makes the temporal network ape-riodic and the graph Laplacians noncommutative. InFig. 6, through direct simulations [33], we show that tem-poral networks still outperform optimal static networksin the presence of these random fluctuations. Here, weuse the same model of optoelectronic oscillators and thediscrete-switching network considered in Fig. 5, exceptthat independent random Gaussian perturbations of zeromean and standard deviation 0 . /n (10% of the averageedge weight) are added to the strength of each edge atevery time step. For temporal activity A = 0 (Fig. 6a),synchronization cannot be sustained at coupling strength σ = 1 .
05. For temporal activity A = 0 .
15 (Fig. 6b), syn-chronization is stabilized at the same coupling strengthby the variation in network structure. The network sizeis set to n = 99 and the switching period to T = 10 inour simulations, although the results do not depend onthese two parameters. FIG. 6.
Improved synchronization in aperiodic andnoncommutative temporal networks.
The temporal net-works are based on the discrete-switching networks given byEq. (22), which are further made aperiodic and noncommuta-tive by applying random Gaussian perturbations of zero meanto the strength of each link independently at every time step t . The standard deviation of the perturbations is fixed at0 . /n (10% of the average link weight). The spacetime plotsshow the evolution of the optoelectronic oscillators on tem-poral networks with a : temporal activity A = 0 and b : tem-poral activity A = 0 .
15. Both systems are initialized close tothe synchronous state. Synchronization persists only in thesecond system, even though the network in the first systemis an optimal static network (complete graph with uniformlink weights). Other parameters are set to n = 99, T = 10, β = 2 .
8, and σ = 1 . VI. DISCUSSION
To summarize, we have designed temporal networksthat synchronize more efficiently than optimal static net-works. These temporal networks are particularly rele-vant when the coupling budget available in a system tomaintain stable synchrony is limited. We provided ana-lytical insight into the synchronizability of commutativetemporal networks by linking it to the curvature of thecorresponding master stability function. In particular,our analysis reveals the subtle relation between the per-formance of a temporal network and its switching rate.The switching rate plays an especially critical role in sys-tems with high-dimensional oscillator dynamics, and net-works with intermediate switching rate often emerge asthe most effective. Our open-loop design has several advantages comparedto closed-loop schemes where the network structure is ad-justed on-the-fly based on feedbacks from the node states(often modeled by adaptive networks [36]). First, ourdesign does not depend sensitively on the node dynam-ics. As we have shown, the same design works for sys-tems with vastly different node dynamics, and it appliesreadily to both continuous-time and discrete-time sys-tems. Second, we do not need to monitor all the nodesconstantly, which also eliminates the possibility of beingdetrimentally influenced by measurement errors. Third,the evolution of the network is highly predictable andwe can easily control the coupling budget allocated tothe system at any given time t , a task that is far moredifficult in adaptive networks. On the other hand, closed-loop schemes have the advantage of being readily adap-tive to the changing environment and can react quickly tounexpected perturbations. A promising future directionwould be to devise hybrid schemes that combine the bestfrom both worlds, which could enable even more efficientand robust synchronization.In this work, for the sake of analytical tractability, wemostly focused on temporal networks whose Laplacianmatrices from different time instants commute. There isevidence that synchronization in temporal networks canbenefit when L ( t ) L ( t (cid:48) ) (cid:54) = L ( t (cid:48) ) L ( t ) [20]. It would there-fore be interesting to see whether our design of temporalnetworks could be further optimized by allowing non-commuting Laplacian matrices. In particular, can ran-dom fluctuations in the network structure (which give riseto noncommuting Laplacian matrices in general) outper-form our designed temporal networks? More generally,do optimal temporal networks exist for the purpose ofsynchronization, just like there are optimal static net-works? And if so, what are their defining characteristics?Finally, we hope our results can serve as an importantstep towards achieving efficient synchronization in com-plex interconnected systems. For example, many tem-poral networks arise naturally in the real world throughmoving agents, whose interactions depend on their spa-tial distance [37–40]. An exciting next step is to under-stand how our design can be implemented in such systemsand how the time-varying connections can be translatedinto the spatial movement of individual agents. ACKNOWLEDGMENTS
Y.Z. acknowledges support from the Schmidt ScienceFellowship. [1] S. H. Strogatz, Exploring complex networks, Nature ,268 (2001).[2] R. Roy and K. S. Thornburg Jr, Experimental synchro-nization of chaotic lasers, Phys. Rev. Lett. , 2009(1994). [3] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa,Spontaneous synchrony in power-grid networks, Nat.Phys. , 191 (2013).[4] R. E. Mirollo and S. H. Strogatz, Synchronization ofpulse-coupled biological oscillators, SIAM J. Appl. Math , 1645 (1990).[5] D. Zhang, Y. Cao, Q. Ouyang, and Y. Tu, The energycost and optimal design for synchronization of coupledmolecular oscillators, Nat. Phys. , 95 (2020).[6] J. Xi, C. Wang, H. Liu, and Z. Wang, Dynamic outputfeedback guaranteed-cost synchronization for multiagentnetworks with given cost budgets, IEEE Access , 28923(2018).[7] T. Nishikawa and A. E. Motter, Maximum performanceat minimum cost in network synchronization, Physica D , 77 (2006).[8] T. Nishikawa and A. E. Motter, Network synchroniza-tion landscape reveals compensatory structures, quanti-zation, and the positive effect of negative interactions,Proc. Natl. Acad. Sci. U.S.A. , 10342 (2010).[9] R. K. Pan and J. Saram¨aki, Path lengths, correlations,and centrality in temporal networks, Phys. Rev. E ,016105 (2011).[10] M. Starnini, A. Baronchelli, A. Barrat, and R. Pastor-Satorras, Random walks on temporal networks, Phys.Rev. E , 056115 (2012).[11] P. Holme and J. Saram¨aki, Temporal networks, Phys.Rep. , 97 (2012).[12] N. Masuda, K. Klemm, and V. M. Egu´ıluz, Temporalnetworks: Slowing down diffusion by long lasting inter-actions, Phys. Rev. Lett. , 188701 (2013).[13] E. Valdano, L. Ferreri, C. Poletto, and V. Colizza, Ana-lytical computation of the epidemic threshold on tempo-ral networks, Phys. Rev. X , 021005 (2015).[14] A. Paranjape, A. R. Benson, and J. Leskovec, Motifs intemporal networks, in Proceedings of the Tenth ACM In-ternational Conference on Web Search and Data Mining (2017) pp. 601–610.[15] A. Li, S. P. Cornelius, Y.-Y. Liu, L. Wang, and A.-L.Barab´asi, The fundamental advantages of temporal net-works, Science , 1042 (2017).[16] I. V. Belykh, V. N. Belykh, and M. Hasler, Blinkingmodel and synchronization in small-world networks witha time-varying coupling, Physica D , 188 (2004).[17] D. J. Stilwell, E. M. Bollt, and D. G. Roberson, Suf-ficient conditions for fast switching synchronization intime-varying network topologies, SIAM J. Appl. Dyn.Syst. , 140 (2006).[18] S. Boccaletti, D.-U. Hwang, M. Chavez, A. Amann,J. Kurths, and L. M. Pecora, Synchronization in dy-namical networks: Evolution along commutative graphs,Phys. Rev. E , 016102 (2006).[19] M. Porfiri, D. J. Stilwell, and E. M. Bollt, Synchroniza-tion in random weighted directed networks, IEEE Trans.Circuits Syst. I, Reg. Papers , 3170 (2008).[20] R. Amritkar and C.-K. Hu, Synchronized state of coupleddynamics on time-varying networks, Chaos , 015117(2006).[21] R. Jeter and I. Belykh, Synchronization in on-off stochas-tic networks: Windows of opportunity, IEEE Trans. Cir-cuits Syst. I, Reg. Papers , 1260 (2015). [22] S. Zhou, Y. Guo, M. Liu, Y.-C. Lai, and W. Lin, Randomtemporal connections promote network synchronization,Phys. Rev. E , 032302 (2019).[23] M. Porfiri, D. J. Stilwell, E. M. Bollt, and J. D. Skufca,Random talk: Random walk and synchronizability ina moving neighborhood network, Physica D , 102(2006).[24] J. Petit, B. Lauwens, D. Fanelli, and T. Carletti, Theoryof turing patterns on time varying networks, Phys. Rev.Lett. , 148301 (2017).[25] L. Chen, C. Qiu, and H. Huang, Synchronization withon-off coupling: Role of time scales in network dynamics,Phys. Rev. E , 045101 (2009).[26] O. Golovneva, R. Jeter, I. Belykh, and M. Porfiri, Win-dows of opportunity for synchronization in stochasticallycoupled maps, Physica D , 1 (2017).[27] C. Pereti and D. Fanelli, Stabilizing Stuart-Landau oscil-lators via time-varying networks, Chaos Solitons Fractals , 109587 (2020).[28] L. M. Pecora and T. L. Carroll, Master stability functionsfor synchronized coupled systems, Phys. Rev. Lett. ,2109 (1998).[29] A. Forrow, F. G. Woodhouse, and J. Dunkel, Functionalcontrol of network dynamics using designed Laplacianspectra, Phys. Rev. X , 041043 (2018).[30] Y. Kuramoto, Chemical Oscillations, Waves, and Turbu-lence , Vol. 19 (Springer Science & Business Media, Berlin,2012).[31] I. S. Aranson and L. Kramer, The world of the com-plex Ginzburg-Landau equation, Rev. Mod. Phys. , 99(2002).[32] P. A. Kuchment, Floquet theory for partial differentialequations , Vol. 60 (Birkh¨auser, Basel, 2012).[33] Code for performing network dynamics simulations andstability calculations are available at https://github.com/y-z-zhang/temporal_sync .[34] S. Blanes, F. Casas, J.-A. Oteo, and J. Ros, The Magnusexpansion and some of its applications, Phys. Rep. ,151 (2009).[35] J. D. Hart, D. C. Schmadel, T. E. Murphy, and R. Roy,Experiments with arbitrary networks in time-multiplexeddelay systems, Chaos , 121103 (2017).[36] T. Gross and B. Blasius, Adaptive coevolutionary net-works: A review, J. R. Soc. Interface , 259 (2008).[37] M. Frasca, A. Buscarino, A. Rizzo, L. Fortuna, andS. Boccaletti, Synchronization of moving chaotic agents,Phys. Rev. Lett. , 044102 (2008).[38] N. Fujiwara, J. Kurths, and A. D´ıaz-Guilera, Synchro-nization in networks of mobile oscillators, Phys. Rev. E , 025101 (2011).[39] K. P. O’Keeffe, H. Hong, and S. H. Strogatz, Oscillatorsthat sync and swarm, Nat. Commun. , 1504 (2017).[40] D. Levis, I. Pagonabarraga, and A. D´ıaz-Guilera, Syn-chronization in dynamical networks of locally coupledself-propelled oscillators, Phys. Rev. X7