The orbit and density of the Jupiter Trojan satellite system Eurybates-Queta

We report observations of the Jupiter Trojan asteroid (3548) Eurybates and its satellite Queta with the Hubble Space Telescope and use these observations to perform an orbital fit to the system. Queta orbits Eurybates with a semimajor axis of $2350\pm11$ km at a period of $82.46\pm0.06$ days and an eccentricity of $0.125\pm0.009$. From this orbit we derive a mass of Eurybates of $1.51\pm0.03 \times 10^{17}$ kg, corresponding to an estimated density of $1.1\pm0.3$ g cm$^{-3}$, broadly consistent with densities measured for other Trojans, C-type asteroids in the outer main asteroid belt, and small icy objects from the Kuiper belt. Eurybates is the parent body of the only major collisional family among the Jupiter Trojans; its low density suggests that it is a typical member of the Trojan population. Detailed study of this system in 2027 with the Lucy spacecraft flyby should allow significant insight into collisional processes among what appear to be the icy bodies of the Trojan belt.


INTRODUCTION
The Jupiter Trojan asteroid (3548) Eurybates is the parent of the largest known collisional family in the Trojan population (Brož & Rozehnal 2011;Vinogradova 2015;Nesvorný gradova 2015;Nesvorný et al. 2015;Rozehnal et al. 2016); the Eurybates family is the only major family, though a few smaller ones have been suggested. Interestingly, Eurybates and its collisional family members are also among the least red objects in the entire Trojan population (Fornasier et al. 2007;). This group is often classified as C-type asteroids in the DeMeo & Carry (2013) taxonomy, for example, compared to the P-and D-types found among the rest of the Jupiter Trojans. For these objects with essentially featureless spectra out to 2.5 µm and poorly measured albedos, however, it is probably more informative to simply consider that the Jupiter Trojan asteroids are bifurcated in color, with a red population (which generally maps to the D-type classification) and a less-red population (which generally maps to the P-type classification) (Emery et al. 2011;Wong et al. 2014) and that the distribution of colors of Eurybates and its family overlap with the less-red colors but skew blueward (De Luise et al. 2010;).
The small number of collisional families among the Trojans could perhaps suggest compositional and strength differences between the Jupiter Trojan asteroids and the main belt asteroid population (Levison et al. 2009). One possibility for the unique family and atypical colors of Eurybates could be that this object is an interloper, with a composition more similar to main belt asteroids than to the potentially ice-rich Jupiter Trojan population.
Like collisional families, known satellites also appear rare among the Jupiter Trojan asteroid population, with Eurybates only the fourth Jupiter Trojan to have had a satellite directly imaged (Noll et al. 2020), though it is clear that some of this paucity is an observational bias owing to their greater distance than main belt asteroids. It seems likely that, for Eurybates, the existence of a satellite and of a collisional family are related, and that both trace back to a catastrophic impact at some earlier time (Durda et al. 2004). The existence of this satellite allows us to measure the density of Eurybates and potentially help understand the the critical question of the relationship of Eurybates to the remainder of the Jupiter Trojans and to the outer solar system.

OBSERVATIONS AND DATA REDUCTION
The Eurybates system was observed 13 times between 12 Sep 2018 and 12 Feb 2021 (Table 1). The satellite was first detected in the initial pair of observations on 12 and 14 Sep 2018 and followup to confirm the existence of the satellite and determine the orbit began on 11 Dec 2019. Details of the observing strategy are given in Noll et al. (2020). The first successful re-aquisition occured on 3 Jan 2020, just before Eurybates became unobservable due to solar avoidance. The satellite was again detected at the first opportunity post-solar avoidance on 19 July 2020 and then was not detected at the next two attempts. The next detection occurred on 12 October 2020, by which time the orbit solution had converged sufficiently that the next 4 attempts were all successful.
Observations were obtained in a variety of passbands, depending on the expected brightness of Eurybates at the time of observation. The initial observations were obtained with F555W. For followup observations we either used F350LP, when Eurybates was more distant from the Earth and thus fainter, or a combination of F606W and F814W when Eurybates was brighter and closer. In all of the followup observations the primary is saturated in order to get better signal-to-noise on the satellite (additionally unsaturated short exposures were also obtained, but Queta is not detected in these so they were not used in this analysis).
For each set of observations, data reduction proceeds identically. All data reduction steps are performed on the flat-fielded images pro- vided by the HST pipeline (the "*_flt.fits" files from the archive); no use is made of the re-sampled geometrically rectified images. First, the saturated regions and any regions affected by bleeding of charge from these saturated regions along the column are automatically masked; in practice we conservatively mask everything at approximately 70% of the expected saturation level and higher. Next, each image is manually inspected and pixels with obvious cosmic ray hits are added to the mask. Low level cosmic rays are always present and likely go unmasked. These events will contribute to our final uncertainties.
We fit the positions of Eurybates and Queta in two steps: first we determine the point-spread function (PSF) of the telescope and the precise position of Eurybates, and then we fix the previous two parameters and fit the brightness of Eurybates and the brightness and position of Queta to a more localized set of data. This two-step process allows the final astrometry of Queta to be unaffected by the details of the interior (or exterior) part of the PSF, where the largest discrepancies between the model and data often exist. In the first step of our fitting, we create a model point-spread function (PSF) using TinyTim (Krist 1993) for the appropriate filter and pixel position. We find that the true PSF of the observations is more extended than any of the model PSFs, so we convolve the PSF with a symmetric two-dimensional Gaussian function. We perform a least-squares minimization of all of the unmasked pixels using the IDL implementation of MPFIT (Markwardt 2009), allowing the position of Eurybates, the brightness of Eurybates, and the width of the Gaussian kernel to vary. Even with saturated 1.50 a The positional uncertainty, assumed to be identical in RA and dec.
b The distance from the predicted position divided by the uncertainty observations, we find that we can fit the position of Eurybates quite reproducibly (as demonstrated by the small scatter in final calculated value of the offset position from Eurybates to Queta). We next mask all pixels either within 12 or beyond 25 pixels of Eurybates. This step leaves us with an annulus that includes Queta and the wings of the PSF of Eurybates. In this second step, we fix the PSF and the position of Eurybates to that previously determined, and perform a least-squares fit allowing the position of Queta to vary, as well as the flux of Eurybates and Queta. We average all of the observations from a visit to determine the astrometric offset at the mean time of the observation. It is difficult to estimate uncertainties from such a non-linear process a priori so we use the dispersion of the results from multiple measurements during a single HST visit as a proxy for uncertainty of the mean. We conservatively assume that the uncertainties on the measurement in both dimensions of the detector are the same and calculate the uncertainty as the combined standard deviation of the mean in both dimension. In two cases (14 Sep 2018, 3 Jan 2020) a cosmic ray landed close enough to Queta on one of the images to affect the measurement, so that individual observation was discarded from the analysis of that visit. All pixel positions and uncertainties are converted to angles on the sky using the distortion lookup tables provided in the image headers. Uncertainties range from 3 mas for observations where the satellite is wellseparated and 8 observations were obtained to 16 mas for the observation where Queta was closest to Eurybates and only 3 non-cosmic ray affected observations were obtained. The final astrometry is presented in Table 2.

ORBIT FIT
We initially had no knowledge that Earth was passing through the orbital plane of the system and thus the viewing geometry was changing  Table 2. In most cases the uncertainty in the measured position is smaller than the size of the symbol. Projections of the orbit at three separate epochs are shown. The red shows the projection at the epoch of the first two detections (also in red) The yellow projection shows the orbit at the moment of the 3 Jan 2020 re-acquisition of the satellite (yellow point). The predicted position of Queta at the times of the two non-detections previous to the re-acquisition can be seen in grey along the projection of the yellow orbit at positions 3 and 4. The blue orbit projection shows the nearly-edge on configuration on 24 Nov 2020 with the three last detections shown as blue points. The purple predicted points are from Jul and Oct 2020.
rapidly. Much effort was thus expended trying to understand the geometry of the orbit during the period from the first re-acquisition on 3 January 2020 until a good solution -accounting for all of the detections and also the non-detections -was finally obtained after the 19 Sep 2020 nondetection, Here, however, we simply detail the final orbit fitting.
The orbit is initially fit with a least-squares minimization, again using MPFIT, where all parameters are allowed to vary. We fit for semimajor axis, a, orbital period, P , inclination (with respect to the ecliptic), i, the eccentricity vector [e sin ω, e cos ω], the longitude of ascending node, Ω, and the mean anomaly at the epoch of the first observation, M . Neither orbital evo-lution nor external perturbations from the Sun are considered, though we revisit this issue later. Figure 1 shows the observed astrometric positions of Queta as well as the positions predicted from the least-squares fit. The least-squares fit has a χ 2 of 8.2 (note that we remove the nondetections from this fit, but verify that all of the predicted positions at the time of non-detection would indeed have been undetectable). Given the 9 detections (with both ∆RA and ∆dec) and the 7 fit parameters, this value corresponds to a reduced χ 2 of 0.75, suggesting that our approach to the uncertainties is not unreasonable.
To better understand the range of uncertainties and correlations among the parameters, we take the initial least-squares fit as the starting position in a Markov Chain Monte Carlo (MCMC) model. We use the Python package emcee (Foreman-Mackey et al. 2013) which implements the Goodman & Weare (2010) affineinvariant MCMC ensemble sampler, taking uniform priors in the same fit parameters as used above.
While the least-squares fit shows that the non-detections came at times when the satellite was too close for detection, to ensure that the MCMC properly penalizes orbits where these points should have been observed, we include the times of the the non-detections in the likelihood calculation and impose a likelihood penalty if the predicted position is more than 0.45 arcseconds from Eurybates (in practice the predicted positions are all well inside this radius so they have little effect on the final results). We run a series of 100 separate chains and collect 10,000 samples, approximately 100 times longer than the autocorrelation timescale of the chain, as determined after-the-fact. For our final results we discard the initial 2000 samples to ensure that our samples are independent of the initial parameters and take every tenth sample from the chains to select a statistically independent set of samples of the distribution. These samples approximate the posterior dis-tribution of the orbital elements. In these samples, all parameters have approximately Gaussian marginalized distributions, with some correlation between different parameters. As the uncertainties are small, however, in Table 3 we simply give the median and 16th and 84th percentiles as the ±1σ range of the parameters. From the semimajor axis and period we calculate a mass of Eurybates of 1.51 ± 0.03 × 10 17 kg (the contribution to the system mass of the satellite is negligible). The volume of Eurybates has large uncertainties. WISE observations and modeling put the effective diameter at 63.9 ± 0.3 km. Correcting the WISEassumed absolute magnitude of H V = 9.8 to the value of 10.01±0.02 measured by  changes the retrieved diameter by an amount smaller than the listed uncertainty (though it does drop the best-fit albedo from 0.052 to 0.043). Systematic errors in the WISE diameters are difficult to assess, but they include difficulties in the color corrections for cold targets such as Jupiter Trojans (Wright et al. 2010), uncertainties owing to the simplistic thermal modeling used (Grav et al. 2012), and the inability to account for concave surfaces. We thus -perhaps conservatively -assign 10% uncertainties to the WISE diameter. An ongoing occultation campaign (Buie et al. 2020) could provide significantly more accurate results prior to the Lucy flyby. With our assumed volume of 1.4 ± 0.4 × 10 14 m 3 , the bulk density of Eurybates is 1.1 ± 0.3 g cm −3 . The results show that Queta is on a low-eccentricity retrograde (with respect to the ecliptic) orbit with a moderately long period around a lowdensity primary.

ORBITAL EVOLUTION
Queta is widely separated from Eurybates, with a semimajor axis that is 0.104 times that of the Hill radius of Eurybates, so the satellite orbit could be subject to orbital evolution induced by the Sun. We investigate this possibility though a long-term dynamical integration in which we use a Wisdom-Holman integrator (Wisdom & Holman 1991) to explicitly integrate the full Sun-Eurybates-Queta system. As an illustration, Figure 2 shows the evolution of Queta's orbit about Eurybates for the next 3000 years.
To better understand the evolution, the orbital parameters are shown with respect to the plane of the orbit of Eurybates, rather than in the ecliptic coordinate system of Table 3. The large oscillations in eccentricity are driven by Kozai forcing driven by the gravitational effects of the Sun (Kozai 1962). In particular, the solar perturbation causes coupled oscillations between eccentricity (e) and incli-nation (i) such that the magnitude of the component of angular momentum in the direction perpendicular to the plane of Eurybates' heliocentric orbit, l z = a √ 1 − e 2 cos (i), is constant. In this orbital integration, the fact that the argument of pericenter (ω) is circulating indicates that Queta is not in the Kozai resonance itself. Nonetheless, the observed evolution is induced by the resonance islands at ω = 90 • and 270 • , which force the eccentricity to become large as ω sweeps by. The temporal evolution of Queta's orbit, shown in a reference system defined by the plane of the heliocentric orbit of Eurybates, so as to illustrate the Kozai-induced evolution. At the top is the argument of pericenter, ω, which would librate near 90 • or 270 • for an object in the Kozai resonance, rather than circulate, as seen here. The middle panel shows eccentricity, e, with Queta's current observed value of 0.125 shown by the horizontal line. The inclination, i, with respect to the heliocentric orbit of Eurybates, is shown on the bottom, with the current value of 139.4 • shown.
Assuming no other perturbations, the current orbit should oscillate in eccentricity from 0.11 to 0.37 and in inclination from 142 to 146 • (in the plane of the heliocentric orbit of Eurybates) with a timescale of about 500 years. The expected orbital evolution over the 29 month span of our observations is smaller than the uncertainties on our orbital elements, so our orbital  The contours shown are the 50%, 90%, and 99% confidence intervals. These data are shown in a coordinate system centered on Eurybates, with the Sun on the −X axis. Lucy will be targeted to cross the Eurybates-Sun line, thus we define a coordinate system where Lucy always remains on the X-Y plane. The blue vector shows the path of the spacecraft during 500 seconds of the encounter.
fits, which do not include the effect of the Sun, remain valid. With continued monitoring of the orbit, however, the orbital elements should get sufficiently precise that a full fit including the gravitational effects of the Sun will have to be performed. Additional perturbations to the orbit could occur if there are other satellites in the system. Indeed, an additional satellite could disrupt the Kozai resonance, thereby stifling the large eccentricity oscillations. Both (87) Sylvia and (107) Camilla are outer main belt asteroids associated with collisional families that have at least two satellites. In the Kuiper belt, Haumea is the parent of the only known collisional family and also has at least two satellites. While an interior satellite would be difficult to directly image prior to the Lucy encounter, if such a satellite exists, its dynamical effects could perhaps eventually be detected by perturbations to the orbit of Queta or it could be detected through an occultation.

LUCY FLYBY
NASA's Lucy mission will perform a close flyby of Eurybates on August 12, 2027. At closest approach, Lucy will fly within ∼ 1000 km of Eurybates, well within the orbit of Queta. In order to determine the location of Queta at the time of the encounter, we take 3000 randomly selected samples from our MCMC chain and propagate the orbits forward until the time of the encounter. The span of this extrapolation is sufficiently large that we include the Sun's effects in this integration, though we find that this step still makes little difference in the final predictions. The probability distribution function of the predicted position of Queta is shown in Figure 3.
The closest approach distances between Queta and Lucy range from 116 to 1262 km. The average value of the close approach distance is 513 km, while the mode of the distribution is 460 km. Unfortunately, the geometry of the en-counter is such that there is only a 3% chance that Queta's phase angle will be less than 90 • at the moment of close approach, meaning that the Lucy spacecraft will most likely pass over the unilluminated hemisphere. Queta will, however, be able to be imaged as the spacecraft approaches the system.

DISCUSSION
The low density of Eurybates is consistent with the 0.8 g cm −3 of the Patroclus-Menoetius binary estimated from the orbit-derived system mass and a volume determined from stellar occultations (Buie et al. 2015). The system mass of Hektor is well known from the orbit of its satellite (Marchis et al. 2014), but estimates of bulk density vary from 1.0 to 2.5 g cm −3 due to differences in volume estimates for the bilobed primary (Descamps 2015). Outside of the Trojans, Eurybates' density falls in the bottom of the range determined for nine C-complex objects in the main asteroid belt where the system mass has been determined from satellite orbits or spacecraft flyby. Bulk densities determined for these systems range from 0.7 -1.8 g cm −3 with a mean of 1.4 g cm −3 (Thomas et al. 1999;Margot et al. 2015). The effects of porosity on the densities of these objects is unknown, but the low density of Eurybates suggests that it and other low-density asteroids could be derived from the icy population of planetesimals that formed beyond the giant planets for which similar low densities are measured (Morbidelli et al. 2005).
The low density of Eurybates suggests that it is consistent with being a typical member of the Jupiter Trojan population, rather than a unique interloper. If the original Eurybates parent body was indeed a typical Jupiter Trojan asteroid then it is plausible that the impact itself has led to the bluer-than-usual colors for Eurybates and the family members. In the hypothesis of Wong & Brown (2016), catastrophic collisions will remove the original irradiation crust of a Jupiter Trojans, and the new irradiation crust formed at 5 AU will be significantly less red than the original. Family members formed in this manner would suffer the same fate and have similarly less-red colors. This hypothesis was originally developed to explain the two color populations in the Jupiter Trojans and the observations that the smaller, thus more collisionally active, objects in the Jupiter Trojan population are predominantly less red (Wong et al. 2014;Wong & Brown 2015); its consistency with the properties of Eurybates and its family is encouraging. The ability to study this system in detail with the Lucy spacecraft will allow unique insights into the catastrophic impacts and their aftermath in the icy bodies of the Jupiter Trojan population.