Pluto's ephemeris from ground-based stellar occultations (1988-2016)

From 1988 to 2016, several stellar occultations have been observed to characterize Pluto's atmosphere and its evolution (Meza et al, 2019). From each stellar occultation, an accurate astrometric position of Pluto at the observation epoch is derived. These positions mainly depend on the position of the occulted star and the precision of the timing. We present Pluto's astrometric positions derived from 19 occultations from 1988 to 2016 (11 from Meza et al. (2019) and 8 from other publications). Using Gaia DR2 for the positions of the occulted stars, the accuracy of these positions is estimated to 2-10~milliarcsec depending on the observation circumstances. From these astrometric positions, we derive an updated ephemeris of Pluto's system barycentre using the NIMA code (Desmars et al., 2015). The astrometric positions are derived by fitting the occultation's light curves by a model of Pluto's atmosphere. The fits provide the observed position of the body's centre for a reference star position. Other publications usually provide circumstances of the occultation such as the coordinates of the stations, the timing, and the impact parameter (i.e. the closest distance between the station and the centre of the shadow). From these parameters, we use a procedure based on the Bessel method to derive an astrometric position. We derive accurate Pluto's astrometric positions from 1988 to 2016. These positions are used to refine the orbit of Pluto'system barycentre providing an ephemeris, accurate to the milliarcsec level, over the period 2000-2020, allowing better predictions for future stellar occultations.


Introduction
Stellar occultation is a unique technique to obtain the physical parameters of distant objects or to probe their atmosphere and surroundings. For instance, Meza et al. (2019) have used 11 stellar occultations by Pluto from 2002 to 2016 to study the evolution of Pluto's atmosphere. Meanwhile, occultations allow an accurate determination of the relative position of the body's centre compared to the position of the occulted star, leading to an accurate astrometric position of Pluto at the time of occultation if the star position is also known accurately.
The accuracy of the body's position mainly depends on the knowledge of the shape and the size of the body, the modelling of the atmosphere, the precision of the timing system, the velocity of the occultation, the exposure time of the camera, the precision of the stellar position, and the magnitude of the occulted star. Since the publication of the Gaia catalogues in September 2016 for the first release (Gaia Collaboration et al. 2016) and moreover with the second release in April 2018 (Gaia Collaboration et al. 2018) including proper motions and trigonometric parallaxes of the stars, the precision of the stellar catalogues can now reach a tenth of a milliarcsec (mas). For comparison, before Gaia catalogues, the precision of stellar catalogues such as UCAC2 (Zacharias et al. 2004) or UCAC4 (Zacharias et al. 2013), was about 50 to 100 mas per star with also zonal errors. With Gaia, the precision of positions deduced from occultations is expected to be around few mas, taking into account the systematic errors. Thanks to the accuracy of the GaiaDR2 catalogue, occultations can provide the most accurate astrometric measurement of a body in the outer solar system.
In this paper, we present the astrometric positions we derived from occultations presented in Meza et al. (2019)  astrometric positions from other publications, knowing the circumstances of occultations: timing and impact parameter (Appendix). Finally, we present in Section 3 a refined ephemeris of Pluto's system barycentre and we discuss the improvement in the predictions of future occultations by Pluto at a mas level accuracy and as well as the geometry of past events (Sect. 4). Beyond the parameters related to Pluto's atmosphere, another product of the occultations is the astrometric position of the body. From the geometry of the event, we determine the position of the Pluto's centre of figure (α c , δ c ). This position corresponds to the observed position of the object at the time of the occultation for a given star position (α s , δ s ). In particular, the position of the body we derive, only depends on the star position. Before Gaia catalogues, we determined the star position with our own astrometry. Table 1 gives the position of Pluto's centre and its precision we derived from the geometry of the occultation and the corresponding star position from our astrometry. With Gaia, the astrometric position of Pluto's centre can be refined by correcting the star position with the Gaia DR2 star position with the relations:

Astrometric positions from occultations
This refined position only depends on the Gaia DR2 position which is much more accurate than previous catalogues or our own astrometry. The associated position of the occulted stars in Gaia DR2 catalogue (α GDR2 , δ GDR2 ) are listed in Table 2. The positions take into account the proper motions and the parallax from Gaia DR2. The table also presents the Gaia source identifier and the estimated precision of the star position in right ascension and declination, at the time of the occultation, taking into account precision of the stellar position and the proper motions as given in GaiaDR2, for all the occultations studied in this paper.
Finally, Table 3 provides the absolute position in right ascension and declination of Pluto's centre derived from the geometry and from stellar positions of Gaia DR2. The residuals related to JPL ephemeris 1 DE436/PLU055 are also indicated as well as the differential positions between Pluto and Pluto's system barycentre used to refine the orbit (see Sect. 3). A flag indicates if the position is used in the NIMAv8 ephemeris determination. Finally, the reconstructed paths of the occultations are presented in Fig. 6.

Astrometric positions from other publications
Several authors have published circumstances of an occultation by Pluto (e.g. Millis et al. 1993;Sicardy et al. 2003;Elliot et al. 1 DE436 is a planetary ephemerides from JPL providing the positions of the barycentre of the planets, including the barycentre of Pluto's system. It is based on DE430 (Folkner et al. 2014). PLU055 is the JPL ephemeris providing the positions of Pluto and its satellites related to the Pluto's system barycentre, developed by R.Jacobson in 2015 and based on an updated ephemeris of Brozović et al. (2015): https://naif.jpl.nasa.gov/pub/naif/generic_ kernels/spk/satellites/plu055.cmt Young et al. 2008;Person et al. 2008;Gulbis et al. 2015;Olkin et al. 2015;Pasachoff et al. 2016;Pasachoff et al. 2017). From these circumstances (coordinates of the observer, mid-time of the occultation and impact parameter), it is possible to derive an offset between the observation deduced from these circumstances and a reference ephemeris. The procedure, based on the Bessel method used to predict stellar occultations, is described in Appendix A and the details of computation for each occultation are presented in Appendix B. The Pluto's positions deduced from occultations published in other articles besides those of Meza et al. (2019) are presented in Table 3.
The positions derived from Pasachoff et al. (2016) involving single chord events and faint occulted stars, are not accurate enough to discriminate North and South solutions, i.e. to decide if Pluto's centre as seen from the observing site passed North or South of the star. Finally, these positions were not used in the orbit determination.

NIMA ephemeris of Pluto
NIMA (Numerical Integration of the Motion of an Asteroid) has been developed in order to refine the orbits of small bodies, in particular TNOs and Centaurs studied by stellar occultations technique (Desmars et al. 2015). It consists of the numerical integration of the equations of motion perturbed by gravitational accelerations of the planets (Mercury to Neptune). The Earth and Moon are considered at their barycentre and the masses and the positions of the planets are from JPL DE436.
The use of other masses and positions for planetary ephemeris produces insignificant changes, for example, the difference between the solution using DE436 and the solution using INPOP17a (Viswanathan et al. 2017) for Pluto, is less than 0.06 mas on the 1985-2025 period. Moreover, there is no need to take into account the gravitational perturbations of the biggest TNOs. For example, by adding the 6 biggest TNOs (Eris, Haumea, 2007 OR10, Makemake, Quaoar and Sedna) in the model, the difference between the solutions with and without the biggest TNOs are about 0.04 mas in right ascension and declination on the 1985-2025 period, which is 100 times smaller than the mas-level accuracy of the astrometric positions.
The state vector (the heliocentric vector of position and velocity of the body at a specific epoch) is refined by fitting to astrometric observations with the least square method. The main advantage of NIMA is allowing for the use of observations published on the Minor Planet Center 2 together with unpublished observations or astrometric positions of occultations. The quality of the observations is taken into account with a specific weighting scheme, in particular, it takes advantages of the high accuracy of occultations. Finally, after fitting to the observations, NIMA can provide an ephemeris through a bsp file format readable by the SPICE library 3 .
As NIMA is representing the motion of the centre of mass of an object, it allows to compute the position of the Pluto's system barycentre and not the position of Pluto's centre itself. To deal with positions derived from occultations, we need an additional ephemeris representing the position of Pluto relative to its system barycentre. For that purpose, we use the most recent Table 1: Date, timing, position of Pluto's centre deduced from the geometry and the precision, coordinates of the occulted star used to derive the astrometric positions of occultations by Pluto studied in Meza et al. (2019).

Reference date
Pluto's centre position position of star right ascension σ α declination σ δ right ascension declination α c (mas) δ c (mas) α s δ s ephemeris PLU055 developed in 2015. The occultation-derived positions are then corrected from the offset between Pluto and the Pluto's system barycentre (see Table 3) to derive the barycentric positions from the occultations, then used in the NIMA fitting procedure.
The precisions of the positions in right ascension and in declination derived from the occultations are provided in Table 1 for occultations presented in Meza et al. (2019) and in Appendix B for other publications. This precision is deduced from a specific model and reduction (for occultations in Meza et al. 2019) and from the precisions of timing and impact parameters (for other publications), without any estimation of systematic errors. For a realistic estimation of the orbit accuracy, the weighting scheme in the orbit fit needs to take into account the systematic errors (see Desmars et al. 2015 for details). The global accuracy for the positions used in the fitting depends on the accuracy of the stellar positions (from 0.1 to 2 mas), the precision of the derived position (from 0.1 mas to 11 mas), and the accuracy of the Pluto body-Pluto system barycentre ephemeris (estimated to 1-5 mas).
The errors on Pluto's centre determination have in fact various sources: the noise present in each occultation light curve and the spatial distribution of the occultation chords across the body. Assuming a normal noise, a formal error on the centre of the planet can then be derived, using a classical least-square fitting and χ 2 estimation. However, other systematic errors may also be present, such as problems in the absolute timing registration, slow sky transparency variations that make the photometric noise non-gaussian. Finally, the particular choice of the atmospheric model may also induce systematic biases in the centre determination. All those systematic errors are difficult to trace back. In that context, it is instructive to compare the reconstructions of the geometry of a given occultation by independent groups that used different chords and different Pluto's atmospheric models. For example, occultations on 21 August 2002, 4 May 2013 and29 Article number, page 3 of 15 A&A proofs: manuscript no. aa34958 Table 3: Right ascension and declination of Pluto deduced from occultations, residuals (O-C) in mas related to JPL DE436/PLU055 ephemeris, and differential coordinates (PLU-BAR) between Pluto and Pluto barycentre system position from PLU055 ephemeris. June 2015 (see Table 4) indicate differences of few mas, which is much higher than the respective internal precisions (order of 0.1 mas). Case by case studies should be undertaken to explain those inconsistencies. This remains out of the scope of this paper. Meanwhile, for the weighting scheme in the orbit fit, we adopt the estimated precision presented in Table 4 taking into account an estimation of systematic errors for each occultation.  Table 4 and Fig. 2 provide the residuals (O-C) of the positions derived from the occultations, compared with the NIMAv8 ephemeris, and the estimated precision of the positions used in the weighting scheme. After 2011, residuals are mostly below the mas level, which is much better than any ground-based astrometric observation of Pluto. In that context, other classical observations of Pluto, such as CCD, appear to be less useful for ephemerides of Pluto during the period covered by the occultations 1988-2016. Figure 3 shows the difference in right ascension and declination between the most recent ephemerides of Pluto system  and EPM2017 (Pitjeva & Pitjev 2014) compared to NIMAv8. These differences are mostly due to data and weights used for the orbit determination. They reveal periodic terms in the orbit of Pluto system barycentre that are differently estimated in orbit determination. As described in Desmars et al. (2015), the one-year period corresponds to the parallax induced by different geocentric distances given by the ephemerides. It is also another good indication of the improvement of the NIMAv8 ephemeris since the differences between these ephemerides reach 50-100 mas whereas the estimated precision of NIMAv8 is 2-20 mas on the same period.

Discussion
The NIMA ephemeris allows very accurate predictions of stellar occultation by Pluto in the forthcoming years within a few mas level. In particular, we have predicted an occultation of a magnitude 13 star 5 by Pluto on August 15, 2018, above North America to the precision of 2.5 mas, representing only 60 km on the shadow path and a precision of 4 s in time. As shown in Meza et al. (2019), the observation of a central flash allows to probe the deepest layers of Pluto's atmosphere. The central flash can be observed in an small band about 50 km around the centrality path. By reaching a precision of tens of km, we were able to gather observing stations along the centrality and to highly increase the probability of observing a central flash.
The prediction of the August 15, 2018 Pluto occultation was used to assess the accuracy of our predictions using the NIMA approach. Figure 4 represents the prediction of the occultation by Pluto on August 15, 2018 using two different ephemerides: JPL DE436/PLU055 and NIMAv8/PLU055. The prediction using JPL ephemerides is shifted by 36.8 s and 8 mas south (repre-senting about 190 km) compared to the prediction with NIMAv8 ephemeris. Several stations detected the occultation, some of them revealing a central flash. For instance, observers at George Observatory (Texas, USA) report a central flash of typical amplitude 20%, compared to the unocculted stellar flux (T. Blank and P. Maley, private communications).
As the amplitude of the flash roughly scales as the inverse of the closest approach (C/A) distance of the station to the shadow centre, the amplitude may serve to estimate the C/A distance. A central flash reported by Sicardy et al. (2016) was observed at a station in New Zealand during the June 29, 2015 occultation. It had an amplitude of 13%, with C/A distance of 42 km. Thus, the flash observed at George Observatory provides an estimated C/A distance of 25 km for that station. This agrees with the value predicted by the NIMAv8/PLU055 ephemeris, to within 3 km, corresponding to 0.12 mas. This is fully consistent with, but smaller than our 2.5 mas error bar quoted above, possibly indicating an overestimation of our prediction errors.
The precision of our predictions remains at few mas up to 2025 (in particular in declination) and it is even more important since the apparent position of Pluto as seen from Earth is moving away from the Galactic centre, making occultations by Pluto more rare.
Another point of interest is to look at past occultations. In particular, for the occultation of August 19, 1985, Brosch & Mendelson (1985 reported a single chord occultation of a magnitude 11.1 star 6 by Pluto, showing a gradual shape possibly due to Pluto's atmosphere. The observation was performed at Wise observatory in Israel under poor conditions (low elevation, flares from passing planes, close to twilight). Thanks to Gaia DR2 providing the proper motion of the star and to NIMAv8, we make a postdiction of the occultation of August 19, 1985 (Fig. 5). The nominal time for the occultation (the time of the closest approach between the geocentre and the centre of the shadow) is 17:58:57.1 (UTC) leading to a predicted mid-time of 17:59:49.8 (UTC) at Wise observatory. Brosch (1995) gave an approximate observed mid-time of the occultation for Wise observatory at 17:59:54 (about 4 s later than the prediction). The predicted shadow of Pluto at the same time is presented on the figure as well as the observatory's place as a green bullet. Taking into account the uncertainties of the NIMAv8 ephemeris and of the star position, the uncertainty in time for this occultation is about 20 s whereas the crosstrack uncertainty on the path is about 10 mas (representing 220 km). This is fully consistent with the fact that the occultation was indeed observed at Wise observatory.

Conclusions
Stellar occultations by Pluto provide accurate astrometric positions thanks to Gaia catalogues, in particular Gaia DR2. We determine 18 astrometric positions of Pluto from 1988 to 2016 with an estimated precision of 2 to 10 mas.
These positions are used to compute an ephemeris of Pluto system's barycentre thanks to NIMA procedure with an unprecedented precision on the 1985-2015 period. This ephemeris NI-MAv8 was used to study the possible occultation of Pluto observed in 1985 as well to predict the recent occultation by Pluto on August 15th, 2018 or the forthcoming occultations 7 with a precision of 2 mas, a result impossible to reach with classical astrometry and previous stellar catalogues. In fact, the presence of the usually unresolved Charon in classical images, causes significant displacements of the photocentre of the system with respect to its barycentre. As a consequence, and even modeling the effect of Charon, as in Benedetti-Rossi et al. (2014), accuracies below the 50 mas level are difficult to reach.
This method can be extended, for instance for Chariklo, with an even better accuracy of the order of 1 mas (Desmars et al. 2017) and illustrates the power of stellar occultations not only for better studying those bodies, but also for improving their orbital elements. In fact, t m , u, v, w are calculated iteratively by replacing λ c by λ c − Ω(t m − T 0 ), where Ω is the rate of Earth's rotation, to take into account the Earth's rotation during t m − T 0 .
If ∆t is the difference between the observed time of the occultation for the observer and the nominal time of the occultation T 0 , the correction to apply to the Besselian elements x 0 , y 0 are : The quantities ∆x, ∆y are determined iteratively and finally transformed into an offset in right ascension and in declination between the observed occultation and the predicted occultation (from the ephemeris).
For single chord occultation, there are two solutions (North and South), meaning that we do not know whether Pluto's center went North or South of the star as seen from the observing site. Conversely, for multi-chord occultation there is a unique solution. In that case, the astrometric position deduced from the occultation is the reference ephemeris plus the average offset deduced from all the observing sites. This is a powerful method to derive astrometric positions from occultations. It only requires local circumstances of the occultation for the observing sites such as the mid-time of the occultation and the impact parameter. If the impact parameter is not provided, one can deduce it from the timing of immersion and emmersion knowing the size of the object and assuming it is spherical. Thus, the method can be used for any object.