CC BY
3Russian Journal of Nonlinear Dynamics, 2022, vol. 18, no. 4, pp. 527-541. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd221001
NONLINEAR PHYSICS AND MECHANICS
MSC 2010: 70F15
Rotation of a Planet in a Three-Body System: a Non-Resonant Case
O. M. Podvigina
We investigate the temporal evolution of the rotation axis of a planet in a system comprised of the planet (which we call an exo-Earth), a star (an exo-Sun) and a satellite (an exo-Moon). The planet is assumed to be rigid and almost spherical, the difference between the largest and the smallest principal moments of inertia being a small parameter of the problem. The orbit of the planet around the star is a Keplerian ellipse. The orbit of the satellite is a Keplerian ellipse with a constant inclination to the ecliptic, involved in two types of slow precessional motion, nodal and apsidal. Applying time averaging over the fast variables associated with the frequencies of the motion of exo-Earth and exo-Moon, we obtain Hamilton's equations for the evolution of the angular momentum axis of the exo-Earth. Using a canonical change of variables, we show that the equations are integrable. Assuming that the exo-Earth is axially symmetric and its symmetry and rotation axes coincide, we identify possible types of motions of the vector of angular momentum on the celestial sphere. Also, we calculate the range of the nutation angle as a function of the initial conditions. (By the range of the nutation angle we mean the difference between its maximal and minimal values.)
Keywords: nutation angle, exoplanet, averaging, Hamiltonian dynamics
1. Introduction
The recently increased interest in the factors affecting the rotation of a planet is related to the dependence of the planet's climate on the evolution of its rotation axis, while stable benign climate is favourable for the development of intelligent life [2, 15, 16, 18, 40]. For example, the obliquity (the angle between the rotation axis and the normal to the orbital plane) of the Earth
Received July 03, 2022 Accepted August 22, 2022
This research was carried out at the Moscow Aviation Institute and supported by grant 22-21-00560 from the Russian Science Foundation.
Olga M. Podvigina [email protected]
Institute of Earthquake Prediction Theory and Mathematical Geophysics, RAS ul. Profsoyuznaya 84/32, Moscow, 117997 Russian Federation
varies just between 22.1° and 24.5° and the eccentricity of its orbit between 0.0002 and 0.07, but even such small variations lead to the regular occurrence of glacial/interglacial cycles [33].
The rotation of the Earth on various time scales has been extensively studied in the literature, see, e.g., [14] and references therein. Numerical simulations of the rotation of the Earth for over 5Gyr [28] revealed that its obliquity remains relatively stable. By contrast, the obliquity of Mars is chaotic and varies between 0° and 85° [28, 41]. Such a difference in the behaviour of obliquities led to a conjecture that a heavy satellite has a stabilising effect. The conjecture was supported by the numerical results of Laskar and co-authors [26, 27]. The stabilising effect of the Moon and its necessity for the development of advanced exo-life was questioned by Lissauer and co-authors [31].
The models explaining the retrograde rotation of Venus were discussed in [12, 13]. The long-term evolution of giant planets, in particular, the scenarios leading to the recent tiltings of their rotation axis, were considered, e. g., in [8-11]. The rotation of Saturn under the influence of the Sun and Neptune, assuming resonant interaction between Saturn and Neptune was studied by Krasilnikov and Amelin [22]. The authors have shown that the direct impact from Jupiter on the rotation of Saturn prevails over the indirect impact from any Solar system planet.
The rotation of a planet in particular exo-systems was studied, e.g., in [4, 37, 39]. The problem in a more general setup was considered in [38], where the role of resonances and the possibility of chaotic behaviour of spin-axis dynamics of an exo-planet were investigated.
In the problems of celestial mechanics averaging is often applied to study the rotation of a rigid body on large time scales, as, e.g., in the classical works of Beletskii [5-7], Lidov [29, 30], Kozai [20] or Kinoshita [19]. Markeev and Krasilnikov [32] have shown that in a restricted three-body problem with prescribed quasi-periodic orbits averaging reduces the rotation of a satellite to an integrable problem and described possible types of motions. Krasilnikov and Zaharova [24] generalised the results to a system comprised of an arbitrary number of bodies. The approach was used for investigation of the temporal evolution of a rotation axis of a planet, in particular, the
dependence of the range of the nutation angle (the difference between its maximal and minimal values) on the parameters of the orbits of celestial bodies [23, 35]. In [36] the impact of an exo-moon on the rotation of a planet was studied. It was shown that the exo-moon may have a stabilising or a destabilising effect.
In these papers the satellite or the planet were assumed to be almost spherical, the difference between the largest and smallest principal moments being a small parameter of the problem. Other bodies were point masses. Celestial bodies comprising the system moved along quasi-periodic orbits with prescribed non-resonant frequencies u. The Hamiltonian describing the rotation of the satellite or planet was obtained by expansion in a power series in the small parameter and application of time averaging over motions associated with the frequencies u.
Here we follow this approach to study analytically the temporal evolution of the rotation axis of a planet on large time scales in a planetary system comprised of a star, the planet and its satellite, moving along quasi-periodic orbits. The planet and its satellite are called exo-Earth and exo-Moon, respectively. The exo-Earth is an almost spherical rigid body and the difference between the largest and smallest principal moments is a small parameter. The star and the satellite are assumed to be point masses. The orbit of the exo-Moon keeps a constant inclination to the ecliptic and undergoes two types of slow precessional motion, nodal and apsidal, with small respective frequencies an and aa. We assume that the orbital frequencies of the exo-Moon and exo-Earth are of the order of one and non-resonant.
We define the range of the nutation angle as
A(Io, h0) = max I(t, IQ, h0) - mm I(t, IQ, h0),
where I0 and h0 are the initial values of the nutation angle and the longitude of the spin axis and I(t, I0, h0) is the nutation angle at time t for these initial values.
The paper has the following structure:
In Section 2 we recall the Hamilton's equations for the rotation of a rigid body and averaging. The averaged equations involve six coefficients, which are constants in the moonless system and become time-periodic as the exo-Moon is added. The coefficients are computed given the masses and orbits of the celestial bodies. In Section 3 we show that in the star-planet-satellite system considered, under a canonical change of variables, the system becomes integrable. In Section 4, for the case of axially symmetric exo-Earths and its rotation axis coinciding with the symmetry axis, we describe three types of motion of the angular momentum axis on the celestial sphere and calculate analytically the range of the nutation angle in each of these three cases. Finally, we briefly summarise our results and indicate possible directions for further studies.
2. Equations of motion
In this section we recall the averaged Hamilton's equations for the rotation of the exo-Earth at large time scales, see [35] or [36]. The temporal evolution of the exo-Earth's angular momentum vector is considered in a system comprised of N + 1 celestial bodies. The orbits of celestial bodies are quasi-periodic with the respective frequencies divided into two groups: the order one frequencies that are non-resonant and other frequencies that are small. Averaging over the fast variables associated with the order-one frequencies implies a Hamiltonian which involves six quasi-periodic functions Dj (t), the respective frequencies being the small frequencies of the original system. For a more detailed discussion of the procedure of averaging see [23, 35, 36] or [21].
2.1. Hamilton's equations
Denote by OXYZ a non-moving inertial reference frame, by a coordinate system
with origin at the center of mass of the exo-Earth and axes parallel to those of the OXYZ, and by Mxyz a coordinate system with the same origin and coordinate axes coinciding with the exo-Earth's principal axes, where Mz is associated with the maximum moment of inertia. We employ the Andoyer variables [1] (as in [19], they are denoted by (G, H, L, g, h, l)) defined as follows: G is the magnitude of the exo-Earth angular momentum vector L; H is the Z-component of L; L is the z-component of L; g is the angle between the intersections of the plane £ with the planes and Mxy; h is the angle between the axis M( and the intersection of the planes £ and M(n; l is the angle between the axis Mx and the intersection of the planes £ and Mxy and £ is the equatorial plane orthogonal to L.
The Hamilton's equations for the rotation of the rigid exo-Earth then are
d dW d dH
1h l) = d(G, H, L) ' M{G> H> L) = ~d(g,h, l)1
where the Hamiltonian is
H
G2 - L2 /sin21 cos2l
+
A
B
N
L2
n=1
(2.1)
(2.2)
A ^ B < C are the principal moments of inertia of the exo-Earth, and Vn is the potential energy of the gravitational interaction with the nth celestial body.
The smallness of the radius of the planet compared to the distances between celestial bodies implies that
t A ,,2 i r. 2 i ^< 2
Vn 2m
iAin ,1 + B^n, 2 + C^n, a), Vn = fmn,
(2.3)
where f is the universal gravitational constant, mn is the mass of the body, Rn is its geocentric distance, and (Yn 1, Yn 2, Yn 3) are the direction cosines of the vector Rn = (RnX, RnY, RnZ) from the geocenter to the nth body.
Assuming the exo-Earth to be almost spherical, we write its moments of inertia
A = J0 + eA1, B = J0 + eBx, C = J0 + eC1, where 4- < 1 and J0 = A+B+c. Hence, the Hamiltonian (2.2)-(2.3) takes the form
n
2 Jn
+ eH1 + o(e),
(2.4)
where, by (2.3),
/sin21 cos21 L Cl + (G — L ) I-— +
B
1
3 N
+ 2 E % (A+ + Ci7,1s)
n=1 n
(2.5)
Here h, l, G, H and L are slow variables, while g is a fast one.
The rotation of the coordinate system Mxyz relative to can be described by three
Euler's angles. The angle between the axes Mz and M( is called the nutation angle. If the vector of the angular momentum is oriented along the positive direction of Mz, which is the case considered in Section 4, then the nutation angle is I = arccos (|r).
2.2. Averaging
Denote by ш = (w, ш1, ..., wK) the K + 1 prescribed order-one frequencies of the motion of the considered N + 1 celestial bodies (w is the rotation frequency of the exo-Earth and (w1, ..., wK) are the orbital frequencies) and by F the average of a function F(u)t, и..., ojKt) over the motions associated with these frequencies. In the absence of resonances between the frequencies,
П
SlL
2Jn
+ £ Щ,
(2.6)
where the mean Hamiltonian is
H1 _ F (G, L, l)Ç(G, H, h) + C.
(2.7)
Here
2 72
HG, = ^
G2 — L2 f sin2 l cos2 Л L2
+ +
A
1
B
2CU
J02 A1 + B1 - 2C1
2eC1
+
G
2
Q{G, I, h) = —g + sin21 (~D1 sin2 h - D2 cos2 h + D3 + DA sin(2/?,)) -27o
— sin 2I(D5 sin h — D6 cos h),
N
С
ßnDnR — G ) , 1
n=1
arccos
N
D
D D
(n) _ 1 _ HnDnX2 ,
(n) 4
VnDn ßnDnXY,
Y,D)n), 1 < j < 6
n=1 (n)
D D
2
(n)
№nDnY2 , ßnDnXZ,
D
D
(n) 3
(n) 6
VnDnZ2, ßnDnYZ 5
(2.8) (2.9)
(2.10) (2.11) (2.12)
and
D
2п 2п
1
D
nr (2n)K+1 1
npv (2n)K+1
DnR =
(2n)K+1
r2
Rnp
dJÖdJÖ! ... d0K,
00 2n 2n
R R
-^LJ^dede, R
de
K,
00 2п 2п
R3
dede1... deK,
(2.13)
where p and v denote X, Y or Z, в = wt and вк = шкt, 1 ^ k ^ K.
nj
Since -q^ = 0, the averaged system has a first integral G = const.
3
1
2
1
0
0
3. Hamilton's equations for the system comprised of the exo-Sun, the exo-Earth and the exo-Moon
3.1. Coefficients Di related to the exo-Sun and the exo-Moon
(2)
In the system comprised of a star, a planet and its satellite, the coefficients Di = D(1) + D( are the sums of D(1) that result from the torque from the Sun and the sums of D(2) from the Moon. Let aE be the semi-major axis and eE the eccentricity of the planet's orbit around the star. Then (see [35] or [36])
D[l) = 41} = 41} = D[l) = 41} = D™ = 0. (3.1)
E
The orbit of the exo-Moon is a Keplerian ellipse with the exo-Earth in one of the focuses. The inclination i of the lunar orbit to the ecliptic plane does not change in time. The longitude of the ascending node and the argument of the periapsis evolve as
Q = Q0 + ant, w = w0 + aat.
According to [36], the coefficients related to the exo-Moon are
d(2) = 2 (cos2 Q sin2 i + cos2 i), D22) =2 — cos2 Qsin2 i + 1), d32) =2 sin2 i, d42) = 2 sin Q cos Q sin2 i, D^2) = —2 sin Q sin i cos i, D^2) = 2 cos Q sin i cos i,
where
(3.2)
fm
M
2a3 fi p2 A3/2' 2aM I1 eMJ
and mM, aM and eM are the mass, the semi-major axis and the eccentricity of the exo-Moon's orbit around the exo-Earth.
3.2. The canonical change of variables
Since Di = D(1) + D(2), for the system comprised of a star, a planet and its satellite, the function G (2.9) takes the form
G = ^T2+ sin21 (- (d[1] + 42)) sin2 h - (41} + 42)) cos2 h+
(2) (2)
+d32) + d42) sin(2h)) — sin 2I (dJ2) sin h — Df] cos h). (3.3) Substitution of (3.1) and (3.2) into (3.3), followed by a series of algebraic transformations, yields
G2
G = +sinz/ +5 ^ - sin2-i - cos2-i ) +
sin2 i cos(2Q — 2h) J + ^ sin 2/sin 2?; cos(Q - /?,). (3.4)
Using the generating function F2(h, H') = H'(h — Q) where Q = Q0 — ant, we find that in the canonical coordinates (H', h') = (H, h — Q0 — ant) the Hamiltonian (3.4) takes the form
H1 = F(G, L, l)g(G, H', h')+ C — anH'. (3.5)
Therefore, the averaged system has a complete set of three first integrals: G = const, F = const and FG — anH' = const.
4. Calculation of the range of the nutation angle
In this section we consider the evolution of the angular momentum and calculate the range of the nutation angle for the system shown in Fig. 1 under the assumptions that the exo-Earth is axially symmetric (A1 = B1) and the rotation axis coincides with the symmetry axis of the body (J = 0). The angular momentum L = (Lg, Ln, L) takes the form (see the definition of Andoyer variables)
Lg = G sin h sin I,
Ln = —G cos h sin I, (4.1)
Lz = G cos I.
From (2.1), (2.7)-(2.10) and (3.5) the temporal evolution of I and h' = h — Q0 — ant satisfies
dh' _ 1 dg' dl _ 1 dg'
lit ~ Gsin/ dl ' ~dt ~ _Gsin/ dh!' ( '
where
g'(I, h') = sin2 I(—D' + a cos 2h') + 3 sin 2I + a' cos I, (4.3)
D> = |(C,-A,) + 2 (l - \ sin2 i) ), a = - Al)S sin2 i,
3^ ( )
(3 = — A1)Ssin2-i and a' = Gan.
Note that Eq. (4.3) is invariant under the symmetry h' — —h'
g'(I, h') = g'(I, —h'). (4.5)
The equation is also invariant under the transformation (a, I, h') — (—a, n—I, h'+n), therefore, without loss of generality we consider non-negative a' only. From (4.2) and (4.3) we have that
dh' 1
— = —--(sin2/(-JD/ + a cos 2//) +2/3 cos 2/cos to' - a'sin/),
dt G sin I (4 g)
— = — (2 sin la sin 2h! + 2 cos 1/3 sin h'). dt G
The assumption that the angle i is small implies that a ^ /. To compare D' with /, we note that D' ~ while /3 = ^^ ~ q^fm ^^ = 3.(0,-+)^ ^ fhe Earth and fhe
re 2 rm 2
Moon we have that. -§- ~ mg Jf & 30, i.e., that. D' is significantly larger than (3.
1—1 ^^ j-i R j-i
Below we derive an approximation for the range of obliquity depending on the parameters of the system and the initial condition. The expressions that we obtain are different depending on the number of the steady states (two, four or six) that exist in the system. Below I(t) and h(t) are solutions to (4.2) and (4.3). Often I is regarded as a function of h and I(h) instead of I(t) is assumed.
The steady states of (4.6) satisfy
= = o dh! dl
which can be re-written as
sin 2I(—D' + a cos 2h') + 2B cos 2I cos h' — a' sin I = 0; sin h (2 sin Ia cos h' + cos IB) = 0.
Hence, the steady states can be found from the following equations:
(4.7)
h' = 0, n, sin 2I(—D' + a) ± 2B cos2I — a' sin I = 0, (4.8) cos IB
cos /?'=---7-, sin2I(—D' + a cos 2/?/) + 2/3 cos 2/cos/?.' — cr'sin / = 0. (4.9)
2 sin Ia
Below we assume that, similarly to the Earth, D' » B. Under this assumption the steady states (4.8) are
2B 2B
Ks i = 0, Iss i-—,—--7; /4 2 = 0, Iss —--T; (4-10)
ss '1 ss '1 D' — a + a' '2 '2 D' — a + a'
a' a'
/?.' o = 0, cos I„„ o ------, /?', d = 7r, cos I„„ d &--—--(4.11)
' 2(D' — a) ' 2(D' — a) v ;
where the latter two1 exist only if \a'\ ^ D' — a.
To find the steady states satisfying (4.9), we note that B ^ a implies that \ cos I\ ^ \ sin I\, i.e., that I ~ §. Hence, we can write that, cos/?' —/3cos Substituting this into the second equation (4.9), we find that
aa' , /3d'
cos 1ss,5 ~ — o„ f n/ i „\ i /52' cosnss,5
2a(D' + a) + B2' '5 2a(D' + a) + B2
(and Iss 6 = Iss 5, h'ss 6 = —h'ss 5). Therefore, these steady states exist whenever
B\a'\
2a(D' + a) + B2
^ 1.
1 To find the steady states (4.11), we re-write the second equation in (4.8) as
sin I(2 cos I—D' + a) — a') ± 2B cos 2I = 0,
which due to the assumption D' ^ B implies that, unless sin I is small (this gives the steady states (4.10)), the solution to the equation satisfies \2cosI(—D' + a) — a'\ C B.
Overall, the system may have 2, 4 or 6 steady states. The above conditions for the existence of the steady states can be summarised as follows:
case I a' > D' -a two steady states (4.10), both centers
case II a' < D' — a and ßa' > 2a(D' + a) + ß2 four steady states: two centers (4.10), one center and one saddle (4.11)
case III ßa' < 2a(D' + a) + ß2 six steady states: two centers (4.10), two saddles (4.11), two centers (4.9)
(4.12)
The evolution of L on the celestial sphere in these three cases is shown in Fig. 2. Note that the motions of types I and II emerge in the study of rotation of a planet in an evolving orbit, in the absence of exo-Moon [7].
C
I
Fig. 2. Motion of L (4.1) on the celestial sphere for I and h solving Eqs. (4.2) and (4.3) in cases I, II and III
The maxima and minima of I for a particular trajectory (I(t), h'(t)) are achieved at = 0, which due to (4.6) takes place at
h = 0, n or at the points where cos h = —
cos Iß 2 sin Ia '
Hence, for a trajectory through (I0, h0) the extreme values of I, which we label by I(0), i(n) and I(ext) achieved at the three points given above, can be found by solving the equations
G'(I(0), 0)= s; G'(I(n),n)= s
gt^ext)} hiext) ) = ^ œg ¡¿ext) =
cos I (ext)ß
(4.13)
2 sin I(extt)a
for s = G'(I0, h0). The last equation in (4.13) cannot be solved for all trajectories. If it can be solved, it has two solutions (I(ext, h(ext)) and (I(ext, -h(extt)), due to (4.5). Below we study how the range of the nutation angle
A(Io, ho) = max I(i) - min I(i)
-<x<t<<x -<x<t<<x
for a trajectory (I(t), h'(t)) with (I(0), h'(0)) = (I0, 0) depends on I0, considering the three cases outlined in (4.12) individually.
Recall that we assume D' » /, which implies that for a particular trajectory the difference I(t) — I0 is small and we can write I(t) = I0 + I1(i). In case I, when no heteroclinic equilibria exist, the function I(t) can be regarded as a function of h'. Writing
s = G'(I0, n) = sin2I0(-D' + a) - ß sin 2I0 + a' cos I0, sin21 & sin210 + I1(h') sin 2I0 + I2(h') cos 2I0, sin 2I & sin 2I0 + Il(h')2 cos 2I0 - I?(h')2 sin 2I0,
cos I pa cos I0 - Ii{h') sin /0/i(/?•') - if(h') cos y
and substituting these into
G'(I, h') = sin2 I(-D' + a cos 2h') + ß sin 2I cos h' + a' cos I = s,
we find that can be found from
(4.14)
(4.15)
(If0)2 + bIn + c = 0,
a I
where
a = —D' cos 210 + 2/3 sin 210 — -a' cos I0,
b = —D' sin 2I0 — 2/ cos 2I0 — ct' sin I0, c = —2/ sin 2I0.
(In the first two lines in (4.17) we omit a using the fact that D' » a. The formula for the roots of the cubic equation implies that
(n)
\b\ - (b2 - 4ac)1/2
2a
(4.16)
(4.17)
(4.18)
with a, b and c given in (4.17). Hence, for a trajectory (I(t), h'(t)) that has two extrema of I(t), which as discussed above are achieved at h' = 0 and h' = n, the range of obliquity is
A(I0, 0) =
r(n)
(4.19)
with
r(n)
given by (4.18).
The two additional extrema (see (4.13)) are achieved at hiext\ where cos =
Therefore, they exist only for trajectories such that
cos I0/3
^ 1. (Here we use the fact
that. /(e-Ti) pa /0.) Since a <C (3, this implies /0 ~ f and cos pa — cosJ°/3. Fr-om (4.16) and (4.17), for such trajectories the expression for simplifies to
An)
4ß cos I0
a'
By the same algebraic transformations as above we obtain
j(ext) _ ~2a-2¡3 cos I0
(4.20)
(4.21)
1
1
1
sin I0a
1
(n)
and
T(n) T(ext) I1 - I1
The range of obliquity is the maximum of I imply that
2max(2|/cosI0\, \a + /cosI0|)
, therefore (4.20) and (4.21)
A(Io, 0) =
a'
(4.22)
In case II (see Fig. 2) heteroclinic trajectories through the steady state (Iss 4, n) split the celestial sphere into three regions, comprised of a center and a set of trajectories around this steady state. Two of the centers are located near poles and we call the respective regions polar, while the remaining one2 (Iss 3, 0), in general, is not. By contrast, we call equatorial the region near (Iss 3, 0). Inside a region the range of obliquity depends continuously on the initial condition, while it is discontinuous when crossing a boundary.
For a trajectory (I(t), h'(t)) inside the equatorial region the maximal and minimal value of I(t), Imax and Imin, are both achieved at h' = 0. Moreover, as we noted above, I(t) does
not differ much from I0 and (Iss 3, 0) is a center, therefore, Imax — I,
ss, 3
I o- I ■ .
ss, 3 min
The
initial (I0, 0) corresponds either to the maximum or to the minimum of I(t), implying that
A(I0, 0) = Imax — Imin ~ 2\Iss,3 — Io|.
(4.23)
The region is bounded by a heteroclinic trajectory through (Iss 4,n), where Iss 4 = Iss 3 (see (4.11)). The initial condition (I0, 0) belongs to this region if Iss 3 — 5het < I0 < Iss 3 + + 5het. Following the same ideas that are used to calculate A, we find that the value of 5het is a solution to the following equation:
G'(Iss,3 + het^ 0) = G'(Iss,4, n).
Solving the equation, we obtain
àhet =
2ß\ sin 2Iss,3\ \1/2
D'
In polar regions the range is calculated similarly to case I. Namely, when a trajectory through (I0, h0) has one minimum and one maximum, the range is given by (4.18) and (4.17). For the
case of four ext.rema and /0 ^ f we have
A
( ext)
4ß2 cos210 + 2a2 - 2aß cos I0
_2DU£ + 2aß + aa,
If there are four ext.rema and /0 ^ f, then
(ext) T(n)
A = I(ext) - I,
(4.24)
(4.25)
2To show that the steady state (Iss 3, 0) is a center, we note that in the coordinates (It = I(t) — Iss 3, hx = h'(t) — hss 3) nearby trajectories satisfy the equation
It(D' + 2ß sin 2Iss n) + -2asin2 Iss, -
/^sin 2Iss,3 2
In case II we have that
, , , ( o ßsin2/„„ o\
(D + 2/3 sin2/ss 3) 2a sin /ss 3 - 2 ",3j > 0,
which implies the statement.
1
where I(ext) is given by (4.24) and
I
(n)
4 ft cos I0
2D' cos I0 + 4ft + a'
(4.26)
In case III (see Fig. 2) there are three steady states inside the equatorial region. One is (Iss 3, 0), which now is a saddle, and the other two, (Iss 5, hss 5) and (Iss 6, hss 6), are centers. The meridian of the initial conditions (I0, 0) does not cross the boundaries of the emerging regions around the latter steady states.
For a trajectory inside the equatorial region that takes extreme values only at h' = 0 the range can be found from (4.23). If additional extrema at h' = hext (and also at h' = -hext) exist, then both maximal and minimal values along a trajectory are taken at this value of h'. The maximal and minimal values of I1(i) are solutions to the cubic equation (4.16). They are
(min ,max) 1
-b ± (b2 - 4ac)1/2
2a
where
Therefore,
a = Df, b = cos IQ ^-2D' + - a', c = -2a + 2/3 cos I0
A IT ^ _ r(max) r(min) _ \b\ _ \ cos I0 -2D'a + 2ft2) - a'a\
¿M - h ~h -777--7UZ-
a|
The results can be summarised as follows:
case A: equations #
I |/0-f| >a/3 (4.19) and (4.18)
ko-i <<*P (4.22)
II \l0 - § > af3 and |/0 - Iss 31 > 6het (4.19) and (4.18)
\l0 - § < ap and I0 < § (4.24)
/0 - § < ap and I0 < § (4.25), (4.24), (4.26)
14) ~~ hs,3 < $het (4.23)
III \I0 ~ Iss,3 > ^het (4.19) and (4.18)
\h ~ Iss,s\ < 6het and lJ0 - Iss,s\ > aP (4.23)
\l0 - § < a/3 (4.29) and (4.28)
(4.27)
(4.28)
(4.29)
(4.30)
The range of the nutation angle calculated from these approximations is in good agreement with the one found by numerical integration of Eqs. (4.6), see Fig. 3.
5. Conclusion
In this paper we have studied the evolution of the rotation axis of a planet at large times in the system comprised of a star, the planet and its satellite, assuming that the planet is a rigid almost spherical body, the star and the satellite are point masses, and that the frequencies of the motions of celestial bodies are non-resonant. As in [35, 36], the evolution of the rotation axis
Io Io
(a) (b)
Fig. 3. The dependence of A(I0,0) on I0 found by integrating Eqs. (4.6) over time (black line) and from approximations (4.30) (gray line) in cases I (a) and II (b)
of the planet is governed by a Hamiltonian involving six parameters. The parameters depend on ant, where an is the precession frequency of the satellite. Upon a canonical change of variables the system becomes integrable and has a complete set of three first integrals.
Assuming an axially symmetric planet, we identify three possible types of motions of the vector of the angular momentum on the celestial sphere and derive inequalities determining which of these types takes place in a particular system. For each of these three cases we calculate analytically the range of the nutation angle as a function of initial conditions.
A natural continuation of this study is to examine planetary systems with resonances: resonances are common in the Solar system [34] and we expect them to occur in other planetary systems as well. In the presence of resonances the behaviour of an averaged system changes substantially due to the emergence of additional variables [3, 17].
Note that in the studied star-planet-satellite system the impact of the satellite is destabilising. Namely, in the system comprised of the star and the planet only the nutation angle does not change in time, implying that A = 0. After the satellite is added, we have A(I0, h0) = 0 for almost all initial conditions. Hence, our study gives a simple example of destabilising impact of the exo-Moon, showing that the conjecture about the stabilising effect of a satellite [26, 27] does not always hold true.
Conflict of interest
The author declares that she has no conflicts of interest.
References
[1] Andoyer, M.H., Cours de Mécaniquee Céleste: Vol. 1, Paris: Gauthier-Villars, 1923.
[2] Armstrong, J. C., Barnes, R., Domagal-Goldman, S., Breiner, J., Quinn, T. R., and Meadows, V. C., Effects of Extreme Obliquity Variations on the Habitability of Exoplanets, Astrobiology, 2014, vol. 14, no. 4, pp. 277-291.
[3] Arnol'd, V. I., Kozlov, V. V., and Nelshtadt, A. I., Mathematical Aspects of Classical and Celestial Mechanics, 3rd ed., Encyclopaedia Math. Sci., vol. 3, Berlin: Springer, 2006.
[4] Atobe, K. and Ida, S., Obliquity Evolution of Extrasolar Terrestrial Planets, Icarus, 2007, vol. 188, no. 1, pp. 1-17.
[5] Beletskii, V. V., Motion of an Artificial Satellite about Its Center of Mass, Jerusalem: Israel Program for Scientific Translations, 1966.
[6] Beletskii, V. V., Resonance Rotation of Celestial Bodies and Cassini's Laws, Celestial Mech., 1972, vol. 6, pp. 356-378.
[7] Beletskii, V. V., Satellite's Motion about Center of Mass in a Gravitational Field, Moscow: MGU, 1975 (Russian).
[8] Boué, G. and Laskar, J., A Collisionless Scenario for Uranus Tilting, Astrophys. J., 2010, vol. 712, no. 1, pp. L44-L47.
[9] Boué, G., Laskar, J., and Kuchynka, P., Speed Limit on Neptune Migration Imposed by Saturn Tilting, Astrophys. J, 2009, vol. 702, no. 1, pp. L19-L22.
[10] Brunini, A., Correction: Origin of the Obliquities of the Giant Planets in Mutual Interaction in the Early Solar System, Nature, 2006, vol. 443, no. 7114, p. 1013.
[11] Brunini, A., Origin of the Obliquities of the Giant Planets in Mutual Interaction in the Early Solar System, Nature, 2006, vol. 440, no. 7088, pp. 1163-1165.
[12] Correia, A. C. M. and Laskar, J., Long-Term Evolution of the Spin of Venus: 2. Numerical Simulations, Icarus, 2003, vol. 163, no. 1, pp. 24-45.
[13] Correia, A. C.M., Laskar, J., and Néron de Surgy, O., Long-Term Evolution of the Spin of Venus: 1. Theory, Icarus, 2003, vol. 163, no. 1, pp. 1-23.
[14] Dehant, V. and Mathews, P. M., Precession, Nutation and Wobble of the Earth, Cambridge: Cambridge Univ. Press, 2015.
[15] Ferreira, D., Marshall, J., O'Gorman, P. A., and Seager, S., Climate at High-Obliquity, Icarus, 2014, vol. 243, pp. 236-248.
[16] Heller, R., Leconte, J., and Barnes, R., Tidal Obliquity Evolution of Potentially Habitable Planets, Astron. Astrophys., 2011, vol. 528, A27, 16 pp.
[17] Henrard, J. and Lemaitre, A., A Second Fundamental Model for Resonance, Celestial Mech., 1983, vol. 30, pp. 197-218.
[18] Kilic, C., Raible, C.C., and Stocker, T. F., Multiple Climate States of Habitable Exoplanets: The Role of Obliquity and Irradiance, Astrophys. J., 2017, vol. 844, no. 2, Art. 147, 13 pp.
[19] Kinoshita, H., Theory of the Rotation of the Rigid Earth, Celestial Mech., 1977, vol. 15, no. 3, pp. 277-326.
[20] Kozai, Y., Secular Perturbations of Asteroids with High Inclination and Eccentricity, Astron. J., 1962, vol. 67, pp. 591-598.
[21] Krasilnikov, P. S., Applied Methods for the Study of Nonlinear Oscillations, Izhevsk: R&C Dynamics, Institute of Computer Science, 2015 (Russian).
[22] Krasilnikov, P. S. and Amelin, R. N., On the Precession of Saturn, Cosmic Research, 2018, vol. 56, no. 4, pp. 306-316; see also: Kosmicheskie Issledovaniya, 2018, vol. 56, no. 4, pp. 326-336.
[23] Krasilnikov, P. S. and Podvigina, O. M., On Evolution of the Planet's Obliquity in a Non-Resonant Planetary System, Vestn. Udmurtsk. Univ. Mat. Mekh. Komp. Nauki, 2018, vol. 28, no. 4, pp. 549564 (Russian).
[24] Krasil'nikov, P. S. and Zakharova, E. E., Nonresonant Rotations of a Satellite about Its Center of Mass in the Restricted N-Body Problem, Cosmic Research, 1993, vol. 31, no. 6, pp. 596-604; see also: Kosmicheskie Issledovaniya, 1993, vol. 31, no. 6, pp. 11-21.
[25] Laskar, J., Correia, A. M.C., Gastineau, M., Joutel, F., Levrard, B., and Robutel, P., Long Term Evolution and Chaotic Diffusion of the Insolation Quantities of Mars, Icarus, 2004, vol. 170, no. 2, pp. 343-364.
[26] Laskar, J., Joutel, F., and Robutel, P., Stabilization of the Earth's Obliquity by the Moon, Nature, 1993, vol. 361, no. 6413, pp. 615-617.
[27] Laskar, J., Joutel, F., and Boudin, F., Orbital, Precessional and Insolation Quantities for the Earth from -20 Myr to +10 Myr, Astron. Astroph, 1993, vol. 270, nos. 1-2, pp. 522-533.
[28] Laskar, J. and Robutel, P., The Chaotic Obliquity of the Planets, Nature, 1993, vol. 361, no. 6413, pp. 608-612.
[29] Lidov, M. L., The Evolution of Orbits of Artificial Satellites of Planets under the Action of Gravitational Perturbations of External Bodies, Planet. Space Sci, 1962, vol. 9, no. 10, pp. 719-759.
[30] Lidov, M.L. and Ziglin, S.L., The Analysis of Restricted Circular Twice-Averaged Three Body Problem in the Case of Close Orbits, Celestial Mech., 1974, vol. 9, no. 2, pp. 151-173.
[31] Lissauer, J. J., Barnes, J. W., and Chambers, J. E., Obliquity Variations of a Moonless Earth, Icarus, 2011, vol. 217, no. 1, pp. 77-87.
[32] Markeev, A. P. and Krasilnikov, P.S., On Motion of a Satellite Relative to the Center of Mass in the Elliptic Restricted Three Body Problem, Kosmicheskie Issledovaniya, 1981, vol. 19, no. 2, pp. 178-190 (Russian).
[33] Milankovitch, M., Canon of Insolation and the Ice-Age Problem, Jerusalem: Israel Program for Scientific Translations, 1969.
[34] Murray, C.D. and Dermott, S.F., Solar System Dynamics, Cambridge: Cambridge Univ. Press, 2000.
[35] Podvigina, O.M. and Krasilnikov, P.S., Evolution of Obliquity of an Exoplanet: A Non-Resonant Case, Icarus, 2020, vol. 335, Art. 113371.
[36] Podvigina, O. M. and Krasilnikov, P. S., Impact of a Moon on the Evolution of a Planet's Rotation Axis: A Non-Resonant Case, Celestial Mech. Dynam. Astronom., 2022, vol. 134, no. 3, Paper No. 21, 24 pp.
[37] Quarles, B., Barnes, J. W., Lissauer, J. J., and Chambers, J., Obliquity Evolution of the Potentially Habitable Exoplanet Kepler-62f, Astrobiology, 2020, vol. 20, no. 1, pp. 73-90.
[38] Saillenfest, M., Laskar, J., and Boue, G., Secular Spin-Axis Dynamics of Exoplanets, Astron. As-troph., 2019, vol. 623, A4, 21 pp.
[39] Shan, Y. and Li, G., Obliquity Variations of Habitable Zone Planets Kepler-62f and Kepler-186f, Astron. J., 2018, vol. 155, no. 6, Art. 237, 15 pp.
[40] Spiegel, D. S., Menou, K., and Scharf, C. A., Habitable Climates: The Influence of Obliquity, Astro-phys. J., 2009, vol. 691, no. 1, pp. 596-610.
[41] Touma, J. and Wisdom, J., The Chaotic Obliquity of Mars, Science, 1993, vol. 259, no. 5099, pp. 1294-1297.
