Nonautonomous driving induces stability in network of identical oscillators

Nonautonomous driving of an oscillator has been shown to enlarge the Arnold tongue in parameter space, but little is known about the analogous effect for a network of oscillators. To test the hypothesis that deterministic nonautonomous perturbation is a good candidate for stabilising complex dynamics, we consider a network of identical phase oscillators driven by an oscillator with a slowly time-varying frequency. We investigate both the short- and long-term stability of the synchronous solutions of this nonautonomous system. For attractive couplings we show that the region of stability grows as the amplitude of the frequency modulation is increased, through the birth of an intermittent synchronisation regime. For repulsive couplings, we propose a control strategy to stabilise the dynamics by altering very slightly the network topology. We also show how, without changing the topology, time-variability in the driving frequency can itself stabilise the dynamics. As a by-product of the analysis, we observe chimera-like states. We conclude that time-variability-induced stability phenomena are also present in networks, reinforcing the idea that this is quite realistic scenario for living systems to use in maintaining their functioning in the face of ongoing perturbations.


I. INTRODUCTION
Many natural and man-made systems consist of interacting oscillatory units [1][2][3]. Such coupled oscillators can yield a plethora of complex dynamical behaviours. Of utmost importance is the case of synchronous dynamics, when all oscillators self-organise and oscillate identically [4,5]. Synchronous dynamics can be crucially needed, such as for the beating of the heart, and can also be severely detrimental, such as in epileptic seizures [8,9]. Phenomena such as ageing [6] and anaestesia [7] are known to alter synchrony of dynamics. These systems can be divided into two classes, based on the nature of their interactions being either attractive, where all oscillators tend to be synchronous, i.e. have the same phase, or repulsive, where the oscillators tend to be not synchronous [10]. The Kuramoto model [11] has served as a paradigmatic model for such systems, and many later works have extended the theory of synchronisation to networks with complex topologies [12][13][14]. Coupled oscillator models have been applied successfully in diverse areas of the sciences [3,15], including circadian rhythms [16], cardio-respiratory dynamics [17][18][19], and metabolic oscillations [20].
Beyond the internal dynamics, many systems are subject to ongoing external influence from their environment. That external influence drives the system, and its action can, in general, change with time. The day-light cycles drive the suprachiasmatic nucleus neurons in our brain, which in turn control all circadian rhythms in our * m.lucas@lancaster.ac.uk † duccio.fanelli@gmail.com ‡ aneta@lancaster.ac.uk bodies, making sure the body's different clocks tick together [21]. Day-light cycles are subject to seasonality. Moreover, the neurons's firing is influenced by intra and extracellular ionic concentrations, which in turn affects synchronisation [22]. In some cases, like for fish communities in Japan [23], the time-variability of the system yields time-varying properties and stability [24,25], which was found crucial to the functioning of the community. Moreover, in real-world systems, short-term properties are sometimes more physically relevant than the asymptotic or long-term ones [26]. Non-driven networks with time-varying parameters have started to attract attention recently [27][28][29][30][31]. Driven networks of identical oscillators have been studied theoretically in specific cases, such as regular topologies [32], or hierarchical networks with only a few oscillators being driven [33,34]. However, to date, those few studies either considered fixed-frequency driving, or modelled the time-variability statistically as just noise. Here, we consider a deterministically but aperiodically time-varying frequency of the external driving, playing the role of the ever changing environment. We bring insight about the effect of the timevariability on the stability -both short-and long-term -of synchronous solutions, analytically and numerically. Surprisingly, we show that adding time-variability in the driving can increase the stability of the driven network. In the attractive coupling case, we generalise results from [26,35] to networks. Weakly coupled non-linear oscillators with slowly varying parameters can be reduced to a phase model [36][37][38], which suggests that our results apply to many other types of oscillators. Finally, our results are valid for networks of any size and any topology.
The paper is organized as follows. In Sec. II, we introduce our driven network of oscillators model. In Sec. III, we provide a general linear stability analysis of the synchronous solutions. First, in Sec. III A, the analysis is performed for the fixed-frequency driving case. Then, in Sec. III B, the analysis is extended to the general nonautonomous case, where the driving frequency is timedependent. These results are then discussed and confirmed numerically in Sec. IV for the two possible cases. For the attractive and repulsive cases, in Sec. IV A and Sec. IV B respectively, we show how the stability region is enlarged and describe short-time properties. Furthermore, for the repulsive case, in Sec. IV B, we evidence a link between topology and dynamics, and suggest a control strategy based on it. Moreover, we show how timevariability can counter-balance the desynchronising effect of topology and enhance synchrony. As a by-product, we observe chimera-like state. Finally, in Sec. V we conclude with a brief summary of the results.

II. MODEL
We consider a driven network of N identical oscillators defined by phases θ i and frequency ω, (1) for i = 1, . . . , N , with coupling constant D, and where A stands for any (undirected) adjacency matrix with elements A ij ∈ {0, 1}. Each oscillator is driven with strength γ ≥ 0 by the same external oscillator with phase θ 0 (t) and time-varying frequencẏ where ω 0 is the non-modulated frequency, f is a bounded function, and k and ω m are the amplitude and frequency of the imposed modulation, respectively. Note that f (ω m t) is a generic function, and need not be periodic; without loss of generality, we bound its image in [−1, 1]. System (1)-(2) is a direct generalisation to networks of the single driven oscillator system presented in [35]. The non-driven network, γ = 0, is an autonomous system. In this case, a (fully) synchronous solution always exists, i.e. a solution where all oscillators are in the same state at all times, θ i = θ j for all i = j. This is as long as the network is connected, meaning that there is a path between any two nodes (oscillators). The synchronous solution is stable (unstable) if the coupling between oscillators is attractive (repulsive), i.e., D < 0 (> 0) [39]. In the repulsive case, the system is attracted to an incoherent state.

III. THEORETICAL ANALYSIS
In the driven system, γ = 0, synchronous solutions also exist. It is convenient to go to the rotating frame of the driving, ψ i = θ i − θ 0 (t), where system (1) is rewritteṅ is the time-dependent frequency mismatch. We refer to the θ i as the phases, and to the ψ i as the phase differences between the phase of the driving and that of the i-th oscillator. System (3) is nonautonomous due to the time-dependent frequency mismatch, and is hence hard to treat in general. In the rest of the study, we assume ω 0 is modulated slowly, i.e., ω m ω 0 , and use an adiabatic approach to study the existence and stability of synchronous solutions.

A. Autonomous case
To that end, we first consider the simpler case of a constant driving frequency, recovered for k = 0, for which ∆ω(t) = ∆ω = ω − ω 0 is constant and system (3) is autonomous. A synchronous solution, ψ i =ψ for all i = 1, . . . , N , always exists, and obeyṡ which is the so-called Adler equation describing a single driven oscillator [4]. From Eq. (4), two types of synchronous solutions can be identified: synchronised (to the driving), or not . First, if γ ≥ |∆ω| is satisfied, there exists a stable fixed point. This fixed point corresponds to the synchronous synchronised (SS) solution, ψ ss = π − arcsin(−∆ω/γ), characterised by a constant phase difference,ψ ss = 0. Second, if the synchronisation condition is not met, γ < |∆ω|, there exists a synchronous not synchronised (SNS) type of solutions, denoted ψ sns , which grows (or decays) monotonically, ψ sns > 0 (or < 0). We now investigate the linear stability against a small heterogeneous perturbation δψ i around those solutions. This is determined by linearising Eq. (3) for each node around a solutionψ(t), which stands for either the SS or SNS, Here, L ij = A ij − K i δ ij denotes the Laplacian matrix of the network, defined in terms of the connectivity of each node K i = N j=1 A ij , and the Kronecker delta δ ij . One can now decouple this N -dimensional problem by projecting it onto the eigenbasis of the Laplacian, defined as where the φ α are the eigenvectors associated to the eigenvalues Λ α , for α = 1, . . . , N . The latter are non-positive and real, since the network is assumed to be undirected, and ordered The perturbation can be decomposed in that basis, and we look for solutions of the form δψ i (t) = N α=1 c α e t 0 λ α (t ) dt φ α i . Plugging into Eq. (5) and solving for each α yields the instantaneous Lyapunov exponent (LE) spectrum which is completely general and also valid in the nonautonomous case, as will be seen later. Now in the autonomous case,ψ(t) is periodic modulo 2π [4] and hence (long-term) Lyapunov exponent spectrum is well defined as the time-average of the instantaneous values An explicit form of formula (7) is obtained by splitting it into two cases For the synchronised case, γ > |∆ω|, the explicit solution ψ ss = π−arcsin(−∆ω/γ), which is constant, was plugged into formula (7). For the not synchronised case, the solution is ψ sns , for which the averaging term in formula (7) vanishes.

B. Nonautonomous case
In general, the driving frequency is time-dependent, k = 0, and synchronous solutions obey a nonautonomous version of Eq. (3)ψ = ∆ω(t) + γ sinψ, (9) which was studied in [35,40]. As previously mentioned, here and throughout the paper, we assume slow modulation of the frequency, i.e. ω m small. Equivalently, ∆ω(t) varies much more slowly than the dynamics of the system. Hence, there is a separation of timescales: ∆ω(t) is the slow variable, andψ(t) the fast one. Thus, over the fast timescale, the frequency mismatch is quasistatic, and the slowly moving point attractor ψ ss (t) = π − arcsin(−∆ω(t)/γ) is followed adiabatically, when it exists [40]. This state corresponds to the SS solution of the autonomous case, and differs with it in being only quasi -constant, i.e.,ψ ss = 0 over the fast timescale, and existing only at times such that γ > |∆ω(t)|. Consequently, in contrast to the autonomous case, Eq. (9) shows three types of synchronous solutions: two of them correspond to the SS and SNS solutions as discussed in the autonomous case, whereas a third type exhibits intermittent synchronisation. Three regions in parameter space can thus be defined, each corresponding to the existence of one of those types solutions [35], as illustrated in Fig. 1. In region I, the condition γ > |∆ω(t)| is not met at any time. The solution grows (or decays) monotonically and we denote it ψ sns (t). In region II, the condition γ > |∆ω(t)| is met at all times. The solution is denoted by ψ ss (t) and has an approximately null time derivative, as described above. In region III, γ > |∆ω(t)| is met only at certain times. The phase difference alternates between times of growth (or decay), and quasiconstant (bounded) epochs. We call it synchronous intermittently synchronised (SIS), and denote it ψ sis (t). The three types of solutions are illustrated in Fig. 2(a) and will be discussed in the next section.
The adiabatic assumption allows us to obtain the following formula for the adiabatic LE spectrum, as a nonautonomous version of formula (8), by retracing the same reasoning Equation (10) enables us to draw conclusions about the stability of the aforementioned states. The stability of the SS solution can be assessed by the first condition in Eq. (10). The stability of the SNS is determined by the second condition in Eq. (10). Finally, the stability of the SIS is determined by the first condition of Eq. (10) at times such that γ > |∆ω(t)|, and the second condition the rest of the time.
Note that this result holds true regardless of the network size, topology, or the shape of the frequency modulation function f (ω m t).

IV. ON THE STABILITY OF THE SYNCHRONOUS SOLUTION
We now examine the stability of the synchronous solutions further via formula (10) in the attractive and repulsive cases.

A. Attractive case
In this case, D < 0, all oscillators tend to display the same phase, when not driven. For a given network, only the largest LE λ max (t) determines the stability. In this case, it is λ max (t) = λ 1 (t) which corresponds to Λ 1 = 0 in formula (10), and reads a condition which is identical to that obtained for the single oscillator case considered in [35], even though we have now an arbitrary network of oscillators.
For the sake of clarity, we set f (ω m t) = sin(ω m t) in numerical examples. Note, however, that the analysis is independent of the explicit form of f (ω m t), which can in general be aperiodic, or even defined only over a finite timespan, as discussed in [26].
Typical dynamics of the ψ i is shown for the three regions in Fig. 2(a). The SS (solid) and SNS (dotted) are bounded and unbounded, respectively, as mentioned in the previous section. For the SS (SNS), γ > |∆ω(t)| (γ < |∆ω(t)|) is always met, and hence the LE is negative (zero), i.e., stable (neutrally stable) at all times, as seen from Eq. (11). However, the SIS (dashed) stays bounded only intermittently. In this example, the timescale that controls the alternation of those epochs is the period of the imposed modulation, T m = 2π/ω m , which could be arbitrarily long. Note that these epochs of growth are not the typical 2π phase slips of the fixed-frequency single oscillator case [4], observed close to the synchronisation border. Here, the "slip" is a drift caused by the temporary neutral stability, and its relative importance depends on the length of that neutrally stable epoch, as seen from Fig. 2(b). Now, the SIS is analysed further as shown in Fig. 2(b)-(d). Figure 2 (b) shows the adiabatic LE (dashed black) (11). Epochs where it is negative (null) correspond to the ψ i growing (staying bounded). Moreover, Fig. 2 (b) shows the agreement between the instantaneous LE (6) (grey) and the adiabatic LE (11). This confirms a posteriori the adequacy of the adiabatic approach. The intermittency can also be seen in the timefrequency representation of sin θ 1 (t) for a trajectory of system (1), as shown in Fig. 2(c). Here, stability epochs correspond to the frequency being entrained by the driver (one frequency mode), whereas neutral stability epochs correspond to the presence of two frequency modes plus harmonics. Finally, a negative maximum LE guarantees convergence of different initial conditions, as long as they are in the basin of attraction of the synchronous solution,  even if the trajectory is only intermittently synchronised to the driving. This is illustrated in Fig. 2(d).
Note that all other non-maximal LEs of the spectrum are negative at all time, λ i (t) < 0, for i = 2, . . . , N , since the corresponding Laplacian eigenvalues are negative. The exponent λ 1 measures the stability against a homogeneous perturbation, whereas the rest of the spectrum corresponds to any heterogeneous perturbation. In other words, any synchronous solution will stay synchronousall oscillators with the same phase -against any perturbation. During stable epochs, even the common phase of the synchronous state is stable. During neutrally stable epochs, however, the common phase of the synchronous solution can be pushed by a homogeneous perturbation, and change without the perturbation decaying or growing.
On average, over a time T , the maximum exponent is Note that the LE in Eq. (12) is strictly zero only in region I where the system does not synchronise to the driving at any time. Indeed, in region III, the adiabatic LE alternates between zero and negative values, and is negative on average. Moreover, region I decreases in size as k is increased, and so the remainder of parameter space, corresponding to stability λ max < 0, grows. In other words, by increasing the amplitude of the nonautonomous modulation, one makes the region of stability larger in parameter space. In this region of stability, different initial conditions converge to one unique trajectory. In Fig. 3, panels (a)-(c) show the enlargement of the negative LE region in parameter space as k increases from 0.0 to 0.4, and 0.8. Panel (d) combines and shows those LE values for all three values of k, but for a single value of the forcing strength γ = 2.5. The region of stability is the union of regions II and III. Region II, where trajectories are always synchronised to the driving, decreases in size as k is increased, but region III grows enough so that their union grows. The phenomenon does not depend on the explicit form of f (ω m t) and was also illustrated nu-merically, in the single driven oscillator case, for different aperiodic f (ω m t), and f (ω m t) defined over a finite time in [26,35].

B. Repulsive case
In this case, D > 0, all oscillators have a tendency towards pairwise asynchrony, when not driven. In the unforced network, the synchronous solution is always unstable, but if forced, the synchronous solution can be stable or unstable. The largest adiabatic LE, in this case λ max (t) = λ N (t) corresponding to Λ N < 0 in formula (10), is given by Firstly, when the driving strength is not large enough, i.e., γ < |∆ω(t)|, synchronous solutions are unstablei.e., always for the SNS and during the non-synchronised epochs for the SIS. This is in contrast with the attractive case, which exhibited neutral stability under the same condition. Secondly, when γ > |∆ω(t)|, two effects compete: the network couplings push oscillators away from each other (positive first term) whereas the external driving brings them back towards a synchronous state (negative second term). This is again in contrast with the attractive case, which only exhibits stability under the same condition. Those two effects exactly compensate each other if which has a minimum when ∆ω = 0 at γ = −DΛ N . This is shown for the autonomous case k = 0 in Fig. 4(a) (dashed white).
In the autonomous case, k = 0, the frequency mismatch is constant in time. Hence, condition (14) divides region II into two subregions based on stability: regions IIa (IIb) where the SS is unstable (stable) corresponding  to γ values smaller (larger) than that of condition (14). The regions are shown in parameter space in Fig. 4(a), and are confirmed by the computation of the λ max , as shown Fig. 4(b)-(d) for different values of Λ N . As the largest-magnitude eigenvalue Λ N of the Laplacian, or alternatively the network coupling strength D, is increased, the region of negative λ max decreases in size.
This observation can be used as a viable control strat-egy. Indeed, note that the eigenvalues of the Laplacian, and in particular Λ N , are determined by the topology of the network considered. Moreover, the inequality |Λ N | ≤ 2K max holds, where K max is the number of connections of the most connected node [41]. In other words, the stability of the dynamics is directly determined by the topology, and in particular by the connectivity of the most connected node. Here, in the repulsive case, more connections between nodes, and in particular to the most connected one, means less stability. So then, one can optimally decrease the absolute value of Λ N -or equivalently increase the region of stability -by cutting the edges of the most central node, which in turn amounts to reducing K max . This is illustrated in Fig. 5 where the SS is stabilised by cutting only 5 chosen links (red) out of 81 (≈ 6%). The original network is a Barabási-Albert one [see Fig. 5(b)] [42], for which this strategy is most effective, since only a few nodes have very high connectivity.
As a by-product of the analysis, in region IIa, we numerically observed partially-locked states [43], or chimera-like states [44,45], where most of the oscillators are phase-locked to the driving in a quasi-synchronous cluster, while the rest of the oscillators drift independently. An example video of such dynamics is provided in the Supplementary Material.
In general, in the nonautonomous case, k = 0, an additional type of intermittent synchronisation can be exhibited, which here alternates between stability and instability. That is, the adiabatic LE for the synchronous solutions alternates between negative and positive values. This happens in a new region IIc which, as k is made progressively larger than 0, appears at the border between regions IIa and IIb and subsequently grows, as illustrated in Fig. 6(a)  gion IIc can be defined as all pairs (γ, ∆ω) such that γ is greater than the value of condition (14) at certain times and smaller at others. Equivalently, it is {(γ, ∆ω) : ∃ t : γ = (DΛ N ) 2 + ∆ω 2 (t)}. The region is constant in time, and the explicit form of its boundaries is uninformative and is hence omitted here. In region IIc, the long-term LE can be either positive or negative. For small enough values of the coupling strength D, the region of negative LE increases with k, just as in the attractive case of Sec. IV A. This effect can be achieved for relatively small values of the coupling strength D.
However, while alternation between neutral stability and stability guaranteed overall stability and convergence in the previous section, alternation between stability and instability does not. Here, in region IIc, for a synchronous initial condition, a negative LE associated to its trajectory does not guarantee the state will stay synchronous and converge to the SS solution, as shown in Fig. 7(b). Indeed, during epochs of instability, any perturbation can push the state far away from its original synchronous solution exponentially fast. However, if no other attractor exists during the epochs of stability, then the system will return to the synchronised synchronous state when the next epoch of stability occurs. Thus it alternates between the two regimes, as shown numerically in Fig. 7(b). This effect is purely due to the time-variability of the driving, controlled by k. Moreover, there is subregion in parameter space which is unstable (region IIa) when k = 0, but turns intermittently stable (region IIc) for k > 0. This can be used as a control strategy as illustrated in Fig. 7: one forces a completely incoherent state [panel (a)] to be synchronous and synchronised indefinitely often and long, by allowing time-variability in the driving frequency. Thus, time-variability can also be used to counter-balance the desynchronising effect which is produced by a highly connected network, as shown in Fig. 5. Finally, this strategy can be useful when the requirement for synchronicity are not so stringent and the system under inspection does not need to be synchronous at all times.

V. SUMMARY AND CONCLUSIONS
In this work, we studied the effect of driving an arbitrary network of identical phase oscillators with an external time-varying-frequency oscillator. This is in contrast with previous studies of forced networks, which do not consider time-variability, a crucial ingredient in realworld systems. Stability -both short-and long-termof the synchronous solutions was assessed by linear stability analysis, and results were confirmed via numerical simulations. The system studied is nonautonomous and hard to treat in general, hence we assumed slow timevariability. Two cases were treated: attractive and repulsive couplings, where oscillators tend to be in phase and out of phase, respectively.
In the attractive case, we showed that increasing the amplitude k of the frequency modulation enlarged the stability region in parameter space, as defined by a negative Lyapunov exponent. In other words, a more variable driving makes the stability of synchronous behaviour more robust against parameter changes. This is the direct generalisation to networks of a result that was obtained for a single driven oscillator. An additional region of intermittency between stability and neutral stability appears, as a result of the time-variability in the frequency.
In the repulsive case, we first demonstrated a control strategy, where the synchronous solution can be made stable by cutting a few chosen links in the network. Moreover, we established the phase diagram of the synchronous solutions. We then showed that one could counter-balance the desynchronising effect of connectivity by allowing time-variability in the driving frequency. Indeed, we showed that the time-variability of the frequency can drive the system back to a synchronous state intermittently, where it was asynchronous with a fixedfrequency driving. This observation can serve as an alternative control protocol. Finally, as a by-product of the analysis, we numerically observed chimera-like states.
Dynamics of the phase and frequency is illustrated in detail in the different cases. Such classification can potentially be of use to experimentalists who can only mea-sure phase and frequency, and may have only limited or no knowledge of an external driving. Moreover, we illustrated diverse and simple control strategies to enhance synchronisability that could explain how living systems maintain stability in a changing environment, and could also be implemented directly.
Finally, many of the interesting and physically relevant features observed in this system were either happening over finite time, or explained by the finite-time analysis of dynamics. Consequently, a solely asymptotic analysis would have missed much of the dynamical intricacies at hand. From this, we conclude that nonautonomous networks of coupled oscillators, such as the present one, exhibit features reminiscent of those observed in living systems, and that a finite-time approach is crucial to their understanding.