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

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 possible sky positions) of $\sigma_M=0.22M_\oplus (d_x/400\,\textrm{AU})^3.$ A definitive $(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 $\approx8000$ 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 useful cross-checks with Juno/Cassini/Mars ranging, but do not significantly improve the best-achievable $\sigma_M$ values until they are $\gtrsim10\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 Planet X uncertainty until the data get substantially better than they are today.


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(Kepler , 1618 and general motion (Newton 1687), the successful prediction of Nep-dchgomes@gmail.com arXiv:2212.09594v2 [astro-ph.EP] 11 May 2023 tune'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.

Searching for Planet X
Although the prediction and discovery of Neptune motivated searches for additional planets (see Tremaine 1990 andBatygin 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(Tombaugh , 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;. 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 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,  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  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 ⊕ , carrying out full ephemeris solutions to all the available data, and comparing the residuals to those without the addition of Planet X,  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  investigation is based on an earlier version of the INPOP ephemeris model that did not include mass in the Kuiper belt aside from Pluto , and the Holman & Payne (2016a) study is built upon the residuals shown in  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 ∼ 25 • from the ecliptic plane, closer than 500 AU if it is 5 M ⊕ or closer than 650 AU if it is 10 M ⊕ .
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  or . 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 ⊕ .
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.

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 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 tidal P 2 ∝ d 4 T (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 ≈ 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 ≈ 20 µ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 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 1 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 ≈ 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 σ min ≈ 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 • nearly 1000 times in its 10-year survey, including for example ∼ 8000 Jovian Trojan asteroids at H < 14.5. As detailed in Section 6, the collective accuracy on a shift of this entire population is ≈ 2.5µ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× closer, and O(10 5 ) of them will be bright enough to approach σ 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 ≥ 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.

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 σ 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 σ 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.

THE BACK OF THE ENVELOPE
Before solving numerically for the sensitivity σ 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 ≈ 6M ⊕ 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 ⊕ , 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 φ X and polar angle θ X . It will also be useful to define the radial, azimuthal, and vertical directions relative to the ecliptic asρ,φ, and 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.

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 . The = 0 and = 1 terms are a potential shift and a constant acceleration, respectively, which are unobservable with differential measurements within the solar system. The = 2 tidal terms are a sufficient description of the observable effects of Planet X, since ≥ 3 terms are suppressed relative to this by factors of d T /d X 0.01. We can thus rephrase our search for Planet X as a search for unexplained = 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 × 3 tidal matrix, or as the real (for m ≥ 0) and imaginary (for m > 0) coefficients A 2m of the tidal potential A single point source like Planet X has only three degrees of freedom (e.g. M X , θ X , φ 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 (θ X , φ X ) can be expressed as a linear superposition of these, and we can infer a bound on M X /d 3 X 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 3 X , we will quantify all our results as the measurement accuracy σ M attained on M X at fixed nominal distance d X = 400 AU. Then we can rescale σ 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 ⊕ at d X = 800 AU at 5σ significance, which is equivalent to obtaining σ M < 0.125M ⊕ 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: 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 (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 .

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 Ω 0 . For a nominal Planet X at an ecliptic pole, this is 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 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 P/2. This would lead to an apparent angular lag that grows linearly with time, equal to ∆φ = 2π 2 If the distance for a body is determined to sub-meter precision by ranging, then the period shift would eventually become detectable with µ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.

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ẑ epicycle, largely undetectable by ranging) and longitude of perihelion. The latter is accessible to ranging viaρ epicycles, and also to astrometry viaφ 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 Ω 0 , which leads to perturbations whose amplitudes grow linearly over multiple orbital cycles, leading to transverse astrometric displacements (angular shifts in any direction) of The signal accessible to ranging is primarily theρ component of apsidal precession 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 ∼ 1200 Trojans with typical projected LSST astrometric accuracy below 5 mas per component per observation, and 1.6 × 10 5 MBAs that would satisfy the same criterium. The 130× 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 ≈ 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× 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.

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 ≈ a tidal cos mφ 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 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 > 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.

MATHEMATICAL METHODS AND MODEL
We estimate the ability of a set of observational data points 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 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.

The Fisher matrix
For a system with likelihood L(d|p) of obtaining the data given the parameters, the available information on our parameter set is expressed by the N × N Fisher information matrix F: Let us consider an unbiased estimator θ = {θ 1 , θ 2 , ..., θ N } of our parameter set p. The multivariate Cramér-Rao inequality imposes the following restriction on the covariance matrix of any such estimator: where the matrix inequality states that the difference matrix Cov . This condition allows us to establish a lower bound for the variance of any individual estimator θ i of the parameter p i . If the data points d have independent gaussian errors σ d = {σ 1 , σ 2 , ..., σ k } and means µ d = {µ 1 (p), µ 2 (p), ..., µ N (p)}, the Fisher information becomes which simplifies to The lower bound for the unmarginalized variance of θ i is The process of marginalizing over all parameters p j for which j = 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 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/σ 2 M to the appropriate diagonal element of 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 p R and nuisance parameters subset 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 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.

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 3 X , 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 × 3 tensor, or by the Hermitian coefficients of the = 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 , . . . , 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 × 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.
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 . 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 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.
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 σ 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 σ 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.06AU and 3.27AU. Each ring's mass is a free parameter, to which Gaussian priors with σ = 3.3 × 10 −5 M ⊕ 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 AU and at 43.0 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 ⊕ , 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 σ 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.

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 × 10 4 years. We sample the five components of the internal quadrupole field in 50 year intervals, and we measure the 5×5 covariance matrix of the quadrupole components from these samples. Variations are of the order of 10 −3 M CFEPS /(40 AU) 3 , where M CFEPS is the total mass of the population, estimated at 1.5 × 10 −2 M ⊕ . 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.

Mass uncertainties from large TNOs
The CFEPS sample reaches magnitudes as bright as H ≈ 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 σ 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.

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. 1 Beginning with this base simulation, we compute variational derivatives of planetary positions with respect to their (8 × 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 ⊕ and circular orbit at d X = 400 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); 2 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 ≤ 10 −8 , so this small inconsistency can be neglected for the purposes of making Fisher-matrix derivatives.  (2009)  b We 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. c Fraser et al. (2013) gives four possible orbit solutions for the Quaoar system. We assume this ambiguity will be solved during the following years. d No dynamical mass measurement. Mass uncertainty assumed to be half the total mass estimate given its radius and ρ ≈ 2 g/cm 3 . Reference is for object radius.
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 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 · 10 −7 (e.g. Pijpers 1998). The Sun's north pole is tilted from the north ecliptic pole by 16.12 • in right ascension and −2.70 • 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ẑ axis aligns with the solar pole, the force from a J 2 can be written as where x p ,y p ,z p are the cartesian distances from the sun along each axis, and r p = x 2 p + y 2 p + z 2 p . 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 • (Beck & Giles 2005).
For a given observer our Cartesian positions can be converted to an observable position on the sky through α = arctan(y(t)/x(t)), δ = arcsin(z(t)/ x(t) 2 + y(t) 2 + z(t) 2 ), where α and δ 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 − τ ) − x obs (t), for some light-time delay τ ). 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.     c N: number of data points. For the astrometric data with different errors on ra and dec, uncertainty listed is the largest one.

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.

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 website 3 , plus the complete set of Cassini data points available on the INPOP Astrometric Planetary Database website 4 .
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 KBO floor c · · · · · · · · · 0.053 0.066 0.071 100.0% 100.0% a Fraction of the sky over which a 5M ⊕ Planet X would be a ≥ 5σ detection.
b Ranging data from Juno and Cassini.
c Limits purely from degeneracies with uncertainties in Kuiper Belt quadrupole. 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.

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 σ 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 σ M across the sky. We also tabulate the percentage of the sky f 400 for which σ M < 1, meaning that the presence of a 5M ⊕ Planet X at a distance of 400 AU would be detected unambiguously, at 5σ 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σ detection of 5M ⊕ at a distance d X = 800 AU (σ M < 0.125M ⊕ on our plots) is satisfied over f 800 = 4.8% of the sky in the nominal case. This is shown for the nominal scenario (solid), one restricted to already-published ranging, and another extrapolating ranging data through 2034. 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→current→predicted data adds 5-15% of sky coverage attaining any particular σ M .

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 σ M constraints as we inflate the measurement errors on either the Juno/Cassini ranging or on the Mars ranging. The σ 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.

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 σ 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 σ 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 ⊕ , d X = 800 AU model across a substantial portion of the sky. The middle panel of Figure 1 shows the sky map of σ M in the Old+J/C case, illustrating the weaker (but still interesting) level of constraint compared to the upper, nominal case.

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.

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.  . Each curve plots the fraction of the sky (y axis) for which the uncertainty σ 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 σ M values that are needed for a definitive 5σ detection of Planet X if it has mass of 5M ⊕ 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. 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 σ 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 σ psf of the point-spread function (PSF), and the filter band Q ∈ {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 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 1 AU d ≈ 11 • for Jovian Trojans, and will not significantly impact the overall M X sensitivity. The photon count is then given by The best achievable astrometric uncertainty for an unresolved source is derived using Fisher matrix algebra in Appendix B. Reproducing Equation (B16), a circular Gaussian PSF with dispersion σ PSF on each axis yields an astrometric uncertainty on each axis of: . Experience with the Dark Energy Survey (DES) suggests that PSF-fitting astrometry can closely approach these theoretical limits.
The uncalibrated portion σ min of the atmospheric turbulence displacements then must be added in quadrature to the σ x value of Equation (22) to yield the final σ 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 σ min ≈ 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.

Trojan catalog
The list of all known Jovian Trojans was downloaded from the Minor Planet Center database 5 (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 ≤ 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.

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 ∈ {g, r, i, z} is ∼ 1.2 × 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 12year survey, with the times of each and the per-axis measurement error σ 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.

Information vs depth
A useful measure of the constraining power of observations of some minor body population is the collective uncertainty σ pop defined through this sum over the astrometric uncertainties σ ij per axis in observations j of bodies i with absolute magnitude H i brighter than some threshold: This quantity gives the RMS size of a perturbation that can be detected (per axis) at 1σ 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 ≈ 3 µ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.
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 σ min = 2 mas subtends ≈ 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 σ ij value. At ground-based seeing of 0. 5-1. 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 σ 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.

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 σ pop level, i.e. sub-µas and sub-meter scale. We next discuss the systematic errors that we expect to encounter, their expected sizes, and mitigation strategies.

Radiation pressure
If the ratio of incident solar radiation pressure to the solar gravitational force on a spherical target of radius R T , density ρ, and heliocentric distance d T is well below , then we can ignore its influence. This condition is 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 signalbut 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 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.

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 ≈ 200 K , then the total acceleration from outgassing or mass loss is bounded by 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× 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.

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 ∆x(φ 0 ) that is purely a function of the Sun-Target-Observer phase angle φ 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 φ 0 is bound to ± sin −1 1 AU d T . In this range we will approximate ∆x(φ 0 ) = kφ 0 R T 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.
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 ≈ 0.4. These values suggest that the angular PCM for even the smallest Trojan targets will be ≈ 4 km × 0.4/5.2 ≈ 300 m and manifest as 1.5 km under-estimates of inferred radial coordinates. This is several orders of magnitude above both σ pop and the nominal Planet X signatures. Indeed our expected σ 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 ≈ 1 mas level from the rotation of (21) Lutetia. PCM values this high would inflate our expected statistical error floor of ≈ 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 Sun kR ϕ kϕR kR R true apparent 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 φ is kφ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 stochastic rotational portion of PCM, and concern ourselves only with the rotation-averaged values.

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 ∆z = k bj c bi tan z i 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 × larger than the target σ 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 σ min . The much greater challenge is in insuring that coherent residuals to the DCR corrections do not impact our inferences being made at µ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.

Timing accuracy
The maximum apparent motion of a typical Trojan is ≈ 10 mas/s. To keep ephemeris errors below σ pop ≈ 2 µas, any systematic errors in the temporal midpoint of the LSST exposures should be 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. 6 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.

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.

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.

Photocenter motion
The photocenter motion displacement is computed from Equation (28), with its direction given by the phase position angle φ P A . This gives us The kR T term is our free parameter, and derivatives with respect to it are included in our Fisher matrix modeling for each Trojan.

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 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.

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.
The exact proportionality depends on the albedo and shape of the object. Assuming a constant density for all Trojans, the proportionality becomes a R ∝ c RP /(a 2 R T ), where c 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 RP are proportional to derivatives with respect to M . We obtain the former by rescaling the latter according to the radius R T of each Trojan, which is inferred from its magnitude H.

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 ≈ 50 meters) the thermal acceleration can be computed using Equation (16) in Vokrouhlický & Farinella (1998) where α is the albedo, a is the semi major axis in AU, and A( ) 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θ acceleration with magnitude of 10 −12 m s −2 , a rough value for a body with a radius of 50 m and a density of ρ = 3500kg 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).

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 kth 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 F is broken into blocks among these parameter subgroups, all of the blocks representing covariances between p k and p k are zero for k = k . This means that Equation (17) can be successively applied to all Trojans to marginalize down to the global parameters: The marginalization over all Trojans' parameters never requires inversion of matrices larger than 11 × 11 and there is never need to store the full Fisher matrix.
7. RESULTS WITH TROJAN TRACKING

Combinations with ranging data
We now calculate the full-sky Fisher uncertainties σ 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 σ 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 σ M from 1.07 to 0.98M ⊕ . 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 σ M on the sky from 12.4 to 2.2M ⊕ . Thus a 10M ⊕ 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σ detection of M X = 5M ⊕ 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.

Effect of σ min
A reduction on the turbulence floor σ 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 σ M near the ecliptic poles is decreased from 0.98 to 0.89M ⊕ , 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.

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 σ 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.

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 σ M , perhaps because it only generates the m = 0 spherical harmonic, and any non-polar Planet X will have components at m = ±1 and m = ±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 ⊕ ). Since the tidal field scales as M X /d 3 X , and the KBOs are at d ≥ 40 AU while we are calibrating σ M for d X = 400 AU, we expect the ≈ 10 −5 M ⊕ uncertainties in the masses of the largest moon-less KBOs to result in an ≈ 10 −2 M ⊕ floor on σ 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 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 σ m ≈ 0.06M ⊕ 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 ⊕ Planet X at d X = 800 AU anywhere on the sky. . The fraction of the sky with uncertainty σ 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 ⊕ Planet X at 5σ 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.
What observations would be needed to bring the achievable σ 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 0.1 those achievable with LSST. At 20× reduced errors, we can confidently detect the most pessimistic version of Planet X, 5M 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× 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) µas for occultations of brighter Gaia stars. We leave a detailed quantitative study of this possibility for future work.

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σ detection of a 5M ⊕ 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 ∆χ 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 ⊕ placed at 400 < d X < 800 AU, at ecliptic latitudes between ±25 • . 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 3 X 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 ⊕ Planet X at 600 AU is near a 3σ detection in most of their footprint, equivalent to σ M ≈ 0.5M in our reference case of d X = 400 AU . A similar conclusion is drawn from their depiction of ≈ 1σ signals from M X = 5M ⊕ at 800 AU. This constraint is 2 times weaker than the median σ M ≥ 0.24M ⊕ 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 W 1 = 16.66 from a search of the WISE spacecraft images. Using the models of , this W 1 limit corresponds roughly to excluding M X ≥ 20M ⊕ at 622 AU. Our forecasted gravitational bounds on M X are tighter than this directimaging 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× 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 ≥ 400 AU away, but the tidal-field constraints are equally applicable to smaller or multiple bodies that are 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 σ M < 0.7M 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σ significance for d X ≤ 260 AU, and a Mars-mass object for d X ≤ 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 ∝ d −3 X , is more favorable than for searches for reflected solar light, which is ∝ d −4 X . The 5 × 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 vector a(θ , φ ) for a given direction in the sky is found by solving Equation (A9), and can be used to contract the 5-dimensional Fisher matrix, This allows us to find the Fisher information on the mass of Planet X at any position using derivatives computed solely for five positions.
B. OPTIMAL ASTROMETRIC ACCURACY When fitting a model of a point-spread function (PSF) to an unresolved image of an asteroid, there are three free parameters: flux f , and the centroids x 0 , y 0 . The expected signal f k in pixel k at location (x k , y k ) subtending solid angle ∆x∆y is where P is the PSF function, we assume that the sky background b is in units of photocarriers per unit solid angle, the flux f is in units of photocarriers. If the measurement f k is distributed as f k ∼ Poisson(f k ), then it is straightforward to show that the Fisher matrix in Equation (11) for f k is and the full Fisher matrix is the sum over the statistically independent pixels k. For a circularly symmetric PSF P (r), we can combine the above 2 equations and use the parameter vector {f, x 0 , y 0 } to obtain F f f = k ∆x ∆y P 2 (r k ) b + f P (r k ) → 2π r dr P 2 (r) b + f P (r) (B13) F xx = F yy = k ∆x ∆y f 2 b + f P (r k ) ∂P (r) ∂r 2 k cos 2 θ k → πf 2 r dr 1 b + f P (r) ∂P (r) ∂r 2 (B14) with the off-diagonal elements being zero. The uncertainty σ x in the x position is equal to that in the y direction and is the inverse square root of F xx . If we take the PSF to be a circular Gaussian then we can change Equation (B14) to be a dimensionless integral The parameter ν is the ratio of background flux to source flux.