We show that a known condition for having rough basin boundaries in bistable 2D maps holds for highdimensional bistable systems that possess a unique nonattracting chaotic set embedded in their basin boundaries. The condition for roughness is that the crossboundary Lyapunov exponent $ \lambda x $ on the nonattracting set is not the maximal one. Furthermore, we provide a formula for the generally noninteger codimension of the rough basin boundary, which can be viewed as a generalization of the Kantz–Grassberger formula. This codimension that can be at most unity can be thought of as a partial codimension, and, so, it can be matched with a Lyapunov exponent. We show in 2D noninvertible and 3D invertibleminimal models, that, formally, it cannot be matched with $ \lambda x $. Rather, the partial dimension $ D 0 ( x ) $ that $ \lambda x $ is associated with in the case of rough boundaries is trivially unity. Further results hint that the latter holds also in higher dimensions. This is a peculiar feature of rough fractals. Yet, $ D 0 ( x ) $ cannot be measured via the uncertainty exponent along a line that traverses the boundary. Consequently, one cannot determine whether the boundary is a rough or a filamentary fractal by measuring fractal dimensions. Instead, one needs to measure both the maximal and crossboundary Lyapunov exponents numerically or experimentally.
Besides chaotic attractors, nonattracting chaotic sets also have practical relevance.^{1} They are associated with, e.g., basins of attraction in multistable systems or longlived chaotic transients. Respective global properties quantified by characteristic numbers determine predictability and the lifetime of trajectories within a vicinity of the nonattracting set. Regarding nonattracting sets, the most relevant concept of predictability is that of the second kind, concerning the outcome of an experiment in terms of the final state of the system out of a few alternatives.^{2} This kind of predictability can be quantified by the uncertainty exponent,^{3,4} expressing the improvement of predictability of the outcome by determining the initial condition (IC) more precisely. To put this into context, we note that attractors, on the other hand, are associated with predictability of the first kind only, when we are concerned with how the error in the prediction of the future state of the system changes as we change the precision in the definition of the initial conditions. Error growth or decline in specific directions^{5} is measured by a spectrum of Lyapunov exponents (LEs); and, concerning an initial condition with small random error, it is the maximal positive LE (MLE) that determines the asymptotic but infinitesimal error growth for trajectories confined to the invariant set that supports a measure. Lyapunov exponents are thus regarded to quantify local instabilities, while the escape rate $\kappa $, the inverse of the expected or characteristic lifetime of trajectories mentioned above,^{1} is regarded to quantify global instability. LEs can be defined for both attractors and nonattracting invariant sets that support a measure, while attractors are clearly globally stable with no escape from them.^{1}
I. INTRODUCTION
Global characteristic numbers are, in fact, not completely independent. First, the uncertainty exponent $\alpha $ has a straightforward onetoone connection with the fractal dimension $ D b , 0 $ of the basin boundary, being simply the codimension $\alpha =D\u2212 D b , 0 $, where $D$ is the dimension of the phase space, implying that the more spacefilling the set, the poorer the predictability of the second kind. Second, a connection (i) of the predictability of the second kind and global instability, and, third, another connection (ii) of the global and local instabilities, can be given as follows. We consider here discretetime bistable systems that are possibly composed of coupling two subsystems (both of which can be multidimensional),
where the subsystems decouple for $ \u03f5 X =x \u03f5 Y =0$, such that $ f X (X,Y; \u03f5 X =0)= f ~ X (X)$ and $ f Y (X,Y; \u03f5 Y =0)= f ~ Y (Y)$. We assume that the unperturbed ( $ \u03f5 X =0$) subsystem (1) has two coexisting attractors, while the unperturbed ( $ \u03f5 Y =0$) subsystem (2) has a unique globally attracting set. That is, the bistability of the coupled system derives from that of (1). Furthermore, we assume that a single nonattracting chaotic set is embedded in the basin boundary of either (1) ( $ \u03f5 X =0$) or the coupled system (1) and (2), and, so, they possess a single unstable dimension in which escape can occur. Such a nonattracting set (either a chaotic saddle or a repeller) is said to be lowdimensional.^{1} That is, the dimensionality of the nonattracting set in this sense is not to be confused with the dimensionality of the phase space of the system. Note that the stable manifold of the unique nonattracting set coincides with the basin boundary. For lowdimensional nonattracting chaotic sets embedded in socalled filamentary fractal boundaries, “locally consisting of a Cantor set of smooth curves or surfaces,”^{1} the Kantz–Grassberger relation^{6} holds,
where $ \lambda x $ and $ D 1 ( x ) $ are the crossboundary LE (i.e., the Lyapunov exponent describing the instability of the motion across the smooth boundary filament) and associated^{1} partial information dimension, respectively. It implies that

upon a change that leaves the LEs (and so the predictability of the first kind) practically unchanged, a longer characteristic lifetime, i.e., a greater global stability, is implied by a poorer predictability of the second kind,^{2} and vice versa; and

the global instability quantified by $\kappa $ is, in general, weaker than the local instability quantified by $ \lambda x $, simply because $0\u2266 D 1 ( x ) \u22661$, as though the intricate folded geometry could trap the trajectory.
As for (i) above the change may be (i.a) with respect to a parameter $p$, while $ \u03f5 X =0$ in Eq. (1) or (i.b) by introducing a perturbation, $ \u03f5 X \u22600$. In the latter case, it is meant that the basin boundary of the coupled system is still filamentary and that the LEs of the uncoupled and coupled system are approximately the same. In the case of such a weak coupling, we use the term “perturbation” for $Y$ from the point of view of $X$. As shown by Wouters and Lucarini in Refs. 7 and 8, switching on the coupling between the subsystems by setting $ \u03f5 X \u22600$ and/or $ \u03f5 Y \u22600$ can be formally treated as a perturbation using Ruelle’s response theory.^{9}
As discussed in Sec. II B in the following, for another type of boundary, a continuous fractal boundary, which is rough, nowhere differentiable, (KG) does not apply, but rather a more generic formula (25), in which $ \lambda x $ and $(1\u2212 D 1 ( x ) )$ are replaced by the MLE $ \lambda U $ and the codimension of the boundary, i.e., $\alpha $, respectively. This implies completely new properties, including the generic invalidity of (i) and (ii), which will be further discussed in Sec. V. It becomes thus possible to have a very poor predictability of the second kind due to a strongly spacefilling boundary, while the trajectory lifetime is fairly short, with a global instability being about the same as the local crossboundary instability, that is, $\kappa \u2248 \lambda x $. The obvious—and perhaps very common—way that this situation can arise is that we have a bistable system ( $X$) that is perturbed by another system ( $Y$) whose behavior is noiselike^{10} compared to that of the bistable one. A weak enough noiselike perturbation does not alter the global instability^{11} but renders the outcome completely unpredictable (although the thickness of the boundary scales with the perturbation strength). That is, the poor predictability of the first kind of the fast system $Y$ will impact the predictability of the second kind of the coupled system. In the extreme, predictability of the second kind is completely lost, $\alpha =0$, with an extreme timescale separation $ \lambda x / \lambda U \u21920$, even if in the uncoupled/unperturbed ( $ \u03f5 X =0$) bistable system $X$ we had perfect predictability of the second kind, $ D 1 ( x ) =0$, $\alpha =1$.
We encountered such situations in the case of a bistable climate model of intermediate complexity.^{12} In this model, the socalled snowballsnow free bistability is created by the icealbedo positive feedback. This effect can be modeled by very simple 0D energy balance models (EBMs) (where “0” in “0D” refers to the spatial extension of variables). In our model studied in Ref. 12 a 1D diffusive heat equation serves the same purpose, which, however, has in common with the 0D EBM a nonchaotic solution. When this oceanice model component ( $X$) was coupled with the chaotic atmosphere ( $Y$), the originally regular basin boundary was found to turn into a practically spacefilling object.
Attaining our motivating objective, we are able to claim here that what we encountered was, in fact, a rough continuous very thick spacefilling fractal, because, first, the condition $ \lambda x < \lambda U $ found for 2D maps^{13,14} we argue in Sec. II A to be applicable to highdimensional systems, and, second, the condition was, in fact, satisfied by the coupled climate model. As for the second point, we did not measure the crossboundary LE of the coupled/perturbed climate model, only the MLE, but we did measure $ \lambda x $ of the 1D EBM in Ref. 15, and clearly it is not altered significantly by a rather weak coupling and the weakened diffusivity applied in Ref. 12, while $ \lambda x $ and $ \lambda U $ are obviously vastly different representing climatic and weather processes, respectively.
Next, in Sec. II, we present our three (ac) new technical results: (a) we reproduce the condition for roughness in a general highdimensional setting (Sec. II A), (b) provide a dimension formula for the codimension of basin boundaries that is applicable even to rough boundaries (II B 1), and (c) show that in 2 and 3D the partial dimension belonging to the crossboundary LE is 1, and suggest that it might be the case in any high dimension (Sec. II B 2). In Sec. III, we provide minimal models that feature rough boundaries, among them a prototypical model for a new kind of mixed filamentaryrough boundary. In Sec. IV, we report on our numerical computations performed to determine the fractal dimension of a rough boundary, revealing that the direct computation of the capacity dimension might not be so robust in the case of rough fractals, which could be due to their peculiar feature listed under point (c) above. Finally, in Sec. V, we discuss our results here and those in Ref. 12, including the practical (ir)relevance of the type of fractality, rough or filamentary, and pose some open questions of geophysical relevance.
II. THEORY
A. Condition for roughness
Grebogi et al.^{13} provided for the first time a condition for the roughness of the basin boundary in a 2D map when this rough boundary can be described as a Weierstarss function in terms of a Fourier series whose derivative is a nonconvergent series. The latter can be viewed as a recipe for creating a map with a rough basin boundary. However, Vollmer et al.^{14} derives the same condition in an alternative way, not requiring the boundary to be described by a Weierstrass function. This is what we reproduce next, pointing out in addition that (1) it applies not only to 2D maps but to any highdimensional one and (2) also to systems with a fractal boundary, which is filamentary when $ \lambda x = \lambda U $.
Let us consider the evolution of small perturbations $\delta x = x \u2212 x ref $ on an invariant set supporting a probability measure generated by the dynamics $ x n + 1 = f ( x n )$, which includes the coupled system (1) and (2). In the infinitesimal limit, the perturbations are governed by a system obtained by linearization $\delta x n + 1 =\u2207 f ( x n )\u22c5\delta x n $, called the variational equations. Note that the Lyapunov exponents are defined by such infinitesimal perturbations (see Sec. 4.4 of Ref. 16). Therefore, the variational equations can be seen as a natural starting point when we are deriving a condition given by the relationship of Lyapunov exponents. To start with, let us consider a 2D system, and assume that $x$ and $y$ denote some local coordinates that describe motions in a coarse sense “across” and “along” a basin boundary, respectively. Initially, let us also consider only the case of additive perturbation, i.e., $ f x (x,y; \u03f5 X =\u03f5)= f ~ x (x)+\u03f5y$ in Eq. (1). Then, we can write out the evolution equations as follows:
Instead of these recursive formulas, we want to relate the farfuture value at $i=n$ to the initial perturbation at $i=0$. This will be simplified by considering that $\delta x n $ is bounded when the initial perturbation is chosen in a special way that the perturbed trajectory stays on the boundary. That is, the degreeoffreedom of choosing such an initial perturbation is one, not two. For large $i$, $\delta y i \u226b\delta x i $, and, so, the second terms on the right hand side of (4) can be neglected. We can thus write
where the $\Lambda $’s are local Lyapunov numbers.
where we again used the fact that $\delta x n $ is bounded. With constant local Lyapunov numbers, $ \Lambda x , i = \Lambda x $ and $ \Lambda y , i = \Lambda y $, to take a simple first step, this equation would simplify to
A finite value for the left hand side would mean in the limit of $n\u2192\u221e$ that the boundary is locally smooth. It turns out that it is only possible if the Lyapunov number or exponent is smaller across the boundary than the other one ( $r<1$), yielding a convergent series. If not ( $r>1$), the boundary is not smooth locally but can be viewed as rough, not differentiable. However, in a generic setting $ \Lambda x , i $, $ \Lambda y , i $ do vary over the nonattracting set, and, so, in time along a trajectory. We observe that in that case, these are the smaller values of $j$ in Eq. (7) that need to be kept in check, corresponding to larger powers of some $r$ of the timeindependent formula (8), and, so, in the limit, these are indeed the average Lyapunov exponents, as arithmetic averages of $ln\u2061 \Lambda x , i $, $ln\u2061 \Lambda y , i $, that determine whether the series is convergent. That is,
[We can denote the (geometric) average Lyapunov numbers simply by $ \Lambda x $, $ \Lambda y $, just like the constant local Lyapunov numbers in the special case.] Besides that, we note that a spatial dependence of $\delta x 0 /\delta y 0 $ arises from the finite terms of the series. Furthermore, we now relax the assumption of an additive coupling, too, i.e., consider the generic form of (1) and (2). The only difference that this makes is that the constant $\u03f5$ is replaced by a timedependent one, $ \u03f5 i $, in Eq. (3). Following the argument that we applied to the case of timedependent $ \Lambda x , i $, $ \Lambda y , i $ above, embodied by (10)–(12), it is easy to see that a timedependent $ \u03f5 i $ will not change the condition for roughness either. Next, we point out two more possible generalizations.
(1) Consider a highdimensional boundary, with state variables $ y d $ “along” this hypersurface, $d=1,\u2026,D\u22121$, $D$ being the phase space dimension. Then, Eqs. (5) and (6) generalize straightforwardly as
where we retain the simplified notation $ \Lambda y $ for the maximal one out of all $ \Lambda y , d $. Note that in Eq. (14), we only state the asymptotic behavior, indicated by the symbol $\u223c$, which also means that we suppress the indication of a constant of proportionality. We also point out that it is the maximal $ \Lambda y $ indeed that appears in all $D\u22121$ components of the subsystem (14); the components differ only with respect to the suppressed constant of proportionality. Therefore, we can express a directional derivative by rearranging (13) for $\delta x 0 $ and dividing it by a normalized linear combination $ \u2211 d D \u2212 1 c d \delta y d , 0 $, where $ \u2211 d D \u2212 1 c d =1$. This will be finite again, clearly, if $r<1$. It should be noted that Eq. (14) apply generically, i.e., with probability one, unless a direction is taken for the directional derivative that aligns with a covariant Lyapunov vector (CLV) of the system.^{5} In that case, Eq. (14) is modified to be
where $ \Lambda y , d \u2217 $ can be any one of the positive $ \Lambda y , d $’s, belonging to the CLV in question, the same value applying to all components, $d=1,\u2026,D\u22121$, of (15). In this case, it is possible that the derivative exists given that $ \Lambda y , d \u2217 / \Lambda x <1$ while we have a rough surface because $r>1$, i.e., the surface is not necessarily rough in every direction.
(2) Roughness can be regarded as a local property because the derivative belongs to a particular locale on the boundary. (Note that the Lyapunov exponents are defined on the invariant set that supports a stationary probability measure, the saddle or repeller as far as the condition for roughness is concerned; therefore, we are characterizing the basin boundary only around the saddle/repeller, while it might have different features further away.) The neighborhood of such a locale can be characterized the same way in the case of a filamentary fractal boundary and a regular boundary. On the other hand, roughness should be considered as a “global” property too (extending to the whole of the saddle/repeller at least) because if the derivative does not exist in one locale, neither does it exist in any other one, given that the condition is based on global properties, the average Lyapunov numbers/exponents. Therefore, a filamentary fractal boundary should turn rough upon some change that changes $ \lambda x = \lambda U $ to $ \lambda x < \lambda U $. This change may be brought about either by changing $p$ of the coupled system (1) and (2) from some $ p 1 $ to $ p 2 $ or by switching on the coupling going from $ \u03f5 X =0$ to $ \u03f5 X \u22600$. However, the result is a new type of rough boundary.
B. Dimension formulas
1. Codimension of the basin boundary
In order to establish dimension formulas applying to rough boundaries, we invoke the following generic relations established in Ref. 17. The information dimensions of the unstable and stable manifolds of a nonattracting invariant set, respectively, are
where $U(S)$ is the number of positive $ \lambda i + $ (negative $ \lambda i \u2212 $) Lyapunov exponents ( $U+S=D$, $\u2212 \lambda S \u2212 \u2264\u2212 \lambda S \u2212 1 \u2212 \u2264\cdots \u2264\u2212 \lambda 1 \u2212 \u22640\u2264 \lambda 1 + \u2264\cdots \u2264 \lambda U \u2212 1 + \u2264 \lambda U + $), $I(J)$ is the largest integer for which the numerators of the fractions in Eqs. (16) and (17) are still positive, and
is the socalled metric entropy. In the latter, $ D 1 + ( j ) $ are the partial dimensions belonging to the positive LEs $ \lambda j + $, while partial dimensions $ D 1 \u2212 ( i ) $ also belong to the negative LEs $ \lambda i \u2212 $. These sum up to the full dimension of the invariant set,
Otherwise, this dimension is
because the nonattracting set is the intersection set of its stable and unstable manifolds, and, in general,^{18} the codimension of the intersection set $S= S 1 \u2229 S 2 $ is the sum of the codimensions of the two sets $ S 1 $ and $ S 2 $,
Escape takes place along directions in which the partial dimensions are less than maximal, as expressed by the following:
For invertible systems only,
Equations (16), (17), (18), (22), and (23) above are Eqs. (8.21), (8.24), (8.8), (8.7), and (8.10) of Ref. 1, respectively. Note that (KG) is only applicable to regular or filamentary fractal boundaries when $ \lambda x = \lambda j = U + = \lambda U $, and it derives rather straightforwardly from Eq. (22) as follows. Even if there was chaotic dynamics within the boundary, it would be a socalled relative attractor in that smooth subspace with no escape from it. This means that $ D 1 + ( j ) =1$, $j<U$, which leaves $ D 1 ( x ) = D 1 + ( j = U ) $ as the only nontrivial partial dimension belonging to a positive LE $ \lambda j = U + = \lambda x $, leaving a single term of the sum in Eq. (22). Nevertheless, when the perturbations are weak, an extra equation (approximation, in fact) can be obtained observing that $\kappa \u2248 \kappa ~ $ and $ \lambda x \u2248 \lambda ~ x $ (with errors quadratic and higherorder in $ \u03f5 X $), where the tilde specifies quantities belonging to the unperturbed ( $ \u03f5 X =0$) system, which are related by (KG); therefore,
Note that $ \lambda x $ is one of the $ \lambda j + $’s—not necessarily $ \lambda U $, in general.
We can now obtain a formula for the dimension $ D b , 1 $ of the basin boundary, applicable also to a rough basin boundary in the weak perturbation limit, making use of the following two facts. First, Eqs. (18) and (22) can be combined to yield $ K 1 = \u2211 j = 1 U \lambda j + \u2212\kappa $, which can be substituted in to Eq. (17) to eliminate $ K 1 $. Second, Eq. (24) and roughness, $ \lambda x < \lambda U $, implies that $\kappa < \lambda U $ and, therefore, $J=U\u22121$ in Eq. (17). Bearing in mind that the basin boundary is identical to the stable manifold of the unique nonattracting set in it—just like the stable manifold of an attractor is space filling^{1}—we have
We note that, as we set it out in Sec. I, we concern the dimension of the basin boundary with regard to the predictability of the second kind, related with the uncertainty exponent $\alpha $, and, therefore, it is the capacity dimension $ D b , 0 $ that is more important to us than the information dimension $ D b , 1 $. It just happens so that it is the latter that can be related by formulas to dynamical characteristics like the LEs. As we have already indicated, (25) is our generalization of (KG), such that (KG) is recovered when $ \lambda U = \lambda x $ in the case of a filamentary fractal basin boundary. See Eq. (5.13) of Ref. 1, identical with our Eq. (25), which was derived for a 2D noninvertible map.
2. Partial dimensions of rough basin boundaries
Next, we indicate that for rough boundaries, $ D 1 ( x ) $ becomes trivially unity, which is a rather counterintuitive new characteristic, conflicting with (KG), as it suggests—considering Eq. (22)—that escape does not take place in the crossboundary direction. Although we remark that the crossboundary direction can actually not be defined for a rough boundary. We will discuss this issue further in Sec. V, where we also propose an operative definition of the crossboundary Lyapunov exponent.
Consider a 2D noninvertible map with two positive LEs featuring a rough boundary such that $ \lambda x < \lambda y = \lambda U $. Since there are no negative LEs, the fraction in Eq. (16) disappears and $I=0$, yielding a fulldimensional unstable manifold $ D u , 1 =D=2$; and it implies also that the nonattracting set (a repeller here) is identical with the whole basin boundary, $ D s , 1 = D 1 $ (see also Sec. 8.3.1 of Ref. 1 treating the model in Ref. 13). Furthermore, Eqs. (19), (20), (22), and (25) imply that
We reiterate that $ D u , 1 =2$, $ D s , 1 = D 1 $ and $ D 1 ( x ) =1$ are completely new features of rough fractals with respect to those of filamentary fractals.
In the case of a regular unperturbed boundary $ D ~ 1 ( x ) =0$, and so
that is, the dimension can be predicted just by the Lyapunov exponents. Also, $ \lambda x \u2248\kappa $, which is indeed conflicting with (KG). In the case of a filamentary fractal unperturbed boundary (an example of which to be provided in Sec. III B), we have to know the nontrivial fractional $ D ~ 1 ( x ) $ too. We can take the point of view that $ D 1 ( y ) = D ~ 1 ( x ) + D ^ 1 $ is made up of a contribution from the filamentary fractality of the unperturbed boundary and a contribution from roughening, respectively, as if the contributions were partial dimensions. This implies that
That is, even if the perturbation is noiselike, i.e., $ \lambda x \u226a \lambda y $, the contribution from roughening can be at most $(1\u2212 D ~ 1 ( x ) )$, which achieves the maximally possible $ D 1 ( y ) =1$, that is, the contribution of roughening is not independent of $ D ~ 1 ( x ) $, and, so, the latter cannot be considered a partial dimension on its own.
Next, consider a dissipative 3D invertible map with two positive LEs and one negative LE, featuring a rough boundary such that $ \lambda x < \lambda y = \lambda U $. We have one more unknown soughtfor variable $ D 1 ( z ) $ with respect to the noninvertible 2D case, but we have one more equation as well: it is Eq. (23). We can obtain (save the algebraic manipulation), therefore, that Eqs. (26) and (27) hold, and we have a further nontrivial partial dimension,
Note that since $ \lambda y + \lambda x < \lambda z $ in the invertible dissipative system and $0\u2264 D ~ 1 ( x ) \u22641$, we have $ D 1 ( z ) <1$, indeed. We emphasize that the seemingly paradoxical $ D 1 ( x ) =1$ can hold even in invertible, possibly physical systems.
Going further, in high dimensions, the partial dimensions seem to be underdetermined by the available system of three linear equations,
where $\kappa $ is also given in terms of a LE as per Eq. (24). In $D=4$, we can easily parameterize the three remaining partial dimensions by possible values of a selected one. Let us consider the example of the 3D “foldedtowel” map of Rössler^{19} for Eq. (2), whose LEs are $ \lambda 2 + =0.430, \lambda 3 + =0.377, \lambda 1 \u2212 =3.299$, and choose $ \lambda 1 + = \lambda x =0.2$ concerning some 1D system for Eq. (1). Figure 1 shows the “possible” values of the partial dimensions. However, the only possible value for $ D 1 ( x ) $ is 1, given the constraints $0\u2264 D 1 \u2212 ( i ) \u22641$, $0\u2264 D 1 + ( j ) \u22641$. This prompts that we cannot exclude the possibility that $ D 1 ( x ) =1$ in any higher dimension.
III. MINIMAL MODELS FOR ROUGH BASIN BOUNDARIES
Next, we present some mathematical toy models of interest, which serves also as a device to make special properties of rough boundaries tangible. These models are, in fact, analytically tractable as far as global characteristic numbers are concerned, which, in fact, makes it possible for us to look for examples of rough basin boundaries. Otherwise, the construction of the basin boundary is possible only numerically. Here, we also verify numerically some analytical predictions of Lyapunov exponents and the escape rate but devote a separate section (Sec. IV) for the numerical verification of the fractal dimension given that it is not similarly straightforward.
A. Simple rough boundary
A minimal model for bistable systems with “uncomplicated” crossboundary and “complicated” inboundary dynamics can be constructed as follows:
Loosely speaking, the crossboundary dynamics is governed (dominated) by a linear equation (34), and the inboundary dynamics is a chaotic process produced by a noninvertible 1D map, the tent map (35) in this example. Since the tent map is noninvertible, the whole 2D ( $D=2$) system is noninvertible. The coupling is bidirectional. The intocrossboundary coupling is additive, of strength $d$, and the crosstoinboundary coupling is multiplicative through $\mu (x)$. The latter coupling is not needed for a truly minimal model that features a rough basin boundary. The model with $c<0$ can be considered a minimal model for a rough boundary that cannot be described by a (Weierstrass) function $W(y)$. The bistability is between two attractors at $\xb1\u221e$ with respect to $x$. For our exercise, we chose $a=1.1$, $d=0.01$, and $c=\u22121$.
Given that the perturbation is weak, the LEs of the nonattracting set are approximately those of the unstable fixed point and attractor of the uncoupled bistable and chaotic subsystems, respectively,
Note that small values of $x$ correspond to the nonattracting set embedded in the boundary, and so $\mu \u22482$ on the nonattracting set.
It can be visually checked in the case of our simple 2D model that the basin boundary is indeed not filamentary, but rather a jagged, rough curve. To this end, we initialize many trajectories randomly sprinkled in a box containing the basin boundary, and color small markers placed in these initial positions differently with respect to the different attractor reached. The result of this is shown in Fig. 2.
This ensemble can be used to check the exponential decaying of the survival rate or probability in the box over time, shown in Fig. 3, by which we can select an exponentially small fraction of the ensemble on the initial configuration that give birth to longlived trajectories. These initial conditions (ICs) are very near the basin boundary, which is the stable manifold of the nonattracting set embedded in it. Therefore, the small fraction of longlived trajectories should approach the nonattracting set as they evolve (if it is different from the boundary itself—unlike in our case) and then leave the box along the unstable manifold of the nonattracting set. That is, with this procedure, called the sprinkler method,^{1} one can construct the nonattracting set and its stable and unstable manifolds. Snapshots of the longlived trajectory ensemble are displayed in Fig. 4, revealing that the stable manifold is identical with the nonattracting set,^{20} and that the unstable manifold is space filling. These were indeed the predictions derived from generic dimension formulas in Sec. II B and in Sec. 8.3.1 of Ref. 1. These properties are, again, very different from properties of nonattracting chaotic saddle sets embedded in filamentary fractal basin boundaries of dissipative invertible systems, which saddles have a signature doublefractal geometry and, so, their stable and unstable manifolds both consist of filaments, which define two noninteger partial fractal dimensions (see Chapter 6 of Ref. 4).
We note that from the survival decay shown in Fig. 3, we can extract the escape rate $\kappa $, which supports our claim that it is just the uncoupled/unperturbed $ \lambda ~ x $ given weak perturbations and a regular unperturbed boundary, $ D ~ 1 ( x ) =0$. On the other hand, the numerically estimated Lyapunov exponents, shown in Fig. 5, also agree with our claim under (36). We estimated the LEs using the Gram–Schmidt procedure. Since the LEs of interest belong to the nonattracting set, we first have to construct that set and its natural measure. It has, in fact, been already done using the sprinkler method (Fig. 4). The estimates of the LEs at any time are taken as averages of onestep LEs over the ensemble of longlived trajectories followed.
The analytical formulas for the Lyapunov exponents given in Ref. 13 for a system similar to (34) and (35) are identical to (36). An analytical formula for the escape rate given in Ref. 21 for a generalized model of that in Ref. 13 implies for the special case of Ref. 13 that the escape rate equals the crossboundary LE; this is in agreement again with our claim based on our heuristic argument. In Ref. 22, the authors derive a formula also for the fractal dimension $ D b , 0 $, which is identical to what we derive from generic relations for the case of weak perturbations as (26) and (28).
In Sec. IV, we check if the said dimension formula applies indeed to our system. We will do this in the invertible 3D ( $D=3$) system,
where the linear equation (37) [identical with Eq. (34)] is perturbed by the invertible chaotic baker map (38) and (39). The behavior with respect to the basin boundary is similar (results not shown) as in the 2D noninvertible system (34) and (35). We will use $c=1/3$. However, note that the fractal dimension of the boundary does not depend on $c$ [see Eq. (25)]. This 3D model is a minimal model for an invertible system featuring a rough basin boundary with a dissipative dynamics on the nonattracting chaotic set. The nonattracting set in the invertible system can, therefore, be called a saddle. In contrast, the dynamics on the nonattracting chaotic set of the 2D minimal model (34) and (35), called a repeller, is not dissipative.
B. Mixed filamentary and rough boundary
To have a filamentary fractal unperturbed boundary (or perturbed but not rough), we replace the linear Eq. (34) of the crossboundary dynamics by a nonlinear one similar to that studied in Ref. 23,
where $\beta = 2 1 \u2212 q $, $\gamma =(1\u2212a)/(a+b)/2$, $a>2$, and $b>a/(a\u22122)$. Following the methodology of Sec. 2.2.3 of Ref. 1 applied to a similar but not bistable model, we can derive analytical formulas for the unperturbed LEs,
and a formula also for the unperturbed escape rate,
Of these, we need $ \kappa ~ $ and $ \lambda ~ y $ to predict the codimension of a rough boundary, on the one hand, by Eq. (27) as
On the other hand, by predicting $ \lambda x $ and $ \lambda y $, we can predict for what parameter choices do we actually get a rough boundary. Note that for $q=1$ Eq. (42) is that of the tent map (35). However, with that choice, no valid choices of $a$ and $b$ exist to yield a rough boundary, which is the very reason why we use the “camping” (or “multitent”) map (42) instead. For reliable numerical verification of the dimension formula (44), we need as large a contribution $ D ^ 1 $ from roughness, given by (29), as possible. Reasonably favorable parameter choices for this objective we have found as follows: $a=4$, $b=\pi /6+a/(a\u22122)\u22482.52$, $c=\u22121$, $d=0.2$, and $q=3$, for which we predict $ \lambda x =1.18$, $ \lambda y =2.08$, $\kappa =0.11$, $ D 1 ( y ) =0.947$, and $ D ^ 1 =0.04$. In Fig. 6, we visualize the boundary, providing successive zooms on the rightmost “filament” in order to expose its roughness. The provided value for $\kappa $ has been checked to apply—to the given approximation—to both the unperturbed and perturbed dynamics, indicating that $d$ is small enough. Yet, it is large enough to make the “filaments” not only rough but seemingly intertwined.
IV. NUMERICAL CALCULATION OF THE FRACTAL DIMENSION
The codimension $D\u2212 D b , 0 $ of the basin boundary is maximum 1, as it has to act as a separator of different regimes of phase space. Because of this, it is possible to determine this codimension as the codimension $1\u2212 D 0 ( c ) $ of the intersection set of the boundary and a straight line traversing it. This technique is easily applicable in any highdimensional system, as demonstrated in Ref. 12. It is important to appreciate that formula (21) is satisfied for almost any choice of the angle that the traversing line (set $ S 1 $) makes with the boundary (set $ S 2 $), i.e., with probability one for a random choice. Therefore, even if we found that the crossboundary partial dimension is trivially $ D 0 ( x ) =1$, applying the said simple method, we will actually measure the nontrivial codimension. Therefore, even if we wanted to utilize this distinguishing characteristic of a rough boundary for a detection purpose, we cannot.
On the technical side, it is not immediately clear how to obtain sample points of this intersection set. Instead, the codimension can supposedly be determined, being equal to the uncertainty exponent $\alpha $, based on data generated as follows. We populate a $y,z=constant$ line by equispaced ICs, and let them evolve under (37) and (39) until they approach one of the attractors, determining thereby the outcome. A visual representation of the outcomes along the line is given in Fig. 7. It is accompanied by a complicated and seemingly discontinuous function of the lifetimes of the trajectories depending on the ICs, shown in Fig. 8.
The ratio $N/ N 0 $ of uncertain boxes, in which we have different outcomes, of linear size $\epsilon $ is shown in Fig. 9. Values of $N/ N 0 $ are plotted on a logarithmic scale and those of $\epsilon $ are on a linear scale. In the diagram, two regimes can be seen, in both of which the decay is exponentially fast, only the scale (the slope in the log–lin diagram) is different. That is, the decay is faster than the usual polynomial decay that defines the uncertainty exponent:^{1} $N/ N 0 \u223c \epsilon \alpha $. Indeed, plotting in a log–log diagram, seen in Fig. 10, the data do not show a scaling of good quality in any range; the line is curved. The significance of this experience as a negative result is that this alone (without the theory presented in Sec. II A) could call into question whether the basin boundary has a fractal geometry as we know it in our case. Next, we proceed with the alternative method by which we could actually verify the fractality of the set and our prediction for the dimension.
One can apply “formally” a boxcounting algorithm to the outcome data to estimate the information dimension $ D 1 ( c ) $, as if the points where the outcome is reversed by a small shift of the IC corresponded to trajectory data points. The information $I(\epsilon )=\u2212\u2211 P i (\epsilon )ln\u2061 P i (\epsilon )$, $ P i $ being normalized box counts,^{4} as a function of the box size, or rather $\u2212ln\u2061(\epsilon )$, is shown in Fig. 11. This features a scaling of good quality, which, first, verifies fractality (beyond visuals), and, second, we estimate that $ D 1 ( c ) =0.84$. However, this estimate of the information dimension is appropriate only if the measure across the line over the supporting set is a Lebesgue measure, i.e., constant, because the outcome data do not contain information on the dynamics, only the geometry of the nonattracting set. We have already seen evidence in Fig. 4 that the measure is constant, namely, that the ensemble does not change with respect to its uniform distribution initially, achieved by random sprinkling; the distribution remains uniform. This backs the intuition that the constant measure of the uncoupled chaotic dynamics is inherited by the coupled dynamics. Therefore, $ D 1 ( c ) = D 0 ( c ) $. This agrees with what was reported in Ref. 21: “numerical computation indicates that for the repeller, the boxcounting dimension $ D 0 $ and the information dimension $ D 1 $ are equal,” as also noted in Sec. 8.3.1 of Ref. 1. This fact has a relevance to the verification of dimension formulas given in Sec. II, where $ D 1 $ appears instead of $ D 0 $, and if the measure is not uniform, with the above technique we can determine only $ D 0 $. What matters here is that we measure $ D 0 ( c ) = D 0 ( y ) =0.84$, and this numeric value agrees very well with our prediction by Eq. (28) being $ D 0 ( y ) =0.86$.
Using the same method, we have computed the partial dimension $ D 0 ( y ) $ of the new type of rough basin boundary featured by (40)–(42), too. Our result of 0.943 with a 0.007 root mean square error of fitting a straight line to data is in rather good agreement with our prediction of $ D 1 ( y ) =0.947$, even when considering that the contribution from roughness was predicted to be $ D ^ 1 =0.04$. Note that the roughness can be resolved only at sufficiently small length scales due to the weakness of the perturbation, as prompted also by Fig. 6.
V. DISCUSSION AND OUTLOOK
This work has been motivated by the finding of a fractal basin boundary in Ref. 12 in a model whose bistable ocean component is not chaotic. Here, we conjectured that the condition for rough fractals established by Grebogi et al.^{13} in a 2D noninvertible map should apply. Here, (a) we show—following the approach of Vollmer et al.^{14}—that it is really the case. Furthermore, (b) we derived a formula (25) for the codimension of this fractal basin boundary, which is a generalization of the Kantz–Grassberger (KG) formula holding only for filamentary fractal boundaries. We also found that this new formula implies unexpected properties of rough fractals, contrasting those of filamentary fractal boundaries, including the new threeway relationship of the predictability of the second kind and local and global instability. We continue discussing these now.
We speculate that the formal result (c) that the crossboundary partial dimension is maximal, i.e., $ D 1 ( x ) =1$, which is not excluded by us to hold also in higher dimensions, owes perhaps to the fact that the boundary as a Weierstrass function is nondifferentiable and continuous in the same time. However, say, in terms of the model (34) and (35), given that it is a function $W(y)$, the dimension $ D 1 ( c ) $ of the intersection set of the boundary and a straight line $y=constant$ is actually 0. We have checked this numerically with $c=0$. With twoway coupling, e.g., $c=\u22121$, we immediately have a $ D 1 ( c ) >0$.
In the 2D invertible map that we analyzed in detail, it is, instead, $ D 1 ( y ) <1$ that is not fulldimensional, suggesting that escape occurs not crossboundary, but along the $y$direction. The intuitive picture that can be attached to this is that a trajectory has to be “pushed” in the $y$direction into the “crevices” of the rough boundary “landscape” (into the gaps of the Cantor set whose fractal dimension is $ D 1 ( y ) $) to be able to escape eventually in the $x$direction. Although one can raise an issue with this interpretation, namely, that the local directionality of the rough boundary is undefined. Our consideration in this regard, however, is that partial dimensions of the basin boundary can be associated more directly to the maximal and crossboundary Lyapunov exponents, both of which can be defined without a reference to directions, in an operational manner. The latter is important also insomuch that we can check if the condition $ \lambda x < \lambda U $ for roughness is satisfied. While the wellknown Gram–Schmidt method would produce all the LEs, the problem is that we might not know which one of them is the crossboundary LE. The crossboundary LE $ \lambda x $, on the one hand, can be defined by the divergence of ensembles initialized according to the natural measure of the nonattracting set, on the two sides very near the boundary, say, in terms of the ensemble means.^{24,31} The MLE $ \lambda U $, on the other hand, can also be measured by the divergence of trajectorypairs (as done in Ref. 12). The difference between these definitions of $ \lambda x $ and $ \lambda U $ is the opposite order of (A) measuring the difference and (B) taking the ensemble averages with respect to the natural measure supported by the nonattracting set.
We are not aware of a method by which the formal results of $ D 1 ( x ) =1$ and $ D 1 ( y ) <1$ can be verified. By traversing the basin boundary with a line and evaluating the uncertainty exponent of the resulting intersection set, owing to Eq. (21), it is not necessarily the partial dimension $ D 1 ( x ) $ that we can measure, but the always nontrivial codimension $D\u2212 D b , 1 $. Therefore, from a practical point of view, as far as the predictability of the second kind alone is concerned, it makes no difference if the basin boundary is a filamentary or a rough fractal. Conversely, we cannot identify the type of fractal by measuring fractal dimensions. We conclude that it can only be identified by measuring both the maximal Lyapunov exponent, as done in Ref. 12, and the crossboundary Lyapunov exponent, upon which the theoretical condition (Sec. II A) can be checked. We have done this for our prototype model (40)–(42) featuring a new type of rough boundary, using the ensemble shown in Fig. 6 and show the results in Fig. 12. Given that the ensembles are of finite size, initially no exponential separation can be seen. The separation rate yielded by large values of $ln\u2061d$, that is, when the ensembles departed substantially from the nonattracting set, is $ln\u2061a$, and it clearly does not pertain to the nonattracting set for which Eq. (42) predicts the correct value to be $ \lambda x \u22481.18$. We drew a line of that slope into the diagram and find that the initial separation rate (at small values of $ln\u2061d$) is clearly consistent with it. Therefore, the technique proposed for identifying the type of fractal seems to be viable. While a visual identification is possible only in 2D or perhaps 3D systems, the computation of the two LEs should not be hindered by large dimensionality.
We also emphasize that our prediction of the dimension (27) assumes small perturbations, but the condition for roughness has no such assumption. Otherwise, our new Eq. (25), the generalization of the Kantz–Grassberger relation (KG), indicates that the ratio of $ \lambda x $ and $ \lambda U $ not only determines the type of the fractal, but has a bearing also on its dimension. These LEs associated respectively with weather and climatic process in our model in Ref. 12 resulted in an almost completely spacefilling basin boundary, i.e., a complete unpredictability of the outcome on fine scales.
Because of the inapplicability of (KG), point (ii) of Sec. I is not universally valid. Indeed, the local and global instabilities can coincide [compare Eqs. (25) and (28)], i.e., the fractality due to roughness does not necessarily weaken global instability (in proportion with the codimension) like in the case of filamentary fractals. However, we note that when, e.g., the roughening happens through perturbation, there could be a weakening of global instability according to the effect first reported in Ref. 11 and also examined, e.g., in Ref. 25, whereby the maximal effect occurs for some nontrivial perturbation strength $ \u03f5 X $. In contrast, (i.a) remains valid when the boundary is readily rough (when the $X$ subsystem itself maybe a coupled/perturbed system). (i.b), on the other hand, is not necessarily valid, as roughening the boundary by perturbation might change only the fractal dimension of the boundary but not the global instability, as already said. When the global instability is altered by perturbation according to Refs. 11 and 25, however, it should contribute to “further” change in the dimension similarly to filamentary fractals.^{2} The latter, of course, falls under the validity of (i.a), when the parameter in question can be taken to be the perturbation strength, $p= \u03f5 X $.
The decoupled ocean dynamics in Ref. 12 was given by a diffusive energy balance model, the Ghil–Sellers model,^{26} studied also in Ref. 15, whose solution is not chaotic, and so the fractality of the basin boundary in the coupled model came only from roughness. In a stateoftheart Earth system model, or even in an intermediate complexity model with a dynamical ocean like the Planet Simulator atmospheric model coupled to the Goldstein oceansea ice model,^{27} we expect a mixed filamentaryrough fractal basin boundary to be present. It is also clear that due to the vast atmospheric vs climatic time scale separation, the basin boundary is practically space filling. However, the answers to two related questions of geophysical relevance are not immediately clear:
In a different approach, atmospheric perturbations of the slow climatic model components can be viewed as noise perturbations, possibly inducing transitions or exits from one persistent state to another. Key questions concern the most likely path of transition as well as the expected residence time in the persistent state.^{28,29} These properties are governed by a potentiallike quantity.^{1} An efficient algorithm is proposed in Ref. 30 to estimate the potential barrier height from controlled exit time data.
ACKNOWLEDGMENTS
This work was supported financially from the EU Horizon 2020 Project CRESCENDO (Grant No. 641816) and Horizon 2020 Project TiPES (Grant No. 820970). V.L. thanks Celso Grebogi, Arkady Pikovsky, James Yorke, and the late Bruno Eckhardt for useful discussions. The authors are grateful for Tamás Tél for his feedback on an earlier version of the manuscript and for discussions during “The Mathematics of Climate and the Environment” programme of the Institut Henri Poincaré. T.B. was also supported by the Institute for Basic Science (IBS), Republic of Korea, under Grant No. IBSR028D1.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.