1 Introduction

The existence of invariant tori of maximal dimension is important in the study of the stability of dynamical systems; when the number of degrees of freedom is low, invariant tori act as barriers and confine the motion in regions of phase space.

In recent papers (Calleja et al. 2022, 2013), we have developed extremely effective methods to compute very accurately quasi-periodic solutions (KAM tori) of several models of the dissipative spin–orbit problem in celestial mechanics.

The spin–orbit problem describes the motion of an oblate satellite whose center of mass moves on a Keplerian orbit around a central planet. The satellite is not rigid, and its motion generates tidal friction. We consider both a model in which the friction is time dependent along the position in the orbit (the tides are different at different places in the orbit) and another one whose friction is constant (an average of the time dependent friction).

Both models depend on three parameters: the oblateness of the satellite, the orbital eccentricity and the dissipative factor. The limit of zero dissipation is a singular limit, since the long-term behaviors of the model change drastically from a conservative system to a dissipative one even if the dissipation is small.

For dissipative systems, to obtain a quasi-periodic orbit of a certain frequency, one has to adjust a parameter called drift. This is a big contrast with the conservative case which admit many quasi-periodic orbits. To obtain orbits with a given frequency in a conservative system, one just needs to choose the initial conditions and one does not need to adjust parameters.

The astronomically relevant values of the friction are typically small, so that the systems of astronomical interest are in the regime of singular perturbations. Frictions have crucial effects in the formation of planets over millions of years, but are very small in calculations over a few thousand years.

Working close to a singular limit, our numerical explorations, especially when we want to study the objects near the breakdown, require several elaborated numerical methods (see Sects. 4 and 5) as well as some theoretical underpinnings, e.g., condition numbers and theoretical arguments that indicate that the numerical calculations are credible. Both the efficient numerical algorithms and the a posteriori results take advantage of several identities that are obtained using the geometric properties of a (conformally symplectic) map.

The goal of this paper is to show empirically that our methods are able to continue very effectively the quasi-periodic solutions w.r.t. the parameters; at the same time, our methods allow to maintain a very high accuracy, even extremely close to parameter values where the solutions cease to exist. Our work gives results in different directions: (i) It provides a very good test of the methodology to maintain accuracy and reliability; (ii) it shows that the a posteriori theorems developed in Calleja et al. (2013) can be used in very practical situations and effectively cooperate with elaborated numerical methods; and (iii) it yields information of mathematical phenomena happening at breakdown of the invariant tori.

One of our main tools to ascertain the breakdown is that the Sobolev seminorms of the embedding functions describing the tori blowup when the parameters approach the breakdown. We also study the behavior of the stable and tangent bundles; in Calleja and Figueras (2012), it was found numerically that at the breakdown of tori in the dissipative standard map, one has simultaneously a loss of regularity of the torus and a bundle collapse between the stable and tangent bundles. (The center and the stable bundles get close in a dense set, even if they have to remain apart in a set of full measure, hence they have to become very oscillatory.) We note that there are rigorous results Calleja and de la Llave (2010), Calleja et al. (2013) which show that the algorithms can compute solutions arbitrarily close to breakdown and that the breakdown is signaled either by a blowup of Sobolev seminorms of the conjugacy or a loss of hyperbolicity. Note that in conformally symplectic systems, the stable Lyapunov exponent is bounded away from zero, so that the loss of hyperbolicity can only happen because the angles between the stable and tangent bundles are not bounded away from zero.

In contrast to previous observations in other models (mainly maps with very few harmonics), we have found that the Sobolev seminorms do not have an asymptotic power law and that the bundles do not oscillate wildly. They remain separated even if they lose regularity. There do not seem to occur standard scaling relations on the properties of the bundles.

The lack of standard scaling relations found in other models with few harmonics can be described in the renormalization language as saying that the spin–orbit problem is so far from previously studied maps that it is in a different universality class. Indeed, MacKay (1983) asked the question of clarifying the limits of the universality class.

It seems quite likely that some of the phenomena discovered here for the spin–orbit problem should be applicable to more models. These phenomena can be explained by certain properties of the dynamics of the renormalization group. Since the goal of this paper is to show how one can maintain accuracy and reliability up to breakdown in a concrete model of astronomical interest, we have not done a systematic exploration of similar models and we postpone the worthwhile goal of understanding dynamics of the renormalization group.

1.1 Some Comments on the Numerical Tools

To reach our results at the limit of breakdown, we have deployed the following numerical tools: a reduction of the spin–orbit problem to a time-1 map, a very high-order numerical integrator, jet transport to compute the derivatives, and extended precision arithmetic. Let us briefly comment on these tools.

1.1.1 Reduction to Time-1 Map

The reduction to a time-1 map encodes the KAM tori by a lower-dimensional object. The computation of the time-1 map decreases the memory usage but, in our case, it requires several integrations of an ODE, which is computationally expensive. However, the resulting code is very parallelizable since the integration of each orbit can be assigned to a thread (we have used this systematically).

Note that the tori—specially near breakdown—become rough along the transversal direction, while they remain smooth along the flow direction. Hence, it is natural to deal with these directions differently and to treat the integration along the flow (which remains always very smooth) differently than the KAM computation.

Another advantage of basing the results on the time-1 map is that it permits a rigorous justification of the comparison between the spin–orbit model (1) and the averaged version (4) of Section 2 (see Remark 1).

1.1.2 Using High-Order (Taylor) Methods, Jet Transport and Extended Precision; the taylor Package

As it turns out, the KAM method requires very high precision in the torus which can be reached by using extended precision in the arithmetic. The only way to achieve very high precision with reasonable time steps is to use a very high-order method (e.g., order 30).

Moreover, we also need to obtain derivatives of the time-1 map that, in principle, can be accomplished by integrating the variational equations.

All these requirements can be obtained simultaneously with a reasonable amount of programming time using, for instance, the public domain taylor package (Jorba and Zou 2005; Gimeno et al. 2022).

The taylor program takes as input the equations written in a simple symbolic form, and it outputs a C/C++ code that implements a Taylor stepper integrator. One of the options of taylor is to implement the arithmetic in an extended precision package by using, for instance, the MPFR package. If the user is not satisfied with the constant step, taylor allows to select and fine tuning several adaptive steps or even implement his/her own time stepping strategies.

We recall that the Taylor integrator (already used by I. Newton) consists in finding recursively the Taylor coefficients of the solution. The differential equation gives a recursion that prescribes the Taylor coefficients of the solution of order n as a function of the coefficients of lower orders. The initial conditions are the starting point of these recursions, which are easy to formulate when the differential equations are formed from elemental functions, such as \(\sin , \cos , \exp ,\log and \text {pow}\), and using composition and arithmetic operations.

Recently, an important feature of taylor described in Gimeno et al. (2022) is that it implements what is called jet transport and it provides high-order derivatives of the time-t map. Jet transport is the notion of combining a numerical integrator scheme with variables being multivariate polynomials instead of real numbers. As it is proved in Gimeno et al. (2023), the numerical solution is exactly the same as considering the original differential equation and its variational equations and therefore a stepsize strategy based on all the polynomial coefficients is equivalent to a strategy on its variational equations.

The upshot is that, using taylor with a programming effort comparable to typing the equations in TeX, one obtains an arbitrary-order solver for the equations and the variational equations, implemented in extended precision.

In our case, we computed the one-time map and its derivatives with a precision which is (conservatively) smaller than \(10^{-35}\) in a reasonable computer time. Clearly, stronger precisions could easily be accomplished given more computer resources.

1.1.3 Computation of KAM Tori Given the Map; Automatic Reducibility

Once we have the time-1 map very carefully computed, we turn on the problem of computing the KAM tori by a continuation scheme.

We use a Newton-like method based on automatic reducibility developed for symplectic mappings in de la Llave (2001), de la Llave et al. (2005) and for conformally symplectic systems in Calleja et al. (2013). The method seeks a parametrization of the torus (which is a one-dimensional object) that satisfies an invariance equation.

The algorithm takes advantage of several geometric identities to obtain a quadratically convergent method. The method has very low storage requirements and operation count (if the torus is discretized with N points, a step requires O(N) storage and \(O(N\log (N))\) operations). The algorithm (listed in Calleja et al. 2013) consists in about a dozen steps on objects. Using modern computer languages, all of those steps can be implemented in one line in a way that is independent of the level of truncation or accommodating extended precision arithmetic. Hence, significant parts of the problem can be reused.

1.1.4 A Posteriori Results

Since we are going to compute very close to breakdown, reliability is not obvious.

An important tool is the rigorous result of Calleja et al. (2013) stated in an a posteriori format; that is, given an approximate solution, if one computes some (explicit) condition numbers, then one can ensure that there is a true solution nearby. Hence, provided that we can compute solutions with good accuracy and keeping track of the condition numbers, we are confident of our calculations to breakdown. Explicit (but not completely rigorous since they ignored round-off) calculations of the condition numbers which give existence of the tori extremely close to the breakdown appear in Calleja et al. (2022).

We also note that the fact that we have a posteriori results of the right form tells us that a continuation method will progress unless the torus ceases to satisfy a few conditions. Hence, the method guarantees that our continuation methods will reach arbitrarily close to breakdown (if we could use enough discretization modes and enough precision). These methods were called accurate in Llave and Rana (1991). To show that a method is accurate, one needs to perform some mathematical analysis that ensures that the functions involved have domains that are completely covered by the discretization (e.g., if one discretizes by power series, one would need to show that the functions involved have domains which are circles).

In this paper, we have gone even further in parameters by relying on the standard methods of numerical analysis (e.g., reruns, increasing accuracy and checking it does not change).

1.1.5 Some Other Remarks

We conclude by mentioning that, in different problems, Figueras and Haro (2012) developed a similar methodology to compute reliably close to breakdown based on efficient algorithms backed up by a posteriori theorems. This was used in Figueras and Haro (2015) to discover different scenarios for breakdown of normal hyperbolicity.

As we will see, we have obtained what seems a new scenario of breakdown. Investigating the domain of universality will take us far from our main goal and will be postponed.

1.2 Organization of the Paper

The spin–orbit model is presented in Sect. 2, and its Poincaré map is described in Sect. 3. Our approach works with any Diophantine frequency, but in the numerical experiments we have considered two frequencies, namely the golden ratio and another irrational number related to the golden mean but closer to 1-the resonance. Both numbers are noble, i.e. the continued fraction expansion is 1 after a finite number of terms.

We have computed the so-called basins of rotation number, which give a qualitative picture of the dynamics as well as information on the value of the drift parameter, for given values of the dissipation (see Sect. 4).

The computation of the tori and their exploration near breakdown is presented in Sect. 5, while Sect. 6 analyzes different breakdown reasons.

The construction of the invariant bundles and an exploration of the angle between the stable and tangent bundles are given in Sect. 7. Moreover, Sect. 8 shows the width of analyticity domains of the tori and the angle between the invariant bundles.

The computation of the breakdown threshold through Sobolev’s criterion is given in Sect. 9.

The links with renormalization theory and the study of scale-invariant properties of the tori are studied in Sect. 10.

Finally, some conclusions are presented in Sect. 11.

2 The Spin–Orbit Problem with Tidal Effects

In this section, we describe the physical motivation for the model we will investigate; even if not strictly needed for the mathematical or numerical study, the physical interpretation motivates our questions and the ranges of the parameters involved. We will quickly review the most important features, but leave a more detailed derivation to the referenced literature (Peale 2005; Celletti 2010).

The spin–orbit problem describes a simplified model for the rotational dynamics of a satellite, say \(\mathcal {S}\), orbiting around a central planet, say \(\mathcal {P}\), and rotating around an internal spin-axis (Beletsky 2001; Celletti 1990b, a; Correia and Laskar 2004; Wisdom et al. 1984). We assume that the satellite is a triaxial ellipsoid with principal moments of inertia \(\mathcal {A}< \mathcal {B} < \mathcal {C}\). Moreover, we make the following simplifying assumptions:

  1. A1.

    The barycenter of the satellite \(\mathcal {S}\) moves on an elliptic Keplerian orbit with the planet \(\mathcal {P}\) in one focus; the Keplerian orbit is characterized by a semimajor axis a and an eccentricity \(e\).

  2. A2.

    The satellite spin-axis coincides with the direction of the smallest physical axis of the ellipsoid.

  3. A3.

    The spin-axis is perpendicular to the orbital plane.

  4. A4.

    The satellite is assumed to be non-rigid, hence subject to a tidal torque.

  5. A5.

    The tidal torque is a linear function of the angular velocity of rotation (Macdonald 1964; Peale 2005, see also Celletti and Lhotka 2014; Correia and Laskar 2004).

The unit of time is chosen to normalize the orbital period \(T_{orb}\) to \(2\pi \), which implies that the mean motion \(n=2\pi /T_{orb}\) is equal to one.

The equation of motion of the spin–orbit problem is introduced as follows. Let x be the rotation angle formed by the largest axis of the ellipsoid (belonging to the orbital plane, due to A2 and A3) and the periapsis line. Then, the spin–orbit problem under the conditions A1–A4 is described by the equation

$$\begin{aligned} \frac{d^2x}{dt^2}(t) + \varepsilon \left( \frac{a}{r(t)}\right) ^3 \sin \bigl (2 x(t) - 2f(t)\bigr ) = \mathcal {T} _d \biggl (\frac{dx}{dt}(t),t \biggr ), \end{aligned}$$
(1)

where the perturbative parameter \(\varepsilon \,{\mathop {=}\limits ^{\textrm{def}}}\,\frac{3}{2} \frac{\mathcal {B}-\mathcal {A}}{\mathcal {C}}\) measures the equatorial ellipticity of the satellite. The terms \(r(t)=r(t;e)\) and \(f(t)=f(t;e)\) are known functions depending on the eccentricity and related to time by means of the Kepler equation. In fact, denoting by \(t _0\) the initial time, Kepler equation \(nt+t _0=u-e\sin u\) gives the eccentric anomaly u as a function of the time t, while r and f are related to u by

$$\begin{aligned} r=a(1-e\cos u),\quad \tan \frac{f}{2} =\sqrt{\frac{1+e}{1-e}} \tan \frac{u}{2}. \end{aligned}$$
(2)

The right-hand side of (1) represents the dissipative effect, which is given by the tidal torque \(\mathcal {T} _d=\mathcal {T} _d(\frac{dx}{dt}(t),t)\); using the model developed in Macdonald (1964), (Peale 2005) (see Assumption A5), \(\mathcal {T} _d\) takes the form

$$\begin{aligned} \mathcal {T} _d \biggl (\frac{dx}{dt}(t), t\biggr ) \,{\mathop {=}\limits ^{\textrm{def}}}\,- \mu \left( \frac{a}{r(t)}\right) ^6 \biggl (\frac{dx}{dt}(t) - \frac{d f}{dt}(t) \biggr ), \end{aligned}$$
(3)

where \(\mu > 0\) is called the dissipative constant, which depends on the physical features of the satellite (density, rigidity, etc.).

We remark that the astronomical values of \(\varepsilon \) are small for many satellites and planets of the Solar system, typically of the order of \(10^{-4}\) for celestial bodies such as the Moon and Mercury, and that the dissipative term is typically small with respect to the conservative term, say \(O(\varepsilon ^2)\).

For \(\mu =0\), Eq. (1) is conservative, and for \(\varepsilon =\mu =0\) it is integrable. Thus, for small values of \(\varepsilon \) and \(\mu \), we obtain a nearly-integrable and nearly-Hamiltonian system.

Taking the average of the tidal torque over one orbital period (see, e.g., Peale 2005; Correia and Laskar 2004), Eq. (1) becomes

$$\begin{aligned} \frac{d^2x}{dt^2}(t) + \varepsilon \left( \frac{a}{r(t)}\right) ^3 \sin \bigl (2x(t)-2f(t) \bigr ) = -\mu \biggl (\bar{L}(e)\frac{dx}{dt}(t)-\bar{N}(e)\biggr ), \end{aligned}$$
(4)

where

$$\begin{aligned} \begin{aligned} \bar{L}(e)&= {\frac{1}{(1-e^2)^{9/2}}} \biggl (1+3e^2+\frac{3}{8}e^4 \biggr ), \\ \bar{N}(e)&= \frac{1}{(1-e^2)^6} \biggl (1+\frac{15}{2}e^2+\frac{45}{8}e^4+ \frac{5}{16}e^6 \biggr ). \end{aligned} \end{aligned}$$
(5)

The averaged model is obtained from the non-averaged model by eliminating the high-frequency terms in the tidal torque; notice that the terms in (5) are exact and they are not truncation of expansions in the eccentricity.

If we write the second-order Eq. (1) as a first-order system in phase space, we obtain:

$$\begin{aligned} \begin{aligned} \dot{x}&= y\\ \dot{y}&=- \varepsilon \left( \frac{a}{r(t)}\right) ^3 \sin \bigl (2 x-2f(t)\bigr )-\mu \left( \frac{a}{r(t)}\right) ^6 (y - \dot{f}(t)) \\ \dot{t}&=1 \end{aligned} \end{aligned}$$
(6)

and similarly for Eq. (4):

$$\begin{aligned} \begin{aligned} \dot{x}&= y \\ \dot{y}&=- \varepsilon \left( \frac{a}{r(t)}\right) ^3 \sin (2 x-2f(t))-\mu \biggl (\bar{L}(e)y-\bar{N}(e)\biggr ) \\ \dot{t}&=1. \end{aligned} \end{aligned}$$
(7)

Equations (1) and (4) are defined over the phase space \([0,2\pi )\times \mathbb {R}\), which can be endowed with the standard scalar product and the symplectic form \(\Omega = dy \wedge dx\). Both Eqs. (1) and (4) are dissipative in the sense that the phase space volume contracts under the evolution of the flow. Hence, if \(J_t\) is the determinant of the linearized flow, Abel’s formula gives

$$\begin{aligned} \det (J_t) = \exp {\int _0^t {{\,\textrm{Tr}\,}}(A(s))\ ds}, \end{aligned}$$

where \({{\,\textrm{Tr}\,}}\) denotes the trace and A(t) is the differential of the vector field at the time t flow.

With reference to Eq. (6), we obtain

$$\begin{aligned} {{\,\textrm{Tr}\,}}(A(t))=-\mu \left( \frac{a}{r(t)}\right) ^6 \end{aligned}$$

and consequently the volume contracts as

$$\begin{aligned} \det (J_t)=\exp \biggl (-\mu \int _0^t \left( \frac{a}{r(t)}\right) ^6\, ds \biggr ). \end{aligned}$$
(8)

For the averaged tidal torque model (7), we obtain

$$\begin{aligned} {{\,\textrm{Tr}\,}}(A(t))=-\mu \bar{L}(e) \end{aligned}$$

and consequently the volume contracts as

$$\begin{aligned} \det (J_t) = \exp \bigl (-\mu \bar{L}(e)\ t\bigr ). \end{aligned}$$
(9)

It is not difficult to see by an explicit computation (compare with Calleja et al. 2022), that at \(t = 2\pi \) both (8) and (9) give the same asymptotic contraction rate of the volume.

Note that when we do continuation in \(\varepsilon \), the value of the dissipation changes. This is somewhat different from previous studies in which the continuation was done for families of constant dissipation.

Remark 1

The relation between the averaged and non-averaged models has been considered puzzling in the literature. The difficulty is that averaging methods are usually justified for short time. On the other hand, attractors involve the evolution over very long times.

As we will see in Sect. 3, one of the advantages of the approach of Calleja et al. (2022) is that the computations of the attractors are based on time-1 maps (for which the standard justification of the averaging applies). The a posteriori format of the result in Calleja et al. (2013) shows that the attractors of the maps are close, when the maps are close. Hence, the a posteriori results justify that the approximations to finite time lead to conclusions in infinite time.

Besides the general argument above, Calleja et al. (2022) also presents a self-contained elementary argument showing that attractors of the two models are close over all times. In this paper, we present some numerical explorations that show that the rigorous results on the similarity of the attractors of the two models can be observed in numerical experiments, see Figs. 1 and 5.

3 The Poincaré Spin–Orbit Map

To perform the computations that will be presented in the next sections, it is convenient to reduce the study of (1) to the discrete system associated with the vector field (6) by a return map. According to the procedure detailed in Calleja et al. (2022), we consider the Poincaré map \(P _e\) associated with (6), which is obtained as follows. We introduce the map \(P_ e\) as

$$\begin{aligned} P _e(x_0,y_0;\varepsilon ) \,{\mathop {=}\limits ^{\textrm{def}}}\,\begin{pmatrix} x(2\pi ;x_0,y_0,\varepsilon ) \\ y(2\pi ;x_0,y_0,\varepsilon ) \end{pmatrix}, \end{aligned}$$
(10)

where \(x(2\pi ;x_0,y_0,\varepsilon ) \) and \(y(2\pi ;x_0,y_0,\varepsilon )\) represent the solution of (6) at time \(t=2\pi \) with initial conditions \((x_0,y_0)\) at \(t=0\). Writing \(P _e\) in components as \(P _e\equiv (P _e^{(1)},P _e^{(2)})\), the spin–orbit Poincaré map is given by

$$\begin{aligned} \begin{aligned} \bar{x}&= P _e^{(1)}(x,y;\varepsilon ),\\ \bar{y}&= P _e^{(2)}(x,y;\varepsilon ). \end{aligned} \end{aligned}$$
(11)

For computational reasons, it is convenient to consider the change of coordinates

$$\begin{aligned} \Psi _e\,{\mathop {=}\limits ^{\textrm{def}}}\,2\pi \begin{pmatrix} 1 &{} 0 \\ 0 &{} 1 - e\end{pmatrix} \end{aligned}$$
(12)

and define the map

$$\begin{aligned} G_e\,{\mathop {=}\limits ^{\textrm{def}}}\,\Psi _e\circ P _e\circ \Psi _e^{-1}, \end{aligned}$$
(13)

which is in fact the Poincaré spin–orbit map for (6) with the change of coordinates given by Kepler’s equation, \(t = u - e\sin u\).

The map (11), equivalently (13), inherits some properties from the differential Eq. (1) (equivalently, (6)). Notably, (11) (or (13)) is a nearly-integrable dissipative map, which is indeed conformally symplectic according to the following definition (see Calleja et al. 2013).

Definition 1

Let \(\mathcal {M} = \mathbb {T}^n\times U\) with \(U\subseteq \mathbb {R}^n\) an open, simply connected domain with smooth boundary and with symplectic form \(\Omega \). A family of diffeomorphisms \(f _\vartheta :\mathcal M\rightarrow \mathcal M\), depending on a parameter \(\vartheta \), is conformally symplectic, if there exists a function \(\lambda :\mathcal {M}\rightarrow \mathbb {R}\) such that

$$\begin{aligned} f _\vartheta ^* \Omega = \lambda \Omega , \end{aligned}$$
(14)

where \(f _\vartheta ^*\) denotes the pullback.

The quantity \(\lambda \) is called the conformal factor. For \(\lambda =1\), we recover the symplectic case. For \(n=1\), any diffeomorphism is conformally symplectic with \(\lambda (x)=\pm |\det (D f _\vartheta (x))|\). For \(n\ge 2\), it is known that \(\lambda \) is constant (see Calleja et al. 2013 for details).

For later use, we denote by J the symplectic matrix representing \(\Omega \) at \({\underline{z}}\), namely \(\Omega _{\underline{z}}({\underline{u}},{\underline{v}})=({\underline{u}},J({\underline{z}}){\underline{v}})\) for \({\underline{u}},{\underline{v}}\in \mathbb {R}^n\) and \((\cdot ,\cdot )\) is the standard inner product in \(\mathbb {R}^n\). For the spin–orbit map, the matrix J is given by

$$\begin{aligned} J = \begin{pmatrix} 0 &{} 1 \\ -1 &{} 0 \end{pmatrix}. \end{aligned}$$
(15)

The Poincaré spin–orbit map (11) (or (13)) is conformally symplectic (Calleja et al. 2022) with the conformal factor given by

$$\begin{aligned} \lambda = \exp \biggr (-\mu \pi \frac{3 e^4+24 e^2+8}{4 \left( 1-e^2\right) ^{9/2}}\biggl ). \end{aligned}$$
(16)

The system contracts the volume for \(\mu >0\), expands the volume for \(\mu <0\), and it is neutral for \(\mu =0\). In the following sections, we will only take \(\mu > 0\).

Note that the conformal factor of the averaged model (4) will have the same \(\lambda \) given in (16), since at time \(t=2\pi \) both contraction rates (8) and (9) coincide precisely due to the \(\bar{L}(e)\) choice. For the calculation leading to (5), we refer to Calleja et al. (2022, p. 15).

4 Basins of Rotation Numbers

In this section, we present exploratory results based on the computation of the rotation number of many orbits of the spin–orbit problem. These computations provide a procedure to obtain an approximate value of the drift parameter associated with a given frequency of an invariant torus. Such information is essential to start the procedure detailed in Calleja et al. (2022) for the computation of a KAM torus with fixed frequency.

The phase space associated with the spin–orbit map (11) is a subset of the cylinder \([0,2\pi ) \times \mathbb {R}\). For numerical reasons, it is convenient to take the eccentric anomaly u as independent variable, instead of time. This can be achieved by taking into account Kepler’s equation, which provides a relation between t and u: \( t = u - e\sin u\). More precisely, let \(s(x;u,e) \,{\mathop {=}\limits ^{\textrm{def}}}\,\sin (2x(t)-2f(t))\), which depends on u and \(e\) through f, that satisfies the following identities:

$$\begin{aligned} \cos f = \frac{\cos u - e}{1 - e\cos u} \quad \text {and} \quad \sin f = \frac{\sqrt{1 - e^2}\sin u }{1 - e\cos u}. \end{aligned}$$
(17)

We use simple trigonometric identities to write the function s as

$$\begin{aligned} s(x; u, e) = \sin (2x) (2\cos ^2 f-1) - \cos (2x)2\cos f \sin f \end{aligned}$$
(18)

and we define the function c as

$$\begin{aligned} c(x; u, e) \,{\mathop {=}\limits ^{\textrm{def}}}\,\cos (2x) (2\cos ^2 f-1) + \sin (2x)2\cos f \sin f. \end{aligned}$$
(19)

Note that both quantities are easier to compute than the one in terms of \(\tan (f/2)\) in (2) which has singularities depending on the value of u. Moreover, (18) and (19) satisfy the following relations, which will be used during the integration with Taylor’s numerical integration method:

$$\begin{aligned} \frac{\partial s}{\partial x}(x; u, e) = 2c(x; u, e), \quad \frac{\partial c}{\partial x}(x; u, e) = -2 s (x; u, e). \end{aligned}$$

As a consequence of the change of coordinates, we deduce that

$$\begin{aligned} \frac{d f}{dt} = \left( \frac{a}{r(u)}\right) ^2 \sqrt{1 - e^2}. \end{aligned}$$

The spin–orbit problem in (6) can then be expressed in terms of the independent variable u as

$$\begin{aligned}{} & {} \frac{d^2 \beta }{du ^2}(u) - \frac{d\beta }{du}(u) \frac{a}{r(u)} e\sin u + \varepsilon \frac{a}{r(u)} s(\beta ;u, e) \nonumber \\{} & {} \quad = - \mu \left( \frac{a}{r(u)}\right) ^5 \biggl (\frac{d\beta }{du}(u) - \frac{a}{r(u)} \sqrt{1 - e^2}\biggr ), \end{aligned}$$
(20)

where \(\beta \) and \(\gamma \) are defined as

$$\begin{aligned} \beta (u) \,{\mathop {=}\limits ^{\textrm{def}}}\,x (u - e\sin u), \quad \gamma (u) \,{\mathop {=}\limits ^{\textrm{def}}}\,\frac{d\beta }{d u}(u) = \frac{r(u)}{a}y(u - e\sin u). \end{aligned}$$

The time-(\(2\pi \)) map associated with (20) is related, up to a factor, to the map defined in (13).

4.1 Computation of the Rotation Numbers

In this section, we present the computation of the rotation numbers obtained by a direct iteration; in Sect. 4.2, we will map the regions covered by different attractors. First, we recall that the definition of the rotation number associated with Eq. (6) (equivalently (7)) is:

$$\begin{aligned} \rho = \lim _{n \rightarrow \infty } \frac{1}{n} \sum _{j = 1}^{n-1} y(2\pi j). \end{aligned}$$
(21)

The rotation number of orbits of two-dimensional maps could fail to exist, but when the orbits lie in an invariant circle (the main object of our interest in this paper), such rotation number exists. Also, all the orbits in the basin of attraction of a circle have the same rotation number.

The convergence of the limit in (21) may be very slow, say of the order O(1/n). As shown in Das et al. (2017), a method alternative to (21) to compute the rotation number is through the following formula:

$$\begin{aligned} \rho = \lim _{n \rightarrow \infty } \biggl [\frac{1}{\sum _{j = 1}^{n-1} \phi (\tfrac{j}{n})} \sum _{j = 1}^{n-1} y (2\pi j) \phi (\tfrac{j}{n})\biggr ], \end{aligned}$$
(22)

where \(\phi \) denotes a weight function, e.g.,

$$\begin{aligned} \phi (z) \,{\mathop {=}\limits ^{\textrm{def}}}\,\exp \biggl (\frac{-1}{z^2(1-z)^2} \biggr ), \quad z \in (0,1). \end{aligned}$$
(23)

It is shown in Das et al. (2017) that if the rotation number is sufficiently irrational (and the corresponding circle is smooth), the limit in (22) is reached very fast. For example, if the circle is analytic and its frequency is Diophantine, the limit is reached super-exponentially fast.

In contrast, if the limiting rotation number is rational (or if the circle is not smooth), the convergence of (22) is slow. Hence, the speed of convergence on (22) can be used as a diagnostic whether the torus is present and smooth. For our purposes, the main use of (22) is to adjust parameters, so that an attractor with the right rotation number is found.

The computation of the rotation number for the system (6), associated with an attractor close to the initial condition \((x _0, y _0) \in [0,2\pi ) \times \mathbb {R}\), can be performed according to Algorithm A.Footnote 1

Algorithm A

Computation of the rotation number for the spin–orbit problem with tidal torque, see Eq. (6).

  • Inputs: Initial condition: \((x _0, y _0) \in [0,2\pi )\times \mathbb {R}\). Spin–orbit parameters: \(\varepsilon > 0\), \(e\in [0,1)\), \(\mu > 0\), and \(\lambda \) as in (16). Positive integers: \(\delta \), \(n_1\), \(n_2\) with \(n _1 < n _2\), and \(n _0 \leftarrow \lceil -14/\log _{10} \lambda \rceil \). Weight function: \(\phi \), for instance, the one given in (23).

  • Outputs: Approximate rotation number \(\rho \) for (6) (equivalently, for (7)).

    1. 1.

      \((\beta _0, \gamma _0) \leftarrow (x _0,y _0 / (1-e))\).

    2. 2.

      \((\beta _0, \gamma _0) \leftarrow G_e^{n _0}(\beta _0, \gamma _0)\) by integrating (20) for \(2\pi n _0\) times in \([0,\pi )\times \mathbb {R}\).

    3. 3.

      Store \(\gamma _j \leftarrow \gamma (2\pi j; \beta _0, \gamma _0)\) for \(j = 1, \ldots , n _2\) by integrating (20).

    4. 4.

      \(n \leftarrow n _1\).

    5. 5.

      \( \rho \leftarrow \bigl (\sum _{j = 1}^{n-1} \phi (\tfrac{j}{n})\bigr )^{-1}\, \sum _{j = 1}^{n-1} \gamma _j \phi (\tfrac{j}{n})\).

    6. 6.

      Iterate steps (4)–(5) with \(n \leftarrow n + \delta \) until \(\rho \) does not vary or \(n > n _2\).

    7. 7.

      If convergence, return \(\rho / (1 - e)\).

To find a reliable rotation number of a randomly chosen orbit, it is convenient to discard a transient number of iterations to ensure that the orbit has converged to the 1-D attractor; the parameter \(n_0\) in step 2 of Algorithm A implements this. Unfortunately, the time that an orbit needs to settle in an attractor is hard to estimate and it may be surprisingly large (see Celletti and Lhotka 2014), because an orbit may be entertained near an attractor for a long time before landing in the final one. This effect is very hard to predict and changes wildly with the initial conditions.

Since the only harm in overestimating the transition is the use of more computer time, we have tried to overestimate it. Note also that iterating different points is very parallelizable.

We report in Fig. 1 the values of the eccentricities that lead to invariant attractors with given rotation number. As sample cases, we are interested to the following two irrational numbers:

$$\begin{aligned} \omega _1&\,{\mathop {=}\limits ^{\textrm{def}}}\,\gamma _g^+ \approx \mathtt {1.618033988749894848}\ldots , \end{aligned}$$
(24)

and

$$\begin{aligned} \omega _2&\,{\mathop {=}\limits ^{\textrm{def}}}\,1 + \frac{1}{2+\gamma _g^-} \approx \mathtt {1.3819660112501051} \ldots \end{aligned}$$
(25)

with \(\gamma _g^\pm \) equal to the golden ratio and its conjugate: \(\gamma _g^\pm \,{\mathop {=}\limits ^{\textrm{def}}}\,\frac{\sqrt{5}\pm 1}{2}\). Note that both \(\omega _1, \omega _2\) are noble numbers (i.e., the continued fraction expansion has only 1 after a certain point). It is known that one of the predictions of the standard renormalization group, see MacKay (1983) is that for maps in the universality domain many of the fine properties at breakdown are similar for all tori with noble rotation numbers.

Figure 1 shows \(\omega _1\) and \(\omega _2\) with two different dashed lines. Moreover, it shows with a dotted-dashed line the term \(\bar{N}(e)/\bar{L}(e)\) (see (5)), which is the predicted rotation number in the model with the averaged tidal torque, see Celletti and Chierchia (2009), Calleja et al. (2022). Extremely close to the dotted-dashed line, the crosses are obtained by computing the rotation number for the non-averaged model (6). To highlight the proximity between the results associated with the averaged and the non-averaged models, Fig. 1 shows also a zoom-in nearby the frequency \(\omega _2\). This shows that the rigorous results in Calleja et al. (2022) justifying the averaging method give very good numerical values.

Fig. 1
figure 1

Rotation numbers (+) associated with (6) from Algorithm A in terms of the eccentricity \(e\) with weight function as in (23), \(n _1 = 4500\), \(n _2 = 4600\), \(\delta = 10\), \(\varepsilon =10^{-4}\), \(\mu =10^{-5}\) and \((x _0, y _0) = (0, 1.38)\). The fixed frequencies \(\omega _1\) and \(\omega _2\) are given in (24) and (25), and the predicted rotation number for the averaged spin–orbit is given by \(\bar{N}(e)/\bar{L}(e)\). Finally, a zoom-in is shown nearby \(\omega _2\) with the + connected by a straight line

Full size image

4.2 Computation of the Basins of Rotation Numbers

A qualitative picture of the dynamics can conveniently be achieved by computing the basins of rotation numbers. These basins are obtained by applying recurrently Algorithm A; given a grid in \((x _0, y _0)\) of size \(n _x \times n _y\) (for some positive integers \(n_x\), \(n_y\)) in a window of the cylinder \([0,2\pi ) \times \mathbb {R}\), we apply the Algorithm A and represent the rotation number as a color map, using interpolation for the elements that are not in the finite grid.

Notice that the iterations of each orbit are completely independent from each other, which makes the procedure completely parallelizable. Each orbit could be assigned to a thread (using, e.g., openMP Dagum and Menon 1998) or, working on different regions, could be distributed using coarse grained distributors such as GNU-parallel or HTCondor (see Tange 2011; Htcondor 1998). We have used these techniques systematically.

Although we have not yet used GPUs, the computation of these basins is indeed suitable for GPUs. For preliminary calculations, the single precision available even in consumer grade GPUs is perfectly acceptable. However, when dealing with the more detailed calculations of tori, even double precision is not enough and we need extended precision.

Figure 2 shows some basins of rotation numbers for different values of the parameters. On average, each of the plots, using the taylor program Jorba and Zou (2005), needs around 12 min using 143 CPUs. The color scale gives the rotation number for every initial condition. We provide the results for two different values of the eccentricity. The plots clarify regions of librational and rotational motion. At a finer inspection, the regions close to the location of the invariant tori with frequencies \(\omega _1\) in the upper plots of Fig. 2 and \(\omega _2\) in the lower plots become more irregular as the perturbing parameter \(\varepsilon \) ranges over values far below from what will be the critical threshold (left column), close to the critical threshold (middle column) or far above the critical threshold (right column). This behavior is consistent with the theory on the existence of KAM tori developed in Calleja et al. (2013).

Fig. 2
figure 2

Basins of rotation numbers with grid size \(500\times 500\) and window \([0,2\pi ) \times [1,2]\) for the spin–orbit problem (6) with \(\mu =10^{-3}\), different values of \(\varepsilon \) per each column, and \(e= 0.316\) in the upper row, \(e=0.247\) in the lower row. The color scale provides the value of the rotation number for each initial condition (x, y), using interpolation for points outside the grid

Full size image

5 Explorations of the Breakdown

A rotational torus defined for a conformally symplectic system (see Definition 1) is introduced as follows.

Definition 2

Let \(f _\vartheta :\mathcal {M} \rightarrow \mathcal {M}\) be a conformally symplectic system, and let \(\omega \) be an irrational number. A map \(K :\mathbb {T}^n \rightarrow \mathcal {M}\) is said to be a rotational torus for \(f _\vartheta \) with frequency \(\omega \), when there exists \(\vartheta _0\) such that

$$\begin{aligned} f _{\vartheta _0} \circ K(\theta ) = K (\theta + \omega ) \quad \text {for all } \theta \in \mathbb {T}^n. \end{aligned}$$

The algorithm in Calleja et al. (2013) deals with tori whose frequency has Diophantine properties and it computes iteratively corrections for the map K and the parameter \(\vartheta _0\) such that Definition 2 is satisfied up to a numerical tolerance. If there is an extra parameter in the problem, it is standard to use the Newton method to implement a continuation method. The solutions of the invariance equation for a parameter value are taken as an initial guess for the problem with a slightly larger parameter value and applying the Newton method, we obtain a solution for the increased parameter value.

In many situations, this continuation method stops because the algorithm becomes unreliable (e.g., the expansions used in the truncation of the problem do not reach the function) or because the torus ceases to exist (as a smooth object).

In the following sections, we present very accurate computations of these tori and we discuss some interesting phenomena near the spin–orbit breakdown.

5.1 Computation of the Tori and Their Continuation Toward the Breakdown

The theorems and algorithms in Calleja et al. (2013) provide a way to efficiently compute KAM tori in conformally symplectic systems. Those results have been adapted for the spin–orbit problem in Calleja et al. (2022). The main difference between the symplectic and the conformally symplectic cases is that, in the dissipative case at each time that the torus is corrected for a fixed rigid rotation, a parameter of the system must be corrected as well. In fact, while the torus with a fixed frequency can exist for each parameter value in a conservative system, in a dissipative one the torus may only exist for certain parameter values. Moreover, in a conformally symplectic system a torus with a given rotational frequency exists for parameter values belonging to a whole interval, the so-called Arnold tongue (Arnol’d 1965; Calleja et al. 2014).

The method detailed in Calleja et al. (2013) is really efficient and accurate in computing the torus and the drift parameter of the system given approximate values. This makes it extremely suitable for a continuation algorithm. Here we have taken the solution for a parameter value as an approximate solution for a slightly bigger value of the parameterFootnote 2.

Furthermore, the continuability argument in Calleja et al. (2013) shows that the continuation method will get as close as desired to the boundary of existence if given enough computer resources. In practice, getting extremely close to the boundary of existence may require using a very large number of Fourier modes, extended precision arithmetic, etc. Note that several KAM proofs and their attendant algorithms do not satisfy this property. This requires that the parametrization used reaches the maximal domain, see Llave and Rana (1991).

In reality, of course, one has limited computer resources (limited memory, computer time and time of the programmer). In what follows, we will explore the limits obtained using standard computers (desktops or small servers) at the time of the writing. We note that the algorithms based on automatic reducibility specified in Calleja et al. (2013) are very efficient both in memory required, in number of operations and in the dimension of the objects studied. The algorithm is also based on very structured operations which make it easy to change the precision of the arithmetic and the number of Fourier modes in a modern programming language.

One of the results of our explorations is that the main bottleneck for the accuracy of the calculation near the breakdown is the precision of the arithmetic. Therefore, in our studies, we use high precision arithmetic and a very high-order ODE integrator. Both of these improvements are doable in a reasonable programming time, thanks to the taylor program Jorba and Zou (2005), which generates with minimal programming effort a Taylor solver (very high order) with extended precision arithmetic. This is crucial to obtain high precision for the Poincaré spin–orbit map.

5.2 Some Implementation Details; Running Times

For the frequency (24), the computation was done in a machine with Intel Xeon Gold 5220 CPU at 2.20 GHz with 18 CPUs with hyperthreading which simulates 36 CPUs, while for (25) a machine with Intel Xeon Gold 6154 CPU at 3.00 GHz with 74 CPUs with hyperthreading which simulates 148 CPUs. In case of parallelism, we have always requested 32 CPUs.

5.2.1 Choice of Continuation Parameters

The procedure in Calleja et al. (2022) involves computing the derivatives of the spin–orbit Poincaré map (11) with respect to (x, y) and with respect to the parameter to be corrected. In our case, we will use the eccentricity as the parameter to be adapted. The eccentricity has a clear physical interpretation. We will perform our calculations for tori of the two frequencies (24) and (25). Thus, we get \((K _{\omega _1}, e_{\omega _1})\) and \((K _{\omega _2}, e_{\omega _2})\) such that

$$\begin{aligned} P _{e_{\omega _j}} \circ K _{\omega _j} (\theta ) = K _{\omega _j} (\theta + \omega _j), \quad j \in \{1,2\}. \end{aligned}$$
(26)

We fix the dissipation parameter to \(\mu =10^{-3}\), and we perform a continuation w.r.t. the parameter \(\varepsilon \). The continuation is based on using a cubic extrapolation from the previous results to provide an initial guess for the next continuation step and then, polishing this guess by running the Newton method as many times as needed to achieve the desired accuracy.

5.2.2 Step Control and Adaptability of the Algorithm

Besides the usual step control in continuation methods, we incorporate some control flags to ensure that the computed embedded torus satisfies the invariance equation with enough accuracy and it is a reasonable function.

One of the flags monitors the error in the invariance Eq. (26). To declare the computation successful, we require that the error is below a certain value.

We note that in this problem, specially near the breakdown, it is possible to obtain spurious solutions that indeed solve very accurately the invariance equation, but which are very unreasonable functions. Hence, we introduce a second flag which checks that the tail of Fourier coefficients has norms between two tolerances \(\epsilon _1 \le \epsilon _2\). If the size of the tail is bigger than \(\epsilon _2\), we increase the number of coefficients used; if the tail is smaller than \(\epsilon _1\), we decrease the number of Fourier coefficients (to increase the speed). Hence, we declare a computation successful when it achieves small residual and when the solution has a small enough tail, so that it can be computed with different levels of truncation.

In case of failure, in any of the control flags, remesh and increase the number of the Fourier coefficients in the last successful computation, run the Newton method to get more accuracy (and check that the new solution has good flags) and repeat the process. Of course, there is a global limit to the number of Fourier modes and the calculation finishes when we go over the limit. In this run, we have used \(L _{\max } = 65536\). Figure 3 shows the different mesh values used during the continuation. We also use standard adaptative step control in continuation, which stops the calculation when a step below the minimum does not achieve that the Newton method converges. If the Newton method succeeds, the initial guess for the next step will be obtained by doing a cubic extrapolation of the previously computed tori with the same mesh size.

5.2.3 Run Times

We have made several runs with different values for these tolerances and checked that the results do not change appreciably by changing the choices of parameters of the algorithm.

Fig. 3
figure 3

Mesh size and execution time (in log10 min) of successful continuation steps for the non-averaged spin–orbit problem (6) and for the frequencies (24) and (25)

Full size image

For the spin–orbit model, we use 70 digits of precision, \(10^{-35}\) of Newton tolerance for the error of the invariance Eq. (26), tail tolerance between \(\epsilon _1=10^{-55}\) and \(\epsilon _2=10^{-28}\), and \(10^{-70}\) as tolerance for the absolute and relative errors in the Taylor integration.

The last values for the perturbing parameter \(\varepsilon \) that we reached using the \(\varepsilon \) values in the continuation strategy described above are (27) for \(\omega _1\) and (28) for \(\omega _2\):

$$\begin{aligned} \varepsilon _{\omega _1} ^{c}&= \texttt {1.265364670507455833901687945901589e-02}, \end{aligned}$$
(27)
$$\begin{aligned} \varepsilon _{\omega _2} ^{c}&= \texttt {1.372784208166277584502189850802202e-02}. \end{aligned}$$
(28)

The corresponding tori are plotted in Fig. 4, and the values of the drift parameter for \(\omega _1\) and \(\omega _2\) are:

$$\begin{aligned} e_{\omega _1} ^{c}&= \texttt {3.1701530650181080344148405746028386e-01}, \end{aligned}$$
(29)
$$\begin{aligned} e_{\omega _2} ^{c}&= \texttt {2.4797274489383016717182991353216456e-01}. \end{aligned}$$
(30)
Fig. 4
figure 4

Tori at the numerical breakdown for the non-averaged (6) spin–orbit problems and for the frequencies (24) (left panel) and (25), respectively

Full size image

The critical perturbative and drift parameters, (27)–(30), do not substantially differ for the averaged and non-averaged models (respectively, (4) and (6)). For example, Fig. 5 shows that the difference in the drift parameter (the eccentricity) of the KAM tori, obtained using the continuation procedure, for the two models is of the order of \(10^{-7}\). Therefore, in single-precision arithmetic, these averaged and non-averaged drifts may look indistinguishable.

The mesh size is the same in the two models and for the studied frequencies, i.e., Fig. 3 is the same for the averaged model. Also, the computational times to compute the tori are similar. In general, the averaged version is faster far from the breakdown. Figure 6 illustrates the time ratios (using in all the cases a maximum of 32 CPUs), while Fig. 3 provides the time of each successful continuation step in \(\varepsilon \). We strongly encourage to not infer general statements from the computational times, since they may differ due to cache sizes, CPU machines, running of other jobs, changes in the continuation strategies, initial conditions, number of Newton iterations, accurate control flags, etc.

The main computational bottleneck is on the Newton steps. For each discretized value of \(\theta \) parametrizing the invariant torus, a numerical integration of the ODE and its variational equations must be performed. Also, different continuation strategies in the stepsize choice and initial guesses can improve the overall time.

Motivated by these remarks, in what follows, we are going to provide only the plots for the non-averaged spin–orbit problem. The plots of the two models are visually identical, since the difference between them in the regimes we study is about 5 orders of magnitude smaller than the main effects.

Fig. 5
figure 5

Difference of the eccentricities for the averaged and non-averaged spin–orbit problems (4) and (6), respectively, and for the frequencies (24) (left panel) and (25) (right panel)

Full size image
Fig. 6
figure 6

Time ratios of each successful continuation step for the averaged and non-averaged spin–orbit problems (4) and (6), respectively, and for the frequencies (24) (left panel) and (25) (right panel). It run with 32 CPUs

Full size image

6 Studies of the Breakdown of Tori

Our results on the behavior at breakdown of the spin–orbit map are described in Sects. 7 and 9, where we have taken the algorithms to their limits of validity. Our study is motivated by different reasons:

  1. (i)

    the spin–orbit problem provides an excellent testing ground for the algorithms, which is far from the academic standard models;

  2. (ii)

    the a posteriori method allows to compute with confidence very close to the breakdown;

  3. (iii)

    the phenomena at the breakdown of validity of KAM theory have a great mathematical interest, also because they are related to renormalization group theory (Östlund et al. 1983; Rand et al. 1982). For us, the phenomena at breakdown is a side issue, since the main goal of this paper is to develop a methodology to maintain accuracy and reliability even up to breakdown in a model of astronomical interest. Clearly, the phenomena happening at the breakdown are of certain interest and deserve another study.

We anticipate that the results that we obtain are that the breakdown of the dissipative spin–orbit problem does not conform to the behaviors that have been previously found in the literature, MacKay (1983), Rand et al. (1982).

6.1 Some Rigorous Results on the Breakdown

The rigorous result in Calleja et al. (2013) has a remarkable consequence concerning the behavior of the torus, and of the stable and tangent bundles at breakdown. Precisely, if:

  1. (i)

    the invariant attractor is smoothly conjugate to a rotation,

  2. (ii)

    the invariant torus has a smooth invariant direction which is:

    • contractive,

    • with an angle bounded below from the tangent direction,

then the torus can be continued as a smooth curve and the iterative method specified in Calleja et al. (2013) will converge for a small enough perturbation.

As a consequence of the above results, at breakdown at least one of the above requirements (i) or (ii) has to fail. In other words, either the hyperbolicity is lost or there is a loss of regularity.

Hence, our numerical explorations focus on the behavior of the bundles (Sect. 7), the regularity of the tori (Sect. 8) and the regularity in Sect. 9. Since previous studies found scaling relations, we have also studied the possibility of scaling invariant ratios (Sect. 10).

There are some rigorous results that further limit the phenomena that can happen.

As for the breakdown by loss of hyperbolicity, in the spin–orbit problem the determinant is the friction, since the dynamics on the circle is a rotation, the stable exponent is given by the friction, so the only way that hyperbolicity can be lost is if the stable and tangent bundles become close at some points. This phenomenon was found numerically in many examples (Haro and de la Llave 2006b, a) and called bundle collapse. In Haro and de la Llave (2006b), it was also found that this phenomenon happens often and that it satisfies remarkable scaling relations. The paper Haro and de la Llave (2006b) contains a proof that this phenomenon indeed happens in some examples. In Bjerklöv and Saprykina (2008), Figueras and Timoudas (2020) there is a proof of the scaling relations of bundle collapse in some models.

As for the breakdown by loss of regularity, we note:

  • If the torus and the bundles are \(C^1\), using that the normal exponent has to be the fraction, the results in Fenichel (1973-1974) show that the torus and the bundle are \(C^r\) for any \(r \in \mathbb {N}\).

  • Using Herman (1979), Yoccoz (1984), Khanin and Sinai (1986), we obtain that the dynamics is \(C^{r-a}\) conjugate to a rotation for some \(a \in \mathbb {R}_+\) related to the Diophantine exponent.

  • If the torus is \(C ^{r-a}\) conjugated to a rotation, then the bundles are \(C ^{r-2a}\).

  • In Calleja et al. (2013), it is shown that if the torus and the bundles are \(C^r\) for r sufficiently high, then the torus and the bundles are analytic.

Therefore, at the breakdown by loss of regularity, then the regularity of the torus or the bundle has to go from analytic to less than \(C^1\).

We also note that, as we will see, the bundles can be obtained from the conjugating function by solving cohomology equations, so that the only way that bundles can lose regularity is if the regularity of the torus breaks down.

This drastic drop of regularity justifies that we talk about the breakdown. In other situations, the regularity seems to decrease gradually (Yao and Llave 2021).

We also note that, even after the breakdown, there could be other invariant objects of lower regularity. For example, Capiński and Kubica (2020) presents a topological extension of normal hyperbolicity beyond \(C^1\) regularity. The objects therein are continua, not manifolds. For dissipative systems, there are topological (Le Calvez 1986, 1988) or variational arguments (Marò and Sorrentino 2017), showing existence of quasi-periodic orbits of a rotation number, even if they are not dense on a circle.

The closest precedents to our study are Haro and de la Llave (2006b), Haro and de la Llave (2006a) and, specially Calleja and Figueras (2012), which studies conformally symplectic systems. In these papers, it was found at the same time a bundle collapse and a loss of regularity as well as remarkable scaling relations.

7 Behavior of Invariant Bundles

In this section, we study the behavior at breakdown of the stable and tangent bundles. To this end, let \((K,e)\) be an invariant torus with frequency \(\omega \) for the map \(P _e\) given in (11). Let \(N(\theta )\), \(M (\theta )\), \(S(\theta )\) be the quantities defined by

$$\begin{aligned} \begin{aligned} N(\theta )&\,{\mathop {=}\limits ^{\textrm{def}}}\,(DK(\theta )^\top DK(\theta ))^{-1}, \\ M(\theta )&\,{\mathop {=}\limits ^{\textrm{def}}}\,[DK(\theta )\ |\ J^{-1}\circ K(\theta )\ DK(\theta ) N(\theta )], \\ S(\theta )&\,{\mathop {=}\limits ^{\textrm{def}}}\,((DK N)\circ T_\omega )^\top (\theta ) DP_{e} \circ K(\theta ) J^{-1}\circ K(\theta )DK(\theta )N(\theta ), \end{aligned} \end{aligned}$$

where \(T_\omega \) denotes the shift by \(\omega \): \(T_\omega (\theta )=\theta +\omega \). Then, one can show (Calleja et al. 2013) that N, M, S and \(P_e\) satisfy the relation

$$\begin{aligned} D P _e\circ K (\theta ) M (\theta ) = M (\theta + \omega ) \begin{pmatrix} 1 &{} S (\theta ) \\ 0 &{} \lambda \end{pmatrix} \end{aligned}$$

with \(\lambda \) defined in (16). Similarly, the stable and tangent bundles, say \(E ^s(\theta )\) and \(E^c(\theta )\), respectively, must satisfy the reducibility equation (see Calleja and Figueras 2012)

$$\begin{aligned} D P _e\circ K(\theta ) W(\theta ) = W(\theta + \omega ) \begin{pmatrix} 1 &{} 0 \\ 0 &{} \lambda \end{pmatrix} \end{aligned}$$

with \(W (\theta ) = \begin{pmatrix} E^c(\theta )&E ^s(\theta ) \end{pmatrix}\).

To reduce the cocycle, we introduce the change of variables

$$\begin{aligned} \widehat{W} (\theta ) = \begin{pmatrix} 1 &{} B (\theta ) \\ 0 &{} 1 \end{pmatrix} \end{aligned}$$

satisfying

$$\begin{aligned} \begin{pmatrix} 1 &{} S (\theta ) \\ 0 &{} \lambda \end{pmatrix} \widehat{W} (\theta ) = \widehat{W} (\theta + \omega ) \begin{pmatrix} 1 &{} 0 \\ 0 &{} \lambda \end{pmatrix}. \end{aligned}$$

Hence, the unknown function \(B(\theta )\) verifies the following cohomology equation:

$$\begin{aligned} B(\theta ) - \lambda B (\theta + \omega ) = - S (\theta ), \end{aligned}$$
(31)

which can be solved by expanding in Fourier series whenever \(| \lambda | \ne 1\). Finally, the tangent and stable bundles can be recovered from the relation \(W(\theta ) = M (\theta ) \widehat{W} (\theta )\). Explicitly, if K has components \((K _1, K _2)\), then

$$\begin{aligned} W(\theta ) = \begin{pmatrix} D K _1(\theta ) &{} D K _1(\theta ) B (\theta ) - D K _2(\theta ) N (\theta ) \\ D K _2(\theta ) &{} D K _2(\theta ) B (\theta ) + D K _1(\theta ) N (\theta ) \end{pmatrix}. \end{aligned}$$

A straightforward computation of the inner products between the bundles leads to

$$\begin{aligned} \begin{aligned} E^c(\theta )^\top E^s(\theta )&= N(\theta ) ^{-1} B(\theta ), \\ E^c(\theta )^\top E^c(\theta )&= N(\theta ) ^{-1}, \\ E^s(\theta )^\top E^s(\theta )&= N(\theta ) ^{-1} B(\theta )^2 + N(\theta ). \end{aligned} \end{aligned}$$

Therefore, if \(\alpha (\theta )\) is the angle between the invariant stable and tangent bundles at the point \(\theta \in \mathbb {T}\), then

$$\begin{aligned} \cos \alpha (\theta ) = \frac{B(\theta )}{\sqrt{N(\theta )^2 + B(\theta )^2}} \end{aligned}$$

or equivalently

$$\begin{aligned} \tan \alpha (\theta ) = \frac{N(\theta )}{B(\theta )}. \end{aligned}$$
(32)

Figure 7 shows the behavior of the angle \(\alpha \) for the averaged and non-averaged models (7) and (6). We notice that the minimum value of the angle of separation does not go to zero as the parameter \(\varepsilon \) approaches the breakdown threshold.

Fig. 7
figure 7

Minimum angle between the stable and tangent bundles of the KAM tori of (7) (left panel) and (6) (right panel) with \(\mu =10^{-3}\), \(e\) as drift, and \(\varepsilon \) as continuation parameter

Full size image

In Calleja and Figueras (2012), it is shown that the angle between the tangent and stable bundles of the dissipative standard map collided in a very particular manner, while we have realized that in the spin–orbit problem the approach to breakdown occurs in a different way. More precisely, let us consider the angles w.r.t. the semi-axis \(\{x>0\}\), which are obtained by computing the arc tangent between the second coordinate and the first one in each of the bundles. Figure 8 shows these angles for the center and stable bundles as well as their difference for the \(\varepsilon \) values given in (27) and (28). We observe that the difference between these angles does not coincide as in the dissipative standard map, yielding a non-collapse bundle near the breakdown in the sense of Calleja and Figueras (2012). Moreover, the minimum value happens in the same value around \(\theta \approx 0.5\) which is coherent with the results in Bjerklöv and Saprykina (2008).

Fig. 8
figure 8

Angles between the center and stable bundles w.r.t. the semi-axis \(\{x>0\}\) and its difference for \(e\) ((29) upper left, (30) lower left) and \(\varepsilon \) (27) upper right, (28) lower right) values corresponding to each of the frequencies (24) and (25) of the KAM tori of the non-averaged spin–orbit (6) with fixed \(\mu =10^{-3}\)

Full size image
Fig. 9
figure 9

Quotient between the difference in the angle of the center and stable bundles for frequencies (24) and (25), respectively, see Fig. 8. As in that figure, the computations refer to KAM tori of the non-averaged spin–orbit (6) with fixed \(\mu =10^{-3}\)

Full size image

In conclusion, the separation of the bundles remains bounded away from zero as we approach the breakdown. On the other hand, the regularity seems to decrease. This indicates that the destruction of the torus does not happen because of destruction of the hyperbolicity, but just because of loss of regularity.

7.1 Similarity of the Bundle Angles for Different Rotation Numbers

One striking feature of Fig. 8 is the similarity of the position of the bundles for the two frequencies at breakdown. Of course, before breakdown they are completely different. This could be explained because in renormalization, many phenomena in the fine scale at breakdown depend only on the tail of the continued fraction of the expansion. To illustrate this similarity, Fig. 9 shows the difference in the angle of center and stable bundles reported in Fig. 8.

The paper MacKay (1983) includes a renormalization procedure that involves not only the mapping, but also the rotation number. Roughly, the basic idea of the transformation is to change the map by an iterate (related to the return map of the rotation sought), and scale the space and the parameter to breakdown. At the same time, one changes the required rotation number by applying to it the Gauss map \(\omega \mapsto 1/\omega - [ 1/\omega ]\), or equivalently, shifting the continued fraction expansion. If this renormalization map has an asymptotic behavior, which is the same for several maps/numbers (even if the dynamics of the renormalization is more complicated than a fixed point), then the properties close to the breakdown would only depend on the tail of the continued fraction expansion. This is consistent with our data.

8 Width of Analyticity Domains

Since the functions we consider are analytic up to the breakdown, a way to ascertain the breakdown is to explore when the analyticity domain goes to zero.

A widely applicable method to approximate the analyticity domain of a periodic function f is to perform a linear regression of the \(\log _{10} |\hat{f}_k| \) values and possibly discard the first terms.

In our problem, applying the above method, we obtain \(\delta _\varepsilon \) for each value \(\varepsilon \) in the continuation. The function \(f_\varepsilon \) is analytic in \(-\delta _\varepsilon \log (10)/(2\pi ) \).

Figure 10 shows the analyticity domain for the difference of the angles in Fig. 8 and also for the function \(u=u(\theta )\) defined as

$$\begin{aligned} u(\theta ) \,{\mathop {=}\limits ^{\textrm{def}}}\,K _1(\theta ) -\theta , \end{aligned}$$
(33)

where K is the spin–orbit embedding of the torus and \(K _1\) its first component.

In both cases, the values tend to zero as \(\varepsilon \) approaches the threshold value.

Fig. 10
figure 10

Maximum size of the analyticity domains of the center minus stable angles in Fig. 8 and u given in (33) for each of the frequencies (24) (left panel) and (25) (right panel) of the model (6) with fixed \(\mu =10^{-3}\), and \(e\) as drift parameter. The small plots are just zooms in

Full size image

We note that since the domain of analyticity has dimensions of length, if there are scaling relations, then the scalings should go to zero like a power of the distance of \(\varepsilon \) to the breakdown value (de la Llave 1992). Note that Fig. 10 does not support the existence of a scaling behavior.

9 Sobolev’s Criterion

The Sobolev criterion was introduced in Calleja and Celletti (2010) for symplectic mappings and extended in Calleja et al. (2013) to conformally symplectic mappings.

The rigorous underpinning is very strong. The a posteriori theorem (Calleja and Celletti 2010) gives an algorithm for the computation of the analyticity breakdown, based on the fact that if the approximate solution is well behaved, there is a true solution near the approximate one.

We start by giving the definition of Sobolev seminorm as follows.

Definition 3

Let f be a map in \(L ^2(\mathbb {T})\), whose Fourier expansion is written as \(f(\theta ) = \sum _{k \in \mathbb {Z}} \hat{f} _k e^{2\pi \texttt{i}k \theta }\), and let \(\Vert f \Vert _{L ^2} \,{\mathop {=}\limits ^{\textrm{def}}}\,\bigl ( \sum _{k \in \mathbb {Z}} | \hat{f} _k | ^2\bigr ) ^{1/2}\). Then, if r is a real number, the r-th Sobolev seminorm is defined by

$$\begin{aligned} \Vert f \Vert _r \,{\mathop {=}\limits ^{\textrm{def}}}\,\Vert \partial _{\theta }^r f \Vert _{L ^2} = \biggl ( \sum _{k \in \mathbb {Z}} |2 \pi k| ^{2r} | \hat{f} _k | ^{2} \biggr ) ^{1/2} \, \end{aligned}$$
(34)

where \(\partial _\theta ^r\) denotes the r-th derivative w.r.t. to \(\theta \) and we have used the Parseval identity to express the \(L^2\) norm of the derivative in terms of Fourier coefficients.

Note that (i) r does not need to be an integer, (ii) Definition 3 is given in terms of a seminorm on a space of (smooth) periodic functions, and (iii) a map f in Definition 3 will numerically be represented by a finite sum, say \(f^{\le L}\) where L denotes the mesh size, see Fig. 3. Therefore, we compute the seminorm of the truncation of the function f up to an integer L as

$$\begin{aligned} \Vert f ^{\le L} \Vert _r = \biggl ( \sum _{|k| \le L} (2 \pi k) ^{2r} | \hat{f} _k | ^{2} \biggr ) ^{1/2}. \end{aligned}$$
(35)

Remark 2

Note that the computation of (35) could be very sensitive to round-off error, since it involves summing terms of different sizes. It is well known that summing first the larger terms and then the smaller ones is very affected by round-off errors.

If the floating point of the computer satisfies some elementary properties (indeed holding in several popular packages such as MPFR), then there is a remarkable algorithm due to W. Kahan, which computes the sums (in a few extra operations) without round-off error (Knuth 1997, p. 244).

We have made sure that our calculations of the Sobolev seminorm are done using this algorithm. (It is implemented by default in some linear algebra packages.)

At the basis of the Sobolev criterion, there is the following remark. The KAM theorem developed in Calleja et al. (2013) for conformally symplectic systems shows that if we obtain a solution with a bounded Sobolev seminorm high enough and the bundles are separated, then one can continue further. The paper Calleja et al. (2013) provides explicit rigorous estimates of what is the meaning of high enough order. In the numerical estimates stage, one can see that the blowup happens even for orders smaller than what the rigorous results ensure.

Therefore, when the parameters approach the threshold, either all the Sobolev seminorms of the conjugacy of sufficiently high order blow up or—inclusive or—we have that the distance between the stable and unstable bundles goes to zero.

It is important to remark that the method and its justification remain valid for any number of dimensions as well as for other problems such as non-twist circles (González et al. 2022; Calleja et al. 2021).

This method was originally implemented in Calleja and de la Llave (2010) for symplectic maps, where the angle of the bundles did not enter. It was found to be very practical, since the continuation algorithms for quasi-periodic tori do not require adjustment and can run unattended. One can run different paths of continuation in different cores and obtain domains of two parameters easily. Independent implementations confirming and extending the results to asymmetric mappings appeared in Fox and Meiss (2016). Note that the asymmetric maps considered in Fox and Meiss (2014) do not have symmetry lines, so that the methods based on periodic orbits such as Greene’s method have great difficulty. Other independent implementations for twist maps appear in Flesher (2021), which explored even a three-dimensional parameter space. Implementations for the breakdown of two-dimensional tori appear in Blass and de la Llave (2013), Fox and Meiss (2016).

In numerical computations, we consider the function u in (33) and we approximate it as \(u^{\le L}\) for a suitable integer L. The truncated seminorm of \(u^{\le L}\) will be denoted as \(H _r \,{\mathop {=}\limits ^{\textrm{def}}}\,\Vert u ^{\le L} \Vert _r\). The function u depends on the frequency of the torus. Figure 11 provides the seminorm values for different indexes r and for each of the two frequencies considered in this paper, namely \(\omega _1\) in (24) and \(\omega _2\) in (25). As \(\varepsilon \) increases, the seminorms blow up, indicating a breakdown when approaching the critical \(\varepsilon _c\) values.

Fig. 11
figure 11

Sobolev seminorms as in (35) of the spin–orbit embedding u in (33) for different values r for the frequencies (24) (left panel) and (25) (right panel) of the non-averaged model (6)

Full size image

9.1 On Other Methods to Compute the Breakdown Values

There have been several methods proposed to study numerically the breakdown of invariant circles, also for higher dimensions Bustamante et al. (2023). A survey for symplectic systems appears in Calleja and de la Llave (2010); some of the methods exclude circles of a given rotation number, and some of them exclude all circles. An interesting question is to study which of these symplectic methods can be adapted to conformally symplectic systems. This adaptation is not trivial since the dissipative systems require to adjust parameters to maintain the frequency. Besides the conceptual problems of formulation and justification, there are of course practical problems.

We have considered adapting the Greene’s criterion Greene (1979) and the obstruction criterion Olvera and Simó (1986). The Greene’s criterion requires some adaptation (e.g., redefining the residue, clarifying the role of parameters (Calleja et al. 2014), while the obstruction criterion is topological and does not need adaptation. Note also that the obstruction criterion is carefully justified, while the Greene’s criterion has only partial justifications.

Both criteria depend on a detailed study of periodic orbits, whose rotation number approximates the frequency of the torus. Finding periodic orbits is not easy, even in standard-like symplectic maps, in particular when they contain several harmonics (Lomelí and Calleja 2006; Falcolini and de la Llave 1992). For dissipative mappings, Greene’s criterion has been successfully implemented in maps with few harmonics (see Calleja and Celletti 2010, which also includes a comparison with Sobolev criterion in those cases).

Certainly, the spin–orbit map contains many harmonics, there are no symmetry lines and it becomes very hard to track periodic orbits. Note also that one has to adjust the drift parameter; periodic orbits—with different characteristics—exist for drifts in intervals (Arnold tongues).

The obstruction criterion requires a careful computation of the stable and unstable manifolds of periodic orbits. We have given some effort to the implementation, but found very hard to follow periodic orbits of high period in the spin–orbit problem.

At the moment, we are not aware of any computational alternative to the Sobolev criterion, which is, moreover, rigorously justified.

The Sobolev criterion has the weak point (shared with Greene’s criterion) that the breakdown involves a limit. Hence, a finite calculation cannot give a rigorous conclusion. One can use extrapolation techniques to compute these limits as it is standard in numerical analysis. The obstruction criterion produces rigorous results of non-existence with a finite calculation (see Celletti and MacKay 2007 on the non-existence of tori of all rotation numbers). Our KAM theorems produce definite results of existence with a finite calculation (Calleja et al. 2022).

10 Fine Properties of the Breakdown: Scaling Relations, Renormalization Group

One of the consequences of our calculations (see Sect. 5) is that we can explore the phenomena that happen at breakdown of KAM theory, which is an area full of open mathematical problems.

As a consequence of our detailed calculations, we find that the breakdown of invariant circles in the spin–orbit problem does not fit with the customary scaling theory, documented in the breakdown of other mappings (Rand et al. 1982; Östlund et al. 1983).

The emphasis of this paper is in the accurate computation of tori in a concrete spin–orbit problem, and the behavior at breakdown is a (welcome but unexpected) by-product. Therefore, in this paper, we will not investigate systematically other properties such as universality (i.e., the behavior of several maps) which are of great importance to the renormalization community, but which require a point of view different from the one considered here.

10.1 Scaling Theory and Renormalization Group

In this section, we present a very quick overview of scaling theory and renormalization group. The presentation is informal, and we omit many important considerations such as domains of definitions and scaling properties of angles.

A very important discovery (Rand et al. 1982; Rand 1992a, b, c; MacKay 1983; Kadanoff 1981; Shenker and Kadanoff 1982; Feigenbaum et al. 1982) in the 1980s was that, for the systems studied, the breakdown of KAM tori had similarities with the theory of phase transitions in statistical mechanics: It satisfies scaling relations with universal exponents (i.e., the exponents appearing in the scaling relations are the same for similar systems). These scaling properties and the universality of the exponents are explained and predicted by properties of the dynamics of a renormalization group (a dynamical system in an infinite-dimensional space of maps).

If we have a breakdown of a KAM torus for \(\varepsilon = \varepsilon _c\), let us denote by \({\tilde{\varepsilon }}\,{\mathop {=}\limits ^{\textrm{def}}}\,\varepsilon - \varepsilon _c \). Scaling theory tells us that scaling the parameter is asymptotically the same as scaling the conjugacy. More precisely, for some numbers \(\delta \) and \(\eta \), we have

$$\begin{aligned} K_{\delta ^{-1} {\tilde{\varepsilon }}}(\theta ) \approx K_{\tilde{\varepsilon }}(\eta \theta ) \end{aligned}$$
(36)

and there is a limit function \(K^*\), such that

$$\begin{aligned} K_{\delta ^{-n} {\tilde{\varepsilon }}}(\eta ^n \theta ) \rightarrow K^*(\theta ), \quad \text {as } n \rightarrow \infty . \end{aligned}$$
(37)

Furthermore, the scaling factors \(\delta \) and \(\eta \) are independent of the family considered; besides, the \(K^*\) that can appear considering different families (or taking different initial \(\varepsilon \)) belong to a one-dimensional family of mappings.

This scaling behavior can be explained by assuming that a renormalization acting on mappings has a non-trivial fixed point with a one-dimensional unstable manifold and a codimension 1 stable manifold. The main object of investigation in our paper is the scaling behavior, which is a consequence of dynamical assumptions of a renormalization operator.

We omit the details of the well-known derivation of scaling behaviors from the renormalization picture, since they are not relevant for our discussion. We just recall that renormalization theory introduces an operator in the spaces of maps, which is just a change of scale in time and in space. Assuming that there is a fixed point of this operator, the scaling factor \(\delta \) is the unstable eigenvalue at the fixed point, the scaling \(\eta \) is given by properties of the fixed point and the set of scaling limits is the unstable manifold. The different maps at the threshold for many families form the stable manifold of the fixed point under the renormalization dynamics.

The renormalization analysis also predicts other scaling relations, such as the scaling properties of the analyticity domain of the solutions (see Sect. 10.2).

The results reported in Sect. 10.2 show that the spin–orbit model does not seem to satisfy the scaling relations (36).

As a consequence of our computations, we conclude that the dynamics of the renormalization group in the neighborhood of the breakdown for the spin–orbit map is more complicated than the simple dynamics observed before.

Indeed, Figs. 12 and 13 are compatible with the renormalization group having some type of non-trivial recurrence (e.g., a limit cycle or a homoclinic orbit).

The results in this paper are completely compatible with the previous papers reporting other behaviors. In renormalization theory, the universality of the exponents holds only on some universality domain (an open domain in the space of maps for which the dynamics of the renormalization is described by a hyperbolic fixed point). Indeed, one of the questions raised in MacKay (1983), MacKay (1993) is precisely to explore the boundary of the universality domain.

To obtain results in maps with many harmonics, it was crucial to have the method based on Sobolev criterion, which is very robust. Other methods, such as Greene’s method, seem to depend very well in having periodic orbits working in a predictable way even up to close to breakdown. This seems to be linked to the renormalization having a simple dynamics (de la Llave and Olvera 2006).

Since the main goal of this paper is to develop methods to compute extremely accurately until breakdown in concrete models (where the computation is challenging), we have not done a systematic study of easy to implement maps in a neighborhood. This study is natural from the point of view of renormalization group, but it goes in a direction different from the main goal of this paper.

10.2 Scale-Invariant Observables

In this section, we look at the scale invariance properties of the Sobolev seminorms. Given the relevance of scalings, it will be important for us to perform measurements on the approach to breakdown which are scale invariant.

Following Calleja and de la Llave (2010), we observe that, by Parseval identity, if \( \beta > 0, \eta \in \mathbb {N}\), and \(u _\eta (\theta ) \,{\mathop {=}\limits ^{\textrm{def}}}\,u(\eta \theta )\), then

$$\begin{aligned} \begin{aligned} \Vert \beta ^{-1} u _\eta \Vert ^2_r&= \beta ^{-2} \int _0^1 |D^r u(\eta \theta ) |^2 \, d\theta = \beta ^{-2}\int _0^1 \eta ^{2r-1} |(D^r u) (\eta \theta ) |^2 \, d (\eta \theta )\\&= \beta ^{-2} \eta ^{2 r -1} \int _0^\eta | D^r u(\sigma )|^2 \, d \sigma = \beta ^{-2} \eta ^{2r} \int _0^1 | D^r u(\sigma )|^2 \, d \sigma \\&= \beta ^{-2} \eta ^{2r} \Vert u \Vert ^2_r. \end{aligned} \end{aligned}$$

Even if the derivation above is mathematically exact only when \(\eta \) is a large number, similar scaling relations are true in the leading order for large \(\eta \).

Therefore, if we consider observables which are the products of N Sobolev seminorms, we see that, under scalings they transform as:

$$\begin{aligned} \prod _{j=1}^N \Vert \beta ^{-1} u _\eta \Vert _{r_j}^{\gamma _j} = \eta ^{ \sum _{j = 1}^N \gamma _j r _j}\cdot \left( \beta ^{-1} \eta ^{-1/2} \right) ^{ \sum _{i = 1}^N \gamma _i}\cdot \prod _{j = 1}^N\Vert u \Vert _{r _j} ^{\gamma _j}. \end{aligned}$$
(38)

Therefore, if we chose \(\gamma _j\), \(r _j\) such that

$$\begin{aligned} \sum _{i = 1}^N \gamma _i r _i = 0 \quad \text {and} \quad \sum _{i = 1}^N \gamma _i = 0, \end{aligned}$$
(39)

then the observables (38) are invariant under the scaling no matter what are the scaling parameters \(\eta ,\beta \). We refer to these observables whose value does not change under scaling as scale-invariant observables.

If scaling theory applies, we must see that the scaling invariant observables converge to a limit when we approach the breakdown.

One can think of the scaling invariant observables as some coordinates in the space of mappings. If the standard renormalization group picture (the breakdown corresponds to the stable manifold of a non-trivial fixed point) applies, one should see that near the breakdown all the scale-invariant variables accumulate near the values corresponding to the scaling limits (the one-dimensional unstable manifold of the renormalization group).

When the breakdown is described by a fixed point of the renormalization, the ratios of the scaling with \(\varrho =0\) tend to a fixed value monotonically. However, that is not the case for the spin–orbit model as it is illustrated in Figs. 12 and 13.

Fig. 12
figure 12

Scale-invariant ratios of the Sobolev’s seminorms \(H_r\equiv \Vert u^{\le L}_r\Vert \) of the non-averaged spin–orbit problem (6) for \(\mu =10^{-3}\) and frequencies (24) and (25) versus the continuation parameters. The dashed lines represent the same continuation with smaller maximum stepsize

Full size image
Fig. 13
figure 13

Scale-invariant ratios of the Sobolev’s seminorms (34) for the non-averaged spin–orbit (6) for \(\mu =10^{-3}\) and frequencies (24) and (25). Each coordinate line represents a ratio, so that a parameter value produces a point in \(\mathbb {R}^3\). Note that the lines joining the points are a product of the plotting program. They indicate the progression of the continuation parameter. The dashed lines represent the same continuation with smaller maximum stepsize

Full size image

11 Conclusions

The recent results in (2022, 2022) show that it is possible to compute irrational KAM tori in the spin–orbit problem using very efficient and accurate algorithms. These algorithms are backed up by rigorous a posteriori theorems, so that the results are convincing. Note that the problem is very challenging because it is in the regime of weak dissipation, which is a singular perturbation of the conservative spin–orbit problem.

The challenge faced in this work is to take the algorithm to the limits of validity, and see how one can maintain the reliability and explore the phenomena that happen.

The main conclusion is that the algorithm maintains the extremely high accuracy and reliability very close to breakdown.

The calculations of the breakdown in this paper have been based on the Sobolev criterion of Calleja and de la Llave (2010) for symplectic mappings (adapted in Calleja et al. 2013 to conformally symplectic mappings). We have found it easy to implement and reliable.

One of the results of our exploration is that the breakdown happens in ways that are qualitatively different to previous explorations. Notably, we do not observe scaling behaviors and the stable and tangent directions of the tori remain separated. The main reason why the tori disappear is just the loss of regularity.

The different behaviors at breakdown can be explained by postulating that the renormalization map has different dynamics in different regions in the space of maps. This paper includes tools such as accurate computation and diagnostics like scale invariants that open the possibility of a more systematic study of the dynamics of renormalization, but this will require studying breakdown in many families of maps even if they are not relevant for astrodynamics.

The (very interesting) goal of exploring the renormalization theory will be postponed, since it is different from the main goal of this paper: maintaining high precision and reliability in astronomically motivated models even very close to breakdown.