Free Access
Issue
A&A
Volume 580, August 2015
Article Number L14
Number of page(s) 7
Section Letters
DOI https://doi.org/10.1051/0004-6361/201526800
Published online 20 August 2015

© ESO, 2015

1. Introduction

Circumbinary bodies are objects that orbit around a more massive binary system. In the solar system, the small satellites of the Pluto-Charon system are the best example (e.g., Brozović et al. 2015). Planets orbiting two stars, often called circumbinary planets, have also been reported (e.g., Correia et al. 2005; Doyle et al. 2011; Welsh et al. 2012). Most of these bodies are close enough to the central binary to undergo tidal dissipation, which slowly modifies the rotation rate until it becomes close to the mean motion (e.g., MacDonald 1964; Correia et al. 2014).

For an object that moves in an eccentric orbit around a single massive body, the rotation rate can be captured in a half-integer commensurability with the mean motion, usually called a spin-orbit resonance (Colombo 1965; Goldreich & Peale 1966). Moreover, for very eccentric orbits or large axial asymmetries, the rotational libration width of the individual resonances may overlap, and the rotation becomes chaotic (Wisdom et al. 1984; Wisdom 1987). However, for nearly circular orbits, all these equilibria disappear and the only possibility for the spin is synchronous rotation (e.g., Correia & Laskar 2009).

When a third body is added to the classic two-body problem, the mutual gravitational perturbations introduce additional spin-orbit resonances at the perturbing frequency (Goldreich & Peale 1967; Correia & Robutel 2013). In the case of circumbinary bodies there is, in addition, a permanent misalignment of the long inertia axis of the rotating body from the radius vector pointing to each inner body (Fig. 1). The resulting torque on the rotating body’s figure induces some rotational libration that may give rise to some unexpected behaviors for the rotation rate (Showalter & Hamilton 2015). In this Letter we investigate all possibilities for the rotation of circumbinary bodies, and apply them to the spin evolution of the small satellites of the Pluto-Charon system.

2. Spin-orbit coupling

We consider a planar three-body hierarchical system composed of a central binary with masses m0 and m1, together with an external circumbinary companion with mass m (Fig. 1), where mm1m0. For the orbits we use Jacobi canonical coordinates, with rb being the position of m1 relative to m0 (inner orbit), and r the position of m relative to the center of mass of m0 and m1 (outer orbit). The body with mass m has principal moments of inertia ABC and angular velocity given by ω. In Appendix A we provide the full equations of motion for the spin. Since tidal dissipation usually damps the obliquity to zero (e.g., Hut 1980; Correia 2009), for simplicity, here we describe the motion for ω normal to the orbit. In absence of tides, the equation for the rotation angle, θ, is then given by (Eq. (A.5)) ¨θ=σ22[(1δ)(ar0)3sin2(θf0)+δ(ar1)3sin2(θf1)],\begin{eqnarray} \ddot \theta = -\frac{\sigma^2}{2} \left[ (1-\delta) \left(\frac{a}{\rp}\right)^3 \sin 2 (\theta - f_0) \!+\! \delta \left(\frac{a}{\rc}\right)^3 \sin 2 (\theta-f_1) \right] , \label{150616a} \end{eqnarray}(1)where a and n are the semi-major axis and the mean motion of the outer orbit, respectively, (ri,fi) are the radial and the angular coordinates of ri, r0=r+δrbandr1=r+(δ1)rb,δ=m1m0+m1,andσ=n3BAC·\begin{eqnarray} \label{150616b} &&\vrp = \vr + \delta \vrb \quad \mathrm{and} \quad \vrc = \vr + (\delta-1) \vrb , \\ \label{150616c} &&\delta = \frac{m_1}{m_0+m_1} , \quad \mathrm{and} \quad \sigma = n\sqrt{3\frac{B-A}{C}}\cdot \end{eqnarray}The parameter δ is the mass ratio of the inner binary, while σ is approximately the frequency of small-amplitude rotational librations in an unperturbed synchronous resonance.

thumbnail Fig. 1

Reference angles and coordinates. The system is coplanar and is composed of an inner binary with masses m0 and m1, and an outer body m that rotates with angular velocity ω.

We can express r0-3\hbox{$\rp^{-3}$}, sinf0, and cosf0 as r0-3=||r+δrb||-3=r-3(1+2δr·rb||r||2+δ2||rb||2||r||2)3/2,\begin{eqnarray} \rp^{-3} = ||\vr + \delta \vrb||^{-3} = r^{-3} \left( 1 + 2\delta \frac{\vr \cdot \vrb}{||\vr||^2} + \delta^2 \frac{||\vrb||^2}{||\vr||^2} \right)^{-3/2}\! , \label{150617z} \end{eqnarray}(4)and eif0=eif(1+δ||rb||||r||ei(fbf))(1+2δr·rb||r||2+δ2||rb||2||r||2)1/2·\begin{eqnarray} \ep^{\ii f_0} = \ep^{\ii f} \left( 1 + \delta \frac{||\vrb||}{||\vr||} \ep^{\ii (\fb-f)} \right) \left( 1 + 2\delta \frac{\vr \cdot \vrb}{||\vr||^2} + \delta^2 \frac{||\vrb||^2}{||\vr||^2} \right)^{-1/2} \cdot \label{150617y} \end{eqnarray}(5)Similar expressions can be found for r1-3\hbox{$\rc^{-3}$}, sinf1, and cosf1, if we replace δ by (δ − 1). Therefore, since | | rb | | < | | r | |, we can develop Eq. (1) in power series of rb/r. In addition, for the hierarchical problem the orbits are approximately ellipses. Thus, r and rb can also be expanded in terms of Hansen coefficients (e.g., Laskar & Boué 2010). We then rewrite Eq. (1) as ¨θ=σ22l,p,qβl,p,q(e,eb,δ,α)sin2(θpMqMb),\begin{eqnarray} \ddot\theta = -\frac{\sigma^2}{2} \sum_{l,p,q} \beta_{l,p,q} (e, \eb, \delta, \alpha) \sin2 (\theta - p M - q \Mb) , \label{eq:rot_gen} \end{eqnarray}(6)where e is the eccentricity; M is the mean anomaly; βl,p,qαl, with α = ab/a; and p, q are half-integers. Spin-orbit resonances occur whenever \hbox{$\dot\theta = p n + q \nb$}. In Appendix B we provide the complete expression for Eq. (6) to the first order in the eccentricities and to the second order in α.

3. Near circular orbits

For near circular orbits (e ≈ 0 and eb ≈ 0), we can neglect both eccentricities, thus ra and rbab. Retaining terms in α3 in Eq. (6) gives

¨γ=σ22[(1+54ρ2)sin2γ+38ρ2sin(2γ2φ)+358ρ2sin(2γ+2φ)+1516ρ3sin(2γφ)+3516ρ3sin(2γ+φ)\begin{eqnarray} \ddot \gamma &=& -\frac{\sigma^2}{2} \left[ ( 1 + \frac54\vo_2 ) \sin 2\gamma \right.\nonumber \\ && + ~\frac38\vo_2\sin (2\gamma - 2\phi ) +\frac{35}{8}\vo_2\sin (2\gamma +2\phi) \nonumber \\ && +~ \frac{15}{16}\vo_3\sin (2\gamma -\phi) + \frac{35}{16}\vo_3\sin (2\gamma +\phi) \nonumber \\ && \left. +~ \frac{5}{16}\vo_3\sin (2\gamma -3\phi ) + \frac{105}{16}\vo_3\sin (2\gamma +3\phi) \right] , \label{150616z} \end{eqnarray}(7)where γ = θM, φ = MbM, ρ2=δ(1δ)α2,andρ3=δ(1δ)(12δ)α3.\begin{eqnarray} \vo_2 = \delta(1-\delta) \, \alpha^2 , \quad \mathrm{and} \quad \vo_3 = \delta(1-\delta)(1-2\delta)\, \alpha^3 . \label{150616y} \end{eqnarray}(8)We see that we have several islands of rotational libration, γ̇=θ̇n=0,±ν,±ν/2,±3ν/2,withν=φ̇=nbn.\begin{eqnarray} \dot\gamma=\dot\theta-n=0,\pm \nu, \pm\nu/2, \pm 3\nu/2 , \quad \mathrm{with} \quad \nu = \dot \phi = \nb - n . \label{150616x} \end{eqnarray}(9)Each term individually behaves like a pendulum where the rotation can be trapped. Therefore, together with the classic synchronous equilibrium at \hbox{$\dot\theta=n$}, there are additional possibilities for the spin at the super- and sub-synchronous resonances \hbox{$\dot\theta=n + k\nu/2$}, with k ∈Z.

Terms with k ≠ 0 depend on σ2, δ, and α2, that is, on the axial asymmetry of the body, on the mass ratio of the central binary, and on the semi-major axis ratio. Thus, for circumbinary bodies far from the inner binary (α ≪ 1), the amplitude of these higher order terms decreases very quickly to zero. For δ ≈ 0 (i.e., m1m0) we have ρ2 = ρ3 = 0, thus only the synchronous resonance subsists, as in the classic circular two-body problem (e.g., Goldreich & Peale 1966). The effect of the binary mass ratio is maximized for δ = 1/2 (i.e., m0 = m1), but only for the second order resonances (terms in ρ2), since third order resonances also vanish for equal masses (ρ3 = 0). Interestingly, the size of the sub-synchronous resonance with k< 0 is always larger than the symmetrical super-synchronous resonance with k> 0.

3.1. Chaotic motion

When the libration amplitudes of some individual resonant islands overlap, the rotational motion can be chaotic (Chirikov 1979; Morbidelli 2002). This means that the rotation exhibits random variations in short periods of time. This phenomenon was first described for the rotation of Hyperion, which is chaotic because its orbit is eccentric (Wisdom et al. 1984).

In the circular approximation, each individual resonance is placed at \hbox{$\dot\gamma = k \nu/2$}, with libration width σβk\hbox{$\!\sqrt{\beta_k}$}, where βk = ∑ lβl,(1 − k)/2,k/2. Overlap between two resonances with k1<k2 then occurs whenever k1ν/2+σβk1k2ν/2σβk2.\begin{eqnarray} k_1 \nu/2 + \sigma \sqrt{\beta_{k_1}} \gtrsim k_2 \nu/2 - \sigma \sqrt{\beta_{k_2}} . \label{150618a} \end{eqnarray}(10)The synchronous resonance (k = 0) has the largest width, so overlap is more likely for this resonance. The third order resonance with k = −1 is the nearest resonance with larger amplitude, thus chaos sets in for σ(3516ρ3+1+54ρ2)σν2=12(nbn).\begin{eqnarray} \sigma \left( \!\sqrt{\frac{35}{16} \vo_3} + \sqrt{1+\frac54 \vo_2} \right) \approx \sigma \gtrsim \frac{\nu}{2} = \frac12 \left( \nb - n \right) . \label{150618c} \end{eqnarray}(11)Rotation in the chaotic zone is also attitude unstable (Wisdom et al. 1984; Wisdom 1987). A non-zero obliquity introduces additional resonant terms to the rotation (Eq. (7)), increasing the chances of overlapping (e.g., Correia & Laskar 2010). Therefore, in order to correctly account for the chaotic behavior, we need to integrate the full equations of motion (Appendix A).

3.2. Global dynamics

thumbnail Fig. 2

Stability analysis of the rotation rate close to the synchronization for δ = 0.1, α ∈ [ 0.25:0.5 ], and (BA) /C ∈ [ 0:0.9 ]. We adopt three different initial obliquity values ψ = 1°, 30°, and 60°. The color index indicates the proportion of chaotic orbits inside the studied domain: from dark blue for fully regular to red for entirely chaotic. The dashed line (Eq. (11)) roughly delimits the zone where chaotic motion is possible. The dots identify two representative cases that we detail in Fig. 3.

The general problem of the spin-orbit coupling for near circular circumbinary bodies can be reduced to the analysis of \hbox{$\dot\gamma/n$}, which depends on the libration frequency σ/n, the mass ratio δ, and the semi-major axis ratio α (Eq. (7)). However, since δ ∈ [ 0,1/2 ] and it is a well-known parameter for a specific binary system, the global dynamics is approximately controlled by the only two free parameters: α (or nb/n) and σ/n (or (BA) /C).

Rescaling the time with n = 1, we can perform a stability analysis of \hbox{$\dot\gamma$} in the plane (α,σ) to quickly identify the rotational regime for any system of quasi-circular circumbinay bodies that is near the synchronous equilibrium. If isolated, the half-width of the synchronous resonant island in the direction of \hbox{$\dot\gamma$} is equal to σ. When perturbed, chaos is expected for σ ~ ν (Eq. (11)). Thus, for a given (α,σ) pair, we fix the initial value of γ at the libration center of the synchronous resonance, and select 400 equispaced initial values of \hbox{$\dot\gamma \in [-\nu{:}\nu]$}. The corresponding solutions are integrated using the equations for the three-body problem together with Eqs. (A.1) and (A.2) for the rotational motion (without tidal dissipation). The dynamical nature (stable/unstable) is deduced from frequency analysis (Laskar 1990, 1993), which gives the fraction of chaotic trajectories.

In Fig. 2 we show the general results for δ = 0.1 using three different initial obliquity values ψ = 1°, 30°, and 60°. All the remaining variables are initially set to zero. The color index indicates the proportion of chaotic orbits inside the studied domain: from dark blue for fully regular to red for entirely chaotic. Depending on the value of α and σ, the rotation can present different behaviors, ranging from non-synchronous equilibria to chaotic motion. Figure 2 remains almost unchanged for δ> 0.1 since the global dynamics is not very sensitive to δ. The only exception would be for δ ≪ 0.1 since in that case the system would behave like the classic circular two-body problem where the rotation can only be synchronous.

A common way of visualizing and understanding the different regimes is to use frequency map analysis (Laskar 1993). We selected two representative pairs (α,σ) and plotted the corresponding cross-section maps for the three initial obliquities in Fig. 3. We are able to identify the stable and the chaotic regions very easily. We observe that stable islands can exist in the middle of chaos, implying a large diversity of possibilities for the final rotation. In the chaotic regions for the rotation, the obliquity is also chaotic and can vary between and 180°. For high obliquities, the size of the stable synchronous island shrinks and the chaotic zone is extended. Therefore, depending on the obliquity, the rotation alternates between more or less chaotic. In times of small obliquity, the rotation can be trapped in a stable resonant island, preventing the obliquity from increasing again, and the spin can become stable (see Fig. 4).

thumbnail Fig. 3

Cross sections representing the spin dynamics of Hydra (top) and Nix (bottom) at different initial obliquities ψ = 1°, 30°, and 60°. The dots show the main frequency of γ(t), denoted by ξ, for different initial values of \hbox{$\dot \gamma(0)$}. It is obtained by numerical integration of the three-body problem together with Eqs. (A.1) and (A.2) for the rotational motion (without dissipation). θ is initially set equal to f so the sections go through the middle of the synchronous resonance (\hbox{$\dot{\gamma}=0$}). All the remaining variables are initially set to zero. The dots are blue for stable rotation (libration or circulation), and red for chaotic motion. The distinction between stable and unstable trajectories is based on the estimate of the second derivative of ξ with respect to the initial condition \hbox{$\dot\gamma(0)$} (see Laskar 1993). The plateaus correspond to the resonance crossings (ξ/ν = k/2). The red dots randomly distributed between the plateaus indicate the overlapping of the associated spin-orbit resonances, which generates chaotic motion.

For small values of α and/or σ, the individual resonances (Eq. (7)), associated with the plateaus in Fig. 3, are well separated. Thus, for dissipative systems, the rotation can be captured in individual spin-orbit resonances and stay there. For rotation rates increasing from lower values, the sub-synchronous resonance with \hbox{$\dot\theta=n-\nu = 2n-\nb$} has the largest amplitude, so this is the most likely spin-orbit resonance to occur. For rotation rates decreasing from higher values, the super-synchronous resonance with \hbox{$\dot\theta=n+\nu=\nb$} becomes the most likely possibility. This case is very interesting as it corresponds to a synchronous rotation with the inner orbit period. However, since the amplitude of this resonance is small, the rotation can easily escape it and subsequently evolve into the classic synchronous rotation with the outer orbit (\hbox{$\dot\theta = n$}).

The amplitude of the resonant terms increases with α and σ. At some point, the individual libration islands merge and chaotic motion can be expected (see Sect. 3.1). The transition between the two regimes is roughly given by the dashed curved, obtained with expression (11). Near the transition regime, chaotic rotation can be observed around the separatrix of the synchronous resonance, but the motion is still regular nearer the center (Fig. 3, left). Therefore, when the rotation of a body is evolving by tidal effect, the rotation becomes chaotic when approaching the synchronous resonance. However, depending on the strength and on the geometry of the tidal torque, it is still possible to find a pathway to the synchronous rotation. The wandering in the chaotic region may also provide a path into another stable super- or sub-synchronous resonance (see Fig. 4).

For α> 0.4 and σ ~ 1, there is a large overlap between several resonant terms, in particular for the negative ones (k< 0). As a consequence, there is a large chaotic zone for the spin (Fig. 3, bottom). Some of the individual resonances may still subsist, including the synchronous one, but they have small stable widths and they can only be reached at periods of small obliquity. The chaotic motion is maximized for σ ≈ 1, which corresponds to (BA) /C ≈ 0.35, because secondary resonances generate strong instabilities inside the synchronous island (Robutel et al. 2012). Finally, for α ≈ 0.49 the system reaches the 3/1 mean motion resonance, which introduces additional forcing to the rotation, and all trajectories become chaotic.

4. Application to the Pluto-Charon system

In 1978, regular series of observations of Pluto showed that the images were consistently elongated revealing the presence of Charon (Christy & Harrington 1978). The orbital parameters hence determined have shown that the two bodies evolved in an almost circular orbit with a 6.4-day period (e.g., Tholen & Buie 1997). Charon has an important fraction of the total mass (about 11%), and therefore the system is considered a binary planet rather than a planet and a moon. Later, it was found that four tiny satellites move around the barycenter of the Pluto-Charon system (Weaver et al. 2006; Brozović et al. 2015), also in nearly circular and coplanar orbits (see Table 1).

Table 1

Mean orbital and physical parameters for Pluto’s moons (Brozović et al. 2015).

Pluto’s brightness also varies with a period of 6.4 days (e.g., Tholen & Tedesco 1994). Since Charon is too dim to account for the amplitude of the variation, this period has been identified as the rotation of Pluto. Therefore, at present, the spin of Pluto is synchronous with the orbit of Charon, a configuration acquired from the action of tidal torques raised on Pluto by Charon (e.g., Farinella et al. 1979; Cheng et al. 2014). Tidal torques raised on Charon and on the remaining satellites by Pluto are even stronger, so all these bodies are presumably also tidally evolved.

Showalter & Hamilton (2015) have measured the brightness variations of the small satellites. They concluded that Hydra and Nix show no obvious pattern, suggesting that their rotation is chaotic. Assuming that these two satellites are uniform triaxial ellipsoids, they additionally estimate their figure semi-axis ratios, which allows us to compute B/C ≈ 0.92 and A/C ≈ 0.44 for Hydra, and B/C ≈ 0.94 and A/C ≈ 0.31 for Nix.

In Fig. 2 we already show the global dynamics for the Pluto-Charon system, for which δ = 0.1085 ≈ 0.1 (Brozović et al. 2015). Therefore, we plot vertical lines at the α values corresponding to the small satellites in the system (Table 1). Moreover, since we have an estimation for the axial asymmetries of Hydra and Nix, we also show a dot for their (BA) /C values in Fig. 2. These two satellites are in a different dynamical regime: the rotation of Hydra can present a stable spin-orbit coupling (Fig. 3, top), while the rotation of Nix is most likely chaotic (Fig. 3, bottom), as discussed in Sect. 3.2. The spin dynamics of Kerberos is probably similar to that of Hydra, while the rotation of Styx can be even more chaotic than that of Nix, depending on their (BA) /C values.

To test the reliability of the dynamical picture described in the previous section, we now add the effect from tides to our model, so that we can follow the long-term evolution of the spin of these satellites. For simplicity, we adopt a constant time-lag linear model, whose contribution to the spin is given by expression (A.13). In Fig. 4 we show some examples for the final evolution of Hydra’s spin with tidal dissipation, starting with slightly different initial values of θ. We adopt an initial retrograde rotation of 4.4 days (\hbox{$\dot\gamma/\nu=-2.3$}), and 30° for the initial obliquity, to force the rotation to cross the large amplitude sub-synchronous resonances. An initial retrograde rotation for such small bodies is as likely as a prograde one (e.g., Dones & Tremaine 1993). For the uncertain parameters we used R = 45 km, k2/Q = 10-4, and C/(mR2) = 0.4.

In one example (Fig. 4a) the rotation of Hydra is trapped in the spin-orbit resonance with \hbox{$\dot\gamma = - 3\nu/2$}, while in the others the rotation reaches the chaotic zone. However, after some wandering in this zone some simulations are able to find a path into a stable spin-orbit resonance with k> 0 (Fig. 4b) or with k< 0 (Fig. 4c). Interestingly, while for resonances with k> 0 the spin axis stabilizes near zero degrees, for resonances with k< 0, the spin axis stabilizes with a high obliquity value. The last example (Fig. 4d) remained chaotic for the length of the integration and therefore can represent the observed present state.

Over 10 Myr, most of our simulations remained chaotic, but several captures in stable non-synchronous resonances also occurred. The final scenario depends on the initial conditions, and also on the tidal model. A constant-Q model would prevent any capture in resonance (Goldreich & Peale 1966), while a visco-elastic model would increase the chances of capture (e.g., Makarov 2012; Correia et al. 2014). Strange attractors can also exist in the chaotic zone, which may prevent the spin from stabilizing (e.g., Batygin & Morbidelli 2011). However, the global dynamics described in Sect. 3.2 is very robust and reliable, since it does not depend on dissipative forces (Fig. 2).

5. Discussion

This work was motivated by the recent observations on the rotation of Hydra and Nix (Showalter & Hamilton 2015), which appear to be chaotic. Our model confirms that chaotic rotation is a likely scenario for both satellites, but stable spin-orbit coupling could also be possible, in particular for Hydra. This model assumes a hierarchical three-body system with coplanar orbits. However, we also integrated the spins of Hydra and Nix using the ephemerides for the full system provided by Brozović et al. (2015) and SPICE routines (Acton 1996), and we have found no evidence of any substantial differences.

In our study, we conclude that stable spin-orbit coupling is also a plausible scenario for near circular circumbinary bodies with small α and/or small σ (Fig. 2). Equilibrium rotation occurs for \hbox{$\dot\theta=n + k\nu/2$}, with ν = nbn and k ∈Z. The largest amplitude non-synchronous resonance corresponds to the sub-synchronous resonance at nν = 2nnb. Bodies captured in this resonance present retrograde rotation if nb/n> 2. This condition is verified for the small satellites of the Pluto-Charon system (Table 1), and likely for any circumbinary system, since large orbital

instabilities are expected for period ratios below the 2/1 mean motion resonance. Therefore, bodies trapped in a non-synchronous resonance are also likely to present retrograde rotation.

Lately, many planets have been detected around binary stars (e.g., Welsh et al. 2012). So far, all these planets are gaseous giants, for which the axial asymmetry is very low, for instance (BA) /C ~ 10-7 for Jupiter (Jacobson 2001). Therefore, although many of these planets are close enough to their stars to undergo tidal dissipation, spin-orbit coupling is very unlikely. However, for smaller mass Earth-like circumbinary planets these states are possible since (BA) /C ~ 10-5 (e.g., Yoder 1995). In particular, non-synchronous rotation is possible, which is an important point to take into account in future habitability studies (e.g., Selsis et al. 2007). Unlike the Pluto-Charon system, circumbinary exoplanets usually present eccentric orbits. As a consequence, the number of spin-orbit resonances drastically increases (Eq. (B.3)). Overlap of the different contributions is then easier, so chaotic rotation can also be more likely in this case.

Online material

thumbnail Fig. 4

Possible evolutions for the spin of Hydra. We show the rotation \hbox{$\dot\gamma/\nu$} (in red) and the cosine of the obliquity (in blue). We numerically integrated the equations of the three-body problem together with Eqs. (A.1), (A.2), and (A.13) for the rotational motion. We adopt an initial retrograde rotation of 4.4 days with different initial values of θ, and an initial obliquity of 30°.

Appendix A: Equations of motion

We let (\hbox{$\vi, \vj, \vk$}) be a non-inertial frame attached to the body’s principal inertial axes, with inertia tensor I = diag(A,B,C). In this frame, ω = (ωa,ωb,ωc) and the angular momentum L = (Aωa,Bωb,Cωc). The equations of motion are given by I·ω˙+ω×L=T0+T1,\appendix \setcounter{section}{1} \begin{eqnarray} \vv{I} \cdot \dot{\vw} + \vw \times \vL = \vv{T}_0 + \vv{T}_1 , \label{150626a} \end{eqnarray}(A.1)where Ti is the gravitational torque on the body’s figure. For second order interactions, we have (e.g., Goldstein 1950) Ti=3Gmiri5[(BA)yiri×+(CA)ziri×],\appendix \setcounter{section}{1} \begin{eqnarray} \vv{T}_i = \frac{3 G m_i}{r_i^5} \left[ (B-A) \, y_i \, \vr_i \times \vj + (C-A) \, z_i \, \vr_i \times \vk \right] , \label{050104b} \end{eqnarray}(A.2)where ri has coordinates (xi,yi,zi) in the body’s frame. Thus, projecting Eq. (A.1) over each axis (\hbox{$\vi, \vj, \vk$}) gives ω̇a=CBA(3Gm0r05y0z0+3Gm1r15y1z1ωbωc),ω̇b=ACB(3Gm0r05z0x0+3Gm1r15z1x1ωcωa),ω̇c=BAC(3Gm0r05x0y0+3Gm1r15x1y1ωaωb).\appendix \setcounter{section}{1} \begin{eqnarray} \label{150703a} &&\dot \omega_a = \frac{C-B}{A} \left( \frac{3 G m_0}{r_0^5} \, y_0 \, z_0 + \frac{3 G m_1}{r_1^5} \, y_1 \, z_1 - \w_b \, \w_c \right) , \\ \label{150703b} &&\dot \omega_b = \frac{A-C}{B} \left( \frac{3 G m_0}{r_0^5} \, z_0 \, x_0 + \frac{3 G m_1}{r_1^5} \, z_1 \, x_1 - \w_c \, \w_a \right) , \\ \label{150626b} &&\dot \omega_c = \frac{B-A}{C} \left( \frac{3 G m_0}{r_0^5} \, x_0 \, y_0 + \frac{3 G m_1}{r_1^5} \, x_1 \, y_1 - \w_a \, \w_b \right) . \end{eqnarray}To solve these equations, a set of generalized coordinates to specify the orientation of the satellite must be chosen. We adopt the modified Euler angles (θ,ϕ,ψ) as defined in Wisdom et al. (1984). Starting the ĉ axis coincident with the normal to the orbital plane, and â along an inertial direction, we first rotate the body about the ĉ axis by an angle θ, then we rotate about â by an angle ϕ, and finally we rotate about the \hbox{$\vj$} by an angle ψ. Then, xi/ri=cos(θfi)cosψsin(θfi)sinψsinϕ,yi/ri=sin(θfi)cosϕ,zi/ri=cos(θfi)sinψ+sin(θfi)cosψsinϕ,\appendix \setcounter{section}{1} \begin{eqnarray} \label{150703c}&& x_i/r_i = \cos(\theta-f_i) \cos\psi - \sin(\theta-f_i) \sin\psi \sin\vp , \\ \label{150703d} && y_i/r_i = - \sin(\theta-f_i) \cos\vp , \\ \label{150626x} && z_i/r_i = \cos(\theta-f_i) \sin\psi + \sin(\theta-f_i) \cos\psi \sin\vp , \end{eqnarray}and θ̇=(ωccosψωasinψ)/cosϕ,ϕ̇=ωcsinψ+ωacosψ,ψ̇=ωb(ωccosψωasinψ)tanϕ.\appendix \setcounter{section}{1} \begin{eqnarray} \label{150627b} && \dot\theta = (\omega_c \cos\psi - \omega_a \sin\psi)/\cos\vp , \\ \label{150703e} && \dot\vp = \omega_c \sin\psi + \omega_a \cos\psi , \\ \label{150703f} && \dot\psi = \omega_b - (\omega_c \cos\psi - \omega_a \sin\psi) \tan\vp . \end{eqnarray}The full equations of motion for the spin are thus described by the set of variables (ωa,ωb,ωc,θ,ϕ,ψ), whose derivatives are given by Eqs. (A.3)(A.5) together with (A.9)(A.11). The angle ε between the ĉ axis and the normal to the orbit, usually called obliquity, can be obtained from cosε=cosϕcosψ.\appendix \setcounter{section}{1} \begin{eqnarray} \cos \ve = \cos \vp \cos \psi . \label{150715a} \end{eqnarray}(A.12)For the dissipative tidal torque, we adopt a constant time-lag linear model, whose contribution to the spin (Eq. (A.1)) is given by (e.g., Mignard 1979) Td=3k2GR5Δti=0,1mi2ri8[(ri·ω)riri2ω+ri×ri˙],\appendix \setcounter{section}{1} \begin{eqnarray} \vv{T}_d = 3 k_2 G R^5 \Delta t \sum_{i\,=\,0,\,1} \frac{m_i^2}{r_i^8} \left[ ( \vr_i \cdot \vw ) \, \vr_i - r_i^2 \vw + \vr_i \times \dot{\vr_i} \right] , \label{150626c} \end{eqnarray}(A.13)

where k2 is the Love number, Δt is the time lag, Q-1nΔt is the dissipation factor, and R is the radius of the rotating body. Also note that ri˙\hbox{$\dot{\vr_i}$} is the derivative of ri in an inertial reference frame.

Appendix B: Series with eccentricity

Equation (1) for the rotational motion can be expanded in power series of the eccentricities and semi-major axis ratio α as given by the general expression (6) ¨θ=σ22Γ,\appendix \setcounter{section}{2} \begin{eqnarray} \ddot \theta = -\frac{\sigma^2}{2} \Gamma , \label{150616az} \end{eqnarray}(B.1)with Γ=l,p,qβl,p,q(e,eb,δ,α)sin2(θpMqMb).\appendix \setcounter{section}{2} \begin{eqnarray} \Gamma = \sum_{l,p,q} \beta_{l,p,q} \left(e, \eb, \delta, \alpha\right) \sin2 \left(\theta - p M - q \Mb\right) . \label{150621a} \end{eqnarray}(B.2)When we truncate the series to the first order in the eccentricities and to the second order in α, we get Γ=(1+54ρ2)sin(2θ2M)+38ρ2sin(2θ2Mb)+358ρ2sin(2θ4M+2Mb)+(72+458ρ2)esin(2θ3M)(1258ρ2)esin(2θM)10516ρ2esin(2θ3M+2Mb)+45516ρ2esin(2θ5M+2Mb)+1516ρ2esin(2θ+M2Mb)+1516ρ2esin(2θM2Mb)+38ρ2ebsin(2θ3Mb)98ρ2ebsin(2θMb)1058ρ2ebsin(2θ4M+Mb)+358ρ2ebsin(2θ4M+3Mb)54ρ2ebsin(2θ2MMb)\appendix \setcounter{section}{2} \begin{eqnarray} \Gamma &=& \left( 1 + \frac54\vo_2 \right) \sin \left(2 \theta - 2 M\right) \nonumber \\ && + \frac38\vo_2\sin \left(2 \theta - 2 \Mb \right) \nonumber \\ && +\frac{35}{8}\vo_2\sin \left(2\theta - 4 M + 2 \Mb\right) \nonumber \\ && + \left(\frac72 + \frac{45}{8} \vo_2 \right) \, e \sin\left( 2\theta -3M\right) \nonumber \\ && - \left(\frac12 - \frac{5}{8} \vo_2 \right) \, e \sin\left( 2\theta -M\right) \nonumber \\ && -\frac{105}{16} \vo_2 \, e \sin\left(2\theta - 3M + 2\Mb\right) \nonumber \\ && +\frac{455}{16} \vo_2 \, e \sin\left(2\theta - 5M + 2\Mb\right) \nonumber \\ && +\frac{15}{16} \vo_2 \, e \sin\left(2\theta + M- 2\Mb\right) \nonumber \\ && +\frac{15}{16} \vo_2 \, e \sin\left(2\theta - M- 2\Mb\right) \nonumber \\ && +\frac{3}{8} \vo_2 \, \eb \sin\left(2\theta - 3\Mb\right) \nonumber\\ && -\frac{9}{8} \vo_2 \, \eb \sin\left(2\theta - \Mb\right) \nonumber \\ && -\frac{105}{8} \vo_2 \, \eb \sin\left(2\theta - 4M + \Mb\right) \nonumber \\ && +\frac{35}{8} \vo_2 \, \eb \sin\left(2\theta - 4M + 3\Mb\right) \nonumber \\ && -\frac{5}{4} \vo_2 \, \eb \sin\left(2\theta - 2M -\Mb\right) \nonumber \\ && -\frac{5}{4} \vo_2 \, \eb \sin\left(2\theta - 2M +\Mb\right) . \label{150616w} \end{eqnarray}(B.3)

Acknowledgments

We thank Doug Hamilton and Mark Showalter for helpful suggestions. We acknowledge support from the “conseil scientifique” of the Observatory of Paris and CIDMA strategic project UID/MAT/04106/2013.

References

  1. Acton, C. H. 1996, Planet. Space Sci., 44, 65 [NASA ADS] [CrossRef] [Google Scholar]
  2. Batygin, K., & Morbidelli, A. 2011, Cel. Mech. Dyn. Astron., 111, 219 [NASA ADS] [CrossRef] [Google Scholar]
  3. Brozović, M., Showalter, M. R., Jacobson, R. A., & Buie, M. W. 2015, Icarus, 246, 317 [NASA ADS] [CrossRef] [Google Scholar]
  4. Cheng, W. H., Lee, M. H., & Peale, S. J. 2014, Icarus, 233, 242 [NASA ADS] [CrossRef] [Google Scholar]
  5. Chirikov, B. V. 1979, Phys. Rep., 52, 263 [Google Scholar]
  6. Christy, J. W., & Harrington, R. S. 1978, AJ, 83, 1005 [NASA ADS] [CrossRef] [Google Scholar]
  7. Colombo, G. 1965, Nature, 208, 575 [NASA ADS] [CrossRef] [Google Scholar]
  8. Correia, A. C. M. 2009, ApJ, 704, L1 [NASA ADS] [CrossRef] [Google Scholar]
  9. Correia, A. C. M., & Laskar, J. 2009, Icarus, 201, 1 [NASA ADS] [CrossRef] [Google Scholar]
  10. Correia, A. C. M., & Laskar, J. 2010, Icarus, 205, 338 [NASA ADS] [CrossRef] [Google Scholar]
  11. Correia, A. C. M., & Robutel, P. 2013, ApJ, 779, 20 [NASA ADS] [CrossRef] [Google Scholar]
  12. Correia, A. C. M., Udry, S., Mayor, M., et al. 2005, A&A, 440, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Dones, L., & Tremaine, S. 1993, Icarus, 103, 67 [NASA ADS] [CrossRef] [Google Scholar]
  15. Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
  16. Farinella, P., Milani, A., Nobili, A. M., & Valsecchi, G. B. 1979, Moon Planets, 20, 415 [NASA ADS] [CrossRef] [Google Scholar]
  17. Goldreich, P., & Peale, S. 1966, AJ, 71, 425 [NASA ADS] [CrossRef] [Google Scholar]
  18. Goldreich, P., & Peale, S. 1967, AJ, 72, 662 [NASA ADS] [CrossRef] [Google Scholar]
  19. Goldstein, H. 1950, Classical mechanics (Reading: Addison-Wesley) [Google Scholar]
  20. Hut, P. 1980, A&A, 92, 167 [NASA ADS] [Google Scholar]
  21. Jacobson, R. A. 2001, BAAS, 33, 1039 [NASA ADS] [Google Scholar]
  22. Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
  23. Laskar, J. 1993, Physica D Nonlinear Phenomena, 67, 257 [Google Scholar]
  24. Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. MacDonald, G. J. F. 1964, Rev. Geophys., 2, 467 [NASA ADS] [CrossRef] [Google Scholar]
  26. Makarov, V. V. 2012, ApJ, 752, 73 [NASA ADS] [CrossRef] [Google Scholar]
  27. Mignard, F. 1979, Moon Planets, 20, 30 [Google Scholar]
  28. Morbidelli, A. 2002, Modern celestial mechanics: aspects of solar system dynamics (London: Taylor & Francis) [Google Scholar]
  29. Robutel, P., Rambaux, N., & El Moutamid, M. 2012, Cel. Mech. Dyn. Astron., 113, 1 [NASA ADS] [CrossRef] [Google Scholar]
  30. Selsis, F., Kasting, J., Levrard, B., et al. 2007, A&A, 476, 1373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Showalter, M. R., & Hamilton, D. P. 2015, Nature, 522, 45 [NASA ADS] [CrossRef] [Google Scholar]
  32. Tholen, D. J., & Buie, M. W. 1997, Icarus, 125, 245 [NASA ADS] [CrossRef] [Google Scholar]
  33. Tholen, D. J., & Tedesco, E. F. 1994, Icarus, 108, 200 [NASA ADS] [CrossRef] [Google Scholar]
  34. Weaver, H. A., Stern, S. A., Mutchler, M. J., et al. 2006, Nature, 439, 943 [NASA ADS] [CrossRef] [Google Scholar]
  35. Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  36. Wisdom, J. 1987, AJ, 94, 1350 [NASA ADS] [CrossRef] [Google Scholar]
  37. Wisdom, J., Peale, S. J., & Mignard, F. 1984, Icarus, 58, 137 [NASA ADS] [CrossRef] [Google Scholar]
  38. Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants (Washington D.C: American Geophysical Union), 1 [Google Scholar]

All Tables

Table 1

Mean orbital and physical parameters for Pluto’s moons (Brozović et al. 2015).

All Figures

thumbnail Fig. 1

Reference angles and coordinates. The system is coplanar and is composed of an inner binary with masses m0 and m1, and an outer body m that rotates with angular velocity ω.

In the text
thumbnail Fig. 2

Stability analysis of the rotation rate close to the synchronization for δ = 0.1, α ∈ [ 0.25:0.5 ], and (BA) /C ∈ [ 0:0.9 ]. We adopt three different initial obliquity values ψ = 1°, 30°, and 60°. The color index indicates the proportion of chaotic orbits inside the studied domain: from dark blue for fully regular to red for entirely chaotic. The dashed line (Eq. (11)) roughly delimits the zone where chaotic motion is possible. The dots identify two representative cases that we detail in Fig. 3.

In the text
thumbnail Fig. 3

Cross sections representing the spin dynamics of Hydra (top) and Nix (bottom) at different initial obliquities ψ = 1°, 30°, and 60°. The dots show the main frequency of γ(t), denoted by ξ, for different initial values of \hbox{$\dot \gamma(0)$}. It is obtained by numerical integration of the three-body problem together with Eqs. (A.1) and (A.2) for the rotational motion (without dissipation). θ is initially set equal to f so the sections go through the middle of the synchronous resonance (\hbox{$\dot{\gamma}=0$}). All the remaining variables are initially set to zero. The dots are blue for stable rotation (libration or circulation), and red for chaotic motion. The distinction between stable and unstable trajectories is based on the estimate of the second derivative of ξ with respect to the initial condition \hbox{$\dot\gamma(0)$} (see Laskar 1993). The plateaus correspond to the resonance crossings (ξ/ν = k/2). The red dots randomly distributed between the plateaus indicate the overlapping of the associated spin-orbit resonances, which generates chaotic motion.

In the text
thumbnail Fig. 4

Possible evolutions for the spin of Hydra. We show the rotation \hbox{$\dot\gamma/\nu$} (in red) and the cosine of the obliquity (in blue). We numerically integrated the equations of the three-body problem together with Eqs. (A.1), (A.2), and (A.13) for the rotational motion. We adopt an initial retrograde rotation of 4.4 days with different initial values of θ, and an initial obliquity of 30°.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.