The Case for a Large-scale Occultation Network

Authors, notably including Trujillo & Sheppard (2014), Batygin & Brown (2016), and Malhotra et al. (2016), have raised the possibility that the solar system contains an as-yet-undetected planet orbiting far beyond Neptune (typical models often invoke an M ∼ 10 M_⊕ body on an eccentric orbit with P ∼ 2 × 10⁴ yr). Evidence for such a "Planet Nine" is indirect, however, and rests primarily on the larger-than-expected degree of orbital alignment among the most distant KBOs. Additional clues are provided by the observed population of Centaur-like objects with high inclinations that are difficult to explain through alternative mechanisms (Batygin & Brown 2016). Numerical simulations by Holman & Payne (2016b), Millholland & Laughlin (2017), Becker et al. (2017), and Batygin & Morbidelli (2017; among others) have steadily refined the dynamical evidence and now point toward a planet with m ∼ 5−10 M_⊕ with semimajor axis a ∼ 400–800 au. Parameter ranges based on current constraints are listed in Table 1 (Batygin et al. 2019), along with the assigned values used throughout this work.

Direct optical searches for Planet Nine in reflected sunlight are rendered difficult by the r⁻⁴ decrease in brightness with distance, and they become extremely challenging for V ≳ 23–24, especially in regions that are crowded with stars. Directed searches at infrared wavelengths (e.g., Meisner et al. 2018) for self-luminous trans-Neptunian planets are stymied by the detection limits of wide-field space-based surveys (e.g., WISE, Wright et al. 2010; NEOWISE, Mainzer et al. 2011) unless a planet exists that is much larger and warmer than expected. Because Planet Nine's emission is expected to peak at submillimeter wavelengths (50–100 μm), cosmic microwave background (CMB) surveys provide an alternative method to search for Planet Nine through its parallactic motion (Cowan et al. 2016). However, only a few current and upcoming CMB surveys have the requisite sensitivity (Cowan et al. 2016), and the effectiveness of this search method relies upon assumptions about Planet Nine's flux while also requiring a specific observing mode that repeatedly maps the same large regions of the sky at a cadence of a few months. Synergies between multiple detection methods may exist (Baxter et al. 2018), but a definitive assessment remains elusive. Furthermore, in the case that a ninth planet does not exist, it is difficult to rule out its presence simply through the lack of a direct detection.

Fortunately, however, indirect detection via the planet's tidal gravity signature may be feasible. Shortly after Batygin & Brown (2016) issued their detailed prediction of the Planet Nine orbit, Fienga et al. (2016) combined a high-resolution solar system dynamical model with telemetry from the Cassini mission to assess whether the data and model are consistent with the presence of Planet Nine at different locations along Batygin & Brown's (2016) fiducial orbit. They found no constraint on the planet near apoastron while simultaneously determining that its presence near periastron is inconsistent with the data. This work has since been revisited by Folkner et al. (2016), who found that, under certain assumptions regarding the mass distribution in the Kuiper Belt, perturbing bodies larger than 10 M_⊕ at 1000 au are disallowed. Holman & Payne (2016a, 2016b) further analyzed the proposed Planet Nine's effects on the known solar system ephemeris, concluding that the models are best fit where Planet Nine is placed at or near apoastron in its orbit.

The instantaneous tidal differential acceleration a_t induced by a perturbing object of mass m at a distance d from the Sun is given by

$$\begin{eqnarray}&&{a}_{t}=\displaystyle \frac{{x}_{B}{Gm}}{{d}^{3}},\end{eqnarray} \tag{ 1 }$$

where x_B is the effective baseline across which the tidal acceleration is applied. For an as-yet-undetected 7.5 M_⊕ object at d = 810 au (the apoastron Sun–Planet Nine distance corresponding to the Planet Nine properties adopted from Table 1) and measured with a baseline x_B ∼ 10 au appropriate to Jupiter's Trojan asteroids, a_t ∼ 3 × 10⁻¹³ cm s⁻², corresponding to a displacement ${dx}\sim \tfrac{1}{2}{a}_{t}{t}^{2}\sim 30\,{\rm{m}}$ over a 5 yr timescale.

A net displacement of order dx ∼ 30 m incurred over 5 yr may, at first glance, appear to be insignificant, but it has the same order of precision to which solar system ephemerides are currently determined. For example, the range residuals to the Cassini orbiter's telemetry reported by Fienga et al. (2016) are of order dx ∼ 10–100 m, whereas the Martian range has been determined to within 3.6 m using just over 1 yr of Mars Express data and the lunar range to within 2 cm using 10 yr of Apollo retroreflector data from the Apache Point Observatory and 2 yr of infrared lunar laser ranging data from the Grasse station (Viswanathan et al. 2017, 2018).

Occultation monitoring can be used to measure the accumulated effects of the tidal acceleration imparted by massive outer solar system objects, including Planet Nine. To fix ideas, we investigate the use of the aggregate Jovian Trojan populations, located at the L4 and L5 Jupiter–Sun Lagrange points, as a probe of tidal gravity. This population provides an abundance of available objects to track—a significant statistical sample that permits triangulation to distinguish perturbations of the Earth's orbit from those of the Trojans' orbit—as well as a relatively long baseline over which the differential tidal acceleration from an external perturber can be measured (∼10 au across a Jovian Trojan orbit). Over 7000 Jovian Trojans are currently cataloged in the JPL Small-Body Database,³ and this number is anticipated to surge with the advent of LSST. Both Jovian Lagrange points are well populated, with number asymmetry N(L₄)/N(L₅) = 1.34 reported from WISE/NEOWISE observations (Marzari et al. 2002; Grav et al. 2012). Hereafter in Section 3 we use the terms "Trojan" and "Jovian Trojan" interchangeably.

An accurate measurement of tidal acceleration must combine precise astrometry with a rigorous dynamical model of the full solar system to integrate the trajectories of occulting objects. Together, the combination of measurements and model enables the prediction of future sky positions that can be compared with future observations. The requisite stellar astrometric precision is provided by the Gaia mission, while solar system ephemeris models (e.g., Giorgini et al. 1996; Fienga et al. 2014; Folkner et al. 2014; Pitjeva & Pitjev 2018) combined with occultation measurements provide ever-stricter constraints on the projected positions of foreground asteroids over time.

The aforementioned Gaia mission is obtaining accurate sky positions and sky motion models for N ∼ 10⁹ stars with G < 20.0, where Gaia's G filter is similar to the standard V-band filter (Gaia Collaboration et al. 2016). For this large aggregate of stars, the original mission specifications mandated that 50% of sources have G < 19 and 2.27% (∼26 million stars) have G < 14.5. The actual astrometric precision to which Gaia's observed stars can be fixed varies inversely with brightness, with design expectations ranging from ω = 7 μas at V = 10 to ω = 12–25 μas at V = 15 and ω = 100–300 μas at V = 20, where the ranges at a given magnitude depend upon the intrinsic brightness of the star (de Bruijne et al. 2005). Gaia Collaboration et al. (2018) recently published Gaia's Data Release 2 (DR2), an interim compendium of measurements based on 22 months of spacecraft operations.

Most small bodies in the solar system have ephemerides determined from photometric measurements that, in the absence of radar range measurements, have accuracy tied to the pixel scale, p, of the employed imaging detector.⁴ For individual observations, this is typically of order p = 1'' pixel⁻¹, corresponding to sky-plane displacements of ∼3700 km at the distance of the Trojan asteroids and ∼30,000 km for objects in the Kuiper Belt. When a minor planet occults a star, however, its instantaneous ephemeris is radically improved. An observed blocking of starlight indicates that a chord traversing the line of sight to the object also traversed the star, thereby establishing a time accuracy Δt ∼ δ/v_T, where δ is the angle subtended by the transit chord and v_T is the angular velocity of the body on the sky. The spatial accuracy on the plane of the sky is of order ω, the astrometric accuracy of the star's position. With the advent of Gaia DR2, the improvement factor from previous Hipparcos astrometric uncertainties is often measured in the thousands (Perryman et al. 1997; Lindegren et al. 2018). A Trojan asteroid observed to occult a bright Gaia star generates a one-time sky-plane positional accuracy of just a few tens of m, comparable to the order-of-magnitude tidal displacement dx ∼ 30 m induced by a Planet Nine–like object over a half-orbit timescale for the Trojan occulter. This startling potential accuracy is the primary motivator for the investigation presented here.

The population of Jovian Trojan asteroids is understood to be extensive, with estimates ranging from N ∼ 1 × 10⁴ to 2 × 10⁵ objects with D > 2 km in the combined L₄ and L₅ swarms (Jewitt et al. 2000; Yoshida & Nakamura 2005; Fernández et al. 2009; Wong & Brown 2015). Variations in N stem from assumptions regarding the underlying albedo distribution. The Gaia DR2, moreover, contains of order 3 × 10⁷ stars of apparent brightness V < 15 whose sky positions have been determined to an accuracy of ω ∼ 20 μas (de Bruijne et al. 2005). At a typical Trojan distance of a = 5.2 au, this is equivalent to a sky-plane positional accuracy of dx ∼ 75 m for the portion of the asteroidal body responsible for the occultation.

To estimate the rate of observable Trojan occultations, we must first assume a Trojan size distribution. We adopt the broken power-law differential magnitude distribution

$$\begin{eqnarray}\displaystyle \frac{{dN}}{{dH}}=\left\{\begin{array}{ll}{10}^{{\alpha }_{1}(H-{H}_{0})} & H\lt {H}_{b}\\ {10}^{{\alpha }_{2}H+({\alpha }_{1}-{\alpha }_{2}){H}_{b}-{\alpha }_{1}{H}_{0}} & H\geqslant {H}_{b}\end{array}\right.\end{eqnarray} \tag{ 2 }$$

and associated best-fitting parameters α₁ = 0.91, α₂ = 0.43, H₀ = 7.22, and H_b = 8.46, obtained in Wong & Brown (2015) by fitting the distribution of bright Trojans. The resultant magnitude distribution is shown in Figure 2. We combine this magnitude distribution with Equation (11) to obtain the corresponding size distribution, where we use the process described in Appendix A to set albedo values.

Zoom In Zoom Out Reset image size

We then sum the diameters of all inferred Trojans with D > 2 km (7.2 < H < 16.42) to obtain a total length of 1.4 × 10⁵ km for all 3.6 × 10⁴ large Trojans lined up side-by-side. Given P ∼ 12 yr orbital periods, each Trojan sweeps out ∼30° yr⁻¹, and the aggregate of Trojans sweeps out a projected area of ∼0.31 deg² yr⁻¹. For a uniform distribution of V < 15 stars on the sky, this corresponds to N ∼ 225 potentially observable occultations per year at each location on Earth. Using a network of N ∼ 2000 small telescopes, approximately 4.5 × 10⁵ events are thus observable annually, some of which may not be unique—i.e., multiple telescopes may observe a single occultation event. Assuming an f = 0.6 duty cycle to account for air mass, solar angle, and overlapping targets, as well as an additional f = 0.4 duty cycle when considering the day–night cycle, one plausibly expects ∼1.1 × 10⁵ unique occultation events per year on the network, of which a subset would be observable by multiple telescopes.

Focusing on the catalog of known Trojans in the JPL Small-Body Database, we provide a more detailed assessment to compare with the foregoing estimates. The occultation rate for an individual Trojan is

$$\begin{eqnarray}&&r={n}_{* }\sigma {v}_{{\rm{T}}},\end{eqnarray} \tag{ 3 }$$

where n_* is the areal number density of background stars, σ is the linear cross section for occultations, and v_T is the angular speed of the Trojan relative to the observer at the time of occultation. For the purpose of assessing rates, rather than planning specific observations, we adopt a simplified solar system model containing only the Earth, the Sun, and the Trojan of interest. We ignore Earth's axial precession and assume that background stars are point sources, which is an appropriate approximation at the scale of this problem. For example, a G2 solar-type dwarf star located at d = 1000 pc has V ∼ 15 and an angular size of δ ∼ 9.3 μas, smaller than the projected end-of-mission ω = 12–25 μas astrometric precision of Gaia.

To study the frequency of observable occultation events, we simulate the Sun–Earth–Trojan system for each of the 7207 Jovian Trojans listed in the JPL Small-Body Database on 2018 November 1. The starting spatial distribution of these Trojans relative to Jupiter and the Sun in our simulations is shown in Figure 3. Combining the REBOUND  orbital integrator (Rein & Liu 2012) with standard coordinate frame transformations, we obtain the R.A./decl. coordinates for each Trojan at 2000 equally spaced time steps over 24 yr—approximately the duration of two full Trojan orbits. We calculate the expected occultation rate for each Trojan at each time step (Figure 4), with each component of the rate equation obtained using the following procedure.

• 1.  
  n_*. We categorize the total number of stars along the ecliptic plane (to rough approximation, the orbital plane of the Jovian Trojans) into 360 R.A. bins that are evenly spaced such that each bin encompasses 1° in R.A. space. We choose to parameterize over R.A. space in order to encompass regions predominantly passing through the Galactic plane within individual bins.To accomplish this, we advance Jupiter's orbit at 360 evenly spaced time steps (corresponding to slightly varying step sizes in R.A. space) using REBOUND. Between each set of adjacent R.A. values, we query the Gaia DR2 database within ±25 from the midpoint decl. value to obtain a star count within each bin. Using these total star counts and the size of each bin on the sky, we calculate the stellar density in each bin. From this point, we linearly interpolate our results to obtain stellar densities in evenly spaced R.A. bins. We obtain n_* at each time step by selecting the R.A. bin within which the Trojan falls at that time, then retrieving the corresponding stellar density.
• 2.  
  σ. The effective occultation cross section for each Trojan can be estimated as the 2545 km latitudinal extent of the continental United States subtended at the instantaneous Earth–Trojan distance. This cross section takes into account the full range of telescopes extending across the country such that nearly all asteroids with D ≳ 2 km passing across the United States in the east–west direction can be seen in occultation by one or more network telescopes.Among the latitudinal distribution of telescopes, gaps larger than the diameter of the occulter may prevent the detection of an occultation, since small asteroids can slip through these gaps. Few large latitudinal gaps (33 gaps larger than 10 km in extent) exist within the cISP network design that we employ, and the full distribution of gaps is displayed in Figure 5. For the initial estimates presented here, we do not explore this in further detail but acknowledge that these gaps may lead to an underestimate of the total occultation rate. If desired, these gaps could be addressed by adding new, strategically placed telescope sites to the current network design, increasing the total number of sites from N = 1913 to ∼2000.
• 3.  
  v_T. This is the effective angular velocity of the Trojan, which is given by the absolute value of the sky-plane velocity of the Trojan relative to that of the Earth $(| {{\boldsymbol{v}}}_{{\rm{E}}}-{{\boldsymbol{v}}}_{{\rm{T}}}| )$ at each time step.

Zoom In Zoom Out Reset image size

Our resulting n_*, σ, and v_T values, as well as the corresponding occultation rates r, are displayed in Figure 4. All four maps order the Trojans by θ_jst, the angle between the Sun–Jupiter and Sun–Trojan vectors, where Figure 3 provides the orientation of the full system. We use the convention that negative angles correspond to Trojans lagging behind Jupiter in its orbit (the L5 swarm), while positive angles correspond to Trojans leading Jupiter (the L4 swarm).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The structure in the top panel of Figure 4 reflects the periodicity with which the two Trojan swarms pass through the Galactic midplane. Bright maxima in the background stellar density n_* of each swarm occur with an ∼12 yr periodicity corresponding to the frequency with which they pass through the Galactic midplane in the direction of the central Milky Way. Dimmer local maxima occur 6 yr apart from each of these brighter maxima, corresponding to passage through the outer region of the Galactic plane with a lower stellar density. The offset between central peaks in L4 and L5 is approximately 4 yr, in accordance with the spatial separation of the two swarms along Jupiter's orbit.

The second panel shows the linear cross section σ of the United States as viewed by each Trojan asteroid as it moves through its orbit. Alternating bright and dark streaks, each with a duration of ∼6 months, reflect the Earth's periodic movement toward and away from the two Trojan swarms during its orbit.

The third panel displays the sky-plane velocity v_T of each Trojan as a function of time, where dark streaks correspond to times during which the Trojan's motion is perpendicular to Earth's sky plane (approximately twice per Earth orbit). Between these dark streaks are alternating high- and medium-velocity local maxima. High-velocity local maxima correspond to times during which the Earth and the Trojan move in opposite directions within their orbits (for example, when they are on opposite sides of the Sun), increasing the Trojan's relative velocity as viewed from Earth. Conversely, medium-velocity local maxima occur where the Earth and the Trojan move in the same direction, such as when they are both on the same side of the Sun, reducing the Trojan's relative sky-plane velocity.

Finally, the fourth panel is produced by multiplying together the first three panels and incorporating an additional condition that all regions with θ_set < 38° have the occultation rate set to zero, where θ_set is the angle between the Earth–Sun and Earth–Trojan vectors. This final condition accounts for solar angle and air-mass considerations, where we require that, for a Trojan to be observable during the night, (1) the solar angle must be −8° or lower relative to the horizon, slightly below civil twilight (Nichols & Powers 1964), when the Trojan is in the sky, and (2) the asteroid must be at air mass z < 2, or placement in the sky at least 30° above the horizon. We assume that all Trojans are observable from the continental United States due to the orientation of their orbits roughly along the solar system ecliptic plane.

The overall structure of occultation rates in Figure 4 follows that of n_*, with overlaid dark streaks reflecting the influence of σ and v_T. Longer-lasting dark streaks are superimposed as well due to our requirement that r = 0 where θ_set < 38°.

From this distribution, we extract a median rate of 7.42 occultations per Trojan per year, with the rate distribution displayed in Figure 6. Extrapolating to the full Jovian Trojan population with D > 2 km (N ∼ 3.6 × 10⁴), this corresponds to ∼2.7 × 10⁵ total occultations per year. Some of these events will occur during the daytime and thus will not be observable. Approximating that, on average, 9/24 hr day^–1 can be spent observing, we estimate a total of 1.0 × 10⁵ observable occultation events per year, in excellent agreement with our previous estimate. This is a sufficiently high occultation rate to regularly track a large number of Trojan asteroids in search for the signatures of a massive perturber. We envision that, in conjunction with the proposed program to track Trojan asteroids, the network will be kept continually busy with observations of various small-body occultations throughout the solar system.

Zoom In Zoom Out Reset image size

Although this estimate does not account for weather conditions, we argue that, for all but the smallest targets, the extended nature of this network will mitigate weather problems. An occulting body will cross numerous network sites as it traverses the continent, making it unlikely for weather conditions to prevent observations by all telescopes along the path. In particular, the characteristic Rossby parameter ($\beta =2\omega \cos \phi /a$, where ϕ is the latitude, ω is Earth's angular rotation speed, and a is Earth's radius) at midlatitudes generally ensures nonuniform weather conditions across the country. For the smallest targets, which occult only one to two sites, however, weather conditions will play a larger role.

A set of three occultation measurements over a time baseline of several months or more permits, in idealized principle, the determination of all six osculating Keplerian orbital elements to the precision of the measured sky-plane positions (Bernstein & Khushalani 2000). To numerically evaluate the precision to which we can pin down the orbit of a Trojan with a given number of occultation measurements, we simulate a series of successive occultations and study their effect on refining the orbital elements of the occulter. We report results for medium-sized Jovian Trojan 3794 Sthenelos from the JPL Small-Body Database, with properties summarized in Table 2.

Table 2.  3794 Sthenelos

Parameter	Value	Uncertainty
a	5.204 au	1.246e−7 au
e	0.147	6.766e−8
i	0106	1163e−7
Ω	5990	9742e−7
ω	6215	1052e−6
M0	5355	5122e−7
D	34.531 km	0.364 km

Download table as:  ASCIITypeset image

We first assume that the fiducial reported orbit from the JPL Small-Body Database is the "true" orbit. We select N_occ points along the orbit at which we simulate occultations, and for each of these points, we obtain a separate coordinate transformation with $\hat{x}$ along the Earth–Trojan line of sight and $\hat{y}$ and $\hat{z}$ in the sky plane as measured from Earth. We set $\hat{z}$ in the direction of the Trojan's sky-plane velocity, while the perpendicular sky-plane direction is set as $\hat{y}$ (see Figure 7).

Zoom In Zoom Out Reset image size

The precision to which an asteroid's position can be determined within the sky plane is dependent upon two factors: (1) the astrometric uncertainty associated with the background star and (2) the network's observing cadence. Occultations do not provide a constraint in the $\hat{x}$ direction for any individual measurement; information from multiple events, however, can be combined to provide constraints in all three directions.

In the $\hat{y}$ direction, the positional uncertainty is equivalent to the occulted star's astrometric uncertainty. In our model, we set all astrometric uncertainties to 20 μas, corresponding to the projected end-of-mission Gaia astrometric performance for a star in the range V = 12–15 mag (de Bruijne et al. 2005; Gaia Collaboration et al. 2016). This is equivalent to a sky-plane positional uncertainty σ_y ∼ 75 m at a distance of 5.2 au; in our simulation, we set uncertainties at each occultation based on the fiducial Earth–Trojan distance at that time.

In the along-track direction within the sky plane, $\hat{z}$, the asteroid's positional uncertainty is limited by the observing cadence needed to obtain the requisite signal-to-noise ratio (S/N) for the occultation event. We aim for an observing cadence sufficiently high to reach the limits of the Gaia astrometric precision in the along-track direction.

Integrating over visible and near-infrared wavelengths from 400 to 800 nm, ∼4900 photons from a V = 15 star strike the aperture of a 16 inch telescope each second. To achieve 20 μas precision at d = 5.2 au, we must observe at 225 frames per second (fps), which results in 21 photons frame^–1 before accounting for inefficiencies in instrumentation. Achievement of this observing cadence is thus permitted by physics and will hinge on the available equipment. Frame rates for CCD and CMOS sensors are limited by the effective number of pixels used by the camera (e.g., El-Desouki et al. 2009). Because high resolution is not required for our purposes, in which each telescope observes a single star at a given time, we propose that the requisite precision can be attained by binning pixels to meet the appropriate frame rates needed for each star. As a result, we set our uncertainty in the $\hat{z}$ direction to 20 μas, as in the $\hat{y}$ direction.

To quantify the orbital element improvement with improved positional constraints, we employ a Markov chain Monte Carlo (MCMC) analysis using an affine-invariant ensemble sampler from the emcee Python package (Foreman-Mackey et al. 2013). We consider a set of five occultations measured in 1 yr intervals, directly reflecting the expected orbital improvement over a 5 yr baseline across which Planet Nine's perturbational signal should be observable. We note that we have adopted a conservative occultation measurement rate of 1 yr⁻¹, and in practice, it therefore should be possible to constrain the orbital elements to even higher precision with more frequent measurements. For comparison purposes, we also consider the case in which five occultation measurements are spread evenly across 12 yr, corresponding to ∼1 full Trojan orbit.

In each case, we use 100 walkers to sample six parameters (the six orbital elements), initialized with Gaussian priors where the fiducial orbital elements and associated uncertainties are set as the mean and standard deviation, respectively. At each step of the MCMC, we transform every proposed set of orbital elements into the coordinate frame of each occultation to find where that set of parameters places the asteroid at the measured occultation time. Then, models are accepted or rejected based on the log-likelihood function in Equation (4), where x_n and σ_n represent the fiducial asteroid position and associated uncertainty at the nth occultation, respectively, while μ_n represents the model position. We sum over the $\hat{y}$ and $\hat{z}$ directions for each occultation using the associated positional uncertainties σ_y and σ_z:

$$\begin{eqnarray}&&{ \mathcal L }(\mu | x,\sigma )=-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{n}\left(\displaystyle \frac{{\left({x}_{n}-{\mu }_{n}\right)}^{2}}{{\sigma }_{n}^{2}}+\mathrm{ln}(2\pi {\sigma }_{n}^{2})\right).\end{eqnarray} \tag{ 4 }$$

We run this process for 1000 iterations, and we find that the walkers stabilize after ∼200 iterations. We discard our first 800 iterations—well beyond the point of walker stability—as our "burn-in" and draw posteriors only from the final 200 iterations. Our results are displayed in Figure 8. The posterior distribution demonstrates that five occultation measurements across 5 yr can constrain each orbital element to a fractional uncertainty σ/μ ∼ 1.5 × 10⁻⁹ or better. With occultation measurements that are more spread out in time, the width of each posterior shrinks further, reflecting the exquisite precision to which orbits can be determined from occultation measurements over an extended temporal baseline.

Zoom In Zoom Out Reset image size

By incorporating a 2D model of the Kuiper Belt into the Ephemerides of the Planets and Moon 2017 version (EPM2017) planetary ephemerides, Pitjeva & Pitjev (2018) found that perturbations to the giant outer solar system planets due to the presence of the Kuiper Belt are comparable in magnitude to those induced by a 10 M_⊕ planet at a solar distance d = 550 au. Here we investigate how the two sources of influence can be disentangled by tracking their respective perturbations on Jupiter's Trojan asteroids.

We examine the change in each orbital element with time due to the presence of an external perturbing body. We complete this study in several steps. First, we compare analytical expressions for the force applied to the Trojan due to a ring and point source (Section 3.4.1). We then follow the formulation detailed in Murray & Dermott (1999; originally derived in Burns 1976), which derives analytical expressions for the derivative of each orbital element due to an external perturber (Section 3.4.2). Lastly, we simulate the effect of a point source versus ring perturber in REBOUND to assess net perturbations over an integrated orbit of the Trojan asteroid (Section 3.4.3). Throughout these analyses, we use the values listed in Table 1 and set Ω = ω = 0 and M = π to place Planet Nine at apoastron in its orbit.

3.4.1. Forces Imparted by an External Perturber

Evaluation of the instantaneous force imparted by each external perturber illuminates the comparative influences from Planet Nine and the Kuiper Belt, which we approximate as a 2D ring of uniform material. Planet Nine can be treated as a point mass (m_P9); thus, the acceleration, or force per unit mass imparted upon another point mass a distance d from Planet Nine, is a_P9 = −Gm_P9/d².

We consider two collinear configurations for the Trojan: one in which it is at the farthest point in its orbit from Planet Nine (Figure 9; top), and another in which it is at its closest approach to Planet Nine (Figure 9; bottom). The acceleration imparted on the Trojan is a_min = 2.01 × 10⁻¹¹ cm s⁻² in the former configuration and a_max = 2.06 × 10⁻¹¹ cm s⁻² in the latter.

Zoom In Zoom Out Reset image size

For comparison, we model the Kuiper Belt as a uniform circular ring with radius extending from R_min to R_max, where the angle between Planet Nine and any point in the Kuiper Belt is given by θ (see Figure 10). The acceleration imparted upon our test Trojan by the 2D Kuiper Belt is given by Equation (5), where we use the law of cosines to integrate the force per unit mass along the full angular and radial extent of the ring. Here σ_KB denotes the 2D mass density of the Kuiper Belt. The location of the Trojan is given as (r, 0) in polar heliocentric coordinates, where we define the axis representing θ = 0 along the line connecting the Sun and the Trojan.

$$\begin{eqnarray}&&{a}_{\mathrm{KB}}=-G{\sigma }_{\mathrm{KB}}{\int }_{{R}_{\min }}^{{R}_{\max }}{\int }_{0}^{2\pi }\displaystyle \frac{R(r-R\cos \theta )}{{\left({r}^{2}+{R}^{2}-2{rR}\cos \theta \right)}^{3/2}}d\theta {dR}\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

We place our test Trojan at r = 5.2 au and apply Equation (5) to study the instantaneous perturbation to the Trojan imparted by the Kuiper Belt. We obtain surface density σ_KB by dividing the total mass of the Kuiper Belt used in Pitjeva & Pitjev (2018), m = 1.97 × 10⁻² M_⊕, by the full surface area of our 2D disk. We integrate Equation (5) with limits R_min = 39.4 au and R_max = 47.8 au to emulate the disk model from Pitjeva & Pitjev (2018), resulting in a_KB = 1.10 × 10⁻¹⁰ cm s⁻². This is approximately an order of magnitude larger than the acceleration obtained from Planet Nine. Due to the symmetry of the Kuiper Belt, the magnitude of this acceleration is unaffected by the starting location of the Trojan, though the direction changes throughout the Trojan orbit. This result underscores the findings of Pitjeva & Pitjev (2018) and motivates further study to disentangle the perturbational effect of Planet Nine from that of the Kuiper Belt.

3.4.2. Analytical Perturbations to Orbital Elements

To explore this problem in greater detail, we also consider the effect of each structure on the orbital element evolution of our fiducial Trojan (Roy 2005). We consider a force $F=\bar{R}\hat{r}\,+\bar{T}\hat{\theta }+\bar{N}\hat{z}$ imparted by a perturber with mass M_p, where $\bar{R}$, $\bar{T}$, and $\bar{N}$ are components of force in the $\hat{r}$, $\hat{\theta }$, and $\hat{z}$ cylindrical coordinate directions, respectively. We set all mutual inclinations to zero and thus have $\bar{N}=0$. Zero force applied in the normal direction necessitates that $\dot{I}$ and $\dot{{\rm{\Omega }}}$ must be equivalently zero. Expressions for all other orbital element derivatives, following the formulation of Murray & Dermott (1999), are provided in Equations (6)–(10). We report changes in τ, the time of perihelion passage, rather than M₀, where the two are related by M₀ = −nτ.

$$\begin{eqnarray}&&\mu ={{GM}}_{p}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\dot{a}=2\displaystyle \frac{{a}^{3/2}}{\sqrt{\mu (1-{e}^{2})}}[\bar{R}e\sin f+\bar{T}(1+e\cos f)]\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\dot{e}=\sqrt{\displaystyle \frac{a(1-{e}^{2})}{\mu }}[\bar{R}\sin f+\bar{T}(\cos f+\cos E)]\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\dot{\omega }=\sqrt{\displaystyle \frac{a(1-{e}^{2})}{\mu {e}^{2}}}\left[-\bar{R}\cos f+\bar{T}\sin f\displaystyle \frac{2+e\cos f}{1+e\cos f}\right]\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\begin{array}{rcl}\dot{\tau } & = & \left[3(\tau -t)\sqrt{\displaystyle \frac{a}{\mu (1-{e}^{2})}}e\sin f\right.\\ & & +\left.{a}^{2}(1-{e}^{2}){\mu }^{-1}\left(\displaystyle \frac{-\cos f}{e}+\displaystyle \frac{2}{1+e\cos f}\right)\right]\bar{R}\\ & & +\left[3(\tau -t)\sqrt{\displaystyle \frac{a}{\mu (1-{e}^{2})}}(1+e\cos f)\right.\\ & & +\left.{a}^{2}{\mu }^{-1}(1-{e}^{2})\left(\displaystyle \frac{\sin f(2+e\cos f)}{e(1+e\cos f)}\right)\right]\bar{T}.\end{array}\end{eqnarray} \tag{ 10 }$$

We first study the effect of a single point-mass perturber to characterize the gravitational impact of Planet Nine. We neglect the tangential and normal components of the force ($\bar{T}=\bar{N}=0$) and focus solely on Planet Nine's radial force. For our purposes, this simplified model is sufficient to obtain an estimate of the expected change in orbital elements over time, since Planet Nine's large predicted separation ensures that its radial force is substantially stronger than the forces in the other two reference directions. We calculate $\bar{R}$ using Planet Nine's apoastron distance from the Sun. We report results for Trojan 3794 Sthenelos, again using the orbital elements from Table 2. The resultant orbital element derivatives are provided in Table 3.

Table 3.  Instantaneous Change in Orbital Elements Induced by Planet Nine and the Kuiper Belt

	Planet Nine	Kuiper Belt
$\dot{a}$	−8.50e−10	1.69e−7
$\dot{e}$	−5.44e−10	1.65e−9
$\dot{\omega }$	−4.77e−9	−3.82e−7
$\dot{\tau }$	−3.47e−9	−1.06e−6

Note. All values are given in units of the orbital element per year, where a is in astronomical units and ω is in radians.

Download table as:  ASCIITypeset image

We also complete this exercise for a circular 2D ring in the Trojan's orbital plane to obtain the instantaneous change in orbital elements induced by the Kuiper Belt. This model includes a tangential force but again has zero normal force. Following the procedure of Pitjeva & Pitjev (2018), our model of the Kuiper Belt includes three rings, where the inner and outer ring are located at R_min = 39.4 au and R_max = 47.8 au, respectively, corresponding to the 3:2 and 2:1 orbital resonances with Neptune. The inner and outer ring each contain 40 discrete, evenly spaced point masses. The central ring contains 80 point masses and is located at r = 44.0 au. The total mass of the Kuiper Belt, set as 1.97 × 10⁻² M_⊕ (Pitjeva & Pitjev 2018), is evenly distributed among these 160 points comprising the Kuiper Belt. The model described here is visualized in Figure 11, and the instantaneous changes in orbital elements resulting from the Kuiper Belt are listed in Table 3.

Zoom In Zoom Out Reset image size

The magnitude of the instantaneous change in each orbital element is significantly smaller in the case of Planet Nine than in the case of the Kuiper Belt. However, a continuous force in a single direction—from the effectively static Planet Nine at apoastron—should have a different effect over time from a ring imparting a symmetric force throughout the Trojan's orbit, where the effect would in part cancel out. We explore this further in Section 3.4.3.

3.4.3. Numerical Integration with REBOUND

We next compound the perturbative effects of both Planet Nine and the Kuiper Belt on a single Trojan using the REBOUND orbital integrator. We consider two test systems and compare the perturbations upon the Trojan in each: (1) a system including only the Sun, the test Trojan, and Planet Nine, and (2) a system including the Sun, the test Trojan, and the same model of the Kuiper Belt described in Section 3.4.2 (Figure 11).

All components of each model are treated as point masses, and we integrate over one full orbit of the test Trojan, tracking perturbations from the initial Keplerian orbit over time. For simplicity, we set all inclinations in both models to zero such that all bodies lie within the same orbital plane. We note that our analogous tests that incorporate the inclination of Planet Nine (∼20°) behave qualitatively in the same way within the plane of perturbations (Δr × Δϕ) that we consider here.

For consistency, we again provide results for Trojan 3794 Sthenelos, with properties listed in Table 2, with the exception that we set the inclination i = 0°. We also run the same REBOUND simulations for several other Trojans of various sizes and confirm that there are no qualitative differences in the results.

The results of our integrations are shown in Figures 12 and 13. We use the convention Δϕ = ϕ_p − ϕ_T and Δr = r_p − r_T, where the p subscript refers to the system with a perturber, while the T subscript refers to the system including only the Sun and the Trojan, with no external perturber. Thus, a positive value for Δr, as is the case for the Kuiper Belt where e_T = 0, indicates that the perturber has pushed the Trojan outward in its orbit.

Zoom In Zoom Out Reset image size

In Figure 12, we display the perturbations induced by Planet Nine for a range of initial angles between the Sun–Planet Nine and Sun–Trojan vectors (θ_pst). The perturbation induced by the Kuiper Belt is not clearly visible on this scale and is thus omitted. At the beginning of each orbit, Planet Nine pushes the Trojan asteroid radially outward (positive Δr) for θ_pst < π/2 and radially inward (negative Δr) for θ_pst > π/2. Here Δϕ is set by the location of the Trojan in its orbit: positive as the Trojan moves toward Planet Nine and negative as it moves away. Thus, all four quadrants of this space are populated by perturbations, with the perturbative direction set by the Planet Nine–Sun–Trojan orientation.

However, the differential perturbation induced by Planet Nine—the measurable perturbation from within the solar system—is substantially smaller than the global perturbation shown in Figure 12. This is because Figure 12 does not incorporate the acceleration also imparted upon the Sun by Planet Nine.

In Figure 13, we subtract the Sun's motion from the system to show the displacement of the Trojan relative to the Sun. To elucidate the behavior of the Trojan in multiple scenarios, we display the perturbations produced in two cases: (1) where the Trojan eccentricity e_T is set to zero (left and middle panels in Figure 13) and (2) where we use the fiducial, nonzero Trojan eccentricity extracted from the JPL Small-Body Database (right panel in Figure 13).

Zoom In Zoom Out Reset image size

For nonzero Trojan eccentricities, the Kuiper Belt can produce perturbations in the +Δr, −Δr, and −Δϕ directions but never in the +Δϕ region of parameter space. In contrast, Planet Nine can produce perturbations in all four quadrants. We obtain similar results when, instead of using the Kuiper Belt model in Figure 10, we model all but the three most massive TNOs listed in the JPL Small-Body Database (Pluto, Eris, and Makemake). Despite the lack of perfect symmetry in this case, we still find only small (∼5 μas) deviations into the +Δϕ space, suggesting that, in the absence of additional undiscovered Pluto-sized objects, the symmetry of the Kuiper Belt is sufficient such that Planet Nine should be gravitationally distinguishable from the Kuiper Belt. Given the presence of an additional Pluto-sized perturber in the Kuiper Belt, this network would provide strong evidence for its existence and constraints upon its location—an interesting result in itself.

It is possible to distinguish the effects of Planet Nine and the Kuiper Belt, therefore, by studying occultations of the Trojans in regions of the sky that correspond to positive Δϕ for a Planet Nine perturber but still result in a negative Δϕ from the Kuiper Belt. These include regions in which the Trojan is moving toward Planet Nine in its orbit and thus is pulled further forward by the perturber, as illustrated in Figure 14. The perturbational signature of the putative Planet Nine will deviate most from the noise produced by the Kuiper Belt during the half of the Trojan orbit in which perturbations are in the +Δϕ direction. The peak of this +Δϕ deviation occurs when θ_pst = π/2 and the Trojan is moving toward Planet Nine. As a result, the angular orientation of Planet Nine can be extracted upon detection of the predicted perturbation from a distribution of Trojans.

Zoom In Zoom Out Reset image size

The perturbation size detectable by our network is limited by the angular precision of the occulted star (ω = 12–25 μas at V = 15 from the Gaia catalog). This is much more precise than standard astrometric measurements, and even astrometry with HST typically results in uncertainties of order ∼1 mas or greater (e.g., Bellini et al. 2011; Porter et al. 2018) and can only reach a precision of ∼0.1 mas for the very brightest targets (Benedict et al. 2002). Thus, the extremely high precision necessary to observe Planet Nine's perturbational effect can be decisively obtained through the occultation method.

By measuring a large number of Trojan positions over time and searching for this +Δϕ signature of a point perturber, we can extract the Planet Nine signal and separate it from that of the Kuiper Belt. For an individual small Trojan, we have shown that ephemeris refinement from successive occultations through the course of a full orbit generates a final positional accuracy of order $\delta {\boldsymbol{x}}=75\,{\rm{m}}.$ Figure 13, on the other hand, suggests that systematic deflections on the Trojan orbit due to Planet Nine are of order 25 m, corresponding to a detection of σ_i = 1/3 σ significance. As a consequence, the tracking of a single Trojan is no more useful for detecting Planet Nine than was the Cassini telemetry. With a large occultation network, however, thousands of bodies can be monitored with a detection significance rising to ${\sigma }_{\mathrm{final}}={\sigma }_{i}\sqrt{N}$, where N is the number of Trojans whose orbits have been improved to the DR2-enabled precision. As an example, refinement of N ∼ 225 Trojan orbits would generate a σ_final = 5σ detection significance.

4.1.1. The Lucy Mission

4.1.2. The Psyche Mission

Trojan asteroids, introduced in the Jovian context in Section 3.1, are more generally defined as small bodies librating at the L4 and L5 Lagrange points of a planet–star system in a 1:1 mean-motion resonance, first theorized by Lagrange (1873) as a stable solution to the restricted three-body problem. Wolf (1906) discovered the first known Trojan asteroid at Jupiter's L4 point, and the vast majority of Trojan asteroids discovered since then are also located in the Jovian system.

More recently, a population of Trojans was also found in the Neptunian system (Chiang et al. 2003; Sheppard & Trujillo 2006, 2010a, 2010b; Becker et al. 2008), with a comparable number of asteroids at each Lagrange point and similar dynamics in both populations. Previous works have explored possible capture mechanisms for the Jovian and Neptunian Trojans (e.g., Marzari et al. 2002; Morbidelli et al. 2005; Sheppard & Trujillo 2006; Lykawka et al. 2009), where most of the recent models incorporate capture from the same population that now comprises the modern-day Kuiper Belt, as in the Nice model (Tsiganis et al. 2005). However, the exact capture mechanism of these two populations remains to be definitively established.

Sheppard & Trujillo (2006) and later Jewitt (2018) observed that the color distribution of the Jovian Trojans is statistically indistinguishable from that of the Neptunian Trojans, consistent with a shared origin for both populations. However, Jewitt (2018) found that this distribution does not match that of any dynamical subpopulation within the Kuiper Belt, where these subpopulations are defined in Peixinho et al. (2015). This discrepancy refutes the hypothesis that the two Trojan populations were captured as scattered KBOs, raising the question of where else the Trojans may have instead originated.

A robust comparison of asteroid size distributions would serve as another useful distinguishing factor, since the parent population should have a similar size distribution to that of the captured Trojans. The H-magnitude distribution of the large Jovian Trojans, which is closely tied to the radius distribution, has been measured on several occasions and matches a steep power law with index α (α = 0.9 ± 0.2, Jewitt et al. 2000; α = 1.0 ± 0.2, Fraser et al. 2014; α = 0.91 ± ^+0.19_−0.16, Wong & Brown 2015); however, the paucity of bright Trojans limits the accuracy to which α can be determined. The translation of an H-magnitude distribution to a radius distribution also requires assumptions about the albedos of the Trojans, which could be circumvented through direct size observations using stellar occultations. Currently, the ranges of asteroid sizes measured for the Trojan and Kuiper Belt populations barely overlap—thus, precise constraints on both α and the size distribution of small KBOs would enable a more convincing direct comparison.

Another appeal of this network lies within its ability to pinpoint the positions of distant minor planets, such as TNOs, to high precision. While challenging, stellar occultation measurements of TNOs have been successfully demonstrated on several occasions (Ortiz et al. 2017; Müller et al. 2018; Sickafoose et al. 2019) and have become increasingly accessible with the advent of the Gaia DR2. Deep sky surveys such as the Dark Energy Survey (Dark Energy Survey Collaboration; Fermilab & Flaugher 2005) can also provide astrometric measurements of TNOs with sufficient accuracy to reliably predict upcoming occultation events (Banda-Huarca et al. 2019). By observing TNO occultations over time, it is possible to constrain the allowed parameter space for undiscovered minor planets in the outer solar system, leading to an effective "acceleration map." As new extreme TNOs with a > 250 au are discovered, such a map will contextualize the distribution of known objects and constrain their abundance.

Furthermore, occultation measurements provide the added benefit of determining TNO shapes and sizes. Stellar occultation measurements are currently the only ground-based observational technique capable of constraining the geometry of TNOs to kilometer accuracy. However, as described in Appendix B, the Fresnel diffraction limit becomes important for small TNOs and must be taken into consideration when planning TNO occultation measurements.

Occultation measurements provide an independent test to verify previously obtained asteroid diameters, as well as tighter constraints to decrease the uncertainties in reported values. A combination of measurements from the WISE (Wright et al. 2010) and NEOWISE (Mainzer et al. 2011) infrared surveys has been modeled using the Near-Earth Asteroid Thermal Model (NEATM; Harris 1998), resulting in a database of measured diameters for ∼150,000 asteroids (Mainzer et al. 2016). The methods and underlying assumptions used by the NEOWISE team to obtain these asteroid diameters were recently challenged (Myhrvold 2018a, 2018b, 2018c), with a corresponding rebuttal in Wright et al. (2018). Independent measurements obtained through separate methods can improve confidence in reported values obtained through thermal modeling, and this could serve as an external test of the accuracy of previously measured values.