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

We show that a known condition for having rough basin boundaries in bistable 2D maps holds for high-dimensional bistable systems that possess a unique nonattracting chaotic set embedded in their basin boundaries. The condition for roughness is that the cross-boundary Lyapunov exponent $\lambda_x$ {\bfac on the nonattracting set} is not the maximal one. Furthermore, we provide a formula for the generally noninteger co-dimension of the rough basin boundary, which can be viewed as a generalization of the Kantz-Grassberger formula. This co-dimension that can be at most unity can be thought of as a partial co-dimension, and, so, it can be matched with a Lyapunov exponent. We show {\bfac in 2D noninvertible- and 3D invertible minimal 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. Indeed, 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 cross-boundary Lyapunov exponents numerically or experimentally.


Introduction
Beside chaotic attractors, nonattracting chaotic sets also have practical relevance [1]. They are associated with e.g. basins of attraction in multistabe systems, or, long-lived chaotic transients. Respective global properties quantified by characteristic numbers determine predictability and the life time 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 more precisely. To put this in to 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 (LE); 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 κ, the inverse of the expected or characteristic life time of trajectories mentioned above, 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] 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 = 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)- (2), and so they possesses a single unstable dimension in which escape can occur. Such a nonattracting set (either a chaotic saddle or repeller) is said to be low-dimensional [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 low-dimensional nonattracting chaotic sets embedded in so-called filamentary fractal boundaries, "locally consisting of a Cantor set of smooth curves or surfaces" [1], the Kantz-Grassberger relation [6] holds: where λ x and D (x) 1 are the cross-boundary 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 (i) upon a change that leaves the LEs (and so the predictability of the first kind) practically unchanged, a longer characteristic life time, i.e. a greater global stability, is implied by a poorer predictability of the second kind [2], and vice-versa; and (ii) the global instability quantified by κ is in general weaker than the local instability quantified by λ x , simply because 0 D (x) 1 1, as though the intricate folded geometry could trap the trajectory.
As for (i) above the change may be (i.a) wrt. 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 [7,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 [9] response theory. As discussed in Sec. 2.2 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 (24), in which λ x and (1 − D (x) 1 ) are replaced by the MLE λ max 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. 5. It becomes thus possible to have a very poor predictability of the second kind due to a strongly spacefilling boundary, while the trajectory life time is fairly short, with a global instability being about the same as the local cross-boundary instability. 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 behaviour is noise-like 1 compared to that of the bistable one. A weak enough noise-like perturbation does not alter the global instability [10] 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 /λ max → 0, even if in the uncoupled/unperturbed ( X = 0) bistable system X we had perfect predictability of the second kind, D (x) 1 = 0, α = 1. We encountered such a situations in the case of a bistable climate model of intermediate complexity [11]. 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 0-D energy balace models (EBM) (where "0" in "0-D" refers to the spatial extension of variables). In our model studied in [11] a 1-D diffusive heat equation serves the same purpose, which, however, has in common with the 0-D 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 fractal, because, first, the condition λ x < λ max found for 2D maps [12,13] we argue in Sec. 2.1 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 1-D EBM in [14], and clearly it is not altered significantly by a rather weak coupling and the weakened diffusivity applied in [11], while λ x and λ max are obviously vastly different representing climatic and weather processes, respectively.
Next, in Sec. 2 we reproduce the condition for roughness in a general setting, and also provide dimension formulae for rough boundaries. In Sec. 3 we provide minimal models that feature rough boundaries, among them a prototypical model for a new kind of mixed filamentary-rough boundary. In Sec. 4 we report on our numerical computations performed to determine the fractal dimension of the boundary. Finally, in Sec. 5 we discuss our results here and those in [11], and pose some open questions of geophysical relevance.

Condition for roughness
Grebogi et al. [12] 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. [13] 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 = λ max .
Assume that x and y denote some local coordinates that describe motions in a coarse sense "across" and "along" a basin boundary, respectively, of a discrete-time dynamical system. Therefore, the small perturbations around the boundary evolve as: where the Λ's are local Lyapunov numbers. Note that eq. (3) derives from the linear evolution equation (being a recursive formula) and from eq. (4) the x-dynamics is readily transformed out, having utilized the fact that the dynamics is constrained to a surface being the basin boundary. Note also that a time-independent constant means that, to start with, we consider the case of additive perturbation, i.e., f x (x, y; X = ) =f x (x) + y in eq. (1) simplified to a 2D situation. In eq. (3) we express δy n−1−j using eq. (4) and rearrange it as: which owes to the fact that δx n is bounded when the perturbation is chosen in a special way that the perturbed trajectory stays on the boundary. That is, the degree-of-freedom of choosing such a perturbation is one, not two. With constant local Lyapunov numbers 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. When Λ x,i , Λ y,i do vary over the nonattracting set, and so in time along a trajectory, in terms of j appearing in eq. (6), these are the smaller values of j that need to be kept in check, corresponding to larger powers of some r of the time-independent formula (7), and so in the limit, these are indeed the average Lyapunov exponents, as arithmetic averages of ln Λ x,y,i 's, that determine whether the series is convergent. That is, Λy Λx (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. We also point out that when the coupling is also nonlinear, such that the generic form of eq.
(1) applies, which results in a time-dependent i (to replace the constant in eq. (5)), the condition on roughness is unchanged. Next, we point out two more possible generalisations. 1) Consider a multi-dimensional boundary, with state variables y d "along" this hyper-surface, d = 1, . . . , D−1, D being the phase space dimension. Then, eqs. (3) and (4) generalise straightforwardly as: where we retain the simplified notation Λ y for the maximal one out of all Λ y,d . Note that in eq. (13) we only state the asymptotic behaviour, 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 (13); the components differ only wrt. the suppressed constant of proportionality. Therefore, we can express a directional derivative by rearranging (12) for δx 0 and dividing it by a normalised linear combination This will be finite again, clearly, if r < 1. It should be noted that eqs. (13) 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 eqs. (13) modify 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 of (14). 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. The neighbourhood 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 because if the derivative does not exists 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, also a filamentary fractal boundary should turn rough upon some change that changes λ x = λ max to λ x < λ max . This change may be brought about either by changing p of the coupled system (1)-(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.

Dimension formulae 2.2.1 Co-dimension of the basin boundary
In order to establish dimension formulae applying to rough boundaries, we invoke the following generic relations established in [15]. The information dimensions of the unstable and stable manifolds of a nonattracting invariant set, respectively, are: where is the largest integer for which the numerator of the fraction in eq. (15) ( (16)) is still positive, and is the so-called metric entropy. In the latter D Otherwise, this dimension is because the nonattracting set is the intersection set of its stable and unstable manifolds, and, in general [16], the co-dimension of the intersection set S 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: Eqs. (15,16,17,21,22) as the only nontrivial partial dimension belonging to a positive LE λ + j = λ x = λ U , leaving a single term of the sum in eq. (21). Nevertheless, when the perturbations are weak, an extra equation (approximation, in fact) can be obtained observing that κ ≈κ and λ x ≈λ 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. We can now obtain a formula for the dimension D b,1 , applicable also to a rough basin boundary in the weak perturbation limit, making use of the following two facts. First, eqs. (17) and (21) can be combined to yield which can be substituted in the fraction of eq. (16). Second, eq. (23) and roughness λ x < λ U implies that κ < λ U , and, therefore, J = U − 1 in eq. (16). 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: As we have already indicated, (24) is our generalisation 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. (24), which was derived for a 2D noninvertible map.

Partial dimensions of rough basin boundaries
Next, we indicate that for rough boundaries D (x) 1 becomes trivially unity, which is a rather counterintuitive new characteristic conflicting with (KG), as it suggests-considering eq. (21)-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.
Consider a 2D noninvertible map with two positive LEs featuring a rough boundary λ x < λ y = λ U . Since there are no negative LEs, the fraction in eq. (15) disappears and I = 0, yielding a full-dimensional unstable manifold D u,1 = D = 2; and also that the nonatracting set is identical with the whole basin boundary, D s,1 = D 1 . (See also Sec. 8.3.1 of [1] treating the model in [12].) Furthermore, eqs. (18,19,21,24) imply that We reiterate that D u,1 = 2, D s,1 = D 1 and D (x) 1 = 1 are completely new features of rough fractals with respect to filemanetray fractals.
In the case of a regular unperturbed boundaryD that is, the dimension can be predicted just by the Lyapunov exponents. In the case of a filamentary fractal unperturbed boundary we have to know the nontrivial fractionalD 1 +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 (x) 1 ), which achieves the maximally possible D (y) 1 = 1. I.e., the contribution of roughening is not independent ofD (x) 1 , and so the latter cannot be considered a partial dimension on its own. Next, consider a 3D invertible map with two positive LEs and one negative LE, featuring a rough boundary λ x < λ y = λ U . We have one more unknown sought-for variable D (z) 1 with respect to the noninvertible 2D case, but we have one more equation as well: it is eq. (22). We can obtain (save the algebraic manipulation), therefore, that eqs. (25)-(26) hold, and we have a further nontrivial partial dimension: Note that since λ y + λ x < λ z in the invertible dissipative system and 0 ≤D  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. (23). In D = 4, we can easily parametrize the three remaining partial dimensions by possible values of a selected one. Let us consider the example of the "folded-towel" map of Rössler [17], whose LEs are: λ + 2 = 0.430, λ + 3 = 0.377, λ − 1 = 3.299 and choose λ + 1 = λ x = 0.2. Fig. 1 shows the "possible" values of the partial dimensions. However, the only possible value for D

Minimal models for rough basin boundaries
In order to clarify some of the previous findings, we present here some mathematical toy models of interest.
Loosely speaking, the cross-boundary dynamics is governed (dominated) by a linear equation (33), and the inboundary dynamics is a chaotic process produced by a noninvertible 1D map, the tent map (34) 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/parametric 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 ±∞ wrt.
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 initialise 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 (IC) 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 nonatracting 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 long-lived trajectory ensemble are displayed in Fig. 4, revealing that the stable manifold is identical with the nonattracting set 2 , and that the unstable manifold is space filling. These were indeed the predictions derived from generic dimension formulae in Sec. 2.2 and in Sec. 8.3.1 of [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 [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 (x) 1 = 0. On the other hand, the numerically estimated Lyapunov exponents, shown in Fig. 5, also agree with our claim under (35). 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 formulae for the Lyapunov exponents given in [12] for a system similar to (33)-(34) are identical to (35). An analytical formula for the escape rate given in [18] for a generalised model of that in [12] implies for the special case of [12] that the escape rate equals the cross-boundary LE; this is in agreement again with our claim based on our heuristic argument. In [19] 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 (25) and (27).
In Sec. 4 we check if the said dimension formula applies indeed to our system. We will do this in the invertible 3D (D = 3) system: y n+1 = mod(2y n + x n , 1) y n > 1/2 mod(1 + 2(y n − 1) + x n , 1) y n ≤ 1/2 (37) z n+1 = cz n y n > 1/2 1 + c(z n − 1) y n ≤ 1/2 (38) where the linear equation (36) (identical with eq. (33)) is perturbed by the invertible chaotic baker map (37)-(38). The behaviour wrt. the basin boundary is similar (results not shown) as in the 2D noninvertible system (33)-(34). We will use c = 1/3. However, note that the fractal dimension of the boundary does not depend on c; see eq. (24). 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 (33)-(34), called a repellor, is not dissipative.

Mixed filamentary and rough boundary
To have a filamentary fractal unperturbed boundary (or perturbed but not rough), we replace the linear eq. (33) of the cross-boundary dynamics by a nonlinear one similar to that studied in [20]: y n+1 = µ/β min[ mod (y n , β), β − mod (y n , β)], µ = 1 + exp(c|x n |), 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. (26) 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. (40) is that of the tent map (34). 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 (40) instead. For reliable numerical verification of the dimension formula (43) we need as large a contribution D 1 from roughness, given by (28), as possible. Reasonably favorable parameter choices for this objective we have found as: a = 4, b = π/6 + a/(a − 2) ≈ 2.52, c = −1, d = 0.2, q = 3, for which we predict λ x = 1.18, λ y = 2.08, κ = 0.11, D (y) 1 = 0.947,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.

Numerical calculation of the fractal dimension
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 (c) 0 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 [11]. It is important to appreciate that the formula (20) 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 (x) 0 = 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 (36)-(38) 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 [1]: N/N 0 ∼ α . Indeed, plotting in a log-log diagram, seen in Fig. 10, the data does 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. 2.1) 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 1 , 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 normalised box counts [4], 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 (c) 1 = 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 does 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 0 . This agrees with the authors of [18] whose "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.2 of [1]. This fact has a relevance to the verification of dimension formulae given in Sec. 2, 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

Discussion and outlook
This work has been motivated by the finding of a fractal basin boundary in Ref. [11] in a model whose bistable ocean component is not chaotic. There we conjectured that the condition for rough fractals established by Grebogi et al. [12] in a 2D noninvertible map should apply. Here we show-following the approach of Vollmer et al. [13]-that it is really the case. Furthermore, we derived a formula (24) 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 of D (x) 1 = 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, given that it is a function W (y), the dimension D (c) 1 of the intersection set of the boundary and a straight line y, z = constant is actually 0. We have checked this numerically with c = 0 in (33)-(34). With two-way coupling, e.g. c = −1, we immediately have a D (c) In the 2D invertible map that we analyzed in detail, it is, instead, D (y) 1 < 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 (y) 1 ) to be able to escape eventually in the x-direction. Although the local directionality of the rough boundary is undefined, its partial dimensions can be associated more directly to the maximal and cross-boundary Laypunov exponents, both of which can be defined without a reference to directions, in an operational manner. The latter is important insomuch that we can check if the condition λ x < λ U for roughness is satisfied. The cross-boundary LE λ x , on the one hand, can be defined by the divergence of ensembles initialised according to the natural measure of the nonattracting set, on the two sides very near the boundary, say, in terms of the ensemble means 3 . The MLE λ U , on the other hand, can also be measured by the divergence of trajectory-pairs (as done in [11]). The difference between these definitions of λ x and λ U is the opposite order of A) measuring the difference and B) taking the ensemble averages wrt. the natural measure supported by the nonatracting set.
We are not aware of a method by which the formal results of D 1 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 [11], and the cross-boundary Lyapunov exponent, upon which the theoretical condition (Sec. 2.1) can be checked. We have done this for our prototype model (39)-(40) featuring a new type of rough boundary, using the ensemble shown in Fig. 4, 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. (41) 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 (26) assumes small perturbations; but the condition for roughness has no such assumption. Otherwise, our new eq. (24), 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 [11] 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. 1 is not universally valid. Indeed the local and global instabilities can coincide [compare eqs. (24) and (27)], 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 [10] and also examined e.g. in [22], 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 [10] and [22], 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 = X .
The decoupled ocean dynamics in [11] was given by a diffusive energy balance model, the Ghil-Sellers model [23], studied also in [14], 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 [24], 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 [10,22,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 [25,26]. These properties are governed by a potential-like quantity [1]. An efficient algorithm is proposed in Ref. [27] to estimate the potential barrier height from controlled exit time data.