A numerical study of rigidity of hyperbolic splittings in simple two-dimensional maps

Chaotic hyperbolic dynamical systems enjoy a surprising degree of rigidity, a fact which is well known in the mathematics community but perhaps less so in theoretical physics circles. Low-dimensional hyperbolic systems are either conjugate to linear automorphisms, that is, dynamically equivalent to the Arnold cat map and its variants, or their hyperbolic structure is not smooth. We illustrate this dichotomy using a family of analytic maps, for which we show by means of numerical simulations that the corresponding hyperbolic structure is not smooth, thereby providing an example for a global mechanism which produces non-smooth phase space structures in an otherwise smooth dynamical system.

For decades, nonlinear dynamics has had a significant impact on various applications in many scientific fields, largely due to the study of seemingly simple model systems that exhibit surprisingly complex and counterintuitive behaviour, often with great relevance for real-world phenomena.Evidence of this impact can be found in the now-classic studies of one-and two-dimensional time-discrete dynamical systems, showing universal transitions from regular to chaotic dynamical behaviour, via the so-called Feigenbaum and quasiperiodic routes to chaos, which allow for fairly rigorous underpinnings through renormalisation group treatments [10,11].More importantly, the relevance of the mathematical findings have been confirmed in a variety of experiments (see, for example, [15,33]).Furthermore, simple low-dimensional chaotic maps are at the core of the foundations of statistical physics, as they provide a rigorous mechanism for the emergence of irreversibility, the approach to thermodynamic equilibrium, and the understanding of transport properties in Hamiltonian time-invariant dynamical systems, without taking recourse to phenomenological coarse graining [16,4,9].
While the case of one-dimensional maps has been studied thoroughly, both empirically and rigorously, for higher-dimensional invertible maps the situation is somewhat different.Despite a wealth of available numerical data, rigorous statements are sparse due to the considerable technical challenges involved.The majority of explicit calculations going beyond numerical simulations are based on a few paradigmatic systems, such as the Smale horseshoe or the Arnold cat map and its variants.These studies, as well as most applied investigations, build on the concept of hyperbolicity, see [22,Ch. 6 and 19] for a comprehensive overview.
In order to keep technical details to a minimum we shall focus on maps T mapping the two-dimensional torus [0, 1) Omitting technical details and looking at the orientation-preserving case only, hyperbolicity amounts to the existence of two normalised transversal vector fields e s (θ) and e u (θ), the so-called stable and unstable directions, which are invariant under the dynamics in the sense that DT (θ)e u (θ) = λ u (θ)e u (T (θ)) Here, DT (θ) denotes the Jacobian of the map, and the so-called local contraction and expansion rates obey 0 < λ s (θ) < 1 < λ u (θ).A diffeomorphism T exhibiting this property globally is referred to as an Anosov diffeomorphism.An explicit calculation of this hyperbolic structure, that is, closed expressions for e u (θ) and e s (θ), have been obtained only in very few cases where the Jacobian does not depend on the phase space coordinates.These cases serve as textbook examples for hyperbolic dynamics (see, for example, [9]).The hyperbolic structure is a key mechanism for chaotic dynamics in time-reversible systems, and hence a key ingredient for irreversibility occurring in otherwise Hamiltonian dynamics.While an explicit computation of the stable and unstable directions poses considerable challenges, a numerical calculation is more or less straightforward.This can be achieved through a suitable forward or backward iteration of the linear variational equations (1), which mimics the simple power method for computing leading eigenvalues and eigenvectors of matrices.In fact, this empirical numerical scheme is also at the heart of rigorously proving the existence of a hyperbolic splitting via an invariant cone field using the Alekseev criterion [1].We now turn to the surprising degree of rigidity enjoyed by hyperbolic systems, a fact well known in the mathematical dynamics community, but quite counterintuitive and perhaps less well known in the applied sciences.Naively one would expect that the degree of smoothness of the dynamical system is reflected by the smoothness of the hyperbolic structure; for example, one might suspect that analytic hyperbolic dynamical systems also have an analytic hyperbolic splitting.This turns out to be far from the truth.Anosov showed in [2] that the hyperbolic splitting of any smooth Anosov diffeomorphism is C α (that is, α-Hölder continuous) for some α ∈ (0, 1), but need not be C1 or even Lipschitz 1 .Here, the regularity of the hyperbolic splitting refers to that of the stable and unstable subspaces with respect to the phase space point.Higher-order differentiability is rare, making the Hölder setting natural for hyperbolic structures of (even smooth) Anosov diffeomorphisms.In the special case of symplectic smooth Anosov diffeomorphisms on the two-torus [0, 1)2 , the natural setting typically becomes C 1+α with α ∈ (0, 1), see, for example, [18].If in this case the hyperbolic splitting is of regularity 2 C 2 , the diffeomorphism is in fact necessarily smoothly conjugated to a linear Anosov diffeomorphism, that is, it is dynamically equivalent to a variant of the cat map.See [21,14,19] for these, and more nuanced related results.
This obstruction to higher regularity in arbitrarily smooth (nonlinear) hyperbolic systems is somewhat unexpected at first sight, and shows that in physical dynamical systems nonsmooth behaviour can be produced by a global mechanism.Although these results are wellestablished, to the best of our knowledge, it is difficult to find explicit examples in the literature where the lack of smoothness has been demonstrated purely through numerical methods.In this article we aim to fill this gap by explicitly demonstrating the lack of smoothness of the hyperbolic splitting for a specific class of systems that are, to a significant extent, amenable to an analytic treatment.
In order to construct the desired examples we shall resort to a class of analytic maps on the torus derived from Blaschke maps, for which full spectral information is available.This means that, in addition to ergodic invariant measures, it is also possible to compute exponential decay rates of correlation functions for analytic observations in closed form.The case of one-dimensional Blaschke maps has been studied for decades [24,7].However, the complete analytic treatment of their spectral properties is a more recent development [30,5], for a precursor see also [23].While formal generalisations to two dimensions have been proposed a while ago [28] the construction of two-dimensional hyperbolic diffeomorphisms with accessible spectra is a very recent finding [31,27,32].
More specifically, we shall investigate the properties of the hyperbolic splitting for the following two-dimensional mixing analytic map on the torus with This map can be considered as a (strong) perturbation of the Arnold cat map with (3) denoting the deformation of the map governed by the parameters µ and α.The unique fixed point of the map is given by with the Jacobian obeying det(DT (θ * )) = 1 and The above formula for the trace implies in particular that for µ > 0 the map T given by ( 2) is not smoothly conjugate to the Arnold cat map (which arises from (2) by choosing µ = 0), since the Jacobians at the unique fixed point of both maps, do not coincide.It turns out that T is a symplectic map for which Lebesgue measure is invariant and mixing.Moreover, the point spectrum of the compact Perron-Frobenius operator defined on a suitable anisotropic Hilbert space is given by {(−µ) |n| exp(inα) : n ∈ Z}, and correlation functions of analytic observables decay exponentially [31].Furthermore, the structure of (2) ensures the existence of an invariant cone field, given by the standard quadrants in the tangent space.Hence by the Alekseev criterion [1], T has a hyperbolic splitting as defined in (1).
While the hyperbolic splitting is hard, if not impossible, to compute by analytic means, it is fairly easy to get some impression of the hyperbolic structure by computing the stable and unstable manifolds of the fixed point θ * numerically, that is, by a forward and backward iteration of the map, see Figure 1 (see also [8, §4.3]).The hyperbolic splitting visible in this figure is not twice continuously differentiable, since the map is not smoothly conjugate to the Arnold cat map.However, this degree of roughness is difficult to discern visually.Results in [21] provide an upper bound on the smoothness of the hyperbolic splitting for properly nonlinear maps and suggest that typically the splitting is, up to logarithmic corrections, almost twice differentiable.We are going to demonstrate this phenomenon for our map T .In our subsequent numerical calculations we will use the numerical values µ = 0.7 and α = 0.3 for the parameters of the map.For this setting one can on the one hand still verify hyperbolicity of the map by an invariant cone condition, while on the other hand a non-vanishing value for α avoids an additional inversion symmetry, which would occur if µ exp(iα) were real valued.Nevertheless, the following numerical results do not seem to depend on these considerations3 .Manifolds of finite length have been computed numerically by forward, respectively backward, iteration (and an additional adaptive bisection scheme) of initial conditions close to the fixed point and located on the stable, respectively unstable, direction.
Since (1) relates the local expansion and contraction rates with the stable and unstable directions, the degree of smoothness of the hyperbolic splitting is mirrored by the smoothness of these rates.Hence, in order to evaluate the degree of smoothness of the hyperbolic structure it is sufficient to examine the smoothness of the local expansion rate λ u (θ).
Using (1) we numerically evaluate the local expansion rate mimicking the power method: given a phase space point θ we compute a backward orbit of finite length and then use forward iteration of (1) along this orbit to obtain a numerical approximation of λ u (θ) from the normalisation of the image vectors.At quadruple precision a backward orbit of length of about 100 is sufficient to obtain numerically accurate results for λ u (θ) up to 30 digits.Results for the local expansion rate are shown in Figure 2(a).In order to evaluate the degree of smoothness of this graph we consider difference quotients to estimate the derivative.Here we just show results for the partial difference quotients in the θ 2 -direction.Results for the symmetric first order difference quotient which estimates the first order partial derivative are shown in Figure 2(b), for a fixed value of the offset h.Again we obtain an apparently smooth function, indicating that the local expansion rate may be continuously differentiable.In order to obtain evidence that the hyperbolic splitting is not twice differentiable we shall now consider the second order difference quotient Figure 2(c) shows a rough surface with large fluctuations, suggesting that the second order difference quotient does not converge to a well-defined second derivative, as expected.In order to underpin this observation, namely that the hyperbolic splitting is not C 2 , we shall now consider the scaling behaviour of the difference quotients ( 6) and ( 7) with respect to the offset h in greater detail.(b) First order difference quotient defined in (6), for h = 10 −4 .(c) Second order difference quotient defined in (7), for h = 10 −4 .
To begin with, we look at the convergence of the first order difference quotient towards the limit value, that is, at the dependence of ∆λ u given in ( 6) on the offset h.In order to capture possible dependencies on the phase space coordinate we perform the numerical evaluation for θ-values on a regular grid of size 40 × 40.Even though one could apply more sophisticated methods this turns out to be sufficient to estimate the potential limit by computing the difference quotient for a very small value of the offset.In our case, the value h = 10 −16 turns out to be sufficient.The dependence of the difference between ∆λ u and the limit value is displayed in Figure 3(a).Clearly one observes convergence of the difference quotient towards a limit at a rate of order O(h).This rate is slower than the rate one would expect for an analytic function, but it is in line with the fact that the local expansion rate cannot be twice differentiable.The picture looks different for the second order difference quotient, see Figure 3(b).The second order difference quotient ∆ 2 λ u , given in (7), apparently fails to converge and is even unbounded as its behaviour is affected by logarithmic or sub-logarithmic corrections.Even though the precise scaling is difficult to estimate from the numerical data, the results shown in Figure 3(b) suggest that ∆ 2 λ u = O(| ln h|).This matches the best known regularity result for the hyperbolic splitting in this setting, namely that it is in the Zygmund class, see [21].Hence the data provide compelling evidence that the hyperbolic splitting has modulus of continuity O(s| log s|), implying it is C 2−ε for any ε > 0, that is, it is C 1 and its derivative is α-Hölder continuous for all α ∈ (0, 1).
In summary, we have presented a numerical illustration of a well-established fact, namely, that, phrased in colloquial terms, for analytic chaotic maps there occurs a dichotomy between maps with a non-smooth hyperbolic splitting and maps that are smoothly conjugate to the cat map or its variants.This observation points to a mechanism creating subtle structures in dynamical systems, and defies the naive belief that nature prefers smooth structures.Our computations suggest that the local expansion rate λ u (θ) is continuously differentiable but the second order derivatives fail to converge due to logarithmic corrections.Even though the precise form of the logarithmic corrections is difficult to assess from numerical simulations, the data shown in Figure 3  Cor 3.6]), which may not be as readily accessible to a general audience as they deserve to be.Furthermore, the logarithmic scaling seems to occur for typical phase space points, at least in our example.Overall, the numerics suggests that the hyperbolic splitting is C 2−ε , that is, it is almost twice differentiable -as asserted by the theory.These fine details are certainly not visible in simple illustrations such as Figure 1 which superficially indicates analytic behaviour.
We have focussed here on the simplest case of two-dimensional symplectic diffeomorphisms of the torus, which allowed for a straightforward formulation of the dichotomy between conjugacy to linear diffeomorphisms and occurrence of a non-smooth hyperbolic splitting, and enabled a simple numerical check in terms of a single scalar function like the local expansion rate.It is tempting to look at higher-dimensional cases, where, however, the dichotomy becomes much more involved.The hyperbolic splitting of a general Anosov diffeomorphism is a priori only Hölder continuous, regardless of the smoothness of the map, and a numerical check of the properties of the splitting may now require a more laborious and numerically demanding computation of stable and unstable eigenspaces (see, for instance, [12,26]).This makes the numerical detection of any kind of potential logarithmic corrections to scaling extremely challenging, as it requires scanning quantities over a large range of magnitudes.
We would also like to remark that the lack of higher-order regularity of statistical quantities pertaining to a smooth dynamical system is also of great import for applications.For example, the lack of differentiability of unstable Jacobians has posed an obstacle for the efficient computation of linear response using Ruelle's formula [29] (see [13,20,6,25] and [3, §5.3] for background and a broader overview of linear response).It is worth noting that our choice of area-preserving Blaschke maps as perturbations is not well-suited to numerically study this particular problem (due to the perturbations all sharing Lebesgue measure as their relevant invariant measure, thus giving rise to vanishing derivatives).While we expect that other area-preserving perturbations of linear toral automorphisms yield comparable numerical performance, our choice of Blaschke maps was motivated by the fact that these allow for explicit analytic computation of dynamical properties such as correlation decay and spectra of transfer operators [31,27,32].The explicit calculation of the spectra of transfer operators in these cases relies on analytic properties of the dynamical system and the introduction of suitably chosen anisotropic Hilbert spaces on which the transfer operator is compact and has a spectrum that is accessible by means of a certain composition operator.There is no obvious relation between these features and topological properties in phase space such as hyperbolic splittings.In fact, we have not been able to produce a map where a non-trivial splitting can be obtained in closed form.It is likely that an example of this type would trigger renewed interest in a striking global dynamical phenomenon that appears to have been settled from a mathematical perspective.

1 Figure 1 .
Figure 1.Stable (blue) and unstable (bronze) manifold of the fixed point (black filled circle) of the map T given in (2) and (3) for µ = 0.7 and α = 0.3.Manifolds of finite length have been computed numerically by forward, respectively backward, iteration (and an additional adaptive bisection scheme) of initial conditions close to the fixed point and located on the stable, respectively unstable, direction.

Figure 3 .
Figure 3. (a) Log-log plot of the difference between ∆λ u and a numerically "exact" estimate for the partial derivative (obtained via the difference quotient at h = 10 −16 ) in dependence on the offset h.Data are shown for different values of θ on a 40 × 40 square grid (blue), and for five selected θ-values (highlighted, bronze).(b) Semi-log plot of the dependence of the second order difference quotient ∆ 2 λ u on the offset h.The second order difference quotient is scaled by | ln h|.Data are shown for different values of θ on a 40×40 regular grid (blue), and for five selected values (highlighted, bronze).