New constraints on the location of P9 obtained with the INPOP19a planetary ephemeris

Context. We used the new released INPOP19a planetary ephemerides benefiting from Jupiter-updated 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.


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 trans-Neptunian 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. This version benefits from the use of nine very accurate Jupiter-normal 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. 2). 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 Notes. Column 1 gives the observed planet and information on the type of observations, and Col. 2 indicates the number of observations. Columns 3 and 4 give the time interval and the a priori uncertainties provided by space agencies or the navigation teams, respectively. Finally, the WRMS for INPOP19a are given in the last column.
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.

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.

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. 2019Iess et al. 2014Iess et al. , 2019. For each arc, we solved for the spacecraft initial position and velocity, corrections for orbital trim and reaction wheel desaturation maneuvers, and RTG-induced 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 ground-station 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 gravity-dedicated flybys of Titan (Durante et al. 2019) and Grand Finale Saturn pericenters ). These range normal points were obtained by considering a pass-through 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).

Nine perijoves of Juno
The Juno spacecraft has been orbiting Jupiter on a highly eccentric polar orbit since July 2016. A radio-science experiment aims at characterizing the gravity field of the gas giant to unprecedented accuracy Folkner et al. 2017;Iess et al. 2018). The extremely accurate radio-tracking system of Juno enables simultaneous two-way Doppler measurements at X and Ka band with accuracies as low as 10 µs −1 in the radial velocity during the gravity-dedicated passes, which are used to reconstruct the spacecraft trajectory with a meter-level 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 Jupiter-normal points spanning the period from the orbital insertion, back in 2016, to the end of 2018. 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).

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ẍ P9 A such thaẗ 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.

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 Cauchy-Lipschitz equation-of-motion system, where X ∈ R 6n , with n being the number of integrated bodies, and P ∈ R 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 trans-Neptunian objects; oblateness of the Sun; and the Earth-Moon 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 O(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 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, cov H(t, P) = J H (t) cov P t J H (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 Earth-Moon 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 Earth-Moon 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 * ∈ R p−1 . Considering only body A, the EMB-centered 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 cross-covariance 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 EMB-centered RTN-frame 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.
Considering either the RTN-frame or the ICRF-frame, 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 i-th P9-perturbed 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.

Method
We considered the INPOP19a χ 2 noted χ 2 r , 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χ 2 , 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χ 2 r = 1 without loss of generality (see Bernus et al. 2020 for the full demonstration). One can then compute that follows a 0-centered 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.  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.

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 T 1), 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 T 1, Z 1 and Z 2 given in Table 2, are very similar whatever considered distances or masses. Only the percentage of T 1 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 T 1 is concerned. For the 10 M ⊕ cases, zones with T 1 > 75% are visible only for r > 650 AU and enlarge when the distance increases: at 700 AU, the zones with T 1 > 75% correspond to 0.55% of the cases, when for 800 AU, these zones correspond to 11% of the cases with a maximum T 1 of 92% (versus T 1 = 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 T 1 greater than 50% . In the cases where M P9 = 5 M ⊕ (as for 10 M ⊕ ), the zones with T 1 > 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 T 1 of about 96.6 % (versus 93.4% for the 600 AU cases). We finally note that for the 5 M ⊕ cases, zones with T 1 > 75% are visible for r > 500 AU.
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 T 1-zones (Z 1 and Z 2 , see Sect. 4.1) and two (Z 3 and Z 4 ) are newly proposed by T 2. 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 T 2, as is shown by the magenta plot of the acceptable zones for the T 1 criterion in Fig. 6. Again, the residuals do not improve because the likelihood remains below 0.5 in Fig. 6.

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

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 low-acceleration 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 T 2 and maximizing T 1 (gray levels in Fig. 8) and two zones that are acceptable for T 2, but with small T 1 (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 least-squares uncertainties are not affected. This can explain why these zones have small T 1 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 T 2 acceptability of a certain area is due to our poor knowledge of these objects and consequently to the possibility that the P9-induced perturbations are masked by a misleading estimate 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 x-axis represents the RA of P9 when the y-axis is its Dec. The z-axis 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 Saturn-orbit-propagated covariance of the ephemerides including perturbations by P9. of their masses. Independent constraints on the distribution of masses beyond Neptune are necessary in this context.

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 T 1), and the other used the χ 2 likelihood (Sect. 3.3, called T 2). 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 T 1, the limit corresponds  On the right-hand 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 (T 2 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 T 1 is greater than 75% (see Fig. 4).  . The dashed levels indicate zones where T 2 is acceptable (probability of being rejected below 0.997): the grey zones are common to zones maximizing T 1 and the dark blue zones are those selected with high T 2 but with low T 1 Table 3. Smallest possible distances between P9 and the SSB in AU deduced from T 1 and T 2.
to the detection of zones with T 1 > 75%, and for T 2, 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. Notes. In zones D1 and D2, T 1 and T 2 are maximized in the ecliptic. The coordinates given here correspond to the zones in which the likelihood is greater than 0.05 (2σ). 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 T 1 is greater than 75% according to the propagated covariance (see Fig. 4).
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 T 1 and T 2, and the region even maximizes T 2 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 A6, page 10 of 11 A. Fienga et al.: INPOP19a constraints on planet 9 Fig. 10. Same as Fig. 9 for two specific D-zones: D1 on the left-hand side, and D2 on the right-hand side. statistics in these areas remain low, with T 1 below 90% and T 2 below 1.5σ (L < 0.1).
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.