Solving the Babylonian Problem of quasiperiodic rotation rates

TL;DR

This paper introduces the Embedding Continuation Method, a new technique based on Takens' and Birkhoff's theorems, for computing rotation rates from trajectories of quasiperiodic systems, applicable in dimensions one and two.

Contribution

The paper presents a novel method for determining rotation rates from trajectory data, filling a gap in existing techniques especially for higher dimensions.

Findings

Successfully computes rotation rates in 1D and 2D cases.

Demonstrates the method's effectiveness with illustrative examples.

Explores the dependence of rotation rates on coordinate choices.

Abstract

A trajectory $u_n := F^n(u_0), n = 0,1,2, \dots $ is quasiperiodic if the trajectory lies on and is dense in some $d$-dimensional torus, and there is a choice of coordinates on the torus $\mathbb{T}$ for which $F$ has the form $F(\theta) = \theta + \rho\bmod1$ for all $\theta\in\mathbb{T}$ and for some $\rho\in\mathbb{T}$. There is an ancient literature on computing three rotation rates $\rho$ for the Moon. %There is a literature on determining the coordinates of the vector $\rho$, called the rotation rates of $F$. (For $d>1$ we always interpret $\bmod1$ as being applied to each coordinate.) However, even in the case $d=1$ there has been no general method for computing $\rho$ given only the trajectory $u_n$, though there is a literature dealing with special cases. Here we present our Embedding Continuation Method for computing some components of $\rho$ from a trajectory. It is based on the Takens Embedding Theorem and the Birkhoff Ergodic Theorem. Rotation rates are often called "rotation numbers" and both refer to a rate of rotation of a circle. However, the coordinates of $\rho$ depend on the choice of coordinates of $\mathbb{T}$. We explore the various sets of possible rotation rates that $\rho$ can yield. We illustrate our ideas with examples in dimensions $d=1$ and $2$.

Figures

Equations

$${\theta }_{n+1}:=F({\theta }_n)={\theta }_n+\rho \mathrm{mod}1$$

$${\rho }_{\varphi }\mathrm{mod}1=a\cdot \rho \mathrm{mod}1.$$

$${\rho }_{\varphi }:=a\cdot \rho \mathrm{mod}1.$$

$$\varphi (\theta )=g(\theta )+a\cdot \theta \mathrm{mod}1.$$

$${\bar{\theta }}_{n+1}={\bar{\theta }}_n+A\rho \mathrm{mod}1.$$

$$\Delta (\theta )$$

$$=g(\theta +\rho )+a\cdot (\theta +\rho )-[g(\theta )+a\cdot (\theta )]\mathrm{mod}1\text{mbox}(fromEq.\text{ref}eqn:g)$$

$$=a\cdot \rho +g(\theta +\rho )-g(\theta )\mathrm{mod}1.$$

$$\hat{\Delta }(\theta ):=a\cdot \rho +g(\theta +\rho )-g(\theta ).$$

$${\rho }_{\varphi }:=\left(\underset{N\to \infty }{\lim }\frac{{\sum }_{n=0}^{N-1}{\hat{\Delta }}_n}{N}\right)\mathrm{mod}1,$$

$${\rho }_{\varphi }=a\cdot \rho \mathrm{mod}1.$$

$${\int }_{{\mathbb{T}}^d}\hat{\Delta }(\theta )d\theta$$

$$=a\cdot \rho +{\int }_{{\mathbb{T}}^d}g(\theta +\rho )d\theta -{\int }_{{\mathbb{T}}^d}g(\theta )d\theta$$

$$=a\cdot \rho ,$$

$$x=k+(x\mathrm{mod}1)$$

$$e^{2\pi i\varphi (\theta )}=\frac{\gamma (\theta )-P}{\Vert \gamma (\theta )-P\Vert }.$$

$$W(P):={\int }_0^1{\varphi }^{\prime }(\theta +s)ds,$$

$$B_N(f)({\theta }_0):=\frac{1}{N}\sum _{n=0}^{N-1}f({\theta }_n).$$

$$\text{mbox}{WB}_N^{[p]}(f)({\theta }_0):=\sum _{n=0}^{N-1}{\hat{w}}_{n,N}^{[p]}f({\theta }_n),\text{mbox}where{\hat{w}}_{n,N}^{[p]}=\frac{w^{[p]}(n/N)}{{\sum }_{j=0}^{N-1}{w}^{[p]}(j/N)},$$

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

$$\Theta(\theta):=\Theta_{K}^{\psi}(\theta):=\bigg{(}\psi(\theta),\psi(F(\theta)),\cdots,\psi(F^{K-1}(\theta))\bigg{)}\mbox{ for }\theta\in{\mathbb{T}^{d}}.$$

$$2d+1\le KD.$$

$$\Gamma (\theta ):={\Gamma }_K^{\psi }(\theta ):=(\Theta (\theta ),\hat{\Delta }(\theta ))\text{mbox}for\theta \in {\mathbb{T}}^d.$$

$${\Gamma }_j(\theta )=(\Theta (\theta ),\hat{\Delta }(\theta )+j),$$

$$ϵ:=\inf \{\Vert p_1-p_2\Vert :p_1,p_2\in \mathbb{U}\text{mbox}andareindifferent{\Gamma }_j\},$$

$$\Vert ({\Theta }_{n_1},{\Delta }_{n_1}+m_{n_1})-({\Theta }_{n_2},{\Delta }_{n_2}+k)\Vert <\delta ,$$

$${\rho }_{\psi }^N:=\frac{{\sum }_{n=0}^{N-1}{\hat{\Delta }}_n}{N}.$$

$$\underset{N\to \infty }{\lim }{\rho }_{\psi }^N={\rho }_{\psi }.$$

$$[[a_1,a_2]]:={\theta }_{a_2{\sigma }_2-a_1{\sigma }_1}\in {\mathbb{T}}^d.$$

$$0\le a_2{\sigma }_2-a_1{\sigma }_1<N.$$

$$\bar{\theta }\mapsto \bar{\theta }+A\rho \mathrm{mod}1$$

$$\gamma (\theta ):={\hat{\gamma }}_{-1}{z}^{-1}+\widehat{{\gamma }_0}+\widehat{{\gamma }_1}z+\widehat{{\gamma }_2}{z}^2,$$

Full text

though there is a literature dealing with special cases. Here we present our Embedding continuation method for computing $\rho_{\phi}$ from the image $\phi(\theta_{n})$ of a trajectory.

It is based on the Takens Embedding Theorem and the Birkhoff Ergodic Theorem.

Keywords: Quasiperiodic, Birkhoff Ergodic Theorem, Rotation number, Rotation rate, Takens Embedding Theorem, Circular Planar Restricted 3-Body Problem, CR3BP

1 Introduction

The goal of this paper is to show how to compute a rotation rate of a quasiperiodic discrete-time trajectory. We begin with a motivating historical example, followed by a broad overview of our approach to determining rotation rates.

Rotation rates and quasiperiodicity have been studied for millennia; namely, the Moon’s orbit has three periods whose approximate values were found 2500 years ago by the Babylonians [1]. Although computation of the periods of the Moon is an easy problem today, we use it to give context to the problems we investigate. The Babylonians found that the periods of the Moon - measured relative to the distant stars - are approximately 27.3 days (the sidereal month), 8.85 years for the rotation of the apogee (the local maximum distance from the Earth), and 18.6 years for the rotation of the intersection of the Earth-Sun plane with the Moon-Earth plane. They also measured the variation in the speed of the Moon through the field of stars, and the speed is inversely correlated with the distance of the Moon. They used their results to predict eclipses of the Moon, which occur only when the Sun, Earth and Moon are sufficiently aligned to allow the Moon to pass through the shadow of the Earth. How they obtained their estimates is not fully understood but it was through observations of the trajectory of the Moon through the distant stars in the sky. In essence they viewed the Moon projected onto the two-dimensional space of distant stars. We too work with quasiperiodic motions which have been projected into one or two dimensions.

The Moon has three periods because the Moon’s orbit is basically three-dimensionally quasiperiodic, traveling on a three-dimensional torus ${\mathbb{T}}^{3}$ that is embedded in six (position+velocity) dimensions. The torus is topologically the product of three circles, and the Moon has an (average) rotation rate – i.e. the reciprocal of the period – along each of these circles. While the Moon’s orbit has many intricacies, one can capture some of the subtleties by approximating the Sun-Earth-Moon system as three point masses using Newtonian gravitational laws. This leads to the study of the Moon’s orbit as a circular restricted three-body problem (CR3BP) in which the Earth travels on a circle about the Sun and the Moon has negligible mass. Using rotating coordinates in which the Earth and Sun are fixed while the Moon moves in three-dimensions, the orbit can thus be approximated by the above mentioned three-dimensional torus ${\mathbb{T}}^{3}$ in $\mathbb{R}^{6}$. Such a model ignores several factors including long-term tidal forces and the small influence of the other planets.

As another motivating example, the direction $\phi$ of Mars from the Earth (viewed against the backdrop of the fixed stars) does not change monotonically. This apparent non-monotonic movement is called “retrograde motion.” Now imagine that exactly once each year the direction $\phi$ is determined. How do we determine the rotation rate of Mars compared with an Earth year from such data?

This kind of problem has been unsolved in full generality even for images of one-dimensional quasiperiodic maps.

This paper considers a setting more general than just the Moon or on Mars, although both give good illustrations of our setting.

For typical discrete-time dynamical systems, it is conjectured that the three kinds of recurrent motions that are likely to be seen in a dynamical system are periodic orbits, chaotic orbits, and quasiperiodic orbits  [2]. Starting with a $d$-dimensional quasiperiodic orbit on a torus ${\mathbb{T}^{d}}$ for some $d$ and a map $\phi:{\mathbb{T}^{d}}\to S^{1}$, we establish a new method for computing rotation rates from a discrete-time quasiperiodic orbit. By discrete time, we mean that the trajectory observations are a discrete sequence $\phi_{n},n=0,1,2,\cdots$, as for example when a Poincaré return map is used for the planar circular restricted three-body problem (CR3BP), and $\phi_{n}$ is the angle of the image of the trajectory as measured from the perspective of some reference point

at the $n^{th}$ time the trajectory crosses some specified Poincaré surface.

In the rest of this introduction, we give a non-technical summary of our results, ending with a comparison to previous work on this topic. We then proceed with a more technical parts of the paper, in which we describe our methods and results in detail and give numerical examples for which we compute rotation rates.

Quasiperiodicity defined. Let ${\mathbb{T}^{d}}$ be a $d$-dimensional torus. A quasiperiodic orbit is an orbit that is dense on a $d$-dimensional torus and such that there exists a choice of coordinates $\theta\in{\mathbb{T}^{d}}:=[0,1]^{d}\bmod 1$ (where $\bmod 1$ is applied to each coordinate) for the torus such that the dynamics on the orbit are given by the map

[TABLE]

for some rotation vector $\rho\in{\mathbb{T}^{d}}$ where the coordinates $\rho_{i}$ of the $\rho$ are irrational and rationally independent, i.e. if $a_{k}$ are rational numbers for $k=1,\cdots,d$ for which $a_{1}\rho_{1}+\cdots+a_{d}\rho_{d}=0$, then $a_{k}=0$ for all $k=1,\cdots,d$. We will say such a rotation vector $\rho$ is irrational.

The Babylonian Problem. One might imagine that our goal would be to compute $\rho$ in Eq. 1 from whatever knowledge we could obtain about the torus ${\mathbb{T}^{d}}$. Although the Babylonians did not know about three-dimensional tori, they none the less obtained three meaningful rotation rates. To abstract their situation, we assume there is a smooth map $\psi:{\mathbb{T}^{d}}\to M$ where $M$ is a manifold, usually of dimension $1$ or $2$. The Babylonian Problem is to compute a rotation rate $\rho_{\psi}$ from knowledge of the projection of a trajectory. We assume we only have the values $\psi_{n}$ of $\psi$ at a sequence of times (though might have a continuous time series instead). We now describe the case where the manifold $M$ is the circle $S^{1}$.

“Projections” of a torus to a circle. Maps $\phi:{\mathbb{T}^{d}}\to S^{1}$ have a nice representation. Let $a=(a_{1},\cdots,a_{d})$ where $a_{1},\cdots,a_{d}$ are integers and let $\theta=(\theta_{1},\cdots,\theta_{d})\in{\mathbb{T}^{d}}$. The simplest $\phi$ has the form $\phi(\theta)=a_{1}\theta_{1}+\cdots+a_{d}\theta_{d}\bmod 1.$ Then $\phi$ is a continuous map of the torus to a circle. For any initial point $\theta_{0}\in{\mathbb{T}^{d}}$, we have $\theta_{n}=\theta_{0}+n(a_{1}\rho_{1}+\cdots+a_{d}\rho_{d})\bmod 1$ and in this very simple case $\theta_{n+1}-\theta_{n}=a\cdot\rho\bmod 1:=a_{1}\rho_{1}+\cdots+a_{d}\rho_{d}\bmod 1$ is constant and in this very special case we obtain a constant rotation rate for $\phi(\theta)$, namely

[TABLE]

See Eq. 9. For $d=1$, Eq. 2 says $\rho_{\phi}=a_{1}\rho$ where $a_{1}$ is an integer. The integer $a_{1}$ depends on the choice of $\phi$, so even when $|a_{1}|=1$ we can get $\rho$ for one choice and $-\rho$ for another choice.

[TABLE]

We note that for every map $\phi$ of a torus to a circle, there are integers $a_{j}$ and a periodic function $g:{\mathbb{T}^{d}}\to{\mathbb{R}}$ such that

[TABLE]

Computing a rotation rate for this map can be difficult. In fact, after we define the rotation rate below, it will turn out that Eq. 3 will still be true, independent of $g$, but this formula will not be very helpful in determining $\rho_{\phi}$ from the image of a trajectory, $\phi(\theta_{n})$.

Changes of variables. Define $\bar{\theta}:=A\theta$ where $A$ is a unimodular transformation, that is, an invertible $d\times d$ matrix with integer coefficients. In these coordinates, Eq. 1 becomes

[TABLE]

Note that $\rho$ is irrational if and only if $A\rho$ is so that concept is well defined. However, as we discuss in Section 2.2, for a given irrational $\rho$ the set of $A\rho$ for all such matrices $A$ is dense in ${\mathbb{T}^{d}}$. If for example we wanted to know the vector $\rho$ with 30-digit precision, every 30-digit vectors in ${\mathbb{T}^{d}}$ would be valid approximations for an appropriate choice of coordinate matrix $A$.

We assume throughout this paper each continuous function such as those denoted by $F,\phi,\gamma,$ and $\psi$, and each manifold is smooth, by which we mean infinitely differentiable (denoted $C^{\infty}$). This assures rapid convergence of our numerical methods.

Defining $\Delta$ and its lift $\hat{\Delta}$ for a projection $\phi:{\mathbb{T}^{d}}\to S^{1}$ to a circle. Rotation rates are key characteristics of any quasiperiodic trajectory. Suppose there exists a continuous map $\phi:{\mathbb{T}^{d}}\to S^{1}$ from the dynamical system to a circle, but we only know the image $\phi_{n}:=\phi(n\rho)$ sequence of a trajectory $F(\theta_{n})=\theta_{n+1}=\theta_{n}+\rho\bmod 1$ on a torus. Define

[TABLE]

We say $\hat{\Delta}$ is a lift of $\Delta:{\mathbb{T}^{d}}\to S^{1}$ if (i)$\hat{\Delta}:{\mathbb{T}^{d}}\to\mathbb{R}$, (ii) $\hat{\Delta}$ is continuous; and (iii) $\hat{\Delta}(\theta)\bmod 1=\Delta(\theta)$. Motivated by Eq. 6, we define

[TABLE]

Then (i),(ii), and (iii) are satisfied so $\hat{\Delta}$ is a lift of $\Delta$.

Define $\hat{\Delta}_{n}=\hat{\Delta}(\theta_{n})$.

Proposition 1.1

Assume $\theta_{n}$ is quasiperiodic. There is a well-defined rotation rate $\rho_{\phi}$,

[TABLE]

and using the notation of Eq. 4,

[TABLE]

The existence of the limit is guaranteed by the Birkhoff Ergodic Theorem (See Theorem 1.2), which says that the limit in Eq. 8 is

[TABLE]

since $\int_{\mathbb{T}^{d}}d\theta=1$ and the two integrals in Eq. 10 are equal. Hence $\rho_{\phi}=a\cdot\rho\bmod 1.$

Different choices of the lift $\hat{\Delta}$ can change $\rho_{\phi}$ by an integer, so $\rho_{\phi}\bmod 1$ is independent of the choice of lift $\hat{\Delta}$. The rotation rate is this $\rho_{\phi}\bmod 1$.

The limit in Eq. 8 exists and is the same for all initial $\theta_{0}$.

A caveat. We note however, that in practice we do not know $a\cdot\rho$ so in practice we need to determine numerically what $\hat{\Delta}$ is and we must numerically evaluate the limit.

Throughout this paper we consider ${\mathbb{T}^{d}}$ to be $[0,1]^{d}\bmod 1$, where each copy of $[0,1]$ is the fraction of revolution around a circle. Furthermore $\theta\in{\mathbb{T}^{d}}$ can be treated as a set of $d$ real numbers in $[0,1)$. This will enable us to write $x\in{\mathbb{R}}$ unambiguously as

[TABLE]

where $k$ is an integer.

There are cases where it is easy to compute the rotation rate $\rho_{\phi}$. If the angle always makes small positive increases, we can convert $\phi_{n+1}-\phi_{n}\bmod 1$ into a small real positive number in $[0,1)$, and we can think of $\Delta_{n}=\phi_{n+1}-\phi_{n}$ as numbers in $(0,\alpha)$, where $0<\alpha<1$. The limit of the average of $\Delta_{n}$ is the rotation rate. The average of two or more angles in $S^{1}$ is not well defined. Hence we must average real numbers, not angles, and making that transition can be difficult.

Numerical determination of a lift $\hat{\Delta}$. The essential problem in computing $\rho_{\phi}$ is the determination of a lift $\hat{\Delta}$ for $\phi$. Given a lift, we can compute $\rho_{\phi}$ using Eq. 8. While we know the fractional part of $\hat{\Delta}$ is $\Delta\in[0,1),$ as we will explain later, we must choose the integer part $k_{n}$ of each $\hat{\Delta}_{n}$ so that all of the points $(\theta_{n},\hat{\Delta}_{n}):=(\theta_{n},k_{n}+\Delta_{n})$ lie on a connected curve in $S^{1}\times{\mathbb{R}}$ (for $d=1)$ or a connected surface in in ${\mathbb{T}^{d}}\times{\mathbb{R}}$ (for $d>1)$. We must choose these integer parts despite the fact that we do not know which $\theta_{n}$ corresponds to $\Delta_{n}$.

Even in that case $d=1$ there has been no general method for computing the lift in order to find $\rho_{\phi}$, though there is a literature dealing with special cases. See for example [3, 4, 5]. We have established a general method for determining the lift $\hat{\Delta}$, as summarized in the Figs. 1-5. Our method is based on the Theorem 1.3, a version of the Embedding Theorems of Whitney and Takens, described in detail in Section 2.

Defining $\phi$ from a planar projection $\gamma$. Assume that we are given a planar projection $\gamma:{\mathbb{T}^{d}}\to{\mathbb{R}}^{2}$ and the images $\gamma(\theta_{n})$. Fix a reference point $P\in{\mathbb{R}}^{2}$ that is not in the image $\gamma({\mathbb{T}^{d}})$. Let ${\mathbb{R}}^{2}$ be the complex plane $\mathbb{C}$, so that we can define $\phi(\theta)\in[0,1)\bmod 1=S^{1}$ by

[TABLE]

The winding number around $P$ is

[TABLE]

where $\phi^{\prime}=\frac{d\phi}{dt}$. Note that $W(P)$ is an integral over the circle so it does not depend on $\theta$. The value of $W$ is piecewise constant and integer-valued. In our examples, it is critical that the projection of our quasiperiodic trajectory into ${\mathbb{R}}^{2}$ is such that there exists a point $P$ in ${\mathbb{R}}^{2}$ with $|W(P)|=1$. That is because the measured rotation rate will be higher by a factor of $|W(P)|$. For degenerate cases, there may be no point for which $|W(P)|=1$, as shown in the next paragraph.

A non-generic map $\gamma$. Consider the map given by $\gamma(z)=z^{2}$ where $z\in\mathbb{C}$. The map $\gamma$ maps the unit circle onto the unit circle and for any value of $P\in\mathbb{C}$, $W(P)=0$ if the reference point $P$ is outside that circle, and $W(P)=2$ if inside, and $W(P)$ is not defined if $P$ is on the unit circle. Thus there is no point $P$ such that $W(P)=1$.

Two illustrative examples of complicated images of a quasiperiodic process. Figure 1 shows the projections maps $\gamma:S^{1}\to\mathbb{R}^{2}$, showing how the winding number differs in different connected components of the figure. On the left panel, every point inside the interior connected region that contains $P_{1}$ can act as a reference point for measuring angles and yields the same value of $\rho_{\phi}$. If the map is sufficiently simple, (i.e., the nonlinearity $g$ in Eq. 4 is sufficiently small), the rotation rate can immediately be computed as the average of these angle differences. However, if the map $\gamma$ is more complicated, measurement of angle is compounded by overlap of lifts of the angle between two iterates, since they can be represented by multiple values (values differing by an integer).

Projections to ${\mathbb{R}}$. Sometimes we are only provided with a scalar-valued function $\gamma:{\mathbb{T}^{d}}\to{\mathbb{R}}$, and yet we can still construct a two-dimensional map and use the methods described for ${\mathbb{R}}^{2}$ projections. For example, Fig. 2 shows how we can recover a planar map from only the first component $Re~{}\gamma_{n}$ of the flower map by considering planar points $(Re~{}\gamma_{n-1},Re~{}\gamma_{n})$. This map still gives same rotation rate as obtained by using the map in Fig. 1.

A similar example occurs with the Moon. The mean time between lunar apogees is 27.53 days, slightly longer than the 27.3-day sidereal month. Suppose we measure the distance $D_{n}$ between the centers of the Earth and Moon once each sidereal month, $n=0,1,2,\cdots$. Then the sequence $D_{n}$ has an oscillation period of 8.85 years and can be measured using our approach by plotting $D_{n-1}$ against $D_{n}$, and the point $(D_{n-1},D_{n})$ oscillates around a point $P=(D_{av},D_{av}),$ where $D_{av}$ is the average of the values $D_{n}$. Small changes in $P$ have no effect on the rotation rate.

Yet another case arises from The Moon’s orbit being tilted about 5 degrees from the Earth-Sun plane. The line of intersection where the Moon’s orbit crosses the Earth-Sun plane precesses with a period of 18.6 years. The plane of the ecliptic is a path in the distant stars through which the planets travel. Measuring the Moon’s angular distance from this plane once each sidereal month gives scalar time series with that period of 18.6 years. This example can be handled like the apogee example above.

As a last example, see also our treatment of the circular planar restricted three body problem in Section 4.2 where we compute two rotation rates of the lunar orbit, the first by plotting the rotation rate around a central point and the second by plotting $(r,dr/dt)$, deriving the rotation rate from a single variable $r(t)$, the distance from a central point, where $t$ is time.

The Birkhoff Ergodic Theorem. This theorem assumes there is an invariant set, which in our case is the set ${\mathbb{T}^{d}}$. Since we are interested here only in quasiperiodic dynamics, we can assume the dynamics are given by Eq. 1 where $\rho$ is irrational. Lebesgue measure is invariant; that is, each measurable set $E$ has the same measure as $F(E)=E+\rho$ and as $F^{-1}(E)=E-\rho$. This map is “ergodic” because if $E$ is a set for which $E=F(E)=E+\rho$, then the measure of $E$ is either [math] or $1$.

The measure $\mu$ enables the computation of the space-average $\int_{\mathbb{T}^{d}}fd\mu$ for any $L^{1}$ function $f:{\mathbb{T}^{d}}\to{\mathbb{R}}$ when a time series is the only information available. Since $\mu$ is Lebesgue measure, we can rewrite that integral as $\int_{\mathbb{T}^{d}}f(\theta)d\theta$. We note that the Lebesgue measure of the entire torus is 1, so Lebesgue measure is a probability measure. Hence $\int_{\mathbb{T}^{d}}d\theta=1$.

For a map $F:{\mathbb{T}^{d}}\to{\mathbb{T}^{d}}$, the Birkhoff average of a function $f:{\mathbb{T}^{d}}\to{\mathbb{R}}$ along the trajectory $\theta_{n}=F^{n}\theta_{0}$ is

[TABLE]

Theorem 1.2** (Quasiperiodic case of the Birkhoff Ergodic Theorem [6])**

Let $F:{\mathbb{T}^{d}}\to{\mathbb{T}^{d}}$ satisfy Eq. 1 where $\rho\in{\mathbb{T}^{d}}$ is irrational. Let $\mu$ be Lebesgue measure on ${\mathbb{T}^{d}}$. Then for every**The ergodic theorem for general ergodic maps replaces “for every” with “for almost every” but for quasiperiodic maps the “almost” can be omitted. initial $\theta_{0}\in{\mathbb{T}^{d}}$, $\lim_{N\to\infty}B_{N}(f)(\theta_{0})$ exists and equals $\int fd\mu$.*

The Weighted Birkhoff Averaging method ($\mbox{WB}^{[p]}_{N}$). We have recently developed a method for speeding up the convergence of the Birkhoff sum in Theorem 1.2 through introducing a $C^{\infty}$ weighting function by orders of magnitude when the process is quasiperiodic and the function $f$ is $C^{\infty}$, a method we describe in  [7, 8, 9]. In [9] it is proved that the limit of using $\mbox{WB}^{[p]}_{N}$ is the same as Birkhoff’s limit.

Weighted Birkhoff ($\mbox{WB}^{[p]}_{N}$) average of $f$ is calculated by

[TABLE]

where the $C^{\infty}$ weighting function $w$ is chosen as

[TABLE]

In our calculations of the rotation rates, we use $p=1$ or $2$. See in particular [7] for details and a discussion of how the method relates to other approaches. Note that essentially the same weight function as for $p=1$ case is discussed by Laskar [10] in the Remark 2 of the Annex, but he does not implement it.

Delay Coordinate Embeddings. For manifolds $M_{1}$ and $M_{2}$, a map $h:M_{1}\to M_{2}$ is an embedding (of $M_{1}$) if $h$ is a diffeomorphism of $M_{1}$ onto its image $h(M_{1})$. In particular the map must be one-to-one.

Let $\psi:{\mathbb{T}^{d}}\to M_{0}$ be $C^{2}$ where $M_{0}$ is a smooth manifold of dimension $D$. In our applications below, $\psi$ is either $\phi:{\mathbb{T}^{d}}\to S^{1}$ or $\gamma:{\mathbb{T}^{d}}\to{\mathbb{R}}^{2}$. While $d$ is the dimension of the domain ${\mathbb{T}^{d}}$ of $\psi$, $D$ is the dimension of the range.

For a positive integer $K$, define $\Theta_{K}^{\psi}:{\mathbb{T}^{d}}\to(M_{0})^{K}$ as

[TABLE]

$K$ is referred to as the delay number and is more precisely the number of coordinates used in defining $\Theta$. See Discussion, Section 5. In the theorem below, if $K=1$, we have a Whitney-type embedding theorem, or if $D=1$, a Takens-like result.

In order to include both of the projection maps $\phi:{\mathbb{T}^{d}}\to S^{1}$ and $\gamma:{\mathbb{T}^{d}}\to{\mathbb{R}}^{2}$, we introduce the more general notation $\psi:{\mathbb{T}^{d}}\to M_{0}$, where the manifold $M_{0}$ is $D$-dimensional. Hence $\phi$ or $\gamma$ can be substituted for $\psi$ with $D=1$ or $2$, respectively.

Theorem 1.3

[Special case of Theorem 2.5 in [11]] Let $M_{0}$ be a smooth $D$-dimensional manifold. Assume $F:{\mathbb{T}^{d}}\to{\mathbb{T}^{d}}$ is quasiperiodic where $F$ is given in Eq. 1 and $\rho$ is irrational. Assume

[TABLE]

Then for almost every $C^{2}$ function $\psi:{\mathbb{T}^{d}}\to M_{0}$, the map $\Theta:{\mathbb{T}^{d}}\to M_{0}^{K}$ is an embedding of ${\mathbb{T}^{d}}$.

While this result gives a lower bound on the delay number $K$, it is often convenient to choose $K$ much larger than required.

Define $\Gamma=\Gamma_{K}^{\psi}:{\mathbb{T}^{d}}\to(M_{0})^{K}\times{\mathbb{R}}$ as

[TABLE]

where $\hat{\Delta}$ is given in Eq. 7. See Fig. 5. The following corollary follows immediately from Theorem 1.3.

Corollary 1.4

Assume the hypotheses of Theorem 1.3. Then for almost every smooth ($C^{2}$) function $\psi:{\mathbb{T}^{d}}\to M_{0}$, the map $\Gamma_{K}^{\psi}:{\mathbb{T}^{d}}\to M_{0}^{K}\times{\mathbb{R}}$ is an embedding of ${\mathbb{T}^{d}}$.

Theorem 2.1 explains how this result is used when we have the image of a trajectory such as $(\gamma(\theta_{n}))_{n=0}^{N-1}$ – when $N$ is sufficiently large.

Comparison to previous work. We have written previously about computation of rotation rate in the papers [7, 8, 9]. A complete streamlined method for the case $d=1$ is provided in Section 2; the Embedding continuation method is announced in [8], but this is the first paper in which it is explained. In addition, this paper is the first time that we have applied our methods to cases where $d>1$. While we used the example (CR3BP) in [7], there we used a Poincaré return map whereas here in Section 4.2 no return map is used. We discuss the connections to our work with [3, 4, 5] in the subsequent sections of the paper. Those papers do investigate the Babylonian Problem, starting with only a set of iterates for a single finite length forward trajectory with the goal of finding a rotation number for some projection of a torus.

The investigation of quasiperiodic orbits is considered in [12, 13, 14, 15]. The approach in these papers assumes access to the full form of the original defining equations. Those papers are not investigating the Babylonian Problem.

Our paper proceeds as follows. We give a detailed description of our Embedding continuation method and an algorithm to implement it, in Section 2. Theorem 2.1 gives a proof of convergence of our method. In Section 3, we illustrate our methods using two one-dimensional examples ($d=1$). We refer to these as the fish map (introduced by Luque and Villanueva [5]) and the flower map, based on the shapes of the graphs. In Section 4 we give two-dimensional ($d=2$) examples of maps for which we explore the difficulty of determining their rotation rates about a reference point. We end in Section 5 with a discussion.

2 Embedding continuation method.

We have established that there is a lift $\hat{\Delta}$ of $\Delta$ and that $\Theta$ and $\Gamma_{0}:=\Gamma$ are embeddings of ${\mathbb{T}^{d}}$ for almost every $\psi$. We will assume in this section that $\psi$ has been chosen so that $\Theta$ and $\Gamma_{0}$ are embeddings.

If we are given the image of a trajectory, either $\phi(\theta_{n})$ or $\gamma(\theta_{n})$, we do not yet know what the corresponding $\hat{\Delta}_{n}$ is. In this section, we describe how we find the lift of a map using our Embedding continuation method. A schematic of these ideas is depicted in Fig. 5.

A major difficulty in evaluating $\rho_{\phi}$ is that $\hat{\Delta}(\theta_{n})$ is not known even though $\hat{\Delta}(\theta)\bmod 1=\Delta(\theta)$. This is because $\hat{\Delta}(\theta)\in\mathbb{R}$ is a lift of $\Delta(\theta)\in S^{1}$; i.e., they differ by an (unknown) integer $m(\theta):=\hat{\Delta}(\theta)-\Delta(\theta)$. The key fact is that from its definition, $\hat{\Delta}(\theta)$ is continuous and since it is defined on a compact set it is uniformly continuous. We describe in Steps 1 and 2 below how to choose the integer part of $\hat{\Delta}(\theta_{n})$ consistently, that is, so that $\hat{\Delta}(\theta_{n})$ is continuous on $S^{1}$. They collectively constitute our Embedding continuation method.

Step 1. The embedding. Let $N$ be given; in practice we usually use $N\sim 10^{5}$ or $10^{6}$ if $d=1$. Choose the delay number $K$ so that $2d+1\leq KD$. Recall that $\psi$ is either $\gamma$ or $\psi$ in our applications. Since $\psi(\theta)\in M_{0}$, we have $\Theta(\theta)\in M_{0}^{K}$. By our version of the Takens Embedding Theorems, Theorem 1.3, if $2d+1\leq KD$, then for almost every smooth function $\phi$, the map $\Theta$ is an embedding. In particular, there are no self intersections i.e., if $\Theta(\theta_{1})=\Theta(\theta_{2})$, then $\theta_{1}=\theta_{2}$. That implies $\Gamma$ defined by Eq. 17 is also an embedding of ${\mathbb{T}^{d}}$. We point out above that having an embedding guarantees that there are no self intersections, but there can be points far apart whose images are close to each other, and we try to avoid that by choosing $K$ large.

Denote $\mathbb{U}$ ={$(\Theta(\theta),\Delta(\theta)+k)$ : for all $\theta\in{\mathbb{T}^{d}}$ and all $k\in\mathbb{Z}$}.

The minimum distance $\epsilon$ between components of $\mathbb{U}$. For each $j\in\mathbb{Z}$, define

[TABLE]

and write $\Gamma_{j}:=\Gamma_{j}({\mathbb{T}^{d}}).$ Of course $\Gamma_{0}=\Gamma({\mathbb{T}^{d}})$. Then $\mathbb{U}$ is the union of all $\Gamma_{j}$. These sets are “vertical” translates of $\Gamma({\mathbb{T}^{d}})$ by an integer $j$, i.e. translates in the second coordinate. These are all disjoint from each other (since $\Theta({\mathbb{T}^{d}})$ is assumed to be an embedding). See Fig. 5 for an illustration.

Define

[TABLE]

where $\|\cdot\|$ is the Euclidean norm on $\mathbb{R}^{2K+1}$.

Then $\epsilon>0$ and $\epsilon$ is the minimum distance between points on different components of $\mathbb{U}$. In general $\epsilon$ is hard to compute from just the time series $\psi_{n}:=\psi(\theta_{n})$, so we have to fix a threshold $\delta>0$, assuming that $\delta<\epsilon$. Then if $p_{1},p_{2}\in\mathbb{U}$ and $\|p_{1}-p_{2}\|<\delta$, it follows that $p_{1}$ and $p_{2}$ are in the same component of $\mathbb{U}.$

The choice of the delay number $K$. It is important to note that this separation distance $\epsilon$ depends on the choice of $K$ and we observe that increasing $K$ increases $\epsilon$, so that while Theorem 1.3 guarantees we have an embedding and therefore $\epsilon>0$, this $\epsilon$ may be small. That might make it necessary to have a very large $N$, so instead we choose $K$ much larger than the theorem requires.

Step 2. Extending by $\delta$-continuation. Write $\Theta_{n}:=\Theta(n\rho)$. The goal is to choose integers $m_{n}$ so that all of the points $(\Theta_{n},\Delta_{n}+m_{n})$ for $n=0,\cdots,N-1$ are in the same component. This may be impossible if $N$ is not large enough. The point $(\Theta_{0},\Delta_{0})$ is in some component and we choose $m_{0}=0$ which determines a component. Let $\mathbb{A}$ be the set of $n\in\{0,\cdots,N-1\}$ for which $m_{n}$ has an assigned value. This set $\mathbb{A}$ changes as the calculation proceeds. Initially $m_{n}$ is assigned only for $n=0$ so at this point in the calculation the set $\mathbb{A}$ contains only [math]. Each time we assign a value to some $m_{n}$, that subscript $n$ becomes an element of $\mathbb{A}$. If there is an $n_{1}\in\mathbb{A}$ and an $n_{2}\notin\mathbb{A}$ and an integer $k$ such that

[TABLE]

then the two points are in the same component and we assign $m_{n_{2}}=k$, which adds one element, $n_{2}$ to the set $\mathbb{A}$. Keep repeating this process (if possible) until all $\{m_{n}\}_{n=0}^{N-1}$ are assigned. (We will make this procedure precise in Prop. 2.2.)

For $N$ sufficiently large, all can be assigned values, in which case we define $\hat{\Delta}_{n}=\Delta_{n}+m_{n}$ for all $n\in\{0,\cdots,N-1\}$. Define

[TABLE]

In the following theorem, we want $\delta<\epsilon$ where $\epsilon$ is in Eq. 18.

Theorem 2.1

For a $d$-quasiperiodic map assume $\Theta$ is an embedding. Given a map $\psi$, for $\delta$ sufficiently small, for all sufficiently large $N$ (depending on $\delta$), the above value $\rho_{\psi}^{N}$ is well defined (since all $m_{n}$ are defined), and

[TABLE]

2.1 Continuation algorithm: long chains of little steps on ${\mathbb{T}^{d}}$.

To determine all $\hat{\Delta}(\theta_{n})$ for all $n\in\{0,\cdots,n_{N-1}\}$, we begin knowing only $\hat{\Delta}(\theta_{0})$. Knowledge of $\hat{\Delta}$ can spread like an infection, transmitted between nearby $\theta_{n}$. The epidemic is spread only in little steps. The goal is to describe a continuation algorithm that identifies chains of $n_{j}$’s starting from $n_{j}=0$ and can reach every $n_{j}\in\{0,\cdots,n_{N-1}\}$.

To define “little step” we need a metric. Let $d(\cdot,\cdot)$ be a metric on ${\mathbb{T}^{d}}$ which is translation invariant, i.e. $d(x,y)=d(x+z,y+z)$ for all $x,y,z\in{\mathbb{T}^{d}}$. Furthermore for all $x=(x_{1},\cdots,x_{d})$ where all $|x_{j}|<0.5$, let $d(x,0)=\sum_{j}|x_{j}|$ (where here $d$ denotes the distance on the “$d$”-dimensional torus).

According to Theorem 1.3, $\Theta$ is almost always an embedding of the (rigid-rotation) torus into a higher dimensional space, so we can reasonably assume the following hypothesis.

$H_{1}$. $\Theta$ is an embedding. (Hence $\Gamma$ is also an embedding by Cor. 1.4.)

In this section we will assume $\epsilon$ is given by Eq. 18. Then $(\Theta,\hat{\Delta})({\mathbb{T}^{d}})$ is a smooth graph over $\Theta({\mathbb{T}^{d}})$. Hence if two points $\theta_{1}$ and $\theta_{2}$ in the are sufficiently close to each other, their images in $(\Theta,\hat{\Delta})({\mathbb{T}^{d}})$ will be less $\delta$ apart. That is given $\delta$ there is a $\delta_{1}>0$ such that ($d(\theta_{n_{1}},\theta_{n_{2}})<\delta_{1}$) implies Ineq. 19 will be satisfied. Hence, if $m_{n_{1}}$ has been assigned, and $m_{n_{2}}$ has not, then we will now be able to assign it a value.

We say $(n_{0},n_{1},\cdots,n_{k})$ is an $N$-$\delta_{1}$-chain from $\theta_{n_{0}}$ to $\theta_{n_{k}}$ if $n_{j}\in\{0,\cdots,N-1\}$ for all $j\in\{0,\dots,k-1\}$ and $d(\theta_{n_{j}},\theta_{n_{j+1}})<\delta_{1}$ for all $j\in\{0,\cdots,k-2\}.$

Proposition 2.2

(Long Chains of Little Steps).* Let $F:{\mathbb{T}^{d}}\to{\mathbb{T}^{d}}$ be the rigid rotation with rotation vector $\rho$ with a dense trajectory. For $\delta_{1}>0$, there is $N>0$ such that for every $n\in\{0,\cdots,N-1\}$ there is a $N$-$\delta_{1}$-${\bf chain}$ from $\theta_{0}$ to $\theta_{n}$.*

The following corollary interprets this proposition in terms of lifts and its proof is immediate.

Corollary 2.3

Assume $H_{1}$. Assume $\delta_{1}>0$ is such that $d(\theta_{n_{1}},\theta_{n_{2}})<\delta_{1}$ implies Ineq. 19. Then, since $m_{0}=0$, all $m_{n}$ can be determined. Write $\hat{\Delta}_{j}=\Delta_{j}+m_{j}$. Then all $\hat{\Delta}_{j}$ are in the same lift of $\Delta$. In other words, $(\Theta^{\phi}_{K,j},\hat{\Delta}_{j})$ are all in the same component of $\mathbb{U}$ where $\mathbb{U}$ is defined in Fig. 5.

To sketch a proof of the Proposition, we need the following fact. It is an elementary fact whose proof we leave to the reader.

Given $\delta_{1}>0$, there exists an $N$ with the following property.

$H_{2}$. There exist integers $0<\sigma_{1}<\sigma_{2}<\cdots<\sigma_{P}$ for some integer $P>1$ that $(i)$ the $\sigma_{j}$ are relatively prime (i.e., the greatest common factor of all $\sigma_{j}$ is $1$) and $(ii)$ $\theta_{\sigma_{j}}$ are within $\delta_{1}$ of $\theta_{0}$. Furthermore, $\sigma_{1}+\sigma_{P}<N$.

It is always possible to choose $N$ sufficiently large that $P=2$ in $H_{2}$; however, we might not want to choose such a large $N$, and we might be satisfied with having $P>2$.

An example of a pair $\theta_{\sigma_{1}}$ and $\theta_{\sigma_{2}}$ with relatively prime subscripts in dimension $d=1$. The algorithm for creating chains does not depend on the dimension $d$. Here we let $d=1$ and $\rho=\pi-3$ and $N=200$ and $\delta_{1}=0.01$ (where $d(0,x)=|x|$ for $x$ close to [math]). Then we can choose $\sigma_{1}=7$ and $\sigma_{2}=113$ since $\pi 7-22\approx 0.008$ and $355-\pi{113}\approx 0.00003$ so $\theta_{7}$ and $\theta_{113}$ are within $\delta_{1}$ of [math] and the subscripts $7$ and $113$ are relatively prime. We can reach every subscript in $\{0,\cdots,N-1\}$ by starting from [math] taking little steps, either increasing the subscript by $113$ or decreasing it by $7$, all the while staying between [math] and $N$, taking steps of size less than $\delta_{1}$.

An example of a pair $\theta_{\sigma_{1}}$ and $\theta_{\sigma_{2}}$ with relatively prime subscripts in dimension $d=2$. See Fig. 6. On the left where $N=100$, a chain is shown from [math] taking only steps of either $+93$ or $-4$. Both are within $\delta_{1}=0.13$ of [math]. It would work equally well to take only steps of $-93$ or $+4$. When $N=20,000$ on the right, there are two relatively prime subscripts $4109$ and $11700$ whose $\theta$ values are within $\delta_{1}=0.011$ of [math].

Proof of Proposition 2.2. We now describe why each $\theta_{n}$ can be reached by a chain starting from $\theta_{0}$.

We assume that for the given $\delta_{1}$, the $N$ and $\sigma_{j}$ have been chosen so that $(i)$ and $(ii)$ in $H_{2}$ are satisfied.

Let $B:=B(\delta_{1})$ denote the $\delta_{1}$ neighborhood of $\theta_{0}=0$. First we assume the number $P$ of $\sigma_{j}$ satisfies $P=2$, so $\theta_{\sigma_{1}}$ and $\theta_{\sigma_{2}}$ are in $B$ and their subscripts are in $\{0,\cdots,N-1\}$ and are relatively prime.

For non-negative integers $a_{1},a_{2}$, write

[TABLE]

Suppose $a_{1}$ and $a_{2}$ are such that

[TABLE]

Since $0\leq\sigma_{1}+\sigma_{2}<N$, we can either increase $a_{1}$ or $a_{2}$ by $1$ (thereby decreasing $a_{2}\sigma_{2}-a_{1}\sigma_{1}$ by $\sigma_{1}$ or increasing it by $\sigma_{2}$, respectively) and still have Condition 20 satisfied.

Notice that the distance from $[[a_{1},a_{2}]]$ to $[[a_{1}\pm 1,a_{2}]]$ or $[[a_{1},a_{2}\pm 1]]$ is less than $\delta_{1}$. That is, changing either $a_{1}$ or $a_{2}$ by $1$ moves $[[a_{1},a_{2}]]$ by less than $\delta_{1}$.

The key step of the proof is the following.

Algorithm. We choose a chain, a finite sequence $(\theta_{j})$ of such points, each of the form $[[a_{1},a_{2}]]$ as follows. Our algorithm begins at $\theta_{0}$ with $a_{1}=a_{2}=0.$

$A_{1}.$ Increase $a_{2}$ by $1$ provided the subscript remains non negative; otherwise increase $a_{1}$ by $1$. Repeat the process. Eventually the subscript returns to [math] (with $a_{1}=\sigma_{2}$ and $a_{2}=\sigma_{1}$. We have thereby created a chain of points on the torus, but we most likely have not encountered all the $\theta_{j}$.

Next,

$A_{2}.$ for each point in that chain, increase $a_{2}$ by $1$, and repeat as long as Condition 20 satisfied. This process yields $\theta_{n}$ for every $n\in\{0,\cdots,N-1\}$.

If $P=2$, we are done. When $P>2$, the greatest common factor, denoted $\Psi_{2}$, of $\sigma_{1}$ and $\sigma_{2}$ is greater than $1$. Then the above procedure reaches all points with subscripts divisible by $\Psi_{2}$ and no others. The next step is essentially the same as $A_{2}$ except that steps are taken by adding $\sigma_{3}$ to the subscript; that is,

$A_{3}$. Repeatedly add $\sigma_{3}$ to the subscript, as long as it remains less than $N$.

Taking all of those points and taking a small step for each by adding or subtracting $\sigma_{3}$ repeatedly will reach all points whose subscript is divisible by $\Psi_{3}:=$ the greatest common divisor of $\sigma_{1}$, $\sigma_{2}$, and $\sigma_{3}$.

$A_{j}$. For each point that has been found so far, repeatedly add $\sigma_{j}$ to the subscript as long as it remains less than $N$.

Eventually all $\theta_{n}$ for $0<n<N$ will be reached.

2.2 A dense set of equivalent representations for each rotation vector

While the definition of quasiperiodicity requires that the map has some coordinate system that turns the map into Eq. 1, that requirement by itself does not determine the coordinates of $\rho$. Fixing a coordinate system allows us to write $\rho=(\rho_{1},\cdots,\rho_{d})$ We have defined $\rho$ in Eq. 1 in terms of a given coordinate system. Let $\bar{\theta}=A\theta$ where $\bar{\theta}\in{\mathbb{T}^{d}}$ and $A$ is a unimodular transformation, that is an integer-entried matrix with determinant $|detA|=1$, then in this new coordinate system Eq. 1 becomes

[TABLE]

which is essentially Eq. 5. Hence $A\rho$ is also a rotation vector for the same torus map. Below we show we have a dense set of rotation vector representations.

Let $\mathcal{S}$ denote the set of integer-entried $d\times d$ matrices with determinant $\pm 1$. Observe that for any $B\in\mathcal{S}$, $B^{-1}\in\mathcal{S}.$ A matrix in $\mathcal{S}$ can be viewed as a change of variables on the torus, since it preserves volume. Therefore we call a vector $\tilde{\rho}\in{\mathbb{R}}^{d}$ a rotation representation of $\rho\in{\mathbb{R}}^{d}$ if $\tilde{\rho}=A\rho$ for some $A\in\mathcal{S}$. We ask: When the vector $\rho$ is irrational, what are all the possible rotation vectors (i.e., rotation representations), assuming $A\in\mathcal{S}$?

Proposition 2.4

Assume dimension $d\geq 2$. For an irrational rotation vector $\rho$, the set of its rotation vector representations is $S\rho$ (i.e., $\{A\rho:A\in S\}$), and $S\rho\bmod 1$ is dense in ${\mathbb{T}^{d}}.$

Proof. To simplify notation we prove only the case of $d=2$. The proof for $d>2$ is analogous. See [16]. Write $\rho=(\rho_{1},\rho_{2})$. Note that the matrices B_{m}:=\left(\begin{array}[]{c}1~{}~{}m\\ 0~{}~{}1\end{array}\right) and C_{k}:=\left(\begin{array}[]{c}1~{}~{}0\\ k~{}~{}1\end{array}\right) are in $\mathcal{S}$ for all integers $m$ and $k$, as is $A:=B_{m}C_{k}$. Then the vectors $(\rho_{1},y_{k})=C_{k}(\rho_{1},\rho_{2})\mod~{}1$ are vertical translates (translates in the direction $(0,1)$) of $(\rho_{1},\rho_{2})\mod~{}1$, where $\{y_{k}\}$ is a dense set in $S^{1}$. When we similarly apply $B_{m}$ for all $m$ to each $(\rho_{1},y_{k})$ we obtain a dense set of horizontal translates of $(\rho_{1},y_{k})$ and thereby obtain a dense set in ${\mathbb{T}^{2}}$. Every coordinate of every point in that dense set is of the form $k_{1}\rho_{1}+k_{2}\rho_{2}\mod~{}1$ where $k_{1}$ and $k_{2}$ are integers.

3 Examples of one-dimensional quasiperiodicity.

In this section, we give a detailed explanation of how we compute rotation rates for quasiperiodic maps on one-dimensional tori. For the first example computation of the rotation number is easy and straight forward while in the second it is sufficiently hard that we need our method. The pair of examples makes it clear when our method should be used.

One advantage of the examples below is that we know $\rho$ and therefore we can compare it with the computed rotation rates.

Example 1. The fish map. Luque and Villanueva [5] addressed the case of a quasiperiodic planar curve $\gamma:S^{1}\to\mathbb{C}$ and introduced what we call the fish map, depicted in the left panel of Fig. 1. Let

[TABLE]

where $z=z(\theta):=e^{i2\pi\theta}$ and $\hat{\gamma}_{-1}:=1.4-2i,$ $\hat{\gamma}_{0}:=4.1+1.34i,$ $\hat{\gamma}_{1}:=-2+2.412i,$ $\hat{\gamma_{2}}:=-2.5-1.752i$. (See Fig. 5 and Eq. 31 in [5]). They chose the rotation rate $\rho=(\sqrt{5}-1)/2\approx 0.618$ for the trajectory $\gamma_{n}=\gamma(n\rho)$ for $n=0,1,\cdots$ so we also use that $\rho$. The method in [5] requires a step of unfolding $\gamma$, which our method bypasses. We measure angles with respect to $P_{1}=8.25+4.4i$, where the winding number $|W(P_{1})|=1$.

Example 2. The flower map. We have created an example, the flower map in Fig. 1, right, to be more challenging than the fish. Let

[TABLE]

We use the same $\rho=(\sqrt{5}-1)/2$ as above. The choice of a reference point $P_{1}$ for which $|W(P_{1})|=1$ is shown in the right panel of Fig. 1. For our computations, we use $P=P_{1}:=0.5+1.5i$. Points $P_{j}$ with $|W(P_{j})|=j$ for $j=1,2,3,6$ are also shown. For $\rho_{\gamma_{6}}$, the rotation rate of $\gamma_{6}$, to yield $\rho$ or $1-\rho\bmod 1$ is essential to choose a point $P$ where $|W(P)|=1$. In this example the values of $\Delta_{n}$ are dense in $S^{1}$, and $\max_{\theta}\hat{\Delta}(\theta)-\min_{\theta}\hat{\Delta}(\theta)\approx 1.2$.

For both examples, Fig. 1 shows two successive iterates $\gamma_{n}$ and $\gamma_{n+1}$, and the angle $\Delta_{n}$ between these two iterates, computed with respect to a reference point $P_{1}$. It was computed by finding $\phi_{n}$, the angle of $\gamma_{n}$ with respect to $P_{1}$ as in Eq. 12. Using this, $\Delta_{n}=\phi_{n+1}-\phi_{n}\in[0,1)\approx S^{1}$. On the left, in the fish map case, if we choose $\hat{\Delta}_{0}:=\Delta_{0}$ (or alternatively $:=\Delta_{0}+m$ for some $m$), then we have selected the component in which all $\hat{\Delta}_{n}$ must lie. This is what is referred to below as the easy case. Choose some $k$, write $J_{k}:=[a,b]$. Choose $m_{n}$ is the integer for which $\Delta_{n}+m_{n}\in J_{k}$. It is not as easy to do this for the flower map on the right. Fig. 3 right shows that the possible lifts when plotted against $\phi$ form a tangled mess which does not resolve into bounded components, while when plotted against $\theta$ we obtain components that are diffeomorphic to $S^{1}$.

Figs.  3 and  4 show the possible lift values $\hat{\Delta}_{n}$ of the angle difference $\Delta_{n}$ plotted with respect to angle $\theta$ in Fig.  3 and $\phi$ in Fig. 4. For the fish map on the left, we see that we can set $[a,b]\approx[0.18+k,1.05+k]$ for any integer $k$. Furthermore, we investigated the rotation rate of the signal viewed from $P_{1}=7+4i$. Using the Weighted Birkhoff Average, we observe that the deviations of the approximate rotation from $\rho$ falls below $10^{-30}$ when the iteration number exceeds $N=20,000$, and since we know the actual rotation rate, we can report that the error in the rotation rate is then below $10^{-30}$. Once we have found a proper lift for the flower map, we can do the same procedure. The next section explains how we go about finding a lift in this more complicated case.

4 Higher-dimensional quasiperiodic examples

We develop a higher-dimensional method to compute the rotation vector $\rho$ purely from knowledge of the sequence $\theta_{n+1}:=F(\theta_{n}).$ The question of how to compute the rotation vector is actually two questions. Question 1: If we compute a rotation vector, what are the possible values? Question 2: How do we compute any of the possible values for the rotation vector in difficult cases? Figs 8, 9, and 10, demonstrate that like in one dimension, in $d$ dimensions we are able to use $d$ independent planar projections combined with a higher-dimensional version of our Embedding continuation method in order to find a lift, each projection leading to one component of a $d$-dimensional rotation vector. In fact, these rotation vectors are not unique. In this section, we give a detailed discussion of our higher-dimensional method, describing the possible values we can achieve in calculating a rotation vector. We then illustrate our method for three examples: the fish torus, the flower torus, and the restricted three-body problem.

4.1 Two examples in a higher dimension: fish and flower tori ${\mathbb{T}^{2}}$.

We use the fish and flower maps from the previous section in order to create 2-dimensional torus maps. We will explore the problem of computing rotation rates for these examples where we know the rotation rates for the quasiperiodic maps. Let $\rho:=(\sqrt{5}-1)/2$ and $\phi:=\sqrt{3}/2$, and define

[TABLE]

Let $\gamma$ be either the fish or the flower map defined in the previous section. Define the torus-version $f_{T}$ of the $\gamma$ map(s) as follows. Let $Re(\cdot)$ and $Im(\cdot)$ denote the real and imaginary components of a complex number, and let $f_{T}:{\mathbb{T}^{2}}\to{\mathbb{R}}^{3}$. Write $f_{T}(\theta_{n},y_{n})=(f_{1},f_{2},f_{3})(\theta_{n},y_{n})$, where

[TABLE]

The “+2” is just for convenience so that the torus image can wrap around the origin rather than having to wrap it around some other point. For each $\gamma$, the map $f_{T}$ takes a quasiperiodic trajectory into ${\mathbb{R}}^{3}$.

Two projections of a torus for two rotation rates. Figure 8 shows two independent projections of $f_{T}$ to ${\mathbb{R}}^{2}$. For the first rotation rate, we project $f_{T}$ to $(f_{1},f_{2})$ in the plane. Then we measure the angle $\phi$ from a reference point $P$ which is not in the image of the torus. In particular, $P=(0,1.5)$ for the fish torus, and $(0,0.1)$ for the flower torus. For both maps, this projection gives a rotation rate of $\phi/2\pi$ (the denominator $2\pi$ comes from the fact that we are measuring angles in $[0,1]$).

For a second rotation rate, let $R_{\alpha}$ be the rotation matrix that tilts by angle $\alpha=0.05\pi$ in the $f_{2}-f_{3}$ plane. Namely

[TABLE]

Set

[TABLE]

Define $r=\sqrt{h_{1}^{2}+h_{2}^{2}}$. Then our projection is to the value $(r,f_{3})$. We measure the angle of this projection relative to the point $(8.25,4.4)$ for the fish torus, and $(2.6,1.4)$ for the flower torus. For both maps, this projection gives rotation rates of $1-\phi/(2\pi)$ and $1-\rho$. Why the tilt by $0.05\pi$ rather than use value of $r$ with respect to the original coordinates? Because without the tilt (i.e. $\alpha=0$), the projection would be a curve rather than a thick strip, which would not give a true test of our Embedding continuation method in two dimensions.

In both cases, we get a map whose image has at least one hole (in which the winding number $=\pm 1$), and we can measure angles $\phi$ and angle differences $\Delta$ compared to a point inside one of the holes, as long as the torus has a winding number $|W(p)|=1$ with respect to points in this hole. Thus just as for the one-dimensional case, we compute the lift, and then compute the rotation rate for these two different projections.

Figures 9 and 10 show the original values of the angle difference, and the computed lift, respectively. Note that fish torus lift is easy to compute while the flower torus requires and embedding.

As mentioned in Section 1, rather than using Birkhoff Averages, we achieve more rapid convergence using our Weighted Birkhoff Average, denoted $\mbox{WB}^{[p]}_{N}$ in Eq.14. Define the $\rho$ approximation $\rho_{N}:=\mbox{WB}^{[1]}_{N}(\hat{\Delta}_{n})$, when $p=1$. Fluctuations in $\rho_{N}$ fall below $10^{-30}$ for $N>20,000$. Since we know the actual rotation rate, we can report that the error $|\rho-\rho_{N}|$ is then below $10^{-30}$.

4.2 The circular planar restricted three-body problem (CR3BP)

CR3BP is an idealized model of the motion of a planet, a moon, and an asteroid governed by Newtonian mechanics Poincaré [17, 18] introduced his method of return maps using this model. In particular, we consider a circular planar three-body problem consisting of two massive bodies (“planet” and a large “moon”) moving in circles about their center of mass and a third body (“asteroid”) whose mass is infinitesimal, having no effect on the dynamics of the other two.

This model can also (simplistically) represent the Sun-Earth-Moon system discussed in the introduction though the parameter $\mu$ has to be changed, and the Moon is the body that is assumed to have negligible mass. All three travel in a plane.

We assume that the moon has mass $\mu$ and the planet mass is $1-\mu$ where $\mu=0.1$, and writing equations in rotating coordinates around the center of mass. Thus the planet remains fixed at $(q_{1},p_{1})=(-0.1,0)$, and the moon is fixed at $(q_{2},p_{2})=(0.9,0)$. In these coordinates, the satellite’s location and velocity are given by the generalized position vector $(q_{1},q_{2})$ and generalized velocity vector $(p_{1},p_{2})$.

Define the distance of the asteroid from the moon and planet are

[TABLE]

[TABLE]

The following function $H$ is a Hamiltonian (see [19] p.59 Eqs. 63-66) for this system

[TABLE]

where $p_{1}=\dot{q_{1}}-q_{2}$ and $p_{2}=\dot{q_{2}}+q_{1}$. We get the equations of motion from

[TABLE]

That is, the equations of motion are as follows:

[TABLE]

We measure angles as a fraction of a full rotation and not in terms of radians. The asteroid’s orbit in rotating coordinates is shown in Fig. 11. Here time is continuous so we can measure the total angle through which a trajectory travels, retaining the integer part. The first rotation rate $\rho^{*}_{\theta}$ of the asteroid’s orbit is its average rate of rotation about the planet, that is, the average rate of change of the angle $\theta$ measured from $(q_{1},q_{2})=(-0.1,0)$. We compute that $\rho^{*}_{\theta}=-2.497823504839344460408394$ rev/sec, that is, about -2.5 $\theta$-revolutions per unit time Fig. 13 (left) shows the error in convergence to the value 0.001$\rho^{*}_{\theta}$ Note that the rotation rate in the fixed coordinate frame is $\rho^{*}_{\theta}+1/(2\pi).$

The second rotation rate $\rho^{*}_{\phi}$ measures the oscillation in the distance $r$ from the planet. In particular, we project to the $(r,r^{\prime})$ plane, where $r^{\prime}:=dr/dt$. That is, define $r=\sqrt{(q_{1}+0.1)^{2}+q^{2}_{2}}$ and $r^{\prime}=\frac{dr}{dt}=((q_{1}+0.1)\frac{dq_{1}}{dt}+q_{2}\frac{dq_{2}}{dt})/r$, as shown in Fig. 12. The angle $\phi$ is measured from $(r,r^{\prime})=(0.15,0)$. The fast convergence to the value 0.001$\rho^{*}_{\phi}$ by the Weighted Birkhoff Average $\mbox{WB}^{[2]}_{N}$ in 14 is seen in Fig. 13 (right), where $\rho^{*}_{\phi}=-2.3380583953388194764236520190142509$ rev/sec. The period of time between perigees is the reciprocal, or about 0.43 time units.

We used the 8th-order Runge-Kutta method in Butcher [20] to compute trajectories of CR3BP with time steps of $h=2\times 10^{-5}$.

The meaning of rotation rates for the CR3BP. In [7], we investigated the same asteroid orbit of the CR3BP as is studied here, but instead of the continuous-time trajectory that lies on a two-dimensional torus as presented above, there we used a Poincaré map. The coordinates of the asteroid were recorded each time the asteroid crossed the line $q_{2}=0$ with $dq_{2}/dt>0$. In the cases we study, the map trajectory is a quasiperiodic trajectory on a closed curve. Hence there is only one rotation rate, a much simpler situation. Choosing a point inside the closed curve, we computed a rotation rate, namely the average angular rotation per iteration of the Poincaré map. The rotation rate $\rho^{*}_{P}$ per Poincaré map on the Poincaré surface $q_{2}=0$ (or equivalently, $\theta=0$) around $(q_{1},p_{1})=(-0.25,0)$ was computed as 0.0639617287574530971640777244014426955. We felt that the issues of rotation rates could be clarified if we computed the trajectory as a continuous orbit as we do here. The two rotation rates computed here $\rho^{*}_{\phi}$ and $\rho^{*}_{\theta}$ and our previous result $\rho^{*}_{P}$ bear the following relation to our previous results:

[TABLE]

[TABLE]

See the caption of Fig. 13. We solved the differential equation using an $8^{th}$-order Runge-Kutta method using quadruple precision. Both approaches are based on rotating coordinates, but there is another approach.

The orbit as a slowly rotating ellipse. The asteroid rotates about the planet at a rate of $\rho^{*}_{\theta}$ revolutions per unit time when viewed in the rotating coordinate in which the moon and planet are fixed. The sidereal rotation rate (as viewed in the coordinates of the fixed stars) is $\rho^{*}_{\theta}\ +1/(2\pi)$. We can think of the orbit as an approximate ellipse whose major axis rotates and even changes eccentricity (being more eccentric when the asteroid apogee is aligned with the planet moon axis).

Without the moon the asteroid orbit would be perfectly elliptical with its major axis fixed in position, but the moon causes the ellipse to rotate slowly. The angle $\phi(t)$ tells where the asteroid is on its roughly elliptical orbit; Fig. 12 shows that the apogee occurs when when the distance from the planet $r\sim 0.27$ and the perigee when $r\sim 0.05$, with some variation. The time between successive perigees averages $1/\rho^{*}_{\phi}$. Note that the difference in these rates satisfies

[TABLE]

Hence relative to the fixed stars, that is, in non-rotating coordinates, the asteroid’s ellipse’s major axis precesses slowly. Its apogee point returns to its original position (in non-rotating coordinates) after the asteroid passes through its apogee approximately $1639$ times.

5 Discussion and conclusions

What does it mean to ask for one or more rotation rates of the $d$-dimensional quasiperiodic map Eq. 1? One might expect that one should find $\rho$ or rather its coordinates. As we explain below and in Section 2.2, this is an ill-posed problem (especially for $d>1$. The Babylonians computed rotation rates for projections of the Moon’s trajectory onto the globe of fixed stars (as we have discussed in the Introduction). So we refer to their approach as the “Babylonian Problem”: computing rotation rates for a projection of a quasiperiodic process.

We have developed our Embedding continuation method for calculating the rotation rate for “almost every” Babylonian Problem, that is, for smooth projection $\psi$ from of a quasiperiodic dynamical system on ${\mathbb{T}^{d}}$. “Almost every” is in the sense of prevalence - and in practice there will be difficult cases especially since the number $N$ of interates needed increases as $d$ increases. Our Weighted Birkhoff Method of computing rotation numbers significantly shortens the computation time for computing rotation numbers, making our approach effective in practice. See the Introduction.

A key motivating difference between $d=1$ and $d>1$ is that in the higher-dimensional case, for the rotation vector $\rho\in{\mathbb{T}^{d}}$ there are infinitely many ways of choosing coordinates on ${\mathbb{T}^{d}}$ for the map Eq. 1. In Section 2.2 we show that the set of resulting coordinate representations $(\rho_{1},\cdots,\rho_{d})$ of $\rho$ are dense in ${\mathbb{T}^{d}}$. Every point $r$ in ${\mathbb{T}^{d}}$ is arbitrarily close to such representations. Hence instead of trying to find the coordinates of $\rho$, we have learned from the Babylonians, and we phrase our goals in terms of finding a rotation number $\rho_{\psi}$ (usually, $\rho_{\phi}$ or $\rho_{\gamma}$) for some projection from ${\mathbb{T}^{d}}$ into a one or two-dimensional space.

Even for $d=1$, there is some uncertainty for obtaining $\rho$ depending on the choice of orientation on $S^{1}$. We can obtain either $\rho~{}\bmod 1$ or $1-\rho~{}\bmod 1$.

In Section 4.2, we apply our method to the quasiperiodic torus occurring for a 4-dimensional circular restricted 3-body problem, depicted in Fig. 11. In particular, we explain the relationship between the two rotation rates obtained from the original differential equation system and the rotation rate which was previously obtained from the Poincaré map. The fact that the rotation rate of an asteroid will be different depending on whether on uses rotating coordinates or sidereal coordinates (in which the distant stars are fixed) is an example of how the rotation rate can depend on the projection.

Notes on delay coordinate embedding theorems. H. Whitney [21] showed that a topologically generic smooth map $\Gamma$ from a $d$-dimensional smooth compact manifold $M$ into ${\mathbb{R}}^{D}$ where $2d+1\leq D$ is a diffeomorphism on $M$; in particular the map $\Gamma:M\to F(M)$ is an embedding of $M$.

Sauer et al [11] modified Takens’ result in two ways. First, it replaced “topologically generic” by “almost every” (in the sense of “prevalence”) in Theorems 2.1 and 2.3 in [22]. See also [23]. For physical purposes “almost every” has significance while residual sets do not seem to. In this paper, in Theorem 1.3, we have adapted the “almost every” approach.

For completeness, we mention the second way [11] generalized Takens’ approach, even though this second way is not used here, because the sets we deal with are manifolds. The second way is that [11] allblack replaced the assumption that $M$ is a manifold by assuming only that $M\subset{\mathbb{R}}^{k}$ for some $k$ is an invariant set of some map and that $M$ has box dimension $boxdim(M)$ and $\Gamma$ is a mapping of a neighborhood of $M$ into ${\mathbb{R}}^{D}$ where $D>2\cdot boxdim(M)$. The great majority of citations to Takens [24] are for the case where $M$ is a chaotic attractor that is not a manifold so that Takens’ Theorem does not apply. Those papers actually use the results in [11], not in Takens’ [24]. One unusual aspect of our current paper is that we actually only need the case that Takens proved. Here $M$ is a quasiperiodic torus so it is a manifold.

The Takens Theorem also has assumptions that the set of periodic points $F:M\to M$ for some smooth map was in some sense small, in our case there are no periodic points so those assumptions are automatically satisfied. Hence we only state it in a special case needed here.

We have demonstrated that in one dimension, a rotation rate can be computed precisely with minimal ambiguity, but higher dimensional cases (${\mathbb{T}^{d}}$ with $d>1$) are more complicated. Projections into the plane yield rotation rates, but there are infinitely many topologically distinct ways to project a higher dimensional torus onto a circle, each of which yields a different rotation rate. This makes it important for the investigator to explain the meaning of any particular rotation rate. In fact, a rotation rate is a rate specifying an average change per unit time, where there can be considerable choice in the time units. To illustrate this point, we more carefully consider the CR3BP example with a focus on what the rotation rates tell us about the trajectories of an asteroid.

Acknowledgments. We would like to thank the referee for many helpful comments. YS was partially supported by JSPS KAKENHI grant 17K05360 and JST PRESTO grant JPMJPR16E5. ES was partially supported by NSF grant DMS-1407087. JY was partially supported by National Research Initiative Competitive grants 2009-35205-05209 and 2008-04049 from the USDA.