Reconciling results of 2019 and 2020 stellar occultations on Pluto’s atmosphere New constraints from both the 5 September 2019 event and consistency analysis

A stellar occultation by Pluto on 5 September 2019 yielded positive detections at two separate stations. Using an approach consistent with comparable studies, we derived a surface pressure of 11 . 478 ± 0 . 55 µ bar for Pluto’s atmosphere from the observations of this event. In addition, to avoid potential method inconsistencies when comparing with historical pressure measurements, we reanalyzed the data for the 15 August 2018 and 17 July 2019 events. All the new measurements provide a bridge between the two di ff erent perspectives on the pressure variation since 2015: a rapid pressure drop from previous studies of the 15 August 2018 and 17 July 2019 events and a plateau phase from that of the 6 June 2020 event. The pressure measurement from the 5 September 2019 event aligns with those from 2016, 2018, and 2020, supporting the latter perspective. While the measurements from the 4 June 2011 and 17 July 2019 events suggest probable V-shaped pressure variations that are unaccounted for by the volatile transport model (VTM), the VTM remains applicable on average. Furthermore, the validity of the V-shaped variations is debatable given the stellar faintness of the 4 June 2011 event and the grazing single-chord geometry of the 17 July 2019 event. To reveal and understand all of the significant pressure variations of Pluto’s atmosphere, it is essential to provide constraints on both the short-term and long-term evolution of the interacting atmosphere and surface by continuous pressure monitoring through occultation observations whenever possible, and to complement these with frequent spectroscopy and photometry of the surface.

A compilation of twelve occultations observed between 1988 and 2016 revealed a three-fold monotonic increase in the atmospheric pressure of Pluto during that period (Meza et al. 2019).This increase can be explained by the volatile transport model (VTM) of the Laboratoire de Météorologie Dynamique (LMD) (Bertrand & Forget 2016;Forget et al. 2017;Bertrand et al. 2018Bertrand et al. , 2019) ) fine-turned by Meza et al. (2019).This model provides a framework for simulating the volatile cycles on Pluto's globe over both seasonal and astronomical timescales, allowing for an understanding of the long-term evolution of Pluto's atmosphere and its response to seasonal variations over its 248-year heliocentric orbital period (Meza et al. 2019).According to the LMD VTM in Meza et al. (2019) (VTM19, hereafter), Pluto's atmospheric pressure is expected to reach its peak around 2020.The pressure increase is attributed to the progression of summer over the northern hemisphere of Pluto, exposing Sputnik Planitia (SP) 1 to solar radiation.The surface of SP, composed of nitrogen (N 2 ), methane (CH 4 ), and carbon monoxide (CO) ices, is believed to sublimate and release volatile gases into the atmosphere during this period, leading to a pressure increase.After reaching its peak, the model predicts a gradual decline in pressure over the next two centuries under the combined effects of Pluto's recession from the Sun and the prevalence of the winter season over SP.
On one hand, the VTM19 remains consistent with Sicardy et al. ( 2021)'s analysis of the 6 June 2020 occultation observed at Devasthal, where two co-located telescopes were used.Their analysis suggests that Pluto's atmosphere has been in a plateau phase since mid-2015.This aligns with the model predictions that the atmospheric pressure reached its peak around 2020.
On the other hand, Arimatsu et al. (2020)'s analysis of the 17 July 2019 occultation observed by a single telescope (TUHO) suggests a rapid pressure decrease between 2016 and 2019.They detected a significant pressure drop at the 2.4σ level.However, it is worth noting that the geometry of this occultation is grazing.This may have introduced larger correlations between the pressure and the geocentric closest approach distance to Pluto's shadow axis, leading to the lack of the necessary precision to confidently support the claim of a large pressure decrease followed by a return to the pressure level close to that of 2015 in 2020 (Sicardy et al. 2021).
These contrasting results highlight the need for occultation observations between 2019 and 2020 to better understand the behavior and evolution of Pluto's atmosphere during this time period.Furthermore, while Young et al. (2021) support the presence of a pressure drop based on their analysis of the 15 August 2018 occultation, Sicardy et al. (2021) suggest that careful comparisons between measurements by independent teams should be made before drawing any conclusions on the pressure evolution.
Observations of the 5 September 2019 occultation which has not been reported by other teams are presented in Section 2, followed by a description of the light-curve fitting methods in Section 3.These unique observations allow us to track the changes in Pluto's atmosphere during the time period between the events studied by Arimatsu et al. (2020) and Sicardy et al. (2021).Results are detailed in Section 4, and the pressure evolution is discussed in Section 5, including comparisons with the reanalyzed 15 August 2018 and 17 July 2019 events.Conclusions and recommendations are provided in Section 6.

Occultation observations
Two observation campaigns were organized in China for occultations in 2019 (see Appendix A).One occurred on 17 July 2019, which has been studied by Arimatsu et al. (2020), and the other on 5 September 2019, which is reported in the present paper for the first time.Due to bad weather conditions in many areas, no effective light curves were observed by our stations for the first occultation, and only two for the second.
Table 1 lists the circumstances of the 5 September 2019 event.Figure 1 presents all the observation stations and the reconstructed path of the shadow of Pluto2 during this event.Table A.2 lists the circumstances of stations with positive detections.Their station codes are DWM and HNU.
To ensure accurate and precise timing in stellar occultations, some stations (e.g., DWM as shown in Table A.2) were equipped with QHY174GPS cameras.These cameras, manufactured by  2022), used to derive (α s , δ s ). (b) International Celestial Reference Frame.
(c) Stern et al. (2015), where G is the gravitational constant. (d) Assumed to be the only constituent in the light-curve model, following, e.g., Sicardy et al. (2021).
(h) Positive (resp.negative) values mean that the shadow center passes north (resp.south) of the geocenter. (i) Timings by QHY174GPS camera are used as time references, considering their reliability and accuracy.
QHYCCD3 , offer precise recording of observation time and location for each frame using a GPS-based function.And, they have been used in many stellar occultation studies (e.g., Buie et al. 2020a,b;Morgado et al. 2021Morgado et al. , 2022;;Pereira et al. 2023).
In the light-curve fitting procedures described in Section 3.2, the time recording offsets of the QHY174GPS cameras are fixed to zero, considering their reliability and accuracy as time references.
All observational data were captured in the FITS format.These data were processed using the Tangra occultation photometric tool4 (Pavlov 2020) and our data reduction code (see Appendix B).It was ensured that the targets and reference stars in all the images we used were not overexposed.The resulting light curves from the observations, after being normalized, are presented in Figure 2.Each data point on the light curves is represented by f i (t) ± σ i (t), where i indicates the quantities associated with a specific station, t represents the recorded timing of each frame, f the normalized total observed flux of the occulted star and the Pluto's system, and σ the measurement error associated with each data point.

Light-curve model
In order to simulate observed light curves, we implemented a light-curve model, ϕ(t; A, s, ∆t, ∆τ, ∆ρ, p 0 ), described in Appendix C and consistent with the DO15 (Dias-Oliveira et al. 2015; Sicardy et al. 2016;Meza et al. 2019;Sicardy et al. 2021).As a function of model parameters, its time-dependent Jacobian matrix was also implemented to represent the sensitivity of the model with respect to the corresponding parameters that will be estimated through fitting procedures.
The light-curve model of a given station can be formally written as: where i indicates the quantities associated with the station, and for details the reader is referred to Appendix C. Here, the reference ephemerides we use are the NIMAv9 5 asteroidal ephemeris (Desmars et al. 2015(Desmars et al. , 2019) ) for the orbit of the Pluto system barycenter with respect to the Sun, the PLU058 6 satellite ephemerides (Brozović et al. 2015;Jacobson et al. 2019) for the orbit of Pluto with respect to the Pluto system barycenter, and the DE440 7 planetary ephemerides (Park et al. 2021) for the orbits of the Earth and the Sun with respect to the Solar system barycenter.The reference star catalog where the data of the occulted star are obtained is the Gaia DR3.

Fitting procedure
The light-curve model was fitted to the normalized observed light curves simultaneously by nonlinear least squares, returning a χ 2 type value of goodness-of-fit.The goal was to minimize the objective function given by: where t i j represents the mid-exposure time of the j-th observation of the station i.
In addition, with the used reference ephemerides and star catalog, a piece of a priori information on ∆ρ can be obtained: where the uncertainty σ ρ is set to 72 km using the positional uncertainties listed in NIMAv9's "orbit quality" table and in Gaia DR3.This σ ρ value corresponds to about 3 mas on the sky at the geocentric distance of Pluto.The a priori information can be treated as an independent observational data and used in the model fitting, with the objective function modified as: The fitting steps are as following: -In order to find all local minima, at some of which a nonlinear least-squares fitting could get stuck, we explored the two parameter space (∆ρ, p 0 ) by generating the variation of χ 2 as a function of them.Figure 3 presents such two χ 2 maps, labeled (a) and (b), which will be analyzed in Section 4. The maps are generated by minimizing χ 2 obs or χ 2 apr at each fixed (∆ρ, p 0 ) point on a regularly spaced grid.The Levenberg-Marquardt (LM) method, implemented in the LMFIT package8 , was used in each fitting procedure.The free parameters to be adjusted are ∆τ9 , ∆t i of any station having no reliable time reference system like QHY174GPS, and s i and A i of each station.
-For a more accurate best-fitting solution for (∆ρ, p 0 ), the LM method is used again, with ∆ρ and p 0 adjusted with initial guesses located at all known local minima of each χ 2 map.-Each χ 2 map, which provides information about the quality of the fit, is used to define confidence limits based on constant χ 2 boundaries (Press et al. 2007).

Results
Figure 3a shows the χ 2 obs map for the 5 September 2019 occultation.Two local minima are observed.However, considering the significant χ 2 difference of 9 between the two local minima, the global minimum is more likely the correct solution.In addition, the ∆ρ value at the global minimum is more consistent with the NIMAv9 solution, ∆ρ = 0 km, at the 0.16 σ ρ level, compared with the other local one at the 2.44 σ ρ level.
In an effort to mitigate or, at least, further weaken the presence of multiple local minima, we calculated the χ 2 apr map by adding the χ 2 type value of the a priori informaiton, (∆ρ/σ ρ ) 2 , into the χ 2 obs map. Figure 3b presents the results, showing two local minima still being present, but with a χ 2 difference of about 14.5, which is larger than that of the χ 2 obs map.So, the global minimum is confidently accepted as the solution for (ρ cag , p surf ), as provided in Table 1.
Moreover, Figure 3 presents the consistency of our derived p surf across the two different local minima.It demonstrates that the specific choice of local minima does not significantly affect the value of p surf , further supporting the reliability of our solution for p surf .

Comparisons and necessary reanalyses of historical events
In Figure 4, the red plot represents our p surf measurement from the 5 September 2019 occultation.We also include other published measurements (Hinson et al. 2017;Meza et al. 2019;Arimatsu et al. 2020;Young et al. 2021;Sicardy et al. 2021) and the pressure evolution predicted by the VTM19 to provide a comprehensive view of the pressure variations on Pluto.
To avoid potential inconsistencies arising from different analysis methods, as discussed by Sicardy et al. (2021), we reanalyzed the 15 August 2018 event studied by Young et al. (2021) using the IXON observational data of Silva-Cabrera et al. (2022).The derived pressure measurement, presented in Appendix D.1, is 12.027 +0.09 −0.08 µbar.In addition, we also reanalyzed the 17 July 2019 event in Appendix D.2, deriving a pressure of p surf = 9.421 +0.68  −0.75 µbar, similar to that of Arimatsu et al. (2020), 9.56 +0.52  −0.34 µbar.This similarity is expected because the same DO15 method is used.Since the method is used by Meza et al. (2019) and Sicardy et al. (2021), their pressure measurements, along with that of Arimatsu et al. (2020), can be fully compared with our new ones.Both the remeasurements are plotted in black in Figure 4.
The pressure measurement from the 5 September 2019 event shows alignments with those from the 19 July 2016, 15 August 2018, and 6 June 2020 events within their combined 1σ levels.Our new measurement from the 15 August 2018 event shows no significant pressure drop previously reported by Young et al. (2021).The previously reported pressure drop between the 19 July 2016 and 17 July 2019 events is still detected at the same level as in Arimatsu et al. (2020).

Discussions on pressure variations
While the VTM19, on average, remains applicable and capable of predicting the main atmospheric behavior during the observed years, there are also two probable V-shaped pressure variations observed from 2010 to 2015 and from 2015 to 2020, especially when considering the measurements from the 4 June 2011 and 17 July 2019 events.These V-shaped variations suggest the presence of additional and unaccounted factors.Specifically, shortterm changes in Pluto's surface ices and their interaction with the atmosphere are likely contributing to the variation.In fact, spectral monitoring of the surface composition has revealed some short-term changes in the ices over several Earth years (e.g., Grundy et al. 2014;Lellouch et al. 2022;Holler et al. 2022).
However, the validity of the V-shaped variations is debatable due to the stellar faintness of the 4 June 2011 event and the grazing single-chord geometry of 17 July 2019 event.If the debatable measurement from the 17 July 2019 event is discarded, no significant changes would be observed between 2016 and 2020.This more likely supports the plateau phase since 2015 predicted by the VTM19.In order to better understand the relationship between these factors, further observations using multiple observa- apr maps, respectively.The χ 2 obs denotes the goodness-of-fit only using observational data, while the χ 2 apr with the additional a priori information represented as the χ 2 type value, (∆ρ/σ ρ ) 2 .These maps are used to derive the best-fitting atmospheric pressure p 0 at the reference radius r 0 of 1215 km and the cross track correction ∆ρ to the ephemerides, of which the a priori uncertainty σ ρ is 72 km.The surface pressure p surf and the geocentric closest approach distance ρ cag to the shadow center are obtained by linear transformations of ∆ρ and p 0 , respectively.The best-fitting χ 2 obs and χ 2 apr values per degree of freedom are 1.160 and 1.155, respectively.tional techniques (occultation, spectroscopy, and photometry), as well as simulations with a refined VTM, are required.

Conclusions
The unique observations of the 5 September 2019 occultation have provided a surface pressure of p surf = 11.478± 0.55 µbar.In order to avoid potential method inconsistencies in comparing with historical pressure measurements (Sicardy et al. 2021), we also reanalyzed the 15 August 2018 and 17 July 2019 events based on publicly available data (Silva-Cabrera et al. 2022;Arimatsu et al. 2020).All measurements are presented in Figure 4.
The VTM19 remains applicable on average.Besides, we also observed unaccounted V-shaped pressure variations with the previously reported pressure drop being a part of these variations, which, however, are debatable.To better understand all significant pressure variations of Pluto, continuous pressure monitoring through occultation observations, whenever possible, is essential.Also, simultaneous and frequent spectroscopic and photometric monitoring of its surface ice changes are important, as  Then, the values of r ± can be numerically solved by using the relationship between z and r ± : z = ±(r ± + D • ω(r ± )), (C.4) where D represents the light travel distance, which can be approximated by the geocentric distance of Pluto, and ω(r) is the total deviation angle at r.As detailed in Dias-Oliveira et al. (2015) and Sicardy (2022), ω(r) can be obtained by a ray-tracing code that follows these steps: set the temperature profile T (r) of Pluto using Equation (4) from Dias-Oliveira et al. (2015), with parameters obtained from Sicardy et al. (2016); set the gas molecular refractivity K corresponding to a stellar wavelength of λ µm ; set the atmospheric pressure p 0 at the reference radius r 0 of 1215 km either directly, or determined from a surface pressure p surf (at R p = 1187 km) with the ratio p surf /p 0 = 1.837 used by previous studies (Meza et al. 2019;Sicardy et al. 2021); Fig. 1: Reconstructed occultation map of the 5 September 2019 event.

Fig. 2 :
Fig. 2: Occultation observations and the best-fitting light-curve model of the 5 September 2019 event.Panel a: observed and simultaneously fitted light curves.Panels b and c: Reconstructed stellar paths seen by DWM and HNU, respectively.

Fig. 3 :
Fig. 3: The χ 2 maps of the 5 September 2019 event.Panels a and b: the χ 2obs and χ 2 apr maps, respectively.The χ 2 obs denotes the goodness-of-fit only using observational data, while the χ 2 apr with the additional a priori information represented as the χ 2 type value, (∆ρ/σ ρ ) 2 .These maps are used to derive the best-fitting atmospheric pressure p 0 at the reference radius r 0 of 1215 km and the cross track correction ∆ρ to the ephemerides, of which the a priori uncertainty σ ρ is 72 km.The surface pressure p surf and the geocentric closest approach distance ρ cag to the shadow center are obtained by linear transformations of ∆ρ and p 0 , respectively.The best-fitting χ 2 obs and χ 2 apr values per degree of freedom are 1.160 and 1.155, respectively.

Fig. 4 :
Fig. 4: The pressure evolution of Pluto over a 248-year heliocentric orbital period predicted by the VTM19, along with the measured pressures with 1σ error bars.The 3σ error bars of our new measurements and of the previous measurements of the 19 July 2016, 17 July 2019, and 6 June 2020 events are also presented, using the published χ 2 maps.A N 2 denotes the albedo of nitrogen ice.

Table 1 :
Circumstances and light-curve fitting results of the 5 September 2019 event.