Accurate computations up to break-down of quasi-periodic attractors in the dissipative spin-orbit problem

We consider a Celestial Mechanics model: the spin-orbit problem with a dissipative tidal torque, which is a singular perturbation of a conservative system. The goal of this paper is to show that it is possible to compute quasi-periodic attractors accurately and reliably for parameter values extremely close to the breakdown. Therefore, it is possible to obtain information on mathematical phenomena at breakdown. The method we use incorporates the same time numerical and rigorous improvements. Among them (i) the formalism is based on studying the time-one map of the spin-orbit problem (which reduces the dimensionality of the problem) and has mathematical advantages; (ii) very accurate integration of the ODE (high order Taylor methods implemented with extended precision) for the map at its jets; (iii) a very efficient KAM method for maps which computes the attractor and its tangent spaces ( quadratically convergent step with low storage requirements, and low operation count); (iv) the algorithms are backed by a rigorous a-posteriori KAM Theorem, which establishes that if the algorithm, produces a very approximate solution of functional equation with reasonable condition numbers. then there is a true solution nearby; and (v) the continuation algorithm is guaranteed to reach arbitrarily close to the border of existence if it is given enough computer resources. As a byproduct of the accuracy that we maintain till breakdown, we study several scale invariant observables of the tori used in the renormalization group of infinite dimensional spaces. In contrast with previously studied simple models, the behavior at breakdown of the spin-orbit problem does not satisfy standard scaling relations which implies that the spin-orbit problem is not described by a hyperbolic fixed point of a renormalization operator.


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 ( [CCGdlL22a,CCGdlL22b]), 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) as well as 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 admits 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). Computing by direct iterations requires long transients.
Note that numerical methods based on contractions (e.g. the numerical methods suggested in [CCFdlL14]) slow down in the limit of small dissipation, whereas the Newton methods we use maintain the speed in the limit of small dissipation.
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 Sections 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 [CCdlL13] can be used in very practical situations and effectively cooperate with elaborated numerical methods; (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 blow up when the parameters approach the breakdown. We also study the behavior of the stable and tangent bundles; in [CF12], 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 [CdlL10b,CCdlL13] which show that the algorithms can compute arbitrarily close to breakdown and that the breakdown is signaled either by a blow up of Sobolev seminorms of the conjugacy or a loss of hyperbolicity (in conformally symplectic systems, the stable Lyapunov exponent is bounded away from zero, so that the loss of hyperbolicity can only happen because the minimum of 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 (with log-periodic corrections) 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 failure of the scaling behaviour found in other models with few harmonics can be described in the renormalization group language as saying that the spinorbit problem is so far away from previously studied models that it is in a different universality class. Indeed, [Mac83] asked the question of clarifying the limits of the universality class.
It seems quite likely that some of the phenomena discovered here for the spinorbit problem should be applicable to more models (they are 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) (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 achieved only by using extended precision in the arithmetic and in the integration step. The only way to achieve sufficient precision in the integration using reasonable size steps is to use a very high order method (e.g. order 30).
Moreover, the KAM algorithm uses the derivatives with respect to initial conditions and parameters of the time one map, which requires integrating the variational equations.
All these requirements (very high precision calculations for the map and the variational equations) can be obtained simultaneously with a reasonable amount of programming time using the public domain taylor package [JZ05,GJZ22].
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, pow, etc. and using composition and arithmetic operations.
An important recently added feature of taylor described in [GJZ22] is that it implements what is called jet transport and it provides high-order derivatives of the time-t map. Jet transport is a technique that considers the phase space variables as multivariate polynomials instead of real numbers. As it is proved in [GJJC + 23] the numerical solution of the ODE with polynomial variables is equivalent to solving the original equation and its variational equations. The Taylor method for polynomials can be obtained by overloading the polynomial types in the arithmetic operations and one can also adapt step control strategies.
The upshot is that, using taylor, with a programming effort comparable to typing the equations in T E X, one obtains an arbitrary order solver for the equations and the variational equations, implemented in extended precision. As we will see, the running times are very competitive.
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. Higher precisions could be accomplished by changing a few parameters in the run at the price of 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.
The method is based on formulating a functional equation for the parameterization of a torus that expresses invariance and the motion is a rotation. The unknowns to be solved for are the parameterization and the parameters of the model that need to be adjusted.
We use a Newton-like method based on automatic reducibility developed for symplectic mappings in [dlL01,dlLGJV05] and for conformally symplectic systems in [CCdlL13].
Note that the algorithm only needs to manipulate functions of as many variables as the dimension of the torus in the return map. This number of variables is smaller than the dimension of the phase space (other methods based on normal forms or expansions use functions in phase space or even families of such functions). Reducing the number of variables of the functions considered is important because the computational cost of handling a function grows extremely rapidly with the number of variables (the curse of dimensionality).
The algorithm we use 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 [CCdlL13]) 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. Using more traditional languages, one can implement a layer of operations among objects and then a sort main program. 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 consideration is the rigorous result of [CCdlL13] stated in an aposteriori 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 [CCGdlL22b]. This has been reconsiderd in [Lin23].
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 [dlLR91]. 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, [FH12] developed a similar methodology to compute reliably close to breakdown based on efficient algorithms backed up by a-posteriori theorems. This was used in [FH15] 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 Section 2 and its Poincaré map is described in Section 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 (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 Section 4).
The computation of the tori and their exploration near breakdown is presented in Section 5, while Section 6 presents the scaling theory of breakdown of tori (a growing part of the theory has been made rigorous by now).
The construction of the invariant bundles and an exploration of the angle between the stable and tangent bundles, is given in Section 7. Moreover, Section 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 Section 9.
The links with renormalization theory and the study of scale-invariant properties of the tori are studied in Section 10.
Finally, some conclusions are presented in Section 11.

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 ( [Pea05,Cel10]).
The spin-orbit problem describes a simplified model for the rotational dynamics of a satellite, say S, orbiting around a central planet, say P, and rotating around an internal spin-axis ([Bel01, Cel90b, Cel90a, CL04, WPM84]). We assume that the satellite is a triaxial ellipsoid with principal moments of inertia A < B < C. Moreover, we make the following simplifying assumptions: A1. The barycenter of the satellite S moves on an elliptic Keplerian orbit with the planet P in one focus; the Keplerian orbit is characterized by a semimajor axis a and an eccentricity e.
A2. The satellite spin-axis coincides with the direction of the smallest physical axis of the ellipsoid.
A3. The spin-axis is perpendicular to the orbital plane.
A4. The satellite is assumed to be non-rigid, hence subject to a tidal torque.
A5. The tidal torque is a linear function of the angular velocity of rotation ( [Mac64,Pea05], see also [CL14,CL04]). The unit of time is chosen to normalize the orbital period T orb to 2π, which implies that the mean motion n = 2π/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 where the perturbative parameter ε 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 (2) The right hand side of (1) represents the dissipative effect, which is given by the tidal torque T d = T d ( dx dt (t), t); using the model developed in [Mac64,Pea05] (see Assumption A5), T d takes the form where µ > 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 ε 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(ε 2 ).
For µ = 0, equation (1) is conservative, and for ε = µ = 0 it is integrable. Thus, for small values of ε and µ, we obtain a nearly-integrable and nearly-Hamiltonian system.
Taking the average of the tidal torque over one orbital period (see, e.g., [Pea05,CL04]), equation (1) becomes whereL (e) = 1 (1 − e 2 ) 9/2 1 + 3e 2 + 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 equation (1) as a first order system in phase space, we obtain:ẋ = ẏ and similarly for equation (4): Equations (1) and (4) are defined over the phase space [0, 2π) × R, which can be endowed with the standard scalar product and the symplectic form Ω = dy ∧ dx. Both equations (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 where Tr denotes the trace and A(t) is the differential of the vector field at the time t flow.
With reference to equation (6), we obtain Tr(A(t)) = −µ a r(t) 6 and consequently the volume contracts as For the averaged tidal torque model (7), we obtain and consequently the volume contracts as It is not difficult to see by an explicit computation (compare with [CCGdlL22a]), that at t = 2π both (8) and (9) give the same asymptotic contraction rate of the volume.
Note that, when we do continuation in ε, 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 Section 3, one of the advantages of the approach of [CCGdlL22a] 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 [CCdlL13] shows that, when the maps are close, the attractors of the maps are close. Hence, the a-posteriori results, that justify an approximations for finite time, justify results for objects over all times.
Besides the general argument above, [CCGdlL22a] 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 Figure 1 and Figure 5.

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 to the vector field (6) by a return map. According to the procedure detailed in [CCGdlL22a], we consider the Poincaré map P e associated to (6), which is obtained as follows. We introduce the map P e as where x(2π; x 0 , y 0 , ε) and y(2π; x 0 , y 0 , ε) represent the solution of (6) at time t = 2π with initial conditions (x 0 , y 0 ) at t = 0. Writing P e in components as P e ≡ (P (1) e , P (2) e ), the spin-orbit Poincaré map is given bȳ e (x, y; ε) .
For computational reasons, it is convenient to consider the change of coordinates and define the map 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.
where f * ϑ denotes the pull-back. The quantity λ is called the conformal factor. For λ = 1 we recover the symplectic case. For n = 1 any diffeomorphism is conformally symplectic with λ(x) = ±| det(Df ϑ (x))|. For n ≥ 2 it is known that λ is constant (see [CCdlL13] for details).
For later use, we denote by J the symplectic matrix representing Ω at z, namely For the spin-orbit map, the matrix J is given by The Poincaré spin-orbit map (11) (or (13)) is conformally symplectic ([CCGdlL22a]) with the conformal factor given by The system contracts the volume for µ > 0, expands the volume for µ < 0, and it is neutral for µ = 0. In the following sections, we will only take µ > 0. Note that the conformal factor of the averaged model (4) will have the same λ given in (16), since at time t = 2π both contraction rates (8) and (9) coincide precisely due to theL(e) choice. For the calculation leading to (5) we refer to [CCGdlL22a,p. 15].

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 to a given frequency of an invariant torus. Such information is essential to start the procedure detailed in [CCGdlL22a] for the computation of a KAM torus with fixed frequency.
The phase space associated to the spin-orbit map (11) is a subset of the cylinder [0, 2π)×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) def = sin(2x(t) − 2f (t)), which depends on u and e through f , that satisfies the following identities: We use simple trigonometric identities to write the function s as and we define the function c as 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: As a consequence of the change of coordinates, we deduce that The spin-orbit problem in (6) can then be expressed in terms of the independent variable u as where β and γ are defined as The time-(2π) map associated to (20) is related, up to a factor, to the map defined in (13).

Computation of the rotation numbers.
In this section, we present the computation of the rotation numbers obtained by a direct iteration; in Section 4.2, we will map the regions covered by different attractors. First, we recall that the definition of the rotation number associated to equation (6) (equivalently (7)) is: 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 [DSSY17], a method alternative to (21) to compute the rotation number is through the following formula: where φ denotes a weight function, e.g., It is shown in [DSSY17] 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 superexponentially 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 to an attractor close to the initial condition (x 0 , y 0 ) ∈ [0, 2π) × R, can be performed according to Algorithm A, described below 1 .
Algorithm A. Computation of the rotation number for the spin-orbit problem with tidal torque, see equation (6).
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 [CL14]), 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 paralellizable.
We report in Figure 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: and with γ ± g equal to the golden ratio and its conjugate: 2 . Note that both ω 1 , ω 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 [Mac83] 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 ω 1 and ω 2 with two different dashed lines. Moreover, it shows with a dotted-dashed line the termN (e)/L(e) (see (5)), which is the predicted rotation number in the model with the averaged tidal torque, see [CC09,CCGdlL22a]. 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 to the averaged and the non-averaged models, Figure 1 shows also a zoom-in nearby the frequency ω 2 . This shows that the rigorous results in [CCGdlL22a] justifying the averaging method give very good numerical values.

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 × n y (for some positive integers n x , n y ) in a window of the cylinder [0, 2π) × 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 [DM98]) or, working on different regions, could be distributed using coarse grained distributors such as GNU-parallel or HTCondor (see [Tan11,htc]). We have used these techniques systematically.
Although we have not yet used GPU's, the computation of these basins are indeed suitable for GPU's. For preliminary calculations, the single precision available even in consumer grade GPU's is perfectly acceptable. However, when dealing with the more detailed calculations of tori, even double precision is not enough and we need extended precision.  (23), n 1 = 4500, n 2 = 4600, δ = 10, ε = 10 −4 , µ = 10 −5 and (x 0 , y 0 ) = (0, 1.38). The fixed frequencies ω 1 and ω 2 are given in (24) and (25), the predicted rotation number for the averaged spin-orbit is given bȳ N(e)/L(e). Finally a zoom-in is shown nearby ω 2 with the + connected by a straight line.
around 12 minutes 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 ω 1 in the upper plots of Figure 2 and ω 2 in the lower plots, become more irregular as the perturbing parameter ε 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 [CCdlL13].

Explorations of the breakdown
A rotational torus defined for a conformally symplectic system (see Definition 1) is introduced as follows.
Definition 2. Let f ϑ : M → M be a conformally symplectic system and let ω be an irrational number. A map K : T n → M is said to be a rotational torus for f ϑ with frequency ω, when there exists ϑ 0 such that The algorithm in [CCdlL13] deals with tori whose frequency has Diophantine properties and it computes iteratively corrections for the map K and the parameter ϑ 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 an 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 towards the breakdown. The theorems and algorithms in [CCdlL13] provide a way to efficiently compute KAM tori in conformally symplectic systems. Those results have been adapted for the spin-orbit problem in [CCGdlL22a]. 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 ( [Arn65,CCFdlL14]).
The method detailed in [CCdlL13] 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 (we take the solution for a parameter value as an approximate solution for a slightly bigger value of the parameter 2 ).
Furthermore the continuability argument in [CCdlL13] shows that the continuation method will get as close as desired to the boundary of existence if given enough computer resources (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 [dlLR91]).
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 [CCdlL13] 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 [JZ05], 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.20GHz with 18 CPUs with hyperthreading which simulates 36 CPUs, while for (25) a machine with Intel Xeon Gold 6154 CPU at 3.00GHz 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 [CCGdlL22a] 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 ω 1 , e ω 1 ) and (K ω 2 , e ω 2 ) such that We fix the dissipation parameter to µ = 10 −3 and we perform a continuation w.r.t. the parameter ε. 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, see below).

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 equation (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 ǫ 1 ≤ ǫ 2 . If the size of the tail is bigger than ǫ 2 , we increase the number of coefficients used; if the tail is smaller than ǫ 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 meshsize.

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.
For the spin-orbit model, we use 70 digits of precision, 10 −35 of Newton tolerance for the error of the invariance equation (26), tail tolerance between ǫ 1 = 10 −55 and ǫ 2 = 10 −28 , and 10 −70 as tolerance for the absolute and relative errors in the Taylor integration.
The corresponding tori are plotted in Figure 4 and the values of the drift parameter for ω 1 and ω 2 are: e c ω 2 = 2.4797274489383016717182991353216456e-01 .  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, Figure 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., Figure 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 Figure 3 provides the time of each successful continuation step in ε. 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 θ 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.

Studies of the breakdown of tori
Our results on the behavior at breakdown of the spin-orbit map are described in Sections 7 and 9, where we have taken the algorithms to their limits of validity. Our study is motivated by different reasons: (i) the spin-orbit problem provides an excellent testing ground for the algorithms, which is far from the academic standard models; (ii) the a-posteriori method allows to compute with confidence very close to the breakdown; (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 ( [ORSS83,ROSS82]). 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, [Mac83,ROSS82].
6.1. Some rigorous results on the breakdown. The rigorous result in [CCdlL13] has a remarkable consequence concerning the behavior of the torus, and of the stable and tangent bundles at breakdown. Precisely, if: (i) the invariant attractor is smoothly conjugate to a rotation, (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 [CCdlL13] 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 (Section 7), the regularity of the tori (Section 8), and the regularity in Sections 9. Since previous studies found scaling relations, we have also studied the possibility of scaling invariant ratios (Section 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 ( [HdlL06b,HdlL06a]) and called bundle collapse. In [HdlL06b] it was also found that this phenomenon happens often and that it satisfies remarkable scaling relations. The paper [HdlL06b] contains a proof that this phenomenon indeed happens in some examples. In [BS08,FOT20] 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 [Fen74] show that the torus and the bundle are C r for any r ∈ N. • Using [Her79,Yoc84,KS86], we obtain that the dynamics is C r−a conjugate to a rotation for some a ∈ 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 [CCdlL13] 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 ( [YdlL21]).
We also note that, even after the breakdown, there could be other invariant objects of lower regularity. For example, [CK20] presents a topological extension of normal hyperbolicity beyond C 1 regularity (the objects produced are continua, not manifolds). For dissipative systems, there are topological ( [LC86,LC88]) or variational arguments ( [MS17]), 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 [HdlL06b,HdlL06a] and, specially [CF12], 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.

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 ω for the map P e given in (11). Let N(θ), M(θ), S(θ) be the quantities defined by where T ω denotes the shift by ω: T ω (θ) = θ + ω. Then, one can show ( [CCdlL13]) that N, M, S, and P e satisfy the relation (16). Similarly, the stable and tangent bundles, say E s (θ) and E c (θ) respectively, must satisfy the reducibility equation (see [CF12]) To reduce the cocycle, we introduce the change of variables Hence, the unknown function B(θ) verifies the following cohomology equation: which can be solved by expanding in Fourier series whenever |λ| = 1. Finally, the tangent and stable bundles can be recovered from the relation W (θ) = M(θ) W (θ). Explicitly, if K has components (K 1 , K 2 ), then A straightforward computation of the inner products between the bundles leads to Therefore, if α(θ) is the angle between the invariant stable and tangent bundles at the point θ ∈ T, then cos α(θ) = B(θ) (32) Figure 7 shows the behavior of the angle α 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 ε approaches the breakdown threshold.  In [CF12] 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 ε 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 [CF12]. Moreover, the minimum value happens in the same value around θ ≈ 0.5 which is coherent with the results in [BS08].
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 Figure 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, Figure 9 shows the difference in the angle of center and stable bundles reported in Figure 8.  Figure 9. Quotient between the difference in the angle of the center and stable bundles for frequencies (24) and (25)  The paper [Mac83] 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 ω → 1/ω − [1/ω], 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.

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 analiticity domain of a periodic function f is to perform a linear regression of the log 10 |f k | values (one has to discard the first terms).
In our problem, applying the above method, we obtain δ ε for each value ε in the continuation. The function f ε is analytic in −δ ε log(10)/(2π). Figure 10 shows the analiticity domain for the difference of the angles in Figure 8 and also for the function u = u(θ) defined as 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 ε approaches the threshold value. 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 ε to the breakdown value [dlL92]. Note that Figure 10 does not support the existence of a scaling behavior.

Sobolev's criterion
The Sobolev criterion was introduced in [CdlL09] for symplectic mappings and extended in [CCdlL13] to conformally symplectic mappings.
The rigorous underpinning is very strong. The a-posteriori theorem ([CdlL09]) gives an algorithm for the computation of the analiticity 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 (T), whose Fourier expansion is written as f (θ) = k∈Zf k e 2πikθ , and let f L 2 def = k∈Z |f k | 2 1/2 . Then, if r is a real number, the r-th Sobolev seminorm is defined by where ∂ r θ denotes the r-th derivative w.r.t. to θ 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) the 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 ≤L where L denotes the meshsize, see Figure 3. Therefore, we compute the seminorm of the truncation of the function f up to an integer L as Remark 2. Note that the computation of (35) could be very sensitive to roundoff 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 ([Knu97, 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 [CCdlL13] 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 [CCdlL13] provides explicit rigorous estimates of what is the meaning of high enough order (in numerical estimates, one can see that the blow up 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 nontwist circles ( [GHdlL22,CCH21]).
This method was originally implemented in [CdlL09] 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 [FM16] (note that the asymmetric maps considered in [FM14] 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 [Fle21], which explored even a 3 dimensional parameter space. Implementations for the breakdown of two dimensional tori appear in [BdlL13,FM16].
In numerical computations, we consider the function u in (33) and we approximate it as u ≤L for a suitable integer L. The truncated seminorm of u ≤L will be denoted as H r def = u ≤L 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 ω 1 in (24) and ω 2 in (25). As ε increases, the seminorms blow-up, indicating a breakdown when approaching the critical ε c values. 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 [BCL23]. A survey for symplectic systems appears in [CdlL10b]; some of the methods exclude circles of a given rotation number, 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 [Gre79] and the obstruction criterion [OS87]. The Greene's criterion requires some adaptation (e.g. redefining the residue, clarifying the role of parameters (these extensions were accomplished in [CCFdlL14]), 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 ( [LC06,FdlL92]). For dissipative mappings, Greene's criterion has been successfully implemented in maps with few harmonics (see [CC10], 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 (calculations involving extrapolations are standard in numerical analysis, though). The obstruction criterion produces rigorous results of non-existence with a finite calculation (see [CM07] on the non existence of tori of all rotation numbers). Our KAM theorems produce definite results of existence with a finite calculation ( [CCGdlL22b]).

Fine properties of the breakdown: scaling relations, renormalization group
One of the consequences of our calculations (see Section 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 [ROSS82,ORSS83].
The emphasis of this paper is in the accurate computation of tori in a concrete spin-orbit problem, the behavior at breakdown is a (welcome but unexpected) byproduct. 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, scaling properties of angles, etc.
A very important discovery ([ROSS82, Ran92a, Ran92b, Ran92c, Mac83, Kad81, SK82, FKS82]) in the 1980's 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 ε = ε c , let us denote byε def = ε − ε c . Scaling theory tells us that scaling the parameter is asymptotically the same as scaling the conjugacy. More precisely, for some numbers δ and η, we have and there is a limit function K * , such that Furthermore, the scaling factors δ and η are independent of the family considered; besides, the K * that can appear considering different families (or taking different initial ε) belong to a 1-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 1-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. As a matter of fact, (36) becomes exact in the unstable manifold of renormalization.
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 δ is the unstable eigenvalue at the fixed point, the scaling η 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 Section 10.2).
The results reported in Section 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, Figures 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 [Mac83,Mac93] 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 [dlLO06].
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 [CdlL10a], we observe that, by Parseval identity, if β > 0, η ∈ N, and u η (θ) Therefore, if we consider observables which are the products of N Sobolev seminorms, we see that, under scalings they transform as: Therefore, if we chose γ j , r j such that N i=1 γ i r i = 0 and then, the observables (38) are invariant under the scaling no matter what are the scaling parameters η, β. 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 are log-periodic in the parameters in the unstable manifold. (Compute a scale invariant observable in the unstable manifold where (36) is exact. We see that they should satisfy φ(δε) = φ(ε). For arbitrary families, the λ-lemma says that the log-periodicity should be asymptotic near 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 ̺ = 0 tend to a fixed value monotonically. However, that is not the case for the spin-orbit model as it is illustrated in Figures 12 and 13 (1 − ε/ε c ω2 ) −1 ω 2 Figure 12. Scale-invariant ratios of the Sobolev's seminorms H r ≡ u ≤L r of the non-averaged spin-orbit problem (6) for µ = 10 −3 and frequencies (24) and (25) versus the continuation parameters. The dashed lines represent the same continuation with smaller maximum stepsize.

Conclusions
The recent results in [CCGdlL22a,CCGdlL22b] 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 [CdlL10b] for symplectic mappings (adapted in [CCdlL13] 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 1 1.5 2 2.5 3 3.5 4 4.5 1 2 3 4 5 6 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 Figure 13. Scale-invariant ratios of the Sobolev's seminorms (34) for the non-averaged spin-orbit (6) for µ = 10 −3 and frequencies (24) and (25). Each coordinate line represents a ratio, so that a parameter value produces a point in 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.
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.