Complete Tidal Evolution of Pluto-Charon

Complete Tidal Evolution of Pluto-Charon

W. H. Cheng, Man Hoi Lee, S. J. Peale Department of Earth Sciences, The University of Hong Kong, Pokfulam Road, Hong Kong Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong Department of Physics, University of California, Santa Barbara, CA 93106
Abstract

Both Pluto and its satellite Charon have rotation rates synchronous with their orbital mean motion. This is the theoretical end point of tidal evolution where transfer of angular momentum has ceased. Here we follow Pluto’s tidal evolution from an initial state having the current total angular momentum of the system but with Charon in an eccentric orbit with semimajor axis (where is the radius of Pluto), consistent with its impact origin. Two tidal models are used, where the tidal dissipation function 1/frequency and constant, where details of the evolution are strongly model dependent. The inclusion of the gravitational harmonic coefficient of both bodies in the analysis allows smooth, self consistent evolution to the dual synchronous state, whereas its omission frustrates successful evolution in some cases. The zonal harmonic can also be included, but does not cause a significant effect on the overall evolution. The ratio of dissipation in Charon to that in Pluto controls the behavior of the orbital eccentricity, where a judicious choice leads to a nearly constant eccentricity until the final approach to dual synchronous rotation. The tidal models are complete in the sense that every nuance of tidal evolution is realized while conserving total angular momentum — including temporary capture into spin-orbit resonances as Charon’s spin decreases and damped librations about the same.

1 Introduction

Pluto has five known satellites: Charon, Nix, Hydra, Keberos, and Styx, with the latter four much smaller than Charon. Listed in Table 1 are the physical and orbital parameters of Pluto-Charon from Buie et al. (2012), unless otherwise specified. The Charon-Pluto mass ratio () is large when compared with others in the Solar System ( for Moon-Earth and for the other satellites and their planets). The barycenter of the Pluto-Charon system lies outside the surface of Pluto. Hence, some astronomers regard the pair as a binary system (Stern, 1992). The total angular momentum of the Pluto-Charon system is so large that the combined pair would be rotationally unstable (Mignard, 1981a; Lin, 1981).

The Pluto-Charon system is currently in a dual synchronous state (Buie et al., 1997, 2010), which is the endpoint of tidal evolution. As such the expected zero orbital eccentricity has been recently verified (with a 1- upper limit of ), after taking into account the effects of surface albedo variations on Pluto (Buie et al. 2012; see Table 1).

As Pluto-Charon is similar to Earth-Moon, the feasible origin of this system may be chosen from the proposed schemes for the origin of the Earth-Moon system. A giant impact of a Mars-sized body is thought to be the only viable origin of the Moon (e.g., Cameron and Ward 1976; Boss and Peale 1986; Canup 2004) to account for the large angular momentum of the system. McKinnon (1984) proposed a similar origin for Charon. If Charon accumulated from a debris disk resulting from such an impact, the initial eccentricity of Charon’s orbit would be near zero. Dobrovolskis et al. (1997, hereafter DPH97) were thereby motivated to determine the tidal evolution of Charon in a circular orbit to the current dual synchronous state in a time short compared to the age of the Solar System (see also Farinella et al. 1979) as the only possible outcome of the dissipative process. In a circular orbit, Charon would reach synchronous rotation very quickly (e.g., DPH97), and this has generally been assumed (e.g., Peale 1999). However, smoothed particle hydrodynamic (SPH) simulations by Canup (2005) showed that the results of a nearly intact capture in a glancing encounter surround the region of the system much more completely than those of disk-forming impacts. Therefore, capture where Charon comes off nearly intact after a glancing impact is favored and non-zero eccentricity would be more probable.

We are not aware of any previous attempts to examine the tidal evolution of Charon’s orbit incorporating finite eccentricity. As we shall see, Charon in an initially eccentric orbit avoids the almost immediate synchronous rotation heretofore assumed, and the varied and interesting evolutionary sequences that were suppressed in the circular orbit evolution are revealed. Depending on the ratios of rigidity and tidal dissipation function between Pluto and Charon, the eccentricity of Charon’s orbit may either grow or decay during most of the evolution (Ward and Canup, 2006). Permanent quadrupole moments of the bodies may also lead to spin-orbit resonance, and such resonances can have a significant effect on the orbital evolution.

In the following we tidally evolve the Pluto-Charon system with two tidal models distinguished by the dependence of the dissipation function on frequency : and constant. The tidal model developed in Section 2.1 has the tidal distortion of a body responding to the perturbing body a short time in the past. Constant leads to , so we call the model the constant model. In Section 2.2 we develop the equations of evolution for the constant model. Although neither of these frequency dependences represent the behavior of real solid materials (e.g., Castillo-Rogez et al. 2011) and although the evolutionary tracks are model dependent, most if not all of the possible routes from probable initial configurations to the current equilibrium state are demonstrated. In Section 2.3 we develop the contributions of rotational flattening and permanent quadrupole moment to the equations of motion. We describe the adopted system parameters and initial conditions in Section 3 and the numerical methods in Section 4. The results from both the constant and constant models with zero for Pluto and zero for both bodies are shown in Section 5.1, and the effects of non-zero and in Section 5.2, respectively. The results are discussed in Section 6, and the conclusions are summarized in Section 7.

2 Tidal Models

Tides are raised on Pluto and Charon by each other. Friction delays the response of the tidal bulge to the tide raising potential and causes tidal lag. The lagged bulge leads to angular momentum exchange between itself and the tide raising body, which leads to rotational and orbital evolution.

2.1 Constant Δt Tidal Model

The idea of approximating tidal evolution with a single bulge that lags by a constant was introduced by Gerstenkorn (1955), and developed and used by Singer (1968), Alexander (1973), Mignard (1979, 1980, 1981b), Hut (1981), and Peale (2005, 2007). The advantage of assuming a single, lagged bulge is that the tidal forces and torques can be calculated in closed form for arbitrary eccentricity and inclination. Either instantaneous or orbit-averaged tidal forces and torques can be used to determine the evolution.

The geometry is illustrated in Fig. 1, where and are the angular displacements of the axes of minimum moment of inertia from the inertial axis for Pluto and Charon, respectively, is the longitude of periapse, is the true anomaly, and and are the azimuthal spherical coordinates appearing in the potentials for Pluto and Charon, respectively. The and coordinates are those of Charon relative to Pluto with the - plane being the Pluto-Charon orbit plane. Both spin axes are assumed to be perpendicular to the orbit plane (see Section 3). The motion is thereby two dimensional, and the coordinate is ignorable.

The tidal contributions to the equations of motion for Charon for this model are found from the gradient of the tidal potential expanded to first order in (Mignard, 1980; Peale, 2007):

 MPC¨x = −3k2PGM2CR5Pr8[x+2r⋅˙rxΔtPr2+(˙ψPy+˙x)ΔtP] −3k2CGM2PR5Cr8[x+2r⋅˙rxΔtCr2+(˙ψCy+˙x)ΔtC], MPC¨y = −3k2PGM2CR5Pr8[y+2r⋅˙ryΔtPr2+(−˙ψPx+˙y)ΔtP] (1) −3k2CGM2PR5Cr8[y+2r⋅˙ryΔtCr2+(−˙ψCx+˙y)ΔtC],

where is the gravitational constant, and are the position and velocity of Charon relative to Pluto, , , , and are the mass, radius, spin angular velocity, and second order potential Love number, respectively, of body ( for Pluto and for Charon), and is the reduced mass. The first term on the right hand side of the first (second) equation in Eq. (1) is the -component (-component) of the force due to the tides raised on Pluto by Charon, and the second term is the force due to the tides raised on Charon by Pluto. The equations of motion for the spins are found from the negative of the torques on the bodies determined from the tidal forces:

 CP¨ψP = −3k2PGM2CR5PΔtPr6[˙ψP+−˙yx+y˙xr2], CC¨ψC = −3k2CGM2PR5CΔtCr6[˙ψC+−˙yx+y˙xr2], (2)

where is the moment of inertia of body about its spin axis.

Eqs. (1) and (2) can be used directly in numerical integration of the equations of motion in Cartesian coordinates. Alternatively, one can average the tidal forces and torques over an orbit to obtain the orbit-averaged equations for the variation of the spin rate , orbital semimajor axis , and eccentricity (Mignard, 1980, 1981b):

 1n⟨d˙ψidt⟩ = −3GCia6k2iΔtiM2jR5i[f1(e)˙ψin−f2(e)], (3) 1a⟨dadt⟩ = 6GMPCa8k2PΔtPM2CR5P[f2(e)(˙ψPn+AΔt˙ψCn)−f3(e)(1+AΔt)], (4) 1e⟨dedt⟩ = 27GMPCa8k2PΔtPM2CR5P[f4(e)1118(˙ψPn+AΔt˙ψCn)−f5(e)(1+AΔt)], (5)

where the subscript Charon if Pluto, and vice versa, denotes averaging over an orbit, is the mean motion, and

 f1(e) = (1+3e2+38e4)/(1−e2)9/2, f2(e) = (1+152e2+458e4+516e6)/(1−e2)6, f3(e) = (1+312e2+2558e4+18516e6+2564e8)/(1−e2)15/2, (6) f4(e) = (1+32e2+18e4)/(1−e2)5, f5(e) = (1+154e2+158e4+564e6)/(1−e2)13/2.

In Eqs. (4) and (5),

 AΔt=k2Ck2PΔtCΔtP(MPMC)2(RCRP)5 (7)

is a measure of the relative rate of tidal dissipation in Charon and Pluto, and the same as defined in Mignard (1980) and in Eq. (13) of Touma and Wisdom (1998). depends on the tidal model and we add subscripts to distinguish them whenever necessary.

The Love number measures the elastic distortion of the body in response to the second order spherical harmonic of the deforming potential. It can be modeled as (Eq. [5.6.2] of Munk and MacDonald 1960; Eq. [40a] of Peale 1973)

 k2=kf1+~μ, (8)

where is the fluid Love number and is the effective rigidity. The fluid Love number for homogeneous sphere and its reduction due to differentiation is sometimes ignored (e.g., DPH97). The effective rigidity is a dimensionless quantity, and

 ~μ=19μ2ρgR (9)

for an incompressible homogeneous sphere of radius , rigidity , density , and surface gravity (Eq. [5.6.1] of Munk and MacDonald 1960; Eq. [4.77] of Murray and Dermott 1999). For small solid body, and the approximation

 k2≈3ρgR19μ (10)

is commonly used (e.g., Peale and Cassen, 1978; Yoder and Peale, 1981; Peale, 1999, 2007). Then

 AΔt≈μPμCΔtCΔtPRCRP. (11)

For non-zero eccentricity, the orbit-averaged tidal torque vanishes at a value of , and goes to an asymptotic spin rate that increases with eccentricity. This asymptotic spin is often called a “pseudo-synchronous” state. The case for Mercury was illustrated by the curve frequency in Fig. 1 of Goldreich and Peale (1966). Pseudo-synchronous spin rate can be obtained by setting in Eq. (3):

 ˙ψpsn=f2(e)f1(e)=1+6e2+38e4+1738e6+O(e8). (12)

For the orbit-averaged equations, we ignore the periapse motion due to the tidal bulges, as it will be small initially compared to that caused by the rotational distortion of Pluto (see Section 3). Apsidal motion does not affect our discussion as far as the tidal evolution of Pluto-Charon is concerned, but it does play an important role in the study of the hypothesis that the small satellites, Nix, Hydra, Keberos, and Styx, were brought to their current orbits by mean-motion resonances with Charon (Cheng, Lee and Peale, in preparation; hereafter paper II).

2.2 Constant Q Tidal Model

A number of authors (e.g., Goldreich 1966; Goldreich and Peale 1966; Yoder and Peale 1981) have developed and applied a tidal model with a constant , based on the work of Kaula (1964). This approach has been used by the previous studies of Pluto-Charon (Farinella et al., 1979; DPH97; Ward and Canup 2006). To develop the equations of evolution for the constant model, expansions in the orbital elements are necessary, and the orbit-averaged effect on the orbit can be derived from the Gauss planetary equations. The truncation of the expansion means the equations are no longer exact and angular momentum is no longer strictly conserved.

The tidal potential acting on one of the bodies can be written as a sum of periodic terms. The frequency of the Fourier component of the tidal potential is

 σlm=ln−m˙ψ, (13)

where is the mean motion and is the spin angular velocity of the body. A phase lag is inserted into the response to each of the periodic terms in the expansion. Zahn (1977) demonstrated the procedure to expand the lagged tidal potential in for the second order spherical harmonic. The tidal evolution equations in the constant model can be derived from Eqs. (3.6)–(3.8) of Zahn (1977), using the substitution

 εlm2=k2Qsgn(σlm), (14)

where the tidal coefficient for equilibrium tide111Zahn (1977) used the notation for the apsidal motion constant, which is smaller than the Love number by a factor 2., and the phase angle is assumed to be small and independent of frequency (except for the sign). They are

 ⟨d˙ψidt⟩ = −3GM2j2Cik2iQiR5ia6[sgn(˙ψi−n)+e2Di+O(e4)], (15) 1a⟨dadt⟩ = 3nk2PQPMCMP(RPa)5[sgn(˙ψP−n)+AQsgn(˙ψC−n) (16) +e2(EP+AQEC)+O(e4)], 1e⟨dedt⟩ = nk2PQPMCMP(RPa)5[FP+AQFC+O(e2)] (17)

where

 AQ = k2Ck2PQPQC(MPMC)2(RCRP)5 (18) ≈ μPμCQPQCRCRP. (19)

Our is the same as in Eq. (5) of Yoder and Peale (1981) and agrees with the definition of in Ward and Canup (2006). The above equations conserve the total angular momentum of the system with an error of .

The coefficients , , and depend on the spin of , as listed in Table 2. Discontinuous dependence on of these coefficients arises from the sign changes of . The coefficients for and have been documented in the literature: in Ferraz-Mello et al. (2008) and and in Peale et al. (1980) and Yoder and Peale (1981). for synchronous rotation differs from the value in the literature, which included the effect of permanent quadrupole moment (see Section 12.1 of Ferraz-Mello et al. 2008 and footnote 6 of Efroimsky and Williams 2009). We treat permanent quadrupole moment separately in the next subsection.

A closer inspection of Eq. (15) reveals that the asymptotic spin rate () of the body discontinuously depends on the orbital eccentricity. Eq. (15) changes sign when increases from below to above. A body in asymptotic spin would then increase its spin from synchronous to in the spin evolution timescale. The discontinuity occurs at if we take higher order terms in into consideration (Goldreich and Peale, 1966). Fig. 1 of Goldreich and Peale (1966) showed the next discontinuity as well. We calculate its position to be at , using coefficients up to in Eq. (80) of Efroimsky and Williams (2009). As the coefficients of Eq. (16) and (17) for higher order terms in are not easily derivable, we restrict our analysis to the current order and note that our results from this model are qualitatively inaccurate for .

Before we turn to the effects of rotation induced oblateness and permanent axial asymmetry in the next subsection, we note that the equations of Zahn (1977) can also be used to derive the evolution equations expanded in eccentricity for the constant model. For small phase lag, if we let and , then Eqs. (3.6)–(3.8) of Zahn (1977) give

 1n⟨d˙ψidt⟩ = −3GCia6k2iΔtiM2jR5i[(1+152e2)˙ψin−(1+272e2)+O(e4)], (20) 1a⟨dadt⟩ = 6GMPCa8k2PΔtPM2CR5P[(1+272e2)(˙ψPn+AΔt˙ψCn) (21) 1e⟨dedt⟩ = 27GMPCa8k2PΔtPM2CR5P[1118(˙ψPn+AΔt˙ψCn)−(1+AΔt)+O(e2)]. (22)

These equations agree with the exact equations (Eqs. [3]–[5]) to , as expected. We will compare the results from these equations and from the exact equations to get an idea how good the results from the equations of the constant model are for .

2.3 Rotational Flattening and Permanent Quadrupole Moment

Rotational flattening and internal uneven mass distribution give non-zero gravitational harmonic coefficients and , where are the principal moments of inertia. The contributions of these terms to the equations of motion of a spinning rigid body orbiting another body were derived by Touma and Wisdom (1994b) and in Chapter 5 of Murray and Dermott (1999), which include the change in the spin rate of the body and the feedback on the orbit. In our aligned configuration with the rigid body rotating about its axis of maximum moment of inertia, the motion of the rigid body is confined on a plane and both equations can be greatly simplified. The spin equations become

 Ci¨ψi=−6GMjr5C22iMiR2i[(x2−y2)sin2ψi−2xycos2ψi], (23)

where the notation is as shown in Fig. 1. The contributions of and to the acceleration of Charon relative to Pluto are

 MPC¨x = GMPMC{−3J2PR2Px2r5 + 3(C22PR2Pcos2ψP+C22CR2Ccos2ψC)[2r5−5(x2−y2)r7]x + 3(C22PR2Psin2ψP+C22CR2Csin2ψC)[2r5−10x2r7]y} MPC¨y = GMPMC{−3J2PR2Py2r5 (24) + 3(C22PR2Pcos2ψP+C22CR2Ccos2ψC)[−2r5−5(x2−y2)r7]y + 3(C22PR2Psin2ψP+C22CR2Csin2ψC)[2r5−10y2r7]x},

where of Charon is omitted.

3 Parameters and Initial Conditions

In all our calculations, Charon is assumed to start in an eccentric orbit with semimajor axis (Pluto radii), consistent with its origin in a nearly intact capture from a glancing impact on Pluto (Canup, 2005). Since a significant portion of the angular momentum of the impactor is transferred to the spin of the target in the collision (see Table 1 of Canup 2005), the spin axis of Pluto should be close to being perpendicular to Charon’s orbit initially. The spin axis of Charon, which could be inclined from the orbit normal initially, would quickly approach a Cassini state with the spin axis close to the orbit normal, on a timescale comparable to the timescale for Charon to reach asymptotic spin rate in at least one of the tidal models (see Eq. [53] of Hut 1981). Moreover, tidal evolution of the orbit and spin rates is unaffected to first order in orbital inclination (Hut, 1981; Ferraz-Mello et al., 2008; Efroimsky and Williams, 2009). Thus, to reduce the complexity and the parameter space of the problem, we only examine the aligned configuration of Pluto-Charon, where the orbit normal aligns with their spin axes.

We specify , , and as our initial conditions, and initial is calculated by assuming the same total angular momentum as the current Pluto-Charon system:

 (25)

where and are the current separation and mean motion of Pluto-Charon, respectively. Charon’s spin angular momentum is always small compared to the total. The numerical value is determined under the assumption that the dimensionless moments of inertia and . The numerical value for follows from a two-layer model (DPH97) with a rocky core (density ) and an icy mantle (density ). Pluto should be so differentiated by the giant impact if it was not earlier (McKinnon, 1989; Canup, 2005). The internal structure of Charon is uncertain (e.g., McKinnon et al., 2008, Section 6.2) and we assume that Charon remains homogeneous, i.e., (McKinnon, 1989; DPH97).

Pluto’s Love number , computed using Eq. (8) with and of water ice. For the constant model, we adopt seconds, same as that for the Earth (Mignard, 1980; Touma and Wisdom, 1994a). For the constant model, we adopt , as typically assumed for solid bodies (DPH97; Tables 4.1 and 4.2 of Murray and Dermott 1999). These parameters are fixed throughout the tidal evolution. Our incomplete knowledge of the physics of tides and of the composition and internal structure of Pluto means that the actual values of these parameters are not well constrained. Yet uncertainties in these parameters only affect the overall timescale of tidal evolution, as far as spin-orbit resonance is not included. If the rigidities () and dissipation ( or ) of Pluto and Charon are comparable, from Eqs. (11) and (19), one would expect and . Alternatively, if the Love numbers () and dissipation ( or ) of Pluto and Charon are comparable, one would expect and . In our integrations, we focus on those values of and that can keep roughly constant until the end of the tidal evolution. If is too large, the orbit circularizes quickly, and the tidal evolution would be similar to that already studied by DPH97. If is too small, can approach , and the system can become unstable (see Section 5). The evolutions with roughly constant throughout most of the tidal evolution are also the most likely ones that allow migration of the small satellites in resonances, since resonances cannot be maintained if is too small and become unstable if is too large (Ward and Canup, 2006; Lithwick and Wu, 2008). We discuss the details in paper II.

We estimate the largest value that of Pluto is likely to be by the hydrostatic value just after the impact that captured Charon. For rotation about the axis of maximum moment of inertia, the changes in the principal components of the inertia tensor from rotation are given by (e.g., Peale 1973)

 ΔA=ΔB = −kfPR5P˙ψ2P9G, ΔC = +2kfPR5P˙ψ2P9G, (26)

where is the fluid Love number of Pluto. Then

 J2P=ΔC−(ΔA+ΔB)/2MPR2P=kfPR3P˙ψ2P3GMP. (27)

For an initial Pluto spin period of hours, compatible with our typical initial conditions of , , and , if (by analogy with the Earth) to (for homogeneous sphere). The large value means that inclusion of higher order terms in the rotational distortion would be appropriate, but we have not tried to obtain a more accurate estimate of . The large reflects the fact that the estimated spin of Pluto for is close to the limit of rotational instability according to the ratio of rotational to gravitational binding energy for a homogeneous Pluto (DPH97 and references therein).

Satellite motion around an oblate body deviates from that of Keplerian. If is negligibly small ( or less), corrections on the mean motion and the rate of periapse precession are given by Eqs. (6.244) and (6.249) of Murray and Dermott (1999):

 n2 = G(MP+MC)a3[1+32J2P(RPa)2+O(RPa)4], (28) ˙ϖ = [G(MP+MC)a3]1/2[32J2P(RPa)2+O(RPa)4]. (29)

We estimate the correction on the mean motion to be less than initially. It decreases further with and hence can be ignored. On the other hand, the precession due to the oblateness of Pluto is much larger than that from the oblateness of Charon and tidal deformation during most of the evolution, and the remnant supported by internal stress is likely to be greater than the hydrostatic value and the tidal value for the current configuration.

The value of initial we use in our integrations is or , representing the two extreme cases with no or very fast precession of the orbit. As Charon moves outward, we assume that decreases with (Eq. [27]). We ignore the smaller effect of Charon’s , and we choose or for the integrations, where the latter value is comparable to the measured values of other nearly spherical solid bodies in the Solar System.

4 Numerical Methods

4.1 Runge-Kutta Codes

For the calculations without and , the most efficient way to evolve Pluto-Charon is to solve the orbit-averaged tidal evolution equations: Eqs. (3)–(5) and Eqs. (15)–(17) in the two tidal models. We use the 4th order Runge-Kutta method with an adaptive time step (e.g., Press et al. 1992). We find that it is simpler and more accurate to treat as variables rather than in these codes, through

Discontinuities of coefficients in the equations of the constant model need special treatment to prevent the adaptive time step algorithm from crashing when the system comes across them. One can smooth the discontinuities by assuming that has a very weak power dependence on frequency using Eq. (86) of Efroimsky and Williams (2009), but we find it convenient to use a simple smoothing function recommended by Rauch and Holman (1999). We modify the smoothing function in their Eq. (29) to give a step from (at ) to (at ) with controllable steepness:

 κ(δ,ϵ)=tanh[4δ/ϵ1−(δ/ϵ)2]. (31)

Here denotes the percentage difference of from discontinuity, and is an adjustable parameter within which smoothing is applied. The advantage of this smoothing function is that all orders of derivatives vanish at both the beginning and end of the transition over the discontinuity. The position and amplitude of the smoothing are transformed to replace the discontinuous step, such that when , the coefficients would be calculated using the above equation instead. The parameter can be set to at most 20% without overlapping for the discontinuities at and 3/2. We use % in most cases. In the -body codes described in the next subsection, smoothing of the discontinuities in the constant model can be turned off. We compare the results from -body calculations with and without smoothing and find very small differences for our typical of .

We perform three classes of tests on the Runge-Kutta codes. We test the implementation of each equation separately in the first class of test. By setting the right hand side of all but one of the tidal evolution equations to zero, analytical solutions are available for each equation, except for Eq. (5). In the constant model, analytical solution of each equation is valid only for no discontinuity crossing. We monitor the angular momentum budget of the system as a second class of test. In the constant model, total angular momentum of the system is conserved to better than the tolerance parameter ( for the results presented in Section 5) of the adaptive time step algorithm throughout the evolution to the current dual synchronous state. In the constant model, Eqs. (15)–(17) do not conserve the total angular momentum of the system. The discrepancy becomes worse when eccentricity is large. Smoothing the discontinuities also contributes to the error. When tolerance is small (including the adopted ), the angular momentum error of the system is dominated by the order of the equations in . The third class of test aims at testing the implementation of Eq. (5). In the constant model, Hut (1981) found that the tidal evolution between a body and a point mass (i.e., tides are raised on one body only) depends on the value of only for given initial and , where is the ratio of orbit to spin angular momentum at the dual synchronous state. Figs. 5–8 of Hut (1981) show the flow lines of tidal evolution in the () space for four different values of , where is the separation in the dual synchronous state. By treating either Pluto or Charon as a point mass, we are able to reproduce the flow lines in all four figures by choosing suitable initial conditions, including those with initial .

4.2 N-body Codes

As the effects of are on suborbital timescale, no analytic, orbit-averaged equations are available for evolving spins. Thus, to study the consequences of non-zero quadrupole moments represented by and on the tidal evolution, it is necessary to perform -body integrations.

For the constant model, the equations of motion in Cartesian coordinates with the instantaneous tidal forces and torques and the effects of and (Eqs. [1], [2], [23], and [24]) can be integrated directly. We have written a code implementing these equations using the Bulirsch-Stoer method, and its accuracy is verified by the conservation of angular momentum to more than 8 significant figures for integration from typical initial conditions to the dual synchronous state. This Bulirsch-Stoer code has the advantage of solving the exact equations of motion. However, it is slow for realistic values of , because the tidal forces and torques are computed many times over an orbit, even though they are weak and affect the evolution only on the tidal evolution timescale. In addition, this approach does not work for the constant model, where the instantaneous tidal forces and torques are not available.

Thus, for both tidal models, we also modify the Wisdom-Holman (1991) integrator in the SWIFT222hal/swift.html. package (Levison and Duncan, 1994) to simulate the tidal, rotational and axial-asymmetry effects. The SWIFT package allows non-zero for the central body (i.e., Pluto in our case). For non-zero initial , we adjust it to decrease throughout the evolution according to Eq. (27). Other effects are imposed following the approach of Lee and Peale (2002):

 Ee(mτ2)Ea(mτ2)E˙ψ(mτ2)EC22(τ2)Erot(τ)EWH(τ)EC22(τ2)m copiesE˙ψ(mτ2)Ea(mτ2)Ee(mτ2). (32)

Here denotes a complete step of time step in the Wisdom-Holman scheme, and the other ’s are evaluations for each effect. In each evaluation, only those variables concerned are evolved and others are kept constant.

We follow the method presented by Touma and Wisdom (1994b) for rotation and the effects of . We represent the pointing directions of the long axes of both bodies by unit vectors. rotates them according to the instantaneous of the bodies. The bodies are treated as axisymmetric in , and hence commutes with . changes the spins and velocities of the bodies according to Eqs. (23) and (24). These two evaluations, and , are the spin analog to the leapfrog integration, except that the feedback on the orbit has to be included. Our -body simulations typically start with the long axis of both Pluto and Charon pointing along the inertial -axis and Charon at periapse on the -axis.

The substeps , , and correspond to the changes in , , and in the orbit-averaged tidal evolution equations. Their sequence is chosen such that the computationally expensive in Eq. (6) are calculated once only in each half-step. Eq. (30) is not used here as and are applied sequentially. In the constant model, solving Eq. (4) analytically either involves complicated expressions or transformations back and forth between and in each step, and we use explicit midpoint method in and and analytical solution in . In the constant model, analytical expressions for all the substeps are available, assuming all coefficients are constant during the step.

The parameter is an integer, and tides should be applied on tidal evolution timescale by using a large . This is done for numerical efficiency and to reduce roundoff error. There is an error introduced by the conversion between the positions and velocities and the osculating orbital elements. In the regime of eccentricity and step size of our problem, this error is tested to be secularly increasing with the number of tidal steps taken, which can be significantly reduced by the use of a large .

We use an initial seconds, which is about 60 steps per orbit for initial , and an initial seconds. Since is kept constant and tidal evolution is proportional to a large negative power of , the relative angular momentum error introduced by our second order solution saturates at soon after starts to increase. We integrate the system up to a point when has increased by a significant factor (), then we increase and . The step size is increased by a factor of 5 to give a similar number of steps per orbit as initially. The tidal step size is increased by a factor of 100, which is smaller than one would use to keep the right hand side of the tidal equations comparable in magnitude as initially. We choose this factor of 100 so that the increase in the step size does not further increase the already saturated relative angular momentum error of the system. Because of the slower evolution rate in the constant model, we increase the step size once more, when and are around their maximum for the smallest (see below). This time is increased by a factor of 2, and the tidal step size is increased by a factor of 10.

We perform several tests on the -body codes. When and are set to zero, results from the Wisdom-Holman and Bulirsch-Stoer codes coincide with those from the Runge-Kutta codes. For uniform rotation, the pointing directions of unit vectors are tested to change at the expected rate. Precession due to an oblate Pluto is tested to agree with the expected rate given by Eq. (29). The effects without tides are tested to conserve total angular momentum of the system. Finally, we compare results from the Wisdom-Holman and Bulirsch-Stoer codes, and they agree in all cases examined.

5 Results

5.1 Tidal Evolution with Zero J2p and C22

In this subsection, we present the tidal evolution of Pluto-Charon for both the constant and constant models, with and for both Pluto and Charon. The results are obtained using the Runge-Kutta codes, unless otherwise specified.

Figs. 2 and 3 show the evolution in both tidal models using our typical initial conditions of and for a range of initial eccentricities. The relative rate of tidal dissipation in Charon and Pluto, and defined in Eqs. (7) and (18), is chosen such that is kept roughly constant throughout most of the evolution (, and and for initial and , respectively). Note that all evolutions shown reach the current dual synchronous state of Pluto-Charon, as predicted by DPH97. For , the spin of Charon drops to synchronous quickly, as estimated by DPH97. However, the assumption that the spin of Charon is synchronous throughout most of the evolution does not necessarily hold for non-zero . For constant , the spin of Charon achieves the pseudo-synchronous state quickly instead, and evolves according to afterwards (Eq. [12]). For constant , the asymptotic spin rate for is no longer synchronous but , as mentioned in Section 2.2. Hence, for larger initial and spin of Charon above , the spin of Charon first reaches and stays at , and falls to synchronous depending on the eccentricity evolution (see, e.g., the evolution with initial in Fig. 3). The angular momentum carried in the rotation of Charon is small throughout the evolution, as the moment of inertia of Charon is much smaller than that of Pluto (a factor of ) and its spin stays within a factor of two of synchronous for mild eccentricity (). Note that rises to (higher for larger eccentricity) before falling to synchronous, even though the rotation rate of Pluto is monotonically decreasing. This initial rise in is due to decreasing faster than .

The tidal evolution can be drastically affected by the relative rate of tidal dissipation in Charon and Pluto ( or ). If the spin angular velocity of Pluto sufficiently exceeds the orbital angular velocity of Charon at periapse, the maximum tide at that point in the orbit gives Charon a kick that tends to increase the eccentricity. Otherwise tides raised on Pluto will decrease the eccentricity (see Eqs. [5] and [17]). Since Charon’s rotation stays within a factor of two of synchronous for mild eccentricity, tides raised on Charon typically damp the eccentricity. Since Pluto will be initially spinning very fast, we expect there will be a tendency for tides raised on Pluto to increase the orbital eccentricity that will be counteracted by tides raised on Charon tending to decrease the eccentricity. Which wins depends on the value of .

Fig. 4 shows the evolution in the constant model with initial and a range of . While can be kept more or less constant as increases if is used, larger (smaller) would result in decreasing (increasing) throughout most of the evolution. If the eccentricity is still large () when reaches the current value (), then by the conservation of angular momentum, it is expected that would overshoot before coming back to the current value when decays. This is seen clearly in the case with in Fig. 4. For initial larger than those used in Figs. 2 and 4 (e.g., ), and can drop initially as declines rapidly if .

Fig. 5 shows the evolution in the constant model with initial , and a range of . We choose this initial spin of Charon so that the discontinuity is avoided. As long as no discontinuity is crossed, the coefficients in the tidal evolution equations (Eqs. [15]–[17]) remain constant and the eccentricity evolution is exponential. Hence the lines in the versus graph look straight before the spin of Pluto drops to below . The spin of Charon stays slightly larger than synchronous, as described by Greenberg and Weidenschilling (1984), with the difference depending on both and the smoothing range (note from Eq. [15] that, without smoothing, for and for and ). Same as in the constant model, a suitable value of can be chosen to keep nearly constant.

In Section 3 we show that and if Pluto and Charon have similar tidal response (in terms of rigidities or Love numbers) and dissipation (in terms of or ). For constant , this range includes the value () required to keep nearly constant. For constant , this range is significantly below the value () required to keep nearly constant. Fig. 6 shows the evolution for . We include integrations for and using the Bulirsch-Stoer code. For the case of axial symmetry for both bodies (), the dashed curves in Fig. 6 show that approaches 1, while the growth in well beyond the current value drives to values far above the 2:1 spin-orbit resonance and to . The system can become unstable if the apoapse distance approaches Pluto’s Hill sphere radius. The existence of the additional satellites Nix, Hydra, Keberos, and Styx preclude even the large values of eccentricity on the way to stable equilibrium in the case of axial asymmetry, if these satellites were in orbit prior to the tidal expansion of Charon’s orbit. For this reason, we have considered larger values of that keep the value of at reasonably low values. The case with permanent quadrupole moments () in Fig. 6 is discussed in the next subsection.

5.2 Tidal Evolution with Non-zero J2p or C22

In this subsection, we examine the effects of and on the tidal evolution of Pluto-Charon using results from -body integrations.

Fig. 7 shows the comparison between initial and evolution in the constant model, with initial and . For the case with initial , we assume that decreases with as in equation (27). Only the early stages of the evolution are affected by when Charon is close, and we see in Fig. 7 that, with the exception of the fluctuations in the osculating eccentricity, the overall evolution is not that different from the case with . In Section 3 we estimate from Eq. (28) that the change of mean motion by is initially. Oscillations of in Fig. 7 has initial amplitude , which is about the order of the change of mean motion due to . Because the effect of is relatively minor compared to that of , we shall usually not include it in the evolutionary calculations.

When are non-zero, Charon can be captured into spin-orbit resonances in both tidal models, which may result in very different evolution when compared with the zero cases. The profound effect of non-zero on the evolution is demonstrated in Fig. 6, where we compare integrations with