© A. Fienga et al. 2020

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 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, 2020). 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. 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.

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 RTG-induced anisotropic acceleration. In addition, stochastic accelerations at the level of 5 × 10⁻¹³ km s⁻² (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 (Iess et al. 2019). 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).

2.1.2. 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 (Durante et al. 2020; 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⁻¹ 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.

2.2. Trans-Neptunian 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 ${\ddot{\mathit{x}}}_A^{\mathrm{P}9}$ such that

$$\begin{array}{c}\hfill {\ddot{\mathit{x}}}_{\mathit{iA}}^{\mathrm{P}9}=\frac{G{M}_{\mathrm{P}9}}{r^3}[3({\mathit{x}}_A·\mathit{u})\mathit{u}-{\mathit{x}}_A],\end{array}$$(1)

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 χ² 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 Cauchy-Lipschitz equation-of-motion system,

$$\begin{array}{c}\hfill \frac{\mathrm{d}X}{\mathrm{d}t}=F(X;P),X(t=J2000)=X_0,\end{array}$$(2)

where X ∈ ℝ⁶ⁿ , 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 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₀. 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

$$\begin{array}{c}\hfill R(t_i,P)=(\mathcal{O}(t_i,P)-C(t_i)),\end{array}$$(3)

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 χ² = ^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

$$\begin{array}{c}\hfill \delta P={({}^{t}JWJ)}^{-1}{}^{t}JWR\end{array}$$

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

$$\begin{array}{c}\hfill \mathrm{cov}P={({}^{t}JWJ)}^{-1}.\end{array}$$(4)

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,

$$\begin{array}{c}\hfill \mathrm{cov}H(t,P)=J_H(t)\mathrm{cov}P{}^t{J}_H(t),\end{array}$$(5)

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¹, which are defined according to the following orthonormal basis for any planet A:

$$\begin{array}{cc}& {\mathit{u}}_A=\frac{{\mathit{x}}_A-{\mathit{x}}_{\mathit{EMB}}}{|{\mathit{x}}_A-{\mathit{x}}_{\mathit{EMB}}|},\hfill \end{array}$$(6)

$$\begin{array}{cc}& {\mathit{w}}_A=\frac{{\mathit{u}}_A\times {\dot{\mathit{x}}}_A}{|{\mathit{u}}_A\times {\dot{\mathit{x}}}_A|},\hfill \end{array}$$(7)

$$\begin{array}{cc}& {\mathit{v}}_A={\mathit{w}}_A\times {\mathit{u}}_A,\hfill \end{array}$$(8)

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

$$\begin{array}{cc}& R_A=({\mathit{x}}_A-{\mathit{x}}_{\mathit{EMB}})·{\mathit{u}}_A\hfill \\ \hfill & T_A=({\mathit{x}}_A-{\mathit{x}}_{\mathit{EMB}})·{\mathit{v}}_A\hfill \\ \hfill & N_A=({\mathit{x}}_A-{\mathit{x}}_{\mathit{EMB}})·{\mathit{w}}_A.\hfill \end{array}$$(9)

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

$$\begin{array}{c}\hfill J=\frac{\partial f(P)}{\partial P}=\frac{\partial (P^{\ast }-X_{\mathit{EMB}})}{\partial P},\end{array}$$(10)

with

$$\begin{array}{c}\hfill X_{\mathit{EMB}}=[\underset{n-1}{\underbrace{{\mathit{x}}_{\mathit{EMB}},\cdots ,{\mathit{x}}_{\mathit{EMB}}}},0,\cdots ]\in {\mathbb{R}}^{p-1}\end{array}$$(11)

and P^* ∈ ℝ^p − 1. Considering only body A, the EMB-centered covariances therefore become with ${\mathit{x}}_A^E={\mathit{x}}_A-{\mathit{x}}_{\mathit{EMB}}$

$$\begin{array}{c}\hfill \mathrm{cov}\left({\mathit{x}}_A^E\right)=\mathrm{cov}\left({\mathit{x}}_A\right)+\mathrm{cov}\left({\mathit{x}}_{\mathit{EMB}}\right)-2\mathrm{cov}({\mathit{x}}_A,{\mathit{x}}_{\mathit{EMB}}).\end{array}$$(12)

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.

Fig. 1. Evolution of the 3σ uncertainty for Saturn positions within the Cassini measurement interval. In particular, the three components of the EMB-centered RTN-frame uncertainties are shown in meters.

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₁ and P₂ 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₁ and P₂ 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²

$$\begin{array}{c}\hfill d(t)=\sqrt{({\mathit{x}}_{A_{\mathrm{ref}}}^E(t)-{\mathit{x}}_{A_i}^E(t))·{\left(\mathrm{cov}{\mathit{x}}_A^E\right)}^{-1}·{}^t({\mathit{x}}_{A_{\mathrm{ref}}}^E(t)-{\mathit{x}}_{A_i}^E(t))}.\end{array}$$(13)

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 χ² noted ${\chi }_r^2$, and the χ² of the fitted ephemeris, noted χ²(P9) obtained by including the perturbations of P9 with a mass M_P9 located at a given position (RA, Dec, r). The χ² is defined as nχ² =  ^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χ² follows a n degrees of freedom χ² law, and when n → ∞

$$\begin{array}{c}\hfill z_{\mathrm{P}9}=\sqrt{\frac{n}{2}}({\chi }^2(P9)-{\chi }_r^2)\underset{n\to \infty }{\to }\mathcal{N}(0,1).\end{array}$$(14)

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 χ², noted ${\stackrel{\sim }{\chi }}^2$, which can be related to the actual χ²s by:

$$\begin{array}{c}\hfill {\chi }^2(P9)-{\chi }_r^2=\frac{\stackrel{\sim }{n}}{n}({\stackrel{\sim }{\chi }}^2(P9)-{\stackrel{\sim }{\chi }}_r^2)\end{array}$$(15)

where $\stackrel{\sim }{n}$ 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 ${\stackrel{\sim }{\chi }}_r^2=1$ without loss of generality (see Bernus et al. 2020 for the full demonstration). One can then compute that

$$\begin{array}{c}\hfill z_{\mathrm{P}9}=\sqrt{\frac{\stackrel{\sim }{n}}{2}}({\stackrel{\sim }{\chi }}^2(P9)-1)\end{array}$$(16)

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:

$$\begin{array}{c}\hfill \mathrm{L}(\mathrm{P}9)=1-\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{z_{\mathrm{P}9}}}exp(-\frac{x^2}{2})\mathrm{d}x.\end{array}$$(17)

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 χ² computation.

Fig. 2. Differences between INPOP19a postfit residuals and postfit residuals of (RA, Dec) solutions obtained with r = 800 AU and MP9 = 10 M⊕. The x-axis represents the RA of P9 when the y-axis is the Dec of P9. The z-axis 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.

Fig. 3. Same as Fig. 2, but for r = 600 AU and MP9 = 5 M⊕.

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₁ and Z₂ 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 MP9 = 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.

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 χ² 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 MP9 = 5, 10 M⊕. The x-axis represents the RA of P9 when the y-axis is its Dec. The z-axis 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 T1-zones (Z₁ and Z₂, see Sect. 4.1) and two (Z₃ and Z₄) 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 χ².

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 MP9 = 10 M⊕: on the left-hand side we show the full map of the likelihood extrapolation. 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 (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 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 χ² likelihood is very similar to be that for INPOP19a.

Fig. 7. Accelerations induced by P9 (aP9) on Saturn (left-hand side), Mars (middle) and Earth-Moon barycenter (EMB, right-hand 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⁻¹¹ UA³ day⁻² (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 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 χ² 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 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 P9-induced 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 χ² 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.

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 MP9 = 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 D-zones: D1 on the left-hand side, and D2 on the right-hand side.

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 co-funded by CNES.

References

All Tables

All Figures

	Fig. 1. Evolution of the 3σ uncertainty for Saturn positions within the Cassini measurement interval. In particular, the three components of the EMB-centered RTN-frame 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 MP9 = 10 M⊕. The x-axis represents the RA of P9 when the y-axis is the Dec of P9. The z-axis 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 MP9 = 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 MP9 = 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.

In the text

	Fig. 5. Likelihood considering r = 400, 500, 600, 650, 700, 750, and 800 AU and MP9 = 5, 10 M⊕. The x-axis represents the RA of P9 when the y-axis is its Dec. The z-axis 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 MP9 = 10 M⊕: on the left-hand side we show the full map of the likelihood extrapolation. 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 (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 (aP9) on Saturn (left-hand side), Mars (middle) and Earth-Moon barycenter (EMB, right-hand 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 MP9 = 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 D-zones: D1 on the left-hand side, and D2 on the right-hand side.
In the text