The Spin–Spin Problem in Celestial Mechanics

We study the dynamics of two homogeneous rigid ellipsoids subject to their mutual gravitational influence. We assume that the spin axis of each ellipsoid coincides with its shortest physical axis and is perpendicular to the orbital plane. Due to such assumptions, the problem is planar and depends on particular parameters of the ellipsoids, most notably, the equatorial oblateness and the flattening with respect to the shortest physical axes. We consider two models for such configuration: while in the full model, there is a coupling between the orbital and rotational motions, in the Keplerian model, the centers of mass of the bodies are constrained to move on coplanar Keplerian ellipses. The Keplerian case, in the approximation that includes the coupling between the spins of the two ellipsoids, is what we call spin–spin problem, that is a generalization of the classical spin–orbit problem. In this paper we continue the investigations of Misquero (Nonlinearity 34:2191–2219, 2021) on the spin–spin problem by comparing it with the spin–orbit problem and also with the full model. Beside detailing the models associated to the spin–orbit and spin–spin problems, we introduce the notions of standard and balanced resonances, which lead us to investigate the existence of periodic and quasi-periodic solutions. We also give a qualitative description of the phase space and provide results on the linear stability of solutions for the spin–orbit and spin–spin problems. We conclude by providing a comparison between the full and the Keplerian models with particular reference to the interaction between the rotational and orbital motions.


Introduction
The dynamics of two rigid bodies orbiting under their mutual gravitational attraction is a classical problem of Celestial Mechanics known as the Full Two-Body problem. In this context, Kinoshita investigated the problem by using Hori-Deprit perturbation theory [Kin72], assuming that one of the bodies is spherical and the other body is triaxial. Later, the problem of two extended rigid bodies was studied in [Mac95] as a Hamiltonian system with respect to a non-canonical structure, which is used to characterize the relative equilibria. A seminal work was performed in [Bou17] to which we refer for an alternative description of the model of the full two rigid body problem using spherical harmonics and Wigner D-matrices. In [Sch09], the problem is restricted to a planar configuration with the potential expanded to order 1/r 3 , where r is the relative distance between the two rigid bodies; under this condition, [Sch09] describes the relative equilibria and their stability properties.
In this paper, we investigate different simplified models of rotational dynamics of celestial bodies, subject to the mutual gravitational attraction. The spin-spin problem was introduced in [Mis21] as a planar version of the Full Two-Body problem for ellipsoids (compare with [BM15,BL09]), by using the expansion of the potential up to order 1/r 5 , which results in the coupling of the spins of both bodies. An equivalent model was studied in [JA16] (see also [HX17]).
Indeed, we consider a hierarchy of models with different complexity. In particular, we start by considering two homogeneous rigid ellipsoidal bodies subject to the following assumptions: Assumption 1. The spin axis of each ellipsoid is perpendicular to the orbital plane.
Assumption 2. The spin axis of each ellipsoid is aligned with the shortest physical axis of the satellite.
Assumptions 1 and 2 imply that the motion takes place on a plane. Following [Mis21], we introduce a Hamiltonian function that includes both the orbital and rotational motions. Using the conservation of the angular momentum, the system is described by a Hamiltonian with 3 degrees of freedom that depends on several parameters of each ellipsoid, among which there are the equatorial oblateness and the flattening with respect to the shortest physical axis. Such parameters are typically small for natural bodies of the solar system. We refer to this model as the full problem, since it includes the coupling between the orbital and rotational motions.
The potential of the problem can be written as V = V 0 + ∞ l=1 V 2l , where V 0 denotes the Keplerian potential and the terms V 2l are proportional to 1/r 2l+1 , where r is the instantaneous distance of the two centers of mass. If we consider the expansion up to order l = 2, say V = V 0 + V 2 + V 4 , we obtain that the model includes the coupling of the spins of the two ellipsoids through the term V 4 , which contains combinations of the rotation angles of the two satellites. When the two spins interact, we refer to the problem as the full spin-spin model. If, instead, we limit the potential to V = V 0 + V 2 , then we obtain two decoupled systems, the full spin-orbit models (see, e.g., [Bel66,Cel90,GP66]). When one of the bodies is spherical, its spin is uniform and the dynamics of the spin-spin model becomes very similar to that of the spin-orbit model, but including the terms in V 4 .
Next, we introduce another assumption, namely: Assumption 3. The orbital motion of the ellipsoids coincides with that of two point masses, so that both centers of mass move on coplanar Keplerian orbits with eccentricity e ∈ [0, 1) and with a common focus at the barycenter of the system.
Assumption 3 implies that the orbit is not affected by the rotational motion. To the full spin-spin and spin-orbit problems, it corresponds the Keplerian spin-spin and spin-orbit models, described by a Hamiltonian function with an explicit periodic time dependence.
Note that we are considering rigid bodies only, which means that dissipative effects due to tidal torques are not considered. We refer to [CC08,CC09,MO20,Mis21,GP66] for a description of the dissipative spin-orbit and spin-spin problems.
The previously described models have some symmetries that are a direct consequence of the mirror symmetries of the ellipsoids. This fact leads us to introduce the following two types of resonances within the spin-orbit problem of a single ellipsoid: (R1) We call standard m : n spin-orbit resonance, for some integers m, n, when the spinning body makes m rotations during n orbital revolutions; (R2) We call balanced m : 2 spin-orbit resonance, for some integers m, when the spinning body makes m/2 of a rotation during one orbital revolution.
We remark that a balanced spin-orbit resonance is also a standard resonance, but the converse is not always true. For example, a balanced 2k : 2 resonance for k ∈ Z is equivalent to a k : 1 spin-orbit resonance, but there is not such an equivalence for order (2k + 1) : 2. Both definitions (R1) and (R2) extend to the spin-spin problem of type (m 1 : n 1 , m 2 : n 2 ) for integers m 1 , m 2 , n 1 , n 2 , when the first ellipsoid is in a m 1 : n 1 spin-orbit resonance and the second ellipsoid in a m 2 : n 2 spin-orbit resonance.
We also stress that spin-orbit resonances find many applications in the solar system; in fact, the Moon is an example of a 1 : 1 spin-orbit resonance 1 , since it makes a rotation in the same period it takes to make an orbit around the Earth. This is also called a synchronous spin-orbit resonance, which is common to many satellites of other planets, including Mars, Jupiter, Saturn, Uranus, Neptune. Among the planets, Mercury is locked in a 3 : 2 spin-orbit resonance around the Sun. On the other hand, the Pluto-Charon system is locked in the double synchronous spin-spin resonance (1 : 1, 1 : 1).
In this work, we study the behavior of the solutions of the spin-orbit and spin-spin problems as the parameters and the initial conditions are varied. In particular, we investigate the boundary conditions that lead to the existence of symmetric periodic orbits. Such results (see Propositions 5 and 10) use some symmetry properties of the equations of motion. We remark that these symmetries are lost if we include dissipation; however, such periodic solutions might be continued to the dissipative setting as shown in [MO20,Mis21]. Beside the study of the periodic orbits, we provide the conditions for the existence of quasi-periodic solutions of the Keplerian version of the spin-spin model.
We also give a qualitative study of the spin-orbit problem as well as the spin-spin problem with spherical and non-spherical companion. Within such investigation, we discover some new features, like the measure synchronization (see, e.g., [HZ99]) for the spin-spin problem with identical bodies. Our study leads to analyze the multiplicity of solutions and the linear stability of the periodic orbits (compare with [CC00]). In general, we find that there is not a unique solution associated to a particular resonance; however, for some values of the parameters such a uniqueness exists. Finally, we provide some results on the comparison between the full and Keplerian models, as well as on the interaction between the spin and the orbital motion, motivated by the fact that the coupling between the rotational and orbital motions has not been much explored in the literature.
This work is organized as follows. In Section 2 we present the spin-orbit and spin-spin models. The definition of resonances and the existence of periodic and quasi-periodic orbits are given in Section 3. A qualitative description of the phase space is given in Section 4. The linear stability of symmetric periodic orbits is investigated in Section 5. Finally, a comparison between the full and Keplerian models is presented in Section 6.

The models for the spin-orbit and spin-spin coupling
The aim of this section is to present the so-called spin-orbit and spin-spin models, that we are going to introduce as follows. The assumptions and the notations are given in Section 2.1; different models, subject to some or all the assumptions listed in Section 2.1, are presented in Sections 2.2 and 2.3.
2.1. General assumptions. Consider two homogeneous rigid ellipsoids, say E 1 and E 2 , with masses M 1 and M 2 , respectively. Let A j < B j < C j , j = 1, 2, be their principal moments of inertia with corresponding principal semi-axes a j > b j > c j . We refer to the Full Two-Body Problem (hereafter F2BP) as the problem of two rigid bodies interacting gravitationally (see, e.g., [Sch02,Sch09]). When the bodies have ellipsoidal shape, we speak of the ellipsoidal F2BP, where we make the assumptions 1, 2, 3 of Section 1: Assumptions 1 and 2 guarantee that the problem we deal with is a planar problem. Additionally, Assumption 3 restricts the problem so that we obtain a model with two degrees of freedom and a periodic time dependence. This assumption is equivalent to say that the spin motion will not influence the orbital motion. Besides, note that we are neglecting the gravitational contribution of other bodies, we are not considering any dissipative effect that might arise, for example, from the non-rigidity of the ellipsoidal bodies and we do not take into account the obliquity, namely the inclination of the spin-axes with respect to the orbital plane.
We are going to work with units adapted to the system. If we call τ the orbital period of the Keplerian orbit, then we will use units such that Recall Kepler's third law for the Two-Body Problem where G is the gravitational constant, a is the semi-major axis associated to the motion of the reduced mass of the system, say µ = M 1 M 2 in our units. In consequence, G = a 3 in our units. Let us now define the parameters for each ellipsoid the quantity d j /C j measures the equatorial oblateness of each ellipsoid with respect to the plane formed by the directions of a j and b j , whereas q j /C j measures the flattening with respect to the direction corresponding to the c j -axis. Note that if A j ≤ B j ≤ C j , then, in our units there are some bounds for the parameters of the system given by The last relation in (1) comes from the fact that the moments of inertia of an ellipsoid hold the identities 2.2. The full models. First, let us derive the equations of the full models of spin-orbit and spin-spin coupling, for which only Assumptions 1 and 2 hold, namely we do not constrain the centers of mass of E 1 and E 2 to move on Keplerian ellipses. The equations of motion are obtained by computing the Hamiltonian function through a Legendre transformation of the Lagrangian, say L = T − V , where T is the kinetic energy and V the potential energy of the system. We split T in two parts, associated respectively to the orbital and rotational motions, say T = T orb + T rot .
Let us identify the orbital plane with the complex plane C, consider the center of mass of the system fixed in the origin and let the position of each ellipsoid be r j ∈ C. Then, by definition of the barycenter we have M 1 r 1 + M 2 r 2 = 0 .
If we define the relative position vector r = r 2 − r 1 , since in our units M 1 + M 2 = 1, then we have r 1 = −M 2 r , r 2 = M 1 r .
(2) Figure 2. Generalized coordinates of the full planar model.
We introduce the Lagrangian generalized coordinates (r, f, θ 1 , θ 2 ), illustrated in Figure 2, in the following way. The distance r > 0 and the angle f define the relative position vector r = r exp(if ) ∈ C with respect to an inertial reference frame with a fixed x-axis. The angles θ 1 , θ 2 provide the orientation of the axes a 1 , a 2 with respect to the x-axis. The variables r and f define the orbital motion of the system and θ j the spin motion of the ellipsoid E j .
To write the kinetic energy, we notice that we can express r in components as r = (r cos f, r sin f ), which givesṙ = (ṙ cos f − rḟ sin f,ṙ sin f + rḟ cos f ). Then, using (2), the total orbital and the rotational kinetic energies are given by where, according to [Mis21,Bou17], the full expansion of the potential energy of the system V = V (r, f, θ 1 , θ 2 ) takes the form where Υ = {(l, m) ∈ Z 2 : 0 ≤ |m| ≤ l} , and the constants Λ l 1 ,m 1 l 2 ,m 2 are defined in Appendix A; we refer to [Mis21] for full details. Let us define the momenta as p r = ∂ṙL = µṙ, p f = ∂ḟ L = µr 2ḟ , p j = ∂θ j L = C jθj ; then, the Hamiltonian of the system becomes H(r, f, θ 1 , θ 2 , p r , p f , p 1 , p 2 ) = p 2 Consequently, the equations of motion arė We remark that with these variables, the problem splits in two parts: equations (5) describe the orbital motion, while equations (6) describe the rotational motion. The evolution of both parts is coupled through the potential V . The potential V , expanded in (3), can be written in a perturbative way as We remark that the term V 0 corresponds to the classical Keplerian form for the potential and the term V per provides the coupling between the spin and the orbital motions.
Moreover, we can expand V per as V per = ∞ l=1 V 2l , where V 2l are suitable terms proportional to 1/r 2l+1 . A truncation of the expansion of V per will result in an approximated dynamics of our system. The explicit expressions of the first two terms of such a expansion are given in [Mis21] and we report them here: From now on, we will refer to (5) and (6) as the full models: full spin-orbit model if we take V per = V 2 , in which the angles θ 1 , θ 2 appear in different trigonometric terms, and full spin-spin model if V per = V 2 + V 4 , which contains trigonometric terms with combinations of the rotation angles θ 1 , θ 2 . These names are motivated by the well-known spin-orbit model, [Cel10], and the spin-spin model from [Mis21]. If we consider the models under Assumption 3 which gives a constraint on the orbit, we speak of Keplerian spin-orbit model and Keplerian spin-spin model (compare with Section 2.3).
Then, the Hamiltonian (4) in this new set of variables is given by where V(r, φ 1 , φ 2 ) = V (r, f, φ 1 + f, φ 2 + f ). Now it is clear that f is an ignorable variable in (10) and that P f is a constant of motion, corresponding to the total angular momentum of the system. In summary, in the evolution of the system, there is a transfer of angular momentum between the spin part, given for each body by p j = C jθj , and the orbital part, given by p f = µr 2ḟ .
2.3. The Keplerian models. In this section, we introduce the Assumption 3 to the model of Section 2.2. From (5) and (6), it is straightforward to constrain the orbit to be Keplerian; only in the orbital part we retain just the term V 0 of the potential V = V 0 + V per in (7), thus obtaining the following equations: A convenient procedure to numerically integrate the equations of motion (12) is presented in Appendix B.
2.3.1. Orbital motion. Note that since ∂ θ j V 0 = 0, the system (11) is now decoupled from (12). Moreover, (11) is the Kepler problem, whose solutions depend on the eccentricity e and the semi-major axis a of the orbit. Here we assume for simplicity that the orbit is a 2π-periodic Keplerian ellipse of eccentricity e ∈ [0, 1) with focus at the origin and with the periapsis on the positive x-axis.
Since the orbital period is 2π, then, we can take the time t to coincide with the mean anomaly. We denote by u the eccentric anomaly, which, in our units, is related to the mean anomaly by Kepler's equation The orbital radius is related to u by r = a(1 − e cos u) .
We can write the vector r ∈ C in terms of the eccentric anomaly also as r exp(if ) = a(cos u − e + i √ 1 − e 2 sin u) .
Note that for t = 0 we assumed, without loss of generality, that f = u = 0, and consequently, f = u = π when t = π. From (15) we obtain the following useful relations between f and u With the previous definitions, the Keplerian orbit of eccentricity e and semi-major axis a is given by the functions r = r(t; a, e), f = f (t; e), p r = p r (t; a, e), that correspond to the solution of equations (11) generated by the initial conditions 2.3.2. Spin motion. The spin motion is described by equations (12) with the Keplerian periodic input (17) given implicitly by eqs. (13) to (15). This motion can be described by the non-autonomous Hamiltonian H K (t, θ 1 , θ 2 , p 1 , p 2 ) = H(r(t; a, e), f (t; e), θ 1 , θ 2 , p r (t; a, e), p f (a, e), p 1 , p 2 ) , where H(r, f, θ 1 , θ 2 , p r , p f , p 1 , p 2 ) is the Hamiltonian of the full model defined in (4). The Hamiltonian (19) is hence of the form where the potential W is 2π-periodic in t and π-periodic in θ 1 and θ 2 . The equations of motion (12) take the formθ Let us define the non-dimensional parameters of the model: where λ j represents the equatorial oblateness of E j ; σ j is the ratio between the moment of inertia of E j and the orbital one; andq j measures the flattening of E j with respect to the size of the orbit. Note that the parameters in (22) are small for bodies that are close to spherical. Besides, not all the parameters defined previously are free, because we have the constraint C 1 σ 2 = C 2 σ 1 . If we take V per = V 2 , then the system (12) becomes that is a system of two uncoupled spin-orbit problems. Each of these problems depends just on two parameters: (e, λ j ). On the other hand, if V per = V 2 + V 4 , from (12) and (8) we obtain the following system for j = 1, 2, 0 =θ j + λ j 2 a r(t; e) 3 sin(2θ j − 2f (t; e))+ + a r(t; e) that we call spin-spin problem. From the previous discussion, this model depends on seven independent parameters 2 (e; C 1 , λ 1 , λ 2 , σ 1 ,q 1 ,q 2 ). Note that in (24), the coupling between the dynamics of θ 1 and θ 2 is given by σ 1 and σ 2 . Moreover, ifq j = σ j = 0, the spin-spin problem (24) is reduced to a pair of spin-orbit problems (23).
Let us now consider (24) in the case that E 2 is a sphere, that is, d 2 = q 2 = 0. Then, λ 2 = σ 2 =q 2 = 0, which implies that E 2 is in uniform rotation θ 2 (t) =θ 2 (0)t + θ 2 (0). The dynamics of θ 1 is uncoupled from θ 2 and is given by that is a spin-orbit problem up to order 1/r 5 . An equivalent system was studied previously in [JNA16]. Here the parameters σ 1 andq 1 perturb the framework of the spin-orbit problem (23).
2.3.3. Reversing symmetries. The equations for the spin motion in the Keplerian models have some reversing symmetries, i.e., transformations in the phase space that keep invariant the equations of motion with a time reversal.
Definition 1. Consider the differential equation and let x(t; x 0 ) be the solution of (26) with initial condition x(0; x 0 ) = x 0 .
One can easily check that each transformation defined by is a reversing symmetry of equation (26) with F as in (27) because The same is true replacing in (28) the quantities θ 1 ,θ 1 by θ 2 ,θ 2 . On the other hand, the spin-spin problem in its Hamiltonian formulation is given by (21), that can be written as (26) with Each transformation defined by is a reversing symmetry of equation (26).
In the following sections we are going to emphasize the study of orbits that are symmetric with respect to the mentioned reversing symmetries using the following lemma.
Lemma 2 (from Theorem 4.1 in [LR98]). An orbit of (26) is R-symmetric if, and only if, it intersects the fixed point set Fix(R).

Periodic and quasi-periodic solutions of the Keplerian models
In this section we provide the definitions of spin-orbit and spin-spin resonances, giving some results on the existence of periodic orbits (Section 3.1) and KAM tori (Section 3.2).
3.1. The spin-orbit and spin-spin resonances. The periodic solutions of the Keplerian models presented in the Section 2.3 correspond to resonances between the orbital motion and the spin motion. The expansion of the potential in (8) contains much interesting information concerning such resonances. First, let us introduce the definition of spin-orbit resonances.
Definition 3. We say that the ellipsoid E 1 is in a standard spin-orbit resonance of order The associated resonant angle ψ m:n 1 (t) = mt − nθ 1 (t) is a periodic function of period 2πn. The same definition holds for E 2 .
Recalling that we have normalized the mean motion to unity, we remark that Definition 3 states that the ratio of the orbital period of E j over its period of rotation is m/n, n = 0. Additionally, according to Definition 3, the resonance m : n, n = 0, is also of order km : kn, k ∈ Z\{0}, but the converse is not true in general. For example, with (31), the resonance 1 : 1 is also of order 2 : 2, but a resonance 2 : 2 may not be of order 1 : 1. We can say that the resonance m : n is of higher order than the resonance m : n if m/n > m /n . This will be denoted using the notation m : n > m : n .
Next, we introduce a different definition of spin-orbit resonance.
(32) In this case, the resonant angle ψ m:2 1 (t) is 2π-periodic. The same definition holds for E 2 . Notice that the two notions (31) and (32) are not equivalent: (32) implies (31) for m : 2, but the converse is not true. Actually, note that a balanced 2k : 2 resonance, with k ∈ Z, is a spin-orbit resonance of order k : 1. This new definition was motivated by [BL75], where the solutions associated to the resonance 3 : 2 of the spin-orbit problem (say, (23) with j = 1) were studied numerically. Basically, they found out that the solutions satisfying (31), but not (32), appear only for large values of λ 1 ( 1), also, for a given point (e, λ 1 ) the solutions appear in multiplets, and finally, the corresponding resonant angles have large amplitudes (|ψ 3:2 1 (t)| 0.75), see Table 1 in [BL75]. On the other hand, solutions that obey (32) exist for any point in the (e, λ 1 )-plane, including large regions of uniqueness of solution and resonant angle with small amplitude. Note that such amplitude is a measure of the deviation of the solution with respect to the uniform rotation of angular velocity 3 2 t. In Definition 4 we generalize these two types of resonances for any order m : 2, since this let us determine the main resonances in the first orbital revolution. From [DeV58], for instance, we know that a good tool to study a differential equation that has some reversing symmetries is by studying the periodic orbits that are invariant under such transformations.
In the next proposition we provide some boundary conditions that characterize the symmetric orbits in spin-orbit resonances.
(2) There are two independent types of R α,β -symmetric orbits representing a balanced m : 2 spin-orbit resonance and are given by: Type Moreover, the corresponding symmetry relations are: The same is true for (23), with j = 2 (ellipsoid E 2 ), and also for (25).
Let us now consider the balanced m : 2 case in item 2. Note that, using the definition (32) instead of (31), we can follow the same procedure as before to prove that a balanced m : 2 solution satisfies the following boundary conditions and the symmetry relation θ 1 (t) = 2β − θ 1 (2α − t). Then, we get that the four types are in this case (34), with α = 0, β = 0; (35), with α = 0, β = π/2; with α = π, β = 0, and with α = π, β = π/2. Since an ellipsoid has a mirror symmetry with respect to any plane containing a pair of semi-axes, the angle θ 1 is equivalent to θ 1 + kπ, k ∈ Z. In consequence, if m = 2k 1 , k 1 ∈ Z, type 0 is equivalent to type 0 and type π 2 is equivalent to type π 2 . Likewise, for m = 2k 2 + 1, k 2 ∈ Z, type 0 is equivalent to type π 2 and type π 2 is equivalent to type 0. Then, for resonances of order m : 2, m ∈ Z, we will take types 0 and π 2 as representatives. With this we have proved item 2. The previous facts rely only on the symmetries and the periodicity of equation (23), including the discussion in [CC00]. Since equation (25) for the spin motion of E 1 , with E 2 spherical, has exactly the same properties, then, the proof above is also valid in such case.
Proposition 5 let us characterize all the balanced spin-orbit resonances in the first half of an orbital revolution, additionally, the corresponding solutions have a certain symmetry relation. In the generalization to spin-spin resonances we want to combine different spin-orbit resonances and we will use the boundary conditions in the same time interval.
Remark 6. Note that solutions of type π 2 can be recovered with the conditions of type 0 by considering negative values of λ 1 . More precisely, solutions of type π 2 of equation (23) with λ 1 = λ * > 0 are equivalent to solutions of (23) with λ 1 = −λ * satisfying conditions of type 0. The same holds for (25).
Remark 7. Proposition 5 gives a way to numerically search for balanced resonances. Indeed, eqs. (34) and (35) can be used to apply a Newton method, that is, to find the initial conditionθ 1 (0) such that the conditions for either Type 0 or Type π 2 are satisfied. Next we introduce the following definition, which deals with the spins of both objects.
Definition 8. We say that the ellipsoids E 1 , E 2 are in a standard spin-spin resonance 3 of order (m 1 : n 1 , m 2 : n 2 ) with m 1 , m 2 ∈ Z, n 1 , n 2 ∈ Z\{0}, if the ellipsoid E j is in a m j : n j spin-orbit resonance. In such case, the resonant angles ψ m 1 :n 1 m 2 :n 2 (t) = ψ m 1 :n 1 1 (t) ± ψ m 2 :n 2 2 (t) are 2πn-periodic functions, where n is the least common multiple of n 1 and n 2 .
An analogous definition holds for resonances of the balanced type.
Definition 9. We say that the ellipsoids E 1 , E 2 are in a balanced spin-spin resonance of order (m 1 : 2, m 2 : 2) with m 1 , m 2 ∈ Z, if the ellipsoid E j is in a m j : 2 spin-orbit resonance for j = 1, 2.
The following proposition generalizes Proposition 5 to the spin-spin problem.
Proof. This proof is based on the fact that the same arguments used to prove Proposition 5 can be generalized in a straightforward way for the spin-spin problem (24). Now we apply Lemma 2 to (26) with F(x) given by (29). The fixed point set of each reversing symmetry R α,β 1 ,β 2 is Then, due to the periodicity of W in (20), the periodic orbits can be found at θ j (α) = β j , where we can take any combination between α ∈ {0, π} and β j ∈ {0, π 2 }. A similar method was used in [Gre79] for the standard map and was generalized for the spin-orbit problem in [CC00] and for a standard map of two degrees of freedom in [CFL04].
The rest of the proof of item i follows analogously to the proof of item (1) of Proposition 5 using the reversing symmetries.
The proof of item ii needs more detail. In the same way as (38), we obtain easily that the balanced resonances (m 1 : 2, m 2 : 2) are given by We see that (39) corresponds to (40) with α = 0 and the positive sign. From the case α = π and the negative sign we have that Note that if, for example, we take j = 1 and m 1 = 2k 1 , with k 1 ∈ Z, then, the conditions (39) and (41) are equivalent. On the other hand, if m 1 = 2k 1 + 1, then, the conditions (39) with β 1 = 0 and π 2 are equivalent respectively to (41) with β 1 = π 2 and 0. The same is true for j = 2. Then, with (39) all the possibilities are covered.
Remark 11. Results in Proposition 10 allow to apply a Newton method to compute resonances for the spin-spin problem just considering as unknownsθ j (0) and correct them by imposing the conditions in (39) or (41).
For circular orbits (e = 0), each of the spin-orbit models (23) is a classical pendulum whose only stable equilibrium point corresponds to a 1 : 1 resonance that is given by θ j (t) = t. Similarly, the spin-spin model (24) consists of two coupled penduli whose only stable solution is θ 1 (t) = θ 2 (t) = t that is a (1 : 1, 1 : 1) resonance. However, for e = 0, f (t; e) does not coincide with t and more stable spin-orbit resonances may appear.
In order to study the spin-orbit and spin-spin resonances, it is useful to compute the expansion of V 2 and V 4 up to some power of the eccentricity. This expansion is obtained solving Kepler's equation (13) up to a finite order in the eccentricity, then inserting the solution u = u(t) in (14), (16), expand them in series of the eccentricity and finally expanding the trigonometric terms appearing in (8).
This procedure leads to the expansions of V 2 and V 4 that, for simplicity, we give up to the order 2 in the eccentricity in the Appendix C. In those expressions, the trigonometric terms with arguments ψ m j :n j j (t) = m j t − n j θ j are associated to m j : n j spin-orbit resonances for E j , whereas the terms with argument ψ m 1 :n 1 m 2 :n 2 (t) = (m 1 ±m 2 )t−n 1 θ 1 (t)∓n 2 θ 2 (t) correspond to spin-spin resonances by combining spin-orbit resonances of orders m 1 : n 1 and m 2 : n 2 . For each order of the expansion in the eccentricity e, there are some resonances appearing. They are shown in Table 1, where we can recognize a hierarchy: the most important spin-orbit resonance is 1 : 1, then we have 3 : 2, 1 : 2 and so on, because they appear for low orders of the eccentricity. Resonances of further orders are relevant only for large eccentricities.
Note that the most relevant spin-orbit resonances are of order m : 2, m ∈ Z. Additionally, spin-spin resonances appearing at order e α in V 4 are obtained by combining spin-orbit resonances appearing at order e α 1 and e α 2 in V 2 such that α 1 + α 2 = α.
Finally, note that for V 2 , the coefficients associated to spin-orbit resonances are of order one in d 1 , d 2 , whereas for V 4 , the spin-orbit coefficients are of order two in d 1 , d 2 , q 1 , q 2 , and the spin-spin coupling coefficients are of order two in d 1 , d 2 .
the frequency vector with The Hessian matrix associated to (20) has determinant different from zero, whenever p 1 = 0, p 2 = 0. This implies that (20) satisfies the Kolmogorov non-degeneracy condition, which is a requirement for the applicability of KAM theorem [Kol54,Arn63,Mos62]. The other essential requirement in KAM theory is the assumption that the frequency satisfies a Diophantine inequality, namely there exist C > 0 and ξ ≥ 2, such that We remark that a possible choice of ω satisfying (43) can be obtained as follows. Let α be an algebraic number of degree 3, namely the solution of a polynomial equation of degree 3 with integer coefficients, not being the root of polynomial equations of lower degree. Let us consider the vector ω = (1, ω 1 , ω 2 ) obtained as where (b 1 , b 2 ) and the matrix A ≡ (a mn ) have rational coefficients and det A = 0. By number theory results (see, e.g., [CFL04]), a vector ω as in (44) satisfies (43) with ξ = 2. Under smallness conditions of the parameters, say λ j in (22), KAM theory ensures the existence of a quasi-periodic torus with Diophantine frequency. We remark that the theory presented in [CCGdlL21c] for the spin-orbit problem (see also [CCGdlL21b,CCGdlL21a]) can be extended to provide explicit estimates for (20) and an explicit algorithm to construct quasi-periodic solutions.

Qualitative description of the spin models
In this section, we give a qualitative description of the phase space associated to the spin-orbit problem (Section 4.1), the spin-spin problem with spherical companion (Section 4.2) and with non-spherical companion (Section 4.3). 4.1. Spin-orbit problem (V per = V 2 ). Since the system (23) corresponds to two uncoupled spin-orbit problems, then the Poincaré map of the whole system can be understood as the direct product of the Poincaré maps for each of the bodies. Consider for example the dynamics of E 1 . Figure 3 is a typical Poincaré map of the spin-orbit problem obtained using a Taylor integrator [JZ05] and using a similar approach as the one explained in Appendix B, in this case with parameters (e, λ 1 ) = (0.06, 0.05); it represents solutions at t = 2πk, k ∈ Z, in the plane (θ 1 ,θ 1 ) restricted to θ 1 ∈ [−π, π]. The main stable resonances are tagged with their corresponding order m : n. The Poincaré map for sufficiently small parameters (e, λ 1 ) has the following features: 1) The main stable spin-orbit resonances are represented by fixed points in the plane (θ 1 ,θ 1 ) surrounded by islands of invariant librational tori. High order resonances appear above low order resonances. 2) It looks that the spin-orbit resonances of order m : 2 ≥ 1 : 1 are balanced: solutions of type 0 (34) are stable and those of type π 2 (35) are unstable. On the contrary, for the 1 : 2 resonance, type 0 is unstable and type π 2 is stable. Concerning other resonances (33), for instance in the 3 : 4 resonance, types with α = 0 are unstable and types with α = π are stable. Exactly the opposite occurs for the case 5 : 4.
3) It is possible to have stable secondary resonances, namely small resonances surrounding other resonances. This is clear for the 1 : 1 resonance. Beyond the librational islands associated to the resonance 1 : 1 there is a chaotic region that includes the unstable resonances and that is larger for large parameters (e, λ 1 ). The chaotic region can appear for other resonances and is limited by rotational tori that also distinguish the domains of resonances of different orders. 4.2. Spin-spin problem (V per = V 2 + V 4 ) with spherical companion. Let us now consider the case when E 2 is a sphere. Then, θ 2 (t) = θ 2 (0) +θ 2 (0)t and the dynamics of θ 1 is given by (25), that depends on the parameters (e, λ 1 ,q 1 , σ 1 ). Here the parameters σ 1 andq 1 perturb the previous framework of the spin-orbit problem, see Section 4.1.
We can see a comparison between both problems in Figure 4: we see that the only qualitative difference between both cases is that the Poincaré map associated to equation (25) is slightly more chaotic. This minor difference is due to the fact that we take σ 1 =q 1 = 0.01, that are small parameters. As we will see in Section 6, the spin-spin model is a good approximation for the dynamics of two ellipsoids for σ j andq 2 up to the order of magnitude of about 10 −2 , because larger values could lead to a collision, see Section 6.2. 4.3. Spin-spin problem (V per = V 2 + V 4 ) with non-spherical companion. Now we deal with the general system (24). Let Ψ(t) = (θ 1 (t), θ 2 (t),θ 1 (t),θ 2 (t)) be a solution of (24) and its respective projections Ψ 1 (t) = (θ 1 (t),θ 1 (t)) and Ψ 2 (t) = (θ 2 (t),θ 2 (t)). From now on we restrict θ j (t) to [−π, π]. Let the Poincaré map associated to such solution be defined by P(Ψ(0)) = Ψ(2π), and its projections by P j (Ψ j (0)) = Ψ j (2π). It is not possible to represent the iteration of the map P in a single plot, because it is 4dimensional, so we will represent the projections P j . That is to say, for a solution Ψ(t), we are interested in the behavior of the two families of points Ψ 1 (2πk) and Ψ 2 (2πk), k ∈ Z, in a single 2-dimensional strip Π = {(x, y) ∈ R 2 : |x| ≤ π}. We recognize the following features: 1) A solution in spin-spin resonance corresponds to a family of isolated recurrent points of P (namely, the successive points on the Poincaré map), whose projections are represented in Π as a pair of families of recurrent points. 2) If the spin-spin resonance is stable, then nearby solutions would librate around such points. While in the uncoupled system, librating solutions belong to 2dimensional tori, here tori can be higher dimensional. As a result, the projected points represented in Π are distributed in two clouds of points surrounding each recurrent point. A cloud of this kind covers an annulus-like region of a certain thickness that is usually thicker for stronger couplings. A similar behavior occurs for rotational solutions, whose corresponding clouds are distributed in strips of certain thickness, see Figure 5. 3) We expect that, for small enough coupling parameters σ j , the spin-spin resonances are located close to the corresponding spin-orbit resonances for each ellipsoid, and would keep the same stability as for the uncoupled problem. However, we have found that, for larger σ j , the stability may change with respect to the uncoupled system, see Figure 6. 4) The coupled system is 5-dimensional, then, invariant tori, if there exist, would not confine solutions in determined regions (as in the uncoupled system), but Arnold diffusion is expected to take place.
A particular behavior occurs only when both bodies are identical (C 1 = C 2 = 0.5, λ = λ 1 = λ 2 , σ = σ 1 = σ 2 andq =q 1 =q 2 ), the so-called measure synchronization. This phenomenon was observed numerically in [HZ99] for an autonomous Hamiltonian system of two degrees of freedom (a pair of identical coupled oscillators): the system librates around a stable periodic solution in a very particular way described as follows for our system (see the phenomenon illustrated in Figure 7). Take a solution Ψ(t) librating around a stable spin-spin resonance of order (m : n, m : n). Consider the two families of points Ψ 1 (2πk) and Ψ 2 (2πk) and their corresponding annulus-like region in the plane Π. There are two possibilities: either both clouds of points are distributed in separated annulus-like regions or both regions coincide. That is to say, either the overlap is empty or there is a complete overlap. Moreover, if we start with a solution with separated regions, we can obtain the complete overlap by increasing the coupling parameter σ (keeping the same Ψ(0)). The phenomenon takes place suddenly for a specific σ = σ * when the outer boundary of the inner ring touches the inner boundary of the outer ring. At that moment, there is a concentration of density of points in the contact region. The phenomenon of measure synchronization disappears when the bodies are not equal. See in Figure 8 how the domains of both ellipsoids can overlap without merging into a single ring.

Linear stability of resonances
In this section we analyze the stability (in the linear approximation) of solutions associated to the resonances in different models, namely the spin-orbit problem (Section 5.1) and the spin-spin problem with spherical (Section 5.2) and non-spherical (Section 5.3) companion.
In all cases, we will only deal with balanced resonances (32), because they appear to be simpler and more relevant (see Section 4) than resonances of the general type (31). Actually, we can establish regions in the space of parameters where solutions associated to balanced resonances are unique or have some low multiplicity. In the case of the spinspin problem, we restrict ourselves to regions of uniqueness. The study of linear stability of such periodic solutions in the space of parameters complements the understanding of the dynamics that we presented in previous sections, especially Section 4. 5.1. Spin-orbit problem (V per = V 2 ). Consider the spin-orbit problem (23) with j = 1, that is, the motion of the ellipsoid E 1 . Let θ 1 = θ * (t) be a solution in a balanced m : 2 spin-orbit resonance, whose associated variational equation is y + λ 1 a r(t; e) 3 cos(2θ * (t) − 2f (t; e))y = 0, y ∈ R.
stability, whereas if | Tr(M )| > 2 we have hyperbolic instability. In the parabolic case, when | Tr(M )| = 2, the system is stable if the Jordan canonical form of M is 1 2 or −1 2 , otherwise the system is unstable. Actually, if the system is parabolic unstable, the instability of the linear system is linear in time, but hyperbolic instability is associated to an exponential divergence in time. In our case we want to distinguish regions of stability and instability in the (e, λ 1 )-plane for a given solution, which is continuous in (e, λ 1 ). For a given point (e, λ 1 ) and a given balanced m : 2 resonance, we want to know how many solutions there are of each type (34) or (35), and their linear stability. First, recall Remark 6: solutions of type π 2 satisfy conditions of type 0 for the equation (23) with j = 1, taking −λ 1 instead of λ 1 . Consequently, for each (e, λ 1 ), with e ∈ [0, 1) and λ 1 ∈ (−3, 3), we can obtain all the solutions corresponding to a balanced resonance by applying the shooting method: take a solution θ 1 (t) with initial conditions 4 θ 1 (π) = mπ/2,θ 1 (π) = γ ∈ R and let γ vary until the boundary condition θ 1 (0) = 0 is reached. Finally, we obtain the linear stability of the solution by computing Φ(t) such that Φ(π) = 1 2 . Actually, for this procedure, we can take any of the boundary conditions in eqs. (34) and (35), we just need one type to generate all solutions.
The results of this method for the main balanced spin-orbit resonances are shown in Figure 9. For these computations we used a Runge-Kutta Verner 8(9) integrator [Ver78], instead of a Taylor integrator [JZ05]; the reason is that, for some parameter values, the solutions are constant or polynomials and the Taylor method suffers in choosing a good step size. Thus, Figure 9 required around 3.5 days with 34 CPUs to be generated with a discretization mesh size of 2000 × 2000 × 2000 for (e, λ 1 , γ).
We can recognize the following characteristics 5 : 1) Each of the balanced resonances is represented in the (e, λ 1 , γ)-space by a continuous surface. In the case 1 : 1, the surface is made of two sheets connected only in one point (e, λ 1 ) = (0, 1). 2) The region of uniqueness in the (e, λ 1 )-plane is quite large. The multiplicity is generated by bifurcation of solutions: the surface folds generating multiple solutions (from one to three, as far as we see). Particularly, the bifurcation in the case 1 : 1 occurs at (e, λ 1 ) = (0, 1) producing a secondary sheet behind the main one. In general, the multiplicity takes place for some regions with |λ 1 | > 1. The resonances of order m : 2 with m > 2, have two characteristic folds, one with a V-like shape in solutions of type 0 for λ 1 > 1, and another small one in solutions of type π 2 for λ 1 ∼ −2 and very large eccentricities. 3) Instability is predominant in the diagrams, especially in solutions of type 0 for the resonance 1:2 and of type π 2 for the rest of the resonances. We see that the main regions of linear stability are continuation of stable solutions from e = 0, much of which are close to small |λ 1 |. Except for the resonance 1:2, the stability region for large eccentricities (e > 0.6) of the other resonances has a similar shape, characterized by a bifurcation with an interchange of stability. It is interesting to note that the folds producing the multiplicity are associated to some stable regions with peculiar shapes. In the resonance 1 : 1 we find two unstable regions bifurcating from the exact solution θ 1 (t) = t for e = 0: one at λ 1 = 1/4 = 0.25 (main sheet) and the other one at λ 1 = 9/4 = 2.25 (secondary sheet). Now let us consider both bodies. Since the system (23) is uncoupled, then, the multiplicity and stability associated to a spin-spin resonance is given by each of the separated problems. For example, take the (1 : 1, 3 : 2) balanced spin-spin resonance of type (0, π 2 ) for (e, λ 1 , λ 2 ) = (0.3, 0.1, 0.5), that is, the red points shown in Figure 9. It has associated a unique solution and it is unstable because it is so for θ 2 .
5.2. Spin-spin problem (V per = V 2 + V 4 ) with spherical companion. In this case we know that E 2 is in uniform rotation θ 2 (t) = θ 2 (0) +θ 2 (0)t, while the dynamics of θ 1 is given by (25), that depends on (e, λ 1 ,q 1 , σ 1 ). For this problem, we can proceed in the same way as for the spin-orbit problem of Section 5.1. On one hand, the variational equation associated to a solution in a spin-orbit resonance is a Hill's equation like (45), so the linear stability of the solution is characterized by the corresponding monodromy matrix. On the other hand, we can find all the solutions associated to a balanced spinorbit resonance using the shooting method for only one type of boundary conditions by including negative values of λ 1 , see Remark 6. Comparing Figures 9 and 10 we can see, for example, how the balanced 3 : 2 resonance is perturbed when we turn on the parameters (q 1 , σ 1 ): Figure 10. Diagrams of linear stability for the 3 : 2 balanced resonance of (25) for different values of the indicated parameters (q 1 , σ 1 ). Blue denotes instability and yellow denotes stability. Each of the columns have required around 15h using 24 CPUs and a mesh size of 512 × 512 × 512.
1) The effect of the new parameters is remarkable for large e and |λ 1 |. Especially whenq 1 and σ 1 are large. 2) At very large eccentricities, the surface has a complicated structure resulting in multiple solutions. The effect ofq is mainly to alter the stability for large e: solutions of type π 2 become always unstable, while stable regions of solutions of type 0 are more concentrated. The growth ofq also increases the multiplicity of solutions of type 0 and very large e. On the other hand, increasing σ 1 has a more dramatic effect on the complexity of the surface and also modifies the stability for large e. Actually, for some values of σ 1 , the existing V-shaped fold connects with the complex structure of large eccentricities in the upper part of the diagram. 5.3. Spin-spin problem (V per = V 2 +V 4 ) with non-spherical companion. We study the linear stability of the resonances of the full coupled spin-spin model in (24). In this case linear stability is determined by a more general theory and we will restrict ourselves to zones in the space of parameters where there is uniqueness.
Since we have two degrees of freedom, the first variation at a particular resonance is not a Hill's equation anymore, but a coupled system of second order. A Hill's equation is a particular case of linear Hamiltonian system with periodic coefficients, hereafter LPH systems.
If we define z = (θ 1 , θ 2 , p 1 , p 2 ) T , where p j = C jθj , the spin-spin problem in the Hamiltonian form (21) can be written asż where J l is the square matrix of order 2l given by with 1 l the unit matrix of order l. The non-autonomous Hamiltonian H K (t, z) is given in (19) for V per = V 2 + V 4 . Let z = z * (t) be a solution of (46) that is in a balanced spin-spin resonance of type (m 1 : 2, m 2 : 2). Then, the first variation at such solution has the formẏ = J 2 ∂ z,z H K (t, z * (t))y, y ∈ R 4 , (48) where ∂ z,z H K denotes the Hessian matrix of the Hamiltonian H K in the 4-dimensional variable z. The system (48) is an LPH system of period 2π. Assume that Φ(t) is a matrix solution of (48) with Φ(t 0 ) = 1 4 . The stability of (48) is determined by the Floquet multipliers, that are the eigenvalues of the monodromy matrix Φ(t 0 + 2π).
An LPH system has particular stability properties, see the general theory in [YS75,Eke90] and an application to the double synchronous spin-spin resonance in [Mis21]. For example, assume that ϕ ∈ C is a Floquet multiplier of an LPH system. Then, its inverse ϕ −1 , its complex conjugatesφ andφ −1 are also multipliers and have the same multiplicity as ϕ. That is, the Floquet multipliers have a symmetric distribution with respect to the real line and the unit circle of the complex plane. In consequence, a necessary condition for stability is that all multipliers have modulus 1. Moreover, if all multipliers have modulus 1 and they are all different, then the system is stable. When all the multipliers have modulus 1 and some of them coincide, the situation is not trivial and the stability depends on further algebraic properties of the multipliers (Krein's theory, [Kre50]). Then, unlike for Hill's equations, here we do not have a quantity like the trace of the monodromy matrix in order to characterize the boundary of stability/instability regions. Instead, we will use the following definition.
Definition 12. We will say that an LPH system is hyperbolic unstable if max k |ϕ k | > 1, where ϕ k , k = 1, 2, . . . , are all the Floquet multipliers of the system.
Assume that the solution z * (t) is continuous in some domain of the parameters of the model (e; C 1 , λ 1 , λ 2 , σ 1 ,q 1 ,q 2 ). The equation max k=1,...,4 |ϕ k | = 1 + ε, with ε > 0, defines a 1-parametric family of hyperbolic unstable manifolds in the space of parameters; then, the boundary of hyperbolic instability will be found if we take the limit of such manifolds as ε → 0.
On the other hand, it is possible to find all the solutions of a given type of a balanced spin-spin resonance using the shooting method as in the previous section. However, since the phase space and the space of parameters have large dimensions, our approach will be to obtain the solutions by continuation. Note that the solutions of (24) for λ j = 0 and any e ∈ [0, 1) are exactly given by θ j (t) = θ j (0) +θ j (0)t. Then, the unique solution of type (β 1 , β 2 ) of a balanced spin-spin resonance (m 1 : 2, m 2 : 2) is just θ j (t) = β j + m j 2 t. Such solution can be continued for |λ j | > 0. For small enough |λ j |, the systems in eqs. (23) to (25) can be regarded as different perturbations of the systemθ j = 0, then, Figures 9 and 10 give us a quite clear idea of a region of uniqueness of a given type associated to a balanced spin-spin resonance of (24). Moreover, such solution can be found by continuation in the space of parameters. Figure 11. Regions of (hyperbolic) linear instability of solutions associated to different spin-spin resonances of the problem (24) (case: equal bodies) for different values of the parameters (q, σ). Blue denotes instability and yellow denotes stability. Each plot has required around 10 minutes using 15 CPUs and a mesh size 750 × 750.
2) The effect of the coupling is different in each case: for the resonance (1 : 2, 3 : 2) of type (0, 0), we only see an additional thin unstable region for e < 0.1 and 0.25 < λ < 1. For the resonance (1 : 2, 3 : 2) of type ( π 2 , 0), we see that the region with small e becomes unstable, whereas for large e the stability is somewhat altered. Finally, without coupling, the resonance (1 : 2, 3 : 2) of type (0, 0) is unstable for almost all the points, but the coupling introduces the stability for small e. Note that this in agreement with Figure 6.

Comparison between the full and the Keplerian models
In this section, we provide some results on the comparison between the full and Keplerian problems (Section 6.1), and we give some numerical results on the interaction between the spin and orbital motions (Section 6.2).
6.1. Hamiltonian approach. It will be useful to write the dynamical equations of the Hamiltonian of the full model (4) in the compact forṁ where z = (r, f, θ 1 , θ 2 , p r , p f , p 1 , p 2 ) T and J 4 is defined by (47). Let us now formulate the Keplerian models (spin-orbit and spin-spin, see Section 2.3) as perturbations of the full model, so we can compare both families of models. Consider a function ζ(t) that satisfies the equatioṅ where the vector function is responsible of subtracting the perturbative part of the potential (see equation (7)) only in the equations forṗ r andṗ f . Thus, equation (50) represents the Keplerian models including the orbital part: the spin-orbit model corresponds to (50) with V per = V 2 and the spin-spin model with V per = V 2 + V 4 . The corresponding equations of motion are (11) and (12), so we can write the solution in the form ζ(t) = (r(t; a, e), f (t; e), θ 1 (t), θ 2 (t), p r (t; a, e), p f (a, e), p 1 (t), p 2 (t)), where a and e are the semimajor axis and the eccentricity of the Keplerian orbit, see (17). Now it is clear that the Keplerian model (50) is not Hamiltonian, even though we can split it into one autonomous Hamiltonian system (orbital part with V = V 0 ) and another non-autonomous one (spin part with V = V 0 + V per ). Moreover, unlike for the full model, in a Keplerian model neither the total angular momentum P f = p f + p 1 + p 2 is conserved 6 , since we can compute thaṫ ∂ θ j V per (r(t; a, e), f (t; e), θ 1 (t), θ 2 (t)).
We want to know if ζ(t) is a good approximation for a solution of the full problem (49) with the same initial conditions. Consider the solution z = z(t) = ζ(t) − δz(t) of (49) such that z(0) = ζ(0). Then, expanding the equation (49) up to first order in δz, we obtain that the function δz = δz(t) satisfies the equation with initial condition δz(0) = 0. In equation (51), ∂ 2 z,z H(z) is the Hessian matrix associated to the Hamiltonian in (4). If the system (51) is stable, then the norm ||ζ(t) − z(t)|| is bounded. Additionally, it is easy to see that the system (51) is Lyapunov stable if and only if the trivial solution of the homogeneous parṫ is Lyapunov stable. A general form of ζ(t) is unknown because, although the orbital part is given by the classical Kepler problem, the spin part is given by a non-autonomous periodic system of nonlinear equations. However, ζ(t) is a periodic solution at spin-spin resonances, with period 2π for balanced resonances. In such cases, equation (52) is an LPH system, some of whose properties were mentioned in Section 5.3. At this point, we wonder if it is possible for the periodic system (52) to be stable. Let us point out an argument supporting a negative answer, even for stable spin-spin resonances. It is known that the periodic solutions of the planar Kepler problem are not linearly stable. This can be easily checked using Poincaré variables (action-angle variables): we obtain a positive eigenvalue for a linear system of constant coefficients. See further related discussions in [BO16,Sch72]. We can see that 1 is the only associated Floquet multiplier and has multiplicity four; then, the instability of the periodic solutions of the Kepler problem is not hyperbolic, according to Definition 12, but rather of parabolic kind. From this discussion, we expect that ζ(t) and z(t) are divergent functions, but we want to know if such divergence is exponential in time, that is, when (52) is hyperbolic unstable.
Suppose that the function ζ(t) is continuous on a domain of the space of parameters (a, e; C 1 , λ 1 , λ 2 , σ 1 ,q 1 ,q 2 ). Note that we add a to the parameters of the spin-spin model because it varies with the orbital initial conditions. On the other hand, the Hamiltonian H depends on the parameters of the full model, we take the independent set (µ, C 1 , d 1 , d 2 , q 1 , q 2 ), where the value of µ informs us about the disparity in the masses of the bodies because µ = M 1 M 2 = M 1 (1 − M 1 ). We can obtain (µ, C 1 , d 1 , d 2 , q 1 , q 2 ) from (a, e; C 1 , λ 1 , λ 2 , σ 1 ,q 1 ,q 2 ) as follows: First, from (22), we see that, assuming σ 1 > 0, 6 Instead, the Keplerian assumption results in the conservation of the orbital angular momentum p f (t) = µr(t; a, e) 2ḟ (t; e) = µa 2 √ 1 − e 2 . from this value we can obtain M 1 because µ = M 1 (1 − M 1 ) with M 1 ∈ (0, 1). This is enough to write M 2 = 1 − M 1 , C 2 = 1 − C 1 and, from the definitions (22), it holds that From these relations, it is obvious that, in order to compare the full and the Keplerian models, the spin-orbit versions are not enough, and we need to take V per = V 2 + V 4 , say, the spin-spin models. Additionally, the initial conditions associated to ζ(t) are given by (18), for the orbital part, and the values θ j (0),θ j (0) are such that the spin part satisfies the boundary conditions for the spin-spin problem in a spin-spin resonance of a certain type as in (39). Recall also that, in our units, the Keplerian orbital period is T = 2π and G = a 3 .
As in Section 5.3, we can find the region (in the parameters space) of hyperbolic instability and its boundary, that is given by the limit as ε → 0 of the manifolds that satisfy max k=1,...,8 |ϕ k | = 1 + ε, where ϕ k are the Floquet multipliers of (52). It is important that we focus on a region of the parameters where the corresponding solution associated to a spin-spin resonance is linearly stable, so we can see what is the effect of putting the orbital and the spin parts altogether. Actually, not even the spin part of (52) is completely equivalent to (48), because all the derivatives of V per (r, f, θ 1 , θ 2 ), evaluated at ζ(t), appear in (52) for both orbital and spin variables. Figure 12. Regions of (hyperbolic) linear instability of (1 : 2, 3 : 2) of type (π/2, 0) of the equation (52) for equal bodies and different values of the parameters. Blue denotes instability and yellow denotes stability. Each plot has required around 10 minutes using 15 CPUs and a mesh size of 750 × 750. Figure 12 shows the regions of hyperbolic instability of (52) in the case of equal bodies for a particular spin-spin resonance with different values of the parameters. We observe the following: (1) We can compare Figure 12 with Figure 11 for similar values ofq and σ. In the case of equal bodies (µ = 0.25 and C 1 = 0.5), the value of a is determined by (53): smaller σ also implies further bodies. Then, the plots with σ = 0.0001 in Figure 12 and be compared approximately with those with σ = 0 in Figure 11. (2) From the comparison, we see that the plots in Figure 12 with Figure 11 witĥ q = 0 are similar for small eccentricities, but, for large e, the unstable regions of plots in Figure 12 are larger. On the other hand, the plots in Figure 12 witĥ q = 0.01 include unstable stripes invading the stable regions for small e (they grow with σ), whereas for large e the diagrams become totally unstable.
(3) As we increase σ (closer bodies), the plots become more and more unstable.
We conclude that the spin-spin model should be a reliable simplification of the full model in the regions of the parameters with small e, large a and close to regions with stable spin-spin resonances. We remark also that instability of spin-spin resonances not necessarily implies hyperbolic instability of (48): there are some regions for which the Floquet multipliers of the spin-spin resonance are not too far from the unit circle, so that when we consider the corresponding system (48), all its Floquet multipliers belong to the unit circle. It is also remarkable that, although the comparison makes sense only between the full and Keplerian spin-spin models, we obtain almost identical plots as in Figure 12 independently if we use V per = V 2 or V per = V 2 + V 4 in the Hamiltonian in (52).
6.2. Quantitative numerical approach. In this section, we want to investigate, from a numerical point of view, two questions that the Hamiltonian approach left unanswered: (Q1) Can we quantify the influence of the spin motion on the orbital one? (Q2) Is it possible that the bodies end up colliding, even though (52) is not hyperbolic unstable?
In Section 6.1, we saw that the solution of a Keplerian model ζ(t) corresponding to a spin-spin resonance should diverge from a solution of the full model z(t) with identical initial conditions. In the region of parameters of hyperbolic instability the divergence should be exponential; however, in the rest of the parametric space, we want to quantify how different both motions are. So, let us take the point (a, e; C 1 , λ 1 , λ 2 , σ 1 ,q 1 ,q 2 ) and compute the functions ζ(t) and z(t) in order to compare them.
The orbital motion of ζ(t) is characterized unambiguously by the set of Keplerian elements (a, e, ω), where ω = 0 is the argument of the periapsis, together with t, the mean anomaly. On the other hand, let the orbital part of the solution z(t) of the full model be given by (r F (t), f F (t)). Then, we can transform the orbital position r F (t) = r F (t) exp(if F (t)) to the osculating Keplerian elements (a F (t), e F (t), ω F (t)) of the twobody problem using the following expressions. From the geometrical identity Now let us define the orbital angular momentum per unit mass h F = r F ∧ṙ F , where ∧ is the vector product, and the eccentricity vector whose modulus is given by Whenever e F = 0, we can define ω F ∈ [−π, π) as the polar angle of e F . The mean anomaly associated to the full model can be defined too, but we will not use it in our study. The balanced spin-spin resonance of order (m 1 : 2, m 2 : 2) in the function ζ(t) is characterized by the modified resonant angles ψ m j :2 j (t) = m j f (t, e) − 2θ j (t). Seemingly, the spin part of the solution z(t) of the full model is given by (θ 1,F (t), θ 2,F (t)), so, in order to compare with ζ(t), let us define ψ where δ a is the relative deviation in semi-major axis, while δ e and δ res,j are the absolute deviations in eccentricity and resonant angles, respectively. Finally, recall from (1) that we have an expression for a j in terms of the parameters of the model. Then, for our purpose, we will say that there is a collision if r F (t) ≤ a 1 + a 2 . Now we are in a position to compare the solutions of the full model with respect with those of the Keplerian models.
In Figure 13 we see, for different values of the parameters and the case of identical bodies, the comparison between z(t) and ζ(t) in a resonance (1 : 1, 3 : 2) of type (0, 0). Particularly, we see the evolution of the δ-functions in (54) and some Kepler elements of z(t) in 100 revolutions. In addition, Table 2 informs us about the corresponding Floquet multipliers of both z(t) and ζ(t). We summarize the following observations: (1) We chose two cases for the comparison. In the first case, e = 0 and σ = 10 −3 , so the bodies are relatively close to each other (a ≈ 39) in circular orbits, whereas in the second case, e = 0.1 and σ = 10 −7 , so the bodies are much further (a ≈ 3900) in moderately elliptic orbits. For both cases, we took values of parameters whose corresponding resonance is not hyperbolic unstable, say, λ = 0.05 andq = 0.01. (2) From Figure 13, we see that the first case has a more regular behavior than the second one. In the first case, the orbit of the full model oscillates regularly very close to the circular orbit of the Keplerian model, with a precession that increases uniformly in average. Actually, in a shorter time scale, the eccentricity vector describes a circles passing through the origin. On the other hand, the orbit of the second case is greatly irregular, with variations of size up to ∼ 25%. Even its precession changes direction at some moment and the eccentricity vector swings chaotically in short time scales. (3) The variation of the angles is quite regular in both cases. Except for the resonant angle of the resonance 1:1 of the first case, which oscillates regularly with small amplitude, the rest of the resonant angles increase/decrease more or less constantly, with larger variations in the second case. (4) The Floquet multipliers in Table 2 gives us more information about both cases.
We remark that we display the multipliers of ζ(t) corresponding to spin and orbit separately, say, the first two rows are the four multipliers of the spin-spin model (24), whereas the third one shows the four coincident multipliers of the Kepler problem. We confirm that the spin of ζ(t) is elliptic stable in the two cases. On the other hand, the multipliers of z(t) are computed with (49). We find out that, although the first case is more regular, the solution is actually hyperbolic unstable, whereas the opposite occurs for the second case: it is irregular but not hyperbolic. This is consistent with Figure 12. e = 0, λ = 0.05,q = 0.01, σ = 10 −3 e = 0.1, λ = 0.05,q = 0.01, σ = 10 −7  Table 2. Modulus and argument of the Floquet multipliers of the functions z(t) and ζ(t) compared in Figure 13.
We conclude that, even when there is hyperbolic instability, the solution of full system remains quite close to the solution of the Keplerian system with circular orbit. However, even with no hyperbolic instability, for e = 0, the behavior of both solutions is remarkably different. This fact should be attenuated with very small eccentricities.
We can make a step further in the comparison when e = 0. Figure 14 shows a quantitative comparison in the orbits of the full and Keplerian models. We vary the parameters σ and λ keepingq = 0 constant. For large enough σ (small a), there is a region of parameters leading to collision. The rest of the diagrams show different orders of magnitude of δ-functions corresponding to a and e. We check that, for very small λ, the orbit of the full model is very close to the Keplerian one.

Conclusions
This research is, primarily, a careful numerical study of the spin-spin model, presented as such in [Mis21]. For this purpose, we have considered several parallel models describing the planar gravitational dynamics of two ellipsoids under the Assumptions 1 to 3. We distinguish between full and Keplerian models, and also, between spin-orbit and spin-spin models. The known dynamics of the Keplerian spin-orbit problem was used as reference to develop our results.
In first place, we realize that the symmetries of the ellipsoids lead to certain symmetries in the equations of the Keplerian models, and so, to symmetric periodic solutions that structure the overall dynamics. We complete this general view by providing conditions for existence of quasi-periodic solutions. Between the periodic solutions, we focus on a special class that comprises the simplest solutions from the physical interpretation, the balanced spin-orbit and spin-spin resonances. We study the high-dimensional phase space of the spin-spin problem by means of projections of Poincaré maps, taking those of the spin-orbit problem as reference.
First, we see that when one body is spherical, the dynamics of the spin-orbit problem is reproduced with small variations. For two ellipsoids, the projections of the Poincaré maps show structures similar to those of the spin-orbit ones. There is a particular behavior when both bodies are identical, the so-called measure synchronization.
After that, we study existence, multiplicity and linear stability of the solutions in balanced resonances as we vary the parameters of the problem. We focus on the qualitative changes in the stability diagrams produced by the variation of each parameter. We use the case of identical bodies to illustrate our observations most of the time.
In the last part of the paper, we compare the solutions of the full and the Keplerian models by means of rewriting the equations in a Hamiltonian-like form. We produce stability diagrams similar to the previous ones to show such a comparison. Representing different aspects with similar diagrams allows one a much more global understanding of a problem with so many variables and parameters. We end up by comparing the evolution of particular solutions for certain parameters. The orbit of the full problem is characterized by Kepler's elements to facilitate the comparison. From the parameters, we observe that the comparison only makes sense for the spin-spin models. Here we find out that, beside the case of circular orbits, for which full and Keplerian models show a close behavior, for eccentric orbits, the trajectories are very dissimilar. We finish the comparison in the case of circular orbits by showing that the space of parameters is divided into regions leading to collision or regions with orbits differing specific orders of magnitude.
The topics and results presented in this work are just a taster of some features of the rich spin-spin dynamics. We also include a comparison with the full models, which is not widely studied in the literature; however, this is crucial, because it shows us the limits of validity of the models for applications. We believe that here is still a lot to explore in this model, which can be used to investigate problems like Arnold's diffusion or the existence of normally hyperbolic manifolds.
Appendix A. Coefficients of the expansion of the potential The full expansion of the potential energy of the Full Two-Body Problem was derived in [Bou17]. Later, [Mis21] detailed that expression to the case of planar motion of two ellipsoids, which is given by 2l 2 ,2m 2 cos(2m 1 (θ 2 − f ) + 2m 2 (θ 2 − f )) , To perform the integration of (12), it is convenient to adopt the eccentric anomaly u as independent variable, according to the following procedure.
Let us consider the change of variables given by the Kepler's equation x j (u) = θ j (u − e sin u) , j = 1, 2 .