Transient chaos in time-delayed systems subjected to parameter drift

External and internal factors may cause a system's parameter to vary with time before it stabilizes. This drift induces a regime shift when the parameter crosses a bifurcation. Here, we study the case of an infinite dimensional system: a time-delayed oscillator whose time delay varies at a small but non-negligible rate. Our research shows that due to this parameter drift, trajectories from a chaotic attractor tip to other states with a certain probability. This causes the appearance of the phenomenon of transient chaos. By using an ensemble approach, we find a gamma distribution of transient lifetimes, unlike in other non-delayed systems where normal distributions have been found to govern the process. Furthermore, we analyze how the parameter change rate influences the tipping probability, and we derive a scaling law relating the parameter value for which the tipping takes place and the lifetime of the transient chaos with the parameter change rate.


Introduction
A vast variety of systems are not static over time, in the sense that they can not be modeled by the same set of equations and parameters as time passes. External factors as the increase in greenhouse gases in the context of climate dynamics [1] are one of the possible causes. Also, this change can be caused by internal factors, i.e., by the nature of the system itself as in the case of some engineering systems that wear-out. Furthermore, some systems may be modified manually by a experimentalist that, for instance, injects chemical substances in a reactor at a slowly varying flux [2]. Mathematically, these processes are reflected by a drift in one of the parameters of the model.
For systems suffering a drift at slow rates, the analysis of the frozen-in system, that is, the model without time dependence in the parameters, is of great importance. A vital part of this analysis consists on the computation of bifurcation diagrams that depict the dynamics of the frozen-in system for different values of a parameter. It has been shown in [3,4] that there is a delay in the regime shift when a parameter drifts, at small but non-negligible rates, and crosses a bifurcation point (p bif ). This transition occurs at a value of the parameter p cr > p bif . This is usually referred to as the delay effect in systems with parameter drift. It was also shown that the value of p cr depends on the parameter change rate. This type of bifurcations that are crossed due to the time dependence of one of the parameters are called dynamic bifurcations [5].
Previous work on dynamic bifurcations has been initially restrained to systems where only regular attractors were involved. More recently, some attention has been focused on chaotic attractors by using either maps [6,7] or flows [8]. In both cases, an ensemble approach is needed and a normal distribution was found for the values of p cr . Here, we aim to broaden systems currently explored, by studying dynamic bifurcations in an infinite dimensional system: a time delayed oscillator. As far as we know, this type of systems have not received much attention in the context of dynamic bifurcations and they are worth a study due to the ubiquity in Nature of time-delayed systems. Physically, the time delay accounts for the finite propagation time of information. In the case of laser arrays, for instance, this is due to the finite speed of light [9] and the time delay determines its stability [10]. Here, we consider that the time delay slowly drifts with time and we explore its implications. This drift may account, for instance, for the seasonal fluctuations to the regeneration time of a resource [11].
Besides the delay effect, the parameter drift provokes another interesting phenomenon: bifurcation-induced tipping [12]. When the system presents multistability trajectories, for p > p cr , have to tip to any of the attractors that are stable with a certain probability. This tipping probability also depends on the parameter change rate. Finally, in a time framework, we may consider that after the tipping, the drifting system is in its steady state. Before that, it is in its transient state. Depending on the parameter change rate, the transient state would last for a longer or shorter period of time. This may be a problem in some engineering systems as a parameter drift would be identified long after its start, possibly causing parameter drift failure [13,14]. This article is organized as follows. In section 2 we present our system: the Duffing oscillator with time delay, and we analyze its dynamics when the parameters are frozenin. In section 3 we introduce the time dependence in the time delay parameter. We explore the tipping probabilities in section 4 and the delay effect in section 5, deriving a scaling law that relates the delay and the parameter change rate. Finally, the time framework is analyzed in section 6, where the arising transient chaos phenomenon is characterized. Discussions and conclusions are drawn at the end.

Time-delayed Duffing oscillator
As a paradigmatic example of a nonlinear oscillator with delay, we study the transient dynamics of the undamped and unforced Duffing oscillator with a delay term of the form γx(t − τ ), where τ is the time delay and γ is the amplitude of the delay. The equation for the system reads: We consider the system in absence of dissipation since we focus our study in the effect of the time delay term. This can be performed experimentally if we consider that the experiments are carried out in the laboratory under very low pressure and therefore in the situation of vacuum. In this physical context, thanks to the time delay and the absence of damping, the oscillator presents a region of parameters for which a chaotic attractor is stable. We choose the same parameter values considered in [15], that is, α = −1, β = 0.1 and γ = −0.3. With these values the potential has the shape of a double well, with an unstable fixed point at the origin and two stable fixed points at the bottom of the wells. Taking x(t) = x(t − τ ) = x * , we find the numerical values for the fixed points as In this section, we present the system dynamics for the frozen-in case, that is, for the system when the time delay does not vary with time. For different but fixed in time values of τ , the system's behavior changes. To uncover this, we explore the attractors that appear for two values of τ that correspond to the starting and ending scenarios of the ramp τ (t) that is explored in the next section.
In the context of time-delayed systems, initial conditions are replaced by history functions: u 0 (t) = (x 0 ,ẋ 0 ), which are a set of initial conditions in the continuous time interval [−τ, 0]. The history functions that we use here are, for simplicity, constant functions. Due to the multistability of the system, different history functions lead to different attractors. In figure 1, we show all the coexisting attractors in phase space for τ = 3 and τ = 4.
In figure 1 (a), we observe two attractors: a limit cycle, L, and a chaotic attractor, A. The bifurcation diagram for τ ∈ (0, 5] and history functions u 0 = 1 and u 0 = −1, which was calculated in [16], showed that the chaotic attractor disappeared at τ = 3.6. Consistently, for τ = 4, in figure 1 (b), it can be seen that the chaotic attractor is no longer stable and that there are two new attractors, corresponding to the previously calculated fixed points, x + and x − from equation (2). Although the limit cycle is still stable in both scenarios, it is important to remark that its amplitude decreases with τ .  (1) for every history function with the algorithm Tsit5 from Julia [17] for delay differential equations, which is a Runge-Kutta method with adaptative stepsize-control. The parameters for the numerical calculation have to be taken cautiously as the results have to be very precise to lead to the correct attractor, specially in the chaotic region. In our case, we took RelTol=10 −9 , AbsTol=10 −12 and maxiters=10 9 . The criterion for convergence to fixed points is that the trajectory enters a square of 0.01 × 0.01. In the case of the limit cycle, the criterion for convergence is that the peaks from the time series become equally spaced, thus periodic. When none of these criteria are met, the chaotic attractor is assigned.
The basin of attraction for τ = 3 is depicted in 2 (a). In blue, the history functions that lead to the limit cycle, L, and in pink the history functions that lead to the chaotic attractor, A. The fixed points x + /x − are also represented as black crosses. As it can be seen, there is a bone-shaped region that corresponds to the chaotic attractor basin. Outside this bone-shaped structure, all the trajectories approach the limit cycle.
Following a similar sketch, the basin of attraction for τ = 4 is depicted in 2 (b). The yellow and green regions correspond to the set of history functions that lead to x + and x − respectively. As it can be seen, around each fixed point there is a circular region Figure 2. Basins of attraction before and after the bifurcation at τ = 3.6. Two black crosses point the x ± attractors. (a) For τ = 3, the pink and bone-shaped region corresponds to the set of history functions that go the the chaotic attractor. On blue, the ones that go to the limit cycle. (b) For τ = 4, the bone shaped structure is replaced by the almost circular basins of x + in yellow and x − in green around the attractors, which are also intermingled with the blue basin maintaining the bone-shaped structure.
of points with smooth boundaries that leads to them, in other words, trajectories that start near x ± = ±( √ 13, 0) end up in the yellow/green attractor. Also, there are other history functions, mainly inside the previous bone-shaped structure, that approach x ± .

Time dependent delay, τ (t)
Once the behavior of the system for fixed values of the time delay (before and after the bifurcation at τ = 3.6) is known, we let this parameter slowly evolve with time between those scenarios. The main difference with the previous section is that now we are interested in the case of the parameter variation during the evolution of the system. In other words, the parameter turns into a slowly varying function of time of the form τ (t) = τ 0 + εt, where ε is a sufficiently small parameter compared to the natural time scale of the system. Dynamic bifurcations, as mentioned before, are the bifurcations which are crossed due to the time dependence of this parameter [5]. Here, we aim to study this phenomenon for systems where not only chaotic attractors are involved, but time delay is present too. This is an interesting area as models that include a time delay cover the unavoidable phenomenon of the finite propagation of information. We let not just a normal multiplicative or additive parameter vary with time, but the time delay itself.
The first difference in this analysis from the one of regular attractors is that single trajectories are no longer representative and do not contain all the possible dynamics. Thus, we follow an ensemble of trajectories starting on the same set of history functions used in the previous section to calculate the basins of attraction.
Furthermore, we ought to redefine our system including the time delay dependence with time. We replace τ in equation (1) where τ 0 is the initial value of τ , t 1 is the time for which the parameter shift starts and t 2 when it ends. We refer to this system as the non-autonomous system, in contrast with the frozen-in system. For simulation purposes, we take τ 0 = 3 and we let t 1 = 100. This way, we let the system evolve to its steady state (remember figure 2 (a)), before the shift in τ starts. The orbits evolve as in the previous section until t = t 1 , when τ starts shifting and finally they reach one of the attractors that are stable after the bifurcation (see figure  2 (b)). For t 2 , we choose sufficient and different values depending on ε so that all the trajectories have the time to reach one of the attractors.

Random tipping
Using the same terminology as in [18], we refer to the basin for the non-autonomous system as scenario-dependent basin, because it depends on the parameters of the parameter shift. Figure 3 depicts the scenario-dependent basin of attraction for ε = 10 −2 . In this case, t 2 = 350. Figure 3. State dependent basin of attraction for ε = 10 −2 . As the chaotic attractor looses stability at τ = 3.6, the trajectories tip to one the attractors that are stable (L, x ± ). This leads to phenomenon of random tipping, and as result we obtain a riddle basin inside the bone-shaped structure. We distinguish three different basins corresponding to the attractors past the bifurcation, see figure 1 (b). Outside the bone-shaped structure, all initial functions approach the limit cycle. If we compare this with figure 2, we observe that the trajectories initially in the limit cycle never jump to a different attractor, not even for faster parameter change rates. Inside the bone-shaped structure, the yellow/green smooth boundary regions present for the frozen-in system disappear now, even if we stop the shift at τ = 4. Due to the presence of the chaotic attractor, there is no pattern, and the history functions that lead to each basin are completely intermingled. This is because the final destination of each trajectory depends on the precise moment that it is caught wandering in the chaotic attractor when it tips. We may say that the chaotic attractor acts as a memory-loss agent. In other words, predictability of individual trajectories is lost because the passage through the chaotic attractor induces fractal basins of attraction. Similar phenomena [18] for non time-delayed systems has been addressed as a random tipping.
We may repeat the same procedure for different values of the parameter rate of change, ε. In table 1 we show how the number of trajectories that approach each attractor varies with ε. This is also referred to as tipping probability as it reflects the probability of a random history function to end in each of the attractors.
Due to the symmetry of the potential well, the percentages for attractors x + and x − are the same. For faster changes in τ , the probability that a trajectory ends up in the fixed points x ± rather than in the limit cycle, L grows. For the sake of clarity, this is also represented in figure 4, where the blue points correspond to the percentage of trajectories that leave the chaotic attractor, A, and end up in the limit cycle, L. In the same way, the yellow/green points correspond to the percentage of trajectories that leave the chaotic attractor, A, and end up in x ± .

Scaling law
Although we have uncovered the scenario-dependent basins of attraction, one question is still lacking. When does the orbit jump from the chaotic attractor to either of the attractors that are stable at the end of the parameter shift? From [16], we know that the bifurcation through which the chaotic attractor looses stability occurs at τ = 3.6 . Dependence on the parameter change rate for the percentage of trajectories that leave the chaotic attractor and approach the limit cycle (blue points) or the fixed points, x ± (yellow/green points). For higher parameter change rates, more trajectories end up in the limit cycle.
for the frozen-in system. However, the picture here is different as our parameter varies during the evolution of the system. Now, we explore if the value of τ for this transition, which we call τ cr , is different from the one in the frozen-in case.
In figure 5, we show two time series corresponding to two different history functions for the non-autonomous system with ε = 10 −2 . The secondary x-axis marks the evolution of the parameter. As it can be seen, in both cases, the behavior at first is chaotic. At t = 100, τ starts increasing. At τ = 3.6 (or t = 160), the chaotic attractor lost stability for the frozen-in case. However, the chaotic behavior is still present for a while after until the system jumps to x − in (a) or L in (b). The vertical dotted lines indicate this moment and the value of τ cr . Between τ = 3.6 and τ cr , the chaotic attractor is a metastable state of the system.
For numerical purposes, we calculated τ cr as the value of τ for which the orbit is at a distance of dx = 0.001 from x ± , or the value of τ for which the period of oscillation stabilizes and the amplitude starts decreasing.
As previously mentioned, due to the chaotic attractor, this analysis has to be performed for a large set of history functions. This is why for each value of the parameter change rate, we do not obtain only one value of τ cr ; instead we obtain a distribution of values. In figure 6 we can see the form of this distribution for the cases A → x ± and A → L when ε = 7 · 10 −3 . Both of them are fitted to a gamma distribution with positive skew. It is interesting to notice that the variance is more pronounced for the trajectories going to the limit cycle, while for the ones going to x ± , the deviation in the values for the transition is smaller. Inside both histograms, we have indicated the median value of τ cr , that we denote asτ cr , since it is more representative in this case than the mean The secondary x-axis shows the time dependence in the τ parameter. It can be seen that the transient lasts for a long period of time and that the transition to the steady state starts at a value of τ past the bifurcation value for the frozen-in system, this is, τ = 3.6. This value (τ cr ) is marked with a vertical dotted line.
values. We conclude that history functions that end up in the limit cycle, on average, leave the chaotic attractor before and for a wider range of values of τ . Figure 6. Gamma distribution of τ cr for ε = 7 ·10 −3 for the trajectories that leave the chaotic attractor and approach (a) the fixed points or (b) the limit cycle. The variance is higher in the later case, but the median value is smaller. Thus, history functions that end up in the limit cycle, on average leave the chaotic attractor before and for a wider range of values of τ .
This type of distribution is in contrast with other analysis performed for non timedelayed systems. In [6][7][8] it was found that the values of the parameter for the transition follow a normal distribution instead of a gamma distribution. This might be one of the effects that delay causes in the system. It implies that there are some trajectories that tip at much larger values than the average.
If we repeat the same process for different values of ε, we may deduce how it affects the value of τ for the transition. Figure 7 depicts the dependence ofτ cr with ε for the trajectories that end up in the fixed points (yellow/green) and the limit cycle (blue). As we can see, for higher rates the delay phenomenon is more severe in both cases. However, we can see that on average the transition occurs before for the trajectories heading to L rather than x ± . Figure 7. Scaling law for the median valueτ cr with the parameter change rate. The points (in blue for trajectories that approach the limit cycle and in yellow/green the trajectories that approach the fixed points) correspond to numerically calculated values and in red we show a power law fit.
The points in figure 7 are fitted to a power law that constitutes the numerical scaling law for the median valueτ cr A → x ± :τ cr = a · ε 4/15 + τ 0 (4) where a, b > 0 are constants and τ 0 = 3.6 in our case, with a R-square: R 2 = 0.9993 and R 2 = 0.9999 respectively. We have added the point for the limit ε → 0 which corresponds to the frozen-in case for which the chaotic attractor looses stability at τ = 3.6. This law indicates that an increasing parameter change rate increases the parameter value for the transition. However, the increase inτ cr is reduced for higher values of ε as the slope of the curves reduces with the parameter change rate.

Transient chaos interpretation
In this section we consider again our system with parameter drift and the previous results from a different perspective. Until now, we have considered that the drift in τ represented a small perturbation to the associated frozen-in system and that the chaotic behavior after τ = 3.6 was a metastable state that ended at τ cr . Also, we calculated the value ofτ cr for which the delayed bifurcation takes place. Now, we study the system without any previous knowledge of the behavior of the frozen-in system. From a experimentalist point of view, sometimes it is the regime shift which indicates that there is a parameter drift and not the other way around. However, if this regime shift occurs later than expected due to the delay effect, the parameter may have reached dangerous values when it is noticed by the experimentalist [13,14]. This is why, in this section, we are interested in measuring the time that the system takes before the tipping.
In this time framework, we deal with a non-autonomous system that behaves chaotically for a finite time before reaching one of the attractors (see figure 5). In this sense, the system presents transient chaos and the scaling law predicts the end of the transient state.
In order to characterize the transient chaos regime, we analyze the decay with time in the number of trajectories that still present a chaotic behavior. In figure 8 we represent N(t) as the normalized number of trajectories in the chaotic attractor for a time t. The blue lines correspond to the trajectories that end in the limit cycle and the yellow/green lines to the ones that end in the x ± attractors. This is calculated for ε = 5 · 10 −3 (dotted lines) and ε = 10 −2 (full lines). Note that the decay in τ starts at t = 100, and that τ and time t are equivalent through equation (3). When the curves decrease to zero, the transient chaos regime ends and every trajectory reaches its steady state.
As it can be seen, the decay with time slows down at the end of the curve. This is related to the gamma distribution for τ cr ( figure 6). In fact, we are representing nothing more than the complementary cumulative distribution function of figure 6 in terms of time. There are some trajectories with a transient lifetime higher than the average and this produces a delay in reaching the steady state for the full set of history functions. Similarly, the variance in figure 6 is reflected now as the time since the curve N(t) starts decreasing until it reaches zero, this is more pronounced for the set of trajectories going to the limit cycle (blue lines) while the others present a more step-like shape.
If we compare between different values of ε, we can assert by looking at figure 8 that the transient lifetime decreases with the change rate of the parameter shift. For faster changes in the parameter the steady state is reached faster, although it occurs for further values of τ from τ 0 .
Finally, we may translate the scaling law deduced in the previous section to a time framework just by using equation (3): A → L :t tr = d · ε −6/7 + t 1 , wheret tr refers to the median lifetime of the transient dynamics which follows the same gamma distribution as τ cr in figure 6. Also, t 1 in our case is 100 and c, d are positive constants. The numerically calculated points and the scaling laws are depicted in figure  9 for a visual guidance. Again, we obtain that the median value for the transient lifetime decreases with the parameter change rate following a power law. This decrease is more pronounced in the case of the trajectories that approach the fixed points as the slope is bigger in that case.
For any change rate, the transient lifetime of the trajectories heading towards the limit cycle is shorter. Of course, when ε → 0, the transient lifetime tends to infinity as in the limit when the parameter does not vary with time, the trajectories stay in the chaotic attractor forever. This is an interesting result as if a parameter is varying very slowly, a experimentalist may not notice that it is varying as the regime would not shift for a long time after the start of the drift. Specially, for the history functions corresponding to the tail of the lifetimes distribution. This may be dangerous in some engineering systems due to the parameter drift failure mechanism [13,14].

Conclusions
In the present work, we have analyzed the dynamics of a time-delayed oscillator whose time delay suffers a drift in time. The time delay increases linearly with time crossing a bifurcation value, after which the system presents multistability. We have found that trajectories initially in the chaotic attractor tip to one of the remaining attractors with a certain probability that depends on the parameter change rate. For faster rates more trajectories tip to the fixed points x ± instead of the limit cycle. However, predictability is lost in the non-autonomous system as the basins are riddle-like. In this sense the chaotic attractor acts as a memory-loss agent.
Also, we have found that the delay effect in the regime shift is present in timedelayed systems. However, it acts in a different way from the previously reported in other type of systems. The distribution of values of τ for which trajectories tip, this is τ cr , follows a gamma distribution instead of a normal distribution. This implies that there is a set of trajectories that tip at larger values of τ .
Furthermore, we derived two scaling laws relating the median value of τ cr , that is, τ cr and the parameter change rate for the cases for which the trajectories tip to the limit cycle and for the ones that tip to the fixed points x ± . We have also shown that trajectories that end up in the limit cycle tip before than the ones that tip to x ± and that for faster parameter change rates,τ cr increases in both cases.
Finally, we have characterized the system behavior in terms of time. As the chaotic dynamics lasts for a finite amount of time, we may say that the system presents transient chaos and we derived a scaling law for the transient lifetime. We conclude that for very small parameter drifts, the transient may last for unexpected long times. This may cause problems in the context of engineering due to parameter drift failure mechanisms. On the other hand, for fixed parameter change rates, trajectories that end up in the limit cycle have a shorter transient regime on average.