Designing temporal networks that synchronize under resource constraints

Being fundamentally a non-equilibrium process, synchronization comes with unavoidable energy costs and has to be maintained under the constraint of limited resources. Such resource constraints are often reflected as a finite coupling budget available in a network to facilitate interaction and communication. Here, we show that introducing temporal variation in the network structure can lead to efficient synchronization even when stable synchrony is impossible in any static network under 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 promoting synchronization for systems with vastly different dynamics, including periodic and chaotic dynamics in both discrete-time and continuous-time models. Furthermore, we link the dynamic stabilization effect of the changing topology to the curvature of the master stability function, which provides analytical insights into synchronization on temporal networks in general. In particular, our results shed light on the effect of network switching rate and explain why certain temporal networks synchronize only for intermediate switching rate.

S ynchronization is critical to the function of many interconnected systems 1 , from physical 2 to technological 3 and biological 4 . Many such systems need to synchronize under the constraint of limited resources. For instance, energy dissipation is required to couple molecular biochemical oscillators through oscillator-oscillator exchange reactions, which are responsible for synchronization in systems such as the cyanobacterial circadian clock 5 . For multiagent networks with distributed control protocols, including robotic swarms, the synchronization performance is limited by the available budget of control energy 6 .
Similarly, for networks of coupled oscillators, one important resource is the total coupling budget 7 , which determines how strongly the oscillators can influence each other. For a typical oscillator network, a minimum coupling strength σ c is needed to overcome transversal instability and maintain synchronization. The network structures that achieve synchronization with the minimum coupling strength are optimal, and they are characterized by a complete degenerate spectrum 8 -all eigenvalues of the Laplacian matrix are identical, except the trivial zero eigenvalue associated with perturbations along the synchronization trajectory. Below σ c , there is no network structure that can maintain synchrony without violating the resource constraint.
The results above, however, are derived assuming the network to be static. That is, the network connections do not change over time. Previous studies have shown that temporal networks [9][10][11][12][13][14][15] can synchronize better than two of their static counterpartsnamely, those obtained either by freezing the network at given time instants [16][17][18][19] or by averaging the network structure over time [20][21][22] . But it remains unclear whether there are temporal networks that can outperform all possible static networks. In particular, can temporal variations synchronize systems beyond the fundamental limit set by the optimal static networks? This question is especially interesting given that past studies have often focused on the fast-switching limit, for which the network structure changes much faster than the node dynamics. These fast-switching networks are equivalent to their static, time-averaged counterparts in terms of synchronization stability 17,[23][24][25] . Thus, no temporal networks can outperform optimal static networks in the fast-switching limit.
In this article, we show that the full potential of temporal networks lies beyond the fast-switching limit, a message echoed by several recent studies 21,26,27 . Importantly, by allowing a network to vary in time at a suitable rate, synchronization can be maintained even when the coupling strength is below σ c for all time t. We also develop a general theory to characterize the synchronizability of commutative temporal networks. The use of commutative graphs in synchronization was pioneered in refs. 18,20 and subsequently adopted in numerous studies 22,26,28 for its potential of generating analytical insights beyond the fastswitching limit. An insight provided by our theory is that the effectiveness of introducing time-varying coupling depends critically on the curvature of the master stability function 29 at its first zero, which extends the results presented in ref. 30 . Moreover, we demonstrate analytically that the condition for improved synchronizability in temporal networks is universally satisfied by coupled one-dimensional maps.

Results
Networks of coupled oscillators. We start by considering systems described by the following dynamical equations: where L = (L ij ) is the normalized Laplacian matrix representing a diffusively coupled network. Here, L ij = δ ij ∑ k A ik − A ij , with δ ij being the Kronecker delta and A ij encoding the edge weight from node j to node i. An overall normalization factor is chosen so that the sum of all entries in A, ∑ 1≤i,j≤n A ij , equals n−1. As a consequence, 1 nÀ1 ∑ n i¼1 L ii ðtÞ ¼ 1 nÀ1 ∑ n i¼2 λ i ðtÞ ¼ 1, where the sum over the eigenvalues λ i (t) starts from i = 2 because the trivial eigenvalue λ 1 associated with the eigenvector v 1 = ð1; 1; ; 1Þ > is always 0. As a result of the normalization, the amount of resources (per node) used to maintain synchronization can be quantified solely by the coupling strength σ for networks of different sizes and densities. The d-dimensional vector x i describes the state of node i, F is the vector field dictating the intrinsic node dynamics, and H is the coupling function mediating interactions between different nodes.
To determine the stability of the synchronization state x 1 (t) = x 2 (t) = ⋯ = x n (t) = s(t), we study the variational equation Here, δ ¼ ðx 1 À s; ; x n À sÞ > is the perturbation vector, 1 n is the n × n identity matrix, ⊗ represents the Kronecker product, and J is the Jacobian operator. When the Laplacian matrices L(t) and Lðt 0 Þ commute for any t and t 0 , following the master stability function formalism 18,29 , we can find an orthogonal matrix Q such that Q ⊤ L(t)Q is diagonal for all time t, thus decoupling Eq. (2) into n independent d-dimensional equations Here, {η i } is linked to the original coordinates through the relation ðη 1 ; ; η n Þ > ¼ ðQ > 1 d Þδ. Each decoupled equation describes the evolution of an independent perturbation mode η i . In order for synchronization to be stable, all perturbation modes transverse to the synchronization manifold (namely, the modes η 2 to η n ) must asymptotically decay to zero. Since the decoupled variational equations are all of the same form and only differ in λ i (t), it is informative to study the maximum Lyapunov exponent of the equation as a function of α. We refer to this function as the master stability function and denote it as Λ(α). As we will show throughout the rest of the paper, if Λ″(α 0 ) < 0 when Λ(α) first becomes negative at α 0 = σ c (Fig. 1), then it is guaranteed that there exist temporal networks that outperform optimal static networks. Intuitively, this is because introducing temporal variation in the network structure allows all nonzero λ i (t) to spend a significant amount of time above 1, the optimal value achievable by static networks. (For static networks, because ∑ n i¼2 λ i ¼ n À 1, there must exist 0 < λ i < 1 unless all nonzero eigenvalues are identical, in which case λ i = 1 for all i ≥ 2 and the network is optimal.) If Λ″(α 0 ) < 0, the synchronization state can gain more stability while λ i (t) > 1 than the stability it loses during the period when λ i (t) < 1.
Temporal networks that outperform optimal static networks.
In order to illustrate a simple scheme for designing temporal networks that synchronize for coupling strength below the critical value σ c , we construct a class of Laplacian matrices that have the following spectrum ( Fig. 2a): The nonzero eigenvalues split into two groups with a timevarying gap between them, whereas their sum remains equal to n−1 for all time t. Intuitively, some of the perturbation modes borrow resources from the others to remain stable and then return the favor at a later time. As a result, this kind of dynamic stabilization achieves global synchronization with very limited resources.
One can design networks with a given spectrum by specifying a set of orthonormal eigenvectors {v i } 18 . For our purpose, any choice of {v i } containing v 1 = (1, 1, …, 1) ⊤ is valid, which gives rise to a whole range of synchronization-boosting temporal networks. Here, for concreteness, we adopt the eigenbasis proposed in ref. 31 : where i ≥ 2. Combining Eqs. (5) and (6) using the formula LðtÞ ¼ ∑ n i¼2 λ i ðtÞv i v > i gives rise to a temporal network described by the following weighted adjacency matrix (Fig. 2b): Substituting Eq. (5) into Eq. (7) shows that edges connecting the first m + 1 nodes have a time-dependent weight of 1 n þ nðnÀ1Þ 2 Àmðmþ1ÞðnÀ1Þ nmðmþ1ÞðnÀmÀ1Þ A sinðωtÞ, whereas the weight of the other edges evolves according to 1 n À ðnÀ1Þ nðnÀmÀ1Þ A sinðωtÞ. The choice of the time-varying term sinðωtÞ is not essential; the sine function can be replaced by any other periodic function p(t) with period T that satisfies R T 0 pðtÞ dt ¼ 0. When assuming n odd and m ¼ nÀ1 2 , we get a particularly simple class of temporal networks whose transverse perturbation modes all have the same stability (analogous to the defining property of optimal static networks): Critical role of the switching rate. To demonstrate the effectiveness of our design, we equip the temporal networks described by Eq. (9) with concrete node dynamics and probe their synchronizability in depth. Here, we choose Stuart-Landau oscillators as our first example, as they represent the canonical dynamics of systems in the vicinity of a Hopf bifurcation 32 . The oscillators evolve according to the following dynamical equation: where Z j ¼ x j þ iy j ¼ r j e iθ j 2 C represents the state of the jth oscillator. Equation (10) is the discrete-space counterpart of the Ginzburg-Landau equation 33 and admits a limit-cycle synchronous state Z j ðtÞ ¼ e Àic 2 t 8j. By writing the perturbations in polar coordinates, we find that the Jacobian terms in Eq. (4) become , both of which are constant matrices. Thus, according to Eq. (4), the master stability function can be obtained by solving a characteristic polynomial equation and has the following form 28 : Figure 3a shows Λ(α) for c 1 = −1.8 and c 2 = 4, which clearly has Λ″(α 0 ) < 0 at its first zero α 0 ≈ 3. For Stuart-Landau oscillators coupled on temporal networks, Eq. (3) dictates the stability of individual perturbation modes and can be written as where B i ðtÞ ¼ À2 À σλ i ðtÞ c 1 σλ i ðtÞ À2c 2 À c 1 σλ i ðtÞ Àσλ i ðtÞ is periodic with period T ¼ 2π ω (henceforth, we drop the subscript i to ease the notation). According to Floquet theory 34 , the solution to Eq. (12) Fig. 2 Designing temporal networks that synchronize better than optimal static networks. a Evolution of the nonzero Laplacian eigenvalues described in Eq. (5), which are split into two degenerate groups. b Temporal network constructed from the Laplacian eigenvalues in a. The weight of each edge is represented by its thickness. In addition, edges whose weight is larger than 1 n are colored orange, whereas those with weight less than 1 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 different parts of the network are being strengthened in an alternating fashion. must be of the form e μt P(t), where P(t) has period T. The Floquet exponents μ 1 and μ 2 can be extracted by finding the principal fundamental matrix, and their real parts are the corresponding Lyapunov exponents 35 . Figure 3b shows the maximum Lyapunov exponent Γ ¼ maxfReðμ 1 Þ; Reðμ 2 Þg as a function of ω for different values of the temporal activity A. (It is clear from Eq. (8) that all transverse perturbation modes have the same Γ. Thus, Γ is also the maximum transverse Lyapunov exponent and determines the synchronization stability.) We set the coupling strength to slightly below σ c at σ = 2.9 so that no static network can synchronize. As the temporal activity A is increased, Γ becomes negative for an increasingly wide range of switching rate ω, signaling that the temporal variation in the network structure is successfully stabilizing synchronization under the given coupling budget.
Since the only difference between Eqs. (3) and (4) is the periodic λ(t) vs. the fixed α, it is natural to expect the stability of the temporal network to be related to the master stability function averaged over a suitable range of α. Specifically, one might reasonably associate Γ with the averaged master stability function 18,22,26,27,30 Λ ¼ where W(λ) is the probability distribution of λ (it follows that R λ max λ min WðλÞ dλ ¼ 1). However, it is clear that Λ cannot be used to predict Γ in general. One immediate observation is that Λ does not depend on the rate in which λ(t) is changing (it only depends on the distribution of λ), whereas the curves representing Γ in Fig. 3b clearly depend on the switching rate ω. Indeed, in order to go from Γ to Λ, we are required to shuffle B(t) temporally in Eq. (12). This operation is forbidden when the matrices fBðtÞjt 2 Rg do not commute (or, equivalently, when fBðtÞjt 2 Rg cannot be simultaneously diagonalized). To see why, we can look at the formal solution to Eq. (12) expressed in terms of the matrix exponential: where Ω(t) is given by the Magnus expansion 36 : Here, BðτÞ; Bðτ 0 Þ ½ ¼ BðτÞBðτ 0 Þ À Bðτ 0 ÞBðτÞ is the matrix commutator. Equation (15) makes it clear that {B(τ)|0 < τ < t} can be shuffled without affecting Ω(t) if and only if BðτÞ; Bðτ 0 Þ ½ ¼ 0 for all τ 0 <τ < t, in which case everything on the right-hand side except the first term vanishes.

However,
Λ is still extremely informative on whether a given temporal network can synchronize or not. In particular, for ω → 0 (i.e., slow-switching networks 30 ), Γ approaches the value of Λ, as demonstrated in Fig. 3b. Intuitively, this can be understood through a process we call "grow and rotate". When the matrices fBðtÞjt 2 Rg commute, η can be decomposed into components that grow independently along the eigendirections of B(t), whose growth rates are dictated by the corresponding eigenvalues. Eventually, the component along the direction with the largest eigenvalue becomes dominant. However, when fBðtÞjt 2 Rg do not commute, the growth along the eigendirections are often "interrupted", as the eigenvectors of B(t) are no longer fixed and will rotate over time. To keep track of the growth of the dominant component, we must project η onto the new dominant eigendirection upon rotation. These frequent projections can significantly influence the asymptotic growth rate (this is also why the maximum Lyapunov exponent is usually not the mean of the maximum local Lyapunov exponents). At the slow-switching limit, η can grow along an eigendirection uninterrupted for long enough that the effect of the projections becomes negligible. In this case, Γ is determined by the average growth 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 oscillator models 26,30 . As a result, Λ″(α 0 ) < 0 is a robust indicator that synchronization in a system can benefit from temporal networks. This observation echoes recent results in ref. 30 , which demonstrates the importance of a master stability function's curvature for synchronization in the special case of networks with fixed topology and time-varying overall coupling strength. To see why curvature has such a critical role, we assume the temporal variation of λ around 1 to be small and Taylor expand Λ(α) around α 0 . Then, the averaged master stability function for coupling strength σ = σ c is Thus, if Λ″(α 0 ) = Λ″(σ c ) < 0, then Λ < 0 at σ = σ c and stability is guaranteed to be improved at the slow-switching limit, where Γ ¼ Λ. This improvement is expected to extend into the intermediate switching rate due to the continuity of Γ as a function of ω. At the other limit, for ω → ∞ (i.e., fast-switching networks), Γ clearly does not match with Λ. In particular, Γ does not depend on the temporal activity A. For the system in Fig. 3b, Γ approaches Λ(σ) as ω → ∞, which is the value expected for an optimal static network at coupling strength σ (in this case the time-averaged network is a complete graph with uniform edge weights). The mapping from a temporal network to its timeaveraged counterpart at the fast-switching limit is intuitive and well established in the literature 17,23,25 .
The results above provide new insights into the intriguing phenomenon that certain temporal networks only synchronize for intermediate switching rate 21,26,27 : when switching is too fast, the temporal network reduces to its static counterpart and one cannot take full advantage of the temporal variation in the connections; when switching is too slow, although the asymptotic stability might be maximized, the system would have desynchronized long before the network experiences any meaningful change. Thus, the sweet spot often emerges at an intermediate switching rate.
In Fig. 4, we show typical trajectories of n = 11 Stuart-Landau oscillators on the temporal networks described by Eq. (9), with the temporal activity set to A = 0.15. Systems in all three panels are initiated close to the synchronous state, and their only difference lies in the switching rate ω, which allows us to compare networks with static, moderate-switching, and fast-switching topologies. By monitoring the synchronization error Δ(t), defined as the standard deviation among Z j (t), we see that only the system with an intermediate switching rate (ω = 1, panel b) can maintain stable synchrony. Interestingly, Δ(t) in that system goes down non-monotonically and is bounded from above by periodic envelopes. The width of each envelope is 2π, which coincides with the period of the changing network topology.
Universal stabilization of low-dimensional maps. The framework developed so far can be readily transferred from differential equations to discrete maps, from continuous variation in the network topology to discrete switching, and from periodic oscillator dynamics to chaotic ones. The discrete-time analog of Eq. (1) can be written as To demonstrate the advantage of temporal networks in these settings, we focus on the following class of coupled onedimensional discrete maps: where F : R ! R is the mapping function. As we show below, this setup allows us to develop an elegant theory that offers new insights. Similar to the continuous-time case, the synchronization stability is determined by the decoupled variational equations For fixed λ, the Lyapunov exponent of Eq. (19) is given by constant. Thus, the master stability function has a universal form (illustrated in Fig. 1) Taking the second derivative with respect to α, we see that Thus, synchronization in any system described by Eq. (18) can benefit from the temporal networks designed in this paper. In particular, this holds for any mapping function F, which encompasses important dynamical systems such as logistic maps, circle maps, and Bernoulli maps. For concreteness, we set FðxÞ ¼ sin 2 ðx þ π=4Þ and β = 2.8 (the corresponding Γ s = − 0.5855), which models the dynamics of coupled optoelectronic oscillators 37 and exhibits chaotic dynamics. The time-discretized version of the temporal networks described by Eq. (7) works out-of-the-box for the optoelectronic oscillators, despite the vastly different node dynamics. Here, to demonstrate the flexibility of our network design, we consider the following slightly modified switching scheme, which is also more natural for discrete-time systems: where ⌊⋅⌋ is the floor function. Basically, the network switches between two configurations every T iterations, with each configuration being the extremal in the continuous scheme described by Eq. (9). Consequently, every nonzero eigenvalue of the temporal Laplacian alternates between 1 + 2A and 1 − 2A with period T. Again, the averaged master stability function Λ accurately predicts the stability of the temporal network at the slowswitching limit. More interestingly, for systems described by Eq. (18), the connection is much stronger: Λ determines the stability of the temporal network for all switching periods T. To see why, we note that the synchronization stability is determined by the limit product Q 1 t¼1 β À σλ½t À Á F 0 ðs½tÞ. Normally, these are matrix products and cannot be reordered. However, since 1 × 1 matrix multiplications commute, for one-dimensional maps we can reorder them to obtain This independence of Γ on T might seem contradictory to the fact that, at the fast-switching limit, temporal networks can be reduced to their static counterparts. But notice that there is usually no fast switching in discrete-time systems-even if the network topology changes at every iteration, it is still evolving at the same timescale as the node dynamics. Moreover, unlike in continuous-time systems 16,17,23,25 , the discrete nature of the dynamics precludes the use of the averaging techniques 16,25 essential for connecting fast-switching networks with their timeaveraged counterparts. Thus, one cannot map a temporal network to its time-averaged counterpart in discrete-time systems even when the network topology changes much more rapidly than the node dynamics. In Fig. 5, we show the maximum transverse Lyapunov exponent Γ of the synchronization state in the optoelectronic system for σ = 1, which is slightly below σ c . The dashed line corresponds to the theoretical prediction of Γ based on the averaged master stability function Λ ¼ 1 2 ln j1 þ 2A À βj þ ln j1 À 2A À βj À Á þ Γ s . As expected, the static network (A = 0), despite being optimal, is unstable. As the temporal activity A is increased, Λ deceases and synchronization is eventually stabilized. On the other hand, the solid lines represent Γ obtained numerically by evolving Eq. (19) for different switching periods T. These lines are shifted vertically by different amounts in Fig. 5, purely as an aid to the eye. The unshifted versions are shown in the inset. Notice that all the lines collapse onto a single curve, demonstrating the excellent agreement between theory and simulations.
An interesting question is what happens when we introduce random fluctuations to the network structure at each time step t, which makes the temporal network aperiodic and the graph Laplacians noncommutative. In Fig. 6, through direct simulations 35 , we show that temporal networks still outperform optimal static networks in the presence of these random fluctuations. Here, we use the same model of optoelectronic oscillators and the discrete-switching network considered in Fig. 5, except that independent random Gaussian perturbations of zero mean and standard deviation 0.1/n (10% of the average edge weight) are added to the strength of each edge at every 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), synchronization is stabilized at the same coupling strength by the variation in network structure. The network size is set to n = 99 and the switching period to T = 10 in our simulations, although the results do not depend sensitively on these two parameters.

Discussion
To summarize, we have designed temporal networks that synchronize more efficiently than optimal static networks. These temporal networks are particularly relevant when the coupling budget available in a system to maintain stable synchrony is limited. We provided analytical insight into the synchronizability of commutative temporal networks by linking it to the curvature of the corresponding master stability function. In particular, our analysis reveals the subtle relation between the performance of a temporal network and its switching rate. The switching rate has an especially critical role in systems with high-dimensional oscillator dynamics, and networks with intermediate switching rates often emerge as the most effective.
Our open-loop design has several advantages compared with closed-loop schemes where the network structure is adjusted Fig. 5 Temporal networks promote synchronization universally in discrete-time systems with low-dimensional node dynamics. The maximum transverse Lyapunov exponent Γ decreases as the temporal activity A is increased. The dashed line represents the theoretical prediction based on Λ, whereas the solid lines (shifted vertically for visibility) are calculated directly from Eq. (19) for different switching periods T. Without the shift, all curves completely overlap (inset), which confirms our prediction that the stabilization provided by temporal networks does not depend on the switching rate for coupled one-dimensional maps. Fig. 6 Improved synchronization in aperiodic and noncommutative temporal networks. The temporal networks are based on the discreteswitching networks given by Eq. (22), which are further made aperiodic and noncommutative by applying random Gaussian perturbations of zero mean to the strength of each edge independently at every time step t. The standard deviation of the perturbations is fixed at 0.1/n (10% of the average edge weight). The spacetime plots show the evolution of the optoelectronic oscillators on temporal networks with a: temporal activity A = 0 and b: temporal activity A = 0.15. Both systems are initialized close to the synchronous state. Synchronization persists only in the second system, even though the network in the first system is an optimal static network (a complete graph with uniform edge weights). Other parameters are set to n = 99, T = 10, β = 2.8, and σ = 1.05.
on-the-fly based on feedback from the node states (often modeled by adaptive networks 38 ). First, our design does not depend sensitively on the node dynamics. As we have shown, the same design works for systems with vastly different node dynamics, and it applies readily to both continuous-time and discrete-time systems. Second, we do not need to monitor all the nodes constantly, which also eliminates the possibility of being detrimentally influenced by measurement errors. Third, the evolution of the network is highly predictable and we can easily control the coupling budget allocated to the system at any given time t, a task that is far more difficult in adaptive networks. On the other hand, closed-loop schemes have the advantage of being readily adaptive to the changing environment and can react quickly to unexpected perturbations 39,40 . A promising future direction would be to devise hybrid schemes that combine the best from both worlds, which could enable even more efficient and robust synchronization.
In this work, for the sake of analytical tractability, we mostly focused on temporal networks whose Laplacian matrices from different time instants commute. There is evidence that synchronization in temporal networks can benefit when LðtÞLðt 0 Þ≠Lðt 0 ÞLðtÞ 20 . It would therefore be interesting to see whether our design of temporal networks could be further optimized by allowing noncommuting Laplacian matrices. In particular, can random fluctuations in the network structure (which give rise to noncommuting Laplacian matrices in general) outperform our designed temporal networks? More generally, do optimal temporal networks exist for the purpose of synchronization, just like there are optimal static networks? And if so, what are their defining characteristics?
Finally, we hope our results can serve as an important step towards achieving efficient synchronization in complex interconnected systems. For example, many temporal networks arise naturally in the real world through moving agents, whose interactions depend on their spatial distance [41][42][43][44] . An exciting next step is to understand how our design can be implemented in such systems and how the time-varying connections can be translated into the spatial movement of individual agents.

Data availability
All data needed to evaluate the conclusions are presented in the paper. Additional data related to this paper may be requested from the authors.

Code availability
Code for performing network dynamics simulations and stability calculations are available at https://github.com/y-z-zhang/temporal_sync. An archived version of the code is also provided 35 . Additional source code may be requested from the authors.