Influence of Energy Dissipation on Spin Pole Precession Trajectories for Synchronous Rotators

We present an analysis of the influence of tidal energy dissipation on the precessional motion of the spin pole of a synchronously rotating body, whose orbit pole is also precessing. In the absence of dissipation, this problem has been frequently analyzed in a Hamiltonian context. The Cassini states, or fixed points of the motion of the spin pole, correspond to tangent intersections of the Hamiltonian function with the unit sphere. When dissipation is added to the dynamics, the locations of the fixed points change. We present a new geometrical characterization of the nullclines of the dissipative flow, or motion of the spin pole, and use it to find the locations of the dissipative fixed points. Using the Poincaré−Hopf index theorem, we show that in the case where there are four fixed points, one is a saddle, and depending on the sign of the dissipation rate parameter, the other three comprise either two sources and one sink or one source and two sinks. We also examine the domains of attraction of the sinks and present an analysis of their stability to initial conditions and model parameters.


Introduction
We examine the role of a dissipative term in the equations of motion for the spin pole precession of a rotating triaxial body, which is moving along an orbit that is uniformly precessing about a fixed pole.
In the usual, nondissipative treatment, the motion conserves a scalar quantity. Trajectories lie on the intersection of the unit sphere and a parabolic cylinder. Depending on the radius of curvature of the parabola, at its vertex, there can be two, three, or four fixed points. These are the Cassini states of Colombo (1966) and Peale (1969).
When dissipation is added, there is no longer a conserved scalar, but we can find fixed points as the intersections of the nullclines, which are curves on the surface of the unit sphere along which the latitudinal and longitudinal components of velocity vanish. As in the nondissipative case, there can be two, three, or four fixed points. When there are two fixed points, one is a source and the other is a sink. When there are four fixed points, there is a saddle and either two sources and one sink or one source and two sinks. The size and shapes of the domains of attraction of the sinks depend on the parameter values.
In planetary or satellite systems, with finite values of either orbital eccentricity or spin pole obliquity, the tides will be time dependent, and the stress and strain associated with those tides will lead to energy dissipation. This dissipation will tend to drive the spin and orbit parameters of the system toward one of a family of end states, in which the dissipation is generally reduced. The primary objective of this paper is to examine how obliquity tides influence spin pole evolution.
In isolated binary systems, eccentricity tides (Peale 1999;Wisdom 2008) will act to remove energy from the system as a whole and transfer angular momentum between the spins of the bodies and the orbit. The usual outcome is synchronous rotation for the smaller body and damping of the orbital eccentricity (Counselman 1973). In the double synchronous end state, such as Pluto-Charon (Farinella et al. 1979;Taylor & Margot 2011), tidal dissipation ceases. In cases like the Galilean satellite system, where resonant interactions between the orbits can pump up the eccentricities, the end state will have synchronous rotation but may still have finite orbital eccentricities (Peale 1976;Brown 1977;Henrard 1983).
Obliquity tides tend to damp the obliquity and bring the spin pole toward alignment with the orbit pole. In this case, the circumstance that can prevent eventual vanishing of obliquity, in the dynamical end state, is precession of the orbit pole. This behavior has been studied in a number of different contexts, but most of them have only considered dissipation in a rather perfunctory manner, assuming that it will drive the system into an end state, but that the configuration of that end state does not depend on dissipation. We present a simple model of spin pole motion, for a synchronously rotating triaxial body, whose spin pole is both damping toward the orbit pole and precessing about the orbit pole, while the orbit pole, in turn, is precessing about an inertially fixed direction.
To simplify the analysis, we assume that the binary orbit that our triaxial body shares with a point-mass partner is circular. If the orbit were actually eccentric, it would make small changes to the spin pole precession and tidal damping rate parameters, but other than that, the theory should remain unchanged.
The functional form that we assume for the dissipative influence on the spin pole motion is admittedly somewhat arbitrary. Other functional forms could be proposed (Peale 1974;Ward 1975;Fabrycky et al. 2007). However, our model has two key features, which it will share with any other plausible model of this process. First is that it moves the unit vector aligned with the spin pole in a way that keeps it a unit vector. Thus, the spin pole unit vector must move in a direction orthogonal to itself. Second is that the motion of the spin pole should decrease the potential energy associated with orientation of the triaxial body.
Thus, any other plausible model will differ from ours only in higher-order terms, involving the rate parameter, the obliquity, or both. Our specific choice for a dissipative contribution to spin pole motion was motivated by finding the simplest model with these two attributes.

Background
The problem of interest in this study is the influence of a dissipative torque on the spin precession of a rigid triaxial body, which is moving along an orbit that is precessing at uniform angular rate about a fixed axis. That problem is clearly related to the so-called "Colombo top" (Colombo 1966;Henrard & Murigande 1987;Haponiak et al. 2020), which is very nearly the same dynamical problem but treated in a Hamiltonian setting, thus ignoring the influence of dissipation.
The early history of that problem started with a description of the motion of the Moon, by Cassini (1693), and was later summarized by Tisserand (1891), in three simple statements: (1) the rotation and orbital periods are equal; (2) the rotation axis has a constant inclination to the ecliptic; and (3) the three axes, related to spin, orbit, and precession, are coplanar. Colombo (1966) showed that the first law is quite independent of the second and third. We retain a formulation for the spin pole precession rate that is applicable in either resonant or nonresonant cases. There is a significant, and rapidly growing, literature on the rotational dynamics of planets and satellites. An important subset of that literature concerns the fixed points of spin precession, for bodies whose orbits are also precessing. It was noted by G. D. Cassini, in 1693, that the Moon is in a fixed spin-orbit precessional state and has adjusted its obliquity, or angular separation between spin and orbit poles, so that the period of spin precession matches the period of orbit precession, even though their angular rates of motion are different. In addition, the lunar phases are also adjusted so that the Moon's spin pole sˆremains coplanar with the orbit pole n and invariable pole kˆ, even as n precesses about kˆand sp recesses about n. For a nice discussion of the relevant lunar dynamics and some history, see Eckhardt (1981).
Cassini's description was purely kinematic. He described what happens but gave little insight into the causal mechanisms. Recent progress on understanding the dynamics of this situation started with Colombo (1966), who showed that the fixed point occupied by the Moon is understandable in terms of the geometry of the Hamiltonian function for the system. The Hamiltonian has the geometric form of a parabolic cylinder, and the spin pole must lie on the unit sphere. Any spin pole trajectory thus lies on the intersection of those two surfaces. If the Hamiltonian has a tangent intersection with the unit sphere, the trajectory collapses to a fixed point. This configuration allows the spin pole precession rate driven by the gravitational torque to match the precession rate of the orbit. It yields a configuration that minimizes energy dissipation. Peale (1969) extended the analysis to triaxial bodies and eccentric orbits and showed that there are either two or four such fixed states, now known as Cassini states, depending on values of the relevant parameters.
When there are four fixed states, three of them are centers, and the other is a saddle. Two homoclinic loops emanate from the saddle. These loops divide the surface of the unit sphere into three disjoint regions, each containing a center. As the system parameters change, one of the loops shrinks, and eventually the saddle merges with one of the centers. Further change leads to a state with only two centers. Henrard & Murigande (1987) showed how to compute the areas of the three enclosed regions and gave an explicit criterion for the transition state, where the number of fixed points changes. Ward (1975) applied this theory to past motion of the Moon and showed that it experienced a dramatic transition in obliquity when the Cassini state it had initially occupied ceased to exist and quickly transitioned to the state it now occupies.
A somewhat curious aspect of past studies of these Cassini states is that most analyses have ignored the role of dissipation. This has the advantage of making the Hamiltonian analysis relevant. However, without dissipation, there is no reason to expect a body to actually occupy a Cassini state. In a nondissipative system, the spin pole will cycle around one of the centers, with an amplitude set by the initial conditions. A notable exception is the work of Gladman et al. (1996), which examined some aspects of capture into synchronous rotation and damped precession. However, even there, the results were mainly numerical integrations of the equation of motion. Fabrycky et al. (2007) also discuss dissipation, but mainly in the context of how much dissipation can be produced by obliquity tides. In addition, there have been several recent analyses in which the influence of dissipation on the orientation of the spin pole, compared to that expected in a classical Cassini state, has been considered, including Baland et al. (2017) and Organowski & Dumberry (2020).
The objective of the present study is to examine aspects of the precessional motion of the spin pole of a synchronously rotating planet or satellite, when obliquity tide energy dissipation is included. We will develop the relevant equations of motion, find formulae for fixed points, and numerically integrate solutions. In attempting to describe the geometry of fixed states when dissipation is included, we have used nullclines, which are the curves on the unit sphere along which either the latitudinal or longitudinal velocity vanishes. The intersections of these two curves yield fixed points. We have found simple geometric characterizations of the nullclines, which apply to both dissipative and nondissipative cases.
In the case without dissipation, the spin pole trajectories are set by system rate parameters and the initial position of the spin pole. The trajectories are, as noted above, simple closed curves on the surface of the unit sphere, along which the sphere intersects a parabolic cylinder. The addition of finite dissipation changes the character of the flow, or motion of the spin pole, considerably. It changes the positions and character of the fixed points. In the nondissipative case, all of the fixed points lie in the plane defined by the orbit pole and invariable pole. Dissipation moves them out of that plane. Instead of having fixed points that are either centers or saddles, dissipation changes the character of the fixed points, with centers being transformed to spiral sinks or spiral sources.
The remainder of the paper is divided into seven additional sections. Section 3 examines the nondissipative equations of motion and recalls much of the relevant theory for that case. Section 4 adds dissipation and examines how the nullclines and fixed points change. Section 5 examines vorticity and divergence of the flow, or motion of the spin pole on the unit sphere. Section 6 considers stability of the fixed points and classifies them based on the eigenvalues of the relevant Jacobian matrix. Section 7 uses the Poincare−Hopf index theorem to consider the number of possible fixed points and examines the basins of attraction for fixed points that are attractors. Section 8 applies our theory to the Moon. Section 9 gives a summary.

Conservative Spin Pole Precession
As a prelude to discussion of dissipative motion, we first examine several aspects of the precessional motion of a synchronous triaxial body moving under the influence of rigidbody gravitational torques. We first write the corresponding differential equations for motion of the spin pole, in several different coordinate frames. Next, we examine a scalar quantity that is conserved along the resulting precessional trajectories. Finally, we examine fixed points of the motion.

Equations of Motion
We first suppose that the orbit pole, specified by unit vector n, is precessing, at fixed angular rate γ and fixed inclination i, about an inertially fixed direction, or invariable pole kˆ, so that Next, we suppose that the spin pole, denoted by unit vector sˆ, is precessing about the instantaneous orientation of the orbit pole n with rate parameters α and β, so that ds dt n s n s • .
In the case of a triaxial synchronous rotator, in a circular orbit, with principal moments of inertia A < B < C, these rate parameters are given by (Varadi et al. 2005;Newman 2014) where n is the orbital mean motion.
In the case of a nonsynchronous rotation, we simply average the potential associated with the triaxial body, over a rotation period. That has the effect of making the body behave, in terms of spin pole precession, as an oblate spheroid, with A = B, and thus β = 0. In addition, the prefactor changes from n to n 2 /s, where s is the spin rate. In the following, we retain the full form, including a possibly nonzero value for β.
The obliquity ε, or angular separation of spin and orbit poles, is given by n s cos . 5 e = [ ]ˆ·ˆ( ) It will be convenient to describe the motion of the spin pole, as seen in a uniformly rotating reference frame that is fixed in the precessing orbit plane. In the orbit-fixed frame, the spin pole motion, which is described by Equation (3) in an inertially fixed frame, is governed by We begin our analysis by defining three orthogonal unit vectors, in the orbit-fixed frame. One of them is just the orbit normal n. The other two lie in the orbit plane. The first of them is m k n k n k n k n . 7 = --ˆ(ˆ·ˆ)|ˆ(ˆ·ˆ)ˆ( ) It lies in the plane defined by kˆand n and is orthogonal to n. The third basis vector is just If we now project the spin pole motion, as given in Equation (6), onto this orthogonal triad, we obtain ds dt a n s b n s c m s where the new coefficients (a, b, c) are related to the old ones (α, β, γ) via This transformation converts four parameters {α, β, γ, i} into three parameters {a, b, c} and is thus clearly not invertible. As will be seen below, the precessional behavior is governed by these three parameters, and it does not matter which set of four input parameters were used to obtain them. We will use the term "inertial precession parameter" for c. It will be convenient to use a spherical polar coordinate system {r, θ, f}, with radius r, measured from the intersection point of the unit vectors n and sˆ; colatitude θ, measured from the direction of the orbit normal n; and longitude f, measured in the orbit plane, from the direction of the unit vector lˆ.
In this l n m , , {ˆˆˆ} frame, the Cartesian coordinates {x, y, z} for the spin pole position sˆare s x y z , , sin cos , sin sin , cos .
The equations of motion in the Cartesian frame can thus be written as The corresponding velocity field in the spherical coordinates frame is given by

Conserved Scalar
We note that solutions to the differential equations (12) conserve a scalar quantity, which is given by Similarly, in spherical form, Equations (13) yield a conserved scalar This is the Hamiltonian discussed by Colombo (1966), Peale (1969), and Henrard & Murigande (1987). In classical treatments of Hamiltonian mechanics, the Hamiltonian function is the total energy of the system, represented by the sum of kinetic and potential energies (Goldstein et al. 2002). If the quantity in either Equation (15) or Equation (16) is multiplied by n C, where n is orbital mean motion and C is the polar moment of inertia of the rotating body, it has the expected dimensions of energy. We note that Henrard & Murigande (1987), in their treatment of the conservative system, use a two-parameter version of this Hamiltonian, written as where their parameters are related to ours via

Fixed Points
A significant advantage of the Hamiltonian approach is that conservation of this scalar quantity along any valid trajectory of precessional motion of the spin pole allows much to be inferred about the dynamics of the system, in a purely geometrical fashion. In particular, the set of points on which this function has a given constant value, K x y z K , , , defines a surface in the Cartesian coordinate frame {x, y, z}. In particular, it is a parabola in the {y, z} plane and has no dependence on x, and thus it is a parabolic cylinder. Any trajectory of the spin pole, generated by the equations of motion, will move along curves where the unit sphere intersects this parabolic cylinder.
This parabola can be written succinctly as The axis of the parabola is at The radius of curvature, at the vertex, is Abs[w]. We can use the control parameters {a, b, c} to generate a family of parabolas, sharing the same axis and having the same radius of curvature (c/a) at the vertex. As the numerical value of the scalar K p is changed, the parabola moves along the central axis. Intersections of this parabolic cylinder with the unit sphere define allowed trajectories of the spin pole. Among this family of allowed trajectories, those whose associated parabolas yield tangent intersections with the unit sphere play a special role. It will be important to distinguish between internal and external tangent intersections. We will call an intersection of the parabola with the unit sphere "internal" if points on the parabola, near the intersection, lie inside the sphere. That will occur if the radius of curvature of the parabola, at the intersection point, is less than 1. Likewise, at an "external" tangent intersection, the nearby points on the parabola are outside of the unit sphere. An osculating intersection is one for which the radius of curvature of the parabola, at the intersection point, exactly matches that of the sphere.
Though there can be two, three, or four tangent intersections of the parabola with the sphere, at most one of them is internal, having radius of curvature, at the tangent intersection point, less than 1. The radius of curvature of that parabola, at that point, is given by (Gray et al. 2006 When the value of R is (less than (internal), equal to (osculating), greater than (external)) the value 1, there are (four, three, two) possible tangent intersections. This result, for the transition in number of fixed points, in terms of parameters equivalent to our {a, b, c} was given in Henrard & Murigande (1987), but their derivation was less geometrical in style.
If the tangent intersection is external, the allowed trajectory will collapse to a point. For internal tangent intersections, the trajectory will divide the surface of the unit sphere into three disjoint regions, each of which will contain an external tangent intersection, or fixed point. Figure 1 illustrates such a case. Also shown in the figure are nullclines of the flow, which are curves along which the θcomponent or f-component of velocity vanishes. At points where the θ and f nullclines intersect, the velocity vector vanishes: these are the fixed points. We will have more to say about nullclines subsequently.
Note that, in this and all subsequent figures, we use latitude rather than colatitude θ.

Precession with Tidal Dissipation
When we add dissipation to the equations of motion, the Hamiltonian is no longer conserved, and we need to find another approach to studying the behavior of a precessing spin pole. The approach we will adopt is to examine the "nullclines" of the equations of motion. In general, nullclines are surfaces in the phase space on which one of the spatial variables has zero time derivative. Intersections of nullcline surfaces yield fixed points. For a general treatment of nullclines in dynamical system analysis, see Strogatz (1994).

Augmented Equations of Motion
For a synchronous rotator, with finite values for orbital eccentricity e and obliquity ε, the rate of tidal energy dissipation can be written as (Wisdom 2004(Wisdom , 2008 where k is the tidal Love number (Munk & MacDonald 1960), Q is the loss factor, n is orbital mean motion, R is the mean radius of the body, and G is the gravitational constant. Finite eccentricity causes the subprimary point to move in longitude, while the finite obliquity causes it to move in latitude. The primary effects of the eccentricity and obliquity tides are to damp the eccentricity and obliquity, respectively. We now seek a formulation for the influence of energy dissipation, associated with obliquity tides, on the motion of the spin pole. This has two main aspects. We expect that the dissipation will, acting alone, bring the spin pole closer to alignment with the orbit pole. Also, the rate of spin pole motion should have a plausible connection to the rate of tidal energy dissipation.
An initial ansatz for the functional form of the motion is where δ is a rate parameter and q is a dimensionless scalar factor. For positive values of q, it does have the desired property of moving the spin pole toward the orbit pole, and it does so at a rate related to the difference in orientation of the two poles. However, unless q n s • , 25 =ˆˆ( ) it also changes the length of the vector sˆ. To see that, consider In order to make that rate vanish, we adopt the functional form We now seek an estimate of the numerical value of the dissipative motion parameter δ, using an argument based on energy dissipation rates. Using the Hamiltonian given by Ward (1975), for the rotation of the Moon, the potential energy is seen to depend on the principal moments of inertia (A < B < C) and the orientation of the spin pole. The lowest value is achieved when the spin and orbit poles are aligned. The potential energy, at finite obliquity ε, can be written as The partial derivative of this potential energy, with respect to obliquity, is then The rate of energy loss, associated with obliquity tides, is given by Equation (23)  which, for small obliquity values, then yields This allows a connection between our spin pole damping rate parameter δ and the commonly used parameterization of tidal dissipation in terms of k and Q, as in Wisdom (2004Wisdom ( , 2008 and Winn & Holman (2005).
We restrict the parameter to be positive so that it tends to move the spin pole into alignment with the orbit pole and thus decreases the obliquity contribution to potential energy of the configuration.
We are not aware of a previous presentation of this particular formulation for the influence of tidal dissipation on orientation of the spin pole. Since this is something of uncharted territory, we provide the simplest mechanism that achieves two key features. The dissipation-driven motion of the spin pole needs to reduce the angular separation between spin and orbit poles and needs to preserve sˆas a unit vector. The form given in Equation (27) has both attributes. Any other valid form for the rate of change in spin pole orientation, due to obliquity tides, will surely be similar to this.
A potential problem with this approach is that it assumes 100% efficiency in the conversion of dissipated energy into reduction of the potential energy associated with spin pole orientation. We expect a relatively low efficiency, with most of the energy instead going into heating of the body. In fact, when we apply this formulation to the Moon, we will see that the mechanical efficiency is quite low.
We also note that the influence of eccentricity tides on the spin state of a planetary body will tend to drive the spin rate toward resonance with the orbital mean motions. It can be either synchronous, like the Moon, or a low-order commensurability, like Mercury. The focus of the present study is rather different: we seek to address the question of how obliquity tides drive the spin pole orientation. For a treatment of the approach to synchronous rotation, as driven by eccentricity tides, see Gladman et al. (1996).
When this dissipative term is added to the vector form of the equations of motion (9) and we let the unit vectors have projections onto the Cartesian coordinates s x y z n m , , 0, 0, 1 0, 1, 0 , 33 the resulting set of scalar differential equations is Thus, the dissipative term only changes the time derivative in the θ direction. However, as will emerge below, the locations of the fixed points can be changed in both θ and f. Note that if a = b = c = 0, so that the dissipation term is the only contribution to moving the spin pole, we obtain the solution so that the parameter δ is just the exponential damping rate of the tangent of half the obliquity.

Fixed Points
There is considerable interest in finding fixed points for the flow defined by this model of dissipative precession. The next several sections of this paper will focus on locating the fixed points and addressing their stability. In principle, finding the fixed points is very simple. It is just a matter of finding solutions of the coupled set of trigonometric equations that result from setting both time derivatives, in the spherical form of the equations of motion (Equations (37) and (38) Since the primary objective of this study is to see how dissipation influences the spin pole precessional dynamics, we need to do more. The approach we follow is to first find the nullcline curves, along which either v θ or v f vanishes. Intersections of these two spatial curves then yield the fixed points, at which both components of the velocity field vanish. These nullcline curves are useful in understanding the flow (Thomas 2006), as will be seen below.
Another point of interest is that the geometry of the flow, including locations of fixed points and nullclines, is independent of the actual values of the control parameters {a, b, c, δ}. That is, if we take a given parameter set and multiply all of them by a real constant q, so that both parameter sets yield identical geometries. Of course, the rate of motion of the spin pole precession will differ, by the factor q, but the trajectories will be identical. As a result, we could arbitrarily choose to normalize our parameter set so that, for example, c = 1. In terms of understanding the geometry of the precessional trajectories, this effectively reduces the dimensions of the parameter space from four to three. For the nondissipative case, this is exactly what Henrard & Murigande (1987) did. In that case, the equations of motion were collapsed from a threedimensional parameter space to one with two dimensions.
Other than in Section 8, where we explicitly consider the Moon, we will use values for the rate parameters {a, b, c, δ} that are all positive and of order unity. As noted above, multiplying all the rate parameters by a constant does not change the trajectories or locations of the fixed points.
A negative value for δ would lead to growth, rather than decay, of the obliquity. As seen in Equation (41), a change in sign of c changes the location of the nullcline, with θ → −θ. Other combinations of sign changes for the rate parameters generally correspond to reflections about the equator and/or prime meridian.

Low-order Solutions
Before commencing our fully numerical approach to finding fixed points, we examine some low-order approximations.
Along the θ nullcline, we have

Latitudinal Nullcline
We now examine the full, nonlinear behavior of the θ nullcline. If we set the v θ component of velocity in Equation (39) The first is just the unit sphere, and the second is a circle in the x-y planes. The circle has radius S c 1 2 , 5 8 d = ( ) and its center is at In that notation, Equation (57) is equivalent to The θ nullcline curves thus lie on the intersections of the unit sphere with a circular cylinder. The cylinder has its axis parallel to the z-axis, with radius and center as noted above. The exact geometry of the nullcline curves depends on the values of the parameters c and δ. There are three different fundamental cases for the pattern of the θ nullcline. The case is parameterized by the ratio δ/c and is thus a measure of the dissipative rate term δ, compared to the inertial rate term c, which is defined in Equation (10). In each of the three cases, the θ nullcline curves divide the surface of the unit sphere into either two (S > 1/2) or three (S 1/2) disjoint regions, with differing signs for the θ derivative of the flow field.
For the equilibrium case Abs[δ/c] = 1, the θ nullcline curves have the appearance shown in Figure 2(a). In this case, there are two separate curves, which divide the surface of the unit sphere into three regions. Figure 2. (a) Equilibrium case. Spatial pattern of θ nullcline when δ/c = 1. In this case, the nullcline divides the surface of the sphere into three disjoint regions, with opposite signs for the θ component of velocity. (b) Obliquity tidal dissipation dominant. Spatial pattern of θ nullcline curves when dissipation rate parameter δ exceeds parameter c. In this case, the nullcline divides the surface of the sphere into three disjoint regions, with differing signs for the θ derivative of the flow. (c) Inertial precession dominant. Spatial pattern of θ nullcline, in cases where dissipation rate parameter δ is less than inertial precession parameter c. As with the case shown above, there are two disjoint lines and regions.
If the dissipation rate δ is larger than the inertial precession parameter c, then the form is as shown in Figure 2(b). There are now, as in the low-dissipation case, two distinct curves and three separate regions.
If the dissipation rate is reduced, we enter the region dominated by the inertial precession parameter c and the θ nullclines change to the form shown in Figure 2(c). Now, there is a single curve that divides the surface into two regions.

Longitudinal Nullcline
We now consider the f nullclines. Setting the velocity component in the f-direction, as given in Equation (42) The general appearance of these nullcline curves changes when we vary the input parameters {a, b, c}. For the "internal" case, when the radius of curvature R of the conserved scalar parabola, as given in Equation (25), has value R < 1, the nullcline curves have the form shown in Figure 3(a).
In this case, there are three distinct curves, one near each pole, and one at low latitude. The low-latitude curve spans all longitudes, and the near-polar curves are limited to the longitude range 0 < f < p. These three curves divide the surface of the unit sphere into four disjoint regions.
If R = 1, for the osculating case, then the two upper (λ > 0) nullcline curves merge, as is shown in Figure 3(b). The three curves still divide the unit sphere into four regions, but three of the regions join at a single point.
In the external case where R > 1, there are only two separate curves, each of which intersects a pole, as is shown in Figure 3(c). These two curves divide the unit sphere into three disjoint regions.
Similar to what was seen with the nullclines in Figure 2, these nullclines separate the sphere into either three (R > 1) or four (R 1) disjoint regions, with differing signs for the fcomponent of the velocity of the flow.
When we convert Equation (61) to Cartesian form, it is equivalent to two simultaneous constraints: The first is just the unit sphere. The second is a cubic curve in the y-z planes, with no dependence on x values. Intersections of that surface with the unit sphere are the curves shown in the figures above.
Because Equation (65) has no dependence on x, we can examine the nullcline in the y-z plane. Figure 4 shows this cubic curve, overlaid on the unit circle for the same values of {a, b, c} used in Figures 3(a)-(c).

Vorticity and Divergence of the Flow Field
In much the same way that three-dimensional vector fields in seismology (Alterman et al. 1959) and geomagnetism (Backus 1986) can be described in terms of spheroidal and toroidal modes, we can usefully examine the curl and divergence of the flow on the surface of the unit sphere associated with our model of damped spin precession. In the study of two-dimensional, nondissipative, incompressible fluid flow on the sphere, which encompasses much of oceanic and atmospheric dynamics, vorticity and its mean square (enstrophy) play important roles complementary to that of energy (Charney 1971;Verkley & Lynch 2009 These components can be written in terms of the control parameters (a, b, c, δ) and unit vectors (s n m , ,ˆˆ) as This partitioning clearly illustrates that the nondissipative parameters (a, b, c) completely control the vorticity, while the single dissipative parameter (δ) controls the divergence of the flow.

Stability of Fixed Points
It is useful to characterize the stability of the fixed points. For that purpose, we will examine how the flow changes in the immediate vicinity of such a point. We give the names F and G to the functions that yield the time derivatives of angular coordinates, as in Equation (13), In terms of these functions, the local behavior, near a point on the unit sphere P , , 79 where the incremental position vector is and the Jacobian, or matrix of partial derivatives, is given by where the derivatives are evaluated at the fixed point X = 0.
Explicit values for these derivatives are The solution to this linearized set of equations can be written as This matrix exponential (Moler & Van Loan 2003) can be written in terms of the eigenvectors and eigenvalues of the Jacobian matrix J. For a stable fixed point, the real parts of both eigenvalues must be negative.
The eigenvalues of the Jacobian matrix J, as given above, are The sum of the eigenvalues is thus The discriminant, which marks the dividing line between real and complex eigenvalues, is given by If the discriminant is positive, both eigenvalues are real; if it is negative, the eigenvalues form a complex conjugate pair. When applied at fixed points, for which F = G = 0, if the real parts of both eigenvalues are negative, the location is a sink. If both are positive, it is a source. If they are of mixed sign, it is a saddle.
In the system under investigation, we previously noted that multiplying all the system parameters {a, b, c, δ} by a real constant q does not change the geometry of the fixed points and nullclines. However, if the constant q is negative, it will change sources into sinks, and vice versa.
We note that the Hartman-Grobman theorem (Grobman 1959;Hartman 1960;Wiggins 1990) effectively states that the behavior of a dynamical system near a hyperbolic fixed point is qualitatively the same as the behavior of its linearization near that fixed point. This implies that the characterization of the types of fixed points (source, sink, saddle) associated with a parameter set {a, b, c, δ} is actually applicable over the entire sphere.

Poincaré−Hopf Index Theorem
Poincaré provided a very important constraint on the number and types of fixed points that can exist for a continuous vector field on a sphere. Consider a continuous vector field V and a closed curve h, both defined on a manifold M. At each point P, on the manifold M, there is a value of the vector field V. If we move the point P around the curve h and return to the starting point, the vector field along the curve will have rotated an integer number of times. That number is the winding number associated with V and h. If the curve h is continuously deformed, the winding number will not change, unless it passes through a critical point of V. The index of a vector field, at a critical point, is the limit of the winding number of the field, on a small circle centered on the point, as the circle shrinks to zero radius.
The Poincaré−Hopf index theorem (Milnor 1965) states that S, the sum of the indices of the critical points of a vector field over a manifold, is equal to χ, the Euler characteristic of the manifold. In the case of a sphere, the Euler characteristic is 2. The index of a source or sink is +1, while the index of a saddle is -1.
As a result of this, when we have two fixed points with dissipation, one is a source and the other is a sink. When a third fixed point appears, it must be a saddle. If there are four fixed points, one is a saddle, and the others comprise either one sink and two sources or two sinks and one source.
We will now examine the situations in which the dissipative spin pole precession parameters {a, b, c, δ} are such that there are two, three, or four fixed points.

Two Fixed Points
Because the time derivatives, given in Equations (41) and (42), depend only on the parameters (a, b, c, δ) and the current position (θ, f), there is only one trajectory passing through each nonsingular point. In the case where there are only two fixed points, we can say more. The single trajectory passing through any nonsingular point has a heteroclinic trajectory passing through it, which starts at one of the fixed points and ends at the other. One of the fixed points is thus an attractor, and the other is a repeller (Milnor 1985;Sprott & Hoover 2017;Li et al. 2020). In that case, the basin of attraction, for the attractor, is the entire sphere.

Three Fixed Points
Locations in the four-dimensional parameter space {a, b, c, δ} at which there are three fixed points for the flow associated with the differential Equations (41) and (42) are of interest mainly because they compose the boundary between the larger domains that have two and four fixed points.
We saw above that, when there is no dissipation (δ = 0), a relatively simple criterion (Equation (25)) marks the boundary between two and four fixed points. When dissipation is included, the geometrical requirement for three fixed points is that one of the intersections between the θ nullclines and the f nullclines is a tangent intersection. Equating the slopes along the θ and f nullclines yields another constraint equation Given the three constraint equations for a tangent intersection of nullclines (44,45,89) and the six relevant parameters, comprising four rates {a, b, c, δ} and two coordinates {θ, f}, we can specify any three parameters and then solve for the other three.
As an example, we can specify a = 2.8, b = 1, c = 1, and solve for the dissipation parameter (δ = 1.0724) and tangent intersection point (f = 3°.48; θ = 111°.45). Figure 5(a) shows a related but slightly different case (δ = 1.07), in which there are four fixed points, but two of them are close to merging. The fixed points are color-coded: green is the sink, black is the saddle, and both red and blue are sources. In this figure and several following, blue is the source on the left-hand side, while red is the source on the right-hand side. In addition to the nullclines and fixed points, we also show the vector field, as given by Equations (39). Figure 5(b) shows another tangent intersection case, with identical values for a, b, and c, but a smaller value for the dissipation parameter (δ = 0.29). Note the differing geometry of the tangent nullcline intersections. That is, the saddle is close to merging with the red (right-hand side) source.
We note that since an actual case of three fixed points has a single source and a single sink, dynamical simulations started at any initial position on the unit sphere will eventually converge on the sink.
The cases shown in Figures 5(a) and (b) are at the upper and lower ends of the range of values of the parameter δ, for the given values of {a, b, c}, which yield four fixed points. Higher or lower values yield only two fixed points. This transition in number of fixed points is similar to that seen in the conservative (δ = 0) cases discussed above.

Four Fixed Points
At some level, the most interesting case is that in which there are four fixed points. From the index theorem, we know that one of them is a saddle and the other three are centers, sources, or sinks. In the nondissipative case, all three of them are centers. In the dissipative case, we expect either two sinks and a single source or the opposite case with two sources and a single sink. In all the cases we have examined, a positive value for δ yields a single sink. However, we have no proof of that conjecture. Figure 6 illustrates various aspects of the behavior of the spin pole precession system with control parameter values a = 2.8, b = 1, c = 1, and δ = 0.6. These values lead to four fixed points, consisting of two sources (blue and red), one saddle (black), and one sink (green).
Figure 6(a) shows the fixed points and nullclines, along with the vector flow field. As was the case in Figure 5, the flow field not only confirms the nullcline information but also interpolates the directions of flow over the entire unit sphere. Figure 6(b) again shows the nullclines and fixed points and also illustrates integrated trajectories that start near the saddle. The red lines started to the right of the saddle, while the blue lines started to the left. The solid lines were integrated backward in time, and both converged on the single source, which upon time reversal acts as a sink. The dotted lines were integrated forward in time, and each converged on one of the sinks. Figure 6(c) illustrates the domains of attraction of the two sinks. In addition to the nullclines, fixed points, and integrated trajectories, which are shown in Figure 6(b), this figure also indicates initial positions for which the trajectories, when integrated forward in time, converge on the two different sinks.
We ran 1923 simulations with initial conditions uniformly spaced over the surface of the sphere, based on the icosahedral grid of Baumgardner & Frederickson (1985), and integrated the trajectories forward in time, to see which of the two sources each of them emerged from. The small red circles indicate initial conditions that emerged from the left-hand source, and the small blue circles are initial conditions that emerged from the right-hand source. In this case, of the 1923 simulations, 1549 are red and 374 are blue.
We note that, with a change in sign of all four control parameters (a, b, c, δ}, or equivalently, a change in the sign of the time variable, we would have two sinks and one source, and these colored circles would represent the basins of attraction for the sinks. Without those changes, they are the basins of repulsion of the positive parameter system. This tabulation of initial conditions, and noting which of them lead from each of the two sources, clearly shows that the trajectories integrated from the vicinity of the saddle cleanly delineate the basins of repulsion. In the vicinity of the sink, the basins of repulsion are tightly interleaved. However, they are not chaotic (Nusse & Yorke 1996).
In the case of nondissipative spin pole precession, Henrard & Murigande (1987) have given analytic expressions for the areas of regions on the unit sphere containing each of the Cassini states. We expect that a similar result can be found for the areas of the basins of attraction or repulsion in the dissipative case, but we do not yet have the answer.

Application to the Moon
We now apply our theory to the Moon. There are several motivations for doing so. The Moon is, of course, the body to which the concept of a Cassini state was first applied. It also has the advantage of very precise measurements of its orbital and rotational motions, via lunar laser ranging (Dickey et al. 1994), and the gravity field, via the GRAIL mission (Williams et al. 2014).
We wish to see how well our simple model does at two tasks. The first task is to give the observed obliquity of the Moon. The second is to explain the small observed offset of the lunar spin In particular, Organowski & Dumberry (2020) have argued that part of the observed displacement is due to friction at the core−mantle boundary and that only 0 15 is due to obliquity tides. Colombo (1966) gave the first treatment of Cassini states in terms of a Hamiltonian formulation but mainly focused on the axisymmetric (A = B) case. Peale (1969) extended the formulation to include triaxial bodies (A < B < C). In application to the Moon, he developed a constraint equation that relates the two key angles of the problem, which are orbital inclination i and spin pole obliquity ε.
To second order in orbital eccentricity, the constraint can be written as As noted by Peale (1969), the contribution of the equatorial moment difference (B − A) is not large. In fact, it only emerges at third order in obliquity. When the equation is linearized in the small parameters e, i, ε, and γ/n, the solution has the form n C C A i 2 3 .

Orbit Precession
Our model assumes that the orbit plane is uniformly precessing, at fixed inclination about an inertially fixed pole. Is that a good approximation, when applied to the Moon?
Orbital elements of the Moon, from a numerically integrated solution, are available from the JPL HORIZONS website, https://ssd.jpl.nasa.gov/horizons.
We obtained such elements for 200 yr (1900-2100), at 1-day time steps. The reference frame is the J2000 ecliptic and equinox.
The inclination has no obvious linear trend, over a time span of a few hundred years. On much longer times, the inclination does have secular variations (Touma & Wisdom 1994). At the present epoch, the mean value is The precessional motion of the orbit, as reflected in motion of the longitude of ascending node, is well approximated by a linear trend The mean rate yields a (retrograde) nodal precession period of 2 18.5996 yr. 96 The next several figures show aspects of the lunar orbit pole motion, over a 20 yr time span, centered on 2020. They cover just over one cycle of nodal precession. Figure 7 shows variations in inclination. Most of the variation occurs at two cycles per year. Figure 8 shows the residual variations in longitude of the node, when the best-fitting linear trend is removed. Figure 9 shows a projection of the orbit pole unit vector onto the ecliptic plane. In this view, the semiannual nutation is quite visible. In addition to variations in the orbit pole orientation, the speed of orbital motion also varies. The largest variations are, of course, due to the finite eccentricity of the orbit. However, even the mean motion has time variations, mainly due to solar

Rate Constants
In the next step, we need values of the nondissipative rate constants {a, b, c}. We can then estimate where the lunar spin pole would be oriented, without the influence of dissipation.
The principal moments of inertia (A < B < C) are related to the gravitational potential spherical harmonic coefficients via (Yoder 1995 These parameter values place the Moon in the regime where there are only two fixed points for the spin pole precession. The radius of curvature of the parabolic Hamiltonian, as given by Equation (25), is 3.31, which is well above the transition point.
However, early in its history, the Moon presumably had access to four possible Cassini states, and the transition from four to two was likely quite dramatic (Ward 1975;Ćuk et al. 2016, 2019.

Dissipation Rate
Finding a robust value for the parameter δ is somewhat problematic, in that we have no direct observations to constrain it. The speed at which the Moon's spin pole was initially damped into a Cassini state is unknown. However, we can use observations of the current spin pole orientation to provide some information. In particular, we can take either of two approaches. In one approach, we solve for the value of δ that provides the observed offset in longitude of the damped spin pole, in the vicinity of the conservative spin pole noted just above. In an alternative approach, we use Equation (32) and values of the relevant parameters, to obtain an estimate of δ.
The value of δ required to produce an angular offset of 0 15, as suggested by Organowski & Dumberry (2020), is 6.04 10 s . 106 This value depends only on the adopted values of the inertial rate parameters {a, b, c}, as given above in Equation (104), and shows that the dissipative rate parameter for the Moon is much smaller. In the figures presented above, we have mainly used values where δ is comparable in magnitude to the other rates. This suggests that the dissipative rate parameter is much smaller than the inertial rate parameters {a, b, c}.
The second approach yields a quite different value. We find that the value for δ obtained from Equation (32), using values of k = 0.030 2 ± 0.001 2 and Q = 26.5 ± 1 from Dickey et al. In either case, values obtained this way are roughly 6000 times larger than the estimate in Equation (106). It is not obvious how well we should expect the two results to agree. Our model makes several assumptions that are not strictly applicable to the Moon, including a circular orbit and a uniform rate of orbit plane precession. However, the obliquity value for the fixed point is recovered with reasonable accuracy.
Instead, it appears that our arguments leading to damping time estimates, in Section 4.1, appreciably overestimate the damping rate. As noted earlier, this suggests that much of the tidal dissipation, associated with the obliquity tide, simply heats the body, rather than reorienting it.

Summary and Conclusions
We have examined the dynamics of spin pole precession of triaxial bodies, moving along a uniformly precessing orbit, under the influence of gravitational torques and tidal dissipation. The inclusion of dissipation in the torque balance changes the geometry and stability of the fixed points.
In the canonical Cassini state analysis, where dissipation is ignored, there are two, three, or four fixed points of the precessional motion, depending on the relative rates of spin pole and orbit pole precession. These fixed points are known as Cassini states (Colombo 1966;Peale 1969).
When there are four fixed points, three of them are centers, and one is a saddle. In that case, two homoclinic trajectories emerge from the saddle, and they divide the surface of the unit