Post-new-horizons Orbits and Masses for the Satellites of Pluto

Charon's astrometry is listed in Table 1. Harrington & Christy (1980) reported on position measurements of Charon with respect to Pluto that predate its discovery in 1978 (Christy & Harrington 1978). Pluto and Charon were not fully resolved in the photographic plates, but it was possible to measure the image elongation and orientation and express Charon's position relative to Pluto in terms of position-angle Δθ and separation Δρ. These early observations had measurement uncertainties of ∼150 miliarcseconds (mas). Soon after its discovery, Charon was also observed via speckle interferometry from several large telescopes around the world (Bonneau & Foy 1980; Baier et al. 1982; Hege et al. 1982; Hetterich & Weigelt 1983; Hege & Drummond 1984; Baier & Weigelt 1987; Beletic et al. 1989). In 1984 and 1985, Beletic et al. (1989) reported 56 speckle interferometric measurements of Charon's separation from Pluto in terms of R.A. projected onto a tangential plane (${\rm{\Delta }}\alpha \cos (\delta )$) and decl. (Δδ).

Table 1. Charon Astrometric Observations

Time Span	Type	Points	rms	Type	Points	rms	References
			(mas)			(mas)
1965 Apr–1979 Apr	ρΔθ	23	105	Δρ	12	146	Harrington & Christy (1980)
1980 Jun–1984 Feb	ρΔθ	2	18	Δρ	2	83	Hege et al. (1982), Hege & Drummond (1984)
1980 Jun	ρΔθ	5	15	Δρ	5	29	Bonneau & Foy (1980)
1981 Apr	ρΔθ	2	11	Δρ	2	22	Hetterich & Weigelt (1983)
1981 Apr–1983 Jul	ρΔθ	8	41	Δρ	8	82	Baier et al. (1982), Baier & Weigelt (1987)
1984 May–1985 May	${\rm{\Delta }}\alpha \cos (\delta )$	56	87	Δδ	56	73	Beletic et al. (1989)
1991 Aug–1993 Aug	Sample	28	11	Line	28	15	Null et al. (1993), Null & Owen (1996)
1992 May–1993 Aug	${\rm{\Delta }}\alpha \cos (\delta )$	60	3	Δδ	60	3	Tholen & Buie (1997), Buie et al. (2012)
1992 Feb–1992 Mar	${\rm{\Delta }}\alpha \cos (\delta )$	80	30	Δδ	80	12	Young et al. (1994)
1998 Mar	α	5	14	δ	5	15	Olkin et al. (2003)
2002 Jun–2003 Jun	${\rm{\Delta }}\alpha \cos (\delta )$	384	3	Δδ	384	2	Buie et al. (2006), Buie et al. (2012)
2005 May	ρΔθ	2	24	Δρ	2	8	Weaver et al. (2006)
2005 Jun	${\rm{\Delta }}\alpha \cos (\delta )$	45	<1	Δδ	45	<1	B. Sicardy et al. (2006, personal communication)
2005 Jul	α	1	4	δ	1	4	Herald et al. (2020)
2006 Apr	${\rm{\Delta }}\alpha \cos (\delta )$	1	11	Δδ	<1	6	Sicardy et al. (2006)
2006 Feb	${\rm{\Delta }}\alpha \cos (\delta )$	10	3	Δδ	10	3	Brozović et al. (2015)
2007 Mar–2007 Jun	${\rm{\Delta }}\alpha \cos (\delta )$	215	3	Δδ	215	3	Buie et al. (2012)
2008 Jun	${\rm{\Delta }}\alpha \cos (\delta )$	1	1	Δδ	1	<1	Sicardy et al. (2011)
2010 Apr–2010 Sep	${\rm{\Delta }}\alpha \cos (\delta )$	95	3	Δδ	95	2	Brozović et al. (2015)
2010 Apr–2010 Sep	${\rm{\Delta }}\alpha \cos (\delta )$	240	3	Δδ	240	2	Buie et al. (2012)
2011 Jun	${\rm{\Delta }}\alpha \cos (\delta )$	2	<1	Δδ	2	<1	B. Sicardy et al. (2015, personal communication),
2011 Jun–2011 Sep	${\rm{\Delta }}\alpha \cos (\delta )$	104	3	Δδ	104	2	Brozović et al. (2015)
2012 Jun–2012 Jul	${\rm{\Delta }}\alpha \cos (\delta )$	548	2	Δδ	548	2	Brozović et al. (2015)
2012 Jul	ρΔθ	1	4	Δρ	1	1	Howell et al. (2012)
2014 Jul	Sample	74	81	Line	74	105	Pelletier et al. (2016), Thompson et al. (2016)
2015 Jan–Feb	Sample	75	39	Line	75	53	Pelletier et al. (2016), Thompson et al. (2016)
2015 Apr–Jul	Sample	379	126	Line	379	158	Pelletier et al. (2016), Thompson et al. (2016)
2014 Sep–2019 Nov	AL	232	1.4	AC	232	142	Tanga et al. (2023)

Note. Summary of Charon's optical astrometry. We list the time of observations, type of measurement, number of astrometry points, and the rms of the residuals. Note that Pelletier et al. (2016) and Thompson et al. (2016) are the New Horizons optical navigation data obtained as the spacecraft was approaching the system.

Download table as:  ASCIITypeset image

Charon's orbit was constrained with lightcurves during the Pluto–Charon mutual event season that lasted from late 1984 until late 1990 (Tholen et al. 1987). There were 64 lightcurves reported: 15 from Palomar (Buratti et al. 1995), 39 from Maunakea Observatory (Tholen 2008), and 10 from McDonald Observatory (Young 1993; Young & Binzel 1993). These lightcurves are available on NASA Planetary Data System Small Bodies Node. ¹

Null et al. (1993) and Null & Owen (1996) used the HST Wide Field Camera to measure absolute positions of Pluto and Charon as "Samples" and "Lines" in camera coordinates. The goal was to detect Pluto's barycentric motion and to estimate the individual masses and bulk densities. We used Gaia DR3 reference stars (Tanga et al. 2023) to re-point these images. Buie et al. (2012) reported more relative measurements (${\rm{\Delta }}\alpha \cos (\delta )$ and Δδ) of Charon. These HST measurements were originally obtained by Tholen & Buie (1997), but Buie et al. (2012) rereduced and corrected them for the center-of-light and center-of-body offset based on an early model for Pluto's albedo.

Young et al. (1994) used Maunakea's Observatory 2.2 m telescope on six nights in 1992 February to demonstrate that it is "... possible to perform precise astrometry of Pluto and Charon from large, ground-based telescope using PSF fitting techniques." Their observations were successful, but the ground-based astrometry remains sparse to this day. During the several decades of followup, Charon's astrometry has been predominantly obtained with the HST. Most of these measurements were relative observations (Buie et al. 2006, 2012; Weaver et al. 2006; Brozović et al. 2015), but Olkin et al. (2003) reported five absolute measurements (α, δ) that were important for early Pluto/Charon mass-ratio determination (0.122 ± 0.008).

B. Sicardy et al. (2006, personal communication) reported relative astrometry of Charon derived from adaptive optics observations (NACO/VLT at Paranal) from 2005. One more point was obtained in 2006 as Hydra was followed up from the same site (Sicardy et al. 2006). The last ground-based measurement with speckle interferometry (ρΔθ, Δρ) was obtained by Howell et al. (2012).

Stellar occultations yield high-accuracy measurements of Charon's position. Sicardy et al. (2011) reported a Pluto and Charon occultation of the same star on 2008 June 22. Two more Charon occultation points (2011 June 4 and 23) were reported by B. Sicardy via personal communication. We also used the 2005 July 11 stellar occultation reported by Herald et al. (2020).

The New Horizons Optical Navigation Team reported Charon's astrometry (Pelletier et al. 2016; Thompson et al. 2016) based on raw images that are publicly available on NASA Planetary Data System (see footnote 1) archive. There were more than five hundred measurements of the absolute position of Charon from 2014 July 19 to 2015 July 15. The close approach occurred on 2015 July 14.

The Gaia space telescope detected Charon in the instances when its location was at least 0.25 mas away from Pluto in the scan direction (AL direction; Tanga et al. 2023). We identified 31 transits with Charon observations in the Gaia FPR. ² The data from different transits are at least 90 minutes apart and the transits span 14 separate dates from 2014 to 2019; 23 out of 31 transits occurred in 2019. Gaia astrometry is reported in bursts of 6–10 α and δ measurements that are 5 s apart. The α and δ uncertainties have values of hundreds of mas and the measurements have high correlations (>0.99999). The correlation of α and δ measurements is a result of coordinate transformation from the original (uncorrelated) AL (along the scan) and AC (across the scan) measurements (Gaia Collaboration et al. 2018). A recent press release by the Gaia project ³ stated that this is only a preliminary data reduction, and that Pluto and Charon data will get fully reduced in Gaia DR5 release.

The measurements of small satellites' positions are listed in Tables 2–5. Most of them were obtained with the HST. There were a few more measurements obtained with NACO/VLT and Magelan I Baade telescopes (Fuentes & Holman 2006; Sicardy et al. 2006). Porter & Canup (2023) reprocessed most of the HST data with respect to the stars in the Gaia DR2 catalog, and these measurements are of especially high quality. Furthermore, Porter & Canup (2023) also reduced the New Horizons LORRI 4 × 4 optical images, and reported hundreds of measurements obtained during two months prior to the closest approach. The measurements in Porter & Canup (2023) were reported relative to the JPL DE430 planetary ephemeris location of the Pluto system barycenter (Folkner et al. 2014). Utilizing the HST or New Horizons astrometry requires knowledge of spacecraft position at each observation time. The reconstructed trajectory files for the spacecraft are available from NASA's Navigation and Ancillary Information Facility website (Acton 1996).

Table 2. Nix Astrometric Observations

Time Span	Type	Points	rms	Type	Points	rms	References
			(mas)			(mas)
2002–2003	${\rm{\Delta }}\alpha \cos (\delta )$	11	17	Δδ	11	12	Buie et al. (2006)
2005	${\rm{\Delta }}\alpha \cos (\delta )$	2	2	Δδ	2	2	Porter & Canup (2023)
2006	${\rm{\Delta }}\alpha \cos (\delta )$	3	<1	Δδ	3	2	Porter & Canup (2023)
2006	${\rm{\Delta }}\alpha \cos (\delta )$	2	136	Δδ	2	121	Fuentes & Holman (2006)
2007	${\rm{\Delta }}\alpha \cos (\delta )$	16	2	Δδ	16	2	Porter & Canup (2023)
2010	${\rm{\Delta }}\alpha \cos (\delta )$	10	1	Δδ	10	2	Porter & Canup (2023)
2011	${\rm{\Delta }}\alpha \cos (\delta )$	6	1	Δδ	6	<1	Porter & Canup (2023)
2012	${\rm{\Delta }}\alpha \cos (\delta )$	6	<1	Δδ	6	<1	Porter & Canup (2023)
2013	${\rm{\Delta }}\alpha \cos (\delta )$	12	<1	Δδ	12	<1	Porter & Canup (2023)
2014	${\rm{\Delta }}\alpha \cos (\delta )$	11	<1	Δδ	11	<1	Porter & Canup (2023)
2015	${\rm{\Delta }}\alpha \cos (\delta )$	40	<1	Δδ	40	<1	Porter & Canup (2023)
2018	${\rm{\Delta }}\alpha \cos (\delta )$	13	<1	Δδ	13	<1	Porter & Canup (2023)
2019	${\rm{\Delta }}\alpha \cos (\delta )$	9	<1	Δδ	9	<1	Porter & Canup (2023)

Note. Summary of Nix's optical astrometry and fit residuals. Note that Fuentes & Holman (2006) from 2006 are ground-based observations. The rest is the HST astrometry (Porter & Canup 2023).

Download table as:  ASCIITypeset image

Table 3. Hydra Astrometric Observations

Time Span	Type	Points	rms	Type	Points	rms	References
			(mas)			(mas)
2002–2003	${\rm{\Delta }}\alpha \cos (\delta )$	11	10	Δδ	11	12	Buie et al. (2006)
2005	${\rm{\Delta }}\alpha \cos (\delta )$	2	<1	Δδ	2	<1	Porter & Canup (2023)
2006	${\rm{\Delta }}\alpha \cos (\delta )$	1	32	Δδ	1	11	Sicardy (2006, personal communication)
2006	${\rm{\Delta }}\alpha \cos (\delta )$	2	135	Δδ	2	88	Fuentes & Holman (2006)
2006	${\rm{\Delta }}\alpha \cos (\delta )$	3	<1	Δδ	3	1	Porter & Canup (2023)
2007	${\rm{\Delta }}\alpha \cos (\delta )$	16	<1	Δδ	16	<1	Porter & Canup (2023)
2010	${\rm{\Delta }}\alpha \cos (\delta )$	10	<1	Δδ	10	<1	Porter & Canup (2023)
2011	${\rm{\Delta }}\alpha \cos (\delta )$	6	<1	Δδ	6	<1	Porter & Canup (2023)
2012	${\rm{\Delta }}\alpha \cos (\delta )$	6	<1	Δδ	6	<1	Porter & Canup (2023)
2013	${\rm{\Delta }}\alpha \cos (\delta )$	12	<1	Δδ	12	<1	Porter & Canup (2023)
2014	${\rm{\Delta }}\alpha \cos (\delta )$	11	<1	Δδ	11	<1	Porter & Canup (2023)
2015	${\rm{\Delta }}\alpha \cos (\delta )$	40	<1	Δδ	40	<1	Porter & Canup (2023)
2018	${\rm{\Delta }}\alpha \cos (\delta )$	13	<1	Δδ	13	<1	Porter & Canup (2023)
2019	${\rm{\Delta }}\alpha \cos (\delta )$	9	<1	Δδ	13	<1	Porter & Canup (2023)

Note. Summary of Hydra's optical astrometry and fit residuals. Note that Fuentes & Holman (2006) and B. Sicardy et al. (2006, personal communication) from 2006 are ground-based observations. The rest is the HST astrometry (Porter & Canup 2023).

Download table as:  ASCIITypeset image

Table 4. Kerberos Astrometric Observations

Time Span	Type	Points	rms	Type	Points	rms	References
			(mas)			(mas)
2005	${\rm{\Delta }}\alpha \cos (\delta )$	2	<1	Δδ	2	1	Porter & Canup (2023)
2006	${\rm{\Delta }}\alpha \cos (\delta )$	3	<1	Δδ	3	<1	Porter & Canup (2023)
2010	${\rm{\Delta }}\alpha \cos (\delta )$	10	2	Δδ	10	<1	Porter & Canup (2023)
2011	${\rm{\Delta }}\alpha \cos (\delta )$	6	1	Δδ	6	1	Porter & Canup (2023)
2012	${\rm{\Delta }}\alpha \cos (\delta )$	6	<1	Δδ	6	<1	Porter & Canup (2023)
2013	${\rm{\Delta }}\alpha \cos (\delta )$	10	<1	Δδ	10	<1	Porter & Canup (2023)
2014	${\rm{\Delta }}\alpha \cos (\delta )$	8	<1	Δδ	8	<1	Porter & Canup (2023)
2015	${\rm{\Delta }}\alpha \cos (\delta )$	32	1	Δδ	32	1	Porter & Canup (2023)
2018	${\rm{\Delta }}\alpha \cos (\delta )$	13	2	Δδ	13	1	Porter & Canup (2023)
2019	${\rm{\Delta }}\alpha \cos (\delta )$	8	2	Δδ	8	1	Porter & Canup (2023)

Note. Summary of the astrometry for Kerberos and fit residuals.

Download table as:  ASCIITypeset image

Table 5. Styx Astrometric Observations

Time Span	Type	Points	rms	Type	Points	rms	References
			(mas)			(mas)
2010	${\rm{\Delta }}\alpha \cos (\delta )$	6	3	Δδ	6	5	Porter & Canup (2023)
2011	${\rm{\Delta }}\alpha \cos (\delta )$	5	2	Δδ	5	5	Porter & Canup (2023)
2012	${\rm{\Delta }}\alpha \cos (\delta )$	4	1	Δδ	4	2	Porter & Canup (2023)
2013	${\rm{\Delta }}\alpha \cos (\delta )$	5	2	Δδ	5	4	Porter & Canup (2023)
2014	${\rm{\Delta }}\alpha \cos (\delta )$	5	2	Δδ	5	2	Porter & Canup (2023)
2015	${\rm{\Delta }}\alpha \cos (\delta )$	16	2	Δδ	16	2	Porter & Canup (2023)
2018	${\rm{\Delta }}\alpha \cos (\delta )$	6	3	Δδ	6	2	Porter & Canup (2023)
2019	${\rm{\Delta }}\alpha \cos (\delta )$	8	3	Δδ	8	2	Porter & Canup (2023)

Note. Summary of the astrometry for Styx's and fit residuals.

Download table as:  ASCIITypeset image

All astrometric observations were processed with models for the observables from Murray (1983), Nautical Almanac Office (1987), Riedel et al. (1990), and Seidelmann (1992). Relative measurements of a satellite's position with respect to Pluto are affected by the planet's surface albedo distribution. Observers most often report center of light as opposed to center-of-body astrometry. This predominantly refers to Charon's astrometry, but it also affects a few measurements of Nix and Hydra reported in 2006 (Fuentes & Holman 2006). Some early HST data on Charon (Null et al. 1993; Null & Owen 1996) were already corrected with an albedo model for the surface of Pluto (Buie et al. 1992), so we kept these measurements in their original format. Otherwise, we used the Buratti et al. (2017) albedo map to correct the astrometry (Young & Binzel 1993; Olkin et al. 2003; Fuentes & Holman 2006; Sicardy et al. 2006; Weaver et al. 2006; Howell et al. 2012; Brozović et al. 2015). This albedo map was also used in Porter & Canup (2023) to remove the scattered light from Pluto in their data reduction.

Our orbit model is based on equations of motion and variational equations described in Peters (1981). This approach has been used repeatedly for the construction of numerically integrated ephemerides for satellites of outer planets, e.g., Jacobson & Owen (2004), Jacobson et al. (2012), and Brozović et al. (2015, 2020). The equations are formulated in the Cartesian coordinates with the Pluto system barycenter at the origin and referenced to the International Celestial Reference Frame (ICRF). We used Barycentric Dynamical Time (TDB) in order to remain consistent with the JPL ephemerides convention.

The model includes external perturbations from the Sun, Jupiter, Saturn, Uranus, and Neptune. JPL planetary ephemeris DE440 (Park et al. 2021) was used to provide positions and the masses of the Sun and planets, see Table 6 for the values and sources of the masses. The integration uses variable order, variable step size Gauss–Jackson method. The integration step is controlled by the maximum acceptable error on the components of satellites' velocity vectors (${v}_{\max }={10}^{-10}$ km s⁻¹). The largest time step was set to 1800 s.

We used a weighted least-squares algorithm in order to estimate our model's dynamical parameters (Lawson & Hanson 1974; Tapley et al. 2004). We estimated the starting epoch Cartesian state vector and GM for each satellite, system GM, and additional observational model parameters. We also applied camera pointing biases to the New Horizons astrometry. We included observer dependent position-angle biases for the early speckle interferometric data. We applied biases in R.A. and decl. in Olkin et al. (2003) to account for the error in the absolute position of the Pluto system barycenter. We used a bias for some relative measurements that had the albedo model applied (Young et al. 1994; Brozović et al. 2015) and for the B. Sicardy et al. (2006, personal communication) points. These additional parameters yield more conservative fit results because they absorb potential systematic errors and help inflate the overall orbital uncertainties.

Table 7 lists ICRF Cartesian coordinates and velocities relative to the Pluto system barycenter for each satellite. We list high-precision state vectors, so they can be used as starting states for any future studies that involve orbit integration. Table 1 contains postfit residual statistics for Charon for ground-based and the HST observations. The root-mean-square (rms) of the residuals is consistent with measurement uncertainties reported by observers. The residuals range from coarse, ∼100 mas for the photographic plates of Harrington & Christy (1980) to <1 mas for the latest Gaia astrometry (Tanga et al. 2023).

Figure 1 shows Dunbar & Tedesco (1986) model fits to the Pluto–Charon mutual events lightcurves reported in Tholen (2008). This model is applicable to complex observing geometries where occultations and eclipses occur at the same time. Figure 1 demonstrates that we were able to model lightcurves throughout the changing observing geometry from 1986 to 1989. Lightcurves continue to be the source of high-quality astrometry for the pre-Hubble era. Figure 2 shows residuals for the New Horizons optical astrometry. The New Horizons LORRI CCD has image size 1024 × 1024 pixels with 102 pixel⁻¹ (Porter & Canup 2023). Charon was unresolved in the early images and all residuals are at subpixel level. The last cluster of data, starting on April 9 and ending on July 14, resolved Charon from ∼570 km pix⁻¹ to ∼2 km pix⁻¹. The rms of the residuals for the last 9 points obtained within 24 hr from the closest approach was on the order of 3 km.

Zoom In Zoom Out Reset image size

The final set of astrometry for Charon was obtained by Gaia. We show the fit residuals in Figure 3. Charon is an extended body ∼50 mas across that has non-uniform albedo distribution. If we do not model any systematic biases, such as center-of-light, center-of figure offsets, the rms of the AL residuals is 1.84 mas, with a mean of 0.16 mas (Figures 3(A) and (B)). The albedo model from Buratti et al. (2017) reduced the rms to 1.39 mas and the mean to 0.04 mas (Figures 3(C) and (D)).

Zoom In Zoom Out Reset image size

We also investigated how the correction to Pluto's planetary ephemeris influenced the residuals. We applied a linear correction to JPL's DE440 using the solution covariance as a priori (R. Park 2024, private communication). This reduced the rms for Gaia residuals to 1.05 mas with a mean of −0.26 mas (Figures 3(E) and (F)). This fit did not involve any albedo or other center-of-light, center-of-figure offsets. The planetary ephemeris correction appeared to compensate for all biases, at least to ∼1 mas level. This result is troubling, because it shows that without high-fidelity albedo model the correction to planetary orbit could be incorrect. DE440 was central to our analysis, because this orbit solution for the Pluto's system barycenter was used for the New Horizons navigation, and any changes to this solution would require a full set of Pluto's astrometry in addition to the correction implied by the Gaia astrometry of Charon. Any overhaul of planetary ephemerides is beyond the scope of this work and it also needs further verification of Gaia data on Pluto and Charon in Gaia DR5. After all, Gaia astrometry of Charon is a small portion of a long data arc and an extensive set of observations that we use in our fit. Our Gaia data set is currently dominated by the 2019 astrometry that was not part of the Gaia DR3 release and it is possible that some of these transits will get rereduced. We can clearly state that Gaia-reported submas level uncertainties in AL direction are overly optimistic. As a result of this investigation, we used only albedo-corrected Gaia astrometry (Figures 3(C) and (D)) in our final orbit fit, achieving 1.4 mas residuals.

Tables 2 through 5 contain postfit residual statistics for Nix, Hydra, Kerberos, and Styx. We subdivided Porter & Canup (2023) astrometry into years, so that we can check for any year-to-year changes in the rms. The HST data were processed to an exceptionally high quality, and we were getting ∼1 mas fit correspondence for all but Styx. This is not surprising considering that Styx is the most challenging to measure among the four satellites.

Zoom In Zoom Out Reset image size

Figure 4 shows New Horizons residuals converted to kilometers for the small satellites. Porter & Canup (2023) noted that physical resolution of the LORRI 4 × 4 (256 × 256 pixels) images is 408 pix⁻¹ that corresponds to few thousands kilometers per pixel resolution in 2015 February and a few hundred kilometers per pixel resolution in early July. This means that the small satellites all fit within a single pixel. The weighted rms of the residuals (reduced χ²) for all New Horizons astrometry are ∼0.6, suggesting that the data were weighted slightly conservatively. We were reluctant to tighten the data weights beyond what is reported in Porter & Canup (2023), because at this level, the residuals are dominated by systematics, such as camera pointing, location of Pluto's barycenter relative to spacecraft, reduction procedure, etc.

Zoom In Zoom Out Reset image size

The system mass and masses of satellites were estimated together with other dynamical parameters. We obtained the following 1σ errors from the covariance of the solution: GM_sys = 975.41 ± 0.09 km³ s⁻² and GM_Charon = 106.10 ± 0.14 km³ s⁻² (Table 8). The masses agree within their statistical errors with Brozović et al. (2015). We note that Brozović et al. (2015) inflated the formal 1σ for the system and mass of Charon from GM_sys = 975.2 ± 0.2 km³ s⁻², GM_Charon = 105.90 ± 0.3 km³ s⁻² to GM_sys = 975.2 ± 1.5 km³ s⁻², GM_Charon = 105.9 ± 1.0 km³ s⁻² in order to account for potential systematic errors in the data. Our new mass estimates, obtained with a significantly more astrometry, remained consistent with Brozović et al. (2015) even within their formal bounds. We tested our results with different data weighing schemes and found that the all mass changes remained well within 3σ. It is still possible that there are some unmodeled sources of error in the data, for which we scaled the solution uncertainties by a factor of two: GM_sys = 975.4 ± 0.2 km³ s⁻² and GM_Charon = 106.1 ± 0.3 km³ s⁻². From here, the uncertainty on the mass of Pluto is GM_Pluto = 869.3 ± 0.4 km³ s⁻² and the ratio of Charon's and Pluto's masses is m_c /m_p = 0.1220 ± 0.0003. We note that Null & Owen (1996) mass ratio of 0.124 ± 0.008 and Olkin et al. (2003) mass ratio of 0.122 ± 0.008 were already very close to the current estimate.

The masses of the small satellites were unconstrained by the data at the time Brozović et al. (2015) was published (JPL ephemeris PLU043). However, this did not affect the quality of the orbit fits, as the New Horizons Optical Navigation Team confirmed (Pelletier et al. 2016). We ran a comparison between PLU043 ephemeris and the orbital solution presented here (PLU060) for the time of the New Horizons flyby on 2015 July 14. The orbit difference for Charon was 2 km in radial and along the track directions and 5 km in out of plane direction. For Styx, Nix, Kerberos, and Hydra the maximum difference was always along the track: 244 km, 47 km, 106 km, and 81 km, respectively.

The new astrometry published in Porter & Canup (2023) allowed for mass determinations of Nix and Hydra, GM_Nix = 0.0017 ± 0.0003 km³ s⁻² and GM_Hydra = 0.0020 ± 0.0002 km³ s⁻². Our fit yields very similar values and error bars: GM_Nix = 0.0014 ± 0.0005 km³ s⁻² and GM_Hydra = 0.0020 ± 0.0003 km³ s⁻². We use a much larger set of parameters in the fit (over 200 including various biases) than Porter & Canup (2023) who fit only 28 parameters. Inclusion of more parameters inflates the covariance of the solution, so it is expected that our formal statistics differs. We also used a completely different orbit-fitting procedure, so it is encouraging to see that two very different methods give consistent results.

Our fit found no sensitivity for the masses of Kerberos and Styx. Section 3.2 shows that the formal fit 1σ uncertainties are ${\sigma }_{{{GM}}_{\mathrm{Kerberos}}}=0.17\times {10}^{-3}\,\,{\mathrm{km}}^{3}\,{{\rm{s}}}^{-2}$ and ${\sigma }_{{{GM}}_{\mathrm{Styx}}}=0.33\,\times {10}^{-3}\,\,{\mathrm{km}}^{3}\,{{\rm{s}}}^{-2}$. These are factors of ∼4 and ∼13 larger than what Porter & Canup (2023) reported for their upper bounds. In order to keep our model physically realistic, we imposed 1 g cm⁻³ density requirement on the satellites, which yielded GM_Kerberos = 6 × 10⁻⁵ km³ s⁻² assuming 6 km radius, and GM_Styx = 4 × 10⁻⁵ km³ s⁻² assuming 5.25 km radius (Porter et al. 2021). We compared the nominal orbit solution with the one where both satellites were kept massless, and we noted that the most significant effect was ∼3 km downtrack change in the orbit of Hydra. This was for the time span for which we have astrometry, 2005–2020. In other words, Kerberos and Styx are so small, that they introduce negligible perturbations on the rest of the system, at least over the period of time for which we have the data.

Table 9 compares current densities with the last two orbit fits for the Pluto system. The values for Pluto and Charon in Brozović et al. (2015) are consistent with the current work. Similarly, densities for Nix and Hydra are consistent with Porter & Canup (2023). The density limits for Nix and Hydra were not well constrained in Brozović et al. (2015) because of the order-of-magnitude masses and sizes estimates. We can tentatively put an upper bound of 2.8 g cm⁻³ (1σ) on the density of Kerberos, but the mass/density of Styx is cannot be constrained with the current data.

We summarize 400 yr (from 1800 January 1 to 2200 January 1) of orbit integration by calculating mean values of the elements a, e, i, and P of the osculating ellipse associated with the integrated orbit. The outline of the procedure can be found in Smart (1960), page 58. The Cartesian states were sampled once per day during this interval. Osculating period is the average sidereal period calculated from instantaneous values of $2\pi \sqrt{(}{GM}/{a}^{3})$. The error bars in Table 10 reflect variations in two-body elements, as opposed to the fit uncertainties. Charon's orbit is nearly circular and its osculating period varies less than 0.2 s. The small satellites are strongly perturbed by Charon and by each other to some extent, hence their osculating elements vary in time. Styx and Nix move in the same plane as Charon, while Kerberos and Hydra have small inclinations, 0.4 deg and 0.3 deg, respectively. Figure 5 shows small variations in the osculating ellipse elements for Kerberos.

Zoom In Zoom Out Reset image size

Table 10. Mean Keplerian Osculating Orbital Elements for the Satellites of Pluto

Element	Charon	Styx	Nix	Kerberos	Hydra
a(km)	${19595.764}_{-0.008}^{+0.007}$	${43163}_{-1085}^{+677}$	${49340}_{-666}^{+419}$	${58272}_{-366}^{+268}$	${65175}_{-247}^{+195}$
e	${1.61}_{-0.01}^{+0.01}\times {10}^{-4}$	${0.025}_{-0.020}^{+0.016}$	${0.015}_{-0.014}^{+0.010}$	${0.010}_{-0.010}^{+0.007}$	${0.009}_{-0.009}^{+0.007}$
i (deg)	⋯	${0.049}_{-0.006}^{+0.006}$	${0.022}_{-0.002}^{+0.002}$	${0.422}_{-0.006}^{+0.006}$	${0.280}_{-0.003}^{+0.003}$
P (days)	${6.387221}_{-0.000003}^{+0.000005}$	${20.89}_{-0.79}^{+0.49}$	${25.51}_{-0.51}^{+0.33}$	${32.75}_{-0.30}^{+0.23}$	${38.74}_{-0.22}^{+0.17}$

Note. We list mean osculating elements where a is the semimajor axis, e is the eccentricity, i is the inclination with respect to Charon's orbit pole (see Section 3.4), and P is the orbital period. Elements for Charon are Plutocentric and elements for the small satellites are barycentric. The upper and lower bounds are based on sampling 400 yr of orbit integration (from 1800 January 1 to 2200 January 1) once per day.

Download table as:  ASCIITypeset image

An alternative way to descriptively summarize integrated orbits is to apply a precessing ellipse theory as in Jacobson & Owen (2004). We sampled 400 yr of orbit integration once per day and used these positions as "measurements." The objective was to find a single set of six Keplerian elements and three rates that minimize rms of the residuals projected along the track, in radial, and out of plane directions. The optimization procedure used the least-squares method. The resulting precessing ellipse captures secular behavior of the integration, while the periodic terms are not modeled (Figure 6). The true position of Kerberos can be ∼500 km different along the track than what a simple analytical theory predicts. However, on average, the ellipse model is good. Table 11 summarizes the results. The nodal regression period for Kerberos is ∼3300 days, and it is ∼5100 days for Hydra, similar to its apsidal precession period.

Zoom In Zoom Out Reset image size

Table 11. Mean Precessing Ellipse Model Orbital Elements for the Satellites of Pluto

Element	Charon	Styx	Nix	Kerberos	Hydra
a (km)	19595.764	42410	48689	57749	64719
e	1.61 × 10−4	⋯	⋯	⋯	0.006
i (deg)	⋯	⋯	⋯	0.42	0.28
n (deg day−1)	56.362530	17.855	14.484	11.191	9.424
P (days)	6.387222	20.16	24.85	32.17	38.20
$\dot{\varpi }$ (deg day−1)	2.5 × 10−6	⋯	⋯	⋯	0.070
$\dot{{\rm{\Omega }}}$ (deg day−1)	⋯	⋯	⋯	−0.11	−0.071

Note. We list mean elements according to a precessing ellipse model. This ellipse represents best fit to 400 yr of orbit integration. We only list eccentricities and inclinations if they are larger than 1 × 10⁻⁵ and 1 × 10⁻³ deg respectively. In addition to parameters in Table 10, we also estimated n, the mean longitude rate, $\dot{\varpi }$, the apsidal precession rate, and $\dot{{\rm{\Omega }}}$ the nodal precession rate. Elements for Charon are Plutocentric and elements for the small satellites are barycentric.

Download table as:  ASCIITypeset image

We were interested in examining the influence of the solar torques on the Pluto–Charon orbit plane. We used Ward (1975) and Ward & Hamilton (2004) analytical expression to estimate rate, $\dot{\psi }$, at which the Charon's orbit normal is precessing about the Pluto's heliocentric orbit normal

$$\begin{eqnarray}&&\dot{\psi }=-\displaystyle \frac{3}{2}\displaystyle \frac{{n}_{\odot }^{2}}{s}\left(\displaystyle \frac{{J}_{2}+q}{\gamma +l}\right)\cos (Q)\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&q=\displaystyle \frac{1}{2}\left(\displaystyle \frac{{{GM}}_{\mathrm{Charon}}}{{{GM}}_{\mathrm{Pluto}}}\right){\left(\displaystyle \frac{{a}_{\mathrm{Charon}}}{{R}_{\mathrm{Pluto}}}\right)}^{2}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&l=\displaystyle \frac{{n}_{\mathrm{Charon}}}{{s}_{\mathrm{Pluto}}}\left(\displaystyle \frac{{{GM}}_{\mathrm{Charon}}}{{{GM}}_{\mathrm{Pluto}}}\right){\left(\displaystyle \frac{{a}_{\mathrm{Charon}}}{{R}_{\mathrm{Pluto}}}\right)}^{2}.\end{eqnarray} \tag{ 3 }$$

In this expression, Q is Charon's orbit obliquity with respect to Pluto's heliocentric orbit, n_Charon is orbital mean motion of Charon, n_☉ is Pluto's heliocentric mean motion, s_Pluto is Pluto's rotation rate (same as n_Charon), R is Pluto's equatorial radius, a_Charon is semimajor axis of the orbit of Charon, GM_Pluto and GM_Charon are masses of Pluto and Charon, J₂ is Pluto's degree 2 zonal gravity harmonic, and γ is Pluto's normalized polar moment of inertia. For the case of Charon, the q and l terms in Equations (2) and (3) dominate over any reasonable physical values of J₂ and γ in the Equation (1). We assumed J₂ = 0, consistent with Pluto being a sphere, and γ = 0.4, consistent with uniform density model.

Values for semimajor axis, mean longitude rate, and Pluto's and Charon's masses can be found in Tables 11 and 8. We used 30,000 yr of JPL's DE441 planetary ephemeris (Park et al. 2021) to find mean longitude rate for Pluto's heliocentric motion (n_☉ = 145.75 deg century⁻¹) and mean obliquity of Charon's orbit (Q = 119.7 deg). For completeness, DE441 yields mean R.A. and decl. of the Pluto's heliocentric pole of 314.087 deg and 66.531 deg, apsidal regression rate of −0.00897 deg century⁻¹, and nodal regression rate of −0.00181 deg century⁻¹. The mean sidereal period of Pluto is 246.99 yr. We obtained a precession rate $\dot{\psi }=3.79\,\times {10}^{-3}\,\deg \,{\mathrm{century}}^{-1}$, yielding a precession period of 9.5 million years. Consequently, Charon's orbit pole precession can be ignored in pole models spanning several thousand years.

We calculated osculating Keplerian R.A. and decl. of Charon's orbit pole for the duration of 5000 yr (Figure 7). The Sun influence is visible as a slope in R.A. and decl. There are also variations in R.A. and decl. that appear consistent with Pluto's orbital period, ∼247 yr. R.A. of the pole can be fit as R.A. = 133.0073 + 0.0036T deg, where T is Julian centuries from the epoch J2000. Decl. of the pole can be fit as decl. = –6.2449 + 2.5 × 10⁻⁵ T deg. These first-order polynomials capture the pole variations to within ∼2 millidegrees over 5000 yr.

Zoom In Zoom Out Reset image size

Jacobson et al. (2012) showed that linear mapping of the covariance matrix obtained from a least-squares fit represents a method for placing a reasonable bound on the satellite's position at some later time. We used this method to evaluate how the positional uncertainties grow for each of the five satellites until 2050 (Figure 8). Unlike the uncertainties listed in Table 10 that represent differences between an analytical theory and true orbits, these uncertainties reflect the fit quality that is based on quantity and quality of measurements.

Zoom In Zoom Out Reset image size

Overall, all orbits are well constrained and none of the satellites is in danger of being lost for the foreseeable future. Styx orbit has the highest uncertainty because this satellite is challenging to measure in the HST images due to its faint magnitude and proximity to the bright Pluto–Charon binary. The dominant uncertainty is always in in-orbit direction and the fit to the slopes yields the mean motion uncertainties: 9 × 10⁻⁷ deg day⁻¹ for Charon, 2 × 10⁻⁵ deg day⁻¹ for Styx, 4 × 10⁻⁶ deg day⁻¹ for Nix, 7 × 10⁻⁶ deg day⁻¹ for Kerberos, and 2 × 10⁻⁶ deg day⁻¹ for Hydra. True uncertainties could be larger than this due to unknown systematic errors.

If we take the average osculating periods from Table 10, the satellites have the following ratios with respect to Charon: 3.27:3.99:5.13:6.07. If we take the periods obtained with the precessing ellipse theory (Table 11), the ratios are: 3.16:3.89:5.04:5.98. This shows that the definition of elements matters for the case of perturbed orbits, and that their interpretation (e.g., proximity to resonance) should be treated with caution.

However, taken at face value, we can note Charon–Kerberos proximity to the fourth-order inclination resonance: ${n}_{\mathrm{Charon}}\,-5{n}_{\mathrm{Kerberos}}+4{\dot{{\rm{\Omega }}}}_{\mathrm{Kerberos}}=-0.03\,\deg \,{\mathrm{day}}^{-1}$. A change of 7% in the nodal precession rate would put Kerberos in a resonance. The current estimate for the nodal precession period is about nine years, and the data arc is only about 12 yr long with sparse astrometry until 2012. Thus, there is a room for improvement on the nodal precession rate. The more likely scenario is that the resonance existed in the past.

We have similar case with Charon and Hydra. Hydra has both small eccentricity and inclination, and the closest combination of the rates to a resonance that we could find is: ${n}_{\mathrm{Charon}}-6{n}_{\mathrm{Hydra}}+4{\dot{\varpi }}_{\mathrm{Hydra}}+{\dot{{\rm{\Omega }}}}_{\mathrm{Hydra}}=0.03\,\deg \,{\mathrm{day}}^{-1}$. The nodal and apsidal precession periods are on the order of 14 yr, and it is possible that these numbers will get revised with more data.

Showalter & Hamilton (2015) reported that Nix, Hydra, and Kerberos reside in a Laplace-like resonance. They identified 3n_Styx − 5n_Nix + 2n_Hydra = − 7 ± 1 × 10⁻³ deg day⁻¹ and suggested that the resonant argument, Φ = 3λ_Styx − 5λ_Nix + 2λ_Hydra, librates around 180 deg. Here, λ refers to mean longitude, measured from the ascending node of the Pluto–Charon orbital plane on the ICRF equator. Showalter & Hamilton (2015) noted that the resonant argument Φ decreased from 191 deg to 184 deg during 2010–2012, and that this was too coincidental to be observed by chance. However, their long-term integration showed that the libration depends on the masses of Nix and Hydra that were not well constrained at that time (Brozović et al. 2015). Showalter & Hamilton (2015) showed in their Figure 2 that the resonant argument was circulating for the case of moderate masses of Nix and Hydra (GM_Hydra = 0.0032 km³ s⁻² and GM_Nix = 0.0044 km³ s⁻²) and that the libration condition needed several tens of percent higher mass of Hydra (Nix was already increased by 50% from Brozović et al. 2015).

We revisited this result in the light of the current orbits and masses. We obtained 3n_Styx − 5n_Nix + 2n_Hydra = − 7.3 × 10⁻³ deg day⁻¹ based on the rates from Table 11. We also calculated resonant argument based on 400 yr of osculating longitudes in Figure 9. The argument shows circulation with ∼140 yr period. The argument indeed goes through 180 deg around 2010–2012 time frame, but this does appear to be a coincidence. The lack of resonance is not surprising given the relatively "light" masses for Nix and Hydra.

Zoom In Zoom Out Reset image size

We have seen in Section 3.2 that the Charon's mass was known to better than 1% even prior to the additional optical navigation measurements from New Horizons. The New Horizons astrometry further reduced this mass uncertainty by a factor of three. Regarding the small satellites of Pluto, we only have constraints on Nix and Hydra, and the principal source of mass-sensitive data are two decades of HST astrometry. The New Horizons astrometry contributed to reducing the GM uncertainty for Nix by ∼15%.

We simulated milliarcsec precision astrometry for all small satellites for the next 20 yr hoping to reveal the masses of Kerberos and Styx, but their masses remained elusive. The uncertainties were still several times larger than their expected values assuming ∼1 g cm⁻³ density. This is not surprising considering that the orbital differences between our nominal solution and the one where Kerberos and Styx were kept massless were small. By 2050, the in-orbit differences grow to 5 km for Styx, 6 km to Nix, 3 km for Kerberos, and 7 km for Hydra. However, more high-precision astrometry would improve the uncertainties on masses of Nix and Hydra by a factor of two, and it would also further improve their state vector uncertainties.