Climate change in mechanical systems: the snapshot view of parallel dynamical evolutions

We argue that typical mechanical systems subjected to a monotonous parameter drift whose timescale is comparable to that of the internal dynamics can be considered to undergo their own climate change. Because of their chaotic dynamics, there are many permitted states at any instant, and their time dependence can be followed—in analogy with the real climate—by monitoring parallel dynamical evolutions originating from different initial conditions. To this end an ensemble view is needed, enabling one to compute ensemble averages characterizing the instantaneous state of the system. We illustrate this on the examples of (i) driven dissipative and (ii) Hamiltonian systems and of (iii) non-driven dissipative ones. We show that in order to find the most transparent view, attention should be paid to the choice of the initial ensemble. While the choice of this ensemble is arbitrary in the case of driven dissipative systems (i), in the Hamiltonian case (ii) either KAM tori or chaotic seas should be taken, and in the third class (iii) the best choice is the KAM tori of the dissipation-free limit. In all cases, the time evolution of the chosen ensemble on snapshots illustrates nicely the geometrical changes occurring in the phase space, including the strengthening, weakening or disappearance of chaos. Furthermore, we show that a Smale horseshoe (a chaotic saddle) that is changing in time is present in all cases. Its disappearance is a geometrical sign of the vanishing of chaos. The so-called ensemble-averaged pairwise distance is found to provide an easily accessible quantitative measure for the strength of chaos in the ensemble. Its slope can be considered as an instantaneous Lyapunov exponent whose zero value signals the vanishing of chaos. Paradigmatic low-dimensional bistable systems are used as illustrative examples whose driving in (i, ii) is chosen to decay in time in order to maintain an analogy with case (iii) where the total energy decreases all the time.

(ii) either KAM tori or chaotic seas should be taken, and in the third class (iii) the best choice is the KAM tori of the dissipation-free limit. In all cases, the time evolution of the chosen ensemble on snapshots illustrates nicely the geometrical changes occurring in the phase space, including the strengthening, weakening or disappearance of chaos. Furthermore, we show that a Smale horseshoe (a chaotic saddle) that is changing in time is present in all cases. Its disappearance is a geometrical sign of the vanishing of chaos. The so-called ensemble-averaged pairwise distance is found to provide an easily accessible quantitative measure for the strength of chaos in the ensemble. Its slope can be considered as an instantaneous Lyapunov exponent whose zero value signals the vanishing of chaos. Paradigmatic low-dimensional bistable systems are used as illustrative examples whose driving in (i, ii) is chosen to decay in time in order to maintain an analogy with case (iii) where the total energy decreases all the time.

Introduction
One of the most relevant phenomena experienced by everyone in our days is climate change. The scientific background is that changes in the atmospheric composition caused by greenhouse gases and aerosols have been leading to a decades-long monotonous increase in the globally averaged surface temperature as detailed in renowned sources, like the IPCC report [75]. The so-called arctic amplification mechanism [14] leads to a stronger temperature increase in the polar region than anywhere else since the major cause of albedo reduction is the loss of sea ice. All this gives rise to a monotonous decrease in the temperature contrast between the pole and the equator.
Atmospheric and oceanic dynamics is turbulent, weather, therefore, exhibits chaos-like properties, and is unpredictable [16,21,40]. The "chaotic attractor" of the climate is, however, high dimensional, one is unable to plot it and consider its size as a measure of unpredictability. For a qualitative view, one can accept that there is a multitude of possible states, and in a stationary case, one can speak of parallel states of the climate. This is a concept that helps interpreting the term "internal variability" widely used in climate science as a synonym of chaos [75]: even if only a single state is observable at a given time instant, many others are also permitted due to the chaotic nature of the dynamics. In a stationary climate, driven by annual periodicity, these states are exactly the same on a given day of the year.
Climate change is, however, the consequence of the aperiodic nature of the driving, e.g., of the temperature contrast between the pole and the equator, or of the CO 2 concentration, whose amplitude increases at a non-negligible rate year by year. A recently spreading new view among researchers is that climate should be characterized by a set of permitted climatic histories, parallel climate realizations [35,78]. A proper description of the time-dependent internal variability is then only possible if one monitors the full plethora of permitted histories. This implies the simultaneous monitoring of many time evolutions within the same physical system. In other words, ensembles of trajectories should be considered which differ in their initial conditions only. Note, that this ensemble becomes independent of these initial conditions after a short time, due to dissipation. An instantaneous state of the system can be characterized by taking a snapshot of the ensemble at a given time instant; the term snapshot view is therefore also in use.
This concept was born in dynamical system theory as part of a progress made in understanding non-autonomous dynamics. The first article in this direction [73] drew attention to an interesting feature: a single trajectory leads to a washed out appearance on a map, while an ensemble of trajectories subjected to the same realization provides a structured pattern. This ensemble-related pattern, the snapshot chaotic attractor, changes its shape all the time in contrast to traditional chaotic attractors which are time-independent if observed on a stroboscopic map for periodic forcings [66]. The concept of snapshot attractors has been used to understand a variety of time-dependent physical phenomena (see, e.g., [28,42,[52][53][54]64,74,80,83]). Meanwhile, a related concept, that of pullback attractors, referring to the union of all snapshot attractors along the time axis, arose in the mathematics literature [13,50]. The spread of these ideas became accelerated by the fact that their relevance in climate dynamics was discovered more than 10 years ago, using the language of pullback attractors [8,24,26,27]. Later, it also became clear [18,43,69,70] that in the context of climate change, not noisy driving, rather a continuous drift of parameters is relevant, and therefore, the time evolution of deterministic snapshot attractors captures the essence of climate change.
An important consequence of these studies is that the traditional "single trajectory" view of a chaotic attractor is not equivalent to the "ensemble" view [19], and out of the two approaches, it is the ensemble-based snapshot view that is appropriate for a faithful statistical description of a changing climate. The reason is that the ensemble also represents a natural probability distribution, supported by the snapshot attractor. Averages should be taken with respect to this distribution at any time instant. The measure supported on the time-dependent attractor can also be reconstructed using response theory as long as one does not hit a tipping point [57,58].
The basic features of the snapshot view can be summarized as follows: -in non-autonomous systems the infinite time limit is meaningless since models leave their range of validity, only finite-time observations are feasible, as emphasized in the context of flows of general time dependence (see, e.g., [30]), -conclusions based on single trajectories are not reliable since such trajectories are unpredictable in the long run, and thus not representative, -ensemble properties, on the contrary, are fully predictable (in harmony with general properties of chaotic systems [77]), including the natural probability distributions, -an instantaneous characterization of the system becomes possible, the use of temporal averages can fully be avoided, -a straightforward way is offered to analyze internal variability in a changing system by means of, e.g., statistical quantifiers of the instantaneous probability distribution of the ensemble.
In the particular example of climate change, the snapshot view can equivalently be formulated as the theory of parallel climate realizations [35,78]. Qualitatively speaking, one imagines many copies of the Earth system moving on different dynamical paths, being subjected to the same physical laws and drivings. It is interesting to note that this is a generalization of an idea of Leith [56] published in 1978, augmented by the theory of snapshot attractors.
The aim of this paper is to emphasize that the knowledge accumulated in climate science can be taken over into classical mechanics, or dynamical systems theory, since many features of non-autonomous systems are yet to be explored. The fact that mechanical systems are often of paradigmatic nature and thus provide great models for more complicated systems, adds further significance to this approach. Climate of our days is nothing but a system subjected, in first approximation, to a time-dependent external non-periodic driving. The timescale on which the amplitude of driving monotonically changes is not much smaller than the basic internal timescales of the system in the lack of parameter drift. In the case of a real climate, this basic internal timescale can be, e.g., a year. Based on this, we can say that a mechanical system undergoes a climate change if it is subjected to a monotonous (or piecewise monotonous) parameter drift whose timescale is comparable to that of the internal dynamics. If the system has multiple time scales, as the real climate, the timescale of the parameter drift should be away from either asymptotics, i.e., from the adiabatically slow change or from the sudden jump. The constraint on monotonicity is prescribed to emphasize the similarity with the climate change we are witnessing nowadays, while in a more general context the constraint of monotonicity can even be omitted. The authors of a recently published paper [67] speak about the evolution of a system's climate in a similar sense, and they do allow for non-monotonous parameter changes, although with a given trend on average, as well as for a parameter jump.
A class of mechanical processes where the analogy with climate change is immediate is that of nonlinear oscillators driven by varying amplitude. Seasonality of the climate corresponds to the presence of a periodic component of the driving; the time-varying driving amplitude is the analog of the change in atmospheric carbon dioxide concentration or the changing temperature contrast between the pole and the equator.
For such systems, we can take over important aspects learned from climate science. The dynamics of a general system is chaotic-like, but also has external time dependence and since the tools of standard chaos science, e.g., periodic orbit theory [12], are only defined for stationary systems, they cannot be used in the presence of parameter drift. Such systems, however, also possess a plethora of permitted states, i.e., internal variability, whose time dependence can be described by following an ensemble, i.e., by accepting the snapshot view. Statistical quantifiers can be based on the natural distribution defined by the ensemble, e.g., averages or momenta, and they do depend on time. If a plausible picture is needed behind the snapshot view, we can speak of the theory of parallel dynamical evolutions of the systems, in analogy with parallel climate realizations. A recent effort in this direction is the application of the snapshot view to epidemics under changing environmental conditions [51], to aperiodically driven Hamiltonian systems [38,39], to tipping phenomena [7,47], and to advection in open flows of changing intensity [81].
The theory of parallel dynamical evolutions can usefully be applied to several branches of dynamical systems. Here, we illustrate this with three classes of mechanics: -(i) Driven dissipative systems subjected to parameter drift (Sect. 2). Since the Earth system dynamics is dissipative, the analogy with climate science is the strongest in this class (although spatial aspects do not show up). We can initiate a large ensemble widely distributed initially in the phase space and follow its evolu-tion. It converges, typically exponentially, to the snapshot attractor which then evolves in time. After convergence took place to the snapshot attractor, statistical averages, e.g., Lyapunov-exponent-like quantities, can be determined. The dynamics can be followed up to nearly vanishing drivings, and thus, the question of how chaos dies [45] can also be investigated here.
-(ii) Driven Hamiltonian systems subjected to nonadiabatic parameter drift (Sect. 3). Attractors do not exist in such systems (the system does not forget its initial conditions), but a lesson of climate science is worth taking over, we need to follow trajectory ensembles. Monitoring an extended ensemble would lead to an unstructured view of the dynamics. To a clear characterization, special sub-ensembles turn out to be useful, corresponding to KAM tori of the stationary phase space. As time goes on, these ensembles take up different shapes in every instant, thus we call them snapshot tori. Outside tori, time-dependent chaotic regions exist, called snapshot chaotic seas [38,39] and they can also be used as sub-ensembles.
-(iii) Undriven dissipative systems (Sect. 4). This class is not subjected to driving and therefore all motion in it stops ultimately. It is the total energy which decreases monotonically in time, and plays thus the role of a kind of driving. All motion, and thus chaos, too, is unavoidably transient. Since, however, the complexity of the transients is itself changing, the term doubly transient chaos has been coined [10,61]. Here, we show that the use of sub-ensembles opens up the skeleton of the motion, when one uses the KAM tori of the frictionless system for this purpose, which are then immediately subjected to the dissipative dynamics.
To provide a treatment easy to follow, examples are taken from the family of double-well potentials. In order to keep similarity with the monotonic energy decrease in the third example, we shall investigate drivings of decreasing strength in the two other cases. The simplicity of the models ensures that some analytic calculations are also possible. In the first two classes, we take advantage of the fact that moving points can be identified, at least in a perturbative approach, which keep their stability character (hyperbolic, nodal or elliptic) during long stretches of time. The temporally changing manifolds of such snapshot hyperbolic points are shown to play important roles. The unstable manifold is interwoven with the snapshot attractor all the time, while the stable one is found to be part of the boundary between time-dependent basins of attractions. Both manifolds are constituents of the snapshot chaotic sea in Hamiltonian cases. Furthermore, we shall demonstrate the presence of Smale horseshoe type structures, snapshot horseshoes [45,46,55], traced out by stable and unstable manifolds of snapshot hyperbolic points, and find evidence for its existence in doubly transient chaos, too. The existence of this object, which can also be called a snapshot chaotic saddle, is the common feature of all three cases (iiii) investigated. This object changes with time, and designates changes in the strength of chaos. In all cases, a geometric signature of the dying out of chaos is the disappearance of the horseshoe, i.e., the disappearance of a large number of intersection points between the manifolds. On a more quantitative level, we show that the pairwise distance of trajectories averaged over the parallel dynamical evolutions, the socalled ensemble-averaged pairwise distance (EAPD), is a quantity which faithfully represents the strength of chaos in all cases by providing instantaneous values of the Lyapunov exponent as well. This method enables us to identify the instant when chaos disappears as the time when this Lyapunov exponent vanishes.

The model
We consider a paradigmatic low-dimensional example, the periodically forced Duffing oscillator (see, e.g., [1,41]): The undriven ε = 0 case corresponds to the standard problem of the motion in a double-well potential with minima at x = ±1. The time unit is set here by the linear term, the length unit by the cubic term, and ω and β represent a dimensionless frequency and damping coefficient, respectively. We perform a stroboscopic mapping by considering the instants t = T n, where T = 2π/ω is the driving period, and n is a positive integer. The periodic term ε cos ωt represents in the climatic context annual driving. With a constant ε, the system's climate is stationary. For usual parameter choices (e.g., β = 0.125, ε = 0.3, ω = 1), the system possesses a typical extended chaotic attractor, and the parameter dependence of the location of the attractor(s) is reflected by a complex Fig. 1 Bifurcation diagram of the driven Duffing system (1) with fixed driving amplitudes ε, and β = 0.01. The x coordinates on the attractors are plotted in the stroboscopic map (omitting the first 450 iterations). The upper and lower straight branches indicate the existence of two period-1 attractors whose locations tend to the potential minima ± 1 for vanishing amplitude. They overlap with the orange lines, which correspond to equation (7) for α = 0. There exist two more period-1 attractors (bent curves) representing large oscillation in each potential well which are born in a saddle-node bifurcation at ε ≈ 0.01. A third one appears at about ε = 0.05 and describes a period-1 oscillation surrounding both wells. The phase space view of these orbits and of a period 5 and 7 attractor existing in relatively long ε intervals are indicated in insets bifurcation diagram [1]. Here, we keep ω = 1, but choose a small damping constant β = 0.01 in order to find interesting dynamics for small driving amplitudes ε < 0.1. This shall enable us to compare the numerics 1 with analytic expressions obtained in the limit of weak driving for a few simple orbits valid even in the presence of parameter drift (see Sect. 2.2). The bifurcation diagram obtained in this range is shown in Fig. 1.
The diagram might appear boring in the lack of dominant chaotic attractors. It would, however, be erroneous to conclude that there will be no chaos observed when the driving amplitude is time-dependent and is scanning through this range. The point is that usual bifurcation diagrams can reflect only the existence of permanent chaos. The picture might be different if the transient form of chaos [55] is not excluded. It was pointed out in [45] that transient chaos characterizing the system with fixed driving amplitudes can give rise  Fig. 1, but only the first 50 iterates of the trajectories are omitted. 900 initial conditions distributed uniformly in −1 ≤ x ≤ 1 have been followed for any fixed driving amplitude. Note how much more empty the traditional diagram in Fig. 1 is. The stroboscopic view of a few chaotic saddles belonging to the ε values 0.03, 0.045, and 0.09 is shown in insets to a chaotic snapshot attractor governing the dynamics in the presence of parameter drift. Therefore, we have produced another bifurcation diagram enabling one to see long lived, and thus perhaps chaotic transients. The result visible in Fig. 2 indicates a cloud of scattered points even far away from the long-term attractors shown in Fig. 1. This supports that there are long-lived transients nearly in the whole interval of the driving amplitude investigated. In order to see if they are chaotic, we determined the underlying chaotic saddle at different fixed driving amplitudes by means of the so-called sprinkler method [55]. Three examples are displayed in insets and suggest that the extension of the saddle decreases with decreasing ε. Our investigations imply that the transients for ε 0.025 are non-chaotic.
We introduce the parameter drift converting the driving amplitude ε to be time dependent ε → ε(t). When the timescale on which ε(t) changes is comparable to that of the original system (1), i.e., to 1/ω or 1/β and is monotonous (or at least piecewise monotonous) we say that the system undergoes its own climate change. In more detail: both climatic and mechanical evolutions are sensitive to periodic (seasonal) changes and to the drift in the amplitude of this external driving.
We consider a linearly time-dependent term representing a decreasing ramp ending in a plateau of zero driving amplitude described by the differential equation where α is the rate of the parameter change.
We define a scenario, starting at time t = 0 as the evolution of the driving amplitude set by parameters ε 0 and α, and the number n of periods up to which dynamics (2) is followed. For the sake of seeing the endpoint of the ramp appearing on one of the stroboscopic sections, we set the rate so that the ramp ends at a multiple integer n 0 of the period, i.e., we set α = ε 0 /(n 0 T ). The number of periods spent on the plateau after the end of the ramp is denoted by n = n − n 0 , for n > n 0 .

Snapshot hyperbolic points and snapshot nodes
Unstable periodic orbits are considered to be the skeleton of traditional chaos, the form of chaos appearing in dissipative systems subjected to constant or periodic temporal driving (see, e.g., [12,66]). An infinite number of hyperbolic cycles exist in such systems, and the motion is imagined to be a random walk among the periodic orbits. Statistical properties, like, e.g., the average Lyapunov exponent, can be expressed by means of the local cycle eigenvalues.
In systems subjected to aperiodic driving periodic orbits do not exist. One can nevertheless imagine that a snapshot attractor contains the union of an infinite number of special aperiodic orbits which are analogs of the periodic ones in the sense that they keep their (in)stability practically unchanged over long stretches of time. That is, they are surrounded by regions that also remain (un)stable for a long time interval. In the simplest and most common case, such stable and unstable orbits grow out of periodic orbits of the parameter-independent case. Such points with hyperbolic character shall be called snapshot hyperbolic points (SHPs), as introduced in [38]. Note that analogous objects have been considered in the literature (see, e.g., [22,37,59,84]), termed as distinguished hyperbolic trajectories, hyperbolic cores, or moving hyperbolic points. We intend to keep the term snapshot to emphasize that there is a full set of snapshot objects in any system exhibiting climate change. In fact, we shall see that snapshot nodal points (SNPs), i.e., timedependent attractor points, also exist. In the Hamiltonian limit, they go over into snapshot elliptic points (SEPs).
The simplicity of our model makes the explicit determination of such snapshot points possible when the dimensionless ε(t) is small during the scenario, an assumption valid in our examples. One of the SHPs is inherited from the unstable fixed point where v =ẋ, of the undriven case. Its behavior can be understood by considering the linear equation obtained from (2) by neglecting the cubic term in it. This is an inhomogeneous linear equation (for details on the solution, see "Appendix A"). Its coordinates in the phase space are found to be: The expressions of the parameters A H , B H , φ H , and θ H are given in Appendix A. This snapshot hyperbolic point moves along a curve spiraling about the origin of the phase space. On the stroboscopic map, we find a point which moves uniformly along a line in the plane of the map. The analytic expression for a pair of snapshot attractor points, snapshot nodal points (SNP) can be obtained in an analogous way by expanding the equation of motion (2) about the stable locations x = ±1, v = 0 of the undriven system (for details see Appendix A). This point again keeps its stability while moving for long times. For the x, v coordinates of the SNPs, we obtain: while on the stroboscopic map we find Parameters A N , B N , φ N , and θ N are also given in Appendix A. In continuous time, these points spiral about x = ±1, and move along a straight line on the map. The accuracy of the approximation used is illustrated by the fact that the nodal point coordinates of the unramped case ((7) evaluated for α = 0) coincide with the numerically determined attractor point positions in the full ε range investigated as Fig. 1 illustrates.

Snapshot attractors
A snapshot attractor is obtained numerically by initializing a large number of trajectories in an extended part of the phase space, typically uniformly distributed, following this ensemble of trajectories up to a certain instant, and plotting their location at this snapshot. Due to dissipation, the ensemble forgets its initial shape and distribution at an exponential rate. Therefore, after a convergence time, the outcome can be considered the snapshot attractor governing the global dynamics. This attractor can of course change its shape, and the distribution on it, all the time.
To follow the dynamics of the dissipative Duffing oscillator subjected to a change of the driving amplitude according to scenario (3), we distribute N = 10 4 points in the range shown in Fig. 3a and show the shape of the ensemble at certain instances in the subsequent panels of Fig. 3. An independent investigation reveals that the convergence time is about 5 periods (about convergence, more can be read in [20]); thus, the subjects visible in panels (b-e) are indeed snapshot attractors. In the language of parallel dynamical evolutions, they represent the union of instantaneous states along all the permitted histories.
The shapes in panels (b) and (c) are fractal-like, as filamentary as usual chaotic attractors. The extended shapes indicate that a huge number of states is permitted at these instances, which can also be interpreted as a sign of chaos. In an analogy with climate science, we can also say that the presence of many possibilities represents the internal variability of the mechanical process. Although there are still some filaments in panel (d), the majority of points accumulates in two patches about (x, v) = (±1, 0). Here, simple motion about fixed points starts to dominate. This process becomes completed by the instant of the last panel, already on the plateau of zero driving, indicating that the snapshot attractor has practically become split into two fixed point attractors. The scatter of the blue points around the two centers is limited; internal variability is diminishing as chaos is in the process of dying out.
The snapshot points of the given instant are also displayed in Fig. 3. It is remarkable that the snapshot hyperbolic point x * H,n marked with red always falls on the snapshot chaotic attractor. Even if other snapshot hyperbolic points are not known, it is reassuring to see that at least the one whose analytic form was found is  Fig. 3. Here, we display the steps n = 10, 20, 40, 60, respectively in (a-d). In a-c, the yellow unstable manifold was originated at n = 5 in the form of a section of length dl = 0.2 along the diagonal crossing x * H.5 , consisting of 10 5 points, while in d it was generated as usual: by iterating a short line segment crossing the origin. The SHP of (5) (red) and the SNPs of (7) (orange) are displayed as well. (Color figure online) part of the attractor. The snapshot nodal points x * N ,n , plotted in orange, are initially not part of the snapshot attractor, but as time proceeds, the snapshot attractor seems to converge toward these simple attractor points (panels (d), (e)). On the plateau of zero driving, x * N corresponds to the minima of the double-well potential at x = ±1.

Snapshot attractors and unstable manifolds of the snapshot hyperbolic point
Usual chaotic attractors are known to be the closures of the unstable manifolds of all the hyperbolic cycles, each of which lie densely on the attractor [77]. As a consequence, even the unstable manifold of only one fixed point provides a good approximant hardly distinguishable by the naked eye from the attractor. To check whether an analog property holds for snapshot attractors, we generate the unstable manifold of the SHP x * H,5 . We do this according to the dynamics (2), (3) by following the images of a short line segment (of length dl = 0.2) lying along the diagonal (corresponding to the direction of the unstable manifold in the unperturbed case) at n = 5. This line segment becomes stretched and folded, and after n − 5 periods it provides an approximant of the unstable manifold at time n. By this time, the snapshot hyperbolic point has moved into x * H,n , but is still found to lie on the snapshot attractor. In Fig. 4, we overlaid this manifold plotted in yellow for n = 10, 20, and 40 with the snapshot attractor (blue dots) of these instances. The agreement is quite con-vincing, and we can conclude that the snapshot attractor at a given instant is approximated by the unstable manifold of that instant. This unstable manifold can be called the snapshot unstable manifold of the SHP x * H,n ; although this is the unstable manifold originating from x * H,5 , it evolves together with the moving SHP x * H,n . The last panel exhibits in yellow the unstable manifold of the origin of the undriven system. The spiraling pattern of the blue snapshot attractor follows the unstable manifold of the origin after arriving at the plateau of zero driving.
We can also obtain snapshot stable manifolds the same way, only now starting from a later time instant than x * H,n from a line segment along the sub-diagonal (the direction of the stable manifold in the unperturbed system), which we then subject to iteration backwards in time according to (2), (3) to arrive at the instant n. The resulting pattern is then said to be the snapshot stable manifold belonging to the SHP x * H,n at instant n.

Time-dependent basin boundaries and snapshot stable manifolds
Figure 5a-c shows the basins of the long-term fixed point attractors (±1, 0) of the plateau using color coding. Same colors indicate the initial conditions that lead to the same attractor. The images show the basins for different starting instants, using the same scenario; that is, they capture different parts of the ramp. We see, again, that approaching the end of the ramp chaos seems to weaken, the boundary between the basins becomes smoother.
In Fig. 5d-f, we compare the basin boundaries with the snapshot stable manifolds of the subsequent time instants. The stable manifold belonging to time instant n is initiated at period n 0 = 50 on a short interval of length dl along the sub-diagonal about the origin. This point is taken because the basin patterns are determined after a long time spent on the plateau, and the hyperbolic point of the undriven case is (0, 0). This short segment is then subjected to backward iteration of length 50 − n. The green lines, indicating points of the instantaneous stable manifold, overlap very well with the basin boundaries. 2

Snapshot horseshoes
Smale horseshoes (or chaotic saddles) are generally thought to be the most essential components of both 2 At the end of Sect. 2.4, the method for creating snapshot stable manifolds was laid out. However, because the large eigenvalues of the SHP (4), points around it are pulled apart very fast (the stretching rate is e 2π ≈ 500 per period), for longer scenarios this has to be augmented with another method in the numerical implementation. That is, one has to initiate new manifolds from the SHP every few steps of the iteration, which are then thought to complement the existing snapshot stable manifold. In our particular cases, we initiate manifolds every 2 periods of the scenario from segments of length dl = 0.02 crossing through the corresponding SHPs. This method was used in Figs. 5 and 6c. permanent and transient chaos [66,77]. The literature on cases of arbitrary time dependence is sparse, but such objects were identified in random maps [55], in time-continuous systems subjected to drivings of changing strength in both dissipative [45,46] and Hamiltonian systems [38]. A system where changing tangles of manifolds can be observed are fluid flows of aperiodic time dependence, as illustrated in a turbulence-related experiment in [59]. Most recently, observational evidence of snapshot chaotic saddles has been reported in solar plasmas [11].
In the horseshoe structure, points of the chaotic saddle are the intersection points of the stable and unstable manifolds of the same hyperbolic point. In traditional cases, pieces of the unstable and stable manifolds can be obtained by iterating small intervals taken along the diagonals of the phase space crossing through an x * H forward and backward, respectively, up to the same number of periods. However, when applying the same method to systems where parameter drift is present, the obtained manifolds would belong to different values of the parameter, in our case, the driving amplitude.
Furthermore, this implies that in systems subjected to parameter drift one should take into account that the hyperbolic point moves in time. The main issue is that the intersection points of the manifolds must belong to the same instant, which for moving hyperbolic points (i.e., SHPs) implies that the points from which the manifolds are initiated must belong to different instances [38,45,46,55]. Observe that at the beginning of the scenario the area where the two manifolds cross is of considerable extension, while near the end of the ramp it shrinks to a few points about the origin. The SHP and the SNPs of (5) and (7), respectively are also displayed (red and orange). (Color figure online) We thus require the intersection points to belong to the same time instant, say n. To this end, we consider the stable manifold initiated at period n +ñ on a short interval along the sub diagonal about x * H,n+ñ and subject it to backward iteration of lengthñ. Analogously, the unstable manifold should be initiated at time n −ñ on a short interval along the diagonal about x * H,n−ñ and be subjected to forward iteration of lengthñ. By overlaying these manifolds, a snapshot horseshoe structure shows up belonging to the time instant n. The use of the term is supported by the fact that for a different time instant n, the visual appearance of the horseshoe is different. Note, that the choice ofñ is not set by any internal constraint of the system. Larger values would result in a larger set of intersection points at the same n, which would, however, contain the ones obtained using a smaller value.
In Fig. 6, we display the snapshot horseshoe structure at different instances. Observe that at the beginning of the scenario (ε = 0.08, n = 10 in panel a) the area where the two manifolds cross is of considerable size, while in the middle of the ramp (panel b) it is smaller, and it all but disappears by ε = 0.02, n = 40 where the two manifolds intersect each other in only a couple of points.
Similarly to how time-independent horseshoes determine the skeleton of chaos in the stationary system, we can say that a snapshot horseshoe signals the spread of chaos in a system subjected to parameter drift. Thus, from the findings of Fig. 6 we can conclude that the strength of chaos decreases during the scenario, to the point where eventually the saddle collapses to a few points. This vanishing of chaos warrants a further, quantifiable investigation, namely by looking at the dynamical instability of the system.

Strength of chaos
There have been different quantities suggested to quantify chaos in systems subjected to parameter drifts. Besides a topological-entropy-like quantity (the expansion entropy) [36], the extension or the standard variation of the snapshot attractor was also suggested [45], as well as the mean normalized distance over a whole time integration [69], and a cross-correlation-based method [70]. As a particularly simple concept, we use here the one suggested in [38] which enables us to derive a quantity for describing chaos in the system, which we shall call the instantaneous Lyapunov exponent. This term has been used in a recent paper [65], too, but the concept differs from the one used here since, e.g., no ensemble average is taken in [65].
The dynamical instability of chaotic systems is typically characterized by the average growth rate of distances between pairs of points [66]. Here, we evaluate the average over the ensemble of the snapshot attractor: where r (t) is the distance of a pair of points at time t which were close to each other initially in the ensemble, and the bracket denotes the average taken over the ensemble.
The results for the ensemble-averaged pairwise distance (EAPD) in the scenario used up to now are shown by the top curve in Fig. 7a. For every 10th point of the ensemble a pair of initial distance r 0 = √ 2 · 10 −10 was assigned (adding 10 −10 in both directions). For the top curve, a red line of slope λ = 0.50 ± 0.01 is fitted to the initial, exponential section indicating chaotic behavior in the first part of the scenario. However, from n ≈ 15, the distance growth starts to slow down until it eventually stops at n ≈ 30. Since the slope of the curve can be considered as an instantaneous Lyapunov exponent, which vanishes here, we conclude that chaos ceases in the system after 30 periods. After this, a slow decrease follows indicating a chaos-free/regular approach toward the ultimate fixed point attractors. Observe that the peak of the curve is still at a considerably small value meaning that, on average, a large number of pairs have remained close during the scenario.
The subsequent curves in panel (a) correspond to the results obtained with ensembles initiated at later times in the scenario. In these cases, the initial shape of the ensemble is the snapshot attractor belonging to starting time n, where the driving amplitude felt by this ensemble is ε n = ε 0 − αnT . The initial slopes (red lines) indicate that the Lyapunov exponent at the starting time decreases, the attractors are less chaotic when feeling a later section of the original scenario. Each curve is bent, and remarkably they all reach a maximum by about n = 30, indicating that the instant when chaos ceases to exist is determined by the system and the scenario, while the time when the observation starts does not influence this. The driving amplitude belonging to the instant when all the curves reach their maxima, 30 periods, is ε = 0.04, which is relatively close to 0.025, the value below which no transient chaos is found in the bifurcation diagram of Fig. 2. This reinforces the observation of [45] according to which snapshot attractors become non-chaotic approximately in the range where not even transient chaos exists in the system without parameter drift. To answer the question of whether the same maxima of the curves is characteristic of this specific scenario or the consequence of some intrinsic trait of the system, we repeat the simulation with the same scenario except changing the damping to β = 0.02. The result can be observed in Fig. 7b. As we can see, the curves again have about the same maxima, the location of which however is now changed to n ≈ 20 where the driving amplitude is ε = 0.06. Note that with the different damping coefficient, the produced chaotic saddles of Fig. 2 would be different as well. From this, we can conclude that in decaying systems such as this, the EAPD curves initiated in the same scenario have their maxima at the same time. The time instant itself is set by the damping coefficient in a given scenario.

The model
We consider here the dissipation-free (β = 0) limit of (2): where ε(t) is given by (3). This system was investigated, e.g., in [29] with α = 0, and in [38] with Note that no structured pattern emerges in this extended ensemble. When following this scenario, we are scanning through the stationary cases of all the panels of Fig. 8 α = 0. System (9) corresponds to the Hamiltonian H ( p, x, t) = p 2 /2 − x 2 /2 + x 4 /4 − xε(t) cos ωt. The total energy is not constant due to the driving; instead, in the presence of this ramped scenario it decreases on average. Without a drift (α = 0), the system possesses a typical divided phase space consisting of a mixture of quasiperiodic tori and chaotic seas. The overall outlook of the phase portraits depends on the driving amplitude ε 0 ; as its value decreases, the chaotic area decreases in size, as Fig. 8 illustrates.

The choice of ensembles
When starting with an extended ensemble, as typically done in dissipative cases (see previous section), one obtains an unstructured pattern, as seen in Fig. 9. The reason for this lack of structures is that the basic components of the phase space, tori, and chaotic seas both become time dependent, as shown in [38], and the extended ensemble is too crude to separate these components. It is clear, however, at this stage already that members of the ensemble are scattered in an extended region for long times, i.e., a strong internal variability is present in the drifting dynamics. This property warrants a refined application of the snapshot view of parallel dynamical evolutions.
It is therefore worth selecting a more specific ensemble, concentrated on one of the basic components of the phase space [38]. Being one-dimensional objects and simpler than chaotic seas, let us take the tori of system (9) with α = 0, i.e., the tori of the stationary phase space, and evolve them under (9) with α = 0. To illustrate a basic effect, we select 4 tori, lying in different regions of the original phase space (light green curves of Fig.10). A comparison with the first panel of Fig. 9 reveals how much these initial ensembles differ from an extended one. These ensembles are followed for a short period of time. As Fig. 10 illustrates, they all remain closed curves by the end of the scenario, we thus see that tori might survive the parameter drift. Because they change in time, the term snapshot tori was coined for them in [38]. The fact that tori remain closed curves shows that we identify by them ensembles of limited internal variability. It was also revealed in [38] that as time goes on, deformed snapshot tori can break up, which is signaled by a sudden stretching in their structure and is completed when the tori become completely embedded into the snapshot chaotic sea. However, this phenomenon cannot be observed in this system, owing to the decreasing trend in the driving amplitude.

Snapshot hyperbolic and elliptic points
The snapshot hyperbolic and nodal points remain well defined in the limit of no dissipation (β = 0). The hyperbolic ones (4) simplify to (10) and v * H =ẋ * H (t), while on the stroboscopic map we find The nodes become converted into elliptic ones, with eigenvalues ±i √ 2. These snapshot elliptic points (SEPs) are, from (6), given by: and v * E =ẋ * E (t), while on the map they are represented by

Evolution of the snapshot phase portrait
In the presence of a parameter drift the evolution of both basic components, tori and chaotic seas are interesting to follow. Therefore, in this section not only tori are chosen as initial ensembles, but the chaotic sea of the initial phase space as well. The latter is an extended set whose shape is, however, defined by the original dynamics, and not artificially chosen as in the first panel of Fig. 9. The instantaneous pattern of the phase space at later times can be called a snapshot phase portrait. In Fig. 11, we display such a case. Panel (a) shows the stationary phase space at ε 0 = 0.1 in which (as in Fig. 8a), besides a few tori, an extended chaotic sea (scattered blue dots) can also be seen. Panel (b) belongs to the very first stroboscopic instant of the ramp when the driving amplitude is 0.08, while panel (c) is the penultimate instance with amplitude 0.02. Observe the deformation of the initial tori, similar to the phenomenon of Fig. 10, and that the dotted region also takes a different form. Because of the extended appearance of the blue dots, such snapshot chaotic seas represent ensembles with pronounced internal variability. The elliptic fixed points now become the SEPs whose coordinates are described by (13) and are plotted in orange. Panel (d) represents the instance of the n 0 = 5th iterate. The driving has now reached zero absolute amplitude and it will be kept that way from here on, meaning that in the following steps (9) will be applied without any driving (ε = 0). Panel (e) shows the evolution of the phase space pattern at n = 15 (10 iterations spent on the plateau). Observe that the snapshot tori started a transition into thin lines orbiting around the two elliptic fixed points x * = ±1, v * = 0. A folding of the tori into approximate orbits of the undriven problem (panel f) can be observed.
The evolution of the set of blue dots has to be interpreted more delicately. They are certainly the images of the initial chaotic sea visible in panel (a) according to the scenario. However, it is unlikely for all of the blue dots to designate regions of chaotic nature since the undriven 2D dynamics at the end does not support chaos at all. On the ramp, e.g., in panels (b), (c) a portion of the blue dots is certainly chaotic, and while the extension and shape of this chaotic region cannot be determined by naked eye, one anticipates that its size decreases with time. The lack of chaos at the end of the ramp (panel e) is also visually clear by observing that these dots start to trace out the well-known trajectories of the double-well system without driving (panel f). From the extended size of the points originating in a chaotic sea, we have to conclude that pronounced internal variability is not necessarily a sign of chaos in Hamiltonian problems. Initial chaos might lead to the appearance of a large number of permitted states in the phase space even after chaos has died out. The last panel exhibits typical phase space trajectories of the undriven system to demonstrate that the stirring pattern visible in panel (e) has started to follow the flow of zero driving. In Sect. 3.5, we show that the size of the instantaneous chaotic sea can be estimated by constructing a snapshot horseshoe belonging to the given instant.

Hamiltonian snapshot horseshoes
The principles of the horseshoe construction are similar to that of the dissipative case. In our Hamiltonian problem, we takeñ = 5 and construct the horseshoe structure for n = 10 and n = 40 in the scenario ε 0 = 0.1, n 0 = 50 as Fig. 12 illustrates. Note that while being further down on the ramp, the snapshot horseshoe structure in Fig. 12b is smaller compared to that in Fig. 12a. Since the horseshoe structure is the approximant of the area where chaos is present, we can conclude that as the scenario sweeps through phase spaces containing less and less chaos (compare Fig. 8), the size of the snapshot horseshoe, and thus the area providing snapshot chaos, shrinks. For this reason snapshot chaotic seas, for example the one seen in Fig. 11a-d, albeit in a 10 times faster scenario, lose their chaotic characteristics for some of their points until the end of the ramp. Figure 12 also provides an explanation as to why tori cannot, or very hardly break up in this system of decreasing amplitude. Two snapshot tori, initially lying practically on the edge of the elliptic islands, are plotted on both images. Note that the snapshot elliptic points  45 for the stable (pink) and x * H, 35 for the unstable (green) manifold. Displayed in both cases are the snapshot tori of initial conditions x 0 = −0.4, 0.85, v 0 = 0 (blue), the SEPs of (13) (orange) and the SHP of (11) (red). Observe that on (a) the structure is larger than that on (b), leading to the conclusion that during the scenario the region where true chaos is present, shrinks. (Color figure online) (13) (orange dots) fall into the middle of the snapshot tori. As one can see, they hardly even deform because the horseshoe, and thus the snapshot chaos, shrinks, making it almost impossible for the tori to be mixed in with the chaotic sea.

Strength of chaos
The EAPD (8) is worth applying here, too, to characterize dynamical instability, just applied to point pairs at time t which were close to each other on either a torus or in a chaotic sea at t = 0, and the bracket in (8) in this case denotes the average taken over such a sub-ensemble. Figure 13 shows this quantity for the time evolution of the initially chaotic ensemble in the scenario ε 0 = 0.1, α ≈ 0.0003 which reaches the undriven case at iteration n 0 = 50. For every 10th point of the ensemble, a pair of initial distance r 0 = √ 2·10 −10 was assigned (adding 10 −10 in both directions). Here, because we are in a chaotic sea, the exponential growth starts right away, in contrast to snapshot tori going through the break-up process (see Fig. 19 and [38]). The remarkable feature here is that while the curve is expected to level off, this happens at a significantly lower ρ value than naively anticipated. The initial width of the chaotic sea is about 3 units (compare Fig. 8a or 11a), thus one would expect the curve to level off at around ρ ≈ ln 3 ≈ 1, since the distance of points cannot grow further than the size of the chaotic sea. In contrast to this, Fig. 13 shows the leveling off at around − 5. This is in part because the second half of the scenario is spent on the plateau of the undriven problem (similar to the scenario of Fig. 11), where there is no chaos at all. Some slight growth can be observed here as well, which is due to the fact that point pairs on nearby tori rotate with slightly different speed as the tori possess different winding numbers, and hence, the pair separates approximately linearly in time.
Note however, that the flat part of the curve starts at around n = 30, well before the start of the plateau at n 0 = 50. The reason for this can be found in the observations of the previous section: although we chose a chaotic sea as our initial ensemble, because of the shrinking of the snapshot horseshoe structure, the area that is actually chaotic in nature (i.e., the region where points are strongly pulled apart) shrinks as well, resulting in less and less pairs being able to separate. We conclude that chaos disappears in the average behavior at about the 30th iterate.
In order to gain more insight, in Fig. 14 we show the point pairs that are followed in the scenario. The ones that are not further apart than r c = 10 −2 are deemed "close" (generally these cannot be distinguished by eyesight) and colored purple. For pairs where r > r c one member is colored red, the other is colored blue and they are connected by a black line. The number of "not close" point pairs is also displayed. The first image of n = 20 is from the exponential part; ten steps later (second image) the pairs with large distance is more than doubled. However, after this instant the number of such pairs is practically constant.

Magnetic pendulum
The Duffing equation, as described above, cannot provide chaotic behavior without driving, because it has only a 2D phase space. To understand how chaos can appear in undriven, damped systems we consider the example of the well-known magnetic pendulum, whose motion can be described in a 4D phase space. The comparison with the previous, driven cases is indirect: it is the energy of the system which changes in time. Its instantaneous value can be considered to be the analog of the instantaneous driving amplitude of the Duffing system. The three-magnet problem is widely used as an example of chaos or transient chaos (see, e.g., [10,61,68,77]). However, in [61], it was shown that this system does not exhibit usual transient chaos since the character of transients changes in time; hence, the term doubly transient chaos was suggested to characterize such cases [10,[61][62][63]. Recently, additional features of doubly transient chaos have been explored [44], and we are applying these observations here to a system particularly well suited to the spirit of this feature article.
To maintain the bistable character of the Duffing equation, we use two identical magnets for our pendulum. An iron bob, attracted by both magnets, is fixed to a rod above the horizontal plane of the magnets. The bob is subject to the influence of gravity, attractive magnetic forces, and drag due to air friction. For simplicity, we further assume that the length of the pendulum rod is long compared to the distance between the magnets, which allows us to describe the dynamics of the bob entirely in terms of the (x, y)-coordinates of the plane and benefit from a small angle approximation: we can approximate the swinging of the pendulum in the x, y plane as that of a harmonic oscillator of some natural frequency ω 0 . The plane of the magnets is assumed to lie at distance d below this plane.
Following Refs. [68,77], we assume, for simplicity, an inverse-square law interaction between the bob and the magnets as if they were point magnetic charges. The resulting dimensionless equations of motion arë where (x i , y i ) = (±1, 0) are the coordinates of the ith magnet (i = 1, 2). The time unit is chosen so that the coefficient of the magnetic force is unity. In (14) and (15), ω 0 is the natural frequency, and β is the damping coefficient; here, D i (x, y) = (x i − x) 2 + y 2 + d 2 is the distances from the pendulum bob to the ith magnet. In our simulations, we set ω 0 = 0.3, and d = 0.3. The locations (x * ± , 0) of the attractors in the plane of the swinging is found to lie very close to the location (±1, 0) of the magnets since d is small compared to unity. The deviation turns out to be proportional to d 3 , and in leading order we find x * ± = (1 − d 3 (ω 2 0 + 1/4)). In a broad range of parameters, the origin of (14) and (15) is a hyperbolic fixed point. The linearized dynamics about the origin decouples in x and y. The resulting conservative force is repulsive in x:ẍ = a 2 x where a 2 = (4 − 2d 2 )(1 + d 2 ) −5/2 − ω 2 0 > 0 in the range of interest. In the y direction, this resulting force is attractive,ÿ = −ω 2 x where ω 2 = 2(1 + d 2 ) −5/2 + ω 2 0 . The eigenvalues of the dissipative linearized dynamics are therefore One of the eigenvalues is real and positive, one is real negative, and two are complex with a negative real part (in underdamped (β < ω) cases we are going to investigate). The hyperbolic point possesses thus a one-dimensional unstable manifold and a threedimensional stable manifold.

Time-dependent basin structures
The dynamics of magnetic pendula is traditionally displayed by showing colored images of basins of attraction. First, therefore, we show the basin of attraction of the two stable fixed points x * ± for a relatively large dissipation, β = 0.05 in Fig. 15. The different colors mark the initial locations associated with the attractor which the motion converges to. The trajectories were launched at a fixed initial energy E = 4.0. We see an intricate dependence of the basins on the initial conditions. Note that in the usual presentation of basins of attractions, the particles are released from rest; hence, the energy is not constant over the configuration space. We choose here the constant energy representation since the monotonous decrease of the energy is a central feature of doubly transient chaos, and it is worth monitoring the energy dependence of the basin structure. At any fixed energy E, the permitted initial states are restricted to a region bounded by an equipotential line of the potential energy arising from the pendulum and the magnets 15 Basins of attraction for the two stable fixed points in the x, y plane at a fixed initial energy E = 4, β = 0.05. White dots indicate the locations of the fixed point attractors at x * ± . The round shape of the accessible region of the phase space follows from the fact that the equipotential lines are nearly circular for high enough energy where constant C is chosen so that the potential energy is zero in the fixed points. The basin structure visible in Fig. 15 belongs to E = 4.0, and the initial velocities are chosen to beẋ = 0,ẏ > 0, with the magnitude of the vertical velocity computed from the set initial energy.
In Fig. 16, we see that the initial energy of the pendulum at launch has a crucial effect on the basin boundaries. For small energy, the pendulum is already in the neighborhood of one of the fixed points, and hence, no fractality occurs. For large energy, however, there is enough time for complicated motion during the decay of the energy, and hence, fractality is present. The dissipation dependence of the basins, at fixed energy, is discussed in Appendix B.
In traditional cases, fractal basin boundaries imply the existence of a chaotic saddle whose stable manifold is part of the basin boundary [55,77]. Combining this with the energy-dependent shape of the boundaries, as illustrated by Fig. 16, the following argument can be applied. Although the set of points in a plane of initial conditions belonging to a given energy will not have the same energy as time goes on, the average energy of the set will monotonously decrease. The change with energy thus resembles the time-evolution of the dissipative system (and of the driving ampli-  Fig. 16 with a temporally changing structure of a certain stable manifold. One can thus speculate about the existence of a snapshot chaotic saddle in (14) and (15) which contains the hyperbolic fixed point of eigenvalues (16) at the origin at any instant. Ultimately, the snapshot chaotic saddle converges to the saddle point at the origin and loses its chaotic character.
The snapshot chaotic saddle lives, of course, in the four-dimensional phase space, and is difficult to visualize. (A three-dimensional cut at a given instant was constructed in the 3-magnet problem in [61].) Because of this difficulty, we turn here to a complementary view, and choose an appropriate form of an ensemble to illustrate, even quantitatively, the features sketched above.

Evolution of initially Hamiltonian tori
Without damping, β = 0, the motion can either be chaotic or regular, depending on the initial conditions. This is illustrated in Fig. 17, where the Poincaré section of this conservative system is shown, containing regular islands in an extended chaotic sea, the usual pattern of divided phase spaces. To unfold new facets of doubly transient chaos, we apply the snapshot view to the undriven dissipative case; we use ensembles to monitor the history of the structural changes of the phase space. In the lack of global attractors (other than the states of rest), and because an extended ensemble would only show the overall shrinking of phase space volume toward two points, we find it useful to start with the dissipation-free case, and define the points of any of the tori there as our ensembles. As has been done before for the driven Hamiltonian Duffing problem in Sect. 3, we select a set of KAM tori embedded within each other. Then, in this case we switch on dissipation, and follow the ensemble of any of the tori to see how they deform and whether and how they break up. In the lack of any driving period, we follow the dissipative dynamics on the same Poincare section as used for the zero dissipation case. The dynamics on this section (because we chooseẋ = 0) is represented by a three-dimensional map, which does not reduce to a planar map since the energy is not constant now. When observing the dynamics on the x, y plane, we see a projection from the full phase space onto the plane in which trajectories might intersect. In Fig. 18, a sequence of Poincaré sections is shown, whose points start from 6 KAM tori in the large island visible in the right-hand side of Fig. 17. The dissipation is relatively small here, β = 0.005, so that the change is relatively weak between two subsequent Poincaré sections.
It is remarkable that images of the original tori keep their closed curve character over several Poincaré steps. We see that the outermost torus starts to deform between Poincaré sections 1 and 2, and it practically disintegrates during the next two steps. In Poincare section 3 (panel d), the internal tori also start to deform before breaking up in step 4. The behavior is very similar to that observed in the driven Hamiltonian case with increasing ramp [38]. In terms of internal variability, the breakup of tori implies that ensembles of limited variability can evolve into ones with pronounced variability which implies chaotic-like behavior in this case as the next section also illustrates.

Strength of chaos
To further characterize the process of the breakup of tori, we measure the EAPD, see (8), for the third torus from outside in Fig. 18a. The time dependence of ρ is shown in Fig. 19 for two different values of the damping coefficient. For small dissipation, β = 0.005, we see that the distance between initially close points, after a short transient, increases exponentially until the distance reaches the size of the available part of the phase space that is of the order of unity, and corresponds to ρ ≈ 0. Then, the average distance starts to decay exponentially as all points relax toward one of the fixed points. Chaos thus ceases to exist in this ensemble at about the 100th iterate. For a large dissipation, β = 0.05, we see that relaxation sets in much faster, before the exponentially increasing range could result in an extension comparable to the distance between the fixed points. This lasts to about the 35th iteration after which chaos disappears. In Fig. 20, we follow a different torus (the outermost in Fig. 18a). However, here we also distinguish between such pairs that finally settle on the same or on different fixed points. We see in Fig. 20 that the two curves are initially very close and they reach the size of the available phase space after an exponential growth. Then, the distance between pairs of points settling to the same fixed point decays due to dissipation. However, this distance converges to the distance |2x * ± | between the two fixed points for those pairs settling on different attractors.
We also measure the ensemble average of the energy of trajectories initiated from the same torus as that used for Fig. 19. Again, we use two different values for the damping coefficient β. The result is shown in Fig. 21.
In this autonomous problem, the time dependence of the energy is the analog of the time dependence of the driving amplitude seen in previous sections. All these quantities decay in time and tend to zero ultimately.

Conclusions and discussions
We have explored methods that enable one to follow the dynamics of systems subjected to parameter drift: the use of the theory of parallel dynamical evolutions (or ensembles in practice) to monitor changes in the phase space, the construction of snapshot Smale horseshoes (chaotic saddles), and the evaluation of the ensembleaveraged pairwise distance to quantitatively characterize the strength of chaos.
The chosen type of scenario leads to a full disappearance of driving in cases (i,ii), which was taken in order to remain in analogy with the undriven dissipative case (iii). From a climate change perspective, this is an extreme scenario. For the real climate this would correspond to a process leading to the full disappearance of the equator-to-pole temperature contrast, which leads to a disappearance of the usual global circulation. A possible mechanism in the background could be a disappearance of the incoming solar radiation, i.e., the dying out of the Sun. One is of course free to stop the drift before the plateau of zero driving is reached, or incorporating a plateau at a nonzero level. A rich dynamics can be explored by allowing an increasing ramp in the driving, too [38].
In the context of dynamical systems, a very special property of our examples was the weakness of the driving. We took advantage of this in finding perturbative results for the analogues of periodic orbits, e.g., snapshot hyperbolic points that are moving in time. Their stable and unstable manifolds have been essential in the horseshoe construction. It might appear to be a limitation of the results that the identification of such orbits is unlikely with stronger drivings. In such cases, however, hyperbolic fixed points or cycle points of the stationary problem can be taken, and their manifolds can be used to approximate the stable and unstable foliations of the drifting problem, as pointed out in [39]. The condition of chaos is then the existence of a large number of intersection points between these foliations, opening a language for the generalization of the manifold-based results of the paper.
The methods presented here can support the construction of snapshot phase portraits in engineering applications, too. These range from ships experiencing wave induced time varying restoring moments [79], through multiple scale dynamics [9], to investigating the effects of a moving mass on structures [76], or the switching between different working regimes [23], or milling stability with time dependent parameters [17]. Tipping events are sudden changes accompanying parameter drifts, and our method is useful to apply in such cases (see, e.g., [47,72]), too, including the avoidance of instability via fast enough rates of the parameter drift [6,85]. In analogy with climate sensitivity, parameter sensitivity or model uncertainty can be defined in mechanical systems, too (see, e.g., [3,49]), concepts to which the ensemble approach proposed here might be natural to extend.
Finally, let us speculate about a possible generalization of the concept of doubly transient chaos. The original terminology stems from the fact that the asymptotic state is an attractor (the state of rest), a chaotic-like dynamics preceding this can then only be transient. But since the character of the transients changes as time goes on, chaos is weakening, the adjective double transient was chosen. Consider now a driven dissipative system, case (i), with a ramp of any sign ending in a plateau of fixed driving amplitude where the attractors existing on the plateau might even be chaotic. The preceding dynamics (on the ramp), if chaos-like, is of different character, because the character of the transients toward the ultimate attractor(s) on the plateau is changing in time due to the change of the driving. Therefore, this type of dynamics can also be called generalized doubly transient chaos.
More generally speaking, as the existence of a snapshot horseshoe (snapshot chaotic saddle) was demonstrated in both types of driven systems (i), (ii), and was made plausible in undriven dissipative cases (iii) too, we can say that the dynamics of any system exhibiting transient chaos, if subjected to parameter drift, can be considered to be examples of generalized doubly transient chaos. to the Erwin Schrödinger Institute, Vienna, and for its program "Mathematical aspects of geophysical flows" held in 2020, which contributed to the initialization of this work.
Funding Open access funding provided by Eötvös Loránd University.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of interest
The authors declare that they have no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.
As a unified treatment, let us consider the equation where ω 2 0 = 2 (ω 2 0 = −1) corresponds to the behavior about the equilibrium points (the behavior about the origin).
As a first step, let us determine a particular solution for this inhomogeneous differential equation. Assuming that this y P inherits the frequency of the driving and keeps the character of its time dependence, let us search for it in the form: y P (t) =A 1 (ε 0 − αt) cos ωt + A 2 (ε 0 − αt) sin ωt + B 1 cos ωt + B 2 sin ωt.
By rearranging (22) as we see that particular solution y P (t) describes a moving point which keeps its eigenvalues, and hence its stability character, for any time, as long as the linear approximation remains valid. Although initially we carried out a linearization in a phase space region, this final result corresponds to a linearization about a distinguished trajectory and indicates its stability property. We can thus say that y P (t) provides the location, while its derivative provides the velocity, of a snapshot point. For ω 2 0 = −1 one of the eigenvalues is positive, this corresponds to a hyperbolic case, and as given in (6). For vanishing damping, β = 0, this goes over into a snapshot elliptic point (SEP), given by (12). We have been able to identify snapshot trajectories in the Duffing system since the driving amplitude ε(t) is assumed to be small. The undriven stationary case (ε = 0) is that of the autonomous double-well problem in which a hyperbolic point (x,ẋ = 0), and a pair of stable points (x = 0,ẋ = ±1) is known to exist out of which snapshot trajectories are inherited in the parameter-dependent case. How analytic forms of snapshot trajectories can be identified in general cases (even in the Duffing system for not small ε(t)s) remains a real challenge.

Appendix B: Dissipation-dependence of the basins in the two-magnet problem
The stable fixed points' basin of attraction taken at a fixed energy depends on the dissipation, too, as  illustrates. We see that for large dissipation, the basin boundary is a relatively simple, non-fractal curve. As the dissipation decreases, the trajectories have more time to relax toward the stable fixed points. This results in more and more complex basins of attraction.