Issue 
A&A
Volume 640, August 2020



Article Number  A6  
Number of page(s)  11  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/202037919  
Published online  28 July 2020 
New constraints on the location of P9 obtained with the INPOP19a planetary ephemeris
^{1}
GéoAzur, CNRS, Observatoire de la Côte d’Azur, Université Côte d’Azur, 250 Av. A. Einstein, Valbonne 06560, France
email: agnes.fienga@geoazur.unice.fr
^{2}
IMCCE, Observatoire de Paris, PSL University, CNRS, Sorbonne Université, 77 Avenue DenfertRochereau, Paris 75014, France
^{3}
Dipartimento di Ingegneria Meccanica e Aerospaziale, Sapienza Universitä di Roma, Via Eudossiana 18, Rome 00184, Italy
Received:
10
March
2020
Accepted:
27
May
2020
Context. We used the new released INPOP19a planetary ephemerides benefiting from Jupiterupdated positions by the Juno mission and reanalyzed Cassini observations.
Aims. We test possible locations of the unknown planet P9. To do this, we used the perturbations it produces on the orbits of the outer planets, more specifically, on the orbit of Saturn.
Methods. Two statistical criteria were used to identify possible acceptable locations of P9 according to (i) the difference in planetary positions when P9 is included compared with the propagated covariance matrix, and (ii) the χ^{2} likelihood of postfit residuals for ephemerides when P9 is included.
Results. No significant improvement of the residuals was found for any of the simulated locations, but we provide zones that induce a significant degradation of the ephemerides.
Conclusions. Based on the INPOP19a planetary ephemerides, we demonstrate that if P9 exists, it cannot be closer than 500 AU with a 5 M_{⊕} and no closer than 650 AU with a 10 M_{⊕}. We also show that there is no clear zone that would indicate the positive existence of planet P9, but there are zones for which the existence of P9 is compatible with the 3σ accuracy of the INPOP planetary ephemerides.
Key words: celestial mechanics / ephemerides / Kuiper belt: general / planets and satellites: detection
© A. Fienga et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
In 2016, Batygin & Brown (2016a) proposed the existence of an unseen massive planet, called P9, to explain the spatial distribution of about 20 Kuiper Belt objects (KBO) that are located away from the influence of Neptune. In a more recent analysis, Batygin et al. (2019) restricted the list of objects showing significant perihelion asymmetry to 11 of the 14 known transNeptunian objects (TNOs) with q > 30 AU, a > 250 AU and i < 40°. Kaib et al. (2019) argued that considering the stability of their orbits, this number can be limited to 9 objects.
Several discussion points arose from Batygin & Brown (2016a). Firstly, the real significance of the orbital clustering for such a limited number of objects is still discussed (Shankman et al. 2017a; Lawler et al. 2017; Brown 2017). Secondly, when the KBO orbit clustering is admitted to exist, a possible planet 9 is not the only viable explanation (Kaib et al. 2019; Madigan & McCourt 2016; Gomes et al. 2006). Finally, other evidence of the existence of P9 has been proposed with more or less positive outcomes (Kaib et al. 2019; Shankman et al. 2017b; Bailey et al. 2016; Batygin & Brown 2016b; Gomes et al. 2017; Millholland & Laughlin 2017; Nesvorný et al. 2017). Several tests of the existence of P9 and its possible location have been proposed. Fienga et al. (2016) have proposed to use planetary ephemerides. The authors showed that a planet of the size of 10 Earth masses (M_{⊕}) at 600 AU and on the orbit required by Batygin & Brown (2016a) would affect the planetary orbits. The accuracy with which the position of Saturn can be deduced from the Cassini tracking data enabled an exclusion zone for P9 to be proposed. From this preliminary work, Holman & Payne (2016) further developed this approach and produced maps of favorable locations for P9. In 2018, Pitjeva & Pitjev (2018) proposed to consider the dynamical effect of the other known KBO objects as potential perturbers of the planetary orbits, which might mask the perturbations of planet 9. More recently, Batygin et al. (2019) reviewed the topic and discussed different evidence that P9 might exist. The authors updated their modeling with new characteristics for the orbit of P9. These latest constraints on the orbit are larger that those that were provided by Batygin & Brown (2016a): the mass can range from 5 to 10 M_{⊕}, the inclination from 15° to 25°, the semimajor axis from 400 to 800 AU, and the eccentricity from 0.2 to 0.5. This corresponds to a distance between P9 and the Solar System barycenter (SSB) from 200 to 640 AU. In 2019, a new version of the INPOP planetary ephemerides has been released (Fienga et al. 2019, 2020). This version benefits from the use of nine very accurate Jupiternormal points deduced from the Juno mission and from a new analysis of the Cassini data (Fienga et al. 2019). INPOP19a also shows great improvement compared to the previous INPOP versions because it includes dynamical modeling of the perturbations that are induced by a ring of TNOs, as proposed by Pitjeva & Pitjev (2018). In this context, INPOP19a appears to be a good tool for testing the hypothesis that P9 exists. We here first recall some important elements related to INPOP19a (Sect. 1). We then implement the P9 perturbations in the planetary orbits (Sect. 3.1) and define acceptable zones for the existence of P9 using two criteria: we match the propagated covariance to INPOP19a (Sect. 4.1), and we determine the likelihood of the solution relative to INPOP19a (Sect. 3.3). In Sect. 4 we give the results that we obtained with 3150 simulations considering P9 distances from 400 AU to 800 AU, two masses (5 and 10 M_{⊕}), and using the two criteria described in Sects. 4.1 and 3.3. We compare and discuss the results in Sect. 5, and we conclude in Sect. 6.
2. INPOP19a planetary ephemerides
2.1. Data sample
The full data sets we used to adjust INPOP19a are presented in Table 1. A detailed description of INPOP19a can be found in Fienga et al. (2019) and in Fienga et al. (2020). Two main inputs were significant for INPOP19a: the reanalysis of some periods of the Cassini navigation data for the orbit of Saturn, including points from the Grand Finale and Titan flyby gravity solutions, and the addition of Juno measurements for Jupiter.
INPOP19a data samples we used to adjust P9.
2.1.1. Independent analysis of Cassini data
The NASA Cassini mission gathered an important amount of scientific data during the 13 year duration of the mission around Saturn. In particular, radio science and navigation measurements represented a unique tool for constraining the position of the Saturn barycenter relative to the Earth. Cassini normal points have been produced from a reanalysis of navigation data for the periods 2006, 2008, 2009, and 2011 (Di Ruscio et al. 2020). The new data analysis relies on the updated knowledge of the Saturnian system as acquired throughout the mission: the better accuracies achieved for the Saturn moon ephemerides, and the last gravity solutions of Saturn and its main satellites produced by the radio science team (Durante et al. 2019, 2020; Iess et al. 2014, 2019). For each arc, we solved for the spacecraft initial position and velocity, corrections for orbital trim and reaction wheel desaturation maneuvers, and RTGinduced anisotropic acceleration. In addition, stochastic accelerations at the level of 5 × 10^{−13} km s^{−2} (updated every 8 h) were included to compensate for any remaining dynamical mismodeling; this is the same approach that the Cassini navigation team followed during the mission. We estimated a correction for range measurements in the form of stochastic biases individually for each tracking pass, with a large a priori uncertainty to account for both station calibration and ephemerides error. The reconstructed Cassini trajectories were thus used to produce normal points, including the estimated range biases in computing the light time for the groundstation to Cassini roundtrip. The uncertainty on the normal points is given by the estimated covariance matrix of the range biases.
We also added additional normal points deduced from the radio science solutions for the gravitydedicated flybys of Titan (Durante et al. 2019) and Grand Finale Saturn pericenters (Iess et al. 2019). These range normal points were obtained by considering a passthrough of range data on the spacecraft orbit reconstructions that were produced for the gravity field solutions using only Doppler data. With these supplementary normal points, the period covered by Cassini was extended up to the end of the mission in 2017, when the spacecraft plunged into the atmosphere of Saturn. In Table 1, the newly analyzed navigation data and the normal points deduced from Titan gravity flybys (TGF) are labeled Nav. + TGF range, and the data set deduced from the Grand Finale is labeled Grand Finale range (see Table 1).
2.1.2. Nine perijoves of Juno
The Juno spacecraft has been orbiting Jupiter on a highly eccentric polar orbit since July 2016. A radioscience experiment aims at characterizing the gravity field of the gas giant to unprecedented accuracy (Durante et al. 2020; Folkner et al. 2017; Iess et al. 2018). The extremely accurate radiotracking system of Juno enables simultaneous twoway Doppler measurements at X and Ka band with accuracies as low as 10 μs^{−1} in the radial velocity during the gravitydedicated passes, which are used to reconstruct the spacecraft trajectory with a meterlevel radial accuracy with respect to the center of mass of Jupiter at perijove (periapsis in orbit around Jupiter). Range data points at X band are collected as well, and the Jovian barycenter positions relative to the Earth can be generated once per perijove pass, provided that we know the position of Juno with respect to the Jovian barycenter. In our fit, we include a total of nine new Jupiternormal points spanning the period from the orbital insertion, back in 2016, to the end of 2018.
2.2. TransNeptunian objects
As described in Di Ruscio et al. (2020) and Fienga et al. (2019), the extension of the Saturn data sets with the positions deduced from the Grand Finale and from the analysis reported by Di Ruscio et al. (2020) required introducing the perturbations induced by the TNOs on the planetary orbits into the INPOP dynamical modeling. These perturbations were introduced first by taking into account individual accelerations produced by the nine most massive binary TNOs, and second, by modeling the remaining TNO perturbations by a set of three rings centered on the SSB, which is located in the ecliptic plane and at distances of 39.4, 44, and 47.5 AU. The sum of the masses of these three rings is estimated by the fit during the construction of INPOP19a. The rings allowed Di Ruscio et al. (2020) to show that the residuals obtained for the Grand Finale are improved by a factor 30 between INPOP17a and INPOP19a and that the mass for the TNO ring is equal to (0.061 ± 0.001) M_{⊕}. This value is higher than the value obtained by Pitjeva & Pitjev (2018). This can be explained by the differences in the dynamical modeling (in addition to the ring, nine masses have been fixed in Di Ruscio et al. 2020, while 31 were included in the fit in Pitjeva & Pitjev 2018) and by the differences in the data that were used (Pitjeva & Pitjev 2018 included neither the Cassini navigation and Grand Finale data nor the Juno observations). When the same data sets and the sum of fixed masses are considered, the two results are consistent at 3σ. For more details, we refer to Di Ruscio et al. (2020) and Fienga et al. (2019).
3. P9 detection method
3.1. Modeling
In order to characterize the regions in which P9 might be located, we built different planetary ephemerides that considered different positions of P9 following the method initially developed by Fienga et al. (2016). However, because the expected semimajor axis of P9 lies between 400 and 800 AU, its orbital circular period (from 8000 years to more than 22 500 years) almost fixes the location of P9 relative to the SSB over the time span of the planetary observations (100 years for the entire data set, and 13 years for the Cassini data sample). We therefore considered here that the dynamical effect of P9 on the planetary orbits can be modeled as a tidal effect depending on its fixed position, that is, right ascension (RA) and declination (Dec), it’s distance to the SSB (r), and its mass (M_{P9}). The acceleration induced by P9 on planet A is such that
where u is the unit vector of coordinates defined by P9 (RA, Dec), x_{A} is the barycentric position vector of planet A, and the dot indicates the scalar product. We then constructed planetary ephemerides by numerically integrating and fitting the orbits to the INPOP19a data sample presented in Table 1, for which we used the same weighting scheme as for INPOP19a (see, e.g., Fienga et al. 2019 for more details). For each r, M_{P9}, RA, and Dec, we obtained a fitted solution. In order to explore more possibilities than those proposed by Batygin et al. (2019), we considered distances to the SSB from 400 AU to 800 AU with two cases: M_{P9} = 5 Earth masses (M_{⊕}), and M_{P9} = 10 M_{⊕} for the farther positions. We made a first RA, Dec grid with a 10° sampling in the ecliptic plane and then added supplementary runs on some specific zones out of the ecliptic plane (see Sect. 5).
After the solutions were iteratively fitted, statistical criteria were used to determine which solution was acceptable with respect to INPOP19a and which should be rejected. On this selection depends the definition of an acceptable zone for P9. We used two complementary criteria: one considering the covariance of each solution propagated over the time coverage of the data sample, and one based on the computation of the χ^{2} likelihood.
3.2. Propagation of the covariance
INPOP was computed by numerically solving the equations of motion. Let X(t) be the state vector in barycentric coordinates that contains positions and velocities of each body whose trajectory is computed. The numerical integrator solves a CauchyLipschitz equationofmotion system,
where X ∈ ℝ^{6n} , with n being the number of integrated bodies, and P ∈ ℝ^{p} with n ≤ p contains all constant parameters of the ephemeris (initial conditions for the planetary orbits; masses of the Sun and of the asteroids, including transNeptunian objects; oblateness of the Sun; and the EarthMoon mass ratio). We note that X and P are not independent variables because P includes the initial condition X_{0}. Modifications of the parameter P modify the trajectories X(t). From this ephemeris, we computed observational simulations in order to compare them to real data. C(t_{i}, P) is the observable quantity at date t_{i} computed with parameters P (we considered that the dependence on X(t_{i}) is included in the dependence on the initial conditions included in P, which are integrated by INPOP). The goal of the ephemeris is to minimize the norm of the residual vector
where 𝒪(t_{i}) is the real observation at date t_{i} (for any matrix A, transposition of matrix A is noted ^{t}A). Usually, and this is what we did here, a linear approximation and Gaussian distribution of the observational noise are assumed. It is then well known that the parameters P minimize χ^{2} = ^{t}RWR, where W is the weighting matrix representing the accuracy of the observational data, and they are given by the algorithm that increments P by iterations by adding
until convergence is reached. Here, R is the residual vector as defined in Eq. (3). J is the Jacobian matrix of R(t_{i}, P) or C(t_{i}, P).
The covariance of P, which represents its uncertainty if the linear approximation and the Gaussian distribution of error are realized, is then
From here, it is possible to linearly propagate the covariance of any variable computed with respect to the ephemeris and its parameters. Let H(t, P)∈R^{h} such a variable. For a linear random Gaussian variation of P characterized by a covariance matrix cov P, we then obtain the covariance of H at date t,
where J_{H}(t) is the Jacobian matrix of H with respect to P at date t. To compute this matrix, we followed the same procedure as for C(t_{i}, P), which is formally equivalent.
In what follows, we compute the linear covariance propagation in RTN geocentric coordinates^{1}, which are defined according to the following orthonormal basis for any planet A:
where x_{A} represents the barycentric coordinates of body A in the International Celestial Reference Frame ICRF, x_{EMB} are those of the EarthMoon barycenter, and the cross represents the vectorial product. We computed the quantities R_{A}, T_{A}, and N_{A} as
Wed define a rotation matrix M from ICRF to RTN and used it to rotate the covariance matrix. It is interesting to analyze the evolution of the RTN components of a specific body because this provides a more intuitive interpretation of the uncertainties along the observed direction.
In order to account for the uncertainty in the position of the EarthMoon barycenter that is now the origin of the new coordinate frame, we used Eq. (5). In this case, the Jacobian J can be written as
with
and P^{*} ∈ ℝ^{p − 1}. Considering only body A, the EMBcentered covariances therefore become with
The terms cov (x) represent the variance of x, or its covariance matrix, and cov (x_{A},x_{EMB}) is the crosscovariance matrix between x_{A} and x_{EMB}.
In Fig. 1, the evolution of the covariance of the position of Saturn as obtained from the INPOP19a solution is shown, transformed into EMBcentered RTNframe using Eqs. (9) and (12). The radial direction, constrained by Cassini range measurements (see Table 1), is estimated much more accurately than the transverse and normal directions.
Fig. 1. Evolution of the 3σ uncertainty for Saturn positions within the Cassini measurement interval. In particular, the three components of the EMBcentered RTNframe uncertainties are shown in meters. 
Considering either the RTNframe or the ICRFframe, we compared the difference between the components and the evolution of the covariance for two sets of parameters P_{1} and P_{2} in order to compute a distance between two ephemerides. We compared this for two different models to determine whether the difference between the two ephemerides represented by the two sets of parameters P_{1} and P_{2} is within the uncertainty ellipsoid that was estimated with the propagation of the covariances.
After we obtained the propagated covariance at time t, we compared the ith P9perturbed solution with respect to the reference solution, INPOP19a. The match between the two propagated covariances of the two solutions at that specific time was therefore assessed using a generalized distance normalized using the covariance metric^{2}
The general match between the two solutions throughout the analyzed period was computed as the percentage of times for which the distance d(t) is within the equivalent 3σ interval of cov P(t). A compatibility of 100% means that the distance d(t) remains within the 3σ interval during the entire considered time interval, whereas a 50% match indicates that during half of the time interval, the distance d(t) is outside of the 3σ covariance interval.
3.3. Likelihood
3.3.1. Method
We considered the INPOP19a χ^{2} noted , and the χ^{2} of the fitted ephemeris, noted χ^{2}(P9) obtained by including the perturbations of P9 with a mass M_{P9} located at a given position (RA, Dec, r). The χ^{2} is defined as nχ^{2} = ^{t}RWR, where n is the total number of observations, and R and W are as defined in Sect. 4.1. It is well known that if the postfit residuals are in a linear vicinity of 0 and follow a Gaussian distribution, nχ^{2} follows a n degrees of freedom χ^{2} law, and when n → ∞
In other words, z_{P9}, as defined previously, tends to follow a normal distribution centered around 0 and of standard dispersion equal to 1. For the sake of efficiency, it is proposed to focus on observations that are the most sensitive to P9 accelerations. In this case, one can compute the reduced χ^{2}, noted , which can be related to the actual χ^{2}s by:
where is the number of sensitive data (see Sect. 3.3.2) and n the total number of data used for the fit. It appears that the standard dispersions of the residuals of the reference solution are very close to the instrumental uncertainties, such that with a rescaling of the instrumental uncertainties, one can set without loss of generality (see Bernus et al. 2020 for the full demonstration). One can then compute that
follows a 0centered normal distribution of dispersion equal to 1. From here we deduce the likelihood of each ephemeris obtained with a given position of P9 and mass value with:
By definition, L(INPOP19a) is equal to 0.5 and any solution with L(P9) ≈ 0.5 is as likely as INPOP19a. If one solution has L(P9) > 0.5, this solution is then more likely than INPOP19a, meaning with smaller residuals. Solutions with L < 0.5 are less likely than INPOP19a. For this work, as for the classical gaussian distributed variable, we take the equivalent of the 3σ criterion: theories for which L(P9) < 0.003 will be rejected with a probability of 0.997. An advantage on the likelihood criterion is that it tells if a tested solution improves significantly the reference solution which is not the case with the matching criterion.
3.3.2. Selection of sensible data sets
Figures 2 and 3 plot the variations of the root mean square (RMS) of the postfit residuals of the most accurate INPOP19a data samples computed for two examples (r = 800 AU, M_{P9} = 10 M_{⊕} and r = 600 AU, M_{P9} = 5 M_{⊕}, respectively). As one can see, the most important variations are induced on the Cassini data sets. In Fig. 2, the variations can be as big as 5.9 m for CassJ, 2.9 m for CassN and 14 m for CassG, representing variations, with respect to the average accuracy, of about 25%, 50% and several hundreds of %, respectively. For Juno, Mars orbiters and Messenger (MSG), the variations represent less than 10% of the average accuracies. The same conclusions can be drawn from Fig. 3. In this context, for the selection of the sensible datasets for the likelihood computation, we select the three Cassini samples. The Mars orbiter observations are also selected as they represent 47% of the full data sample and have then an important role in the likelihood computation despite their weak sensitivity to P9. We also include the Juno and the Messenger data samples even if they do not contribute much to the χ^{2} computation.
Fig. 2. Differences between INPOP19a postfit residuals and postfit residuals of (RA, Dec) solutions obtained with r = 800 AU and M_{P9} = 10 M_{⊕}. The xaxis represents the RA of P9 when the yaxis is the Dec of P9. The zaxis gives the differences in meters for different observational data sets between INPOP19a residuals and the ephemerides including the P9 perturbations. Several observational subsamples are considered: CassN corresponds to the data sample presented in Sect. 2.1.1 that is labeled Nav+TGF range in Table 1, CassJ is the data sample analyzed and distributed by the JPL (see Hees et al. 2014), and finally, CassG corresponds to the Saturn positions deduced from the Grand Finale as described in Sect. 2.1.1. Juno, Mars, and MSG indicate the variations in postfit residuals for the Jupiter observations by Juno, described in Sect. 2.1.2; Mars means the observations of Mars, and MSG those of Mercury provided by Verma & Margot (2016). The Juno residuals are improved relative to INPOP19a in the gray zones. 
4. Results
4.1. Propagation of the covariance
In Fig. 4, are plotted the maps of the matching criterion based on the propagated covariance of Saturn (as defined in Sect. 3.2, hereafter called T1), obtained for r = 400, 500, 600, 650, 700, 750, 800 AU, M_{P9} = 5, 10 M_{⊕} and different (RA, Dec) positions of P9. We chose the Saturn orbit as a marker because it is the most sensible to P9 perturbations, as shown in Figs. 2 and 3. In Fig. 4, it is visible that the zones maximizing T1, Z_{1} and Z_{2} given in Table 2, are very similar whatever considered distances or masses. Only the percentage of T1 changes, increasing when the mass of P9 decreases or its distance to the SSB. For the two zones given in Table 2, the existence of P9 cannot be statistically rejected as far as T1 is concerned.
Fig. 4. Percentage of matches based on the propagated covariance for Saturn, considering r = 400, 500, 600, 650, 700, 750, and 800 AU and M_{P9} = 5 and 10 M_{⊕}. The xaxis represents the RA of P9 when the yaxis is its Dec. The zaxis gives the matches (see Sect. 4.1 for a defintion) in percent between the INPOP19a propagated covariance of the Saturn orbit over the time coverage of the data and the Saturnorbitpropagated covariance of the ephemerides including perturbations by P9. 
Zones maximizing criteria T1 and T2.
For the 10 M_{⊕} cases, zones with T1 > 75% are visible only for r > 650 AU and enlarge when the distance increases: at 700 AU, the zones with T1 > 75% correspond to 0.55% of the cases, when for 800 AU, these zones correspond to 11% of the cases with a maximum T1 of 92% (versus T1 = 89.5% for 700 AU). If one considers the results of the Batygin et al. (2019) simulations with a 10 M_{⊕} P9, with a maximum possible P9 distance from the SSB of 640 AU, we can see that less than 0.7% of the cases have a T1 greater than 50% .
In the cases where M_{P9} = 5 M_{⊕} (as for 10 M_{⊕}), the zones with T1 > 75% increase with the distances: at 600 AU, these zones correspond to 3.7% of the cases when they represent 23.5% of the 800 AU cases with a maximum T1 of about 96.6 % (versus 93.4% for the 600 AU cases). We finally note that for the 5 M_{⊕} cases, zones with T1 > 75% are visible for r > 500 AU.
4.2. Likelihood
In Fig. 5 are shown maps of the χ^{2} likelihood (hereafter called criterion T2) obtained for r = 400, 500, 600, 650, 700, 800 AU, M_{P9} = 5, 10 M_{⊕} and different P9 (RA, Dec) positions.
Fig. 5. Likelihood considering r = 400, 500, 600, 650, 700, 750, and 800 AU and M_{P9} = 5, 10 M_{⊕}. The xaxis represents the RA of P9 when the yaxis is its Dec. The zaxis gives the likelihood (see Sect. 3.3 for a definition) of the ephemerides that include perturbations by P9 with respect to INPOP19a. The white portions of the maps correspond to solutions rejected at 3σ (L < 0.003). 
First thing to say is that whatever the considered distances or masses, there is no improvement of the ephemerides by adding P9. The likelihood stays indeed below 0.5 for all tested configurations. The relevant question is then if P9 does not improve the residuals, where does it not degrade them significantly. We can use the definition of acceptable solutions as defined in Sect. 3.3: a region where the likelihood of the solution is at 3σ from the reference solution, meaning that L > 0.003.
For a 10 M_{⊕} P9, in using such a definition, no acceptable regions are possible for distances below 700 AU. Above 700 AU, acceptable regions are possible. They are given in Table 2 with r = 800 AU. Among the four identified regions, two correspond to T1zones (Z_{1} and Z_{2}, see Sect. 4.1) and two (Z_{3} and Z_{4}) are newly proposed by T2. We stress that for these four zones, the likelihoods are still smaller than 0.3 (corresponding to 1σ probability). One can also see that by increasing the P9 distance from the SSB, its dynamical influence is getting weaker and so the surface of the acceptable zones increases.
For the cases where the mass of P9 is equal to 5M_{⊕} and with r = 600 and 800 AU, one can see that the acceptable zones get reduced. At 600 AU and 5 M_{⊕}, some regions are acceptable but with a probability greater than 75% to be rejected (L < 0.25). These regions are localized in zones very similar to the one noticed at 700 AU with 10 M_{⊕}. At r = 800 AU, the 5 M_{⊕} hypothesis gives obviously larger acceptable zones than with M_{P9} = 10 M_{⊕} as the acceleration induced by P9 is proportional to its mass and so the impact of the P9 mass on the χ^{2}.
Finally, let notice that some features can be noticed at the edge of the ecliptic plane for a P9 mass equal to 10 M_{⊕}. In order to investigate these regions where the likelihood seems to increase, we densified the simulations to regions out from the ecliptic plane for r = 800 AU and M_{P9} = 10 M_{⊕}. The results are presented in Fig. 6. This figure shows that two regions out of the ecliptic indeed present a likelihood criterion greater than 0.4. For these zones, P9 solutions are almost as likely as INPOP19a and therefore cannot be strictly ruled out. Three other zones are labeled A, B, and C in Fig. 6. Their likelihoods are slightly higher (up to 0.2) than in the remainder of the acceptable zones. These zones are only noticeable with T2, as is shown by the magenta plot of the acceptable zones for the T1 criterion in Fig. 6. Again, the residuals do not improve because the likelihood remains below 0.5 in Fig. 6.
Fig. 6. Likelihood obtained for r = 800 AU and M_{P9} = 10 M_{⊕}: on the lefthand side we show the full map of the likelihood extrapolation. On the righthand side we show the same map, but for a more restricted color map: the white portion of the map corresponds to solutions that have a probability greater than 0.997 (3σ) to be rejected (T2 criterion). The labels A, B, and C indicate the most acceptable zones (with higher likelihoods) in the ecliptic plane. For comparison, the dotted magenta lines indicate the zones for which T1 is greater than 75% (see Fig. 4). 
5. Discussion
5.1. Likelihood and covariance versus induced acceleration
In Fig. 7, we plot for each P9 (RA, Dec) position the induced accelerations on Saturn, Mars, and the EarthMoon barycenter (EMB) orbits. As expected, the P9 acceleration on the Saturn orbit has the largest amplitude compared to those on the other planets, and the maximum of the P9 dynamical effect occurs in the ecliptic plane (represented in black). When we compare Fig. 7 with Fig. 6, the zones with a likelihood close to 0.5 (red and yellow regions in Fig. 6 and levels in Fig. 7) , where solutions with P9 are almost equally probable (L > 0.4) as in INPOP19a, correspond to zones where P9 induces very low accelerations (blue and green regions in Fig. 7) on the most sensitive planetary orbits. These zones are out of the ecliptic plane and are excluded by the theoretical modeling of Batygin et al. (2019). However, it is interesting to note that for these regions, we cannot conclude whether the presence of P9 is possible because the induced dynamical effect on planetary orbits is very weak. Consequently, the effect of the fitted residuals is also very weak; the χ^{2} likelihood is very similar to be that for INPOP19a.
Fig. 7. Accelerations induced by P9 (a_{P9}) on Saturn (lefthand side), Mars (middle) and EarthMoon barycenter (EMB, righthand side) orbits averaged over 20 years. The darkblue dotted levels indicate zones where the likelihood criterion T2 is greater than 0.4. 
5.2. Likelihood and covariance versus fitted TNO ring mass
As explained in Sect. 2.2, the addition of a TNO ring was an important update for fitting Cassini observations. Another important aspect of the adjustment of the mass of the TNO ring is its ability to absorb a part of the dynamical contribution of P9 (Pitjeva & Pitjev 2018). In this context, it was interesting to study the change in fitted value of the mass of the TNO ring according to the different P9 effects on the planetary orbits. In Fig. 8, we plot the differences between the INPOP19a gravitational mass (GM) of the TNO ring, which is equal to 5.32 × 10^{−11} UA^{3} day^{−2} (Di Ruscio et al. 2020; Fienga et al. 2019), and the fitted values for each solution that includes contributions of P9 for different positions (RA, Dec). When we compare Figs. 4 and 6 and do not consider the lowacceleration zones discussed previously (Sect. 5.1), there are two zones (one for RA < 50° and one for RA ≈ 200°) that are acceptable for the likelihood criterion T2 and maximizing T1 (gray levels in Fig. 8) and two zones that are acceptable for T2, but with small T1 (dark blue levels in Fig. 8). These latter correspond to zones where the values of GM_{TNO} decrease significantly (see Fig. 8) with respect to the INPOP19a fitted value. In these cases, the perturbation by P9 is partially compensated for in the fit by a decrease in TNO ring contribution (and therefore of its mass) that leads to a χ^{2} that is still similar to that of INPOP19a, and consequently is similar to an acceptable likelihood. On the other hand, the propagation of the covariance is not sensitive to this mechanism because the leastsquares uncertainties are not affected. This can explain why these zones have small T1 in Fig. 4. The consequence of this mechanism is that the detection of the gravitational signature of P9 on planetary orbits appears to be correlated with the estimation of the TNO masses. As explained previously, the T2 acceptability of a certain area is due to our poor knowledge of these objects and consequently to the possibility that the P9induced perturbations are masked by a misleading estimate of their masses. Independent constraints on the distribution of masses beyond Neptune are necessary in this context.
Fig. 8. Differences of the gravitational mass of the TNO ring induced by P9 relative to the INPOP19a value (Di Ruscio et al. 2020); 2019NSTIM.109.....V. The dashed levels indicate zones where T2 is acceptable (probability of being rejected below 0.997): the grey zones are common to zones maximizing T1 and the dark blue zones are those selected with high T2 but with low T1 
6. Conclusions
We have presented the results of 3156 simulations of planetary ephemerides built using the INPOP19a dynamical modeling and data sampling (Fienga et al. 2019), but adding the perturbations induced by the unknown planet P9, which we considered at different locations in the Solar System and with two different masses (5 M_{⊕} and 10 M_{⊕}). Based on these simulations, we performed two statistical tests: one used the propagation of the covariance matrices of planetary orbits and related parameters (Sect. 3.2, called T1), and the other used the χ^{2} likelihood (Sect. 3.3, called T2). Two main conclusions can be drawn from Figs. 4 and 5. First, the planetary ephemerides are not a positive marker for the existence of P9 beause P9 does not improve the planetary residuals regardless of the configurations that are considered. Second, we can deduce from Figs. 4 and 5 thresholds in the distances between P9 and the SSB below which the existence of P9 is ruled out by the planetary ephemerides and above which it is acceptable in some regions. For T1, the limit corresponds to the detection of zones with T1 > 75%, and for T2, the limit corresponds to the occurrence of acceptable zones as defined in Sect. 3.3. Table 3 presents these limits, which can be compared with the distances proposed by Batygin et al. (2019). In their simulations, the largest possible distance is 640 AU, considering a minimum eccentricity of 0.2 and a maximum semimajor axis of 800 AU. When the largest distance of P9 is compared to the SSB proposed by Batygin et al. (2019) with the limits given in Table 3 and considering as possible zones those acceptable for both criteria, it appears that our work rules out the possibility of a 10 M_{⊕} P9. Only a small P9 with a mass of 5 M_{⊕} might been accepted at a distance smaller than 640 AU.
Smallest possible distances between P9 and the SSB in AU deduced from T1 and T2.
In addition to these two main conclusions, we showed in Sect. 4 that possible zones exist for different distances and masses (see Table 2). Focusing on the case (r = 600 AU and M_{P9} = 5 M_{⊕}), which best agrees with Batygin et al. (2019), we show in Figs. 9 and 10 that two specific regions, denoted D1 and D2 (coordinates given in Table 4), are positive for T1 and T2, and the region even maximizes T2 with a likelihood of 0.09 (below 1.5σ). These regions can then be proposed as interesting zones for further investigations, even if the acceptability statistics in these areas remain low, with T1 below 90% and T2 below 1.5σ (L < 0.1).
Fig. 9. Likelihood obtained for r = 600 AU and M_{P9} = 5 M_{⊕}. The white portion of the map corresponds to solutions that have a probability greater than 0.997 (3σ) to be rejected. The label D indicates the most acceptable zones (with higher likelihoods) in the ecliptic plane. For comparison, the dotted magenta lines indicate the zones for which T1 is greater than 75% according to the propagated covariance (see Fig. 4). 
Fig. 10. Same as Fig. 9 for two specific Dzones: D1 on the lefthand side, and D2 on the righthand side. 
Possible zones for a search for P9 given in intervals of equatorial (RA, Dec) in degrees.
Finally, as we discussed in Sect. 5.2, the presence of P9 induces a dynamical effect that can be absorbed by reducing the mass of the TNO ring. An estimate of this mass obtained independently from the planetary ephemerides will be a crucial help for efficiently distinguishing perturbations caused by P9 and by TNO accelerations on the planetary orbits.
Acknowledgments
This work has been funded by the French Space agency (CNES APR) and the Université Côte d’Azur (UCA Academie 3). AF thanks A. Morbidelli and K. Batygin for very helpful discussions.This work was also supported by the Programme National GRAM of CNRS/INSU with INP and IN2P3 cofunded by CNES.
References
 Bailey, E., Batygin, K., & Brown, M. E. 2016, AJ, 152, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Brown, M. E. 2016a, AJ, 151, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Brown, M. E. 2016b, ApJ, 833, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., Adams, F. C., Brown, M. E., & Becker, J. C. 2019, Phys. Rep., 805, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Bernus, L., Minazzoli, O., Fienga, A., et al. 2020, Phys. Rev. Lett., submitted [Google Scholar]
 Brown, M. E. 2017, AJ, 154, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Di Ruscio, A., Fienga, A., Durante, D., et al. 2020, A&A, 640, A7 [CrossRef] [EDP Sciences] [Google Scholar]
 Durante, D., Hemingway, D. J., Racioppa, P., Iess, L., & Stevenson, D. J. 2019, Icarus, 326, 123 [CrossRef] [Google Scholar]
 Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572 [CrossRef] [Google Scholar]
 Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2016, A&A, 587, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fienga, A., Deram, P., Viswanathan, V., et al. 2019, Notes Scientifiques et Techniques de l’Institut de Mecanique Celeste, 109 [Google Scholar]
 Fienga, A., Avdellidou, C., & Hanuš, J. 2020, MNRAS, 492, 589 [CrossRef] [Google Scholar]
 Folkner, W. M., Iess, L., Anderson, J. D., et al. 2017, Geophys. Res. Lett., 44, 4694 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R. S., Matese, J. J., & Lissauer, J. J. 2006, Icarus, 184, 589 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R., Deienno, R., & Morbidelli, A. 2017, AJ, 153, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Hees, A., Folkner, W. M., Jacobson, R. A., & Park, R. S. 2014, Phys. Rev. D, 89, 102002 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Payne, M. J. 2016, AJ, 152, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Iess, L., Stevenson, D. J., Parisi, M., et al. 2014, Science, 344, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Iess, L., Folkner, W. M., Durante, D., et al. 2018, Nature, 555, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science, 364, eeat2965 [CrossRef] [Google Scholar]
 Kaib, N. A., Pike, R., Lawler, S., et al. 2019, AJ, 158, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Lawler, S. M., Shankman, C., Kaib, N., et al. 2017, AJ, 153, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Madigan, A.M., & McCourt, M. 2016, MNRAS, 457, L89 [NASA ADS] [CrossRef] [Google Scholar]
 Millholland, S., & Laughlin, G. 2017, AJ, 153, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvorný, D., Vokrouhlický, D., Dones, L., et al. 2017, ApJ, 845, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Pitjeva, E. V., & Pitjev, N. P. 2018, Celestial Mech. Dyn. Astron., 130, 57 [Google Scholar]
 Shankman, C., Kavelaars, J. J., Bannister, M. T., et al. 2017a, AJ, 154, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Shankman, C., Kavelaars, J. J., Lawler, S. M., Gladman, B. J., & Bannister, M. T. 2017b, AJ, 153, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Verma, A. K., & Margot, J.L. 2016, J. Geophys. Res. (Planets), 121, 1627 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Smallest possible distances between P9 and the SSB in AU deduced from T1 and T2.
Possible zones for a search for P9 given in intervals of equatorial (RA, Dec) in degrees.
All Figures
Fig. 1. Evolution of the 3σ uncertainty for Saturn positions within the Cassini measurement interval. In particular, the three components of the EMBcentered RTNframe uncertainties are shown in meters. 

In the text 
Fig. 2. Differences between INPOP19a postfit residuals and postfit residuals of (RA, Dec) solutions obtained with r = 800 AU and M_{P9} = 10 M_{⊕}. The xaxis represents the RA of P9 when the yaxis is the Dec of P9. The zaxis gives the differences in meters for different observational data sets between INPOP19a residuals and the ephemerides including the P9 perturbations. Several observational subsamples are considered: CassN corresponds to the data sample presented in Sect. 2.1.1 that is labeled Nav+TGF range in Table 1, CassJ is the data sample analyzed and distributed by the JPL (see Hees et al. 2014), and finally, CassG corresponds to the Saturn positions deduced from the Grand Finale as described in Sect. 2.1.1. Juno, Mars, and MSG indicate the variations in postfit residuals for the Jupiter observations by Juno, described in Sect. 2.1.2; Mars means the observations of Mars, and MSG those of Mercury provided by Verma & Margot (2016). The Juno residuals are improved relative to INPOP19a in the gray zones. 

In the text 
Fig. 3. Same as Fig. 2, but for r = 600 AU and M_{P9} = 5 M_{⊕}. 

In the text 
Fig. 4. Percentage of matches based on the propagated covariance for Saturn, considering r = 400, 500, 600, 650, 700, 750, and 800 AU and M_{P9} = 5 and 10 M_{⊕}. The xaxis represents the RA of P9 when the yaxis is its Dec. The zaxis gives the matches (see Sect. 4.1 for a defintion) in percent between the INPOP19a propagated covariance of the Saturn orbit over the time coverage of the data and the Saturnorbitpropagated covariance of the ephemerides including perturbations by P9. 

In the text 
Fig. 5. Likelihood considering r = 400, 500, 600, 650, 700, 750, and 800 AU and M_{P9} = 5, 10 M_{⊕}. The xaxis represents the RA of P9 when the yaxis is its Dec. The zaxis gives the likelihood (see Sect. 3.3 for a definition) of the ephemerides that include perturbations by P9 with respect to INPOP19a. The white portions of the maps correspond to solutions rejected at 3σ (L < 0.003). 

In the text 
Fig. 6. Likelihood obtained for r = 800 AU and M_{P9} = 10 M_{⊕}: on the lefthand side we show the full map of the likelihood extrapolation. On the righthand side we show the same map, but for a more restricted color map: the white portion of the map corresponds to solutions that have a probability greater than 0.997 (3σ) to be rejected (T2 criterion). The labels A, B, and C indicate the most acceptable zones (with higher likelihoods) in the ecliptic plane. For comparison, the dotted magenta lines indicate the zones for which T1 is greater than 75% (see Fig. 4). 

In the text 
Fig. 7. Accelerations induced by P9 (a_{P9}) on Saturn (lefthand side), Mars (middle) and EarthMoon barycenter (EMB, righthand side) orbits averaged over 20 years. The darkblue dotted levels indicate zones where the likelihood criterion T2 is greater than 0.4. 

In the text 
Fig. 8. Differences of the gravitational mass of the TNO ring induced by P9 relative to the INPOP19a value (Di Ruscio et al. 2020); 2019NSTIM.109.....V. The dashed levels indicate zones where T2 is acceptable (probability of being rejected below 0.997): the grey zones are common to zones maximizing T1 and the dark blue zones are those selected with high T2 but with low T1 

In the text 
Fig. 9. Likelihood obtained for r = 600 AU and M_{P9} = 5 M_{⊕}. The white portion of the map corresponds to solutions that have a probability greater than 0.997 (3σ) to be rejected. The label D indicates the most acceptable zones (with higher likelihoods) in the ecliptic plane. For comparison, the dotted magenta lines indicate the zones for which T1 is greater than 75% according to the propagated covariance (see Fig. 4). 

In the text 
Fig. 10. Same as Fig. 9 for two specific Dzones: D1 on the lefthand side, and D2 on the righthand side. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.