# Can the gravitational effect of Planet X be detected in current-era tracking
of the known major and minor planets?

Daniel C. H. Gomes<EMAIL_ADDRESS>Department of Physics & Astronomy,
University of Pennsylvania, 209 S. 33rd St., Philadelphia, PA 19104  Zachary
Murray Center for Astrophysics — Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138  Rafael C. H. Gomes Department of Physics & Astronomy,
University of Pennsylvania, 209 S. 33rd St., Philadelphia, PA 19104  Matthew
J. Holman Center for Astrophysics — Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138  Gary M. Bernstein Department of Physics & Astronomy,
University of Pennsylvania, 209 S. 33rd St., Philadelphia, PA 19104

###### Abstract

Using Fisher information matrices, we forecast the uncertainties $\sigma_{M}$
on the measurement of a “Planet X” at heliocentric distance $d_{X}$ via its
tidal gravitational field’s action on the known planets. Using planetary
measurements currently in hand, including ranging from the Juno, Cassini, and
Mars-orbiting spacecraft, we forecast a median uncertainty (over all sky
positions) of $\sigma_{M}=0.22M_{\oplus}(d_{x}/400\,\textrm{AU})^{3}.$ A
$5\sigma$ detection of a $5M_{\oplus}$ Planet X at $d_{X}=400$ AU should be
possible over the full sky but over only 5% of the sky at $d_{X}=800$ AU. The
gravity of an undiscovered Earth- or Mars-mass object should be detectable
over 90% of the sky to a distance of 260 or 120 AU, respectively. Upcoming
Mars ranging improves these limits only slightly. We also investigate the
power of high-precision astrometry of $\approx 8000$ Jovian Trojans over the
2023–2035 period from the upcoming Legacy Survey of Space and Time (LSST). We
find that the dominant systematic errors in optical Trojan astrometry
(photocenter motion, non-gravitational forces, and differential chromatic
refraction) can be solved internally with minimal loss of information. The
Trojan data allow cross-checks with Juno/Cassini/Mars ranging, but do not
significantly improve the best-achievable $\sigma_{M}$ values until they are
$\gtrsim 10\times$ more accurate than expected from LSST. The ultimate
limiting factor in searches for a Planet X tidal field is confusion with the
tidal field created by the fluctuating quadrupole moment of the Kuiper Belt as
its members orbit. This background will not, however, become the dominant
source of uncertainty until the data get substantially better than they are
today.

Planetary dynamics, Fisher’s Information, Kuiper Belt, Trans-neptunian
objects, Astrometry, Gravitation

##  1 Introduction

State-of-the-art ranging and astrometry of solar system bodies has yielded
several of the greatest discoveries in the history of astronomy and physics,
including the development of laws of planetary motion (Kepler, 1609, 1618) and
general motion (Newton, 1687), the successful prediction of Neptune’s
existence from anomalies in Uranus’s orbit (Le Verrier, 1846), the precise
measurement of the astronomical unit (Muhleman et al., 1962; Smith et al.,
1966; Ash et al., 1966), and the validation of General Relativity (Einstein,
1915; Shapiro et al., 1968). Spectacular mis-interpretations of the data
include the prediction of Planet Vulcan (Le Verrier, 1859). We use a Fisher
matrix formulation to estimate the sensitivity of historical
positional/astrometric measurements of the major planets, plus data expected
in the next decade, to discover the gravitational signature of unknown bodies
in the solar system, which we will generically refer to as “Planet X,” an
appellation first used by Dallet (1901) for a hypothesized trans-Neptunian
planet. There has in particular been attention in recent years to a potential
Neptune-mass Planet X lurking hundreds of AU from the Sun, sometimes referred
to as “Planet 9.” But there is also some theoretical motivation for, and weak
observational limits on, Earth- or Mars-mass bodies at closer distances (e.g.
Volk & Malhotra, 2017; Silsbee & Tremaine, 2018), and the gravitational
constraint forecasts in this work are equally applicable to these as Planet X.

###  1.1 Searching for Planet X

Although the prediction and discovery of Neptune motivated searches for
additional planets (see Tremaine 1990 and Batygin et al. 2019 for reviews),
our census of the solar system remained relatively static for many decades.
Tombaugh’s 1930 discovery of Pluto (Tombaugh, 1946, 1996) marked the beginning
of a continuing, almost fantastical era of outer solar system discovery,
including Oort’s hypothesis of a nearly spherical shell of comets orbiting our
sun at great distances (Oort, 1950); the discovery of the first Centaurs
(Kowal et al., 1979) and the realization that there must be additional
reservoirs supplying such a dynamically short-lived population (Duncan et al.,
1988; Holman & Wisdom, 1993); the discovery of 1992 QB1, the first trans-
neptunian object (TNO) by Jewitt and Luu (Jewitt & Luu, 1993); and the
recognition of the intricate dynamical structure of the small body population
beyond Neptune and what it implies for the formation of the solar system (see
Gladman & Volk 2021 for a review). Of course, this era includes the
unbelievably successful exploration of our solar system by spacecraft
encountering, orbiting, landing on, and crashing into members of our solar
system, large and small.

All the while, the hypothesis of additional massive planets, either at an
earlier era or still residing in the outer solar system, continues to be
compelling. Supporting evidence primarily comprises features in the orbital
structure in the outer solar system that can be explained by the long-term
dynamical influence of an additional planet (Gladman et al., 2002; Brunini &
Melita, 2002; Melita et al., 2004; Gladman & Chan, 2006; Lykawka & Mukai,
2008; Huang et al., 2022). Most recently, the orbital distribution of trans-
neptunian objects (TNOs) with highly eccentric orbits and large perihelion
distances has motivated new searches (Trujillo & Sheppard, 2014; Sheppard &
Trujillo, 2016; Batygin & Brown, 2016; Brown & Batygin, 2016). Attempts to
discover Planet X have followed two primary approaches: searching for the
light reflected or emitted by the planet itself; and searching for the
signatures of the short-term gravitational influence of Planet X on other
bodies in the solar system.

Discovering Planet X directly requires large portions of the sky to be
carefully surveyed to faint magnitudes, unless dynamical or other evidence can
narrow the search region (Brown & Batygin, 2016; Malhotra et al., 2016). Only
a few telescopes have the combination of aperture and field of view required
to support an efficient search over a large area of sky. Furthermore, the
search algorithms and survey cadence used must be sensitive to the very slow
rates of apparent motion for distant solar system bodies observed from Earth.
Nearly every platform with which a meaningful search could be conducted has
been used or considered, including Spacewatch (Larsen et al., 2007), the
Catalina Sky Survey (Brown et al., 2015), Palomar (Brown & Batygin, 2022),
Pan-STARRS (Holman et al., 2018), the Canada-France-Hawaii telescope
(Bannister et al., 2018), DECam (Bernardinelli et al., 2022a, b), NASA’s IRAS
(Beichman, 1988), WISE (Luhman, 2014) and TESS missions (Holman et al., 2019;
Rice & Laughlin, 2020), and even CMB experiments (Baxter et al., 2018; Naess
et al., 2021).

Many claims of the existence of additional unseen planets, based on trends in
the residuals of the observations against the available ephemerides, have been
made over the years. However, Hogg et al. (1991) demonstrated in a series of
numerical experiments that such trends are to be expected when an orbital
solution is fit to observations spanning less than a full period and then
extrapolated to the times of new observations without incorporating the new
data. The trends vanish when the orbit is re-fit using the full set of
observations. Indeed, Standish (1993) dispatched all claims existing at the
time by showing that the trends seen at the time disappear when improved
masses of the giant planets (based on Voyager data) are used and the orbits
re-fit. Hogg et al. (1991) presciently concluded in their investigation that
continued searches for Planet X would not be well motivated until improved
data on the outer planets, from Cassini or Galileo or pulsar timing, became
available.

In an inspiring investigation, Fienga et al. (2016) used just such data, a
long span of Cassini ranging observations, to constrain the mass and location
of Planet X. The orbit for Planet X in Batygin & Brown (2016) included all the
elements except for the true anomaly, or its equivalent, because the dynamical
influence on the extreme TNOs was argued to be primarily secular. By inserting
a synthetic Planet X with an assumed mass of $10~{}M_{\oplus}$, carrying out
full ephemeris solutions to all the available data, and comparing the
residuals to those without the addition of Planet X, Fienga et al. (2016) were
able to show that a broad range of initial longitudes is excluded because the
residuals would be excessive. Furthermore, the residuals would be
significantly improved by the addition of a large planet at other initial
longitudes. Holman & Payne (2016a) extended this work to include a wide range
of masses and initial positions of the planet across the entire sky. The
Fienga et al. (2016) investigation is based on an earlier version of the INPOP
ephemeris model that did not include mass in the Kuiper belt aside from Pluto
(Fienga et al., 2016), and the Holman & Payne (2016a) study is built upon the
residuals shown in Fienga et al. (2016) and therefore inherits most of its
dynamical properties. As demonstrated by Pitjeva & Pitjev (2018), however,
incorporating the mass of the Kuiper belt, through a combination of objects
with known masses and orbits and massive rings to account the gravitational
influence of a distribution of smaller bodies, reduces the residuals in the
ephemeris fits. In particular, it reduces the residuals in the Cassini ranging
observations. Fienga et al. (2020) reprised the investigation using the
INPOP19a model, which includes the gravitational influence of the Kuiper belt,
and a tidal model as in Holman & Payne (2016a). Their results exclude an
additional planet within $\sim 25\arcdeg$ from the ecliptic plane, closer than
$500~{}AU$ if it is $5~{}M_{\oplus}$ or closer than $650~{}AU$ if it is
$10~{}M_{\oplus}.$

Following the approach of Hogg et al. (1991), Holman & Payne (2016b) tested
whether fits to the orbits of Pluto and a set of well-observed TNOs would be
improved by the inclusion of a Planet X. They focused on the data set
developed and used by Buie & Folkner (2015) to refine the orbit of the Pluto-
Charon in preparation for the New Horizons encounter. The results of Holman &
Payne (2016b) suggest a planet that is either more massive or closer than the
best-first results of Batygin & Brown (2016) or Fienga et al. (2016). This
result was driven by a clear trend in the residuals over the course of two
decades which may be due to systematic errors. Alternatively, it could also be
explained by perturbations from a different body, aside from Planet X, that is
closer to Pluto but less massive than $5-10~{}M_{\oplus}$.

In this work, we carefully examine the power of extant and imminent
observations (astrometry, ranging, and masses derived from satellites) to
constrain the mass of Planet X, while allowing for freedom in the masses and
orbits of all known significant gravitating bodies in the solar system. We do
this as a function of all possible Planet X locations in the sky.

###  1.2 A new era for ground-based astrometry

In the past 50 years, active radio ranging to spacecraft orbiting or flying by
planets and radar ranging to the inner planets have largely superseded ground-
based optical astrometry for high-precision ephemeris constraints. Radio
ranging to the terrestrial planets routinely reaches sub-meter accuracy, and
lunar laser ranging attains mm-level accuracy. In the search for a
gravitational signature for Planet X, however, the tidal acceleration $a_{\rm
tidal}$ increases linearly with heliocentric distance $d_{T}$ of the test
body, and the time over which the perturbations oscillate is the period $P,$
leading some signals to scale as $a_{\rm tidal}P^{2}\propto d_{T}^{4}$ (see
Section 2.1.3). Hence, there is a signal-to-noise ($S/N$) premium on
observations of bodies in the outer solar system, as well as the advantage of
reducing the influence of potentially confounding signals from uncertainties
in the mass distribution of the asteroid belt. But spacecraft ranging to the
outer planets is much sparser, with meter-level precision to Jupiter (Juno)
and Saturn (Cassini) and only single Voyager flybys to Uranus and Neptune. Our
first focus will be to assess the power of past and ongoing ranging
observations (primarily of Mars, Jupiter, and Saturn) to constrain Planet X.

Our second focus is to assess the power of traditional ground-based optical
astrometry to augment state-of-the-art ranging data in constraining Planet X.
Ground-based astrometry is currently undergoing a revolution in both _quality_
and _quantity_ of available data. The per-exposure astrometric accuracy of
routine ground-based imaging has been improved to milliarcsecond (mas) level
by three developments:

  * •

Bernstein et al. (2017) demonstrate $\approx 2$ mas accuracy in mapping the
500-Mpix focal plane of the Dark Energy Camera to (relative) sky coordinates.
This is done primarily by solving for internal consistency of stellar
positions on multiple shifted exposures of the same star field.

  * •

The availability of a dense global network of stars with $\approx 20$
$\mu\textrm{as}$ absolute astrometric accuracy from the Gaia spacecraft allows
one to tie any ground-based exposure to the absolute reference frame with very
high accuracy.

  * •

For bright point sources ($m_{r}\lesssim 20,$ depending on the observatory),
the dominant astrometric error from the ground is then stochastic refraction
from the turbulent atmosphere. Fortino et al. (2021a) demonstrate that
turbulence components on scales $\gtrsim 1\arcmin$ can be mapped and removed
from each exposure by fitting a Gaussian Process model to the deviations of
Gaia stars from their published positions. This reduces the RMS turbulence
signal to $\approx 2$ mas per axis.

The revolution in data _quantity_ is best exemplified by the upcoming 10-year
Legacy Survey of Space and Time (LSST) to be conducted from the 8-meter
telescope at the Vera Rubin Observatory. There are good reasons to expect the
performance of the LSST camera to be better than that of the Dark Energy
Survey. We can therefore plausibly expect the floor on astrometric errors of
LSST observations of unresolved minor planets to be $\sigma_{\rm min}\approx
2$ mas. This is better than the Gaia single-epoch uncertainty for sources with
$G>18$ (Tanga et al., 2022).

But LSST will observe essentially _every_ minor planet at declination
$<30\arcdeg$ nearly 1000 times in its 10-year survey, including for example
$\sim 8000$ Jovian Trojan asteroids at $H<14.5$. As detailed in Section 6, the
collective accuracy on a shift of this entire population is $\approx
2.5\mbox{$\mu\textrm{as}$}$, which is 8 meters at the distance of the Trojans.
This is comparable to the accuracy obtained by the Cassini and Juno
spacecraft. LSST’s collective precision on the main-belt asteroid (MBA)
population will be at the sub-meter level, given that they are $2\times$
closer, and $O(10^{5})$ of them will be bright enough to approach $\sigma_{\rm
min}.$

The tracking of many minor planets can have advantages over the Cassini and
Juno ranging. First, ranging is only sensitive to displacements in the radial
direction (namely the invariable plane), whereas astrometry measures 2
transverse displacements. Second, the minor planets have larger eccentricities
and a wide range of inclinations, which enhances the types of perturbations
(e.g. a polar Planet X) that are detectable. Third, the Trojans and MBAs are
distributed in space and we can track each for $\geq 1$ full orbit, which
substantially reduces degeneracies between potential sources of gravity.
Cassini data are available for only 45% of Saturn’s orbital period.

Note that ranging to Mars has much of these latter two advantages while also
having accuracy at cm levels—so it is interesting to ask which combinations of
Mars ranging, Jupiter$+$Saturn ranging, and LSST Trojan tracking yield the
strongest constraints.

This work is inspired by the proposal by Rice & Laughlin (2019) to measure the
Planet X gravitational perturbation via a USA-spanning network of several
thousand small telescopes monitoring occultations of Gaia stars by Jovian
Trojans. The occultation method offers higher precision and reduced systematic
errors relative to traditional astrometry. It is of great interest, however,
to ask what can be achieved—at zero marginal cost in instrumentation—from the
observations to be obtained by LSST, and to investigate more rigorously the
potential confusion between Planet X and other solar system gravitational
anomalies.

Ground-based astrometry of minor planets has some systematic errors that do
not exist for spacecraft tracking, e.g. differential chromatic refraction
(DCR) in the atmosphere, photocenter motion (PCM), and non-gravitational
forces. We will investigate how much these systematic errors degrade the
potential statistical accuracy of the LSST Trojan tracking.

###  1.3 Plan of this paper

We begin in Section 2 by reviewing some analytic expressions for the orbital
perturbations induced by the gravity of a Planet X, to understand the size and
scaling of the signals of interest. In Section 3 we describe how the Fisher
matrix is used to forecast the expected uncertainties $\sigma_{M}$ in the mass
$M_{X}$ of a putative Planet X. We describe the model of our solar system used
to derive these limits, including some unresolvable degeneracies between
Planet X and the tidal fields generated by Kuiper Belt members. In Section 4
we present the numerical methods used to acquire the derivatives of
observations with respect to parameters required for the Fisher method. In
Section 5 we present the results of the Fisher analysis in terms of expected
accuracy $\sigma_{M}$ on the mass of Planet X vs its location in the sky,
using completed and planned observations of the major planets. We quantify how
the addition of LSST observations of Trojans affect the constraints on Planet
X mass, including a formulation of the major systematic errors, in Section 6.
In Section 7 we apply our Fisher matrix results to scenarios including LSST
Trojans, quantifying the accuracy on Planet X and our ability to distinguish
systematic errors in the data from the signature of Planet X. We conclude in
Section 8. The appendices give the details of some aspects of the calculation.

##  2 The back of the envelope

Before solving numerically for the sensitivity $\sigma_{M}$ to the mass of
Planet X, it is useful to have a rough estimate of the sizes of displacements
that the presence of Planet X can cause on decade time scales. We will assume
that we are tracking some test body at heliocentric distance $d_{T}$ with
radius $R_{T}$, and that this is being perturbed by an effectively stationary
(over 50 years) point mass $M_{X}$ at distance $d_{X}.$ Brown & Batygin (2021)
use the dynamics of extreme TNOs to suggest $M_{X}\approx 6M_{\oplus}$ in an
orbit ranging from roughly 300–500 AU heliocentric distance. For nominal
values we will take $d_{T}=5.2$ AU, $M_{X}=5M_{\oplus},$ and $d_{X}=400$ AU,
corresponding to the effect of a Planet X in the center of their range as
influencing Jupiter or its Trojans. The choice of nominal value is simply to
scale the constraints into a useful range, not a statement on the likelihood
of or constraints on a Planet X. We describe Planet X’s position by its
ecliptic azimuthal angle $\phi_{X}$ and polar angle $\theta_{X}$. It will also
be useful to define the radial, azimuthal, and vertical directions relative to
the ecliptic as $\hat{\rho}$, $\hat{\phi}$, and $\hat{z}$, respectively. Brown
& Batygin (2021) derive favored regions of the sky, but we remain agnostic
about Planet X’s location and derive constraints for all sky locations, so
that our forecasts remain generally applicable to any Planet X hypothesis.

###  2.1 The signatures of Planet X

We can decompose the gravitational field of Planet X into spherical harmonics
about the solar system barycenter with radial dependences $r^{\ell}$. The
$\ell=0$ and $\ell=1$ terms are a potential shift and a constant acceleration,
respectively, which are unobservable with differential measurements within the
solar system. The $\ell=2$ tidal terms are a sufficient description of the
observable effects of Planet X, since $\ell\geq 3$ terms are suppressed
relative to this by factors of $d_{T}/d_{X}\lesssim 0.01.$ We can thus
rephrase our search for Planet X as a search for unexplained $\ell=2$
potential terms. There are five real degrees of freedom in the tidal field,
which we can express as the five terms in the traceless, symmetric $3\times 3$
tidal matrix, or as the real (for $m\geq 0$) and imaginary (for $m>0$)
coefficients $A_{2m}$ of the tidal potential

$\psi(r,\theta,\phi)=\frac{4\pi
GM_{X}r^{2}}{5d_{X}^{3}}\sum_{m=-2}^{2}A_{2m}(\theta_{X},\phi_{X})Y_{2m}(\theta,\phi).$
(1)

A single point source like Planet X has only three degrees of freedom (e.g.
$M_{X},\theta_{X},\phi_{X}$) because the principal components of its tidal
matrix must be in the ratio $-1:-1:2.$

Our method will be to construct the Fisher matrix for the 5 terms of the tidal
field. Then the perturbation for a Planet X in any posited direction
$(\theta_{X},\phi_{X})$ can be expressed as a linear superposition of these,
and we can infer a bound on $M_{X}/d_{X}^{3}$ for this position. The Fisher
matrix of tidal spherical harmonics is equivalent to that derived for Planet
Xs placed at five distinct positions on the sky, so we opt for computing the
latter because it is more convenient to obtain derivatives with respect to
point masses than quadrupole coefficients within the dynamical codes. In
Appendix A we show how this matrix is used to place limits on $M_{X}$ over the
full sky. Because all of the gravitational effects of Planet X scale as
$M_{X}/d_{X}^{3},$ we will quantify all our results as the measurement
accuracy $\sigma_{M}$ attained on $M_{X}$ at fixed nominal distance
$d_{X}=400$ AU. Then we can rescale $\sigma_{M}$ to any distance. In
particular, we will consider a measurement to be “conclusive” for the presence
of Planet X if it can detect $M_{X}=5M_{\oplus}$ at $d_{X}=800$ AU at
$5\sigma$ significance, which is equivalent to obtaining
$\sigma_{M}<0.125M_{\oplus}$ at the nominal distance.

The critical parameter describing the shifts, precessions, and extraneous
oscillations that the tidal field generates is the ratio of the tidal
acceleration to the solar acceleration:

$\epsilon\equiv\frac{M_{X}d_{T}^{3}}{M_{\odot}d_{X}^{3}}\approx 3\times
10^{-11}\left(\frac{d_{T}}{5.2\,\textrm{AU}}\right)^{3}\left(\frac{M_{X}}{5M_{\oplus}}\right)\left(\frac{400\,\textrm{AU}}{d_{X}}\right)^{3}.$
(2)

We are clearly justified in calculating our derivatives as first-order
perturbations under this parameter. To give a sense of scale, the total
displacement over an observing period $t$ for constant tidal acceleration for
the nominal Planet X is

$\Delta x\approx\frac{1}{2}a_{\rm tidal}t^{2}=\frac{1}{2}\epsilon
a_{\odot}t^{2}\approx
500\,\textrm{m}\,\left(\frac{d_{T}}{5.2\,\textrm{AU}}\right)\left(\frac{t}{12\,\textrm{yr}}\right)^{2}$
(3)

and the observed angular shift if this is transverse to the line of sight from
the barycenter is independent of the tracer’s distance:

$\Delta\theta=\frac{\Delta x}{d_{T}}\approx
130\,\mbox{$\mu\textrm{as}$}\left(\frac{t}{12\,\textrm{yr}}\right)^{2}.$ (4)

These maximal signals should be detectable at $S/N=$10–100 from either radio
ranging observations of planets (if radial) or the expected minor-planet
optical astrometry.

The constant-acceleration model is, however, only applicable to distant bodies
which traverse a small fraction of an orbit during the LSST campaigns, e.g.
trans-Neptunian objects (TNOs). It is only roughly applicable to the Cassini
data, which span about 45% of one Saturnian period. But it is very difficult
to disambiguate a Planet X tidal signal from other gravitational influences
using such localized measurements. Using more bodies improves our ability to
constrain the spatial/temporal variation, and hence the source, of a
gravitational anomaly.

For data spanning one or more periods of the tracer particles, it is easier to
understand the influence of the tidal fields as a combination of several
effects enumerated below. These grow over time more slowly than $t^{2}$ and
have different scalings with target heliocentric distance $d_{T}.$

####  2.1.1 Radial displacement

One manifestation of the $A_{20}$ field is a perturbation to Kepler’s 3rd Law,
i.e. an offset in radius at fixed orbital period $P$ or mean angular motion
$\Omega_{0}.$ For a nominal Planet X at an ecliptic pole, this is

$\delta r=\frac{\epsilon d_{T}}{3}\approx
9\,\textrm{m}\,\left(\frac{d_{T}}{5.2\,\textrm{AU}}\right).$ (5)

While this distance shift is well above the uncertainties in ranging for Mars
or the outer planets, it requires determining the period of the monitored
object to equally high precision. Radial shifts are observable from astrometry
only through the reflex of Earth’s motion, which suppresses the signal by
${(1\,\textrm{AU})/d_{T}},$ making this likely undetectable from LSST
astrometry. We can alternatively view this effect as changing the period
associated with a fixed semi-major axis, by an amount $\epsilon P/2.$ This
would lead to an apparent angular lag that grows linearly with time, equal to

$\Delta\phi=2\pi\,\frac{\epsilon}{2}\frac{t}{P}\approx
20\mbox{$\mu\textrm{as}$}\left(\frac{tP}{(12\,\textrm{yr})^{2}}\right)\left(\frac{M_{X}}{5M_{\oplus}}\right)\left(\frac{400\,\textrm{AU}}{d_{X}}\right)^{3}.$
(6)

If the distance for a body is determined to sub-meter precision by ranging,
then the period shift would eventually become detectable with
$\mu\textrm{as}$-level astrometry. But we are not aware of any impending
observations that would yield both ranging and astrometric data of the
required precision on the same body.

####  2.1.2 Precession

Another manifestation of the tidal field is a shift of the epicyclic
frequencies away from the orbital frequency, generating a precession of the
ascending node (for the $\hat{z}$ epicycle, largely undetectable by ranging)
and longitude of perihelion. The latter is accessible to ranging via
$\hat{\rho}$ epicycles, and also to astrometry via $\hat{\phi}$ epicycles,
which have amplitudes typically twice as large (in meters) as the radial
oscillations (an advantage for astrometry over ranging). The typical
precession rate is $\epsilon\Omega_{0},$ which leads to perturbations whose
amplitudes grow linearly over multiple orbital cycles, leading to transverse
astrometric displacements (angular shifts in any direction) of

$\Delta\theta\approx\sin(\Omega_{0}t)\left\\{\begin{array}[]{llll}\sin{i}\,\epsilon\Omega_{0}t&\approx
43\,\mbox{$\mu\textrm{as}$}\,\sin{i}&\left(\frac{d_{T}}{5.2\,\textrm{AU}}\right)^{3/2}\left(\frac{t}{12\,\textrm{yr}}\right)&\textrm{(nodal
precession)}\\\ 2e\,\epsilon\Omega_{0}t&\approx
86\,\mbox{$\mu\textrm{as}$}\,e&\left(\frac{d_{T}}{5.2\,\textrm{AU}}\right)^{3/2}\left(\frac{t}{12\,\textrm{yr}}\right)&\textrm{(apsidal
precession)}\\\ \end{array}\right.$ (7)

The signal accessible to ranging is primarily the $\hat{\rho}$ component of
apsidal precession

$\Delta\rho\approx\sin(\Omega_{0}t)\,e\,\epsilon\Omega_{0}td_{T}\approx
160\,\textrm{m}\,e\,\left(\frac{d_{T}}{5.2\,\textrm{AU}}\right)^{5/2}\left(\frac{t}{12\,\textrm{yr}}\right).$
(8)

These precession signals all grow in amplitude toward the outer solar system;
the choice of which minor planet populations are most useful for detecting
Planet X then becomes a compromise between the increasing signal and the
smaller number of bright targets (higher measurement noise) as we move
outwards. There are $\sim 1200$ Trojans with typical projected LSSTastrometric
accuracy below 5 mas per component per observation, and $1.6\times 10^{5}$
MBAs that would satisfy the same criterium. The $130\times$ larger number of
targets among the MBAs is partly offset by the $2^{3/2}$ larger precession
signal of the Jovian Trojans. We will investigate the value of the Jovian
Trojans for constraints on Planet X, but it is possible that the MBAs will
offer stronger constraints if the systematic errors associated with smaller
bodies—namely radiation pressure—can be controlled. We defer investigation of
the MBAs to a future publication.

Jupiter and Saturn have very regular orbits ($e\approx 0.05,$ $\sin i<0.02$
from the invariable plane) whereas the Trojans have an RMS eccentricity of
0.08 and an RMS $\sin i$ of 0.28, more than $10\times$ larger than Jupiter and
Saturn. This makes the precession signals from a polar Planet X much larger in
this minor body population. Mars’ eccentricity of 0.093 is also more favorable
to observations of apsidal precession than those of Jupiter and Saturn.

####  2.1.3 Quadrupole excitations

For non-polar Planet X locations, the $A_{2m}$ terms are non-zero at $m=1$
and/or 2 and produce acceleration components $\approx a_{\rm tidal}\cos m\phi$
around the tracer orbit. These drive observable oscillating perturbations. The
quadrupole in particular is a distinctive signature of the tidal field. Its
characteristic amplitude, expressed as a physical or angular displacement is,
for a tracer with period $P$,

$\displaystyle\Delta x$ $\displaystyle\approx a_{\rm
tidal}\Omega_{0}^{2}/4=\frac{\epsilon d_{T}}{4}\approx
6\,\textrm{m}\,\left(\frac{d_{T}}{5.2\,\textrm{AU}}\right)^{4}$ (9)
$\displaystyle\Delta\theta$ $\displaystyle\approx\Delta
x/d_{T}=\frac{\epsilon}{4}\approx
1.7\,\mbox{$\mu\textrm{as}$}\,\left(\frac{d_{T}}{5.2\,\textrm{AU}}\right)^{3}.$
(10)

These stronger $d_{T}$ dependences favor the use of the Trojans over the MBAs,
and are of similar order to the collective accuracy expected from LSST
observations of Trojans, and much larger than spacecraft ranging accuracies.

Quadrupolar oscillations are also excited by less exotic gravitational
sources, namely the quadrupole moments of the collective TNO and MBA
populations, or the most massive individual bodies in these populations. Main-
belt quadrupoles should be distinguishable from Planet X quadrupoles by their
different radial dependence across the giant-planet region ($r^{-3}$ vs
$r^{1}$ force scaling) and the rapid orbital periods of the MBAs. A TNO
quadrupole is, however, very similar to one from Planet X, and must be limited
either by prior constraint on its amplitude or by measuring some $\ell>2$
component that is known to be associated with the perturber’s mass
distribution.

From these rough calculations, we conclude that Planet X induces displacements
that are potentially detectable at high significance through several channels,
justifying the more thorough numerical investigation of the available
information. We also note that the signals from precession and from a shift in
the period-radius relation continue to grow linearly with time, so the $S/N$
from a monitoring program operated continuously over time $t$ would improve as
$t^{3/2}.$ The oscillatory quadrupole signal would gain $S/N$ more weakly, as
$t^{1/2},$ over multiple orbital periods of the tracer.

##  3 Mathematical methods and model

We estimate the ability of a set of observational data points
$\textbf{d}=\\{d_{1},d_{2},...,d_{k}\\}$—in our case the measurements of
angular positions or ranges to solar system bodies—to constrain a set of
parameters $\textbf{p}=\\{p_{1},p_{2},...,p_{N}\\}$ of a model. In our case
the parameters are the masses of a set of active solar system bodies, and the
phase space positions of active and tracer bodies at a reference epoch. The
parameters also include those governing various systematic effects in the
measurements, such as differential chromatic refraction for ground-based
astrometry. We are interested in constraints on a single parameter—the mass of
a putative PX—after marginalization over all the other “nuisance” parameters.
We will use the Cramér-Rao theorem (Cramér, 1946; Rao & Das Gupta, 1989),
which states that no unbiased estimator of the parameters can yield more
precise constraints than calculated from the Fisher matrix, as described
below.

###  3.1 The Fisher matrix

For a system with likelihood $\mathcal{L}(\textbf{d}|\textbf{p})$ of obtaining
the data given the parameters, the available information on our parameter set
is expressed by the $N\times N$ Fisher information matrix $\mathcal{F}$:

$\mathcal{F}_{ij}=\left<-\frac{\partial^{2}\log{\mathcal{L}}}{\partial
p_{i}\partial p_{j}}\right>$ (11)

Let us consider an unbiased estimator
$\theta=\\{\theta_{1},\theta_{2},...,\theta_{N}\\}$ of our parameter set p.
The multivariate Cramér-Rao inequality imposes the following restriction on
the covariance matrix of any such estimator:

$\text{Cov}(\theta)\geq\mathcal{F}^{-1}$ (12)

where the matrix inequality states that the difference matrix
$\text{Cov}(\theta)-\mathcal{F}^{-1}$ is positive semidefinite (e.g.
$[\textbf{v}^{T}(\text{Cov}(\theta)-\mathcal{F}^{-1})\textbf{v}]\geq 0$ for
any vector $\textbf{v}=\\{v_{1},v_{2},...,v_{N}\\}$). This condition allows us
to establish a lower bound for the variance of any individual estimator
$\theta_{i}$ of the parameter $p_{i}$.

If the data points d have independent gaussian errors
$\sigma_{\textbf{d}}=\\{\sigma_{1},\sigma_{2},...,\sigma_{k}\\}$ and means
$\mu_{\textbf{d}}=\\{\mu_{1}(\textbf{p}),\mu_{2}(\textbf{p}),...,\mu_{N}(\textbf{p})\\}$,
the Fisher information becomes

$\mathcal{F}_{ij}=\sum_{l=0}^{k}\left<-\frac{\partial^{2}}{\partial
p_{i}\partial
p_{j}}\left(-\frac{(d_{l}-\mu_{l}(\textbf{p}))^{2}}{2\sigma_{l}^{2}}\right)\right>$
(13)

which simplifies to

$\mathcal{F}_{ij}=\sum_{l=0}^{k}\frac{\partial\mu_{l}(\textbf{p})}{\partial
p_{i}}\frac{\partial\mu_{l}(\textbf{p})}{\partial
p_{j}}\frac{1}{\sigma_{l}^{2}}$ (14)

The lower bound for the unmarginalized variance of $\theta_{i}$ is

$\text{Var}(\theta_{i})\geq(\mathcal{F}_{ii})^{-1}$ (15)

The process of marginalizing over all parameters $p_{j}$ for which $j\neq i$
weakens the derivable constraints $p_{j}$. Instead of simply inverting the
$ii$ element, the whole matrix is inverted, and then the $ii$ element of the
inverse Fisher matrix is selected. The bound becomes

$\text{Var}(\theta_{i})\geq(\mathcal{F}^{-1})_{ii}.$ (16)

Two additional properties of the Fisher matrix will be exploited here. First,
the Fisher matrices of statistically independent experiments can be summed to
express the information in their joint constraints. For example, we can
incorporate knowledge of a planet’s mass from satellite tracking by adding
$1/\sigma^{2}_{M}$ to the appropriate diagonal element of $\mathcal{F}$.
Second, we can generalize the equation for marginalization in Equation (16) to
the case where we partition the parameters into a retained subset of interest
$\textbf{p}_{R}$ and nuisance parameters subset $\textbf{p}_{N}.$ The
marginalization over nuisance parameters transforms the Fisher matrix for the
full parameter set into a smaller one over the retained parameters via

$\mathcal{F}=\left(\begin{array}[]{cc}\mathcal{F}_{RR}&\mathcal{F}_{RN}\\\
\mathcal{F}_{NR}&\mathcal{F}_{NN}\end{array}\right)\longrightarrow\mathcal{F}_{RR}-\mathcal{F}_{RN}\mathcal{F}_{NN}^{-1}\mathcal{F}_{NR}.$
(17)

We construct a model for solar system dynamics and use it to compute the
derivatives of solar system observables with respect to our model parameters.
From those derivatives and from the uncertainties of the observational data
points, we compute the Fisher matrix with Equation (14). The inference of
solar system parameters from these observations should closely match the
Fisher bound, because (a) the uncertainties on the observations are generally
well described by Gaussians, and (b) the dynamical model is very well
constrained: any departures from a nominal model (i.e. with $M_{X}=0$) are in
the regime where linear perturbation theory will suffice. Under these
conditions, the likelihood is formally Gaussian, and the Cramér-Rao bound is
known to be saturated, i.e. the attained uncertainty equals the lower bound.

The greater threat to the accuracy of Fisher forecasts is the possible
omission from the model of physical or observational effects that are
covariant with $M_{X}$ in the posterior likelihood. We will aim to incorporate
any such effects into our model.

###  3.2 Dynamical model and parameters

Our model for solar-system observations is a numerical integration of the
Newtonian equations of motion, detailed in Section 4.1. We are at liberty to
ignore relativistic effects since they will not significantly alter the
derivatives needed for the Fisher matrix in Equation (14). Only if the Fisher
matrix is near singularity will small perturbations to the derivatives be
important. As long as we assume standard General Relativity (GR), there is no
additional freedom to the model introduced by relativistic terms. Likewise we
can ignore the effects of gravitational lensing, stellar aberration, or other
small but known perturbations to the observed positions.

The presence of a Planet X can be modeled through a free mass parameter
$M_{X}$ at a fixed distance $d_{X}$ from the Solar System barycenter. The
effect of Planet X on the observable members of the solar system over the
scale of a few decades can be decomposed into a constant acceleration (which
is unobservable because it does not create any differential motion between
observers and targets), a tidal acceleration scaling as $M_{X}/d_{X}^{3}$, and
higher-order multipoles that are too small to be detectable with the data
under consideration. A general tidal field is described by a traceless
symmetric $3\times 3$ tensor, or by the Hermitian coefficients of the $\ell=2$
spherical harmonics—either way, 5 real degrees of freedom. It is convenient to
introduce five nominal positions for Planet X, assigning to each a different
mass parameter $M_{X1},\ldots,M_{X5}$, each located $d_{X}=400$ AU from the
barycenter, at selected angular positions listed in Table 1. Since any tidal
field—including the tidal field of a Planet X at any position on the sky—can
be expressed as a linear combination of the tidal fields generated by these 5
masses, the $5\times 5$ covariance matrix of these five free parameters can be
used to derive the expected uncertainty on any posited distant masses. The
mathematics for this operation are derived in Appendix A. Our choice of the 5
point-mass locations is not unique; any 5 locations for which the
transformation matrix between point masses and quadrupole coefficients is
well-conditioned would work.

Table 1: Positions of canonical Planet X’s

| $X$ | $Y$ | $Z$
---|---|---|---
$\boldsymbol{\hat{n_{1}}}$ | 1 | 0 | 0
$\boldsymbol{\hat{n_{2}}}$ | 0 | 1 | 0
$\boldsymbol{\hat{n_{3}}}$ | $+\sqrt{3}^{-1}$ | $+\sqrt{3}^{-1}$ | $-\sqrt{3}^{-1}$
$\boldsymbol{\hat{n_{4}}}$ | $+\sqrt{3}^{-1}$ | $-\sqrt{3}^{-1}$ | $+\sqrt{3}^{-1}$
$\boldsymbol{\hat{n_{5}}}$ | $-\sqrt{3}^{-1}$ | $+\sqrt{3}^{-1}$ | $+\sqrt{3}^{-1}$

Note. — The unit vectors (in ecliptic coordinates) pointing to the positions
of the five fiducial Planet X’s introduced to our dynamical model. The tidal
field of a perturber at any location in the sky can be expressed as a linear
sum over these five, as detailed in Appendix A.

All other parameters of the model are ultimately considered nuisances in the
pursuit of $M_{X}.$ The baseline model contains the following:

  * •

For each of the eight known planets we consider seven free parameters: the six
state vector components at the reference epoch $t_{0}$, taken to be MJD
$=60000$ (25.0 Feb 2023), and the mass $M$, all given for the barycenter of
the planet plus satellites. The masses of the known planets (and of some main
belt asteroids and TNOs) are constrained through close encounters with
spacecraft or tracking of natural or artificial satellites. More precisely, it
is $GM$ that is well constrained, but we will take $G$ as fixed and place all
uncertainties in the masses, which is mathematically equivalent. We take this
information into account by adding independent gaussian priors on each
planet’s mass to our Fisher matrix. The adopted priors are listed in Table 2.

  * •

The Sun has a free mass $M_{\odot}$. There is a degeneracy in our Fisher
matrix in that any overall shift or velocity of the solar system barycenter is
unobservable. Marginalization over the state vector of the barycenter is
mathematically equivalent to treating the initial state vector of any single
body as fixed. We therefore choose to treat the state vector of the Sun as
known.

  * •

We include the $J_{2}$ gravitational moment of the Sun as a free parameter. To
establish an adequate prior, we consider the range of $J_{2}$ values derived
by Mecheri et al. (2004) through helioseismology. We then adopt
$\sigma_{J2}=0.05$.

  * •

We ascribe a free mass parameter to each of the 343 largest main-belt
asteroids. These are the same asteroids that are active in the JPL ephemerides
DE430 and DE431 (Folkner et al., 2014). We fix their initial state vectors to
the values in these ephemerides. For 20 asteroids which have precise masses
determined from close encounters, we introduce individual Gaussian priors to
the mass parameters. For Ceres and Vesta, the priors are listed in Table 2.
The other 18 asteroid priors are taken from the uncertainty of the weighted
averages from Baer & Chesley (2008), from which we only leave out 189 Phthia,
that is not included in our 343-asteroid set. For the remaining 323 asteroids,
we assign a nominal mass based on their $H$ values, and introduce parameters
for deviations from the nominal mass. One general rescaling parameter
multiplies all the masses simultaneously, and an additional parameter for each
asteroid modulates its individual mass. The global parameter has a prior with
standard deviation of 20% from unity. The individual masses have independent
priors with $\sigma$ at 50% of their value.

  * •

In addition to individually modeling the DE431 asteroids, we also model a set
of two rings placed in the main belt to account for mass contributed by
smaller asteroids. These rings lie in the ecliptic plane with radii
$2.06\mathrm{AU}$ and $3.27\mathrm{AU}$. Each ring’s mass is a free parameter,
to which Gaussian priors with $\sigma=3.3\times 10^{-5}M_{\oplus}$ are
assigned. The two radii mark the inner and outer bounds of the bulk of the
main belt. Using two rings introduces our uncertainty in the effective mass-
weighted mean radius of the smaller bodies. This value is chosen as a
compromise between the different belt masses proposed in the literature (e.g.
Carry, 2012; Pitjeva & Pitjev, 2018), an exact value is not needed as the
derivatives are insensitive to the mass assumed.

  * •

We include a free parameter for the mass of the Pluto system, leaving its
state vector fixed due to its low mass. The adopted prior for this parameter
is included in Table 2.

  * •

We model the azimuthally symmetric portion of the Kuiper belt as the sum of
two uniform rings of mass in the ecliptic plane at $39.5\,\mathrm{AU}$ and at
$43.0\,\mathrm{AU}$ from the solar system barycenter, with distances chosen to
correspond with the Plutinos at the 3:2 mean motion resonance with Neptune,
and the classical belt’s rough center. Each is given a nominal mass of
$10^{-2}M_{\oplus}$, based on estimates by Pitjeva & Pitjev (2018). Since the
true Kuiper belt mass is quite uncertain, we attribute to each mass a Gaussian
prior with a $\sigma$ of 50% of its value. Again, the use of two rings
introduces uncertainty in the effective radius of the TNOs. In particular this
means that ratio of the $r^{3}$ and $r^{4}$ potentials that it induces on
bodies interior to it is left free to vary.

###  3.3 Degeneracies with trans-Neptunian quadrupoles

Our assumed ring model for the Kuiper belt ignores the graininess and
asymmetries in the TNO distribution. These features can generate stochastic
tidal fields that are indistinguishable on decade timescales from a Planet X
tidal field on observations interior to Neptune. In order to place priors on
such quadrupoles and incorporate them into our results, we use the L7
synthetic debiased model from the CFEPS project (Kavelaars et al., 2009; Petit
et al., 2011; Gladman et al., 2012) to recreate the evolution of TNO
distribution features in the Neptune co-rotating frame over a time range of
$2\times 10^{4}$ years. We sample the five components of the internal
quadrupole field in 50 year intervals, and we measure the $5\times 5$
covariance matrix of the quadrupole components from these samples. Variations
are of the order of $10^{-3}M_{\text{CFEPS}}/(40\,\textrm{AU})^{3}$, where
$M_{\text{CFEPS}}$ is the total mass of the population, estimated at
$1.5\times 10^{-2}M_{\oplus}$. The covariance matrix is then transformed into
our basis of five point masses at 400 AU, and it is added to the covariance
matrix obtained through the Fisher matrix method.

###  3.4 Mass uncertainties from large TNOs

The CFEPS sample reaches magnitudes as bright as $H\approx 3.4$. We can
reasonably suppose that all Kuiper Belt objects brighter than this limit are
known or will be shortly. Thus, the additional uncertainty in the internal
gravitational field depends on how well their masses are determined. We add to
the quadrupole covariance matrix the contribution of the 15 brightest TNOs,
except Pluto which was modeled separately. (The potential existence of larger
and more distant TNOs of gravitational significance can be considered the
subject of this work.) For Kuiper belt objects (KBOs) with mass determination
from moons, we adopt the published system-mass uncertainties from such data.
For the others, we assign an uncertainty equal to half the mass estimated from
their radius measurement and a presumed density. These uncertainties are shown
in Table 2.

Over periods of time much shorter than a KBO period, the quadrupole
gravitational field generated by KBOs is essentially indistinguishable from a
Planet X field from any observations of test bodies interior to Neptune. This
leads to an irreducible floor on the sensitivity of searches for gravity from
Planet X. We will refer to the combined uncertainties from individual large-
KBO masses, plus the stochastic tidal fields of smaller bodies described in
the previous subsection, as the collective “KBO floor” on $\sigma_{M}$. The
contributions of stochastic fields dominate this uncertainty, and are likely
to continue to set a limit on detectability of Planet X tidal fields for many
decades.

Table 2: Priors on $GM$ of planets Planet | $\sigma_{GM}/GM_{\oplus}$ | Reference
---|---|---
Mercury | $3.0\times 10^{-9}$ | Konopliv et al. (2020)
Venus | $1.5\times 10^{-8}$ | Konopliv et al. (1999)
EarthaaDerived from constraints on the Earth-Moon mass ratio. | $3.3\times 10^{-8}$ | Pitjeva & Standish (2009)
Mars | $1.3\times 10^{-8}$ | Konopliv et al. (2011)
Jupiter | $2.1\times 10^{-5}$ | Durante et al. (2020)
Saturn | $2.8\times 10^{-6}$ | Jacobson et al. (2006)
Uranus | $1.1\times 10^{-5}$ | Jacobson (2014)
Neptune | $2.5\times 10^{-5}$ | Jacobson (2009)
Ceres | $9.9\times 10^{-10}$ | Konopliv et al. (2018)
Vesta | $3.0\times 10^{-11}$ | Konopliv et al. (2014)
Pluto | $5.2\times 10^{-6}$ | Stern et al. (2015)
Eris | $1.4\times 10^{-5}$ | Holler et al. (2021)
Makemake | $1.0\times 10^{-5}$ | presumedbbWe assume that future data on the Makemake system, such as the tracking of its known moon MK2, will constrain its total mass to a similar level of uncertainty as current constraints on Gonggong and Haumea.
Haumea | $6.7\times 10^{-6}$ | Ragozzine & Brown (2009)
Gonggong | $1.2\times 10^{-5}$ | Kiss et al. (2019)
Quaoar | $1.5\times 10^{-5}$ |  Fraser et al. (2013)ccFraser et al. (2013) gives four possible orbit solutions for the Quaoar system. We assume this ambiguity will be solved during the following years.
Orcus | $3.2\times 10^{-7}$ | Grundy et al. (2019b)
Sedna | $8.6\times 10^{-5}$ |  Pál et al. (2012)ddNo dynamical mass measurement. Mass uncertainty assumed to be half the total mass estimate given its radius and $\rho\approx 2$ g/cm3. Reference is for object radius.
2013 FY27 | $3.9\times 10^{-5}$ |  Sheppard et al. (2018)ddNo dynamical mass measurement. Mass uncertainty assumed to be half the total mass estimate given its radius and $\rho\approx 2$ g/cm3. Reference is for object radius.
Varda | $1.1\times 10^{-6}$ | Grundy et al. (2015)
2014 UZ224 | $2.2\times 10^{-5}$ |  Gerdes et al. (2017)ddNo dynamical mass measurement. Mass uncertainty assumed to be half the total mass estimate given its radius and $\rho\approx 2$ g/cm3. Reference is for object radius.
2002 AW197 | $4.0\times 10^{-5}$ |  Vilenius et al. (2014)ddNo dynamical mass measurement. Mass uncertainty assumed to be half the total mass estimate given its radius and $\rho\approx 2$ g/cm3. Reference is for object radius.
G!kún$||$’hòmdímà | $5.5\times 10^{-7}$ | Grundy et al. (2019a)
2002 TX300 | $2.1\times 10^{-6}$ | Elliot et al. (2010)
Ixion | $3.1\times 10^{-5}$ |  Levine et al. (2021)ddNo dynamical mass measurement. Mass uncertainty assumed to be half the total mass estimate given its radius and $\rho\approx 2$ g/cm3. Reference is for object radius.

##  4 Numerical methods

###  4.1 Determination of the derivatives

Variational derivatives for each body were determined by adding bodies and
effects to a base simulation consisting of 9 massive bodies (the sun and 8
planets). Our simulation begins at MJD 60000.0, with the initial states and
masses of the 9 massive bodies queried from the JPL Horizons
service.111https://ssd.jpl.nasa.gov/horizons/app.html

Beginning with this base simulation, we compute variational derivatives of
planetary positions with respect to their $(8\times 7)+1=57$ parameters using
REBOUND’s variational equations module (Rein & Tamayo, 2016). All simulations
are integrated from the starting epoch for 13 years, using the 15th order
adaptive ias15 integrator Rein & Spiegel (2015); Everhart (1985) implemented
in REBOUND (Rein, H. & Liu, S.-F., 2012). We choose ias15 over symplectic
integrators like WHFast since our integrated time span is relatively short (no
longer than a few orbital periods of the Trojan asteroids). A forward
integration through 2036 covers all future observations considered in this
paper, and a backward integration to 1965 covers historical data we use.

We add relevant bodies and forces to the system to calculate derivatives with
respect to other sources of gravity. First, we compute derivatives with
respect to the Planet X parameters by adding a fiducial Planet X along each
one of the five axes described by the unit vectors in Table 1, with nominal
mass of $M_{X}=5M_{\oplus}$ and circular orbit at $d_{X}=400\,\mathrm{AU}.$
Motion of Planet X during the integration has negligible effect, so the
initial state vectors of the Planet X’s are held fixed, leaving five free
parameters to represent Planet X (equivalent to the five free parameters of a
general tidal field). The nominal mass $M_{X}$ is chosen at the low end of the
range suggested by Brown & Batygin (2021).

Next we augment our base simulation by adding 343 massive main belt asteroids.
We obtain the masses of these from Table 12 in Folkner et al. (2014);222Note
that there are 344 columns in table in (Folkner et al., 2014), as the asteroid
499-Hamburga was listed twice. these asteroids were used in the computation of
the DE430 and DE431 JPL ephemerides. The states of these asteroids at the
starting epoch were queried from JPL Horizons. The size of these simulations
(more than $350$ massive bodies, in addition to several thousand test
particles once we examine the Trojan asteroids in Section 7) necessitates
parallelization of the integration. We perform this task by dividing the
Trojan asteroids and massive main belt asteroids into $10$ and $30$ groups
respectively. We then compute variational derivatives for each combination of
every subgroup of main belt asteroids and Trojans, at all Planet X positions,
and compiled the results. Since each individual simulation only contains a
subset of massive asteroids, this means that the trajectories of the bodies in
each subgroup are slightly inconsistent with each other. Fortunately, the main
belt asteroids contribute only a small perturbation to the base simulation of
Trojans and planets. The addition of a typical single asteroid causes a
relative change in the variational derivative of another asteroid of $\leq
10^{-8}$, so this small inconsistency can be neglected for the purposes of
making Fisher-matrix derivatives.

The potential from each ring modeled in Section 3.2 felt by a particle of mass
$m_{p}$ at radius $r_{x}$ can be written as

$\phi=-\frac{GM}{2\pi}\int_{0}^{2\pi}\frac{1}{r_{0}^{2}+r_{x}^{2}+z^{2}-2r_{0}r_{x}\cos{\theta_{0}}}d\theta_{0}.$
(18)

where $r_{0}$ is the radius of the ring and $M$ is its mass. The force due to
this potential is added to the integrator using the REBOUNDx package (Tamayo
et al., 2019) and the derivatives are computed by finite differencing two
simulations with slightly different ring masses.

Finally, we compute derivatives with respect to the $J_{2}$ of the Sun again
by finite differencing two simulations with a slightly different value of
$J_{2}$ about an initial value taken of $2.2\cdot 10^{-7}$ (e.g. Pijpers,
1998). The Sun’s north pole is tilted from the north ecliptic pole by
$16.12^{\circ}$ in right ascension and $-2.70^{\circ}$ in declination, which
can be computed using coordinates of the solar north pole (Archinal et al.,
2018) and the obliquity of the earth (e.g. Capitaine et al., 2003).

In the solar frame, where the $\hat{z}$ axis aligns with the solar pole, the
force from a $J_{2}$ can be written as

$\vec{F}=\frac{3GM_{\odot}J_{2}R_{\odot}}{2r_{p}^{5}}\left[\left(\frac{5z_{p}^{2}}{r_{p}^{2}}-1\right)(x_{p}\hat{x}+y_{p}\hat{y})+\left(\frac{5z_{p}^{2}}{r_{p}^{2}}-3\right)(z_{p}\hat{z})\right].$
(19)

where $x_{p}$,$y_{p}$,$z_{p}$ are the cartesian distances from the sun along
each axis, and $r_{p}=\sqrt{x_{p}^{2}+y_{p}^{2}+z_{p}^{2}}.$ We add this force
(after rotating into the ecliptic plane) into REBOUND using the REBOUNDx
package (Tamayo et al., 2019). Only a single free parameter is involved, since
the orientation of the Sun’s rotation pole is known to within about
$0.1^{\circ}$ (Beck & Giles, 2005).

For a given observer our Cartesian positions can be converted to an observable
position on the sky through $\alpha=\arctan(y(t)/x(t))$,
$\delta=\arcsin(z(t)/\sqrt{x(t)^{2}+y(t)^{2}+z(t)^{2}})$, where $\alpha$ and
$\delta$ are the right ascension and declination of the body and $x(t)$,
$y(t)$ and $z(t)$ are all Cartesian differences between the coordinates of the
body in question at the time the signal is emitted and the coordinates of the
observatory at the time of observation (e.g. $x(t)=x_{p}(t-\tau)-x_{obs}(t)$,
for some light-time delay $\tau$). Conversion between our computed Cartesian
state derivatives and the astrometric derivatives can be computed with the
Jacobian. From the cartesian state positions and its derivatives, we also
compute derivatives of round-trip light time observables, in order to
incorporate spacecraft ranging into our Fisher analysis.

Table 3: Planetary observations assumed for fitting Planet | Observatory/SpacecraftaaMGS: Mars Global Surveyor; MRO: Mars Reconnaissance Orbiter | TypebbA: astrometry; R: range | Time Span | Typical uncertainty $/\sqrt{N}$ ccN: number of data points. For the astrometric data with different errors on ra and dec, uncertainty listed is the largest one.
---|---|---|---|---
Mercury | Messenger | A | 2008–2009 | 0.5 mas
| Messenger | R | 2008–2009 | 20 m
| Messenger | R | 2011–2015 | 0.02 m
Venus | Magellan/VLBI | A | 1990–1994 | 0.4 mas
| Cassini | A | 1998–1999 | 3 mas
| Cassini | R | 1998–1999 | 3 m
| Venus Express/VLBI | A | 2007–2013 | 0.3 mas
Mars | Mariner 9 | R | 1971–1972 | 3 m
| Viking | R | 1976–1982 | 0.5 m
| Phobos 2 /VLBI | A | 1989 | 3 mas
| Pathfinder | R | 1997 | 4 m
| MGS | R | 1999–2006 | 0.002 m
| MGS/VLBI | A | 2001–2003 | 0.3 mas
| Odyssey | R | 2002–2013 | 0.008 m
| Odyssey/VLBI | A | 2002–2013 | 0.03 mas
| MRO | R | 2006–2014 | 0.02 m
| MRO/VLBI | A | 2006–2013 | 0.03 mas
| Odyssey (presumed) | R | 2014–2034 | 0.006 m
| MRO (presumed) | R | 2014–2034 | 0.009 m
Jupiter | Pioneer 10/VLBI | A | 1973 | 100 mas
| Pioneer 10 | R | 1973 | 4000000 m
| Pioneer 11/VLBI | A | 1974 | 30 mas
| Pioneer 11 | R | 1974 | 3000 m
| Voyager 1/VLBI | A | 1979 | 30 mas
| Voyager 1 | R | 1979 | 3000 m
| Voyager 2/VLBI | A | 1979 | 3 mas
| Voyager 2 | R | 1979 | 3000 m
| Ulysses/VLBI | A | 1992 | 5 mas
| Ulysses | R | 1992 | 20 m
| Galileo/VLBI | A | 1996–1997 | 2 mas
| Cassini | A | 2000 | 80 mas
| Cassini | R | 2000 | 4000 m
| Juno | R | 2016–2017 | 7 m
| Juno (presumed) | R | 2018–2025 | 3 m
Saturn | Pioneer 11 | R | 1979 | 2000000 m
| Pioneer 11 | A | 1979 | 600 mas
| Voyager 1 | A | 1980 | 100 mas
| Voyager 1 | R | 1980 | 50 m
| Voyager 2 | A | 1981 | 200 mas
| Voyager 2 | R | 1981 | 100 m
| Cassini/VLBI | A | 2004–2011 | 0.6 mas
| Cassini | R | 2004–2017 | 0.2 m
Uranus | Voyager 2 | A | 1986 | 10 mas
| Voyager 2 | R | 1986 | 200 m
Neptune | Voyager 2 | A | 1989 | 4 mas
| Voyager 2 | R | 1989 | 30 m
Outer | Nikolaev Obs. | A | 1966–1998 | 4 mas
Outer | La Palma Obs. | A | 1984–1998 | 3 mas
Outer | Bordeaux Obs. | A | 1985–1996 | 10 mas
Outer | Flagstaff Obs. | A | 1995–2015 | 0.9 mas
Outer | Table Mountain Obs. | A | 1997–2013 | 1 mas

##  5 Results from major-planet tracking

In this section we forecast the Planet X constraints available from the set of
published historical observations, with some extension to presumed future
spacecraft data. In Section 7, we will consider the addition of a set of
predicted LSST observations of 7664 Trojan asteroids.

###  5.1 Assumed observations

The available historical data consists of astrometric measurements, which give
the angular position of the planet or one of its moons, and radar ranging
measurements, which give the light time between the observing station on Earth
and a spacecraft at the planet (typically referenced to the planetary
barycenter). The observations we use are presented in Table 3, and are
comprised of similar data to those compiled by Park et al. (2021). We restrict
ourselves, however, to the data currently available on the JPL
website333https://ssd.jpl.nasa.gov/planets/obs_data.html, plus the complete
set of Cassini data points available on the INPOP Astrometric Planetary
Database website444http://www.geoazur.fr/astrogeo/?href=observations/base.

We posit that Juno observations of Jupiter’s range will continue up to 2025,
simulating data points with the same errors and spacing in time as the actual
Juno data from 2016 and 2017. This corresponds to an addition of $21$ data
points with spacing of $4-5$ months. For Mars we assume that two probes will
be operating at all times between 2014 and the end of LSST observations,
taking range measurements with similar spacing and precision as those from
Odyssey and Mars Reconnaissance Orbiter. These observations are marked as
presumed in Table 3. The extension adds a total of $5943$ data points for MRO
and $16065$ for Odyssey. Both numbers are found by multiplying the number of
data points taken in $2013$ by $21$ years.

The assumed priors on masses of individual major and minor planets, summarized
in Table 2, are independent of the positional data in Table 3, since they are
derived solely from tracking of natural or artificial satellites around each
mass.

###  5.2 Nominal Case

We start with our “nominal” case, where we use the full set of published
historical planetary data, augmented by our presumed data for Mars and Jupiter
only up to 2021. This simulates data which should already be in hand. The top
panel of Figure 1 plots the measurement noise $\sigma_{M}$ of a Planet X at
400 AU vs position in the sky for this nominal case.

For each observational scenario we consider, we construct such an all-sky map,
and from it extract several summary statistics given in Table 4. We list the
minimum, median, and maximum values of $\sigma_{M}$ across the sky. We also
tabulate the percentage of the sky $f_{400}$ for which $\sigma_{M}<1,$ meaning
that the presence of a $5M_{\oplus}$ Planet X at a distance of 400 AU would be
detected unambiguously, at $5\sigma$ significance. This is 99.2% for the
nominal case, i.e. a Planet X at the near end of its expected range would be
found anywhere on the sky. A more stringent criterion, demanding $5\sigma$
detection of $5M_{\oplus}$ at a distance $d_{X}=800$ AU
($\sigma_{M}<0.125M_{\oplus}$ on our plots) is satisfied over $f_{800}=4.8\%$
of the sky in the nominal case.

Figure 1: Maps of forecasted uncertainty $\sigma_{M}$ on the mass of Planet X at 400 AU as a function of its location in the sky. Plots are in ecliptic coordinates with the Vernal equinox at center. The Mollweide projection is used. The upper panel shows the _nominal_ case using all observational constraints to date, including Juno/Cassini ranging and Mars ranging. “Old” data are the existing measurements without Juno, Cassini, or Mars ranging. In the middle panel, only the old data and the Juno/Cassini data are used (Mars ranging is omitted), as described in Section 5.3. In the lower panel, forecasted LSST observations of Jovian Trojans are included, which makes up some of the precision lost by omission of Mars (see Section 7.1). Table 4: Limits on Planet X mass in multiple scenarios Comment | Data included |  $\sigma_{M}$ for Planet X ($M_{\oplus}$) |  $f_{\rm sky}$ detecting $5M_{\oplus}$aaFraction of the sky over which a $5M_{\oplus}$ Planet X would be a $\geq 5\sigma$ detection.
---|---|---|---
| J/CbbRanging data from Juno and Cassini. | Mars | Trojans | Min | Median | Max | 400 AU | 800 AU
Nominal | $\times$ | $\times$ |  | 0.12 | 0.22 | 1.07 | 99.2% | 4.8%
…published ranging | $\times$ | $\times$ |  | 0.12 | 0.24 | 1.27 | 94.0% | 1.3%
…ranging through 2034 | $\times$ | $\times$ |  | 0.11 | 0.20 | 0.80 | 100.0% | 8.1%
Old only |  |  |  | 9.11 | 12.38 | 43.41 | 0.0% | 0.0%
…$+$Juno/Cassini | $\times$ |  |  | 0.16 | 0.35 | 2.87 | 87.3% | 0.0%
…$+$Mars |  | $\times$ |  | 1.16 | 1.54 | 4.23 | 0.0% | 0.0%
Nominal$+$Trojans | $\times$ | $\times$ | $\times$ | 0.12 | 0.22 | 0.98 | 100.0% | 5.0%
…Trojan errors$\div 5$ | $\times$ | $\times$ | $\times$ | 0.11 | 0.20 | 0.58 | 100.0% | 7.3%
…Trojan errors$\div 20$ | $\times$ | $\times$ | $\times$ | 0.11 | 0.14 | 0.27 | 100.0% | 29.2%
Old$+$Trojans, no sys |  |  | $\times$ | 1.66 | 1.90 | 3.40 | 0.0% | 0.0%
Old$+$Trojans |  |  | $\times$ | 1.84 | 2.15 | 3.88 | 0.0% | 0.0%
…$+$Juno/Cassini | $\times$ |  | $\times$ | 0.15 | 0.33 | 2.07 | 89.2% | 0.0%
…$+$Mars |  | $\times$ | $\times$ | 1.02 | 1.22 | 2.39 | 0.0% | 0.0%
KBO floorccLimits purely from degeneracies with uncertainties in Kuiper Belt quadrupole. |  $\cdots$ |  $\cdots$ |  $\cdots$ | 0.053 | 0.066 | 0.071 | 100.0% | 100.0%

####  5.2.1 Extended ranging

Figure 2 shows the effect of restricting the Mars, Jupiter and Saturn ranging
data to observations already published (dotted line), or extending the ranging
with all presumed data up to 2025 for Juno, and up to 2034 for two Mars
probes. We compare, for each, scenario, the fraction of the sky for which each
level of detection is possible. Progressing from
published$\rightarrow$current$\rightarrow$predicted data adds 5–15% of sky
coverage attaining any particular $\sigma_{M}.$

Figure 2: Fraction of the sky for which the $\sigma_{M}$ of a Planet X at 400
AU is below the value on the $x$ axis. This is shown for the nominal scenario
(solid), one restricted to already-published ranging, and another
extrapolating ranging data through 2034.

####  5.2.2 Degraded ranging

We also consider the scenario in which the uncertainties of the most precise
ranging data are currently underestimated due to systematic errors. In Figure
3 we plot the ranges of $\sigma_{M}$ constraints as we inflate the measurement
errors on either the Juno/Cassini ranging or on the Mars ranging. The
$\sigma_{M}$ constraints are seen to be only very weakly dependent on the
accuracy of the Mars ranging, but are more sensitive to the Juno/Cassini
uncertainties.

Figure 3: As we increase the measurement errors on the Juno/Cassini ranging
(red) or Mars ranging (blue), the forecasted uncertainties $\sigma_{M}$ on the
Planet X mass at 400 AU are degraded. The solid lines show the median value of
$\sigma_{M}$ and the hatched regions bound the 5–95% percentiles of
$\sigma_{M}$ across the sky. The full power of the J/C data are more critical
to the Planet X constraints than that of the Mars ranging.

###  5.3 Data subsets

To determine the relative value of different planets’ data to the Planet X
constraints, and to explore whether complementary sets of data could confirm
each other’s detections, we split the nominal data into three distinct
subsets: First, the Juno$+$Cassini ranging data; second, the Mars ranging
data; and third, the remaining elements of the compilation summarized in Table
3. We will refer to this last subset as the “old” observations (even though
some of the Mercury and Venus data are recent).

We find that the quality of $\sigma_{M}$ constraints improves as we move from
Old, to Old$+$Mars, to Old$+$J/C, to the full Old$+$Mars$+$J/C nominal data.
This is illustrated by the solid curves in Figure 4, which plots the fraction
of the full sky achieving a given level of $\sigma_{M}$ in each scenario.
Numerical data on these scenarios are also given in Table 4. Here it is clear
that the J/C data are significantly more useful than the Mars data, however
their combination is needed if we wish to approach secure detections of the
most pessimistic $5M_{\oplus},$ $d_{X}=800$ AU model across a substantial
portion of the sky.

The middle panel of Figure 1 shows the sky map of $\sigma_{M}$ in the
Old$+$J/C case, illustrating the weaker (but still interesting) level of
constraint compared to the upper, nominal case.

Figure 4: Each curve plots the fraction of the sky ($y$ axis) for which the
uncertainty $\sigma_{M}$ on Planet X at 400 AU is at or below the value on the
$x$ axis. The solid curves are for cases without LSST Trojan data, and are
labelled by the subsets of the observational data they use. All curves make
use of the “old” observations. The dashed curves add the use of LSST Trojan
observations to the scenario of the solid curve of matching color. The
leftmost “KBO floor” curve plots the constraints that would be attained if the
only limitation were confusion with the uncertain tidal field produced by the
members of the Kuiper belt. This floor is up to an order of magnitude below
constraints with the best data expected in the next decade. The two vertical
lines mark the $\sigma_{M}$ values that are needed for a definitive $5\sigma$
detection of Planet X if it has mass of $5M_{\oplus}$ and is at distance of
800 AU (left) or 400 AU (right). Ideally we would be able to detect the former
at 100% of locations on the sky.

###  5.4 Kuiper Belt quadrupole limit

The portion of the Kuiper Belt which was treated separately from the Fisher
matrix, that is, the uncertainty on the tidal field caused by asymmetries in
population distribution and specific point sources whose mass are uncertain,
amount to a fundamental limit, computed as an additive contribution to the
covariance matrix. If this additive term is considered by itself, we achieve
the results labeled as “KBO floor” in Figure 4 and in Table 4. Figure 5 is a
sky map of the $M_{X}$ uncertainties at 400 AU induced purely by our ignorance
of the Kuiper Belt mass distribution. Note that a polar Planet X is better
constrained in this case than one in the ecliptic plane, unlike our other
cases where the pole is the weak spot.

Figure 5: A map of $\sigma_{M},$ the lower limit on achievable $M_{X}$
uncertainty for Planet X at 400 AU, set by ignorance of the Kuiper Belt’s
tidal field. Note the change in color scale from Figure 1.

##  6 Trojans as test particles

In this and following sections we investigate the value of LSST tracking of
Jovian Trojans as constraints on Planet X. A similar analysis could
incorporate other minor body populations, especially the main belt (MBAs), but
we limit our current scope to the Trojans. Future work could also extend this
analysis to gauge the utility of LSST-scale minor planet precision astrometry
for other gravitational signals, e.g. determination of masses of members of
these populations.

###  6.1 Measurement uncertainties

####  6.1.1 Astrometric precision per observation

Our model for uncertainties in individual astrometric measurements of minor
planets by LSST assumes optimal extraction in the face of three sources of
random error. Faint sources are dominated by uniform shot noise from diffuse
foregrounds and backgrounds. Brighter sources become dominated by the shot
noise from their own photons. The brightest sources are limited by uncorrected
stochastic atmospheric deflections, which we will place at a nominal level
$\sigma_{\rm min}=2$ mas per axis.

The LSST survey simulation (Ivezić, 2022) predicts a set of characteristics of
each observation of a run-through of the entire survey. We make use of the
tabulated values of the sky brightness (converted to counts per exposure per
unit sky area $b$), the equivalent Gaussian width $\sigma_{\rm psf}$ of the
point-spread function (PSF), and the filter band $Q\in\\{g,r,i,z\\}$ of each
exposure (we ignore the $u$ and $Y$ bands as unlikely to add useful
constraints). Estimated values of the magnitude zeropoints $m_{Q,0}$ for each
filter—i.e. the magnitude generating a single photocarrier during the
exposure—are taken from Jones (2016). We simplify the survey simulator’s
weather predictions by assuming that all nights with $<1\%$ cloud cover have
no cloud extinction, and discarding observations with higher cloud cover.

The assumed absolute magnitude $H$ of the target is combined with the its
geocentric and heliocentric distances $d_{g}$ and $d_{h}$ to yield a true
apparent magnitude in filter $Q$ of

$m_{Q}=H+5\log_{10}{(d_{g}d_{h})}+c_{Q}$ (20)

where $c_{Q}$ is the typical $Q-V$ color of the population. We have ignored
the illumination-phase term, since the phase angle never exceeds
$\sin^{-1}\frac{1\,\textrm{AU}}{d}\approx 11^{\circ}$ for Jovian Trojans, and
will not significantly impact the overall $M_{X}$ sensitivity. The photon
count is then given by

$f=10^{0.4(m_{Q,0}-m_{Q})}.$ (21)

The best achievable astrometric uncertainty for an unresolved source is
derived using Fisher matrix algebra in Appendix B. Reproducing Equation (B6),
a circular Gaussian PSF with dispersion $\sigma_{\rm PSF}$ on each axis yields
an astrometric uncertainty on each axis of:

$\sigma_{x}^{-2}=\frac{f}{\sigma_{\rm{PSF}}^{2}}\int_{0}^{\infty}u\frac{e^{-2u}}{e^{-u}+\nu}du$
(22)

where $\nu=\frac{2\pi\sigma_{\rm{PSF}}^{2}b}{f}.$ Experience with the Dark
Energy Survey (DES) suggests that PSF-fitting astrometry can closely approach
these theoretical limits.

The uncalibrated portion $\sigma_{\rm min}$ of the atmospheric turbulence
displacements then must be added in quadrature to the $\sigma_{x}$ value of
Equation (22) to yield the final $\sigma_{x}$ of the observation. Fortino et
al. (2021b) demonstrate that calibration of DES images based on the position
of Gaia stars can bring this error to $\sigma_{\rm min}\approx 2$ mas. We
consider this as our nominal value, though we also evaluate the effect that
this value has on our final results in Section 7.2.

####  6.1.2 Trojan catalog

The list of all known Jovian Trojans was downloaded from the Minor Planet
Center database555Date of download: 07-Oct-2021. (Minor Planet Center, 2021),
giving orbital elements and $H$ values for each. The known sources exhibit a
power-law size distribution, i.e. exponential in $H$, for $H\leq 14.5$, and we
presume the observations are nearly complete to this magnitude. For $H>14.5$,
where the known counts depart from the exponential, we generate artificial
Trojans by replicating the orbital elements of randomly selected known Trojans
and assigning them new $H$ values until we match the extrapolation of the
exponential distribution.

This extrapolation ends up being irrelevant for the final $M_{X}$
uncertainties because, as we will show below, most of the information is
carried by Trojans in the regime where the observational catalog is nearly
complete.

####  6.1.3 LSST Observing Plan

We use the Baseline 2.0 simulated 10-year LSST survey produced by the Survey
Cadence Optimization Committee (Ivezić, 2022) which lists the times,
pointings, and observing conditions of every exposure during the simulated
survey. The effects of varying airmass are reflected in the values of expected
sky brightness and PSF size in the table, so the simulation is realistic in
these regards, and the simulation also includes cloud conditions based on
weather history at the site. The total amount of cloudless exposures with
$Q\in\\{g,r,i,z\\}$ is $\sim 1.2\times 10^{6}$.

Because the orbital period of the Trojans is 12 years, we assume with modest
optimism that LSST tracking of Trojans will continue for 2 years past the
nominal 10-year LSST duration. This also makes the results insensitive to the
actual start date of the survey. To implement this, we duplicate the pointings
and weather from years 9 and 10 to create years 11 and 12 of our nominal
survey.

At this point we have a catalog of all of the useful observations $j$ of
Trojan target $i$ over the full 12-year survey, with the times of each and the
per-axis measurement error $\sigma_{ij}$ for the object’s position. We will
assume that these errors are statistically independent, which will be true for
errors from shot noise and atmospheric turbulence.

####  6.1.4 Information vs depth

A useful measure of the constraining power of observations of some minor body
population is the collective uncertainty $\sigma_{\rm pop}$ defined through
this sum over the astrometric uncertainties $\sigma_{ij}$ per axis in
observations $j$ of bodies $i$ with absolute magnitude $H_{i}$ brighter than
some threshold:

$\sigma_{\rm pop}^{-2}(<H)=\sum_{i:H_{i}<H}\sum_{j}\sigma^{-2}_{ij}.$ (23)

This quantity gives the RMS size of a perturbation that can be detected (per
axis) at $1\sigma$ using astrometry of the full population, in the absence of
systematic measurement errors, or other free parameters that are covariant
with the perturbation.

Figure 6 shows this quantity vs the limiting magnitude $H$ for the Trojan
population. This sum asymptotes to $\approx 3\,\mbox{$\mu\textrm{as}$}$ at
$H<14.5$ for the Trojans. For the remainder of our analysis, we will restrict
our analysis to the real Trojan sample (with $H<14.5$) leaving 7664 Trojans,
which will be treated as massless test particles in further dynamical
analyses.

Figure 6: Cumulative information from predicted LSST observations of the
Trojans as a function of the magnitude of the faintest Trojan considered. Top
axis shows object radius of objects with albedo of 0.04. Trojans fainter than
$H=14.5$ (vertical dashed line) do not add significant information—this is
roughly the point at which their astrometry becomes limited by background
noise.

Taking nominal distances and albedos of 5.2 AU and 0.04 for Trojans, the lower
limit to radius $R_{T}$ for Trojans that would contribute significantly to the
full-population constraints is roughly 4 km. An error of $\sigma_{\rm min}=2$
mas subtends $\approx 8$ km, twice the minimum radii of targets of interest.
At the faint end of the useful $H$ range, however, shot noise is dominant over
atmospheric turbulence, so the per-visit centroiding errors are larger than
the angular sizes of the objects. For brighter Trojan targets, the measurement
errors per epoch are smaller than their radii. In either case, any source of
_stochastic_ noise (uncorrelated between observations) becomes unimportant if
it is well below the $\sigma_{ij}$ value. At ground-based seeing of
0$\farcs$5–1$\farcs$0, all of the minor-planet targets will be unresolved, and
trailing in $<30$ s LSST exposures will negligibly degrade the astrometric
accuracy, so we can assume point-source values for measurement errors.

The $\sigma_{\rm pop}$ errors are far below the sizes of the individual
targets (just as in the case of radio ranging to the major planets). Both
astrometry and ranging require careful determination of the location of the
observatory relative to the Earth-Moon barycenter. An advantage of conducting
tests in the outer solar system is that the desired gravitational signals from
Planet X are meter-scale or greater, which is well above geodetic errors.
Signals at cm scales place greater stress on our knowledge of the Earth-Moon
system.

###  6.2 Systematic errors

Other forms of measurement error or un-modelled sources of acceleration can
potentially degrade the Planet X constraints. If these errors are _systematic_
, i.e. are correlated with the model’s derivatives with respect to a free
model parameter, then they must be controlled at the $\sigma_{\rm pop}$ level,
i.e. sub-$\mu\textrm{as}$ and sub-meter scale. We next discuss the systematic
errors that we expect to encounter, their expected sizes, and mitigation
strategies.

####  6.2.1 Radiation pressure

If the ratio of incident solar radiation pressure to the solar gravitational
force on a spherical target of radius $R_{T}$, density $\rho$, and
heliocentric distance $d_{T}$ is well below $\epsilon,$ then we can ignore its
influence. This condition is

$\displaystyle\epsilon$ $\displaystyle\gg\frac{\pi R_{T}^{2}L_{\odot}/4\pi
c}{4\pi GM_{\odot}\rho R_{T}^{3}/3}=\frac{3L_{\odot}}{16\pi GM_{\odot}\rho
c}R_{T}^{-1}$ (24) $\displaystyle\Leftrightarrow\quad R_{T}$ $\displaystyle\gg
17\,\textrm{km}\times\left(\frac{5.2\,\textrm{AU}}{d_{T}}\right)^{3}\left(\frac{1000\,\textrm{kg}\,\textrm{m}^{-3}}{\rho}\right).$
(25)

The $H<14.5$ limit of useful LSST Trojans corresponds to $R_{T}>4$ km at a
typical albedo of 0.04 for the population (Fernández et al., 2003), suggesting
the smallest useful Trojans have incident radiation pressure a few times
larger than the nominal tidal effect. Intrinsic variations in albedo (and
hence $R_{T}$) and density at fixed $H$ will make it difficult to predict this
on a per-target basis. Departures from isotropy in the reflected or re-
radiated energy likewise cannot be predicted accurately, i.e. the Yarkovsky
effect. The radiation-pressure perturbation is too small, however, to be
detectable for any single Trojan, and therefore we need only be concerned
about the average effect over the whole population, which could indeed be
statistically significant at the fainter $H$ values of the proposed sample.
When averaging over the population and over rotation phases, we would expect
the dominant component to be a radial inverse-square force, which does not
induce precession or quadrupole signals that would mimic Planet X. It could,
however, potentially confound the radial displacement signal—but we do not
expect that this signal is strong enough to contribute to Planet X constraints
from LSST. We will test this conclusion numerically by allowing the
population-averaged radiation pressure to have a free nuisance parameter
included in the Fisher analysis.

The angular momentum vector of a body’s orbit or rotation breaks the symmetry
of the system around the sun-target axis, potentially producing non-radial
forces—the Yarkovsky effect. For (101955) Bennu, the asteroid with the best-
constrained Yarkovsky models, the results of Farnocchia et al. (2021) are that
the tangential Yarkovsky acceleration is roughly $g_{Y}=10\%$ of the incident
radiation-pressure acceleration for this body with radius 0.25 km. If we
assume a similar $g_{Y}$, i.e. anisotropy of thermal emission, then the ratio
of non-radial radiation-pressure acceleration $a_{Y}$ to the nominal Planet X
tidal acceleration becomes

$\frac{a_{Y}}{a_{\rm tidal}}=\frac{3g_{Y}L_{\odot}}{16\pi R_{T}\rho
c}\frac{d_{X}^{3}}{GM_{X}d^{3}_{T}}=\left(\frac{1.8\,\textrm{km}}{R_{T}}\right)\left(\frac{g_{Y}}{0.1}\right)\left(\frac{\rho}{1000\,\textrm{kg}\,\textrm{m}^{2}}\right)\left(\frac{5.2\,\textrm{AU}}{d_{T}}\right)^{3}.$
(26)

For the smallest members of the Trojan sample, the Yarkovsky force is only a
few times below the Planet X signal, and must be considered in an analysis on
a population-averaged basis. Note that the population’s vectorial average of
$a_{Y}$ will be reduced because of the variation of spin axes, so the
Yarkovsky signal of the Trojans could be too small to be detectable. Two
characteristics of Yarkovsky acceleration will distinguish it from Planet X
perturbations. First, if there is a tangential force, it has non-vanishing
curl, which is an effect that gravity cannot produce. More importantly, it
will vary systematically with $H$ (as a proxy for $R_{T}$). We will be able to
regress the inferred tidal field against $H$ to determine and correct for the
trend. Introducing a nuisance parameter for the Yarkovsky effect will quantify
the effectiveness of this technique.

The non-gravitational forces on free-flying spacecraft are orders of magnitude
too large for use as gravity probes at our desired level. Only spacecraft that
can be referenced to the barycenters of planets (either in orbit or on fly-
bys) are useful, until such time as ranging of drag-free spacecraft becomes
available.

####  6.2.2 Outgassing

Another potential non-gravitational force is reflex from anisotropic mass loss
(the “rocket effect” Whipple 1950). Visible comae are very rare in both the
Trojan and MBA populations. We may safely assume that the average member of
either population has existed for $>4$ Gyr and is currently losing mass at a
rate below its own mass per 4 Gyr. If we assume that the velocity of ejected
mass is below 500 m/s, roughly the mean thermal velocity of water at $\approx
200$ K , then the total acceleration from outgassing or mass loss is bounded
by

$\frac{g_{R}}{M}\frac{dp}{dt}<g_{R}\frac{v_{\rm
thermal}}{4\,\textrm{Gyr}}\approx 4g_{R}\times
10^{-15}\,\textrm{m}\,\textrm{s}^{-2}.$ (27)

Here $g_{R}$ is a measure of the anisotropy of the mass loss for the rocket
effect, in analogy to the $g_{Y}$ for the Yarkovsky effect. There will again
be averaging due to the random orientation of the spin/rocket vectors among
the population. A single fixed-orientation planar surface outgassing uniformly
into a hemisphere would have $g_{R}=0.5.$ Once we vectorially average over the
different orientations of patches on a given body, the rotation of the body,
and the multiplicity of bodies, an average of $g_{R}<0.2$ seems likely, at
which point the maximal mass-loss force is $10\times$ below the nominal tidal
force. We are probably safe ignoring mass-loss forces, especially if any
objects showing coma in deep LSST coadded images are excluded.

####  6.2.3 Photocenter motion

The displacement between the photocenter and barycenter of a tracer becomes an
error in the ephemeris calculation. We can decompose the photocenter motion
(PCM) into some apparent transverse displacement $\Delta x(\phi_{0})$ that is
purely a function of the Sun-Target-Observer phase angle $\phi_{0}$ and lies
in the plane of this angle; plus a zero-mean component that fluctuates due to
the rotation of a body lacking rotation symmetry. The latter will be
uncorrelated with any gravitational signal and can be viewed as additional
astrometric noise. The former will be coherent with an oscillation at the
body’s synoptic period. The phase angle $\phi_{0}$ is bound to
$\pm\sin^{-1}\frac{1\,\textrm{AU}}{d_{T}}.$ In this range we will approximate

$\Delta x(\phi_{0})=k\phi_{0}R_{T}$ (28)

for some constant $k$ that will depend on the geometry and scattering
properties of the tracer body. As illustrated in Figure 7, the effect of this
synoptic PCM is to move the apparent center of the object toward the Sun by
$kR_{T}.$ PCM is a more severe issue for larger bodies, in contrast to
radiation-pressure systematics.

Figure 7: If the body at right with radius $R$ is illuminated along the
horizontal axis, and the photocenter displacement when viewed (in red) at an
illumination angle of $\phi$ is $k\phi R,$ the naive use of the photocenter
instead of the barycenter is equivalent to a shift of the body toward the Sun
by $kR,$ at the location of the black dot.

The PCM shifts can be quite large. For a specular surface, the value of $k$ is
$1/2.$ Appendix C shows that a Lambertian sphere has $k=3/8.$ Mallama (1993)
estimates the PCM for the Galilean satellites using photometric models of
their surfaces, with results equivalent to $k\approx 0.4.$ These values
suggest that the angular PCM for even the smallest Trojan targets will be
$\approx 4\,\textrm{km}\,\times 0.4/5.2\approx 300$ m and manifest as $\gtrsim
1.5$ km under-estimates of inferred radial coordinates. This is several orders
of magnitude above both $\sigma_{\rm pop}$ and the nominal Planet X
signatures. Indeed our expected $\sigma_{\rm pop}$ for the Trojans is roughly
$10^{-3}$ of the diameter of the smallest targets (Figure 6).

We have reason, however, to be optimistic that this systematic error will not
degrade the Planet X constraints. As noted, PCM generates an apparent radial
displacement, but not the precession or multipole signals that carry most of
the tidal information. Furthermore in Section 6.3 we will examine the impact
of allowing every Trojan test particle to have a free parameter for its $k$
value.

Tanga et al. (2022) demonstrate direct detection of the PCM wobble at $\approx
1$ mas level from the rotation of (21) Lutetia. PCM values this high would
inflate our expected statistical error floor of $\approx 2$ mas. But at
$H=7.3,$ and distance 2.4 AU, Lutetia should have a much larger angular
stochastic PCM than the typical LSST Trojan or MBA target. For present
purposes we can ignore the stochastic rotational portion of PCM, and concern
ourselves only with the rotation-averaged values.

####  6.2.4 Differential chromatic refraction

Images of celestial objects appear closer to the zenith because of refraction
in the atmosphere. Because the refraction is wavelength-dependent and the
observing bandwidth is finite, there is a differential chromatic refraction
(DCR) that depends on the precise spectrum of the source. To first order, the
spectral dependence can be quantified by some broadband color measure $c_{bi}$
for each source $i$ spanning band $b.$ For a horizontally stratified
atmosphere, the DCR in exposure $j$ at zenith angle $z_{j}$ can be quantified
as an apparent shift

$\Delta z=k_{bj}c_{bi}\tan z_{i}$ (29)

where $k$ scales the amplitude of the effect. Values of $k$ measured in the
$griz$ bands by Bernstein et al. (2017) are 45, 8.4, 3.2, and 1.4 mas per
magnitude of $g-i$ color, respectively. Chromaticity of the glasses in the
prime-focus corrector cause an effect (“lateral color”) of similar size and
wavelength dependence, which is oriented radially from the telescope axis.
These effects are clearly large, i.e. $10^{4}\times$ larger than the target
$\sigma_{\rm pop}$ for the LSST Trojans. Exquisite removal of the effect will
be needed.

The DCR is easily estimated to precision $<1$ mas on every exposure by using
color data for each individual target, combined with knowledge of the $k_{bj}$
gained by monitoring reference stars of varying color. In this way we can
insure that uncorrected DCR does not increase the stochastic astrometric error
floor $\sigma_{\rm min}.$ The much greater challenge is in insuring that
coherent residuals to the DCR corrections do not impact our inferences being
made at $\mu\textrm{as}$ level from the full population. At this level, an
accurate prediction of DCR would need both precise measurement of the time-
varying behavior of the atmosphere (the $k_{bi}$) and of the full spectral
energy distribution of each minor planet. It should be possible to infer the
former by using the tens of thousands of Gaia stars that will be imaged at
high $S/N$ in every LSST exposure—indeed we should be able to average many
exposures, since the atmospheric constituents and pressure do not change
rapidly. The $z_{i}$ values are of course known to very high accuracy. The
difficulty would be in knowing the relevant color factor for each body. We
will therefore include $c_{bi}$ as free parameters for every band and target
when making our forecasts (Section 6.3). The DCR effects can be distinguished
from gravitational perturbations because the former are modulated by the
pointing patterns of the telescope, i.e. we always know that the atmospheric
DCR points toward the zenith and the lateral color toward the telescope axis.
The Fisher matrix analysis will tell us how effectively this decoupling works.

####  6.2.5 Timing accuracy

The maximum apparent motion of a typical Trojan is $\approx 10$ mas/s. To keep
ephemeris errors below $\sigma_{\rm pop}\approx 2\,\mbox{$\mu\textrm{as}$},$
any systematic errors in the temporal midpoint of the LSST exposures should be
$\ll 200$ microseconds. The time to sweep the shutter across the LSST focal
plane is $O(1)$ s, so a very precise map of the shutter blade motion is
needed. Fortunately the LSST shutter motion is tracked at $O(10)$ microsecond
precision, but effort will be needed to map this onto the focal plane at this
accuracy.666A. Roodman, private communication. The tolerance for random errors
in the shutter timing is much looser; but one could imagine systematic
variations that could correlate with gravitational signals. For example, the
gravity vector on the shutter blades could be correlated with the Trojan’s
position on the sky and season, since the parallactic angle between celestial
north and altitude vectors at a given time of night depends on these factors.
The faster-moving MBAs could potentially be used as calibrators of the shutter
trajectory, or perhaps the rarer but even faster near-Earth asteroids.

###  6.3 Full model and parameters

In this subsection we detail the model and parameters for the observations of
the Trojans. The techniques would be applicable to other small-body
populations as well in future analyses.

For each of the Trojan asteroids, a.k.a. the targets or tracers, we begin with
six free parameters: its position components $x$, $y$, and $z$, and its
velocity components $v_{x}$, $v_{y}$, and $v_{z}$ at a reference epoch
$t_{0}$, taken to be MJD$=60000$ (25.0 Feb 2023), near the start of the
nominal LSST. We neglect the effect of individual Trojans’ gravity fields on
each other and on the remaining bodies. Therefore, their masses are not
included as free parameters.

####  6.3.1 Trojan observations and derivatives

The Trojans function as test particles in the modeled solar system
gravitational field. The computation of the derivatives of their observable
positions is done as described in Section 4.1, that is, converting the state
vector derivatives given by the REBOUND simulation into derivatives of sky
coordinates.

####  6.3.2 Photocenter motion

The photocenter motion displacement is computed from Equation (28), with its
direction given by the phase position angle $\phi_{PA}$. This gives us

$\Delta\delta=kR_{T}\phi_{0}\cos{\phi_{PA}};\qquad\Delta{(\alpha\cos{\delta})}=kR_{T}\phi_{0}\sin{\phi_{PA}}$
(30)

The $kR_{T}$ term is our free parameter, and derivatives with respect to it
are included in our Fisher matrix modeling for each Trojan.

####  6.3.3 Differential chromatic refraction

We model the differential chromatic refraction according to Equation (29),
with values of $k_{b}$ taken from Bernstein et al. (2017), since we assume
correction over atmospheric conditions can remove the systematic effect of
different observational conditions. The direction of the displacement is given
by the parallactic angle $p$, such that its components are

$\Delta\delta=k_{b}c_{bi}\tan{z_{i}}\cos{p};\qquad\Delta{(\alpha\cos{\delta})}=k_{b}c_{bi}\tan{z_{i}}\sin{p}$
(31)

The zenith and parallactic angles are computed from the simulated
observations. We end up with four color parameters $c_{bi}$ (one per filter)
for each Trojan $i$, each of them with derivatives to be included in the
Fisher matrix.

####  6.3.4 Radiation pressure

The radial acceleration $a_{R}$ due to solar radiation pressure is considered
to be proportional to $A_{r}/a^{2}$, with $A_{r}$ defined in Equation 32 and
$a$ being the semi-major axis.

$A_{r}=\frac{L_{\odot}R_{T}^{2}}{4Mc}$ (32)

The exact proportionality depends on the albedo and shape of the object.
Assuming a constant density for all Trojans, the proportionality becomes
$a_{R}\propto c_{\rm RP}/(a^{2}R_{T})$, where $c_{\rm RP}$ is now a nuisance
parameter that scales the whole population’s radiation pressure force. The
derivatives of the asteroid state vectors with respect to $c_{\rm RP}$ are
proportional to derivatives with respect to $M_{\odot}$. We obtain the former
by rescaling the latter according to the radius $R_{T}$ of each Trojan, which
is inferred from its magnitude $H$.

####  6.3.5 Yarkovsky effect

Finally, we consider the influence of the Yarkovsky effect on the asteroids.
As summarized in Bottke et al. (2006), the Yarkovsky effect can be separated
into two components, known as the seasonal effect (which acts along the spin
axis) and diurnal effect (which acts in a plane orthogonal to the spin axis).
For a population of asteroids with randomly oriented spin axes, there will be
no net effect from the diurnal effect. For the same population the seasonal
effect produces a net acceleration tangent to the orbit. The exact
acceleration experienced by an asteroid is a complicated function of its
shape, albedo, and orbit, but since our derivatives are insensitive to exact
value assumed, an order of magnitude estimate will suffice. If the bodies are
large (with diameters larger than $\approx 50$ meters) the thermal
acceleration can be computed using Equation (16) in Vokrouhlický & Farinella
(1998)

$a_{T}=-\frac{4\alpha}{3}\frac{A_{Y}}{a^{2}}A(\ell),$ (33)

where $\alpha$ is the albedo, $a$ is the semi major axis in AU, and $A(\ell)$
a function of the temperature distribution. $A_{Y}$ sets the scale of the
effect, and is a function of the mass and radius of the object being
irradiated. The overall scaling $A_{Y}$ is taken as a nuisance parameter. The
derivative with respect to $A_{Y}$ of each Trojan’s position is determined
numerically by adding to each test particle in the simulation a small
$\hat{\theta}$ acceleration with magnitude of
$10^{-12}\mathrm{m}\,\mathrm{s}^{-2}$, a rough value for a body with a radius
of $50$ m and a density of $\rho=3500\mathrm{kg}\,\mathrm{m}^{-3}$.
Derivatives of each body’s state vector with respect to this nominal
acceleration are rescaled to derivatives with respect to $A_{Y}$ using
Equation (33).

####  6.3.6 Parameter summary

There are 11 parameters for each Trojan—six phase space components, four
coefficients of DCR $c_{b}$, and a coefficient $k$ for the photocenter motion.
At a global level, there remain the masses of the gravitationally active
bodies, along with the solar $J_{2}$ moment and the state vector parameters of
the major planets, plus scaling parameters $c_{RP}$ and $A_{Y}$ for the radial
and Yarkovsky components of radiation pressure.

The total dimensionality of the Fisher matrix is thus $O(10^{5})$ once all the
Trojans are included. The matrix is, however, highly structured, which allows
a huge acceleration of the marginalization over nuisance parameters. We can
divide the parameter set into a group $p_{k}$ of 11 for the $k$th test body
(Trojan), and then a group $p_{G}$ of “global” parameters comprising the
specifications of all the active bodies, plus the global systematic parameters
for the Trojans. When ${\mathcal{F}}$ is broken into blocks among these
parameter subgroups, all of the blocks representing covariances between
$p_{k}$ and $p_{k^{\prime}}$ are zero for $k\neq k^{\prime}$. This means that
Equation (17) can be successively applied to all Trojans to marginalize down
to the global parameters:

${\mathcal{F}}\rightarrow{\mathcal{F}}_{GG}-\sum_{k}{\mathcal{F}}_{Gk}{\mathcal{F}}_{kk}^{-1}{\mathcal{F}}_{kG}.$
(34)

The marginalization over all Trojans’ parameters never requires inversion of
matrices larger than $11\times 11$ and there is never need to store the full
Fisher matrix.

##  7 Results with Trojan tracking

###  7.1 Combinations with ranging data

We now calculate the full-sky Fisher uncertainties $\sigma_{M}$ on Planet X
mass in scenarios including predicted LSST observations of Jovian Trojans. The
dashed lines in Figure 4 plot the sky fractions achieving a given $\sigma_{M}$
level when the LSST data are combined with the nominal major-planet data, in
the four possible choices of whether the Juno/Cassini and Mars data are
included. The summary statistics for these four scenarios are in the lower
half of Table 4.

The addition of the LSST astrometry for Trojans to the full Old$+$J/C$+$Mars
nominal observation set yields very minor improvements in Planet X
detectability. The biggest gain occurs near the ecliptic poles, where addition
of the Trojans lowers $\sigma_{M}$ from 1.07 to $0.98M_{\oplus}.$ The Trojan
information is, however, more useful if we have dropped one or both of the J/C
and Mars ranging sets from the constraints. Thus the LSST Trojans could be
valuable in building confidence in a detection from the nominal data by
confirming that the detection is robust to deletion of Mars or J/C data.

Adding the Trojans to the Old data reduces the median $\sigma_{M}$ on the sky
from 12.4 to 2.2$M_{\oplus}.$ Thus a $10M_{\oplus}$ Planet X would become
securely detectable over 20% of the sky completely independent of the modern
ranging data.

Adding LSST Trojan data to the Old$+$Juno/Cassini data raises the portion of
the sky amenable to $5\sigma$ detection of $M_{X}=5M_{\oplus}$ at 400 AU from
87% to 89%. The lower panel of Figure 1 shows this improvement relative to the
middle panel. One sees that the Trojans are particularly valuable in filling
in the weak coverage near the ecliptic poles. When Trojan data are added to
Old$+$Mars data, they improve the sensitivity to a polar Planet X by a factor
of 2.

###  7.2 Effect of $\sigma_{\rm min}$

A reduction on the turbulence floor $\sigma_{\rm min}$ of the LSST
observations from 2 mas to 0.5 mas does not yield drastic improvement on
constraining power in the nominal case with Trojans. The $\sigma_{M}$ near the
ecliptic poles is decreased from 0.98 to $0.89M_{\oplus}$, with lesser gains
elsewhere on the sky. This suggests that an astrometric limit of 2 mas already
extracts most of the Planet X information from the LSST exposures.

###  7.3 Suppression of systematic errors

The observational and physical Trojan systematic error parameters can be
crossed out of the Fisher matrix in order to generate an ideal case where
these parameters are fully known, i.e. there are no systematic errors. With
this test, we can examine the level to which the systematic effects are
degenerate with a Planet X signal. For this test we compare the Old$+$Trojans
scenario with and without the systematic error parameters free to vary,
omitting the J/C$+$Mars ranging to maximize the potential harm from Trojan
systematics. The rightmost two curves in Figure 8 plot this comparison, where
we find a reduction of 12% in the median value of $\sigma_{M}$ when these
effects are suppressed. Even though some of the systematic signals are much
larger than the Planet X perturbations, we can successfully exploit their
known dependence on target size or observing geometry to yield only minor
degradations in Planet X precision relative to the no-systematics case.

Figure 8: The fraction of the sky with uncertainty $\sigma_{M}$ on Planet X
mass below a given level is shown for observation scenarios including LSST
observations of Trojans. At the right end we see the relatively small
degradation in the combination of Trojans with old data that results from
marginalizing over parameters of systematic-error models (as we move from the
yellow to the dark curve on its right). At left are shown scenarios in which
we combine all of the proposed observational data, and arbitrarily reduce all
uncertainties on the Trojans by factors of 5 (green) and 20 (magenta) from the
actual forecast (black). It takes the 20-fold reduction for the Trojans to
make a substantial improvement in the Planet X constraints, and gain the
ability to detect a $5M_{\oplus}$ Planet X at $5\sigma$ significance over a
significant portion of the sky if it is at a distance of 800 AU. At far left,
below the 800 AU signal, is the floor set by stochastic KBO quadrupoles.

###  7.4 Avenues for stronger constraints

A fundamental limit in searching for the tidal field from unknown solar system
masses is the degeneracy with the tidal field produced by _known_ elements of
the solar system, namely Kuiper Belt members whose masses are poorly known,
and whose periods are so long that we cannot use time dependence on decade
scales to distinguish them from Planet X. In our Fisher matrix we include
these uncertainties in three forms. First, there is an uncertainty in the
azimuthally symmetric “ring” of the time-averaged Kuiper belt. This we find to
have relatively little effect on $\sigma_{M},$ perhaps because it only
generates the $m=0$ spherical harmonic, and any non-polar Planet X will have
components at $m=\pm 1$ and $m=\pm 2$. Second, we have included the tidal
fields from the largest individual KBOs. These are substantially reduced by
the fact that most large TNOs have moons, and hence we can obtain
uncertainties on their masses at $O(10^{-6}M_{\oplus}).$ Since the tidal field
scales as $M_{X}/d_{X}^{3},$ and the KBOs are at $d\geq 40$ AU while we are
calibrating $\sigma_{M}$ for $d_{X}=400$ AU, we expect the $\approx
10^{-5}M_{\oplus}$ uncertainties in the masses of the largest moon-less KBOs
to result in an $\approx 10^{-2}M_{\oplus}$ floor on $\sigma_{M}.$ We also
note that the tidal fields scale as $d^{-3}$ from their source, so additional
large TNOs that may exist beyond the Kuiper Belt ($d\gtrsim 50$ AU) are of
relatively little importance. We find the dominant contribution to the KBO
floor, plotted in the leftmost curve in Figure 4, is the third component, from
stochastic fluctuations in the quadrupole moment of the large numbers of
smaller KBOs that occur as individual elements orbit the Sun. Reducing this
uncertainty will unfortunately prove more difficult than quantifying only the
few largest KBOs, since it would require discovering and making reasonable
mass estimates for many objects all over the sky, something even LSST will not
do.

The good news is that this KBO floor will not be a problem until the data get
substantially better. The KBO floor at $\sigma_{m}\approx 0.06M_{\oplus}$ is
about 3 times smaller than the uncertainties forecast to be available from
observations in the coming decade even in the most favorable directions, so
there is room for improvement with better data. Furthermore the KBO floor lies
below the accuracy required to find the tidal field of a $5M_{\oplus}$ Planet
X at $d_{X}=800$ AU anywhere on the sky.

What observations would be needed to bring the achievable $\sigma_{M}$
constraints significantly closer to the KBO floor? The three middle curves of
Figure 8 (also in Table 4) consider cases in which we arbitrarily reduce all
of the measurement errors on LSST Trojans. It is apparent that the overall
sensitivity to Planet X begins to significantly improve when these errors are
reduced to $\lesssim 0.1$ those achievable with LSST. At $20\times$ reduced
errors, we can confidently detect the most pessimistic version of Planet X,
$5M_{\odot}$ at 800 AU, suggested by dynamical studies, over 29% of the sky.
The measurement uncertainties would still be significantly above the KBO
floor.

These 10–20$\times$ reductions over LSST capabilities are potentially
achievable over a few decades’ observations with an LSST-like facility, if its
observing time is more strongly concentrated on the Trojans instead of the
full sky, and if we note that the size of the tidally-induced precessions grow
linearly with time span. Ranging to future spacecraft missions to the giant
planets would obviously be of substantial gain, as would the advent of ranging
to drag-free spacecraft, e.g. with cm-scale accuracy to a test body beyond
Jupiter.

A less expensive route to higher accuracy on minor-planet positions is to
implement the stellar occultation detection array proposed by Rice & Laughlin
(2019), which can potentially reduce per-epoch positional uncertainties from
milliarcseconds to $O(10)$ $\mu\textrm{as}$ for occultations of brighter Gaia
stars. We leave a detailed quantitative study of this possibility for future
work.

##  8 Conclusions

Using Fisher matrix analysis, we find that current tracking of major planets
through range and astrometric measurements has achieved a precision such that
the tidal effect of an unknown Planet X can be distinguished from spurious
signals throughout most of the sky. Ranging data from Juno, Cassini and Mars
orbiters contribute most to the total Fisher information. In full-sky
searches, the fraction of the sky where the Cramér-Rao bound allows for a
$5\sigma$ detection of a $5M_{\oplus}$ Planet X is $f_{400}=99.2\%$ for Planet
X at 400 AU and $f_{800}=4.8\%$ at $800AU$.

How do actual attempts to detect Planet X in the data compare to our derived
bounds? Previous gravitational searches are still less conclusive than this
prediction. The strongest published bounds are from Fienga et al. (2020), who
perform a gravitational search for Planet X by perturbing the INPOP19a
ephemeris. Their most readily interpretable results are likelihoods of the
$\Delta\chi^{2}$ of ephemeris best fits with an added Planet X tidal field to
those from a null hypothesis with no Planet X present. Fienga et al. (2020)
test cases of $M_{X}=5$ and $10M_{\oplus}$ placed at $400<d_{X}<800$ AU, at
ecliptic latitudes between $\pm 25\arcdeg.$ Comparison with our forecasts is a
little clouded because they do not search the full sky, and their results do
not display the inversion symmetry and the $M_{X}/d_{X}^{3}$ scaling expected
for a tidal force [our Equation (1) and their Equation (1)], implying that
some numerical noise is present. Their Figure 5 and text suggest that a
$5M_{\oplus}$ Planet X at 600 AU is near a $3\sigma$ detection in most of
their footprint, equivalent to $\sigma_{M}\approx 0.5M_{\odot}$ in our
reference case of $d_{X}=400~{}AU$. A similar conclusion is drawn from their
depiction of $\approx 1\sigma$ signals from $M_{X}=5M_{\oplus}$ at 800 AU.
This constraint is 2 times weaker than the median $\sigma_{M}\geq
0.24M_{\oplus}$ bound that we derive for the comparable case of using
published ranging in Table 4. Thus their results do not violate our bound, and
it appears that stronger Planet X constraints should be possible with the
available data. They remark that Planet X perturbations were often compensated
in the fitting process by changes to the TNO ring mass, which points to a
degeneracy between these parameters. Our analysis shows, however, that a
Planet X signal can be sufficiently distinguished from a symmetric KBO ring
with current major planet data—suggesting that this is not currently a
limiting degeneracy. We find that quadrupolar asymmetries within the ring are
the stronger degeneracy with Planet X, but such perturbations are not present
in the Fienga et al. (2020) model. Also of note is the bound from Meisner et
al. (2017) which ruled out a Planet X brighter than $W1=16.66$ from a search
of the WISE spacecraft images. Using the models of Fienga et al. (2016), this
$W1$ limit corresponds roughly to excluding $M_{X}\geq 20M_{\oplus}$ at $622$
AU. Our forecasted gravitational bounds on $M_{X}$ are tighter than this
direct-imaging result. From this somewhat confusing state of affairs, we
conclude that further dynamical fitting and modelling is of substantial value.

We further investigate whether the predicted observations of Jovian Trojans
with the upcoming LSST survey would improve these constraints. The answer is
that, when all data on the major planets are considered, not much information
would be gained from these additional measurements, and a 10–20$\times$
improvement in their predicted uncertainty would be needed before they are
useful in this scenario. However, when different subsets of the major planet
data are ignored, there is a considerable gain of information from LSST Trojan
tracking. Therefore, these measurements would be useful as cross-checks for
the results obtained with specific major planet data sets. Our calculations
also indicate that astrometric data for a large population of small bodies can
internally solve for (and remove) the major sources of systematic error,
attaining the full meter-scale total statistical accuracy for true
gravitational signals.

Our discussion has been framed around the search for a single sub-Neptune-mass
planet lurking $\geq 400$ AU away, but the tidal-field constraints are equally
applicable to smaller or multiple bodies that are $\gtrsim 60$ AU away such
that their gravitational influence is dominated by a tidal field. For our
nominal case of using extant data to constrain the tidal field’s components,
Figure 2 indicates that $\sigma_{M}<0.7M_{\odot}$ for a single mass at 400 AU
over 90% of the sky. Using the scaling of Equation (1), this indicates that an
unknown Earth-mass planet should be detectable at $5\sigma$ significance for
$d_{X}\leq 260$ AU, and a Mars-mass object for $d_{X}\leq 120$ AU. There is
currently no evidence to exclude the existence of such bodies. The scaling of
a gravitational search for new masses, with signal $\propto d_{X}^{-3},$ is
more favorable than for searches for reflected solar light, which is $\propto
d_{X}^{-4}.$

The $5\times 5$ Fisher matrix we have derived for the tidal field is equally
applicable to cases in which a tidal field is generated from _multiple_
undiscovered bodies. Investigating the bounds on cases of multiple distant
planets would be a straightforward extension. Note that only a subset of the
space of 5-dimensional tidal fields is compatible with generation from a
single point mass, so there is a potential internal test for the multiplicity
of sources. Likewise one can estimate from the Fisher matrix our future
capability to localize any single source on the sky given its mass and
distance.

Performing these localizations will be computationally intensive in practice.
However, since the computation of the variational derivatives can be
parallelized, we anticipate no major computational bottlenecks associated with
processing a large future dataset. The principle challenge will be
incorporating physics with sufficient accuracy to compare simulations with
observations. Necessary improvements would include the $J_{2}$ geopotential
terms for several planets and a precise treatment of general relativistic
perturbations. As far as we know, no open source pipeline currently exists
that could preform this integration, so one will have to be developed for this
purpose.

Our results indicate which steps can be taken towards a better detection of
the Planet X signal. Fits of the major-planet data to integrated ephemerides
should be capable of stronger constraints than have been obtained, motivating
investigation of integrator accuracy and fitting techniques. Improved
efficiency for such integration and fitting processes could be needed. Further
characterization of the potential roles of other minor planets (e.g. MBAs or
Hildas) as test bodies is warranted. In the longer term, there is significant
improvement possible in our sensitivity to Planet X from spacecraft visiting
the giant-planet systems, ranging to drag-free spacecraft, and campaigns to
measure many thousands or millions of stellar occultations by minor planets.

The authors thank the NSF Astronomy & Astrophysics Program for supporting this
work: grant AST-2206194 for MJH and ZM, and grants AST-2205808 and AST-2009210
for GMB, DCH, and RCH. MJH is also grateful for support from the Smithsonian
Scholarly Studies Program and from the NASA YORPD Program (80NSSC22K0239). We
greatly thank the members of the LSST Survey Strategy Optimization Committee,
in particular Lynne Jones, for assembling the simulated LSST surveys and
helping us to use them properly.

## References

  * Aksnes et al. (1986) Aksnes, K., Franklin, F., & Magnusson, P. 1986, AJ, 92, 1436, doi: 10.1086/114280
  * Archinal et al. (2018) Archinal, B. A., Acton, C. H., A’Hearn, M. F., et al. 2018, Celestial Mechanics and Dynamical Astronomy, 130, 22, doi: 10.1007/s10569-017-9805-5
  * Ash et al. (1966) Ash, M. E., Smith, W. B., & Shapiro, I. I. 1966, AJ, 71, 844
  * Baer & Chesley (2008) Baer, J., & Chesley, S. R. 2008, Celestial Mechanics and Dynamical Astronomy, 100, 27, doi: 10.1007/s10569-007-9103-8
  * Bannister et al. (2018) Bannister, M. T., Gladman, B. J., Kavelaars, J. J., et al. 2018, ApJS, 236, 18, doi: 10.3847/1538-4365/aab77a
  * Batygin et al. (2019) Batygin, K., Adams, F. C., Brown, M. E., & Becker, J. C. 2019, Phys. Rep., 805, 1, doi: 10.1016/j.physrep.2019.01.009
  * Batygin & Brown (2016) Batygin, K., & Brown, M. E. 2016, AJ, 151, 22, doi: 10.3847/0004-6256/151/2/22
  * Baxter et al. (2018) Baxter, E. J., Jain, B., Blake, C., et al. 2018, arXiv e-prints, arXiv:1812.08701.  https://arxiv.org/abs/1812.08701
  * Beck & Giles (2005) Beck, J. G., & Giles, P. 2005, ApJ, 621, L153, doi: 10.1086/429224
  * Beichman (1988) Beichman, C. A. 1988, Astrophysical Letters and Communications, 27, 67
  * Bernardinelli et al. (2022a) Bernardinelli, P. H., Bernstein, G. M., Sako, M., et al. 2022a, ApJS, 258, 41, doi: 10.3847/1538-4365/ac3914
  * Bernardinelli et al. (2022b) —. 2022b, ApJS, 258, 41, doi: 10.3847/1538-4365/ac3914
  * Bernstein et al. (2017) Bernstein, G. M., Armstrong, R., Plazas, A. A., et al. 2017, PASP, 129, 074503, doi: 10.1088/1538-3873/aa6c55
  * Bottke et al. (2006) Bottke, William F., J., Vokrouhlický, D., Rubincam, D. P., & Nesvorný, D. 2006, Annual Review of Earth and Planetary Sciences, 34, 157, doi: 10.1146/annurev.earth.34.031405.125154
  * Brown & Batygin (2016) Brown, M. E., & Batygin, K. 2016, ApJ, 824, L23, doi: 10.3847/2041-8205/824/2/L23
  * Brown & Batygin (2021) —. 2021, AJ, 162, 219, doi: 10.3847/1538-3881/ac2056
  * Brown & Batygin (2022) —. 2022, AJ, 163, 102, doi: 10.3847/1538-3881/ac32dd
  * Brown et al. (2015) Brown, M. E., Bannister, M. T., Schmidt, B. P., et al. 2015, AJ, 149, 69, doi: 10.1088/0004-6256/149/2/69
  * Brunini & Melita (2002) Brunini, A., & Melita, M. D. 2002, Icarus, 160, 32, doi: 10.1006/icar.2002.6935
  * Buie & Folkner (2015) Buie, M. W., & Folkner, W. M. 2015, AJ, 149, 22, doi: 10.1088/0004-6256/149/1/22
  * Capitaine et al. (2003) Capitaine, N., Wallace, P. T., & Chapront, J. 2003, A&A, 412, 567, doi: 10.1051/0004-6361:20031539
  * Carry (2012) Carry, B. 2012, Planet. Space Sci., 73, 98, doi: 10.1016/j.pss.2012.03.009
  * Cramér (1946) Cramér, H. 1946, Mathematical Methods of Statistics, Goldstine Printed Materials (Princeton University Press).  https://books.google.com/books?id=_db1jwEACAAJ
  * Dallet (1901) Dallet, G. 1901, Bulletin de la Socieété astronomique de France, 16, 31
  * Duncan et al. (1988) Duncan, M., Quinn, T., & Tremaine, S. 1988, ApJ, 328, L69, doi: 10.1086/185162
  * Durante et al. (2020) Durante, D., Parisi, M., Serra, D., et al. 2020, Geophysical Research Letters, 47, e2019GL086572, doi: https://doi.org/10.1029/2019GL086572
  * Einstein (1915) Einstein, A. 1915, Sitzungsber. preuss.Akad. Wiss, 47, 831
  * Elliot et al. (2010) Elliot, J. L., Person, M. J., Zuluaga, C. A., et al. 2010, Nature, 465, 897
  * Everhart (1985) Everhart, E. 1985, in Astrophysics and Space Science Library, Vol. 115, IAU Colloq. 83: Dynamics of Comets: Their Origin and Evolution, ed. A. Carusi & G. B. Valsecchi, 185, doi: 10.1007/978-94-009-5400-7_17
  * Farnocchia et al. (2021) Farnocchia, D., Chesley, S. R., Takahashi, Y., et al. 2021, Icarus, 369, 114594, doi: 10.1016/j.icarus.2021.114594
  * Fernández et al. (2003) Fernández, Y. R., Sheppard, S. S., & Jewitt, D. C. 2003, AJ, 126, 1563, doi: 10.1086/377015
  * Fienga et al. (2016) Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2016, A&A, 587, L8, doi: 10.1051/0004-6361/201628227
  * Fienga et al. (2016) Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2016, arXiv e-prints, arXiv:1601.00947.  https://arxiv.org/abs/1601.00947
  * Fienga et al. (2020) Fienga, A., Ruscio, A. D., Bernus, L., et al. 2020, A&A, 640, A6, doi: 10.1051/0004-6361/202037919
  * Folkner et al. (2014) Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. 2014, Interplanetary Network Progress Report, 42-196, 1
  * Fortino et al. (2021a) Fortino, W. F., Bernstein, G. M., Bernardinelli, P. H., et al. 2021a, AJ, 162, 106, doi: 10.3847/1538-3881/ac0722
  * Fortino et al. (2021b) —. 2021b, AJ, 162, 106, doi: 10.3847/1538-3881/ac0722
  * Fraser et al. (2013) Fraser, W. C., Batygin, K., Brown, M. E., & Bouchez, A. 2013, Icarus, 222, 357, doi: https://doi.org/10.1016/j.icarus.2012.11.004
  * Gerdes et al. (2017) Gerdes, D. W., Sako, M., Hamilton, S., et al. 2017, The Astrophysical Journal Letters, 839, L15, doi: 10.3847/2041-8213/aa64d8
  * Gladman & Chan (2006) Gladman, B., & Chan, C. 2006, ApJ, 643, L135, doi: 10.1086/505214
  * Gladman et al. (2002) Gladman, B., Holman, M., Grav, T., et al. 2002, Icarus, 157, 269, doi: 10.1006/icar.2002.6860
  * Gladman & Volk (2021) Gladman, B., & Volk, K. 2021, ARA&A, 59, 203, doi: 10.1146/annurev-astro-120920-010005
  * Gladman et al. (2012) Gladman, B., Lawler, S. M., Petit, J.-M., et al. 2012, The Astronomical Journal, 144, 23, doi: 10.1088/0004-6256/144/1/23
  * Grundy et al. (2019a) Grundy, W., Noll, K., Buie, M., et al. 2019a, Icarus, 334, 30, doi: https://doi.org/10.1016/j.icarus.2018.12.037
  * Grundy et al. (2015) Grundy, W., Porter, S., Benecchi, S., et al. 2015, Icarus, 257, 130, doi: https://doi.org/10.1016/j.icarus.2015.04.036
  * Grundy et al. (2019b) Grundy, W., Noll, K., Roe, H., et al. 2019b, Icarus, 334, 62, doi: https://doi.org/10.1016/j.icarus.2019.03.035
  * Hogg et al. (1991) Hogg, D. W., Quinlan, G. D., & Tremaine, S. 1991, AJ, 101, 2274, doi: 10.1086/115849
  * Holler et al. (2021) Holler, B., Grundy, W., Buie, M., & Noll, K. 2021, Icarus, 355, 114130, doi: https://doi.org/10.1016/j.icarus.2020.114130
  * Holman & Payne (2016a) Holman, M. J., & Payne, M. J. 2016a, AJ, 152, 94, doi: 10.3847/0004-6256/152/4/94
  * Holman & Payne (2016b) —. 2016b, AJ, 152, 80, doi: 10.3847/0004-6256/152/4/80
  * Holman et al. (2019) Holman, M. J., Payne, M. J., & Pál, A. 2019, Research Notes of the American Astronomical Society, 3, 160, doi: 10.3847/2515-5172/ab4ea6
  * Holman & Wisdom (1993) Holman, M. J., & Wisdom, J. 1993, AJ, 105, 1987, doi: 10.1086/116574
  * Holman et al. (2018) Holman, M. J., Payne, M. J., Fraser, W., et al. 2018, ApJ, 855, L6, doi: 10.3847/2041-8213/aaadb3
  * Huang et al. (2022) Huang, Y., Gladman, B., Beaudoin, M., & Zhang, K. 2022, ApJ, 938, L23, doi: 10.3847/2041-8213/ac9480
  * Ivezić (2022) Ivezić, Z. 2022, Survey Cadence Optimization Committee’s Phase 1 recommendation, LSST Corporation.  https://pstn-053.lsst.io/
  * Jacobson (2009) Jacobson, R. A. 2009, The Astronomical Journal, 137, 4322, doi: 10.1088/0004-6256/137/5/4322
  * Jacobson (2014) —. 2014, The Astronomical Journal, 148, 76, doi: 10.1088/0004-6256/148/5/76
  * Jacobson et al. (2006) Jacobson, R. A., Antreasian, P. G., Bordi, J. J., et al. 2006, The Astronomical Journal, 132, 2520, doi: 10.1086/508812
  * Jewitt & Luu (1993) Jewitt, D., & Luu, J. 1993, Nature, 362, 730, doi: 10.1038/362730a0
  * Jones (2016) Jones, L. 2016, SMTN-002: Calculating LSST limiting magnitudes and SNR, LSST Corporation.  https://smtn-002.lsst.io/
  * Kavelaars et al. (2009) Kavelaars, J. J., Jones, R. L., Gladman, B. J., et al. 2009, The Astronomical Journal, 137, 4917, doi: 10.1088/0004-6256/137/6/4917
  * Kepler (1609) Kepler, J. 1609, Astronomia nova …, seu physica coelestis, tradita commentariis de motibus stellae martis, doi: 10.3931/e-rara-558
  * Kepler (1618) —. 1618, Epitome astronomiae copernicanae usitata forma quaestionum & responsionum conscripta, inque VII. libros digesta, doi: 10.3931/e-rara-3122
  * Kiss et al. (2019) Kiss, C., Marton, G., Parker, A. H., et al. 2019, Icarus, 334, 3, doi: https://doi.org/10.1016/j.icarus.2019.03.013
  * Konopliv et al. (1999) Konopliv, A., Banerdt, W., & Sjogren, W. 1999, Icarus, 139, 3, doi: https://doi.org/10.1006/icar.1999.6086
  * Konopliv et al. (2020) Konopliv, A., Park, R., & Ermakov, A. 2020, Icarus, 335, 113386, doi: https://doi.org/10.1016/j.icarus.2019.07.020
  * Konopliv et al. (2014) Konopliv, A., Asmar, S., Park, R., et al. 2014, Icarus, 240, 103, doi: https://doi.org/10.1016/j.icarus.2013.09.005
  * Konopliv et al. (2011) Konopliv, A. S., Asmar, S. W., Folkner, W. M., et al. 2011, Icarus, 211, 401, doi: https://doi.org/10.1016/j.icarus.2010.10.004
  * Konopliv et al. (2018) Konopliv, A. S., Park, R. S., Vaughan, A. T., et al. 2018, Icarus, 299, 411, doi: 10.1016/j.icarus.2017.08.005
  * Kowal et al. (1979) Kowal, C. T., Liller, W., & Marsden, B. G. 1979, in Dynamics of the Solar System, ed. R. L. Duncombe, Vol. 81, 245
  * Larsen et al. (2007) Larsen, J. A., Roe, E. S., Albert, C. E., et al. 2007, AJ, 133, 1247, doi: 10.1086/511155
  * Le Verrier (1846) Le Verrier, U. 1846, Comptes Rendus Hebdomadaires des Séances de l’Académie des Sciences, 22, 907
  * Le Verrier (1859) —. 1859, Comptes Rendus Hebdomadaires des Séances de l’Académie des Sciences, 49, 379
  * Levine et al. (2021) Levine, S. E., Zuluaga, C. A., Person, M. J., et al. 2021, AJ, 161, 210, doi: 10.3847/1538-3881/abe76d
  * Luhman (2014) Luhman, K. L. 2014, ApJ, 781, 4, doi: 10.1088/0004-637X/781/1/4
  * Lykawka & Mukai (2008) Lykawka, P. S., & Mukai, T. 2008, AJ, 135, 1161, doi: 10.1088/0004-6256/135/4/1161
  * Malhotra et al. (2016) Malhotra, R., Volk, K., & Wang, X. 2016, ApJ, 824, L22, doi: 10.3847/2041-8205/824/2/L22
  * Mallama (1993) Mallama, A. 1993, J. Geophys. Res., 98, 18873, doi: 10.1029/93JE01333
  * Mecheri et al. (2004) Mecheri, R., Abdelatif, T., Irbah, A., et al. 2004, Solar Physics, 222, 191, doi: 10.1023/B:SOLA.0000043563.96766.21
  * Meisner et al. (2017) Meisner, A. M., Bromley, B. C., Nugent, P. E., et al. 2017, AJ, 153, 65, doi: 10.3847/1538-3881/153/2/65
  * Melita et al. (2004) Melita, M. D., Williams, I. P., Collander-Brown, S. J., & Fitzsimmons, A. 2004, Icarus, 171, 516, doi: 10.1016/j.icarus.2004.05.014
  * Minor Planet Center (2021) Minor Planet Center. 2021, Data available from the Minor Planet Center.  https://minorplanetcenter.net/data
  * Muhleman et al. (1962) Muhleman, D. O., Holdridge, D. B., & Block, N. 1962, AJ, 67, 191, doi: 10.1086/108693
  * Naess et al. (2021) Naess, S., Aiola, S., Battaglia, N., et al. 2021, ApJ, 923, 224, doi: 10.3847/1538-4357/ac2307
  * Newton (1687) Newton, I. 1687, Philosophiae Naturalis Principia Mathematica. Auctore Js. Newton, doi: 10.3931/e-rara-440
  * Oort (1950) Oort, J. H. 1950, Bull. Astron. Inst. Netherlands, 11, 91
  * Park et al. (2021) Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. 2021, The Astronomical Journal, 161, 105, doi: 10.3847/1538-3881/abd414
  * Petit et al. (2011) Petit, J.-M., Kavelaars, J. J., Gladman, B. J., et al. 2011, The Astronomical Journal, 142, 131, doi: 10.1088/0004-6256/142/4/131
  * Pijpers (1998) Pijpers, F. P. 1998, MNRAS, 297, L76, doi: 10.1046/j.1365-8711.1998.01801.x
  * Pitjeva & Pitjev (2018) Pitjeva, E. V., & Pitjev, N. P. 2018, Celestial Mechanics and Dynamical Astronomy, 130, 57, doi: 10.1007/s10569-018-9853-5
  * Pitjeva & Standish (2009) Pitjeva, E. V., & Standish, E. M. 2009, Celestial Mechanics and Dynamical Astronomy, 103, 365, doi: 10.1007/s10569-009-9203-8
  * Pál et al. (2012) Pál, A., Kiss, C., Müller, T. G., et al. 2012, A&A, 541, L6, doi: 10.1051/0004-6361/201218874
  * Ragozzine & Brown (2009) Ragozzine, D., & Brown, M. E. 2009, AJ, 137, 4766, doi: 10.1088/0004-6256/137/6/4766
  * Rao & Das Gupta (1989) Rao, C., & Das Gupta, S. 1989, Selected Papers of C.R. Rao, Selected Papers of C.R. Rao No. v. 5 (Indian Statistical Institute).  https://books.google.com/books?id=0AlL7uE2H9YC
  * Rein & Spiegel (2015) Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424, doi: 10.1093/mnras/stu2164
  * Rein & Tamayo (2016) Rein, H., & Tamayo, D. 2016, MNRAS, 459, 2275, doi: 10.1093/mnras/stw644
  * Rein, H. & Liu, S.-F. (2012) Rein, H., & Liu, S.-F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rice & Laughlin (2019) Rice, M., & Laughlin, G. 2019, AJ, 158, 19, doi: 10.3847/1538-3881/ab21df
  * Rice & Laughlin (2020) —. 2020, , 1, 81, doi: 10.3847/PSJ/abc42c
  * Shapiro et al. (1968) Shapiro, I. I., Pettengill, G. H., Ash, M. E., et al. 1968, Phys. Rev. Lett., 21, 266, doi: 10.1103/PhysRevLett.21.266.3
  * Sheppard et al. (2018) Sheppard, S. S., Fernandez, Y. R., & Moullet, A. 2018, The Astronomical Journal, 156, 270, doi: 10.3847/1538-3881/aae92a
  * Sheppard & Trujillo (2016) Sheppard, S. S., & Trujillo, C. 2016, AJ, 152, 221, doi: 10.3847/1538-3881/152/6/221
  * Silsbee & Tremaine (2018) Silsbee, K., & Tremaine, S. 2018, AJ, 155, 75, doi: 10.3847/1538-3881/aaa19b
  * Smith et al. (1966) Smith, W. B., Shapiro, I. I., & Ash, M. E. 1966, AJ, 71, 871
  * Standish (1993) Standish, E. M. 1993, AJ, 105, 2000, doi: 10.1086/116575
  * Stern et al. (2015) Stern, S. A., Bagenal, F., Ennico, K., et al. 2015, Science, 350, aad1815, doi: 10.1126/science.aad1815
  * Tamayo et al. (2019) Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. 2019, Monthly Notices of the Royal Astronomical Society, 491, 2885, doi: 10.1093/mnras/stz2870
  * Tanga et al. (2022) Tanga, P., Pauwels, T., Mignard, F., et al. 2022, arXiv e-prints, arXiv:2206.05561.  https://arxiv.org/abs/2206.05561
  * Tombaugh (1946) Tombaugh, C. W. 1946, Leaflet of the Astronomical Society of the Pacific, 5, 73
  * Tombaugh (1996) Tombaugh, C. W. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 107, Completing the Inventory of the Solar System, ed. T. Rettig & J. M. Hahn, 157–162
  * Tremaine (1990) Tremaine, S. 1990, in NATO Advanced Study Institute (ASI) Series C, Vol. 306, Neutron Stars and Their Birth Events, ed. D. Lynden-Bell & G. Gilmore, 37
  * Trujillo & Sheppard (2014) Trujillo, C. A., & Sheppard, S. S. 2014, Nature, 507, 471, doi: 10.1038/nature13156
  * Vilenius et al. (2014) Vilenius, E., Kiss, C., Müller, T., et al. 2014, A&A, 564, A35, doi: 10.1051/0004-6361/201322416
  * Vokrouhlický & Farinella (1998) Vokrouhlický, D., & Farinella, P. 1998, AJ, 116, 2032, doi: 10.1086/300565