Rough basin boundaries in high dimension: Can we classify them experimentally?

Global characteristic numbers are, in fact, not completely independent. First, the uncertainty exponent $α$ has a straightforward one-to-one connection with the fractal dimension $D b , 0$ of the basin boundary, being simply the co-dimension $α=D− D b , 0$⁠, where $D$ is the dimension of the phase space, implying that the more space-filling 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 discrete-time bistable systems that are possibly composed of coupling two subsystems (both of which can be multi-dimensional),

where the subsystems decouple for $ϵ X =x ϵ Y =0$⁠, such that $f X (X,Y; ϵ X =0)= f ~ X (X)$ and $f Y (X,Y; ϵ Y =0)= f ~ Y (Y)$⁠. We assume that the unperturbed (⁠ $ϵ X =0$⁠) subsystem (1) has two co-existing attractors, while the unperturbed (⁠ $ϵ 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) (⁠ $ϵ 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 low-dimensional.¹ 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 low-dimensional nonattracting chaotic sets embedded in so-called filamentary fractal boundaries, “locally consisting of a Cantor set of smooth curves or surfaces,”¹ the Kantz–Grassberger relation⁶ holds,

where $λ x$ and $D 1 ( x )$ are the cross-boundary LE (i.e., the Lyapunov exponent describing the instability of the motion across the smooth boundary filament) and associated¹ 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,² and vice versa; and

• the global instability quantified by $κ$ is, in general, weaker than the local instability quantified by $λ x$⁠, simply because $0≦ D 1 ( x ) ≦1$⁠, 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 $ϵ X =0$ in Eq. (1) or (i.b) by introducing a perturbation, $ϵ X ≠0$⁠. 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 $ϵ X ≠0$ and/or $ϵ Y ≠0$ can be formally treated as a perturbation using Ruelle’s response theory.⁹ 

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 $λ x$ and $(1− D 1 ( x ) )$ are replaced by the MLE $λ U$ and the co-dimension of the boundary, i.e., $α$⁠, 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 space-filling boundary, while the trajectory lifetime is fairly short, with a global instability being about the same as the local cross-boundary instability, that is, $κ≈ λ 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 noise-like¹⁰ compared to that of the bistable one. A weak enough noise-like perturbation does not alter the global instability¹¹ 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, $α=0$⁠, with an extreme time-scale separation $λ x / λ U →0$⁠, even if in the uncoupled/unperturbed (⁠ $ϵ X =0$⁠) bistable system $X$ we had perfect predictability of the second kind, $D 1 ( x ) =0$⁠, $α=1$⁠.

We encountered such situations in the case of a bistable climate model of intermediate complexity.¹² In this model, the so-called snowball-snow free bistability is created by the ice-albedo 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 ocean-ice model component (⁠ $X$⁠) was coupled with the chaotic atmosphere (⁠ $Y$⁠), the originally regular basin boundary was found to turn into a practically space-filling object.

Attaining our motivating objective, we are able to claim here that what we encountered was, in fact, a rough continuous very thick space-filling fractal, because, first, the condition $λ x < λ U$ found for 2D maps^13,14 we argue in Sec. II A to be applicable to high-dimensional systems, and, second, the condition was, in fact, satisfied by the coupled climate model. As for the second point, we did not measure the cross-boundary LE of the coupled/perturbed climate model, only the MLE, but we did measure $λ 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 $λ x$ and $λ U$ are obviously vastly different representing climatic and weather processes, respectively.

Next, in Sec. II, we present our three (a-c) new technical results: (a) we reproduce the condition for roughness in a general high-dimensional setting (Sec. II A), (b) provide a dimension formula for the co-dimension 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 cross-boundary 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 filamentary-rough 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.

Grebogi et al.¹³ 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.¹⁴ 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 high-dimensional one and (2) also to systems with a fractal boundary, which is filamentary when $λ x = λ U$⁠.

Let us consider the evolution of small perturbations $δ x = x − 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 $δ x n + 1 =∇ f ( x n )⋅δ 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; ϵ X =ϵ)= f ~ x (x)+ϵy$ in Eq. (1). Then, we can write out the evolution equations as follows:

Instead of these recursive formulas, we want to relate the far-future value at $i=n$ to the initial perturbation at $i=0$⁠. This will be simplified by considering that $δ 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 degree-of-freedom of choosing such an initial perturbation is one, not two. For large $i$⁠, $δ y i ≫δ x i$⁠, and, so, the second terms on the right hand side of (4) can be neglected. We can thus write

where the $Λ$’s are local Lyapunov numbers.

In Eq. (5), we express $δ y n − 1 − j$ using Eq. (6) and rearrange it as

where we again used the fact that $δ x n$ is bounded. With constant local Lyapunov numbers, $Λ x , i = Λ x$ and $Λ y , i = Λ 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→∞$ 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 $Λ x , i$⁠, $Λ 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 time-independent formula (8), and, so, in the limit, these are indeed the average Lyapunov exponents, as arithmetic averages of $ln⁡ Λ x , i$⁠, $ln⁡ Λ y , i$⁠, that determine whether the series is convergent. That is,

[We can denote the (geometric) average Lyapunov numbers simply by $Λ x$⁠, $Λ y$⁠, just like the constant local Lyapunov numbers in the special case.] Besides that, we note that a spatial dependence of $δ x 0 /δ 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 $ϵ$ is replaced by a time-dependent one, $ϵ i$⁠, in Eq. (3). Following the argument that we applied to the case of time-dependent $Λ x , i$⁠, $Λ y , i$ above, embodied by (10)–(12), it is easy to see that a time-dependent $ϵ i$ will not change the condition for roughness either. Next, we point out two more possible generalizations.

(1) Consider a high-dimensional boundary, with state variables $y d$ “along” this hyper-surface, $d=1,…,D−1$⁠, $D$ being the phase space dimension. Then, Eqs. (5) and (6) generalize straightforwardly as

where we retain the simplified notation $Λ y$ for the maximal one out of all $Λ y , d$⁠. Note that in Eq. (14), we only state the asymptotic behavior, indicated by the symbol $∼$⁠, which also means that we suppress the indication of a constant of proportionality. We also point out that it is the maximal $Λ y$ indeed that appears in all $D−1$ 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 $δ x 0$ and dividing it by a normalized linear combination $∑ d D − 1 c d δ y d , 0$⁠, where $∑ d D − 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.⁵ In that case, Eq. (14) is modified to be

where $Λ y , d ∗$ can be any one of the positive $Λ y , d$’s, belonging to the CLV in question, the same value applying to all components, $d=1,…,D−1$⁠, of (15). In this case, it is possible that the derivative exists given that $Λ y , d ∗ / Λ 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 $λ x = λ U$ to $λ x < λ 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 $ϵ X =0$ to $ϵ X ≠0$⁠. However, the result is a new type of rough 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 $λ i +$ (negative $λ i −$⁠) Lyapunov exponents (⁠ $U+S=D$⁠, $− λ S − ≤− λ S − 1 − ≤⋯≤− λ 1 − ≤0≤ λ 1 + ≤⋯≤ λ U − 1 + ≤ λ 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 so-called metric entropy. In the latter, $D 1 + ( j )$ are the partial dimensions belonging to the positive LEs $λ j +$⁠, while partial dimensions $D 1 − ( i )$ also belong to the negative LEs $λ i −$⁠. 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,¹⁸ the co-dimension of the intersection set $S= S 1 ∩ S 2$ is the sum of the co-dimensions 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 $λ x = λ j = U + = λ U$⁠, and it derives rather straightforwardly from Eq. (22) as follows. Even if there was chaotic dynamics within the boundary, it would be a so-called 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 $λ j = U + = λ 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 $κ≈ κ ~$ and $λ x ≈ λ ~ x$ (with errors quadratic and higher-order in $ϵ X$⁠), where the tilde specifies quantities belonging to the unperturbed (⁠ $ϵ X =0$⁠) system, which are related by (KG); therefore,

Note that $λ x$ is one of the $λ j +$’s—not necessarily $λ 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 = ∑ j = 1 U λ j + −κ$⁠, which can be substituted in to Eq. (17) to eliminate $K 1$⁠. Second, Eq. (24) and roughness, $λ x < λ U$⁠, implies that $κ< λ U$ and, therefore, $J=U−1$ 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¹—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 $α$⁠, 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 $λ U = λ 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.

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 cross-boundary direction. Although we remark that the cross-boundary 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 cross-boundary Lyapunov exponent.

Consider a 2D noninvertible map with two positive LEs featuring a rough boundary such that $λ x < λ y = λ U$⁠. Since there are no negative LEs, the fraction in Eq. (16) disappears and $I=0$⁠, yielding a full-dimensional 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, $λ x ≈κ$⁠, 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 noise-like, i.e., $λ x ≪ λ y$⁠, the contribution from roughening can be at most $(1− 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 $λ x < λ y = λ U$⁠. We have one more unknown sought-for 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 $λ y + λ x < λ z$ in the invertible dissipative system and $0≤ D ~ 1 ( x ) ≤1$⁠, 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 under-determined by the available system of three linear equations,

where $κ$ 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 “folded-towel” map of Rössler¹⁹ for Eq. (2), whose LEs are $λ 2 + =0.430, λ 3 + =0.377, λ 1 − =3.299$⁠, and choose $λ 1 + = λ 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≤ D 1 − ( i ) ≤1$⁠, $0≤ D 1 + ( j ) ≤1$⁠. This prompts that we cannot exclude the possibility that $D 1 ( x ) =1$ in any higher dimension.

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 minimal model for bistable systems with “uncomplicated” cross-boundary and “complicated” in-boundary dynamics can be constructed as follows:

Loosely speaking, the cross-boundary dynamics is governed (dominated) by a linear equation (34), and the in-boundary 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 in-to-cross-boundary coupling is additive, of strength $d$⁠, and the cross-to-in-boundary coupling is multiplicative through $μ(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 $±∞$ with respect to $x$⁠. For our exercise, we chose $a=1.1$⁠, $d=0.01$⁠, and $c=−1$⁠.

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 $μ≈2$ 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 long-lived 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 long-lived 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,¹ one can construct the nonattracting set and its stable and unstable manifolds. Snapshots of the long-lived trajectory ensemble are displayed in Fig. 4, revealing that the stable manifold is identical with the nonattracting set,²⁰ 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 double-fractal 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 $κ$⁠, which supports our claim that it is just the uncoupled/unperturbed $λ ~ 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 one-step LEs over the ensemble of long-lived 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 cross-boundary 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.

To have a filamentary fractal unperturbed boundary (or perturbed but not rough), we replace the linear Eq. (34) of the cross-boundary dynamics by a nonlinear one similar to that studied in Ref. 23,

where $β= 2 1 − q$⁠, $γ=(1−a)/(a+b)/2$⁠, $a>2$⁠, and $b>a/(a−2)$⁠. 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 $κ ~$ and $λ ~ y$ to predict the co-dimension of a rough boundary, on the one hand, by Eq. (27) as

On the other hand, by predicting $λ x$ and $λ 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 “multi-tent”) 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=π/6+a/(a−2)≈2.52$⁠, $c=−1$⁠, $d=0.2$⁠, and $q=3$⁠, for which we predict $λ x =1.18$⁠, $λ y =2.08$⁠, $κ=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 $κ$ 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.

The co-dimension $D− 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 co-dimension as the co-dimension $1− D 0 ( c )$ of the intersection set of the boundary and a straight line traversing it. This technique is easily applicable in any high-dimensional 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 cross-boundary partial dimension is trivially $D 0 ( x ) =1$⁠, applying the said simple method, we will actually measure the nontrivial co-dimension. 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 co-dimension can supposedly be determined, being equal to the uncertainty exponent $α$⁠, 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 $ε$ is shown in Fig. 9. Values of $N/ N 0$ are plotted on a logarithmic scale and those of $ε$ 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:¹  $N/ N 0 ∼ ε α$⁠. 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 box-counting 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(ε)=−∑ P i (ε)ln⁡ P i (ε)$⁠, $P i$ being normalized box counts,⁴ as a function of the box size, or rather $−ln⁡(ε)$⁠, 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 box-counting 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.

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.¹³ in a 2D noninvertible map should apply. Here, (a) we show—following the approach of Vollmer et al.¹⁴—that it is really the case. Furthermore, (b) we derived a formula (25) for the co-dimension 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 three-way 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 cross-boundary 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 two-way coupling, e.g., $c=−1$⁠, 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 full-dimensional, suggesting that escape occurs not cross-boundary, 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 cross-boundary 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 $λ x < λ U$ for roughness is satisfied. While the well-known Gram–Schmidt method would produce all the LEs, the problem is that we might not know which one of them is the cross-boundary LE. The cross-boundary LE $λ 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 $λ U$⁠, on the other hand, can also be measured by the divergence of trajectory-pairs (as done in Ref. 12). The difference between these definitions of $λ x$ and $λ 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 co-dimension $D− 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 cross-boundary 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⁡d$⁠, that is, when the ensembles departed substantially from the nonattracting set, is $ln⁡a$⁠, and it clearly does not pertain to the nonattracting set for which Eq. (42) predicts the correct value to be $λ x ≈1.18$⁠. We drew a line of that slope into the diagram and find that the initial separation rate (at small values of $ln⁡d$⁠) 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 $λ x$ and $λ 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 space-filling 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 co-dimension) 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 $ϵ 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.² 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= ϵ X$⁠.

The decoupled ocean dynamics in Ref. 12 was given by a diffusive energy balance model, the Ghil–Sellers model,²⁶ 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 state-of-the-art Earth system model, or even in an intermediate complexity model with a dynamical ocean like the Planet Simulator atmospheric model coupled to the Goldstein ocean-sea ice model,²⁷ we expect a mixed filamentary-rough 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:

• With a realistic Earth-like setup, at what length scale (in phase space) can we resolve the space-filling roughness?

• Is the atmospheric perturbation strong enough to significantly alter the global instability as in Refs. 11, 25, and 2?

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 potential-like quantity.¹ An efficient algorithm is proposed in Ref. 30 to estimate the potential barrier height from controlled exit time data.

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. IBS-R028-D1.

The data that support the findings of this study are available from the corresponding author upon reasonable request.