Timing of Transients: Quantifying Reaching Times and Transient Behavior in Complex Systems

When quantifying the time spent in the transient of a complex dynamical system, the fundamental problem is that for a large class of systems the actual time for reaching an attractor is infinite. Common methods for dealing with this problem usually introduce three additional problems: non-invariance, physical interpretation, and discontinuities, calling for carefully designed methods for quantifying transients. In this article, we discuss how the aforementioned problems emerge and propose two novel metrics, Regularized Reaching Time ($T_{RR}$) and Area under Distance Curve (AUDIC), to solve them, capturing two complementary aspects of the transient dynamics. $T_{RR}$ quantifies the additional time (positive or negative) that a trajectory starting at a chosen initial condition needs to reach the attractor after a reference trajectory has already arrived there. A positive or negative value means that it arrives by this much earlier or later than the reference. Because $T_{RR}$ is an analysis of return times after shocks, it is a systematic approach to the concept of critical slowing down [1]; hence it is naturally an early-warning signal [2] for bifurcations when central statistics over distributions of initial conditions are used. AUDIC is the distance of the trajectory to the attractor integrated over time. Complementary to $T_{RR}$, it measures which trajectories are reluctant, i.e. stay away from the attractor for long, or eager to approach it right away. (... shortened for arxiv listing, full abstract in paper ...) New features in these models can be uncovered, including the surprising regularity of the Roessler system's basin of attraction even in the regime of a chaotic attractor. Additionally, we demonstrate the critical slowing down interpretation by presenting the metrics' sensitivity to prebifurcational change and thus how they act as early-warning signals.


Introduction
In complex dynamical systems, the importance of a trajectory's transient, i.e. the part of the trajectory "away" from the attractor, plays an important role in physics research, e.g. for lasers [7,8], the dynamical Ising model [9] and other parts of statistical physics [10][11][12] as well as in various other fields, including ecology [13,14], biology [15], economics [16], medicine [17] and climate change [1,3] with specific focus on long transients in [3,13,18]. Hastings [14] stresses the importance of different time scales and points out how the transient dynamics can be very different and much more interesting than the asymptotic behavior. In addition, he points out how saddles play a central role by inducing long transients.
In this article, we devise novel metrics that measure how long it takes to reach the system's attractor to foster the study of transients.
Even though common methods for that exist, they are confronted with four essential problems. (I) divergences: the attractor is not reached in finite time for a large class of physically relevant systems; (II) physical interpretation: when using -neighborhoods the results depend strongly on the choice of and similarly for other parametrized methods; (III) discontinuities: small changes in the parameter often have a large (noncontinuous) effect on the measured time; and (IV) non-invariance: the results depend on the choice of variables. Problem (IV) is particularly important, as a result should be a property of the dynamical system and thus independent of the choice of variables, i.e. invariant (or correctly transforming) under smooth transformations of the state space (cf. "smoothly equivalent" in [19]).
As these problems are fundamental and have not yet been solved, we do not aim to just improve a quantitative measure but actually provide new solutions to general properties and consistency requirements.
The two novel metrics proposed in this article, the Regularized Reaching Time (T RR ) and the Area under Distance Curve (audic), have been developed in order to treat the aforementioned problems. The first is defined by the difference of the reaching times with respect to a reference trajectory and thus actually measures a time. The idea is that even though the actual reaching times diverge (problem (I)), the difference converges and we can compare which trajectories reach the attractor earlier or later. The second one, audic, is the distance to the attractor integrated along the trajectory. This means that it takes a different point of view and measures which trajectories are reluctant, i.e. stay away from the attractor for long, or eager, i.e. approach it right away.
Both metrics are shown to be Lyapunov functions and this property of audic is used for the computation. In the outlook, we even suggest this property to be the basis of an improved definition of T RR .
When applying these metrics to the global carbon cycle model [3] and the chaotic Rössler oscillator [5], their potential as early-warning signals [1,2] becomes apparent. Statistics of their distributions in state space represent the system's critical slowing down (CSD) [1,2,20] after a shock, i.e. an instantaneous and noninfinitesimal perturbation, uncovering prebifurcational changes in the transient behavior. In contrast, CSD is usually done with (local) noise only. The usage of shocks has been developed in the context of Basin Stability [21,22] and its extensions [18,[23][24][25].
With this new approach, we have been able to uncover new features of the systems: the basin of attraction in the chaotic Rössler system is unexpectedly regular and the basin separation in the carbon cycle is due to the strong stable manifold acting as a separatrix induced by a saddle, demonstrating the idea how saddles can lead to long transients. Additionally, we show how the metrics work well as early-warning signals by detecting the prebifurcational changes.
The remainder of this article is structured as follows.
After stating the fundamental problems of reaching time definitions in Sec. 2, we present two complementary solutions in Sec. 3 and apply them to examples in Sec. 4 before concluding with a discussion and some outlook in Sec. 5. The Appendix contains more technical comments, calculations and additional information. Moreover, Appendix I contains mathematical definitions and proofs, putting the ideas presented in the main text on solid footing.

The problems of reaching time definitions
A basic property of a large class of complex systems is that trajectories reach the attractor in infinite time only, even for steady states or limit cycles and generally most systems of ordinary differential equations with smooth right hand side functions. Figure 1: (color online) The distance (x-axis) over time (y-axis) for an example system with a stable, spiraling fixed point and a saddle, chosen to show the occurrence of problems (I)-(IV) for commondependent definitions of reaching times, t 1 (x 0 ) and t L (x 0 ), as discussed in the text.
Two common metrics that work around this simply measure when the trajectory starting at x 0 enters an -neighborhood around the attractor the first time (t F (x 0 )) or the last time (t L (x 0 )), i.e. when the neighborhood is not left anymore. In Fig. 1, the (Euclidean) distance (x-axis, dotted light blue line) vs. the time t (y-axis) of a basic example system's trajectory have been plotted (more details and explanations in Appendix A); t F (x 0 ) and t L (x 0 ) have been added.
Firstly, we observe the divergence of t for → 0, corresponding to the actual time until reaching the attractor (problem (I)).
Because t diverges, t F (x 0 ) and t L (x 0 ) depend heavily on the choice of so that a proper physical interpretation is rather difficult (problem (II)). It is far from being obvious what "close" or "when the transient is over" means. Thirdly, the strong discontinuities for t F (x 0 ) and t L (x 0 ) in Fig. 1 make the choice of a proper even harder (problem (III)). Problem (IV) from the introduction is noninvariance: Using a different set of variables, i.e. smoothly transforming the system, gives different values for t F (x 0 ) and t L (x 0 ), because the Euclidean distance is not invariant. This means the results are not a genuine property of the dynamical system but just of the representation.
Other metrics are based on characteristic times [26] and Lyapunov exponents [27]. Though common, the former suffer the same problems as the approaches above and are constant for a 1-dimensional linear system which is counter-intuitive when thinking about reaching times. The latter share these problems but are invariant under changes of variables, i.e. they are physical in that sense. Note that Lyapunov exponents do not capture the transient at all but are an asymptotic feature of the system only.
An extended discussion of these problems including an exemplary model is given in Appendix A.

Two novel, complementary solutions
To overcome the above problems, we introduce two metrics: Area under Distance Curve (abbreviated as audic, D) and Regularized Reaching Time (T RR ), and show that they naturally lead to a transient analysis from separate points of view.
(i) Area under Distance Curve (D) comes from the idea that a trajectory stays far away from the attractor during the transient while it is close in the asymptotics. A distance function d(·, ·) is needed to have notions of "far" and "close" and we define audic as where A is the attractor with the basin B A and x(t) the trajectory. Hence, we look at the cumulative distance and remove the influence of the asymptotics. Note the strong difference of D to t F (x 0 ) and t L (x 0 ). Both of the latter are very sensitive to small changes of the distance function around the attractor and to ; and it is difficult to even find a sensible notion of distance. In the audic measure D instead, choosing a tailor-made distance function d(·, ·) allows to adapt the measure to specific research questions, e.g., by letting d(x, A) represent some form of costs or damages due to being away from the attractor. This approach solves problem (IV), too, because the distance function is transformed correctly.
Initial conditions with high values of audic are called reluctant and low values eager. This terminology is used to emphasize that reluctant states go through large transients far away from the attractor, while eager states approach it "right away".
In the case where A is a hyperbolic fixed point of an ODEẋ = f (x), we show under certain mild conditions on d(·, ·) that audic is a Lyapunov function [6] uniquely defined by D(A) = 0 and d dt is the orbital derivative. Thus, the level sets of D foliate the basin of attraction and are forward invariant under the flow, which will be used for the definition of T RR .
Proofs and further properties are given in Appendix D.1, Appendix G and Appendix I.1, incl. a convergence discussion of Eq. (1).
(ii) Regularized Reaching Time (T RR ) is based on time differences between trajectories. It can be interpreted as the additional time (positive or negative) that a trajectory starting at a point of interest x 0 needs to reach the attractor after a reference trajectory has already arrived there. A positive or negative value means that the trajectory at hand arrives by this much later or earlier, respectively, at the attractor than the reference trajectory does. Since the actual reaching times are both infinite, T RR is formally defined as the limit for −→ 0 of the difference between how long the trajectory at hand and the reference trajectory need to reach the audic level set with value .
This idea is put in equations as follows. First we define the time t audic (x 0 , ) := t where = D (x (t)) and x(t) = x 0 . Note that the forward invariance of audic provides uniqueness of t audic . Next, we choose the initial condition x ref of the reference trajectory. And finally, we define the Regularized Reaching Time by taking the limit While the existence of this limit is far from obvious, we can prove it under mild conditions for an important class of systems: namely hyperbolic fixed points (see Appendix I.2).
Problem (IV), noninvariance, is circumvented because for any smooth invertible coordinate transformation Φ the equation where T RR is the regularized return time T RR in the system with changed variables.
Our numerical results indicate that this idea is sensible for more complex attractors, too, particularly the limit cycle discussed in Sec. 4.3 and the chaotic Rössler system below.
T RR represents the actual time by how much a trajectory reaches the attractor later or earlier than the one starting at the reference point. Note that we used the audic level sets in the definition for T RR as they are parametrized, bounded, forward invariant foliations. Foliations like this are usually hard find or compute but for these particular ones an efficient algorithm was developed (see Appendix F) and their computation can be done even for chaotic systems. Also, the usage of audic avoids local geometric measures that can easily induce problem (IV) (non-invariance).
We find also that T RR is a Lyapunov function, this time with constant (negative) orbital derivative −1; and in Appendix I.2 we show that for hyperbolic fixed-point attractors, T RR is a time-parametrization of the strong stable foliation F ss in B A (see Theorem Appendix I.7) with respect to a reference leaf. These properties could be used in order to find an improved definition of T RR (as discussed later). It follows that T RR diverges to −∞ as it approaches A or its strong stable manifold W ss (A) (if it exists) which is the manifold associated to the smaller Lyapunov exponents (precise definition in [42]). This implies that during the phase space estimation W ss (A) becomes visible as e.g. in the carbon cycle example below. (See Appendix D.1, Appendix H and Appendix I.2 for precise definitions, proofs and further properties.) In contrast to the convergence of the distance between trajectories in isochrons [28,29] this work focuses on the distance to the attractor, giving rise to this timing of transients.

Examples
In order to demonstrate the applicability of the metrics, we selected four examples with differing properties. They have been chosen with increasing levels of complexity and to show different properties of the new metrics. Note that 1-dimensional systems can be solved analytically and we discuss an additional example in Appendix D.1.

Linear system with two different time scales
The first example is the linear system in Eq. (3). We chose dimension 2 in order to show the basic features while still being able to map the phase space. But the results can be applied in any dimension.
While in general cases, T RR and audic can only be tackled numerically, we can solve the linear system analytically. The full details of this calculation are in Appendix B and the main results are (with D(x) = 1 2 While the result for audic seems intuitive, T RR might be more surprising. The result depends only on the x 1axis. This is because the strong stable manifold, i.e. the one corresponding to the smaller Lyapunov exponent, is in the x 2 -direction. But for −→ 0 in Eq. (2), which implies t −→ ∞ for a trajectory, the contribution from the smaller Lyapunov exponent disappears. Thus, only the orthogonal part x 1 is relevant. In order to get a better feeling for these novels metrics, we have chosen two exemplary initial conditions, an early-eager one and a late-eager one, and plotted their trajectories' distance to the attractor over time in Fig. 2. Thus, the blue-shaded area corresponds to the audic value which is the same in both cases of our particular choice. In order to demonstrate the intuition that T RR can be interpreted as the time-shift between the original trajectory and The initial conditions have been chosen such that the audic value, which corresponds to the blue-shaded area, is the same for both trajectories, D(x a 0 ) = D(x b 0 ) = 3.8. But the trajectory starting at x a 0 arrives earlier than the reference trajectory (green in (a) and (b)), which in turn is earlier than the one from 34. In order to show this, the example trajectories (blue) have been shifted in each plot by the value of T RR with respect to the reference trajectory. This demonstrates the intuition behind T RR : it describes by how much one has to shift one trajectory so it matches the asymptotics of the reference trajectory. The proofs in Appendix I provide that this is always possible for a generically chosen reference point x ref . Interestingly the manifold (red dashed line) where audic increases the strongest is tilted with respect to the center manifold. As the Rössler system is 3-dimensional, the above plot depicts only a slice at fixed z = 0.6. Furthermore, the dashed, red line marks the boundary of the attractor's projection to this plane.  The results from Eqs. (4) and (5) can also be seen in the numerical simulations in Figs. 3a and 3c. The coloring describes the values of the metrics (cmp. the colorbar in the right of the figures) and green star represents the reference point for the T RR computation. Note how the manifold where audic increases the slowest (red dashed line in Fig. 3a) is tilted with respect to the center manifold. This really proves that we are looking at transient phenomena and we have to take more into account than only the slowest dynamics.
The exponential lower bound that comes up in the correlation diagram Fig. 3b can be calculated analytically (see Appendix B)

Global carbon cycle
The second example is a conceptual model of the global carbon cycle proposed by Anderies et al. [3], where we used the pre-industrialization version. It consists of three dynamical variables, the terrestrial, maritime and atmospheric carbon stocks denoted c t , c m and c a respectively, and the constraint C = c t + c m + c a = const. Thus, we can reduce the system to 2 state variables c t and c m and rescale the units such that C = 1 arriving aṫ where N EP is the net eco-system production, p photosynthesis, r respiration, α harvesting parameter and I diffusion; indirect dependencies have been omitted and more details are in [3,30]. The whole phase space of Eqs. (7a) and (7b) as depicted in Fig. 3f is the basin of the attraction of a fixed point in the middle marked by a blue dot; the dynamics is drawn as streams. Note that the trajectories starting in the lower part have to pass by a "desert-like" saddle (with c t = 0) at the left (green dot).
The color in this graph depicts T RR and the first finding is the splitting of the basin of attraction. Furthermore, the strong stable manifold becomes visible as a light beige line due to their low values of T RR , i.e. as very early states because T RR diverges to −∞. This proves it being the separatrix for the observed splitting and it will merge to an arm of the stable manifold corresponding to the saddle arising after the subcritical pitchfork bifurcation mentioned below. Also, the expected smooth increase of the return times when distancing (along the trajectories) from the attractor can be observed.
When applying audic to this model (Fig. 3d) the splitting of the basin can be observed again. In contrast to T RR , the stable manifold is not visible because audic can be seen as a (by distance) weighted time and the contributions for the asymptotic part where the difference in the Lyapunov spectrum matters are negligible. Furthermore, we see a clear linear correlation of both metrics in Fig. 3e because all trajectories starting in the lower part have to pass by at the saddle on the left and spend a long time there.
This example shows how saddles can induce long transients, as stressed by Hastings [14], and that our metrics react appropriately.
It is important to note that the metrics are early-warning signals, too.
When increasing α, corresponding to the harvest of terrestrial carbon, the system passes through a subcritical pitchfork bifurcation where the saddle becomes stable and the lower-left part of the phase space splits off. The divergences of the two metrics' statistics as seen in Fig. 4 prove their prebifurcational sensitivity, while other systemic indicators like basin stability [21] do not change (up to numerical fluctuations, see Fig. 4).

Generator in a power grid
As an example of intermediate complexity, we chose the swing equations [4,22] in Eqs. (8a) and (8b), a model describing the dynamics of a single generator connected to a large power grid. It consists of two dynamical variables, the phase θ and angular frequency ω, both in a reference frame rotating at the grid's rated frequency. The parameters of the system correspond to the net power production P = 1 (at the node), the capacity of the transmission line K = 6 and dampening The stable fixed point at ω 0 = 0, φ 0 = arcsin P K describes a state of synchronization. For the chosen set of parameters, the system exhibits another attractor: a limit cycle at larger positive values of ω. For negative values, the two basins of attraction are interleaved. A more detailed introduction and analysis can be found in [22,31] Calculating T RR inside the basin of the stable fixed pointed (ω 0 , θ 0 ) yields Fig. 3i. There is basically no color change away from the attractor, so we can see that a trajectory will barely spend any time in the transient and goes quickly to the attractor. Analogously, Fig. 3g for audic leads to a similar conclusion.
Comparing both metrics, see Fig. 3h, shows that they are closely linked. Note that D is presented on a logarithmic scale, so the relation is exponential and can be explained using the calculations for a linear focus in Appendix B. So what we see here is actually the influence of the linearized part of the system.
Note that roughly 30% of the lower part of the phase space, where the nonlinearities actually have an influence, follow the exponential relation, too, because the transient is very fast and thus its influence is rather low.
The aforementioned limit cycle corresponds to the system being far away from synchrony and generators would usually switch off before reaching it. As it is not so relevant, the treatment of this attractor is in Appendix C.

Chaotic Rössler oscillator
Although we have proven the convergence for fixed points only, we show with the chaotic Rössler system [5,32] that our metrics are applicable to higherdimensional and more complex attractors alsȯ 3l shows a slice of the phase space with the standard parameters a = 0.2, b = 0.2, c = 5.7 for T RR and the expected sensitivity to initial conditions for chaos is observed: early and late trajectories lie closely together and the metric T RR has low spatial correlation.
In contrast, audic shows in Fig. 3j surprisingly smooth changes of an embryo-like shape. Because the focus of this article is on transient dynamics a new Figure 5: (color online) The bifurcation diagram (green) of the Rössler system for varying the parameter a in Eqs. (9a) to (9c) was computed from the local maxima in z of the attractor and T RR (orange) shows a strong sensitivity to these qualitative changes. The gray background is used so the reader can more easily connect the peaks in T RR to the corresponding parts in the bifurcation diagram.
feature of the chaotic Rössler system is uncovered: while the attractor is chaotic, the basin of attraction is very regular. audic focuses on the initial transient and the chaotic asymptotics is filtered out. For comparison, the boundaries of the attractor's projection have been added with dashed red lines in Fig. 3j and depictions of the attractor are in Appendix E.
We can deduce that even though the system is chaotic the strong sensitivity to initial conditions happens rather late in the transient when the trajectory is already close to the attractor, because T RR focuses more on the intermediate transient.
However, this implies that T RR can be successfully applied as an early-warning signal in this case, too. In order to demonstrate this, we chose to vary a as it has a crucial influence on the system's dynamics (see the bifurcation diagram in Fig. 5 (green)). While for values of a < 0.006 (cf. [33]) there is only a single stable fixed point, at a ≈ 0.006 a limit cycle emerges due to a Hopf bifurcation [33]. For a > 0.11, several period doublings are observed, finally leading to chaos for a > 0.155. Even in the chaotic regime, further bifurcations can be observed.
In Fig. 5, the standard deviation of the T RR distribution from randomly chosen initial conditions inside the basin of attraction is given. Due to the sensitive dependence on initial conditions, the reference value varies a lot and hence introduce shifts in the distribution that do not describe actual changes in the system's dynamics. To remove this effect, it is crucial to use central statistics like the standard deviation.
T RR is strongly sensitive to any qualitative changes in the dynamics of the system, incl. even chaoschaos transitions. Closely observing Fig. 5 uncovers that there is a base-line with a little noise at T RR ≈ 10 complemented with strong peaks. In the chaotic regime, the peaks correspond directly to qualitative changes. Also, we observe sensible changes during the period-doubling phase and a strong increase before the Hopf bifurcation at a ≈ 0.006, proving the usefulness as an early-warning signal.
The abrupt downward peak at a ≈ 0.11 is unexpected and more details are needed to clarify it.
Note that sometimes the result of estimating the limit numerically fluctuated slightly, suggesting that an improved definition of T RR could be sensible. The defining property could be the constant negative is the leaf of the strong stable foliation containing the reference point x ref . This turns out to be highly nontrivial and is an issue for future research.

Discussion and outlook
We have treated problems (I)-(IV) arising from common methods of analyzing the transient time needed to reach an attractor in nonlinear systems by introducing two complementing metrics: Area under Distance Curve and Regularized Reaching Time, and applied their properties, in particular them being Lyapunov functions, to the development of an efficient estimation algorithm. Furthermore, in Appendix I we prove under mild conditions the existence of these functions.
In order to show their applicability and usefulness, we chose 4 example systems and analyzed them with respect to these novel metrics.
In the linear system everything could be solved analytically and we found the tilt of the slowest increasing audic axis, demonstrating that even this simple system already can have a rich transient. Furthermore, we explained why the strong-stable manifold turns out to be a singularity of the system's T RR function but could argue that this is irrelevant for the estimation of statistics as the integral is still finite.
The global carbon cycle demonstrated the importance of the transient analysis, as the desert state is only a saddle but nevertheless passing by there would lead to an extinction of humanity. The artificially estimated splitting of the basin in [3] arises naturally here and T RR uncovers that the separatrix is the strongstable manifold of the system. Furthermore, the saddle induces an unexpected linear relation between audic and T RR . Particularly interesting is how the (central) statistics of our metrics are a systemic approach to the concept of critical slowing down [1,2,20] and react strongly to prebifurcational changes. Hence they are early-warning signals for fundamental changes of the system.
The generator in a power grid displays how even a very nonlinear system can have rather unimportant transients while the asymptotics is mostly relevant. Thus our measures show the exponential relation expected from a linear focus. In order to prove the applicability to more complex dynamics, we used our metrics on the Rössler system, too, and found the smoothness of the attractor's basin with audic. As the attractor itself is rather chaotic, this smoothness is surprising. Even though T RR reacts strongly to the sensitivity to initial conditions of the chaotic system its worth is displayed when varying the a parameter. This parameter has strong influence on the Rössler system's dynamics and T RR reacts strongly to the different bifurcations and even the chaos-chaos transitions, proving again its worth as early-warningsignal.
Furthermore, using audic as a cost or damage functions can be applied in the context of earth system analysis and climate impacts. In the future, applying this metric for estimating viable pathways without transgressing e.g. the Planetary Boundaries [34,35] is one goal of this work.
We did not perform any comparative analysis with the mentioned " −→ 0"-approaches because these behave inconsistently and their quantitative results are arbitrary, as discussed in length in Sec. 2 and Appendix A.
Four directions of immediate future research are due: (1) Improving the definition of T RR further using the Lyapunov function properties as described above. In order to generalize these definitions to even more complex systems this step seems crucial.
(2) Applying the metrics to more complex systems to understand them and their properties to topological structures, e.g. in complex networks [36], in more detail.
(3) Introducing more sophisticated methods of Lyapunov function estimations [37]. The curse of dimensionality is going to be a problem for network systems, hence methods for estimation these metrics statistics in such kinds of systems will need novel algorithms.
(4) Comparison of the timing of transients in model output and observation data.
With the observable time, new possibilities for comparison with observation data are available and should be used.

Acknowledgments
This paper was developed within the scope of the IRTG 1740/TRP 2011/50151-0, funded by the DFG/FAPESP. This work was conducted in the framework of PIK's flagship project on coevolutionary pathways (copan). The authors thank CoNDyNet (FKZ 03SF0472A) for their cooperation.
The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Land Brandenburg for supporting this project by providing resources on the high performance computer system at the Potsdam Institute for Climate Impact Research. The authors thank the developers of the used software: Python [38], Numerical Python [39] and Scientific Python [40].

Appendix A. Details on the problems of reaching time definitions
To understand the problems of known reaching time definitions in more detail, we introduce a small example system (Eqs. (A.1a) and (A.1b)) containing a saddle-node-bifurcation at a = 0 followed by a nodefocus transition at a = b 4 64 . Hence, for a > b 4 64 , a flow like the one depicted in Fig. A1a with a stable fixed point at x s = (x 1,s , x 2,s ) = (− √ a, 2(1 + b √ a)) and a saddle at x u = (x 1,u , x 2,u ) = ( √ a, 2(1 − b √ a)) is obtained. Note that this system was mainly chosen in order to present all occurring problems with just a single example.ẋ Problem (I), the divergences, can be observed when plotting the distance d of the example trajectory starting at some initial condition x 0 in Fig. A1a to the stable fixed point x s over time to get Fig. A1b. While d approaches 0 because x s is an attractor, it will reach the attractor in infinite time only, i.e. the reaching time diverges.
These problems are often addressed by measuring the time until the trajectory is "close" to the attractor. The relevant metrics would be the times when the trajectory enters an -neighborhood around x s the first time (t F (x 0 )) and the last time (t L (x 0 )), i.e. the neighborhood is not left anymore. To illustrate this, we turned Fig. A1b around and added these two metrics in Fig. A1c.
While the values for each of these metrics are now finite, their strong dependence on the choice of makes the physical interpretation (problem (II)) hard. Particularly, for the relevant small values of the times can become become as large as one wants due to their divergence for −→ 0. Furthermore, it is even hard to state when the trajectory is "close" or still in the "transient". From Fig. A1c, the discontinuities of the metrics (problem (III)) are observed, too. They go hand in hand with the physical interpretation problem as small changes in the choice of can induce jumps in the values of the time.
Problem (IV), non-invariance, can be understood by the fact, that a distance function is not invariant under a change of variables. This means, we would get different values for t F (x 0 ) and t L (x 0 ) if we change the variables. This principle can also be understood from a different point of view: Which is the correct function that one should choose to measure the distance to the attractor? This might sound trivial, but it is crucial for the results. As only values close to the attractor matter, but the metrics are very sensible to changes (problem (II)), the impact of this choice may become large.

Appendix B. Analytical solutions of a linear system
Appendix B.1. Area under Distance Curve For a linear system, we chose dimension 2 in order to show the basic features while still being able to map the phase space. But the general procedure can be applied in any dimension.
The systems isẋ where x,ẋ ∈ R 2 and A ∈ R 2×2 . This implies a fixed point at x * = (0, 0) T (that we will abbreviate with 0 from now on) and we want A to be (complex) diagonalizable and with negative real parts of all eigenvalues. Hence 0 is exponentially stable.
In this part we discuss the case of a real leading eigenvalue and the case with complex leading eigenvalues can be treated similarly while keeping Proposition Appendix I.5 in mind. For reasons of simplicity we use as distance function d(x, 0) = x 1 2 + x 2 2 . (Note that this is not a distance in mathematical terms but still captures how close on is to the attractor. This has been discussed in the main text and in Appendix G.) An ansatz for audic is D(x) = 1 2 x T Cx with a symmetric matrix C ∈ R 2×2 . Plugging this in the orbital derivative equation for audic leads to This implies that the asymmetric part of CA equals negative unity. As this should be true for all A (including non-symmetric ones) the negative inverse of A is in general not the solution for C. Instead we note that for any asymmetric matrix Q ∈ R 2×2 , d(x, 0) can be written as Thus, we shifted the problem to finding an asymmetric matrix Q s.t. (1 + Q)A −1 is symmetric so we have To solve this, we introduce d(d − 1) parameters for the upper-right triangle of Q (because the lower-left one is then defined by the constraint Q ij = −Q ji ).
In the case of a 2-dimensional system, this is one parameter q: x. (B.8)

Appendix B.2. Regularized Reaching Time
Calculating T RR in the linear case can be done, too. The first step is to chose an eigenvector basis {v i } i=0,...d−1 corresponding to the eigenvalues ordered as Reλ 0 > Reλ 1 > . . . of A and find the (unique) decomposition of x in that basis: We can plug this in the solution of the linear system: (B.10) Taking the limit as in Eq. (H.1) looks at large values of t s.t. the terms for i ≥ 1 can be omitted, because the λ i have been ordered in an appropriate way: For a fixed distance value = x(t) this equation can be solved for t yielding the reaching time: Plugging this in the definition for T RR (using a reference point Note that we did not use the audic level sets as described in the main text. This could be done here but only lengthens the calculation while the result stays the same. For the example in Eq. (3) the matrix has been used, which has as eigenvalues and -vectors The decomposition of Eq. (B.9) is In order to understand the exponential lower bound seen in Fig. 3b, we expand Eq. (B.8).
Because T RR in Eq. (B.18) depends only on x 1 we need to find the lower bound of D for fixed x 1 , i.e. the minimum of the polynomial in x 2 . We find  Fig. C1c). This attractor corresponds to the generator being far from synchrony with the grid and usually it would have been switched off long before reaching it, so it is rather irrelevant and we analyze it here for completeness only. Again, as with the fixed point attractor discussed in Sec. 4.3, we observe in Figs. C1a and C1c that the influence of the transient is rather low and the system approaches the attractor exponentially. Looking at Poincaré maps for fixed φ would give asymptotically exponential behavior. Thus the exponential relation in Fig. C1b can be explained with the calculations for the linear system in Appendix B.

Appendix D. One-dimensional Systems
Appendix D. 1

. Closed Formulas
In one-dimensional systemṡ with a stable fixed point at x * , closed formulae for the T RR and audic can be written down. By separation of variables we get where x 0 is initial state, the difference to the fixed point, ± needs to be chosen depending on the side of x * where the initial state x 0 is and t(x 0 , ) the time. Then the Regularized Reaching Time is In the same manner, the closed expression for audic in one dimension can be derived: where d(·, ·) is the chosen distance function.

Appendix D.2. Quadratic Correction
In order to demonstrate how the derived equations can be applied, we analyze this system which has a linear term plus a quadratic correction: The attractor in this system is a fixed point at x * = 0 and the corresponding basin of attraction is B(0) = (−∞, 1/b). T RR can be calculated in a straightforward manner using Eq. (D.4) and yields in Eq. (D.6b). This result is depicted in Fig. D1 (with different reference states x ref and compared with the linear result) and one can observe a different behavior one each side of x * . Particularly relevant are the negative divergence at x * = 0 and the positive one at x b = 1/b. The latter fixed point is the boundary of the basin of attraction and hence never reaches x * . So we expect T RR in the limit to x b to diverge and this can be seen in Fig. D1.
Furthermore, the curve intersects the x-axis at x ref as expected.
audic can be computed with the closed expression from Eq. (D.5), too: This result is depicted in Fig. D1. Particularly the divergence at x −→ 1 b is visible and due to reaching the basin boundary, analogously to T RR .
At the attractor x * = 0 audic reaches 0 because it never deviates from there, thus the cumulative distance vanishes.   The first problem is the estimation of the attractor itself. Finally, this could be solved starting at various different points and numerically integrate for a very long time. Removing the transient and then sampling the trajectories lead to a good estimate that could still fit in the available memory.

Appendix E. Rössler attractor
The second problem was the distance estimation. Fortunately, KD-Trees are exactly made for this and implemented in Scientific Python [40].
With these ingredients, audic could be calculated. In order to estimate T RR , the cumulative distances for the points along the trajectory were calculated backwards. This gives us the corresponding levelset of audic for each point on the trajectory. With this, the times for entering different audic levelsets could be retrieved and compared with the reference trajectory. Thus several values for T RR were obtained and a limit could be estimated.

Appendix G. Convergence of AuDiC
The convergence of audic depends on two elements: the distance function d and the asymptotic approaching behavior of the trajectories. A common case with a mathematical distance function, i.e. a function fulfilling Eqs. (G.1a) to (G.1d) [41] and an exponentially stable attractor, the convergence can be proven right away as the bigger integral over the exponential envelope converges.
Note the usage of the word mathematical distance (fulfilling the four properties). As we use the distance function only to measure how far a point is away from the attractor, more general functions could be used as well as long as they converge to 0 when a trajectory approaches the attractor. Particularly, cost or damage functions that could be well motivated from the system's context are unlikely to always fulfill the requirements of a mathematical distance function. Even assuming a mathematical distance function, convergence is not necessarily given. Systems that converge slower than exponentially could lead to a divergence in audic. A simple example of such a case isẋ = f (x) = −x 3 2 . The solutions are ± 1 √ c+t where the constant c is fixed by the initial condition. Using the absolute values as the distance function gives a divergence for audic, as The simplest case is systems with finite reaching times because the RHS of Eq. (2) can be split in two limits that converge separately. The results is the difference of the actual reaching times and is expected from the approach. Still, as the reaching times are finite anyway the complex approach with T RR is not necessary and is just to show that it is reasonable in these cases, too.
For infinite reaching times, the values for t audic will become increasingly large in the limit → 0. So the limit in Eq. (H.1) exists only, if for small (i.e. large times) the changes in t audic (x 0 , ) and t audic (x ref , ) will be about the same. This means, that for two different, small 1 > 2 Turning this interpretation around, it means that the trajectories have to behave "similarly" in the asymptotic limit, i.e. close to the attractor. The simplest case is a system with a hyperbolic fixed point where the larges eigenvalue of the corresponding Jacobian is real and of multiplicity 1. If that is the case, the calculation in Appendix B can be used locally around the attractor to understand why it converges and the precise proofs are in Appendix I.2. Having a multiplicity larger than one might be mathematically interesting but is physically rather unlikely because some slight differences in the modeling of the system would usually change these. If this is a persistent property of the system, a precise understanding of the meaning is needed. Note that T RR still converges but will depend on the underlying distance function. The latter is equivalent to problem (IV), non-invariance, and hence the result for these systems should be interpreted with care.
In the case of the largest eigenvalue being complex, the convergence can still be proven but there is a need for the choice of a specific distance function, as shown in Appendix I.2. This is not problematic as the result is also invariant under change of variables and hence, simply the unique result that can be taken.
On the other hand, this suggest that our current definition might have to be improved as for more complex and higher-dimensional attractors we can currently rely only on numerics. The results for the Rössler attractor, particularly the strong sensitivity to the dynamics of the system as shown in Fig. 5 provide the numerical support for our current approach.
An improved version of the definition could be done using the properties as Lyapunov functions with constant negative orbital derivate. Even though we can give a rough outline of how to do that in the Discussion of the main paper, there are many subtle technicalities to be addressed in order to define an improved T RR precisely.

Appendix I. Precise Definitions and Theorems
We consider a deterministic dynamical system of the formẋ and assume that the system contains an exponentially asymptotically stable equilibrium A. That is, f (A) = 0 and all the eigenvalues of Df (A) have negative real part. We denote the basin of attraction of A by B A . Also we denote the flow operator of (I.1) as ϕ(t, ·). Let the spectrum of Df (A) be λ s ∪ σ ss , where λ s may be real or complex, but we assume has multiplicity one. It will also be useful later to define constants α s and α ss such that where we also require 2|α s | > |λ s |. Then the following result holds. We refer to [6, Theorem 2.46] for a proof.
Proposition Appendix I.2. Let d(·, A) be a C 1 function and suppose that there is a class K function That is, the AuDiC function is a Lyapunov function with orbital derivative equal to −d(x, A).

Appendix I.2. Regularized return time (RRT) function
In the following, we further assume that the function d(·, A) from equation (I.2) has been chosen such that d(x, A) = ||x − A|| N for some norm || · || N . Then it is clear from (I.2) that the AuDiC function defines a norm || · || D given by D( For a given initial condition x r ∈ B A , we denote the time taken for an initial condition to enter and remain inside a D-ball of radius as t audic (x 0 , ) := inf{T : ||ϕ(t, x 0 )−A|| D < for all t ≥ T }.
The T RR function is then defined as follows.
Definition Appendix I.3 (T RR function). For a given reference point x r ∈ B A , the T RR function is defined as where the limit exists.
A natural question is under what conditions the limit in (I. 3) exists. To answer this question we distinguish between two cases according to whether λ s is real or complex. Our first result is regarding the existence of the T RR function in both cases, and the dependence on the choice of norm || · || D . In order to state the result for λ s complex, we make the following definition.
Definition Appendix I.4. We define the following equivalence class on norms defined on B A : ) where E s is the invariant subspace corresponding to the leading eigenvalue of Df (A).
Clearly, elements in the above equivalence class are defined by the norm of elements in E s . Proposition Appendix I.5. Let the T RR function be defined as in Definition Appendix I.3 for the system (I.1), and assume x 0 ∈ W ss (A). Then we have the following (i) When λ s is real, the limit (I.3) exists for all choices of norm || · || N . Moreover, the limit is independent of the choice of norm.
We will also show that the T RR function is closely related to the strong stable foliation F ss in the basin of attraction of the equilibrium A. We first recall the following definitions.
Definition Appendix I. 6. A foliation F of an d-dimensional manifold M is a partition of M into a disjoint collection of k-dimensional injectively immersed connected submanifolds (called leaves) such that for each x ∈ M , there is a neighborhood V ⊂ M and a chart φ : V → R k × R d−k , such that each connected component of the intersection of a leaf of F with V is mapped to the set R k × {y}, for some y ∈ R d−k . We call F a C r foliation if each local chart is C r . A continuous foliation whose leaves are C r is called a C r lamination.
We denote the leaf of a foliation through a point x as F(x). A foliation F is invariant under the flow of (I.1) if ϕ(t, F(x)) = F(ϕ(t, x)) for sufficiently small |t|.
Theorem Appendix I.7 ( [42]). Consider the system (I.1) and let R d = E s ⊕ E ss be the direct sum decomposition into the invariant stable and strong stable subspaces for the linear systemẋ = Df (A)x. Then there exists a unique invariant C r lamination F ss in B A , called strong stable foliation, such that each leaf of F ss has dimension equal to dim E s and F ss (A) = W ss (A), where W ss (A) is the strong stable manifold of A.
Solutions x(t) and y(t) that belong to the same leaf of F ss for all time are characterized by strong asymptotic convergence to each other: ||x(t)−y(t)|| D ≤ Ce −α ss t for t sufficiently large.
We will also prove the following result which provides an important characterization of the T RR function.
Proposition Appendix I.8. Let the T RR function be defined as in Definition Appendix I.3 for the system (I.1), and assume x r ∈ W ss (A). In the case λ s is complex, we assume || · || D ∼ || · || P as in Proposition Appendix I.5. Then the level sets of T RR (x 0 ; x r ) are equal to the leaves of F ss . That is, T RR (x 0 ; x r ) = T RR (y 0 ; x r ) ⇔ F ss (x 0 ) = F ss (y 0 ) Furthermore, T RR (x 0 ; x r ) → −∞ as x 0 approaches F ss (A) (= W ss (A)).
The proof of Propositions Appendix I.5 and Appendix I.8 rely on the following result regarding the behavior of solutions in the approach to equilibrium. We refer to [43] for a proof.
Theorem Appendix I.9. Consider the system (I.1), and define λ s , α s and α ss as before. Then there exists κ > 0 such that for all solutions x(t) of (I.1) in B A with ||x(0) − A|| D < κ, the limit η(x(0)) := lim Note that since Φ(0, t) leaves E s invariant and E s is closed, we have η(x(0)) ∈ E s . It also follows from the proof in [43] that η : B κ (A) → E s is continuous.
Lemma Appendix I. 10. Let x(t), y(t) be solutions to (I.1) in B A . Then x(t), y(t) belong to the same leaf of F ss for all t if and only if η(x(s)) = η(y(s)) for s sufficiently large. Furthermore, η(x(s)) = 0 if and only if x(t) ∈ W ss (A).
Proof. By Theorem Appendix I.7, the solutions x(t) and y(t) belong to the same leaf of F ss for all t ≥ 0 if and only if ||x(t) − y(t)|| D ≤ Ce −α ss t . Let s > 0 be large enough so that ||x(s) − A|| D , ||y(s) − A|| D < κ. Now from Theorem Appendix I.9 we have η(x(s)), η(y(s)) ∈ E s , and equation (I.6) implies that this is possible if and only if η(x(s)) = η(y(s)). The last statement follows directly from (I.6) and the theory of stable/unstable manifolds.
Proof of Proposition Appendix I. 5. Let x 0 and x r be as in Definition Appendix I. 3 and let x(t),x(t) be the solutions to (I.1) with x(0) = x 0 andx(0) = x r . Let s > 0 be large enough so that ||x(s) − A|| D , ||x(s) − A|| D < κ, and let > 0 be small. Then from Theorem   The final statement of Proposition Appendix I.8 follows from η(x(s)) = 0 ⇔ x(t) ∈ W ss (A) (see Lemma Appendix I.10) and the fact that η : B κ (A) → E s is continuous.