Breakdown of hyperbolicity for quasiperiodic attracting invariant circles in a family of three-dimensional Henon-like maps

TL;DR

This paper investigates the breakdown of hyperbolicity in quasiperiodic attracting invariant circles within a family of three-dimensional Henon-like maps, revealing collision and loss of smoothness at parameter edges.

Contribution

It introduces a numerical approach using the Weighted Birkhoff method to analyze hyperbolicity breakdown in 3D Henon-like maps, extending understanding beyond skew product systems.

Findings

Hyperbolicity breaks down at parameter edges.

Collision and loss of smoothness of Lyapunov bundles observed.

Lyapunov exponents remain distinct during breakdown.

Abstract

We numerically study quasiperiodic normally hyperbolic attracting invariant circles that appear for certain parameter values in a family of three-dimensional Henon-like maps. These parameter values make up contour segments in the parameter space where the invariant circles have constant irrational rotation number. At the edges of these segments we find a breakdown of the hyperbolicity of the invariant circle. We observe the collision and loss of smoothness of two of the invariant Lyapunov bundles while the Lyapunov exponents all remain distinct. This is very similar to the breakdown of quasiperiodic normally hyperbolic invariant circles studied in previous works that have mostly focused on skew product type systems along with a few other special types of systems. The numerical tools we use for finding the invariant circles and calculating rotation numbers, Lyapunov exponents and bundles are based on the recently developed Weighted Birkhoff method. To apply all of these tools we need for the invariant circles to be attracting (or repelling) and for the system to be invertible. This is a severe restriction compared to alternative methods, but it is very numerically efficient and allows us to study even highly irregular circles.

Equations

$$F_{a}\left(\begin{array}[]{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}[]{c}y\\ z\\ Bx+M_{1}+M_{2}y-z^{2}\end{array}\right),\hskip 28.45274pta=(B,M_{1},M_{2})\in\mathbb{R}^{3},$$

$$\rho =\underset{n\to \infty }{\lim }\frac{{\tilde{f}}^n(x)-x}{n},$$

$$\left|\rho -\frac{m}{n}\right|>D\frac{1}{|n{|}^{\nu }}\text{for all }m,n\in \mathbb{Z},n\ne 0.$$

$$\Vert D{F}_a^n(x)v\Vert \le C_{-}({\lambda }_{-}{)}^n\Vert v\Vert \text{for all }n\ge 0.$$

$$e^{2\pi i\phi (x)}=\frac{\gamma \circ \pi (x)-P}{|\gamma \circ \pi (x)-P|}$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{k=0}^{n-1}{\Delta }_k={\int }_{S^1}\Delta dm={\int }_0^1\tilde{\Delta }(x)dx={\int }_0^1(\phi (x+\rho )-\phi (x))dx$$

$$={\int }_0^1{\int }_0^{\rho }{\phi }^{\prime }(x+s)dsdx={\int }_0^{\rho }{\int }_0^1{\phi }^{\prime }(x+s)dxds=W\rho .$$

$$B_n^{\Delta }=\frac{1}{n}\sum _{k=0}^{n-1}{\Delta }_k.$$

$$W{B}_n^{\Delta }=\frac{1}{A_n}\sum _{k=1}^{n-1}w\left(\frac{k}{n}\right){\Delta }_k,A_n=\sum _{k=1}^{n-1}w\left(\frac{k}{n}\right),$$

$$w(t)=\left\{\begin{array}[]{ll}\exp\left(\frac{1}{t(1-t)}\right),&t\in[0,1]\\ 0,&t\notin[0,1],\end{array}\right.$$

$$\left|W{B}_n^{\Delta }-{\int }_{S^1}\Delta dm\right|=\mathcal{O}\left(\frac{1}{n^k}\right)$$

$$\tilde{\Delta }(x)=\phi (x+\rho )-\phi (x)=\rho +\sum _{m\in \mathbb{Z}}{a}_m(e^{2\pi im\rho }-1)e^{2\pi imx}$$

$$E_n:=\frac{1}{A_n}\sum _{k=0}^{n-1}w\left(\frac{k}{n}\right){\Delta }_k-\rho =\frac{1}{A_n}\sum _{k=0}^{n-1}w\left(\frac{k}{n}\right)\sum _{m\in \mathbb{Z}}{a}_m(e^{2\pi im\rho }-1)e^{2\pi imk\rho }$$

$$=\frac{1}{A_n}\sum _{m\in \mathbb{Z}}{a}_m(e^{2\pi im\rho }-1)\sum _{k\in \mathbb{Z}}w\left(\frac{k}{n}\right)e^{2\pi imk\rho }$$

$$=\frac{1}{A_n}\sum _{m\in \mathbb{Z}}{a}_m(e^{2\pi im\rho }-1)\sum _{k\in \mathbb{Z}}{\int }_{\mathbb{R}}w\left(\frac{t}{n}\right)e^{2\pi imt\rho }{e}^{-2\pi ikt}dt$$

$$E_n=\frac{n}{A_n}\sum _{m\in {\mathbb{Z}}^{(\rho )}}{a}_m(e^{2\pi im\rho }-1)\sum _{k\in \mathbb{Z}}{\int }_0^{1}w(s)e^{2\pi i(m\rho -k)ns}ds$$

$$=\frac{n}{A_n}\sum _{m\in {\mathbb{Z}}^{(\rho )}}{a}_m(e^{2\pi im\rho }-1)\sum _{k\in \mathbb{Z}}\frac{1}{(2\pi i(m\rho -k)n{)}^d}{\int }_0^1{w}^{(d)}(s)e^{2\pi i(m\rho -k)ns}ds.$$

$$c_k=\frac{a_{k}f(b_k)-b_{k}f(a_k)}{f(b_k)-f(a_k)}$$

$$T_x{\mathbb{R}}^3={\mathbb{R}}^3=\underset{i=0}{\overset{s(x)}{⨁}}{H}_x^{i}x\in {\mathbb{R}}^3,\text{ such that }D{F}_a(x)H_x^i=H_{F_a(x)}^i,i=0,1,\dots ,s(x),$$

$$\underset{n\to \pm \infty }{\lim }\frac{1}{n}\log \Vert D{F}_a^n(x)v\Vert ={\lambda }_i(x),v\in H_x^i\setminus \{0\},i=0,1,\dots ,s(x),$$

$$\left\{\begin{array}[]{l}u_{k+1}=\|DF_{a}(x_{k})u_{k}\|^{-1}DF_{a}(x_{k})u_{k},\\ v_{k+1}=\frac{DF_{a}(x_{k})v_{k}-\mathrm{Proj}_{u_{k+1}}DF_{a}(x_{k})v_{k}}{\left\|DF_{a}(x_{k})v_{k}-\mathrm{Proj}_{u_{k+1}}DF_{a}(x_{k})v_{k}\right\|}\end{array}\right.$$

$$\left\{\begin{array}[]{l}\tilde{u}_{k-1}=\|DF_{a}(x_{k-1})^{-1}\tilde{u}_{k}\|^{-1}DF_{a}(x_{k-1})^{-1}\tilde{u}_{k},\\ \tilde{v}_{k-1}=\frac{DF_{a}(x_{k-1})^{-1}\tilde{v}_{k}-\mathrm{Proj}_{\tilde{u}_{k-1}}DF_{a}(x_{k-1})^{-1}\tilde{v}_{k}}{\left\|DF_{a}(x_{k-1})^{-1}\tilde{v}_{k}-\mathrm{Proj}_{\tilde{u}_{k-1}}DF_{a}(x_{k-1})^{-1}\tilde{v}_{k}\right\|}\end{array}\right.$$

$$\left\{\begin{array}[]{l}H_{x_{k}}^{0}\approx\hfil\lx@intercol\mathrm{span}(u_{k}),\\ H_{x_{k}}^{1}\approx\hfil\lx@intercol\mathrm{span}(u_{k},v_{k})\cap\mathrm{span}(\tilde{u}_{k},\tilde{v}_{k}),\\ H_{x_{k}}^{2}\approx\hfil\lx@intercol\mathrm{span}(\tilde{u}_{k}).\end{array}\right.$$

$${\lambda }_i=\underset{n\to \infty }{\lim }\frac{1}{n}\log \Vert D{F}_a^n(x_1)w_1\Vert =\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{k=1}^n\log \Vert D{F}_a(x_k)w_k\Vert$$

$${\lambda }_i=\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{k=1}^n\log \Vert D{F}_a(x_k)w(x_k)\Vert ={\int }_A\log \Vert D{F}_a(x)w(x)\Vert d\mu (x)$$

$$\left|\phi -\frac{m}{n}\right|>D\frac{1}{n^2}\text{for all }m,n\in \mathbb{Z},n\ne 0.$$

Full text

The numerical tools we use for finding the invariant circles and calculating rotation numbers, Lyapunov exponents and bundles are based on the recently developed Weighted Birkhoff method. To apply all of these tools we need for the invariant circles to be attracting (or repelling) and for the system to be invertible. This is a severe restriction compared to alternative methods, but it is very numerically efficient and allows us to study even highly irregular circles.

1 Introduction

The study of invariant manifolds has been at the very heart of dynamical systems since it’s inception, because the existence or nonexistence of them gives great insight into the long term behavior of a dynamical system. Along with the question of existence, researchers have also extensively studied the question of persistence, i.e. for what systems does the invariant manifold persist under small perturbations. Part of this research yielded the notion of a normally hyperbolic manifold, which informally is an invariant $C^{1}$ submanifold where the growth rate of the transverse vectors dominate the growth rate of the tangent vectors. The key results from [HPS70, Fen72] is then that a compact normally hyperbolic manifold without boundary that is invariant under a diffeomorphism persists under small $C^{1}$ perturbations. Later it was shown in [Mañ78] that the converse is also true.

With the persistence of compact normally hyperbolic manifolds well understood, the question then becomes to identify the scenarios under which they can break down. Naturally this is a very broad subject so we restrict ourselves to the case of normally hyperbolic invariant circles, the simplest form of compact manifold. This case has been studied extensively in the context of quasiperiodic skew product systems, see e.g. [HdlL06a]. In this paper we will examine the breakdown of normally hyperbolic invariant circles in a non-skew product type system induced by a diffeomorphism. These types of system has not been studied as extensively as the skew product case but some research has been done, see e.g. [CF12, CH17a].

When studying invariant circles we separate the cases where the internal dynamics has rational rotation number and irrational rotational number. The breakdown of normal hyperbolicity in the case of a rational rotation number is well understood. Rational rotation number and normal hyperbolicity implies the existence of an even number of periodic orbits connected by stable manifolds which together make up the invariant circle. Here the circle will break down exactly when some of the Lyapunov exponents of the transversal dynamics equals the exponent of the internal dynamics, see [Fen72]. The case of irrational rotation number is not is more subtle. Assuming that the invariant circle is $C^{1+\varepsilon}$ and then it follows by a classical result of Denjoy that the internal dynamics is topologically conjugate to a rigid rotation. Like in the periodic case it is possible for the circle to break down when some of the Lyapunov exponents of transversal dynamics equals the internal exponent which is [math]. But the quasiperiodic internal dynamics also introduce new possibilities for breakdown scenarios. In [HdlL06a, FH15] one such scenario is identified. There they present some families of skew product systems where the transversal Lyapunov exponents remain bounded away from the internal exponent but hyperbolicity is lost due to transversal directions losing regularity and colliding with the tangent direction on measure zero set. A similar scenario has been observed in the conformally symplectic setting, see [CF12]. These are very similar to the phenomenon observed in this paper.

In this paper we use the efficient computational tools developed in [DSSY17, DY18] to numerically study the breakdown of attracting invariant curves in non-skew product systems. Of course more general tools for computing normally hyperbolic invariant manifolds could be used, see e.g. [HCF*+*16], but the computational gain of the tools chosen show a great advantage in time and offer all the observables needed for characterizing the breakdown such as invariant bundles, Lyapunov exponents and regularity of invariant curve and so on.

2 Theoretical background/Setting

Throughout this paper we will be considering the family of quadratic maps $F_{a}:\mathbb{R}^{3}\to\mathbb{R}^{3}$ defined by

[TABLE]

which is one possible generalization of the classic Hénon family to three dimensions. We note that the maps all have constant Jacobian $\det(DF_{a})=B\neq 0$ and are invertible. This family has also been considered in [GOST05], where they establish more properties of the maps. Throughout this paper we will assume $B$ to be fixed and $|B|<1$, i.e. $F_{a}$ is a dissipative system.

We will be interested in identifying $F_{a}$-invariant circles, these can be identified as curves $\mathcal{C}\subseteq\mathbb{R}^{3}$ parametrized by an injection $K:S^{1}\to\mathbb{R}^{3}$ such that $K(S^{1})=\mathcal{C}$ and $F_{a}\circ K=K\circ f$. Since $F_{a}$ is invertible the induced dynamics $f:S^{1}\to S^{1}$ satisfying $F\circ K=K\circ f$ is invertible so $f$ is a circle homeomorphism whose dynamics are well understood. We will assume $f$ is orientation preserving, since if it is not we can look at $f\circ f$. Then we know that $f$ has a rotation number $\rho$ defined by

[TABLE]

where $\tilde{f}$ is a lift of $f$, i.e. any map from $\mathbb{R}$ to $\mathbb{R}$ satisfying $f\circ\pi=\pi\circ\tilde{f}$, where $\pi:\mathbb{R}\to\mathbb{R}/\mathbb{Z}=S^{1}$ is the projection $\pi(x)=x+\mathbb{Z}$.

The question now is when does these invariant circles exist? In [GOST05] they identify a curve in the $(M_{1},M_{2})$ space where the map has multipliers $e^{\pm i\varphi}$, $0<\varphi<\pi$ with $2\cos(\varphi)=-(M_{2}+1)/B$. They indicate that the conditions for both subcritical and supercritical Neimark-Sacker bifurcations are fulfilled. In particular we are interested in the supercritical bifurcation that gives birth to an attracting invariant circle. Near the Neimark-Sacker bifurcation there will be Arnold tongues with rational rotation numbers corresponding to values of $\varphi$. We are interested in invariant circles with irrational rotational number, in particular invariant circles that are persistent along some contour in the subspace. These do indeed exist between the Arnold tongues, and some justifications for this can be found in [CH17b] which establishes existence of quasiperiodic normally hyperbolic invariant tori in smooth families of analytic systems. We say more about normal hyperbolicity further down. The proof uses methods from Kolmogorov-Arnold-Moser (KAM) theory and as such demands we restrict ourselves to rotation numbers that are Diophantine.

A real number $\rho$ is Diophantine of type $(D,\nu)$ where $D>0$ and $\nu>2$ if

[TABLE]

A number is Diophantine if it is Diophantine of some type $(D,\nu)$, and we note that the set of Diophantine numbers has full Lebesgue measure in $\mathbb{R}$.

We consider only invariant circles with quasiperiodic inner dynamics that are attracting and normally hyperbolic. This means that there is a splitting $T_{\mathcal{C}}\mathbb{R}^{3}=N_{-}\oplus L$, such that $\pi_{N_{-}}:N_{-}\to\mathcal{C}$ and $\pi_{L}:L\to\mathcal{C}$ are invariant subbundles of $T_{\mathcal{C}}\mathbb{R}^{3}$, $L=DK(TS^{1})$ and there are numbers $0<\lambda_{-}<1$ and $C_{-}>0$ such that for any $x\in\mathcal{C}$ and any vector $v\in(N_{-})_{x}=\pi_{N_{-}}^{-1}(x)$

[TABLE]

Here $\|\cdot\|$ is the regular Euclidian norm on $T_{\mathcal{C}}\mathbb{R}^{3}$. For a general definition of normal hyperbolicity see [HPS70]. If $\mathcal{C}$ is a normally hyperbolic and attracting invariant circle with quasiperiodic inner dynamics we say it is completely reducible if there exists everywhere nonzero sections $e_{0}\in\Gamma(L)$, $e_{1},e_{2}\in\Gamma(N_{-})$ such that $F_{a}e_{0}=e_{0}$, $F_{a}e_{1}=\lambda_{1}e_{1}$ and $F_{a}e_{2}=\lambda_{2}e_{2}$ for real $\lambda_{1},\lambda_{2}$ such that $0<|\lambda_{2}|<|\lambda_{1}|<1$.

Lastly we note that using a method referred to as bootstrap of regularity it can be proved that if $\mathcal{C}$ is normally hyperbolic and the parametrization $K$ is $C^{0}$ then $K$ also has the same regularity as the system, in our case analytic. For more details see [HdlL06b, Section 3.2].

3 Numerical methods

There are several ways to approach the problem of studying invariant circles in a smooth dynamical system and a complete list of all methods is probably yet to be compiled. Some of the history calculating invariant tori and many of the important developments can be found in [HCF*+*16, Chapter 5.1]. One of the possible methods to be considered is the parametrization method. The parametrization method is a very robust method for computing high order approximations and can be used for more general invariant manifolds than just tori, however we chose to forego the parametrization method for a couple of reasons. One reason is that the circles we seek to study are always attracting and quasiperiodic and so it is very simple and computationally efficient to approximate the torus by iteration of a point in its basin of attraction. We will informally refer to this as the ’iteration method’ from now on. Another reason for not using the parametrization method is that it can be sensitive to irregularity of the invariant circle. This is because the method is based on Fourier expansion of the circle, and if the torus becomes highly irregular for some parameter values it can require a very large number of terms in the Fourier expansion to achieve satisfactory accuracy. The iteration method also suffers from a highly irregular circle. As we will see later the Lyapunov exponents associated to the attracting directions stay bounded away from zero but these directions can approach the invariant direction of the internal dynamics as parameters change. As explained in [HdlL06a] the attracting dynamics remains exponential but the constant terms grow when the directions approach. So while a point will converge slower toward an irregular torus it is still an exponential convergence and only adds a bit of computational time and no extra requirements on storage.

With the iteration method we can thus achieve a satisfactory approximation of an orbit on the invariant circle. It then remains to find methods for computing desired quantities for the circle. In our case this will mainly be the rotation number of the internal dynamics, and the Lyapunov exponents and their associated invariant sections on the circle. The invariant sections can be computed by iteration of tangent vectors in a similar manner of the iteration method itself, and has exponential convergence for the same reason. The rest of the quantities can be reformulated as integral problems where we integrate on the torus with respect to its invariant measure. For this we will use a recently developed method referred to as the Weighted Birkhoff method.

3.1 Calculating rotation numbers

Let us begin with an explanation of how to compute the rotation number of an invariant circle given only the points of an orbit on the circle.We base our approach on the methods described in [DDS*+*16]. First choose an appropriate plane to project the orbit to, along with a coordinate system in that plane. For simplicity of presentation we can view these coordinates as real and imaginary parts and the plane as the complex plane $\mathbb{C}$. Let $\{\gamma_{n}\}_{n}$ be the projected orbit and let $\gamma:S^{1}\to\mathbb{C}$ be a map such that $\gamma_{n}=\gamma(n\rho)$, where $\rho\in[0,1)$ is the rotation number we are looking to compute. Such a map will exist as long as the original dynamical system is of high enough differentiability.

Now pick a point $P\in\mathbb{C}$ such that $P\notin\gamma(S^{1})$, and choose a continuous function $\varphi:\mathbb{R}\to\mathbb{R}$ such that

[TABLE]

where $\pi:\mathbb{R}\to S^{1}=\mathbb{R}/\mathbb{Z}$ is the projection map. We note that $\varphi(x+1)=\varphi(x)+W$ where $W\in\mathbb{Z}$ is the winding number around $P$ and let $\tilde{\Delta}:\mathbb{R}\to\mathbb{R}$ be defined by $\tilde{\Delta}(x)=\varphi(x+\rho)-\varphi(x)$. Since $\tilde{\Delta}(x+1)=\varphi(x+\rho+1)-\varphi(x+1)=\varphi(x+\rho)+W-\varphi(x)-W=\tilde{\Delta}(x)$ we see there is $\Delta:S^{1}\to\mathbb{R}$ such that $\tilde{\Delta}=\Delta\circ\pi$. We know that the sequence $\{k\rho\}_{k=0}^{\infty}$ on $S^{1}$ is equidistributed with respect to Lebesgue measure $m$ on $S^{1}$, so with $\Delta_{k}=\Delta(k\rho)$, Birkhoff’s ergodic theorem gives that

[TABLE]

[TABLE]

In practice one usually tries to choose plane, projection and point $P$ in such a way that $W=1$. Computing the numbers $\{\Delta_{k}\}$ in practice requires a bit of work since we don’t immediately have the function $\varphi$. To do this one first finds numbers $\{\Delta_{k}^{*}\}_{k=0}^{N}\subset[0,1)$ such that $\Delta_{k}^{*}=(2\pi)^{-1}(\arg(\gamma_{k+1})-\arg(\gamma_{k}))$ for some choice of $\arg$. Then $\Delta_{k}^{*}$ and $\Delta_{k}$ will only differ by a integers. We can modify $\{\Delta_{k}^{*}\}$ by integers such that if $k\rho$ and $m\rho$ are ’close’ in $S^{1}$ then $\Delta_{k}^{*}$ and $\Delta_{m}^{*}$ should be close as well. What constitutes ’close’ will depend on how the curve behaves, and for larger $N$ we can find closer points. If this process succeeds then $\{\Delta_{k}\}$ and $\{\Delta_{k}^{*}\}$ will differ by the same integer everywhere and we can find it by noting that we want $\rho\in[0,1)$.

One problem with determining when $k\rho$ and $m\rho$ are close in $S^{1}$ is that we do not have these points but only their images under $\gamma$. Since we are looking at a projection of an invariant circle it is possible that the projected circle is self intersecting, i.e. two different points map to the same point. The solution proposed in [DDS*+*16] is to use a delay coordinate embedding $\Gamma:S^{1}\to\mathbb{R}^{2L}$ defined by $\Gamma(x):=(\gamma(x),\gamma(x+\rho),\ldots,\gamma(x+(L-1)\rho)))$, for some $L\geq 2$. This gives a new orbit $\Gamma_{n}=(\gamma_{n},\gamma_{n+1},\ldots,\gamma_{n+L-1})$ in $\mathbb{R}^{2L}$, and when we want to determine if $k\rho$ and $m\rho$ are close we compare $\Gamma_{k}$ and $\Gamma_{m}$. The theoretical justification given in [DDS*+*16] for why this works is that by the Whitney and Takens Embedding theorems the map $\Gamma$ will be an embedding for almost every curve $\gamma$, so this mapping should be a reasonably effective method to remove self intersections.

With this we have reduced the problem of computing the rotation number to the problem of approximating the limit of the Birkhoff sums

[TABLE]

The problem now is that $\{B_{n}^{\Delta}\}$ generally only have $\mathcal{O}(1/n)$-convergence, which is unsatisfactory. A very elegant and simple solution to this is to instead look at weighted Birkhoff sums

[TABLE]

where the weight $w:R\to[0,\infty)$ is a smooth bump function with support $[0,1]$ and $\int_{\mathbb{R}}w(x)dx\neq 0$. The particular weight function we use is

[TABLE]

but any choice satisfying the conditions works. If $\Delta$ is a function and $\rho$ is Diophantine, then

[TABLE]

for any $k\in\mathbb{N}$. In other words $\{W\!B_{n}^{\Delta}\}$ converges super-polynomially. Here the constants in the bounds will depend on $k$, $\rho$, $\Delta$ and $w$. The proof for this is given in [DY18]. We should note here that the constants for $\mathcal{O}\left(1/n^{k}\right)$ derived in the proof contains $\|w^{(k)}\|$ as a factor, and these grow very fast in $k$. Indeed they grow faster than $c^{k}k!$ for all $c>0$. We also note that in the case of non-Diophantine $\rho$, $\{W\!B_{n}^{\Delta}\}$ still converges but only $\mathcal{O}(1/n)$.

Another thing worth mentioning is that in the particular case of calculating rotation numbers, the Weighted Birkhoff method also converges super-polynomially for rational rotation numbers, which can be of relevance for numerical application. We are not aware of this result being published anywhere else so we give an outline of the proof. Assume for simplicity that $W=1$ and let $\varphi(x)=x+\sum_{m\in\mathbb{Z}}a_{m}e^{2\pi imx}$. Then we see that

[TABLE]

Since $\Delta_{k}=\Delta(k\rho)=\tilde{\Delta}(k\rho)$ we get that

[TABLE]

[TABLE]

[TABLE]

where in the last step we used the Poisson summation formula. If $\rho$ is rational then some of the terms are zero, so let \mathbb{Z}^{(\rho)}:=\left\{m\in\mathbb{Z}\,\big{|}\,m\rho\notin\mathbb{Z}\right\}. Then doing a change of variables and integrating by parts $d$ times we get that

[TABLE]

[TABLE]

Estimating all the terms we get that $E_{n}=\mathcal{O}(1/n^{d})$, for any $d\geq 1$. For more details we refer to [DY18] where they do the proof for more general functions and Diophantine rotation numbers. The reason it also works for our case with a rational rotation number is that our function has the special form $\varphi(x+\rho)-\varphi(x)$ which leads to cancelation of terms that otherwise would cause problems.

3.2 Tracing contours in parameter space

Once we have a reliable way of calculating rotation numbers we turn to the problem of finding the parameter values which give a certain rotation number. In other words if $R\subseteq\mathbb{R}^{2}$ is the subset of the parameter space which gives systems with attracting invariant circles and we let $r:R\to[0,1)$ be the function which gives the rotation number, then we are looking for the contours $r^{-1}(\rho)=\{x\in R\mid r(x)=\rho\}$ for some given $\rho\in[0,1)$. In our case we know that the contours we are looking for will be $1$-parametric curves. If we can locate a candidate area for the search we choose a line segment $\ell:[0,1]\to R$ that transversally crosses the contours we are interested in. If done properly the function $r\circ\ell$ will be Devil’s staircase type function and in particular it will be monotone. So by checking if $r\circ\ell(0)\leq\rho\leq r\circ\ell(1)$ or $r\circ\ell(1)\leq\rho\leq r\circ\ell(0)$ we know if the line segment crosses our desired contour.

Here there is a large variety of numerical methods for finding where a function on an interval attains a certain value, but we note that due to the Devil’s staircase nature of our function we can not utilize any method based on taking derivatives. Instead we have chosen to use the False Position Method (FPM) which gives good performance for $\rho$ with good Diophantine properties since near such a number the ’steps’ of the staircase are relatively small. The FPM is based on constructing a nested sequence of intervals $[a_{k},b_{k}]$ enclosing the solution. So for $f(t)=r\circ\ell(t)-\rho$ and $a_{0},b_{0}$ equal to $\ell(0),\ell(1)$ such that $a_{0}\leq b_{0}$. Then given $[a_{k},b_{k}]$ we compute

[TABLE]

and depending on the sign of $f(c_{k})$ we get either $[a_{k+1},b_{k+1}]=[a_{k},c_{k}]$ or $[c_{k},b_{k}]$. There are some problems when analyzing the convergence of the FPM, in particular there are examples of functions for which the method does not converge. However for our application we did not encounter any of these problems, and if one does encounter these problems there are modified versions of the method one can use.

With a good approximation $p_{0}\in R$ of one point on the contour we turn to the problem of finding more. Since the contour continues in two directions we must decide on an initial direction to search. With this we start to examine the value of $r$ on points of a circle around $p_{0}$ with some prescribed small radius. Taking two distinct points on the circle in the chosen initial direction, we evaluate $r$ there and if these values enclose $\rho$ we are done. If not we start stepping along the circle until we find two points with values of $r$ that do enclose $\rho$. Which way we start stepping depends on the values of $r$ we got and which inclination of $r$ on the circle these values suggest. When we have an enclosure we take the corresponding points and look at the line segment between them. Here we can reuse the previously explained FPM to find another good approximation $p_{1}$ of a point on the contour.

From here we repeat the same procedure for the new points using the old points to guess the initial direction to search. Since the contour does not continue indefinitely the procedure will eventually fail, in which case we shrink the radius of the circle used to search in order to get good approximations near the point of breakdown.

3.3 Lyapunov exponents and associated sections

In order to calculate the Lyapunov exponents, we see that applying Oseledet’s ergodic theorem to our (invertible) dynamical system with the derivative cocycle gives us a decomposition of the tangent spaces

[TABLE]

and

[TABLE]

where $\lambda_{0}(x)>\ldots>\lambda_{s(x)}(x)$ are the Lyapunov multipliers. We will restrict ourselves to the basin of attraction of a normally hyperbolic attractor where the subspaces $\{H_{x}^{i}\}$ constitute an orientable subbundle $H^{i}$ of the Tangent bundle. In this case the Lyapunov exponents will not depend on the point $x$ and we will write $\lambda_{0}>\ldots>\lambda_{s}$. Note that in the case of quasiperiodic attractor we get that $\lambda_{0}=0$ and we have either $s=1$ or $s=2$. To simplify exposition a bit we assume $s=2$, but the following methods works for the case $s=1$ too.

We are interested in numerically approximating $\lambda_{i}$ and $H_{x}^{i}$ for $x$ on the attractor, $i=0,1,2$. We begin by finding a procedure for approximating vectors spanning $\{H_{x}^{i}\}$. If we take a vector $u\in\mathbb{R}$ and view it as a tangent vector for some point $x$ we see that $DF_{a}^{n}(x)u$ will be very close to $H_{F_{a}^{n}(x)}^{0}$ for large $n$, as long as $u$ was not in the subspace $H_{x}^{1}\oplus H_{x}^{2}$ (and even if it was, it is possible for numerical inaccuracies to push it out of the subspace). With this in mind we look at the sequence of vectors defined by the Gramm-Schmidt scheme

[TABLE]

where $x_{k+1}=F_{a}(x_{k})$ is the orbit and $\mathrm{Proj}_{u}v=\|u\|^{-1}(u\cdot v)u$. We see then that $\mathrm{span}(u_{k})\to H^{0}_{x_{k}}$ and $\mathrm{span}(u_{k},v_{k})\to H_{x_{k}}^{0}\oplus H_{x_{k}}^{1}$ as $k\to\infty$. Since the dynamical system is invertible we can perform the same scheme backwards

[TABLE]

where $x_{k-1}=F_{a}^{-1}(x_{k})$. Then we get that $\mathrm{span}(\tilde{u}_{k})\to H_{x_{k}}^{2}$ and $\mathrm{span}(\tilde{u}_{k},\tilde{v}_{k})\to H_{x_{k}}^{1}\oplus H_{x_{k}}^{2}$ as $k\to\infty$.

In practice we choose some appropriate $N_{1},N_{2},N_{3}>0$ and $N=N_{1}+N_{2}+N_{3}$ and some point $x_{1}$ that is approximately on the attractor, as well as some choice of vectors $u_{1},v_{1},\tilde{u}_{N},\tilde{v}_{N}$. Then we compute the orbit $\{x_{k}\}_{k=1}^{N}$ and iterate the vectors $u_{1},v_{1}$ forward starting at $x_{1}$ and the vectors $\tilde{u}_{N},\tilde{v}_{N}$ backward starting at $x_{N}$. If $N_{1}$ and $N_{3}$ were large enough then the vectors $\{u_{k},v_{k},\tilde{u}_{k},\tilde{v}_{k}\}_{k=N_{1}}^{N_{1}+N_{2}}$ should give us a good approximation of $H^{0},H^{1},H^{2}$ on the orbit segment $\{x_{k}\}_{k=N_{1}}^{N_{1}+N_{2}}$ and

[TABLE]

For the Lyapunov exponents we note that if $w_{1}\in H_{x_{1}}^{i}$ and $\|w_{1}\|=1$, then we see that

[TABLE]

where $w_{k+1}=\|DF_{a}(x_{k})w_{k}\|^{-1}DF_{a}(x_{k})w_{k}$ for $k\geq 1$. Since $H_{x}^{i}$ are one-dimensional and orientable we see that there is a section $w$ on the attractor such that $w(x_{k})=w_{k}$, and $\|w(x)\|=1$ for all $x$. Then by Birkhoff’s theorem we get that

[TABLE]

where $A$ is the attractor and $\mu$ is the invariant measure induced by the dynamical system. When we have this type of integral over an invariant measure we see that we are in a setting where the Weighted Birkhoff method is applicable, and this is how we use our approximations of the sections spanning the invariant subspaces to approximate the Lyapunov exponents.

4 Results

The methods described in the previous section requires us to have a good idea of where the attractors we are looking for are, both in terms of the physical space that the dynamical system is operating on and i terms of the parameter space. Locating the attractor in physical space is not too hard and we can visually inspect it by plotting it for whichever parameter value we choose. To get a good idea of which parameters have which type of attractors we first choose some grid of values in the parameter space and for each value on the grid take a grid of points in the physical space. For each of these we iterate the point forward and test for periodic behavior. The results can be seen in Figure 1. First thing to note is that the picture is a bit different depending on if we plot the maximum period found or minimum period found. This indicates that there can be more than one attractor for each parameter value. We note that the pattern on the right that is very similar to the one found when studying the parameter space for the classic family of Arnold circle maps. Because of this we have good reason to believe there are curves of parameter values for which we have quasiperiodic attractors.

We begin with searching for attractors with rotation number equal to the golden mean $\varphi$, as this is the number with the “best” Diophantine properties and hence should give us good performance for our numerical methods. In fact the golden mean is what is called a badly approximable number, meaning there is some $D>0$ such that $\varphi$ satisfies

[TABLE]

Using the methods described in Section 3.2 we locate a parameter point giving the desired rotation number and then trace the contour as it continues up to the right. These points are plotted in Figure 2 which is a zoomed in version of Figure 1. As we are tracing the contour we are also calculating and saving the Lyapunov exponents as well as the minimum angles between the Lyapunov bundles of the attractors. These values are plotted in Figure 3. We note that one of the Lyapunov exponents stay constant around zero which is expected since quasiperiodic motion has Lyapunov exponent zero, so this is the exponent for the tangential direction. The other one or two exponents are negative which we expect from an attractor, and we note that the stay away from zero. So if there is a loss of hyperbolicity at the end of the contour it is not due to Lyapunov exponents going to zero. We can also note that there seems to be three distinct exponents at the end, indicating that attractor is completely reducible at the breakdown. Looking at the minimal angles we can see that minimum between the tangential direction and the slow contracting direction is getting very small towards the end, while the other angles all stay away from zero.

In Figure 4 we see two plots of the attractor at the last step. When looking at the whole attractor it appears to have some amount of differentiability, although if we zoom in far enough we find small ’folds’ along the curve so it is not quite as regular as it first appears. But even with these folds it does not appear the curve loses differentiability as the parameters approach the point of breakdown. Since we know the rotation number it is relatively simple to construct the conjugation. In Figure 5 we show the projection of this conjugation on one of the coordinate axes. To the left we have the projection some distance away from the breakdown and it does not give us any good reason to believe it is anything less than analytic. To the right we have the projection at the last step and it appears as if differentiability is about to become lost. In Figure 6 we show the smallest angles between the tangent direction and the slow contracting direction and here we see indeed there appears to be a collision between the two.

We note here that if it is the case that the conjugation looses differentiability at the breakdown then this places distinct limits on how regular the curve itself can be. If the curve is $k>2$ times differentiable then the conjugation is $k-1-\varepsilon$ times differentiable for any $\varepsilon>0$, c.f. [KO89, KGS87]. So if the conjugation is not differentiable then the curve can not be more than $2$ times differentiable.

Acknowledgments.

This work was supported by the Swedish Research Council Grant 621-2011-3629. I want to thank S. V. Gonchenko, Á. Haro, J. Yorke and J.-L. Figueras for fruitful discussions.