Stabilisation of dynamics of oscillatory systems by non-autonomous perturbation

Synchronisation and stability under periodic oscillatory driving are well-understood, but little is known about the effects of aperiodic driving, despite its abundance in nature. Here, we consider oscillators subject to driving with slowly varying frequency, and investigate both short-term and long-term stability properties. For a phase oscillator, we find that, counter-intuitively, such variation is guaranteed to enlarge the Arnold tongue in parameter space. Using analytical and numerical methods that provide information on time-variable dynamical properties, we find that the growth of the Arnold tongue is specifically due to the growth of a region of intermittent synchronisation where trajectories alternate between short-term stability and short-term neutral stability, giving rise to stability on average. We also present examples of higher-dimensional nonlinear oscillators where a similar stabilisation phenomenon is numerically observed. Our findings help support the case that in general, deterministic non-autonomous perturbation is a very good candidate for stabilising complex dynamics.

Synchronisation and stability under periodic oscillatory driving are well-understood, but little is known about the effects of aperiodic driving, despite its abundance in nature. Here, we consider oscillators subject to driving with slowly varying frequency, and investigate both short-term and long-term stability properties. For a phase oscillator, we find that, counter-intuitively, such variation is guaranteed to enlarge the Arnold tongue in parameter space. Using analytical and numerical methods that provide information on time-variable dynamical properties, we find that the growth of the Arnold tongue is specifically due to the growth of a region of intermittent synchronisation where trajectories alternate between short-term stability and short-term neutral stability, giving rise to stability on average. We also present examples of higher-dimensional nonlinear oscillators where a similar stabilisation phenomenon is numerically observed. Our findings help support the case that in general, deterministic non-autonomous perturbation is a very good candidate for stabilising complex dynamics.

I. INTRODUCTION
Complex oscillatory dynamics abounds in nature. Despite many real-life examples exhibiting stable oscillations with a time-varying frequency (e.g., [1][2][3]), little is known theoretically about the properties of this type of behaviour. This kind of oscillation requires aperiodic external driving, making the system non-autonomous by nature [4], such that most of the traditional analytical methods are unusable or insufficient. The case of periodic forcing with a constant frequency, which has been extensively investigated to date, is often too simplistic to account for reality.
Closely linked to the concept of stability is the concept of synchronisation. Synchronisation phenomena in the sciences, including phase synchronisation of complex oscillators, have drawn much attention over the last decades [5]. Different mechanisms for achieving phase synchronisation, such as synchronisation by periodic forcing [6], by noise [7], and by quasi-periodic forcing [8], as well as synchronisation of chaotic oscillators [9,10], and networks of oscillators [11], have been considered. At the same time that theoretical interest in synchronisation phenomena has been growing, real-life applications have been found to be prevalent in many diverse aspects of nature [12][13][14], including circadian rhythms [15,16], cardio-respiratory dynamics [1,17,18], metabolic oscillations [3], the brain [19,20], and climate dynamics [21]. But once again, little is understood theoretically about synchronisation in the context of variable-frequency oscillations.
So then, there is a need for extensive and ongoing study of dynamical behaviour under deterministic oscillatory driving with a time-varying frequency, to fill the gap between the existing theory of deterministically driven systems where constant frequency is typically assumed and the statistical theory of systems driven by noise. In this manuscript, we present three major contributions to the field of interacting non-autonomous systems: Firstly, we present notions of stability, synchronisation and instantaneous frequency entrainment in the non-autonomous setting, and the relationships between these concepts; and we investigate these concepts for the simplest example of a phase oscillator subject to driving with slowly time-varying frequency. In so doing, we enable the notion of chronotaxicity [18,22,23] to be broadened beyond its current description, and we compare the stability properties in this setting with the traditional settings of fixed-frequency driving on the one side and driving by stationary noise on the other. Secondly, we introduce an approach to analysing time-dependent dynamical stability from a time-series consisting of time-localised Lyapunov exponents (LE), that is, finite-time Lyapunov exponents (FTLE) taken over a time-window with a moving centre. By contrast, typically, dynamical stability is assessed only in terms of time-averaged stability, for example by the asymptotic LE [24]. Thirdly, numerically and analytically, we show enlargement of the stability region in parameter space for the phase oscillator subject to driving with slowly varying frequency, and we show that this growth is specifically due to the growth of a subregion characterised by intermittent synchronisation where the time-localised dynamical stability is varying. While we show that slow modulation of the driving frequency guarantees enlargement of the stability region for one-dimensional phase oscillators, we also show numerically that the same phenomenon can readily occur in more general oscillatory systems. This mathematical phenomenon of stabilisation has two major practical implications: (i) deterministically varying the frequency of external driving could be implemented as a means of inducing stability in complex systems, and (ii) dynamical systems where stability is induced by deterministic frequency variation are an excellent candidate for modelling living systems, which are highly complex and yet usually operate stably within a time-varying environment.
The few existing pioneering studies of stability and synchronisation with time-varying frequency of oscillations have considered a simple case of linearly growing frequency [25], more general frequency that is slowly varied (as in our present paper) with particular application to both linearly growing frequency and low-passfiltered noise [26], networks of coupled oscillators [27], and a case of two interacting oscillators, each with the same form of time-varying frequency [18,23]. In the lastmentioned studies, the notion of chronotaxicity was introduced, to indicate a type of synchronisation specifically characteristic of non-autonomous interacting oscillators. The central idea in [26] (also relevant to the analytical approach in the work presented here) is that under driving of sufficiently slowly varying frequency, the phenomenon of synchronisation by common driving can be investigated in terms of the presence of a stable equilibrium for the instantaneous vector field; the goal of [26] is to identify what qualifies as slow variation. Time-series analysis methods have also been developed to resolve in time the dynamical characteristics of time-varyingfrequency oscillators (e.g., wavelet-based spectrum, coherence and bispectrum [28] as well as Bayesian inference of coupling functions [29]) rather than analyse them in a statistical sense (e.g. calculating power-spectrum density) and thereby miss noteworthy time-dependent dynamical features. Normal forms of different types of non-autonomous bifurcations have also been investigated [30][31][32], and safety criteria were derived for aperiodically forced systems [33].
However, this present work is the first investigation of non-autonomously driven oscillators showing growth of the stability region in parameter space, where this growth is due to time-variability without the need for statistical phenomena as in noisy models.
The paper is organised as follows. In Sec. II, we introduce a simple one-dimensional phase oscillator model. We then provide an explanation of notions of synchronisation and stability for non-autonomous systems, followed by a theoretical analysis of the one-dimensional model, showing the enlargement of the stability region, as well as the birth of an intermediate region of intermittent synchronisation. We illustrate these phenomena with numerical results for both long-time and short-time behaviour. In particular, in Sec. II E, we discuss the relationship between the deterministic system considered here and the analogous case with noisy driving as considered in previous works. In Sec. III, we illustrate the stabilisation phenomenon numerically in higher-dimensional systems, and argue that it is of general importance. Finally, in Sec. IV we discuss the results, and in Sec. V we provide a brief summary.

A. Model
The system we consider is a driven phase oscillator of the formθ where the driving has strength γ, phase θ 1 (t), and a timevarying frequencyθ where ω 1 is the non-modulated driving frequency, f is a bounded function, and ω m and k are the modulation frequency and relative amplitude, respectively. The phase oscillator θ may represent, for example, the phase on the stable limit cycle of an oscillator satisfying where r p is the amplitude of the limit cycle and is the restoring constant.
The unforced system (1) with γ = 0 is a typical autonomous phase oscillator [34], and hence its phase is neutrally stable (zero LE). In the forced system, i.e. γ = 0, the traditional constant-frequency forcing case is recovered for k = 0. In this case, depending on the parameters, the system lies in one of two regimes: either 1 : 1 synchronisation (negative LE), or neutral stability (zero LE). The condition for synchronisation, γ > |∆ω|, with the frequency mismatch ∆ω = ω 0 − ω 1 , is derived analytically [5] by requiring that the equation for the phase differencė has a stable fixed point. This condition for synchronisation corresponds to a so-called Arnold-tongue [35] in (γ, ω 0 )-parameter space. This Arnold tongue can be seen in Fig. 2(a) as the region appearing in shades of blue, corresponding to negative values of the numerically computed LE. Since (4) is an autonomous differential equation, the numerical LE computed over a long time wil approximate well the asymptotic LE, except possibly when the parameters lie extremely close to the border of the Arnold tongue. For k = 0, the equation for the phase difference is now the non-autonomous equatioṅ ψ = ∆ω(t) + γ sin ψ, with frequency mismatch ∆ω(t) = ω 0 − ω 1 (1 + kf (ω m t)). Throughout this paper, we assume that the modulation is much slower than the dynamics of the system, i.e. ω m is very small.

B. Synchronisation in autonomous and non-autonomous systems
Suppose an oscillatory system θ with no internal timedependence is subject to driving from another oscillator θ 1 with fixed frequency different from the natural frequency of θ. We say that θ is synchronised to θ 1 if over time, the trajectory of θ loses memory of its precise initial phase and instead follows a periodic behaviour whose period is a rational multiple n m of the period of θ 1 . In this case, we say that θ 1 entrains the frequency of θ, and we describe the synchronisation as n : m synchronisation. This implies in particular that the difference in unwrapped phase between an n-fold cycle of θ 1 and an m-fold cycle of θ stays bounded over all time-a phenomenon referred to as phase-locking between θ and θ 1 .
The particular phenomenon that θ(t) loses memory of its initial phase is called phase stability. If θ is a phase oscillator (as in our model), then phase stability can be assessed in terms of the sign of Lyapunov exponent associated to the trajectory θ(t): a negative LE indicates stability. Stability inherently implies resilience against the effects of other possible perturbations not accounted for in the model. When the phase of a driven oscillator is stabilised by a fixed-frequency driving oscillator, typically this implies n : m synchronisation for some integers n and m; we emphasise that this statement is specific to fixed-frequency driving. For n : m synchronisation where (in lowest terms) n ≥ 2, it is not possible for θ(t) to lose all memory of its initial phase: for any 1 ≤ i ≤ n − 1, delaying the initial phase by a suitable amount will delay the phase of the eventual periodic motion by i n . However, in the case of 1 : m synchronisation, it is possible for θ eventually to lose all memory of its initial phase; in this case, we say that the phase is globally stable.
Synchronisation has also been investigated in the context of systems driven by noise, such as zero-mean Gaussian white noise or a pulse train with i.i.d. consecutive waiting times [5,36,37]. For an oscillator θ driven by such noise, one does not have a notion of n : m synchronisation between this driven oscillator and the noise ξ driving the oscillator. This is because, even if the noise is stationary noise, any one realisation of the noise does not have a deterministic periodic behaviour. As described in [5,Section 15.2], instead of defining synchronisation in terms of "phase locking", one can think of synchronisation here as meaning that over time, θ loses memory of its precise initial state and instead follows some path that is determined by the realisation of the noise-but since the noise itself has no deterministic regular behaviour, this phenomenon can only be physically manifested as synchronisation by common noise between copies of θ.
Synchronisation by common noise is a particular case of the phenomenon of synchronisation by common external driving (which may be noisy or deterministic): Suppose we have an array θ 1 , . . . , θ n of self-sustained oscillators whose internal dynamics are described by exactly the same systemθ i = f (θ i ), where f does not depend on i; and no direct coupling is introduced between these oscillators, but instead all these oscillators are simultaneously subject to driving from the same external driver p(t) (which could be noisy or deterministic). Thus the n driven oscillators are now indirectly coupled, and it may happen that as a result of this indirect coupling, over time the trajectories of θ i lose memory of their initial states and instead follow the same path as each other. This may be viewed as a kind of perfectly instantaneous 1 : 1 synchronisation between the driven oscillators.
The relationship between the above concepts is as follows. For a self-sustained phase oscillator θ subject to driving by an external driver p(t), the following statements are equivalent: (i) the trajectory of θ is globally stable; (ii) the trajectory of θ eventually follows some path that is determined by p(t) independently of the initial phase of θ; (iii) any array of identical copies of θ is synchronised when the indirect coupling of common driving by p(t) is introduced; and in the case that p(t) is a fixed-frequency deterministic oscillator, these typically imply that: (iv) there is 1 : m synchronisation between p and θ for some m.
The physical interpretation of the implication (ii)⇒(i) is that the driving p(t) causes θ to become resilient in its course of following the path laid out by p(t), although as we shall see, this resilience may only be intermittent. Such driving-induced resilience may play an important role in many real-world systems that exhibit remarkable stability in the face of continuous environmental perturbations.
In our model, if k = 0 then the driving is a deterministic oscillator θ 1 with non-fixed frequency. Hence, it will be useful for us to discuss notions of synchronisation for oscillators subject to deterministic oscillatory driving with time-varying frequency. Such driving shares in common with fixed-frequency driving that it is deterministic and oscillatory, and it shares in common with noisy driving that it does not possess a phase which proceeds in cycles of a fixed period. Therefore, on the one hand, as with noisy driving, it is not clear that one can correctly define a notion of n : m frequency entrainment, though the slightly weaker phenomenon of n : m phaselocking can still occur; nonetheless, as in [26], one can still consider the question of whether identical copies of the driven oscillator are caused to synchronise by simultaneous driving from the driving oscillator.
Having stated that n : m frequency entrainment is difficult to define in our setting, let us now highlight our slow variation assumption. Under this assumption, one can define a notion of instantaneous frequency entrainment. In general, if a pair of phase oscillators θ, θ 1 is governed by a non-autonomous differential equation and it is assumed that f 1 (t, ·), f 2 (t, ·, ·) vary slowly with t, then we can say that there is frequency entrainment at time t if the solution of the associated autonomous differential equation exhibits frequency entrainment. In the case of our model, at any time t, there is instantaneous 1 : 1 frequency entrainment between θ and θ 1 if and only if the differential Just as negative Lyapunov exponents are connected with the presence of frequency entrainment for fixedfrequency driving, so likewise instantaneous frequency entrainment will typically be connected with negative finite-time Lyapunov exponents defined over a suitable time-window. Finite-time Lyapunov exponents are a measure of stability over finite time-scales. We will use the term time-localised Lyapunov exponent to refer to FTLE taken over a sliding time-window [t, t + τ ] which slides along with time t. By contrast, we will use the term long-term Lyapunov exponent to refer to a Lyapunov exponent taken over a long time-interval [0, T ]; when clear from the context, we will sometimes drop the word "long-term". Technically, a long-term LE is still a finite-time Lyapunov exponent, but it plays a similar role to asymptotic LE for autonomous systems. Asymptotic LE need not exist for non-autonomous systems-indeed, non-autonomous systems need not even be well-defined over infinite-time. But moreover, even if they can be defined, asymptotic LE may not necessarily be physically relevant for the limited time-scales on which a system is considered in practice.

C. Theoretical analysis
In contrast to the autonomous case, the existence of an attracting equilibrium point for the vector field on the right-hand side of Eq. (5) (regarded as a function of ψ) can change with time t; as shown in Fig. 1(a), for a sinusoidal modulation f (·) = sin(·), if k is large enough then within each modulation period the vector field undergoes two saddle-node bifurcations. Since we assume that the modulation is much slower than the dynamics of the system, the system adiabatically follows the moving attracting point ψ slow (t) = π − arcsin[−∆ω(t)/γ] for Eq. (5), when it exists. (A more technically precise description of how ω m needs to compare with the values of other parameters in order to qualify as "small" for the purposes of this adiabatic approach can be found in . The (∆ω(t))-dependent curve ofψ against ψ moves up and down over time, as indicated by the two solid lines representing where the curve could be at two different instants in time. When this curve lies between the dashed lines, the system has an attracting point, and otherwise, not. (b) Phase diagram showing three regimes. Light (region I in the text), medium (region III), and dark grey (region II) show where the system is never synchronising, intermittently synchronising, and always synchronising, respectively. Solid white curves show the border between synchronisation and non-synchronisation for if f (·) is set to 0; and white dashed curves show the border between synchronisation and nonsynchronisation for if f (·) is set to ±1. When k is increased, regions I and II decrease while region III increases; hence in particular, the Arnold tongue consisting of the union of regions II and III increases. [26].) On faster time-scales, one could view ∆ω(t) as approximately constant and consider the stable point in the quasi-stationary limit.
Following the idea that solutions follow the moving attracting point ψ slow when it exists, we derive three regions, with qualitative features corresponding to the following conditions on Eq. (5): (I) no existence of the attracting point at any time t, (II) existence of the attracting point for all t, and (III) alternation over time between the existence and non-existence of the attracting point. If we assume that f (·) varies throughout the interval [−1, 1], then these conditions on the parameters are precisely as illustrated in Fig. 1(b). In region I, the slow variation assumption imlies that solutions behave similarly to the neutrally stable regime of the fixed-frequency-driving system; solutions of (1) or (5) will exhibit neutral stability, with a long-term Lyapunov exponent that is essentially zero. In region II, the attracting point exists at all times, and attracts solutions starting from throughout the circle to itself; thus, the driven oscillator θ is globally stable, losing memory of its initial state and following the motion of θ 1 (t) + ψ slow (t). In particular, long-term LEs will be negative. There is instantaneous 1 : 1 frequency entrainment at all times; moreover, the attracting point moves within a bounded arc of the circle, and thus we have 1 : 1 phase-locking between θ and θ 1 .
In region III, the attracting point exists at some times but not other times. We refer to the epochs during which the attracting point exists as stable epochs; the remain epochs are epochs of neutrally stable dynamics. During the stable epochs, solutions from throughout the circle are attracted to the attracting point. While following the attracting point, these solutions pick up a negative contribution to the Lyapunov exponent, due to the gradient of the instantaneous vector field being itself negative at the attracting point; and then during each of the epochs of neutral stability, the solutions receive zero net contribution to the Lyapunov exponent, meaning that overall, as in [26], long-term LE are negative and the solutions remain synchronised with each other over all time. There is instantaneous 1 : 1 frequency entrainment during the stable epochs but no instantaneous frequency entrainment during the epochs of neutral stability; we will refer to this phenomenon as intermittent synchronisation. Overall, we do not have phase-locking between θ and θ 1 . However, unlike in the case of fixed-frequency driving, synchrony of an array of identical copies of θ [represented as different simultaneous solutions of Eq. (1)] is achieved and endures (even through the epochs of neutral stability) in the absence of a phase-locking mechanism. In other words, there does not need to be a phase-locking mechanism in place in order for the driving θ 1 (t) to cause θ to lose all memory of its initial condition and follow a path determined by the evolution of θ 1 (t).
Let us mention that there will be some very small subregions of region III where in theory, if one waits long enough, θ will come close to the instantaneous repeller around the start of a stable epoch [38] and thus receive a positive contribution to the LE, such that the reasoning here and in [26] can eventually break down and the asymptotic LE (if f (·) is defined ad infinitum) could even be zero. However, this theoretical phenomenon is unlikely to manifest in practice, due to the precision of fine-tuning of f (·) required for the phenomenon to occur, combined with the unphysical length of time that one is likely to have to wait for the phenomenon to take place. Indeed, no such phenomenon is ever observed in our numerics.
So then, in analogy to the case of fixed-frequency driving, we define the Arnold tongue as being the union of region II and region III, that is, the total region where the long-term LE will be negative.
From Eq. (10), the role of the modulation amplitude k here is clear: as k increases from 0, regions I and II decrease in size (although still extending infinitely), being symmetrically pushed back by the appearance and growth of region III, such that overall, the Arnold tongue is enlarged. In other words, increasing modulation amplitude induces stability.

D. Numerical results
For our numerics in this section, we take θ 1 (t) = ω 1 (t− k ωm cos(ω m t)); so the frequency modulation f (ω m t) is a sine wave f (ω m t) = sin(ω m t). Nonetheless, the results presented are just as valid for more general aperiodic slow modulation, and are demonstrated numerically to be true for aperiodic modulations in Sec. II F. We set ω 1 = 4 and ω m = 0.05, except where stated otherwise, and we investigate the effect of the remaining free parameters γ, ω 0 , and k. We integrate the two-dimensional system (3) with r p = 1 and = 5, except that for Fig. 3 (showing synchronisation between solutions of (1) with different initial conditions) and Fig. 6 (showing long-term LE together with average frequency entrainment), we simply integrate (1). All Lyapunov exponents, both longterm and time-localised, are computed following Benettin's canonical algorithm [39,40]; for the time-localised LE, we use a moving average of the expansion coefficient.
In (3), the radial LE at the limit cycle is equal to −5; therefore, since the maximum LE is greater than −5 in all our numerical experiments, it follows that this maximum LE corresponds to the phase dynamics defined by (1). The same is also true of time-localised LE, at least after the first few moments needed for the trajectory to approach the limit cycle.
First, we investigate the long-term stability of the system, by means of the numerically computed maximum LE defined over a long time-window. Stability is indicated by a negative value for the LE. In Fig. 2, we see that there is an Arnold tongue (shades of blue) similar to that shown in Fig. 1(b). As illustrated in Fig. 3, solutions of Eq. (1) synchronise with each other when the parameters lie in the Arnold tongue, but not when the parameters do not lie in the Arnold tongue. As shown in Fig. 2, the Arnold tongue is enlarged as the amplitude k of the frequency modulation is increased.
In other words, stability is induced by varying the frequency of the forcing over time. Quantitatively, we observe that the width of the Arnold tongue grows linearly with k.
While Fig. 2 shows the long-time stability, region III can only be distinguished and understood from the point of view of time-localised stability. The dynamics of Eq. (1) is illustrated over time for the three regions in Figs. 4(a-f), by time-frequency representation and by time-localised LE-namely, maximum LE defined over the time-window [t, t + τ ] where τ is a fixed number. Here, we take τ = 0.1 s.
Region III is a region of intermittent synchronisation where trajectories alternate between epochs of timelocalised stability and epochs of time-localised neutral stability; indeed, as the time t evolves, the time-localised LE alternates between epochs where it is negative, and epochs where it oscillates with high frequency around an average value of zero, as is seen in Fig. 4(e). Averaging over the total time yields a negative LE, meaning overall stability on average, even though the short-term stability is time-varying. Region I is thus the only region with a long-term LE of zero, and this region decreases in size, which means that the region of stability increases.
The distinction between the three regions can be seen by looking at the time-frequency representation of a trajectory in each of these regions, as shown in Figs. 4(a-c). In all three cases, the changing frequency of the driver is reflected in the frequency content of the driven oscillator. In Fig. 4(a), representing region II, the driving frequency is the only frequency present, as the frequency of the driven oscillator is entrained by that of the driving at all times. In Fig. 4(c), representing region I, we also see the natural frequency of the driven oscillator (cream, representing the highest amplitude), though slightly modulated by the driving. The fact that these two frequency modes are distinct shows that the driven oscillator's frequency is not entrained by the driving at any time. In Fig. 4(b), representing region III, we see the maximumamplitude frequency mode overlapping the driving frequency at some times, but not at other times. The times of overlap are when the frequency of the driven oscillator is entrained by that of the driving, and the other times are when there is no frequency entrainment. Thus, in this region, we have intermittent frequency entrainment. Comparing (a, b, c) with (d, e, f) in Fig. 4, we can see that in all three regions, absence of frequency entrainment coincides with time-localised LE that oscillate about 0, while the occurrence of frequency entrainment coincides with time-localised LE that stay negative over a longer time-interval. Now when investigating numerically the time-evolution of the time-localised LE, as in Fig. 4(e) one can clearly distinguish between those time-intervals where the timelocalised LE oscillates with high frequency around zero, and those time-intervals of length much greater than the periods of these aforementioned high-frequency oscillations during which the time-localised LE remains negative; and hence, one can numerically distinguish between the three regions. The proportion P t of time taken up by time-intervals where the time-localised LE remains negative is plotted in Figs. 5(a, c), across different parameter values. As in Figs. 4(d-f), we expect P t = 1 in region II, 0 < P t < 1 in region III, and P t = 0 in region I. We also plot in Figs. 5(b, d) the analytically obtained proportion P t of time for which the instantaneous vector field has a stable equilibrium. This is given by where ∆ω ± = ω 0 − ω 1 (1 ∓ k). The close resemblance between (a, c) and (b, d) helps confirm the validity of the numerical approach to distinguishing between the three regions. Fig. 5(b) provides a quantitative picture for the qualitative skeleton shown in Fig. 1(b). The average frequency difference Ω ψ = ψ = θ − ω 1 (1), with initial conditions θ i (0) = i−1 5 .2π, subject to the same driving θ1(t). Here, γ = 2.5 and k = 0.4. In (a), ω0 = 4, and so the system is in region II according to Eq. (9); in (b), ω0 = 6, and so the system is in region III according to Eq. (10); in (c), ω0 = 9, and so the system is in region I according to Eq. (8). In each case, the upper plot shows the first 8 seconds of the sine of the five trajectories, while the lower plot shows the distance between θ 1 (t) and θ 2 (t) over about the first 630 seconds (more precisely, 5 cycles of the frequency modulation); in (a) and (b), the inner graph shows the same information on a logarithmic scale. In (a) and (b), the system lies within the Arnold tongue as described in Sec. II C, and the five trajectories are observed to synchronise and to remain in synchrony; by contrast, in (c), the system does not lie within the Arnold tongue, and no synchronisation is observed.  Fig. 6(b), for (g) region III, and (h) region I. Interestingly, in region I, the main observed frequency oscillates in anti-phase with those of the driving frequency (dashed). Note that we here only predict the main frequency, and not the higher harmonics observed in (b, c). is a measure of the "average frequency entrainment" of the system [5]. In the traditional autonomous case k = 0 where the driving frequency is constant, nullity of Ω ψ is equivalent to actual frequency entrainment. The quan- tity Ω ψ is shown in Fig. 6(b) for γ = 2.5 across different values of k. On Fig. 6(a), the corresponding curves for the long-term Lyapunov exponent are displayed. Curves for Ω ψ for k > 0 are extremely similar to the case of driving with bounded noise ξ(t),ψ = ∆ω+γ sin(ψ)+ξ(t) (see [5]). This will be explored in further detail in Sec. II E. The similarity is due to the fact that only averages are considered, and time is forgotten. However, investigation of finite-time dynamics reveals that in region III, the frequency difference alternates between epochs where it is zero and epochs where it is non-zero [as in Fig. 4(g)].
To obtain this, we calculate the main observed frequencẏ θ main /2π of trajectories using the slow modulation assumption, by takinġ where Ω ψ (t) is the Ω ψ -value associated with the autonomous differential equation d ds ψ(s) = ∆ω(t) + γ sin ψ(s); results are plotted in Figs. 4(g) and 4(h).

E. Comparison between non-autonomous and noisy systems
Experimental science essentially seeks to understand the underlying mechanics of a system that gives rise to the observed behaviour. Since the study of timehomogeneous dynamics (deterministic or noisy) is very well developed in comparison to the study of nonautonomous dynamics, there is a tendency to assume For k = 0, the region of phase stability (as given by λ1 < 0) and the region of permanent frequency entrainment (as given by Ω ψ = 0) coincide; but for k > 0, the regions do not coincide: as k is increased, the region with λ1 < 0 is increased while the region with Ω ψ = 0 is decreased. The graphs look very similar to those in the case of harmonic driving with bounded noise [5]; therefore, investigation of time-variable dynamics for non-autonomously driven oscillators is necessary for an accurate understanding of the dynamical nature of the system. that for modelling purposes, the dynamics of a real-world system may be treated as statistically time-homogeneous. In this section we will illustrate, using our aboveidentified phenomenon of intermittent time-localised stability, how such a tendency may lead to the complete misidentification of some key aspect of the internal mechanics of a system.
There are various methods for analysing experimentally obtained time-series that are based on timeaveraged properties of the time-series, such as power spectra. The theory of both deterministic autonomous dynamical systems and autonomous systems perturbed by stationary noise is well developed, and in particular it is well-known that adding noise to a system can create stability which was not present in the absence of noise, e.g. [5,36,37,41,42]. Indeed, a one-dimensional phase oscillator model will almost invariably exhibit asymptotic stability of solutions when driven by stationary white noise [43]. Therefore, when seeking to understand the mechanism by which a real-world system behaves robustly against unpredictable external perturbations, if one observes in a time-series of measurements from the system a power spectrum similar to that of some noisy model, and if moreover this noisy model is known to ex-hibit stability with negative Lyapunov exponents as a consequence of the noise, then naturally one may come to the conclusion that the real-world system under investigation is subject to a significant level of noise and that this noise plays the key role in causing stability.
However, our results for the deterministic system (1)-(2) demonstrate that such a conclusion may be profoundly erroneous. The frequency modulation in (2) may be an entirely deterministic process that is not subject to any significant levels of noise. This gives rise to the deterministic non-autonomous equation (5), and we will illustrate that the time-averaged properties of (5) [with f (·) = sin(·)] are very similar to those of a noisy counterpartψ where ξ(t) is bounded noise. Physically, Eq. (13) represents the phase difference under a model in which the driving frequency modulation is assumed to be noisy. The similarity that we shall illustrate between the timeaveraged properties of (5) and (13) proves an important point: Since real-world systems are open and therefore subject to time-variability, one must examine temporally evolving dynamical properties of a system rather than just time-averaged properties, in order to account for the possibility that the mechanisms behind features of the observed behaviour are due to non-autonomicity. In the case of the system (1)-(2), in region III, the mechanism behind stability is not stationary noise but deterministic intermittent frequency entrainment between driving and driven oscillators, arising from the slow variation of the driving frequency. For simulations, dichotomous Markov noise ξ(t), which switches between ±D at rate µ (with ξ(0) = +D), was used [44]. Nonetheless, it is expected that any bounded noise will show similar behaviour to that presented here [5]; moreover, although the asymptotic properties of unbounded noise models exhibit a slightly different behaviour [5], any noise will effectively serve as bounded noise over typical physically relevant finite time-scales.
On average, the noisy and the non-autonomous systems will have properties as illustrated in Fig. 7 that look essentially the same. Indeed, both can be made to have an increased region for negative LE [see Fig. 7(a)] and a smaller plateau for average frequency entrainment [see Fig. 7(b)], as compared to the autonomous case given by k = 0 or D = 0.
The noisy and the non-autonomous systems can, however, be distinguised based on their dynamics over time. This is illustrated in Fig. 8 by trajectories and their timefrequency representation. In the non-autonomous case, one can see the regularly intermittent frequency entrainment between driving and driven phases in Fig. 8(c), where frequency entrainment corresponds to those times where the instantaneous power-frequency spectrum has only a single peak, and in Fig. 8(b), where frequency entrainment corresponds to the regular plateaus in the phase difference. By contrast, in the noisy case,  4) and (5) respectively were integrated, with ψ(0) = 0; for the noisy case, one sample path of ξ(t) was generated, and Eq. (13) was integrated with ψ(0) = 0, using the same sample path ξ(t) for all ω0 values. Parameters are set to ω1 = 4 and γ = 1; for the non-autonomous case, f = sin(·), k = 0.1 and ωm = 0.05, and for the noisy case D = 1.6 and µ = 10. (a) Lyapunov exponent λ1, and (b) average frequency difference Ω ψ . The nonautonomous and noisy cases, observed on average, present the same enlarging of the negative Lyapunov exponent region, and their Ω ψ is almost exactly identical, including the plateau.
the instantaneous power-frequency spectrum shown in Fig. 8(e) is significantly more bumpy around a peak that stays roughly fixed over time, and the phase difference in Fig. 8(b) looks like it is essentially drifting at all times. Despite these starkly visible differences in timevariable properties, the average power spectra as shown in Fig. 8(d, f) are reasonably similar to each other.
As expected, and shown in Fig. 9, the enlargement of the Arnold tongue holds. Moreover, more quantitatively, the results shown in Fig. 9(b, c) are almost identical to those shown in Fig. 2(c, d) respectively. This is because f (·) oscillates throughout the interval [−1, 1] in both cases, and therefore Eqs. (1) with θ1(t) = ω1(t − k ωm cos(ωmt)) in the non-autonomous case, and θ1(t) = ω1t − t 0 ξ(s) ds in the noisy case, and θ(0) = θ1(0) in both cases. Parameters are set to ω1 = 4, γ = 1, ω0 = 3; for the non-autonomous case k = 0.1 and ωm = 0.05, and for the noisy case D = 1, µ = 4. First ψ(t) is numerically obtained by integrating Eq. (5) or (13) as appropriate, with ψ(0) = 0, and then θ(t) is obtained by θ(t) = ψ(t) + θ1(t). (a) Sine of the driven phase θ, in the non-autonomous setting. (b) Phase difference ψ between the driving and driven oscillators, over time. The non-autonomous case presents epochs of phase-locking as seen by the regular plateaus, whereas the noisy phases' difference drifts without ever phase-locking to the driving, with an average velocity that is close to that of the non-autonomous curve. (c, e) Time-frequency representation (showing power, i.e. square of the amplitude) for θ extracted using continuous Morlet wavelet transform (with p = 1) with central frequency 3, and (d, e) the associated time-averaged power. The main difference is the presence of intermittent frequency entrainment in the non-autonomous case. The average-power spectra are very similar, and do not clearly distinguish the two cases. Long-term LEs were also found to be negative in both cases: for the non-autonomous Eq. (5) with ψ(0) = 0, the LE over 500 s was about −0.32, and for the noisy Eq. (13) with ψ(0) = 0, the LE over 500 s was about −0.24.

III. HIGHER-DIMENSIONAL CASES
In the above section, we showed analytically, and confirmed numerically, that enlargement of the stability region will always occur in the simple one-dimensional case. Nonetheless, the phenomenon of stabilisation by slow variation of the driving frequency, and the phenomenon of intermittent synchronisation under such variation of the driving frequency, may be found in a broader class of systems. To illustrate the more general scope of the stabilisation phenomenon, we illustrate it numerically in non-linear driven oscillators. We consider three cases: first, a typically forced van der Pol (vdP) oscillator; second, a vdP oscillator with the phase driven via diffusive coupling; and finally, a typically forced Duffing oscillator. All three cases are investigated with non-autonomous driving θ 1 (t) = ω 1 (t − k ωm cos(ω m t)), with ω m = 0.02. Long-term LE are computed over 10 cycles of the frequency modulation (about 3140 s).

A. Typically forced van der Pol
We consider a vdP oscillator that is directly forced by the external phase oscillator θ 1 (t), so that the vdP oscillator satisfies the differential equation For investigation of LE, we treat this as a first-order equation in (x, y)-space withẋ = y (and numerically integrate it as such). The long-term maximum LE is shown over parameter space in Fig. 10. The region with negative LE increases with the amplitude k of the frequency modulation, showing that the enlargement of the negative LE region still holds in this non-linear case.
Moreover, consideration of time-localised LE, as shown in Fig. 11(b), suggests that for some parameter values we have intermittency between epochs of stable dynamics and epochs of neutrally stable dynamics. Unlike in Fig. 4(e), the stable epochs are not characterised by negativity of LE defined over a very short sliding window, but rather over a suitably longer sliding window. This is because, if one were to freeze the driving frequencyθ 1 at any moment in time during such an epoch of stable dynamics, the solution of the resulting periodic differential equation (15) would not converge to a fixed point but most likely to a stable periodic orbit, for parts of which the vector field is locally contractive and other parts not, with contraction on average over each period. Hence, the intermittency is demonstrated most clearly by taking time-localised LE over a wider moving time-window, whose width is likely to incorporate several periods of the aforementioned stable periodic orbit. In Fig. 11(b), we see quite clearly (in red) the alternation between plateaus of zero time-localised LE and epochs where the time-localised LE dips to become negative. In the time-frequency representation shown in Fig. 11(a), during the epochs of zero time-localised LE, the power is shared mostly between two distinct peaks in the instantaneous spectrum, but during the epochs of negative time-localised LE, virtually all the power is concentrated around a single peak. As in Fig. 4(e), this suggests that instantaneous 1 : 1 frequency entrainment is taking place during the epochs of negative time-localised LE, but not the epochs of zero time-localised LE, and so overall the system exhibits intermittent synchronisation.

C. Forced and coupled Duffing oscillator
Here, we consider two Duffing oscillators, x and x 1 , unidirectionally diffuively coupled with strength g d so that x 1 drives x. Additionally, the driven Duffing oscillator is directly forced with external non-autonomous driving: We fix parameters δ = 0.3, β = 0.1, ω 1 = 1.2, ω m = 0.02, and coupling strength g d = 0.5. For investigation of LE, we treat the system as a first-order equation in (x, y, x 1 , y 1 )-space withẋ = y,ẋ 1 = y 1 (and numerically integrate it as such). The long-term maximum LE is shown over parameter space in Fig. 13. As k is increased, the chaotic region essentially decreases, giving way to either stability or neutral stability. The stability region does not strictly increase, as parts of the stability region become neutrally stable as k is increased; nonetheless, the phenomenon is still observed that for various fixed values of all the parameters other than k, increasing k has the effect of turning chaos into stability. From a control point of view, if the region of interest in parameter space is the chaotic region, one can stabilise the dynamics by adding time-variation to the forcing frequency.
FIG. 13. Numerically obtained long-term maximum Lyapunov exponent λ1 over parameter space for the coupled and forced Duffing oscillator (17), for different amplitude of frequency modulation k. In each case, 20 random initial conditions (x(0),ẋ(0), x1(0),ẋ1(0)) were taken from [−1, 1] 4 , and the average maximum LE over these trajectories is plotted. The region of chaotic behaviour (red) is reduced as many points are turned stable (shades of blue) as k is increased. Grey represents zero values.

IV. DISCUSSION
The work was motivated by real systems that exhibit dynamics with time-varying frequencies and are stable against external perturbation [1][2][3][17][18][19][20][21]. Surprisingly, not much analytical work has been carried out on such systems, and most of the work that has been carried out has used noisy driving [7,41,42] as the foundation of the model, or noise consisting of impulses at random times [5,37]. In these studies, it was shown that noise can create and increase stability. Outside of a stochastic approach, the only other way to incorporate timevariability is to model the system as a deterministic nonautonomous dynamical system. However, not much analytic theory of non-autonomous systems has been developed yet. The problem is additionally complicated by the fact that an asymptotic approach does not give the full picture as the evolving dynamics over shorter time-scales is missed. As illustrated in this work, changes in dynamical behaviour over shorter time-scales are of crucial importance in the types of systems considered. For if they had not been considered in this work, the phenomenon of intermittent synchronisation would have been missed. The work in this paper has provided a key new insight into systems subject to time-varying influences, by identifying the phenomenon of intermittent synchronisation and the region in parameter space where it occurs, and thereby showing the enlargement of the Arnold tongue. This insight also has potential for being the foundation of future methods to induce stability in complex or other systems; the fact that non-autonomous driving allows for average stability to be achieved without the need to maintain frequency entrainment at all times may be of significant advantage.
The basic adiabatic reasoning underlying our analytical approach is the same as that employed and investigated in [26]. It is this reasoning that has led us to our new discovery that increasing time-variability inherently induces stability in phase oscillators. We have also employed numerical tools to visualise time-localised dynamics as derived by this adiabatic reasoning, namely time-localised LEs as in Figs. 4(d-f) and time-frequency representation as in Figs. 4(a-c). We hypothesise that for a time-frequency representation applied to experimental data, a result resembling Fig. 4(b) could be a signature of intermittent synchronisation. We also investigated the slow variation assumption in higher-dimensional systems, and numerically illustrated the creation of stability as the amplitude of variation is increased, as well as the occurrence of intermittent synchronisation. In this way, we showed that the phenomena of stabilisation and intermittent synchronisation under slow variation of the driving frequency occur more broadly than just in the case of phase oscillators.
Chronotaxic systems have been introduced to model the distinctive feature of real-life oscillatory systems, that they are able to keep their time-varying dynamics resistant to external perturbations [18,22,23]. Chronotaxic systems were defined in previous works by the necessary condition that a time-dependent attractor exists, and that trajectories in its close vicinity always move closer to it [18], or alternatively just by the existence of a positively invariant time-varying region in which the dynamics is always contracting [23]. However, it seems reasonable to expect that in real life, there exist stable oscillatory systems for which no trajectory is always instantaneously locally attractive. Instead, in contrast to currently existing definitions of chronotaxicity, an intermittent synchronisation phenomenon such as identified in this work would give rise to stability on average. Thus, we have broadened the definition of chronotaxic systems, increasing its potential for effectively modelling and understanding real-life systems, which are non-isolated and therefore continuously subjected to time-varying external influences.

V. SUMMARY
We have shown that driving a phase oscillator with an arbitrary slowly varying frequency always induces stability: the larger the amplitude of the frequency modulation, the larger the stabilty region. We have furthermore shown numerically that this phenomenon occurs in more complex cases where the driven oscillator is higher-dimensional and non-linear, hinting at the wider scope and importance of the effect at hand. We have even shown numerically that chaotic regions in parameter space can be made stable by the same mechanism. If only the quantities λ 1 and Ω ψ , which describe time-averaged properties of the system, are considered, the system looks the same as in the case of driving with noise. However, in reality, our fully deterministic example exhibits some time-localised frequency entrainment, whereas none is exhibited in the case of driving with bounded noise. It is therefore clear that the non-autonomous deterministic system could be misinterpreted as a noisy one if only time-averaged quantities are considered.
The enlarged stability region makes time-variable driving very suitable for real-world modelling and for engineering, where a controlled adjustment of the frequencies is often of key importance. We believe that this type of model will find applications in many fields, including physics, biology, medicine, and climate dynamics. ACKNOWLEDGMENT We thank Arkady Pikovsky and Martin Rasmussen for valuable discussions, as well as Phil Clemson, Bastian Pietras and Tomislav Stankovski, for useful comments on the manuscript. We especially thank Yevhen Suprunenko for his invaluable help and for the useful discussions that we had. This work has been funded by the EU's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 642563, and the EPSRC grant EP/M006298/1 A device to detect and measure the progression of dementia by quantifying the interactions between neuronal and cardiovascular oscillations.