Influence of a second satellite on the rotational dynamics of an oblate moon

The gravitational influence of a second satellite on the rotation of an oblate moon is numerically examined. A simplified model, assuming the axis of rotation perpendicular to the (Keplerian) orbit plane, is derived. The differences between the two models, i.e. in the absence and presence of the second satellite, are investigated via bifurcation diagrams and by evolving compact sets of initial conditions in the phase space. It turns out that the presence of another satellite causes some trajectories, that were regular in its absence, to become chaotic. Moreover, the highly structured picture revealed by the bifurcation diagrams in dependence on the eccentricity of the oblate body’s orbit is destroyed when the gravitational influence is included, and the periodicities and critical curves are destroyed as well. For demonstrative purposes, focus is laid on parameters of the Saturn–Titan–Hyperion system, and on oblate satellites on low-eccentric orbits, i.e. e≈0.005\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e\approx 0.005$$\end{document}.


Introduction
Saturn's seventh moon, Hyperion (also known as Saturn VII), was discovered in the XIX century by Bond (1848) and Lassel (1848), but only due to Voyager 2 (Smith et al. 1982) and Cassini (Thomas 2010) missions it became apparent that it is the biggest known highly aspherical celestial body in the Solar System, with a highly elongated shape and dimensions 360 × 266 × 205 km. Since the rotational state of Hyperion was predicted to remain in the chaotic zone (Wisdom et al. 1984) based on the spin-orbit coupling theory (Goldreich and Peale 1966), further analyses and observations, regarding Hyperion as well as other Solar System satellites, were conducted on a regular basis. Hyperion's long-term observations were carried out twice in the post Voyager 2 era. In 1987, Klavetter (1989a, b) performed photometric R band observations over a timespan of more than 50 days, resulting in 38 high-quality data points. In 1999 and 2000, Devyatkin et al. (2002) conducted C (integral), B, V and R band observations. The objective of both analyses was to determine whether Hyperion's rotation is chaotic and to fit a solution of the equation of motion to the observations. To the best of the author's knowledge (Melnikov, priv. comm.) there were no other long-term observations that resulted in a lightcurve allowing the determination of Hyperion's rotational state (see also Strugnell and Taylor 1990;Dourneau 1993 for a list of earlier observations). Although, shortly after the Cassini 2005 passage a ground-based BV R photometry was conducted (Hicks et al. 2008), resulting in 6 nights of measurements (and additional 3 nights of R photometry alone) over a month-long period. Unfortunately, this data was greatly undersampled and period fitting procedures yielded several plausible solutions.
The theoretical and numerical treatment of the rotational dynamics of an oblate satellite have been performed widely. After the seminal paper of Wisdom et al. (1984), Boyd et al. (1994) applied the method of close returns to a sparse and short-term simulated observations of Hyperion's lightcurve. Black et al. (1995) performed numerical experiments using the full set of Euler equations to model long-term dynamical evolution. Beletskii et al. (1996) considered a number of models, including the gravitational, magnetic and tidal moments as well as rotation in gravitational field of two centers. A model with a tidal torque was examined analytically using Melnikov's integrals and assymptotic methods (Khan et al. 1998). The stability of resonances with application to the Solar System satellites was inferred based on a series expansion of the terms in the equation of rotational motion Chierchia 1998, 2000). The Lyapunov exponents and spectra were exhaustively examined for a number of satellites 1 (Shevchenko 2002;Shevchenko and Kouprianov 2002;Shevchenko 2003, 2005). A model of an oblate satellite with dissipation was used to examine the basins of attraction in case of low eccentricities, especially with application to the Moon (Celletti and Chierchia 2008). The dynamical stability was examined for all known satellites by Melnikov and Shevchenko (2010). Again the dynamical modeling using the full Euler equations was conducted by Harbison et al. (2011), who also analyzed the moments of inertia in light of the precessional period. Finally, Tarnopolski (2015a) argued that in order to extract a maximal Lyapunov exponent from the photometric lightcurve of Hyperion, at least one year of dense data is required.
The orbital dynamics of Hyperion in the Saturn-Titan-Hyperion system (see Table 1 for some physical parameters) have been exhaustively examined due to the interesting 4:3 mean motion resonance between Hyperion and Titan (Peale 1976;Taylor 1992;Stellmacher 1999;Rein et al. 2012). While the impact of Titan's gravitation on Hyperion's orbit has been established (Taylor et al. 1987) and the stability of the resonance has been considered in great detail (Colombo et al. 1974;Bevilacqua et al. 1980), introduction of the gravitational impact of a secondary body on the rotation of an oblate satellite has been done before for nearly spherical bodies such as Venus (Beletskii and Levin 1981) or low-eccentric orbits in the Pluto-Charon system (Correia et al. 2015). Herein, numerical integrations will be performed within the chaotic zone of the Saturn-Titan-Hyperion system with parameters ω 2 and e such that the perturbation techniques are not valid (Maciejewski 1995), which to the best of the author's knowlege has not yet been examined and hence is the aim of this work, which is general enough to be applicable to moons other than Hyperion. To focus attention, Hyperion's orbital period T 21.3 d Thomas et al. (2007) throughout the analysis the parameters are set to those of the Saturn-Titan-Hyperion system unless otherwise stated, but low-eccentricity and low-oblateness cases are also investigated for comparison. This paper is organized in the following manner. In Sect. 2 the rotational models in case of the absence and presence of a second satellite's gravitation are derived. In Sect. 3 the phase space is briefly described. Section 4 presents the methods used: the correlation dimension and its benchmark testing, and the bifurcation diagrams. The results are presented in Sect. 5, which is followed by discussion and conclusions gathered in Sect. 6. The computer algebra system mathematica ® is applied throughout this paper.

Rotational model of an oblate moon
The equation of rotational motion can be derived based on the following assumptions (Greiner 2010): 1. the orbit of the satellite around the planet is Keplerian with eccentricity e and major semi-axis a: where f H is the true anomaly given bẏ with M the mass of the planet and the overdot denotes differentation with respect to time; 2. in general, the physical model of the satellite is a triaxial ellipsoid; however, to simplify calculations, the satellite is simulated by a double dumbbell with four mass points 1-4 Fig. 1 Left Rotational model of an oblate moon. Right Geometry of the model including the orbital motion of a second satellite. S, T and H stand exemplary for Saturn, Titan and Hyperion (center of mass), respectively. See text for explanation of the remaining symbols (see Fig. 1) with equal mass m arranged in the orbital plane. The principal moments of inertia are A > B > C; 3. the satellite's spin axis is fixed and perpendicular to the orbit plane; the spin axis is aligned with the shortest physical axis, i.e. the one corresponding to the largest principal moment of inertia.
In case of Hyperion, the first assumption is not precisely valid, as it is well known that due to gravitational interaction with Titan the eccentricity of Hyperion oscillates from ∼ 0.08 to ∼ 0.12 with an 18.8-year period (Taylor et al. 1987). However, as the analysis herein is performed on a time span much shorter than this period (i.e., < 1 yr), the effect of this interaction will be negligible and as such is omitted (Black et al. 1995;Shevchenko and Kouprianov 2002). The second assumption, while might look like an oversimplification at first, does not affect the final equation of motion, which is the same as the one obtained directly from the Euler equations (Danby 1962; see also Appendix 1 for a remark on the moments of inertia in both models). The third assumption is justified for most satellites as the angular momentum is assumed to be constant with great accuracy. However, it should be noted that Wisdom et al. (1984) showed that the chaotic state is attitude unstable, and also the analysis of Voyager 2 images showed that the axis of rotation was far from being perpendicular to the orbital plane. Therefore, the models derived herein are a first approximation that will be expected to give initial insight into the dynamics of the satellite, and the parameters corresponding of the Saturn-Titan-Hyperion system are used for demonstrative reasons.
Defining the oblateness as ω 2 = 3(B−A) C , and choosing the units so that the orbital period T is equal to 2π and the major semi-axis a = 1 (which implies through Kepler's third law G M = 1, and that the orbital mean motion n = 1), eventually the equation of motion takes the formθ where time is dimensionless andθ is measured in units of n, with Eq. (1) for the orbit in the form and Eq. (2) for the true anomaly yieldṡ Moreover, transforming Eq. (3) so that f H is the independent variable leads to the famous Beletskii equation (Beletskii 1963), which was shown to be non-integrable (Maciejewski 1995).

Introducing a second satellite
Herein, a second satellite is assumed to revolve around the planet on a circular orbit, with radius a 0 , in the same plane as the oblate moon. Based on Eq. (5), the true anomaly depends linearly on time: From the triangle T SH (see Fig. 1) one obtains that the distance between the two satellites is equal to The angle α 1 is also required. Using again the triangle T SH one finds The angle x can be found by applying the law of cosines to the same triangle T SH, what gives Inserting Eq. (9) into (8) one arrives at the formula for α 1 . Finally (see Appendix 2), one obtains the following equation of motion: The initial conditions (ICs) for the true anomalies will be assumed throughout to be The backward differentiation formula (BDF) is employed for numerical integrations (Ascher and Petzold 1998).

Phase space properties
In this Section the structure of the phase space of the dynamical system given by Eqs.
(3)-(5) is briefly described. This will allow an insight into how does the gravitational interaction with the second satellite influence the oblate moon's rotation. The phase space is 3-dimensional: But f H is a regular, 2π-periodic function and does not carry much information. Moreover, in dimensionless units the orbital period of the oblate moon is also equal to 2π. Hence, a Poincaré surface of section, constructed by taking the values of (θ,θ) with a time step of 2π, i.e. employing stroboscopic variables, provides insight into the rotational dynamics. Furthermore, the rotation of the satellite by 180 • (i.e., θ → θ + π) gives an equivalent configuration, hence θ can be confined to the interval [0, π). Such surfaces of section are shown in Fig. 2.
As is common in Hamiltonian systems, the phase space is divided into regions occupied with chaotic trajectories, and regions of regular (periodic or quasiperiodic) motion. There are also motions called sticky orbits, when the trajectory initially behaves in a regular manner and diverges into the chaotic zone after some time. (See also Melnikov 2014 for the emergence of strange attractors when dissipation is introduced.) The phase space in Fig. 2a is dominated by a large chaotic sea, formed by merging the chaotic zones surrounding spin-orbit resonances from p=1:2 to p=2:1 when ω 2 increases (Wisdom et al. 1984). Quasiperiodic motions are represented by smooth curves and by closed curves, e.g. the ones surrounding the synchronous p = 1:1 resonance (see Table 1 in Black et al. 1995 for locations of the surviving resonances). When the IC is located near the boundary between regular motion and the chaotic sea, sticky motion occurs. A narrow chaotic zone is present also in Fig. 2b, c obtained for smaller values of the eccentricity and oblateness, respectively. In fact, every Solar System satellite has a chaotic period in its past (Spohn et al. 2014).

Correlation dimension
The algorithm and programme for computing the correlation dimension are briefly described in Sect. 4.1.1, which is followed by the description of methodology and discussion of the results of the benchmark testing in Sect. 4.1.2.

Algorithm
A fractal dimension (or, more precisely, the Hausdorff dimension; Hausdorff 1919;Theiler 1990;Ott 2002) is often measured with the correlation dimension, d C (Grassberger and Procaccia 1983;Grassberger 1986;Theiler 1990;Alligood et al. 2000;Ott 2002), which takes into account the local densities of the points in the examined dataset. For usual 1D, 2D or 3D cases the d C is equal to 1, 2 and 3, respectively. Typically, a fractional correlation dimension is obtained for fractals (Mandelbrot 1983).
The correlation dimension is defined as with the estimate for the correlation function C(R) as where the Heaviside step function Θ adds to C(R) only points x i in a distance smaller than R from x j and vice versa. The total number of points is denoted by N , and the usual Euclidean distance, ||.||, is employed. The limit in Eq. (11) is attained by using multiple values of R and fitting a straight line to the linear part of the obtained dependency. The correlation dimension is estimated as the slope of the linear regression. The computations in this work were performed using the python code from (Tarnopolski 2014), with a slight modification so that ln R, instead of R, is uniformly sampled with a step Δ ln R. Throughout this paper, when the d C is considered, N is set to be 10, 201 (see Sect. 5.2).

Benchmark testing
In order to verify the reliability of the correlation dimension algorithm, benchmark testing is performed on 128 realisations of uniform sampling with N points in each of the regions I and II defined as follows: region I is a unit square, and region II is a unit square without a circle of radius 1/4 placed at the center of the square. Next, the d C is computed as described in Sect. 4.1.1, with ln R ∈ [−8, −2] and Δ ln R = 0.5. The results are gathered in Table 2, from which it follows that the estimated d C is very close to the correct dimension expected for an Euclidean 2-dimensional space.
To validate the performance when clustering is introduced, regions I and II are uniformly sampled with 9000 points, and the remaining 1201 points are introduced in a circular region with radius equal to 1/8 (which is overlaid with points in the unit square), randomly chosen at each of the 128 realisations and lying entirely within region I or II. The results are gathered in Table 3. The correlation dimensions are close to the expected value of 2, but slightly smaller than they were when there was no clustering.
Finally, 128 realisations of uniform distributions of 2101 points within a unit square, which was overlaid with 8000 points distributed uniformly in a 0.25 × 0.25 square located in one corner of the unit one, were generated. The resulting mean d C was 1.937, and σ = 0.008.
Hence it was shown that clustering might result in a correlation dimension systematically lower than the expected one. Figure 3 shows the ln C(R) versus ln R relations for a uniformly sampled unit square and for the last numerical experiment on clustering.

Bifurcation diagrams
In dynamical systems theory, a bifurcation occurs when an infinitesimal change of a (nonlinear) parameter governing the system leads to a topological change in its behaviour. Generally speaking, at a bifurcation the stability of equilibria, periodic orbits or other invariant sets is changed. The theory of bifurcation is a vast field (Crawford 1991;Lichtenberg and Lieberman 1992;Baker and Gollub 1996;Alligood et al. 2000;Ott 2002;Kuznetsov 2004;Peitgen et al. 2004). Herein focus will be laid on the pitchfork bifurcations, that are present e.g. in the logistic map (May 1976;Feigenbaum 1979) and that constitute one of the routes to chaos. Let us consider a system of differential equations in the forṁ x = f (x; α), where α is a parameter, and assume that given x 0 as an IC, for α < α 1 the orbit is 1-periodic. A bifurcation at α = α 1 is a point where the trajectory begins to be 2-periodic and maintains its periodicity up to α = α 2 . Similarly, at α = α 2 a bifurcation occurs on each of the two branches, hence the orbit becomes 4-periodic. This scheme, called a period-doubling (pitchfork) bifurcation cascade, continues until at α = α ∞ < ∞ the orbit becomes chaotic. However, in the chaotic zone, α > α ∞ , windows of periodic motion with arbitrary period occur. E.g., when a 3-periodic trajectory emerges from the chaotic zone it also undergoes the period-doubling, hence produces orbits that are 6periodic, 12-periodic, and so on. A bifurcation diagram is a diagram illustrating this complex mechanism with respect to the nonlinear parameter α. Finally, bifurcations may also occur when α is decreasing (period-halving bifurcations) as well as when |α| is decreasing or increasing. Fig. 4 The phase space evolution after k revolutions of (from top to bottom): IC1 in case of the second satellite's absence, IC1 but including additional gravitational influence, IC2 in the absence of the second satellite, and IC2 in its presence. The sets IC1 and IC2 are depicted in red in each panel Usually, on the bifurcation diagrams there appear to be some curves running through the plot in the chaotic region. These are the so called critical curves (Peitgen et al. 2004) defined by x = f n (x 0 ; α).
IC1 was chosen so that it is located on the edge of the chaotic sea and the domain of quasiperiodic motion (according to Black et al. 1995, the 1:2 resonance is located at (π/2, 0.9); compare also with Fig. 2a), and IC2 was chosen so that it lies entirely within the chaotic zone (according to Wisdom et al. (1984), the 3:2 resonance does not exist).
The equations of motion (3) and (10) are solved numerically for every IC in the sets IC1 and IC2. The integration time is equal to 20 orbital periods of the oblate satellite, i.e. t max = 20 × 2π. After each revolution, starting from t = 0 (corresponding to the sets IC1 and IC2), the value of θ andθ is recorded, and the sets to which IC1 and IC2 evolved after k orbital periods are displayed in Fig. 4. Next, the correlation dimensions of these sets are computed according to Sect. 4.1. The linear regression was performed for ln C(R) vs. ln R relation in the region ln R ∈ (−7, ln R th ), where R th was chosen in each case so that the fit was done only in the linear part of the plot. The results, in graphical form, are displayed in Fig. 5. Because all orbits corresponding to IC1 and the absence of the second satellite are quasiperiodic, the d C , as expected, plateaus to a value near 1 (Fig. 5a). However, when its influence is taken into account, the d C initially tends to a slightly higher value, approximately 1.1, but then starts to rise suddenly (Fig. 5b). As can be seen in the second row of Fig. 4, some of the orbits appear to remain quasiperiodic when the other satellite's gravitation is switched on, but some become chaotic. On the other hand, based on the behaviour of d C in Fig. 5b, one might suspect that sticky chaos is also encountered. While the analysis of different types of rotation-periodic, quasiperiodic, and chaotic, including sticky chaos-is beyond the scope of this work, it is remarkable to note a clear impact of a second satellite on the rotation of an oblate moon: some ICs that would lead to quasiperiodic motion become chaotic.
The evolution of IC2 is very similar in case of both the absence and presence of the second satellite's gravity (third and fourth rows in Fig. 4). Also, the d C behaves in the same manner for both models, reaching a plateau of d C ≈ 1.75 after about 16 orbital periods, as shown in Fig. 5c, d. It is important to note that the value 1.75 cannot be a reliable estimate of the supposed fractality of the attained structure due to the following reasons: • both models given by Eqs. (3) and (10) are Hamiltonian and hence cannot posses a strange attractor (Greiner 2010) that could be characterized by a fractional correlation dimension; assymptotically any set of ICs leading to chaotic orbits should occupy a 2-dimensional subset in the phase space; • the number of points, N = 10, 201, used to estimate the d C is relatively small and hence might bias the outcome (Tarnopolski 2014); • local densities exceeding the average density (clustering) affect the correlation dimension such that it is lower than the dimension of the embedding space (see Sect. 4.1.2 and Tarnopolski 2015b).

Bifurcation diagrams
First, Eq. (3), describing the rotation in the absence of a second satellite, was integrated using the initial conditions (θ 0 ,θ 0 ) = (π/2, 3/2) in the time range t ∈ (0, 10 4 ), and the values ofθ were recorded every revolution, i.e. with a time step of 2π. To obtain the bifurcation diagrams in dependence on the oblateness, the eccentricity e was set to 0.1 and ω 2 was varied. Next, the same procedure was undertaken to obtain the bifurcation diagrams in dependence on the eccentricity, i.e. the oblateness was set to ω 2 = 0.79 and e was varied. The results are presented in Fig. 6a, c. For relatively small values of ω 2 and e, the magnifications in Fig. 6b, d show a very complex structure with alternating quasiperiodic and periodic orbits with a wide range of periods, e.g. a 4-period at ω 2 ≈ 0.08 or a 2-period at e ≈ 0.005 in Fig. 6b, d, respectively. Note that global chaos occurs at ω 2 ≈ 0.083, which corresponds to ω ≈ 0.29, what is in good agreement with the critical value obtained by Wisdom et al. (1984) using the resonance overlap criterion (Chirikov 1979;Lichtenberg and Lieberman 1992), ω RO = 1/ 2 + √ 14e ≈ 0.31. Hence the bifurcation diagrams confirm the applicability of this criterion to the rotation of an oblate moon.
In the same way the bifurcation diagrams were obtained using Eq. (10) including the gravitational influence of a second satellite. In this case the motion is governed by a third parameter, m 1 /M, in addition to ω 2 and e. The diagram in Fig. 7a (where m 1 /M was set to 0.00024) is qualitatively similar to the corresponding one in Fig. 6a, and the dependence on the ratio m 1 /M in Fig. 7b is mainly structureless. However, there is a significant difference when the oblateness and mass ratio are set to the values from Table 1, and the eccentricity e is varied. Note that the range of e in Fig. 7c is about four times smaller than in the corresponding Fig. 6c, as the computational complexity of the problem was much higher when Eq. (10) was applied. However, using a more sparse grid it was confirmed that the range ofθ also increases with the eccentricity (not shown). A remarkable difference becomes apparent in Fig. 7d: the image drawn by the bifurcation diagram is less structured compared to Fig. 6d, and the critical curves are much more tangled. Hence, the impact of a second satellite on the rotation is that the additional body destroys a number of periodic orbits that could occur for low eccentricity values, but for the parameters (i.e., oblateness and eccentricity) of the Saturn-Titan-Hyperion system, its dynamics should also be expected to remain in a chaotic state.
As was pointed out in Sect. 2.1, the model of planar rotation is not applicable to the real Saturn-Titan-Hyperion system as it is attitude unstable. Hence, the bifurcation diagrams are also computed in dependence of ω 2 and the ratio m 1 /M with e = 0.005, which is the mean eccentricity of all Solar System satellites. 2 Figure 8 shows that for a range of ω 2 the rotation is highly structured in the absence of the second satellite. Interestingly, a distinct 2periodic window emerges at ω 2 ≈ 0.79, which surprisingly coincides with the oblateness of Hyperion. Figure 9 shows that the picture becomes much more complex also for such small eccentricity when the gravitational influence of a second satellite is introduced. In particular, the 2-periodic window from Fig. 8 has disappeared altogether. Finally, when ω 2 was set to 0.04, the bifurcation diagrams revealed a highly regular structure (not shown). The aim of this paper was to investigate how does the gravitational interaction with a second satellite influence the rotational dynamics of an oblate moon. A simplified model was designed, resulting in the equation on motion given in Eq. (10), being basically a perturbation of the well known Eq.
(3). The derived equation of motion introduces a third parameter, the mass ratio m 1 /M, additional to the oblateness ω 2 and eccentricity e. To allow comparison, two sets of ICs distributed uniformly in a 0.1 × 0.1 square in the phase space were evolved, in case of the absence and presence of the additional source of gravitation. In case of the set IC1 [centered at (π/2, 0.55)] the difference between the two models is qualitative in nature: when the second satellite was absent all trajectories were quasiperiodic (first row in Fig. 4), as indicated by the d C = 1 in Fig. 5a. Interestingly, when its presence was taken into account, one could observe leaking of the orbits into the chaotic sea (second row in Fig. 4). This phenomenon manifests itself also through a higher d C attained in Fig. 5b. Hence, it turns out that an additional satellite has the ability to change quasiperiodic orbits into chaotic   ones, i.e. it enlarges the chaotic domain. On the other hand, when the set IC2, located in the center of the chaotic region [an 0.1 × 0.1 square centered at (π/2, 1.5)], was considered, no long term (assymptotic) differences could be observed (third and fourth rows in Fig. 4), and the correlation dimension for both models reached a plateau at d C ≈ 1.75 < 2 (Fig. 5c,  d), likely due to clustering. However, this is not that surprising, given that the gravitational influence under investigation was three orders of magnitude smaller than the planet's, and that the rotational model in absence of the second satellite is dominated by the chaotic zone, hence it would be highly unexpected for it to have the ability to change chaotic motion into a regular one. The bifurcation diagrams, especially interesting when m 1 /M and ω 2 were fixed and e was varied, when small values of the eccentricity (i.e., e < 0.03) are considered lead to a conclusion that the regular and highly structured picture (Fig. 6d) becomes much more messy, and the transition to chaos occurs for smaller eccentricities than in the case when additional gravitation source is neglected (Fig. 7d). This is consistent with the results of the other method (i.e., evolving the sets IC1 and IC2) in the sense that the second satellite   changes regular motion into chaotic. The differences in case when ω 2 was varied was not that much remarkable (Figs. 6a, 7a), and both models lead to chaotic motion when larger e are considered (Fig. 6c, 7c). The bifurcation diagram in dependence on the ratio m 1 /M was mostly structureless (Fig. 7b). Eventually, the destruction of regular rotation caused by the second satellite might be ascribed to the destruction of the invariant tori (Tabor 1989; see also Celletti and Chierchia 1998 and references therein). Finally, when e was set to the mean eccentricity of all Solar System satellites (i.e., e = 0.005), the highly structured bifurcation diagram displayed in Fig. 8 also got destroyed and became much more tangled, as shown in Fig. 9. To conclude, the derived simplified model of a second satellite's influence on rotational dynamics of an oblate satellite implies that: 1. the additional source of gravitation can change some regular orbits into chaotic ones, and 2. it destroys the regularity, particularly the periodicities and critical curves, in the bifurcation diagram for small eccentricities e < 0.03.

Introducing the second satellite
The torques coming from an additional satellite and acting on the oblate moon are given by where m 1 is the mass of the second satellite. The total torque acting on an oblate body is then