Slim Fractals: The Geometry of Doubly Transient Chaos

Traditional studies of chaos in conservative and driven dissipative systems have established a correspondence between sensitive dependence on initial conditions and fractal basin boundaries, but much less is known about the relation between geometry and dynamics in undriven dissipative systems. These systems can exhibit a prevalent form of complex dynamics, dubbed doubly transient chaos because not only typical trajectories but also the (otherwise invariant) chaotic saddles are transient. This property, along with a manifest lack of scale invariance, has hindered the study of the geometric properties of basin boundaries in these systems--most remarkably, the very question of whether they are fractal across all scales has yet to be answered. Here we derive a general dynamical condition that answers this question, which we use to demonstrate that the basin boundaries can indeed form a true fractal; in fact, they do so generically in a broad class of transiently chaotic undriven dissipative systems. Using physical examples, we demonstrate that the boundaries typically form a slim fractal, which we define as a set whose dimension at a given resolution decreases when the resolution is increased. To properly characterize such sets, we introduce the notion of equivalent dimension for quantifying their relation with sensitive dependence on initial conditions at all scales. We show that slim fractal boundaries can exhibit complex geometry even when they do not form a true fractal and fractal scaling is observed only above a certain length scale at each boundary point. Thus, our results reveal slim fractals as a geometrical hallmark of transient chaos in undriven dissipative systems.


I. INTRODUCTION
Physicists often relate chaos with fractal basin boundaries and sensitive dependence on initial conditions [1][2][3][4][5]. While the former is a geometrical concept and the latter is inherently dynamical, the correspondence between the two has been established for conservative systems and driven dissipative systems. For example, in driven dissipative systems, the geometry and dynamics of a chaotic attractor are explicitly related through the Kaplan-Yorke formula [6], which connects the information dimension of the attractor with its Lyapunov exponents. A generalization of this formula to chaotic saddles is the Kantz-Grassberger relation [7], which connects the information dimensions along unstable directions with the associated Lyapunov exponents and the overall rate of escape from the saddle. While some fundamental open problems remain subjects of active research (e.g., the properties and applications of transient chaos [8][9][10][11][12], as well as the robustness [13], the classification [14], and the very definition [15] of chaos), studies of chaos in such systems are relatively mature [16].
In contrast, much less is understood about the relation between dynamics and geometry in the large class of physical processes categorized as dissipative but undriven, in which energy dissipated is not balanced by energy injected into the system. Examples of such systems abound, including coalescing binary systems in astrophysics, interacting vortices in viscous flows, chemical reactions approaching equilibrium, and many forms of self-organization. It also includes various arcade games (e.g., pinball) and games of chance (e.g., coin flipping and dice throwing) as well as cue and throwing sports (e.g., billiards and bowling). Due to the monotonic decrease of energy to its minima in such systems, all trajectories in a compact phase space will eventually settle to one of the fixed points, and the fixed points are the only invariant sets. Yet, for a transient period of time the dynamics can be very complicated and demonstrate sensitive dependence on initial conditions. A recent paper by a collaboration involving one of us [17] studied the nature of the dynamics of such systems. It was demonstrated that these systems show fundamentally different properties when compared to driven dissipative systems. In particular, they exhibit doubly transient chaos: system trajectories transiently follow a chaotic saddle which is itself transient. Moreover, the fraction of unsettled trajectories follows a doubly exponential function of time, which corresponds to an exponential settling rate rather than the constant settling rate observed in driven dissipative systems. However, the geometry of the attraction basins has not been characterized, and has been generally perceived as a very hard problem to address because these systems do not enjoy scale invariance (i.e., the basin boundaries do not exhibit any form of self-similarity, not even statistically). While it is known [2,5,17] that the attraction basins are intertwined and appear fractal-like, the absence of invariant chaotic saddles suggests that the basin boundaries may be simple at sufficiently small scales. Hence, the question remains whether the boundaries are true fractals. If the boundaries are fractals, what leads to the fractality despite the lack of invariant chaotic saddles? If they are not fractals, is there a characteristic length scale for the system that defines the resolution at which the boundaries become simple? How can we quantify the sensitive dependence on initial conditions in terms of their geometry? What roles do the observation length scale and computational precision play in one's ability to measure and simulate the dynamics of the system?
In this article, we investigate the geometry of attraction basins to address the questions posed above. We derive the condition under which the boundaries form a true fractal set (i.e., successive magnifications of the boundaries reveal new structures at arbitrarily small scales) and have the Wada property [18] (i.e., any boundary point between two basins is also a boundary point between all basins) for a general class of undriven dissipative systems. We show that this condition is satisfied generically, indicating that true fractal basin boundaries and the associated sensitive dependence on initial conditions are not only possible but are in fact common. The boundaries can also form a finite-scale fractal, characterized at each point by a finite length scale above which the fractal property is observed and below which the boundaries are simple around that point. Through extensive, high-precision numerical simulations on physical examples-the dynamics of a roulette of different shapes-we show that this fractality length scale can be smaller than the resolution typically used in simulations, making such basin boundaries practically indistinguishable from true fractals. We also find that, as a function of phase-space position, the fractality length scale can vary across many orders of magnitude. A common feature shared by the observed fractal and finite-scale fractal basin boundaries is that (at a given phase-space position) the fractal dimension for a given length scale decreases with the decrease of that length scale. Since this property implies that the boundaries would appear to cover less space when observed at higher resolution, we call such sets slim fractals. For characterizing the complex geometry of such boundaries, the existing fractal dimensions are not adequate, whether they are defined asymptotically at zero length scale or defined at a given finite length scale. Thus, to capture the cumulative effect of fractal scaling across all scales, we define the notion of equivalent dimension based on the process of increasing the initial-state accuracy to reduce the final-state uncertainty.
In the following, we first introduce the class of systems we consider and derive the condition for the fractality of their basin boundaries (Sec. II). We then apply the condition to the roulette systems and numerically validate the results (Sec. III). This is followed by the introduction of the equivalent dimension and its application to the roulette systems (Sec. IV). We provide concluding remarks in the final section (Sec. V).

II. FRACTALITY CONDITION FOR BASIN BOUNDARIES
For concreteness, here we focus on the class of twodimensional potential systems with frictional dissipation having n stable equilibria symmetrically located around an unstable equilibrium and separated by "hills" in the potential function. The equations of motion for such a system arë where µ is the dissipation constant and U (x, y) is the potential function. The dynamics of this system can be regarded as a scattering process, in which a trajectory entering the neighborhood of the unstable equilibrium swings back and forth chaotically between the hills before approaching one of the stable equilibria. The dynamics is thus dominated by the shape of the potential in this scattering region near the unstable equilibrium, which we define to be the origin. Writing in polar coordinates, the shape of the potential function near the origin is determined by the leading term in the expansion if U (r, θ) is smooth with respect to r. The symmetry of the system implies that the coefficients are n-fold periodic functions: a i (θ + 2πj n ) = a i (θ) for each integer j. The coefficient a 2 (θ) additionally satisfies a 2 ( 2πj n ) ≤ 0 and a 2 ( 2πj n ) = 0 for each j, because the attracting equilibria can be assumed to be located along the lines θ = 2πj n without loss of generality. We establish that the fractality of the basin boundaries is determined by system trajectories that move down a hill in the potential and approach the neighborhood of the origin. Specifically, we show that the basin boundaries are: (1) fractal if all such trajectories pass through the neighborhood, and (2) not fractal if some of them can asymptotically approach the origin without passing through it. Case (1) includes the generic situations in which a 2 (θ)r 2 is the leading term in Eq. (2), coefficient a 2 (θ) takes both positive and negative values depending on θ [with positive a 2 (θ) in the direction of the hills], and the dissipation is sufficiently weak. Case (2) includes the non-generic situation in which a 2 (θ) is identically zero [thus making the leading term in Eq. (2) cubic or higher] and the leading coefficient a j (θ) takes both positive and negative values.
As an example for case (1), consider the potential [i.e., n = 3, a 2 (θ) = − cos 3θ, and no higher-order terms]. Although this potential makes Eq. (1) an open scattering system with no attractors, it can be regarded as an approximation of a system that has n = 3 attractors far away from the scattering region. There are three possible ways [denoted E 1 , E 2 , and E 3 ; see Fig. 1(a)] for a trajectory to exit the scattering region. We can show that, between any two trajectories starting on the vertical line segment labeled A in Fig. 1(a) with velocity zero and eventually leaving the region through two different exits, we can find another trajectory that goes to the third exit (see Appendix A). Such a situation is illustrated in Fig. 1(a) by the red and green trajectories starting near the color boundary on segment A, which turn around near curve B and exit the region through E 1 and E 2 , respectively. We indeed see that the orange trajectory starting between them turns around near B, passes the neighborhood of the origin, and exits through E 3 . Since the same situation can occur after an arbitrary number of oscillations between the hills (e.g., after bouncing off B once and reaching C), this translates to the following property of the basins on A: between any two segments of different colors, we can always find a segment of the third color. These geometrical properties are verified numerically by successive magnifications near a boundary point in  Fig. 1(b). We note that our argument for segment A (on which the initial velocity is zero) can be extended to an arbitrary line segment in the full four-dimensional phase space connecting points from different basins (see Appendix A). This implies that any cross section of the neighborhood of any boundary point has a similar Cantor-set structure and has the Wada property, establishing that the entire set of basin boundaries is fractal.
As an example for case (2), consider the potential [i.e., n = 3, a 2 (θ) = 0, a 3 (θ) = − cos 3θ, and no higherorder terms]. With this potential, Eq. (1) is also an open scattering system that approximates one with three attractors. In this case, we can show that there exists a finite-length line segment {−r s ≤ x ≤ 0, y = 0} from which all trajectories approach the origin asymptotically (see Appendix A) and that this segment is a simple boundary between the basins of E 2 and E 3 , which does not belong to the boundary of the basin of E 1 . This is because any trajectory starting above (below) this segment with zero initial velocity, no matter how close it is to the segment, moves toward the origin initially but soon curves away and exits through E 2 (E 3 ). The trajectories starting exactly on the segment do not exit the region at all. The green, orange, and black trajectories starting from A in Fig. 1(c) illustrate this situation. Thus, every point on this segment is a boundary point between basins of E 2 and E 3 , and hence is a non-Wada point, implying that successive magnifications around this segment would not reveal any finer structures. We can further show that the segment splits into two branches forming simple boundaries, each of which in turn splits into two branches forming simple boundaries [see branching points indicated by blue arrows in Fig. 1(c)], and so on, composing a binary tree of simple boundary segments. Thus, the boundaries are not fractal [as numerically verified by successive magnifications in Fig. 1(d)]; however, since they have a Cantor set structure down to finite length scales (which are different for different branches), we say that such boundaries form a finite-scale fractal. We now generalize this result to lift the zero initial velocity assumption. Our argument is based on applying the center manifold reduction [19] to the equilibrium at the origin. Transforming Eq. (1) with U 3 (r, θ) into a suitable coordinate system (x,ỹ,ũ,ṽ), we determine the local center manifold to be and the dynamics on that manifold to bė up to second order inx andỹ. These are both visualized in Eq. (6) to the global phase space, we establish that the full set of basin boundaries is a finite-scale fractal (see Appendix A for details).
It is interesting to note that the stable manifold of just one equilibrium (the origin) is responsible for the full complexity of the basin boundaries-whether they are fractal or finitescale fractal-for the class of systems we consider. To see this, note that the basin boundaries consist of all points from which the trajectories never leave the scattering region. Since the only possible asymptotic state in this region is the unstable equilibrium at the origin, any trajectory starting from a boundary point must approach the equilibrium. Conversely, any point from which the trajectory converges to the equilibrium is a boundary point. This is because one can always find an arbitrarily small change to the initial point that would make the trajectory steer left or right just before converging to the equilibrium, and eventually leave the scattering region through one exit or another. Thus, the set of boundary points is the stable manifold of the equilibrium.
In addition to case (2) discussed above, finite-scale fractals can arise when the origin is a local maximum of the potential [e.g., when a 2 (θ) < 0 for all θ], if the higher-order terms in Eq. (2) create unstable saddle points that play a role similar to that played by the origin in our argument above. We will see an example of this situation below. Also, the transition between fractal and finite-scale fractal boundaries can be studied using the class of potentials with arbitrary real parameter α. Indeed, we can fully characterize this fractality transition: the boundaries are fractal if α ≤ 2 and finite-scale fractal if α > 2 [see Appendix B for the analysis and Appendix C for numerical verification]. Finally, we note that our arguments above do not rely on the linearity of the dissipative term in Eq. (1) and can also be applied to systems with nonlinear dissipation (i.e., when µ is not constant and instead depends on the position, such as in electric circuits with nonlinear resistors [20] and in nanomechanical resonators [21]). In particular, our fractality condition based on the behavior of the trajectories approaching the origin remains valid for any nonnegative function µ = µ(r, θ), and the condition can be expressed in terms of µ(r, θ) (see Appendix B). For instance, if the dissipation is of the form µ = µ 0 r q , this condition reads as follows: the boundaries form a true fractal if α ≤ 2(1 + q) and a finite-scale fractal if α > 2(1 + q).

III. ROULETTE AS A MODEL SYSTEM
As a physical example that can be described using a potential of the form (2), consider a roulette system. When the game is played in reality, a ball is released to a spinning roulette with 38 slots labeled with different numbers. The ball collides multiple times with bumps on the surface of the roulette and eventually falls into one of the slots. In our study, we simplify this system by assuming that the roulette is still, has a smooth surface, and has three slots (thus n = 3). We consider three different shapes of the roulette surface, shown in Fig. 3 and given by the following functions: Note that these functions serve also as the (gravitational) potential of the system, and the three slots correspond to three fixed-point attractors A 1 , A 2 , and A 3 located at (r, θ) = (1, 0), (1, 2π 3 ), and (1, 4π 3 ), respectively. This means that the results established above apply to this system, implying that the basin boundaries are fractal for S 1 [for which a 2 (θ) = − cos 3θ takes both positive and negative values], while the boundaries are finite-scale fractal for S 2 [for which Three shapes of the roulette surface we consider, given by the functions S1, S2, and S3 defined in the text. In each panel, the white dot indicates the unstable fixed point at the origin, the green dots the attractors (A1, A2, and A3), and the red dots the saddle points away from the origin (only present for S3). The part of each surface corresponding to Si(r, θ) < 0.5 is shown in the bottom row. Surface colors indicate the value of the function, and a common color scheme is used in all six panels. a 2 (θ) = 0] and for S 3 [for which a 2 (θ) = −(2 + cos 3θ) < 0 for all θ, and the surface has three additional saddle points, as indicated by the red dots in Fig. 3]. To compensate for the fact that our simplified roulette is still, we consider initial conditions in which the ball is placed on the circle r = 2 and has a velocity tangent to the circle. Friction and gravity dominate the motion of the ball. In order to prevent the ball from moving too far from the center of the roulette, we impose a maximum v max (θ 0 ) on the initial speed v 0 , where v max (θ 0 ) is defined as the value of v 0 corresponding to zero centrifugal acceleration when the ball's initial position is (2, θ 0 ) in polar coordinates. The ball experiences a dragging force proportional to its velocity with coefficient µ [representing dissipation, as in Eq. (1)], and here we use µ = 0.2.

Figures 4(a)-4(c)
show that, despite the difference in the fractality resulting from the three shapes, the numerically estimated boundaries between the basins of the three attractors in the phase space show highly convoluted, fractal-like structures in all three cases. Comparing Figs. 4(a) and 4(b), we observe that the basin boundaries appear more complex for S 2 than for S 1 . However, a closer look at the structure around the points P 1 and show that the basin boundary is simple below a certain finite length scale (on the order of 10 −15 ). To systematically quantify this fractality length scale, consider applying the bisection algorithm to a small vertical line segment of length ∆ in the (v 0 , θ 0 )-space, which can be used to estimate the location of a boundary point (to a given numerical resolution).
We define (v 0 , θ 0 ) to be the length of the interval used in the last occurrence of the following situation in the bisection process: the midpoint belongs to a basin that differs from those to which the two end points belong. For example, the quadrupleprecision bisection procedure used to generate the magnification plots in Fig. 5 for P 1 and P 2 gives ≈ 2.58 × 10 −27 and = 1.42 × 10 −15 , respectively (with ∆ = 0.1 and resolution on the order of 10 −27 , see Appendix D for details). Note that the fractality length scale at P 2 is at the limit of doubleprecision calculation and thus could not be clearly resolved without using higher precision. This illustrates the fact that a finite-scale fractal can be numerically indistinguishable from true fractals. The fractality length scale can also be seen as a quantitative measure of the Wada property at a given point (see Ref. [22] for a different numerical approach to quantify this property).
The fractality length scale (v 0 , θ 0 ) can generally depend on phase-space location (v 0 , θ 0 ), and its spatial distribution is quite different for the three example shapes [see Figs. 4(d)-4(f)]. For S 1 , the computed length scale is at the chosen precision (= 10 −13 ) uniformly over the boundary set (although the exact number depends slightly on the details of each bisection sequence), which is consistent with the true fractality of the boundaries. For both S 2 and S 3 , the boundaries are finite-scale fractals, and for S 3 the length scale is indeed well above the scale of the chosen precision across the boundary set. In contrast, shows a mixed behavior for S 2 , where is close to the scale of the chosen precision for the most part, but is well above that scale in certain locations. In this sense, the finite-scale fractal for S 2 is closer to a true fractal than for S 3 . Further analysis of the probability distribution of , as well as of a quantitative measure of the Wada property, corroborates these observations (see Appendix E).  (v0, θ0). The red, green, and beige regions indicate the basins of the attractors A1, A2, and A3 (marked in Fig. 3), respectively. (d)-(f) Spatial distribution of the (color-coded) fractality length scale (v0, θ0) on the boundaries of the basins shown in (a)-(c). Note that v0 and θ0 are normalized by vmax and 2π/3, respectively, only for the axes of the plots and not for the computation of (v0, θ0). We compute (v0, θ0) using double precision and the bisection resolution of 10 −13 for each of the 1,024 × 1,024 grid points [corresponding to ∆ = 2 −10 · vmax(θ0), which ranges from 3.38×10 −3 to 6.48×10 −3 depending on θ0 and the roulette shape].
We expect to see similar geometry of the basin boundaries if we consider the more realistic case of a roulette rotating at a constant angular velocity with zero initial velocity for the ball. Rewriting Eq. (1) in the frame co-rotating with the roulette, we gain two additional terms representing the centrifugal and Coriolis forces. The former effectively adds a constant to the coefficient a 2 (θ) in Eq. (2), while the latter simply shifts the location of the basin boundaries without altering the fractality of the boundaries.

IV. EQUIVALENT DIMENSION FOR SLIM FRACTALS
The fractality of the basin boundaries can be quantitatively characterized also by their dimension, which can be defined through a scaling relation between initial-state accuracy and final-state uncertainty [1,23]. For a self-similar system and an N -dimensional region of its phase space, the scaling is f (ε) ∼ ε N −D , where the constant D is defined as the fractal dimension of the boundaries, and f (ε) is the final-state uncertainty, defined as the fraction of pairs of points belonging to different basins among all pairs that are within the region and apart from each other. In contrast, the scaling exponent is resolution dependent for the systems studied here [as shown in Figs. 6(a)-6(c) for the roulette system], which motivates us to adopt a finite-scale measure of the dimension. With that in mind, we first consider using the effective fractal dimension [24][25][26] given by which is a strictly local measure of how uncertainty changes with resolution. Specifically, the effective dimension describes the relation between small improvement in initial-state accuracy and the resulting reduction in final-state uncertainty at the finite scale ε. The usual (asymptotic) definition of frac- tal dimension is recovered in the limit ε → 0.
In general, for slim fractals-which we define as having D eff that decreases with decreasing ε-the effective dimension at a given scale fails to capture the complexity of the basin boundaries observed at larger scales and its impact on the dynamics. To see this, consider the case of finite-scale fractals, for which we have D eff = N − 1 below the fractality length scale δ > 0. In this case, the final-state uncertainty scales as f (ε) ∼ ε N −D = ε, which is the same as that of a system without sensitive dependence on initial conditions. This means that the improvement in the accuracy of initial conditions (i.e., the amount by which ε is reduced) required to achieve a given level of uncertainty can be much less compared to the case of fractal boundaries with N − D < 1. However, a prerequisite for benefiting from this linear scaling is that ε < δ, which is itself a requirement on the accuracy of initial conditions. A similar argument applies to the case of true (but slim) fractals, since benefiting from smaller D eff (thus larger scaling exponents) requires the initial condition accuracy to be high in the first place.
To characterize the finite-scale sensitive dependence on initial conditions, we define a new dimension D eq (ε) to be the dimension of an equivalent self-similar system, whose finalstate uncertainty is the same as the system being studied at two different scales: ε and a larger reference scale L. We term this quantity equivalent dimension and show that it can be expressed as which, as an integral quantity, properly accounts for the cumulative impact of the effective dimension on the relation between initial-state accuracy and final-state uncertainty in the systems we consider. The equivalent dimension in Eq. (12) can be derived as follows. First, writing the final-state uncertainty of the equivalent self-similar system asf (ε ) = C · (ε ) N −Deq , where C is a constant, we have f (L) = CL N −Deq and f (ε) = Cε N −Deq . Next, we eliminate C from these two equations and obtain D eq = N − ln f (L)−ln f (ε) ln L−ln ε . Since this can also be obtained by using Eq. (11) and rewriting Eq. (12), we see that the equivalent dimension is indeed given by Eq. (12). Thus, we have a more intuitive and direct definition of fractal dimension that considers the entire process of decreasing ε to improve the accuracy of predicting the final state.
For the case of finite-scale fractals, which have fractal dimension D = N − 1, the dependence of the equivalent dimension on ε is given by the general formula for ε < δ [which follows directly from Eqs. (11) and (12)]. When D eq (δ) > D, we see that D eq (ε) slowly (and continuously) decreases from D eq (δ) and approaches D as ε → 0. Thus, the equivalent dimension for scales below δ "feels" the effect of large D eq (δ) (and hence of D eff at scales larger than δ), which reflects the sensitivity to initial conditions observed at scales above δ. While we do not expect Eq. (13) to be followed exactly in practice, as the scaling of f (ε) is never perfect, we do expect D eq (ε) to start decreasing at the fractality length scale and approach the asymptotic dimension D = N − 1. This is indeed observed in Fig. 6(d) for onedimensional cross sections (thus D = N − 1 = 0) of the basin boundaries in our roulette system with shapes S 2 and S 3 . For S 1 , with the basin boundaries forming a true fractal, the equivalent dimension seems to approach D eq ≈ 0.14. The uncertainty-based calculations for all three cases are consistent with the results from another numerical approach (valid for N = 1), based on the fractal dimension estimate, where i is the length of the ith interval identified as part of the Cantor set structure of the basin boundaries (see Appendix F for details, where we account for intervals as small as i = 1.1 × 10 −27 ). Interestingly, Fig. 6(d) shows that the equivalent dimension of the finite-scale fractal for S 2 is significantly larger than that of the true fractal for S 1 for scales above 10 −20 . This, however, is actually consistent with the more complex basin boundaries observed for S 2 [ Fig. 4(b)] than for S 1 [ Fig. 4(a)].
The equivalent dimension fills a gap between classes of systems that can be suitably characterized with existing definitions. For self-similar systems, D eff is constant, as illustrated in Fig. 6(e), which corresponds to a straight line for the graph of ln f (ε) vs. ln(1/ε), as illustrated in Fig. 6(f). In this case, the complexity of the basin boundaries is captured well by the usual asymptotic definition of fractal dimension D (and by D eff at any finite ε). For non-hyperbolic systems (such as Hamiltonian systems with mixed phase space [26]), D eff increases as a function of ln(1/ε) [27], as shown in Fig. 6(e), and this corresponds to a convex curve in Fig. 6(f). In this case, the asymptotic dimension D reflects the complex geometry of the basin boundaries [and is lower bounded by D eff (ε) for finite ε]. In contrast, in the class of undriven dissipative systems we consider here, D eff decreases as a function of ln(1/ε), as shown in Fig. 6(e) [which is directly associated with the decrease of D eq as a function of 1/ε observed in Fig. 6(d)], and this corresponds to the concave curve in Fig. 6(f). This behavior of D eff is the defining characteristic of slim fractals and reflects their structure, which appears sparser at smaller length scales. Since D eff (ε) ≥ D in this case, D is only a "lower bound" for the finite-scale geometrical complexity reflected in D eff (ε), and can in fact indicate no complexity at all (e.g., the case of finite-scale fractals with asymptotic dimension D = N − 1, which equals the dimension of simple boundaries). Figure 6(f) illustrates that the shape of the graph of ln f (ε) vs. ln(1/ε) determines the initial condition accuracy required to achieve a given level of uncertainty f (ε) = f * . The concavity of this graph for slim fractals implies that the required initial condition accuracy ε * SF can be orders of magnitude smaller than the corresponding numbers for the other types of fractals, even when the asymptotic dimension [and thus the asymptotic slope of the curves in Fig. 6(f)] is the same. By design, D eq integrates the finite-scale complexity over a range of different scales, and is therefore suitable for studying such systems. As an integral of D eff , the equivalent dimension also enjoys the benefit of having less numerical errors than D eff .

V. DISCUSSION
We have demonstrated that the basin boundaries in systems exhibiting doubly transient chaos are generically true fractals, with both Cantor set structure and the Wada property observed at arbitrarily small length scales. It is instructive to compare this with the most previously studied forms of transient chaos (i.e., those in driven or conservative systems). In all cases, the basin boundaries correspond to the stable manifolds of an unstable invariant set. However, this set consists of an uncountable number of trajectories in previous cases but of only one unstable fixed point (the origin) in the systems considered here. Accordingly, the basin boundaries consist of one or a few manifolds in our case, as opposed to a bundle of uncountably many manifolds as in previously studied cases. But can a finite number of manifolds really define a fractal? The answer has long been known to be yes; the Koch snowflake is an immediate example-though the curve is non-differentiable and constructed ad hoc-but there are also known examples of a dynamically generated manifold forming a fractal, such as the invariant manifolds in homoclinic tangles [19]. Therefore, our result that such boundaries are true fractals is not the first demonstration of fractal geometry arising from a finite number of manifolds. However, an interesting aspect of the fundamental problem studied here is that, contrary to the case of homoclinic tangles, which embed Smale horseshoes with (permanent) chaotic trajectories, our dissipative systems cannot exhibit any sustained oscillations (chaotic or otherwise): every system trajectory must converge to an equilibrium. This underlies the fact that the stable manifold of a single equilibrium is fully responsible for the complexity of the fractal basin boundaries in the systems we consider.
We have also demonstrated that, even when the boundaries do not form a true fractal, they can give rise to a form of sensitive dependence on initial conditions, which nevertheless is not properly characterized by existing notions of dimension. These results challenge us to think differently about the definition of fractals. In many natural systems, geometric structures similar to fractals are observed, but they disappear at sufficiently small scales due to finite resolution or the nature of physics at that length scale. Nonetheless, those systems are likely to exhibit sensitive dependence on initial conditions at physically relevant length scales (e.g., those relevant for measuring the initial state). For example, games of chance, such as a dice roll, are undriven dissipative systems for which the basin boundaries can be simple at sufficiently high resolution [28,29], but there are no practical methods to measure initial conditions at that resolution and predict the outcomes. Moreover, our results show that the resolution below which boundaries become simple can be highly dependent on the phase-space location. An immediate option for studying such systems is to use an existing notion of scaledependent dimension, such as the effective dimension. However, for being a local measure of uncertainty versus length scale, the effective dimension alone cannot capture the physically observable sensitive dependence on initial conditions. Our integral-based definition of the equivalent dimension addresses this issue and, together with the fractality length scale, offers an analysis framework for studying undriven dissipative systems.
Our findings have profound implications for the physics of undriven dissipative systems. Prominent examples include the following: 1. Astrophysical systems. When two compact objects-e.g., neutron stars, white dwarf stars, or black holes-orbit each other emitting gravitational waves, we have an undriven dissipative system (since energy is lost due to gravitational radiation) [30]. Such coalescing binary systems serve as candidate sources of detectable gravitational waves. Characterizing the dynamics and geometry of these systems has been controversial, with arguments both for [31] and against [32] the existence of chaos and fractal basin boundaries. This issue is significant because sensitive dependence on initial conditions would lead to an explosion in the number of possible theoretical templates of gravitational waves against which the observational data would have to be matched, necessitating alternative detection methods.
2. Fluid systems. Interacting vortices in an otherwise still viscous fluid form undriven dissipative systems whose characterization of chaos is relevant and to which existing tools do not apply. Typically, scenarios involving three or more vortices are considered to allow for chaotic dynamics. In part because of the lack of adequate tools, previous studies of chaotic dynamics in such systems focused primarily on potential flows and other solutions of the Euler equations (in which dissipation due to viscosity is neglected) [33].
Our results established here open the possibility of a selfconsistent study of chaos in solutions that properly account for viscous dissipation.
3. Chemical systems. Nonlinear chemical reactions in thermodynamically closed systems can exhibit chaotic dynamics in the absence of any driving [34,35]. Previous studies of such systems, of which the Belousov-Zhabotinsky reaction is an example, have focused primarily on the far-fromequilibrium regime of strong chaotic oscillations. This regime is nevertheless transient, as dissipation unavoidably makes the system approach thermodynamic equilibrium. Our results can allow the complete characterization of this transition to equilibrium, which thus far could be only partially understood using the tools of conservative and driven dissipative systems.
Ultimately we note that our derivation of the fractality condition and the measures introduced here to quantify slim fractals do not rely on the specifics of the systems considered. Thus, we expect these results to be generalizable to undriven dissipative systems exhibiting doubly transient chaos in higher dimensions and with an arbitrary number of basins.
branches forming simple boundaries, consider a vertical segment at x < −r s , such as segment B at x = −0.5 in Fig. 1(c). Because x < −r s , there is a part of this segment from which trajectories eventually exit through E 1 , but the boundaries between different basins are simple. This is due to the presence of a trajectory that is deflected by the hill at θ = π/3 before approaching the origin as t → ∞ (black curve), similarly to the one starting on segment A and approaching the origin asymptotically. The same argument as above applied to this trajectory shows that the boundary between basins of E 1 and E 3 is simple. Thus, the simple segment of the boundary touching the origin splits into two branches forming simple boundaries (the blue arrow indicates the branching point at x = −r s , y = 0). Repeating this argument with segment C and other similar segments of initial positions, we see that the basin boundaries form a binary tree of simple segments. We observe that, as one moves away from the origin along the branching tree, the gaps between branches narrow, thus making the fractality length scale smaller. Since we have assumed zero initial velocities, the binary tree we just established is a two-dimensional cross section of the basin boundaries in the full four-dimensional phase space.
To see that this full set of boundaries is also not truly fractal, we apply the center manifold reduction [19] to the equilibrium at the origin. The Jacobian matrix at the origin has eigenvalues 0 and −µ, each with multiplicity 2. This implies that there exists a two-dimensional center manifold and a two-dimensional stable manifold in a neighborhood of the origin. Writing Eq. (1) in terms of the eigenvector coordinates (x,ỹ,ũ,ṽ) ≡ (x +ẋ/µ, y +ẏ/µ, −ẋ/µ, −ẏ/µ), we determine the center manifold and the dynamics on it to be given by Eqs. (5) and (6), respectively, up to second order iñ x andỹ. Figure 2(b) shows that the region is divided into three basins (corresponding to exits E 1 , E 2 , and E 3 ) by three segments of simple boundaries: the half lines θ = π/3, θ = π, and θ = 5π/3. Since the stable manifold is two dimensional, these boundaries on the center manifold extend to three pieces of simple, smooth, and thus non-fractal three-dimensional boundaries dividing a four-dimensional neighborhood of the origin. These boundaries, when extended as much as possible, intersect with the subspaceẋ =ẏ = 0 in the line segments 0 ≤ r ≤ r s , θ = π/3, π, 5π/3, in Fig. 1(c). The full set of basin boundaries can then be expressed as the set of all points whose trajectory ultimately falls on one of these local boundaries. This is because approaching the origin is the only asymptotic behavior possible for the system besides leaving the scattering region. Thus, in a sufficiently small neighborhood of any basin boundary point, the boundary is a threedimensional smooth manifold, since it is a pre-image of part of the local boundaries near the origin. Therefore, the global basin boundaries, whose two-dimensional cross section is the binary tree we established above, are not fractal but form a finite-scale fractal inheriting the branching structure. We show that the arguments in Appendix A are also valid for the class of potential functions U α (r, θ) ≡ −r α cos(3θ), which then implies that the basin boundaries form a fractal for α ≤ 2 and a finite-scale fractal for α > 2. In other words, as α decreases through the transition point α = 2, the basin boundaries transform from a branching tree structure to a shape similar to a Cantor fan that exhibits fine basin structures at any resolution. To see why α = 2 is the transition point between the two regimes, note that we can write U α (r, θ) = −β(r)r 2 cos(3θ), where β(r) ≡ r α−2 can be interpreted as an r-dependent prefactor for the quadratic potential U 2 . According to this interpretation, the dynamics on line θ = π/3, r ≥ 0 [governed byr + µṙ + 2β(r)r = 0] would be critically damped if ξ ≡ µ 2 √ 2β(r) = 1, over-damped if ξ > 1, and under-damped if ξ < 1. For α > 2, the same argument we used in the main text for the case α = 3 can be used to show that the basin boundaries form a finite-scale fractal, since the dynamics is also effectively over-damped in this case as long as r < µ 2 8 1 α−2 . In comparison, for α ≤ 2, there is a neighborhood of the origin in which the dynamics is effectively under-damped. We can thus use the same argument used above for the case α = 2 to establish the fractality of the basin boundaries. This fractality transition at α = 2 is numerically verified in Appendix C. In the more general case of nonlinear µ = µ(r, θ) ≥ 0, the condition for the boundaries to be a finite-scale fractal (true fractal) is that ξ(r) = µ(r, π/3) 2 √ 2β(r) > 1 (< 1) for all r sufficiently small. If the dissipation is of the form µ(r, θ) = µ 0 r q , q > 0, for example, the fractality transition occurs at α = 2(1 + q).   In each iteration, we determine the basin to which the midpoint of the interval belongs by integrating the system with quadruple precision and relative accuracy of 10 −4 (with respect to the length of the bisection interval for that iteration). We iterate until the interval length becomes equal to 2 −86 · ∆ ≈ 1.29 × 10 −27 and 2 −84 · ∆ ≈ 5.17 × 10 −27 for P 1 and P 2 , respectively. Numerical integration on these intervals is thus performed with absolute accuracy of 1.29 × 10 −31 and 5.17 × 10 −31 , respectively. In Fig. 5 we show basins on every fourth bisection interval, so two consecutive plots represent magnification by a factor of 2 4 = 16. We show only those intervals with length ≥ 2 −80 · ∆ ≈ 8.27 × 10 −26 . The magnification plots for P 1 demonstrate the existence of fine structure down to the smallest scale resolvable at the limit of our quadruple-precision numerics. Upon magnification of the narrow beige strip on the third-to-last interval (by a factor of 16), we find even narrower green and red strips around it. The green strip is identified by the bisection process only after five more bisection iterations beyond the last interval shown in Fig. 5, when the bisection interval is of length 2.58 × 10 −27 . Thus, the fractality length scale for P 1 (with ∆ = 0.1 and resolution 1.29 × 10 −27 ) is = 2.58 × 10 −27 . In contrast, for the cross section through P 2 , the plots indicate that the boundary becomes simple at a scale well above the numerical resolution, with the narrowest observed strip of basin found on the bisection interval of length 5.68 × 10 −15 (the green part in the middle of the 12th plot in Fig. 5). With two more iterations applied to this interval, we have an interval of length 1.42 × 10 −15 (not shown), and the midpoint of that interval belongs to the green strip. Since this is the last time this situation occurs, the fractality length scale for P 2 is = 1.42 × 10 −15 (with ∆ = 0.1 and resolution 5.17 × 10 −27 ).
Appendix E: Distribution of fractality measures Figure 8(a) shows the probability density functions for the fractality length scale estimated using the same set of line segments used to generate Figs. 4(d)-4(f). Note that is finite and broadly distributed above the numerical precision even for the true fractal in the case of S 1 , since the next smaller scale at which finer basin structure is observed can be below the level of numerical resolution. For S 1 , we have < 10 −10 for most boundary points (more than 95% of the approximately 10 6 line segments used), which is consistent with the true fractality of the boundaries. For both S 2 and S 3 , the boundaries are finite-scale fractals; however, for S 3 is larger than 10 −9 in more than 90% of the boundary points, indicating that the simple boundaries can almost surely be observed after zooming in a few times, while for S 2 is < 10 −10 in about 93% of the cases, suggesting that the simple boundaries at small scales are mostly hidden behind numerical round-off errors.
As a measure to quantify the extent to which the boundaries exhibit the Wada property, we define the construction level N level through the same bisection process we used to define . Rather than using interval length, however, N level is defined as the number of times the same situation (i.e., the end points and the midpoint all belonging to different basins) occurs in the process. Note that being able to continue the bisection process indefinitely implies that points belonging to all three basins can be found in an arbitrarily small interval, indicating the Wada property. Thus, N level can be interpreted as the depth of the Cantor-set construction levels observed by the bisection procedure, and hence a quantitative measure of . Ratio i+1/ i of consecutive basin interval lengths for the roulette system as a function of i. The inset shows a magnification of the region shaded in gray. For each shape of the roulette surface (S1, S2, or S3), ten realizations of the (random) process described in Appendix F are superimposed. The curves are to guide the eyes. the Wada property. Figure 8(b) shows the probability distributions of N level estimated from the set of line segments used for Figs. 4(d)-4(f). The construction levels for S 3 are relatively small, as expected for finite-scale fractals. However, N level for S 2 is significantly larger than for S 1 on average, which is the opposite of what one might expect, since the boundaries form a finite-scale fractal for S 2 and a true fractal for S 1 ; however, this is consistent with the observation that the complexity of the boundaries for S 2 in Fig. 4 appears to be higher than that for S 1 .