Pluto's lower atmosphere and pressure evolution from ground-based stellar occultations, 1988-2016

Context. Pluto's tenuous nitrogen (N2) atmosphere undergoes strong seasonal effects due to high obliquity and orbital eccentricity, and has been recently (July 2015) observed by the New Horizons spacecraft. Goals are (i) construct a well calibrated record of the seasonal evolution of surface pressure on Pluto and (ii) constrain the structure of the lower atmosphere using a central flash observed in 2015. Method: eleven stellar occultations by Pluto observed between 2002 and 2016 are used to retrieve atmospheric profiles (density, pressure, temperature) between $\sim$5 km and $\sim$380 km altitude levels (i.e. pressures from about 10 microbar to 10 nanobar). Results: (i) Pressure has suffered a monotonic increase from 1988 to 2016, that is compared to a seasonal volatile transport model, from which tight constraints on a combination of albedo and emissivity of N2 ice are derived; (ii) A central flash observed on 2015 June 29 is consistent with New Horizons REX profiles, provided that (a) large diurnal temperature variations (not expected by current models) occur over Sputnik Planitia and/or (b) hazes with tangential optical depth of about 0.3 are present at 4-7 km altitude levels and/or (c) the nominal REX density values are overestimated by an implausibly large factor of about 20% and/or (d) higher terrains block part of the flash in the Charon facing hemisphere.


Introduction
The tenuous atmosphere of Pluto was glimpsed during a groundbased stellar occultation observed on 1985 August 19 (Brosch 1995), and fully confirmed on 1988 June 9 during another occultation (Hubbard et al. 1988;Elliot et al. 1989;Millis et al. 1993) that provided the main features of its structure: temperature, composition, pressure, and density; see the review by Yelle & Elliot (1997).
Since then, Earth-based stellar occultations have been a highly efficient method to study the atmosphere of Pluto. They yield, in the best cases, information from a few kilometers above the surface (pressure ∼10 µbar) up to 380 km in altitude (∼10 nbar). As Pluto moved in front of the Galactic center, the yearly rate of stellar occultations dramatically increased during the 2002-2016 period, yielding a few events per year that greatly improved our knowledge of the atmospheric structure and evolution of this planet.
Ground-based occultations also provided a decadal monitoring of the atmosphere. Pluto has a large obliquity (∼120 • , the axial inclination to its orbital plane) and high orbital eccentricity (0.25) that takes the dwarf planet from 29.7 to 49.3 AU during half of its 248-year orbital period. Northern spring equinox occurred in January 1988 and perihelion occurred soon after in September 1989. Consequently, our survey monitored Pluto as it receded from the Sun while exposing more and more of its northern hemisphere to sunlight. More precisely, as of 2016 July 19 (the date of the most recent occultation reported here), the heliocentric distance of Pluto has increased by a factor of 1.12 since perihelion, corresponding to a decrease of about 25% in average insolation. Meanwhile, the subsolar latitude has gone from 0 • at equinox to 54 • north in July 2016. In this context, dramatic seasonal effects are expected, and observed.
Another important aspect of ground-based occultations is that they set the scene for the NASA New Horizons mission (NH hereafter) that flew by the dwarf planet in July 2015 (Stern et al. 2015). A fruitful and complementary comparison between the ground-based and NH results ensued -another facet of this work.
Here we report results derived from eleven Pluto stellar occultations observed between 2002 and 2016, five of them yet unpublished, as mentioned below. We analyze them in a unique and consistent way. Including the 1988 June 9 occultation results, and using the recent surface ice inventory provided by NH, we constrain current seasonal models of the dwarf planet. Moreover, a central flash observed during the 2015 June 29 occultation is used to compare the lower atmosphere structure of Pluto derived from the flash with profiles obtained by the Radio Science EXperiment instrument (REX) on board NH below an altitude of about 115 km.
the present work. In a second part of the table, we list other campaigns that were not used because the occultation light curves had insufficient signal-to-noise ratios (S/Ns) and/or because of deficiencies in the configuration of the occulting chords (grazing chords or single chord) and as such, do not provide relevant measurements of the atmospheric pressure. Details on the prediction procedures can be found in Assafin et al. (2010Assafin et al. ( , 2012 and Benedetti-Rossi et al. (2014). Some of those campaigns are already documented and analyzed in previous publications, namely the 2002July 20, 2002August 21, 2007June 14, 2008June 22, 2012July 18, 2013May 4, and 2015 June 29 events. They were used to constrain Pluto's global atmospheric structure and evolution (Sicardy et al. 2003(Sicardy et al. , 2016Dias-Oliveira et al. 2015;French et al. 2015;Olkin et al. 2015), the structure and composition (CH 4 , CO and HCN abundances) of the lower atmosphere by combination with spectroscopic IR and sub-millimeter data (Lellouch et al. 2009(Lellouch et al. , 2017, the presence of gravity waves (Toigo et al. 2010;French et al. 2015), and Charon's orbit (Sicardy et al. 2011). Finally, one campaign that we organized is absent from Table A .1 (2006 April 10). It did not provide any chord on Pluto, but was used to put an upper limit on its rings (Boissel et al. 2014).
We note that five more (yet unpublished) data sets are included here that were obtained on the following dates June 24, 2010February 14, 2010June 4, 2011June 4, and 2016 July 19.

Light-curve fitting
For all the eleven data sets used here, we used the same procedure as in Dias-Oliveira et al. (2015;DO15 hereafter) and in Sicardy et al. (2016). The procedure consists of simultaneously fitting the refractive occultation light curves by synthetic profiles generated by a ray-tracing code that uses the Snell-Descartes law. The physical parameters adopted in this code are listed in Table 1.
We note in particular that our adopted radius for Pluto is taken from Stern et al. (2015), who use a global fit to full-disk images provided by the Long-Range Reconnaissance Imager (LORRI) of NH to obtain R P = 1187 ± 4 km. Nimmo et al. (2017) improve that value to R P = 1188.3 ± 1.6 km. However, we kept the 1187 km value because it is very close to the deepest level reached by the REX experiment, near the depression Sputnik Planitia; see Sect. 4. Consequently, it is physically more relevant here when discussing the lower atmospheric structure.
We assume a pure N 2 atmosphere, which is justified by the fact that the next most important species (CH 4 ) has an abundance of about 0.5% (Lellouch et al. 2009Gladstone et al. 2016), resulting in negligible effects on refractive occultations.
We also assume a transparent atmosphere, which is supported by the NH findings. As discussed in Sect. 4, the tangential (line-of-sight) optical depth of hazes found by NH for the rays that graze the surface is τ T ∼ 0.24, with a scale height of ∼50 km (Gladstone et al. 2016;Cheng et al. 2017). As our fits are mainly sensitive to levels around 110 km (see below), this means that haze absorption may be neglected in our ray-tracing approach. We return to this topic in Sect. 4.3, which considers the effect of haze absorption on the central flash, possibly caused by the deepest layers accessible using occultations.
Moreover, we take a global spherically symmetric atmosphere, which is again supported by the NH results, at least above the altitude ∼35 km; see Hinson et al. (2017) and Sect. 4. This is in line with global climate models (GCMs), which predict that wind velocities in the lower atmosphere should not exceed v ∼1-10 m s −1 (Forget et al. 2017). If uniform, this wind would create an equator-to-pole radius difference of the corresponding isobar level of at most ∆r ∼ (R P v) 2 /4 GM P < 0.1 km, using Eq. (7) of Sicardy et al. (2006) and the values in Table 1. This expected distortion is too small to significantly affect our synthetic profiles.
Finally, the temperature profile T (r) is taken to be constant. Here, the radius r is counted from the center of Pluto, while the radius found by NH is 1187 km (Table 1). This is the reference radius from which we calculate altitudes. Fixing the pressure at a prescribed level (e.g., the surface) then entirely defines the density profile n(r) to within a uniform scaling factor for all radii r, using the ideal gas equation, hydrostatic equilibrium assumption, and accounting for the variation of gravity with altitude.
Taking T (r) to be constant with time is justified by the fact that the pressure is far more sensitive to surface temperaturethrough the vapor pressure equilibrium equation -than is the profile T (r) to seasonal effects and heliocentric distance, at least from a global point of view. For instance, an increase of 1 K of the free N 2 ice at the surface is enough to multiply the equilibrium pressure by a factor of 1.7 (Fray & Schmitt 2009). We note that this is not inconsistent with our assumption that T (r) is time-independent. In fact, the overall atmospheric pressure is controlled by the temperature a few kilometers above the surface, while our fits use a global profile T (r) well above the surface.
Pluto ground-based stellar occultations probe, for the best data sets, altitudes from ∼5 (pressure level ∼10 µbar) to ∼380 km (∼10 nbar level); see DO15. Rays coming from below ∼5 km are detectable only near the shadow center (typically within 50 km) where the central flash can be detected. The analysis is then complicated by the fact that double (or multiple) stellar images contribute to the flux. Moreover, the possible presence of hazes and/or topographic features can reduce the flux; see Sect. 4.
Conversely, rays coming from above 380 km cause stellar drops that are too small (<∼1%) to be of any use under usual ground-based observing conditions. That said, our ray-tracing method is mainly sensitive to the half-light level, where the star flux has been reduced by 50%. This currently corresponds to a radius of about 1295 km (or an altitude ∼110 km and pressure ∼1.6 µbar).

Primary results
The ray-tracing code returns the best fitting parameters, in particular the pressure at a prescribed radius (e.g., the pressure p surf at the surface, at radius R P = 1187 km) and the ephemeris offset of Pluto perpendicular to its apparent motion, ∆ρ. The ects and heliocentric distance, at least from a global point of view. For instance, an increase of 1 K Fig. 1. Example of a χ 2 (∆ρ, p surf ) map derived from the simultaneous fit to the light curves obtained during the 2016 July 19 occultation. The quantity ∆ρ is the ephemeris offset of Pluto (expressed in kilometers) perpendicular to its apparent motion, as projected in the sky plane. The other parameter (p surf ) is the surface pressure of the DO15 atmospheric model. The white dot marks the best fit, where the minimum value χ 2 min of χ 2 is reached. The value χ 2 min = 4716, using 4432 data points, indicates a satisfactory fit with a χ 2 per degree of freedom of χ 2 dof ∼ 4716/4432 ∼ 1.06. The best fit corresponds to p surf = 12.04 ± 0.41 µbar (1-σ level). The error bar is derived from the 1-σ curve that delineates the χ 2 min + 1 level. The 3-σ level curve (corresponding to the χ 2 min + 9 level) is also shown. ephemeris offset along the motion is treated separately; see DO15 for details. Error bars are obtained from the classical function χ 2 = N 1 [(φ i,obs − φ i,syn )/σ i ] 2 which reflects the noise level σ i of each of the N data points, where φ i,obs and φ i,syn are the observed and synthetic fluxes, respectively. An example of a χ 2 (∆ρ, p surf ) map is displayed in Fig. 1, using a simultaneous fit to the 2015 June 29 occultation light curves. It shows a satisfactory fit for that event, χ 2 dof ∼1.06. Table 2  Two main consequences of those results are now discussed in turn: (1) the temporal evolution of atmospheric pressure on Pluto; (2) the structure of its lower atmosphere using the central flash of 2015 June 29. A third product of these results is the update of the ephemeris using the occultation geometries between 2002 and 2016. This latter will be presented in a separate paper (Desmars et al. 2019).

Constraints from occultations
In 2002, a ground-based stellar occultation revealed that the atmospheric pressure on Pluto had increased by a factor of almost two compared to its value in 1988 (Elliot et al. 2003;Sicardy et al. 2003), although Pluto had receded from the Sun, thus globally cooling down. Models using global volatile Fig. 4. Typical modeled annual evolution of surface pressure obtained with LMD Pluto volatile transport model, assuming permanent deposits of N 2 ice inside Sputnik Planitia and in the depression of mid-northern latitudes, a uniform soil seasonal thermal inertia of 800 J s −1/2 m −2 K −1 , an emissivity N2 = 0.8 and albedo range A N2 = 0.72-0.73 for N 2 ice, chosen to yield a surface pressure near 10-11 µbar in July 2015. The black dots with error bars show the surface pressure (p surf ) inferred from stellar occultation pressure measurements (see Table 2). The curve in magenta corresponds to a similar simulation but assuming a permanent N 2 ice reservoir in the south hemisphere between 52.5 and 67.5 • S, which leads to a pressure peak in 1990. transport predicted this seasonal effect, among different possible scenarios (Binzel 1990;Hansen & Paige 1996).
Those models explored nitrogen cycles, and have subsequently been improved (Young 2012(Young , 2013Hansen et al. 2015). Meanwhile, new models have been developed to simulate possible scenarios for changes over seasonal (248 yr) and astronomical (30 Myr) time scales, accounting for topography and ice viscous flow, as revealed by the NH flyby in July 2015 (Bertrand & Forget 2016;Forget et al. 2017;Bertrand et al. 2018).
The measurements obtained here provide new values of pressure versus time, and are obtained using a unique light curve fitting model (taken from DO15), except for the 1988 occultation; see Table 2. This model may introduce systematic biases, but it can nevertheless be used to derive the relative evolution of pressure from date to date, and thus discriminates the various models of the current seasonal cycle of Pluto. In any case, the DO15 light-curve-fitting model appears to be close to the results derived from NH (see Hinson et al. 2017 and Sect. 4) meaning that those biases remain small. We note that other authors also used stellar occultations to constrain the pressure evolution since 1988 (Young et al. 2008;Bosh et al. 2015;Olkin et al. 2015), but with less comprehensive data sets. We do not include their results here, as they were obtained with different models that might introduce systematic biases in the pressure values. Table 2 provides the pressure derived at each date, at the reference radius r = 1215 km (altitude 28 km), their scaled values at the surface using the DO15 model, as well as the pressure previously derived from the 1988 June 09 occultation. Figure 4 displays the resulting pressure evolution during the time span 1988-2016. As discussed in the previous section, even if the use of the DO15 model induces biases on p surf , it should be a good proxy for the global evolution of the atmosphere, and as such, provides relevant constraints on seasonal models.

Pressure evolution versus a volatile transport model
We interpret our occultation results in the frame of the Pluto volatile transport model developed at the Laboratoire de Météorologie Dynamique (LMD). The model is designed to simulate the volatile cycles over seasonal and astronomical timescales on the whole planetary sphere (Bertrand & Forget 2016;Forget et al. 2017;Bertrand et al. 2018). We use the latest, most realistic version of the model featuring the topography map of Pluto (Schenk et al. 2018a) and large ice reservoirs (Bertrand et al. 2018). In particular, we place permanent reservoirs of nitrogen ice in the Sputnik Planitia basin and in the depressions at mid-northern latitudes (30 • N, 60 • N), as detected by NH (Schmitt et al. 2017) and modeled in Bertrand et al. (2018). Figure 4 shows the annual evolution of surface pressure obtained with the model compared to the data. This evolution is consistent with the continuous increase of pressure observed since equinox in 1988, reaching an overall factor of almost three in 2016. This results from the progressive heating of the nitrogen ice in Sputnik Planitia and in the northern mid-latitudes, when those areas were exposed to the Sun just after the northern spring equinox in 1988, and close in time to the perihelion of 1989, as detailed in Bertrand & Forget (2016).
The model predicts that the pressure will reach its peak value and then drop in the next few years, due to: the orbitally driven decline of insolation over Sputnik Planitia and the northern mid-latitude deposits; and the fact that nitrogen condenses more intensely in the colder southern part of Sputnik Planitia, thus precipitating and hastening the pressure drop. The climate model has several free parameters: the distribution of nitrogen ice, its Bond albedo and emissivity, and the thermal inertia of the subsurface (soil). However, the large number of observation points and the recent NH observations provide strong constraints for those parameters, leading to an almost unique solution.
First, our observations restrict the possible N 2 ice surface distribution. Indeed, the southern hemisphere of Pluto is not expected to be significantly covered by nitrogen ice at the present time, because otherwise the peak of surface pressure would have occurred much earlier than 2015, as suggested by the model simulations (Fig. 4). With our model, we obtain a peak of pressure after 2015 only when considering small mid-latitudinal nitrogen deposits (or no deposit at all) in the southern hemisphere.
In our simulation, nitrogen does not condense significantly in the polar night (outside Sputnik Planitia), in spite of the length of the southern fall and winter. This is because in Pluto conditions, depending of the subsurface thermal inertia, the heat stored in the southern hemisphere during the previous southern hemisphere summer can keep the surface temperature above the nitrogen frost point throughout the cold season, or at least strongly limit the nitrogen condensation.
Consequently, the data points provide us with a second constraint, which is a relatively high subsurface thermal inertia preventing massive condensation in the southern polar night. Using a thermal inertia between 700 and 900 J s −1/2 m −2 K −1 permits us to obtain a surface pressure ratio (p surf,2015 /p surf,1988 ) of around 2.5-3, as observed. Higher (resp. lower) thermal inertia tends to lower (resp. increase) this ratio, as shown in Fig. 2a of Bertrand & Forget (2016).
Finally, the nitrogen cycle is very sensitive to the nitrogen ice Bond albedo A N2 and emissivity N2 , and only a small range for these parameters allows for a satisfactory match to the observations. Figure 4 illustrates that point. To understand it, one can do the thought experiment of imagining Pluto with a flat and isothermal surface at vapor pressure equilibrium. A rough estimate of the equilibrium temperature is provided by the classical equation: where F is the solar constant at Pluto and σ = 5.67 × 10 −8 W m −2 K 4 is the Stefan-Boltzmann constant. The surface pressure p surf is then estimated from the surface temperature T surf assuming N 2 vapor pressure equilibrium (Fray & Schmitt 2009). Consequently, the surface pressure data set inferred from stellar occultations provide us with a constraint on (1 − A N2 )/ N2 . In practice, in the model, we assume large grains for N 2 ice and we fix the emissivity at a relatively high value N2 = 0.8 (Lellouch et al. 2011). Taking F = 1.26 W m −2 (in 2015) and assuming A N2 = 0.72, we find T surf = 37.3 K, and a corresponding vapor pressure p surf = 14.8 µbar for the N 2 ice at the surface. With A N2 = 0.73, we obtain T surf = 37.0 K and p surf = 12.0 µbar. Thus, the simple equation above provides pressure values that are consistent with the volatile transport model displayed in Fig. 4. It can then be used to show that decreasing the nitrogen ice albedo by only 0.01 leads to an increase of surface pressure in 2015 by the large amount of 25%.

The 2015 June 29 occultation
The 2015 June 29 event provided seven chords across the atmosphere of Pluto; see Table A Celestial north is at the top and celestial east to the left; see labels N and E. The equator and prime meridian (facing Charon) are drawn as thicker lines. The direction of planet rotation is along the gray arrow. In the two panels, the stellar motion relative to Pluto is shown as black solid lines as seen from the Bootes-3 and Dunedin stations, with direction of motion marked by the black arrow. The shaded region at the center roughly indicates the zone where a central flash could be detected. In the upper panel, the red and blue lines are the trajectories of the primary and secondary stellar images, respectively, as seen from Bootes-3. In the lower panel, the same is shown as in the upper panel but for the stellar images as seen from Dunedin. For a spherical atmosphere, the position of the star in the sky plane, the center of Pluto and the two images are aligned, as shown in the upper panel (see the dotted line connecting the star symbols). region (Fig. 5). This was a unique opportunity to study the lower atmosphere a mere fortnight before the NH flyby ( 2015 July 14). During this short time lapse, we may assume that the atmosphere did not suffer significant global changes.
For a spherical atmosphere, there are two stellar images at any moment, a primary (near limb) image and a secondary (far limb) image that are aligned with the center of the planet and the star position, as projected in the sky plane; see Fig. 5. Since the ray-tracing code provides the refraction angle corresponding to each image, their positions along the limb of Pluto can be determined at any time (Fig. 5) and then projected onto its surface (Fig. 6).

Comparison with the REX results
The REX instrument recorded an uplinked 4.2 cm radio signal sent from Earth. The phase shift due to the neutral atmosphere The arrows indicate the direction of motion. "Ingress" (resp. "egress") refers to the disappearance (resp. re-appearance) of the images into the atmosphere of Pluto. The diamond-shaped symbols mark the positions of the image at the peak of the flash, corresponding to the time of closest approach of the respective station to the shadow center. In total, the primary image scanned longitudes from 120 to 270 • , while the secondary image scanned longitudes from 310 to 360 • and then from 0 to 70 • . The brace indicates the total duration of the primary flash (∼15 s, see Sect. 4) at Bootes-3, covering a relatively large region of more than 120 • in longitude. A similar extension applies to the secondary flash, but the brace has not been drawn for sake of clarity. The black bullets are the locations of the REX measurements at entry and exit (Hinson et al. 2017). We note the casual proximity of the REX points and the 2015 June 29 flash peaks. Right panel: same as in left panel but for the Dunedin station, where the brace has not been repeated. We note that the tracks and motions of the primary and secondary images are essentially swapped between the two stations. Notes. (a) For the ground-based observations, this is the time of closest approach to shadow center (Sicardy et al. 2016), for the REX experiment, this the beginning and end of occultation by the solid body (Hinson et al. 2017). (b) One "hour" corresponds to a rotation of Pluto of 15 • . A local time before (resp. after) 12.0 h means morning (resp. evening) limb.
was then used to retrieve the n(r), p(r), and T (r) profiles through an inversion method and the usual ideal gas and hydrostatic assumptions (Hinson et al. 2017). The REX radio occultation probed two opposite points of Pluto as the signal disappeared behind the limb (entry) and re-appeared (exit); see Fig. 6. We note that the REX entry point is at the southeast margin of Sputnik Planitia, a depression that is typically 4 km below the surrounding terrains; see Hinson et al. (2017) for details. We also note the (serendipitous) proximity of the regions scanned by the 2015 June 29 central flash and the two zones probed by REX at entry and exit. This permits relevant tests of the REX profiles against the central flash structure. The local circumstances on Pluto for the central flash and the REX occultation are summarized in Table 3. However, the local times are swapped between our observations and REX suboccultation points: the sunrise regions of one being the sunset places of the other and vice versa; see the discussion below.
The REX profiles are in good general agreement with those derived by Sicardy et al. (2016) -based itself on the DO15 procedure -between the altitudes of 5 and 115 km (Figs. 7 and 8), thus validating our approach. However, we see discrepancies at altitudes below ∼25 km (r < 1212 km), in the region where the REX entry and exit profiles diverge from one another.
Part of those differences may stem from the swapping of the sunrise and sunset limbs between the REX measurements and our observations, and from the fact that a diurnal sublimation/condensation cycle of N 2 occurs over Sputnik Planitia. In fact, lower temperatures just above the surface are expected at the end of the afternoon in that region, after an entire day of sublimation (Hinson et al. 2017). Conversely, a warmer profile could prevail at sunrise, after an entire night of condensation. Radius r (km) Fig. 7. Red and blue squares: the REX radio occultation N 2 density profiles, with the shaded area indicating the 1-σ error bar domain (Hinson et al. 2017). Below 1220 km, the errors decrease and become unnoticeable in this plot. The entry (resp. exit) profile is given from r = 1188.4 km (resp. 1193.4 km), up to 1302.4 km, where the error bars become too large for a reliable profile to be retrieved. We note that by construction, the REX entry and exit profiles are identical for r > 1220 km. Below that radius, the two profiles diverge significantly, due to different physical conditions of the boundary layer just above the surface (Fig. 8). The solid red and blue lines connecting the squares are spline interpolations of the REX profiles that are used in our ray-tracing code; see text. The REX profile is extended above r = 1302.4 km as a thin solid line, by adopting a scaled version of the 2015 June 29 profile (i.e., a mere translation of the thick solid line in this (log 10 (n), r) plot), while ensuring continuity with the REX profile. Thick solid line: the profile derived by Sicardy et al. (2016) using the DO15 light-curvefitting model. The formal 1-σ error bar of this profile is smaller than the thickness of the line, but does not account for possible biases; see text.
This warmer profile would then be more in agreement with the DO15 temperature profile. However, the difference between the REX (red) and DO15 (black) profiles in Fig. 8 remains large (more than 20 K at a given radius). This is much larger than expected from current GCMs (e.g. Forget et al. 2017, Fig. 7), which predict diurnal variations of less than 5 K at altitude levels 1-2 km above Sputnik Planitia, and less than 1 K in the ∼4-7 km region that causes the flash (Sicardy et al. 2016). In practice, Forget et al. (2017) predict that above 5 km, the temperature should be uniform over the entire planet at a given radius. This is in contrast to REX observations that reveal different temperature profiles below 25 km (Fig. 8). Thus, ingredients are still missing to fully understand REX observations, for instance the radiative impact of organic hazes, an issue that remains out of the scope of this paper.
We note that the entry REX profile goes deeper than the exit profile. This reflects the fact that the nominal radii of Pluto are at 1187.4 ± 3.6 km at entry and 1192.4 ± 3.6 km at exit (Hinson et al. 2017). This discrepancy is not significant considering the uncertainties on each radius. However, the examination of Fig. 9 shows that the most probable explanation of this mismatch is that REX probed higher terrains at exit than at entry, then providing the same pressure at a given planetocentric radius. This is the hypothesis that we adopt here, which is furthermore supported by the fact that the REX entry point is actually near the depressed region Sputnik Planitia. More precisely, the REX solution for the radius at entry (1187.4 ± 3.6 km) is fully consistent with the radius derived from NH stereo images at the same location, 1186.5 ± 1.6 km (Hinson et al. 2017). That said, we note that our data do not have enough sensitivity to constrain the absolute 35, for rays that went at about 8 km Fig. 8. As in Fig. 7 but for the temperature profiles T (r). By construction, the REX profile uses a boundary condition T b = 95.5 K at the reference radius r b = 1302.4 km in order to connect it to the DO15 profile (solid black line). Thus, the intersection of the REX and DO15 profiles at r b is a mere result of the choice of T b , and is not a measurement. There are no formal error bars on the temperature profile of Sicardy et al. (2016), as most of the errors come in this case from biases; see text. vertical scale of the density profiles at a better level than the REX solution (±3.6 km); see following section.

The 2015 June 29 central flash
The REX profiles extend from the surface (with pressures of 12.8 ± 0.7 and 10.2 ± 0.7 µbar at entry and exit, respectively) up to about 115 km, where the pressure drops to ∼1.2 µbar. Meanwhile, Sicardy et al. (2016) derive a consistent surface pressure of 12.7 µbar, with error domains that are discussed later. Nevertheless, the DO15-type thermal profile for the stratosphere (also called inversion layer) that extends between the surface and the temperature maximum at r = 1215 km is assumed to have a hyperbolic shape. The DO15 profile stops at its bottom at the point where it crosses the vapor pressure equilibrium line, thus defining the surface (assuming no troposphere). While the adopted functional form captures the gross structure of the thermal profile, it remains arbitrary. In fact, as the error bars of the REX profiles decrease with decreasing altitude, it becomes clear that the DO15 profile overestimates the temperature by tens of degrees (compared to REX) in the stratosphere as one  (Sicardy et al. 2016); see also Figs. 7 and 8. (b) As in panels a but using the nominal REX density profile. We note that the synthetic flashes are too high at both stations. (c) As in (a) and (b) but after multiplying the REX density profiles by a factor f = 0.805 and moving the shadow of Pluto 17 km north of the solution of Sicardy et al. (2016). (d) As in panels c using the nominal REX profiles, but with a topographic feature of height h = 1.35 km that blocks the stellar image during part of its motion along the southern Pluto limb (Fig. 5); the shadow has now been moved by 19.5 km north of the solution of Sicardy et al. (2016). In each panel, the value of the χ 2 function per degree of freedom (χ 2 dof ) provides an estimation of the quality of the fit; see text for discussion.
approaches the surface. Also, it ends up at the surface with a thermal gradient (16 K km −1 , see Fig. 8) that is much stronger than in the REX profiles, where it is always less that 10 K km −1 in the stratosphere. As discussed in the previous section however the N 2 diurnal cycle might induce a warmer temperature profile (after nighttime condensation) at an altitude of a few kilometres above Sputnik Planitia. This would result in a larger thermal gradient that would be closer to the DO15 profile, but still too far away from it according to GCM models, as discussed previously.
In this context, we have tested the REX profiles after modifying our ray-tracing procedure to generate new synthetic central flashes. We now account for the fact that the two stellar images that travel along the limb of Pluto probe different density profiles. To simplify the problem as much as possible, we assume that the stellar images that follow the northern and southern limbs probe an atmosphere that, respectively, has the entry and exit REX density profiles, in conformity with the geometry described in Fig. 6. This is an oversimplified approach as the stellar images actually scan relatively large portions of the limb, not just the REX entry and exit points (Fig. 6). However, this exercise allows us to assess how different density profiles may affect the shape of the central flash. To ensure smooth synthetic profiles, the discrete REX points have been interpolated by spline functions, using a vertical sampling of 25 m. Finally, above the radius r = 1302.4 km, the REX profiles have been extrapolated using a scaled version of the DO15 profile (see details in Fig. 7).
Because we want to test the shape of the central flash only, we restrict the generation of the synthetic light curves to the bottom parts of the occultation. We also include in the fit two intervals that bracket the event outside the occultation, where we know that the flux must be unity (Fig. 10). Those external parts do not discriminate the various models, but serve to properly scale the general stellar drop. Thus, the steep descents and ascents of the occultation light curves are avoided, as they would provide too much weight to the fits. Finally, since no calibrations of the light curves are available to assess the contribution φ P of Pluto to the observed flux, a linear least-square fit of the synthetic flux to the data was performed before calculating the residuals. This introduces a supplementary adjustable parameter, φ P to the fits.
Four simple scenarios are considered: (1) We first use the original model of Sicardy et al. (2016) to generate the light curves. (2) We take the REX density profiles at face value and use the modified ray-tracing model described above, fixing the ephemeris offset of Pluto as determined in Case (1). (3) We apply an adjustable, uniform scaling factor f to the two REX density profiles (which also applies to the pressure profile since the temperature is fixed), and we adjust the ephemeris offset accordingly. (4) Turning back to the REX density profiles of Case (2), we assume that a topographic feature of height h (on top of the REX exit radius, 1192.4 km) blocks the stellar image generated by the REX exit profile, that is, that the stellar image that travels along the southern limb (Fig. 5) is turned off below a planetocentric radius 1192.4 + h km.
It should be noted that the amplitude of the synthetic flash is insensitive to the absolute altitude scale that we use for the REX density profiles, to within the ±3.6 km uncertainty discussed in the previous section. For instance, displacing the REX entry profile downward by 1 km, while displacing the exit profile upward by the same amount (because the two errors and anticorrelated; see Hinson et al. 2017) changes the relative amplitude of the flash by a mere 10 −3 , well below the noise level of our observations (Fig. 10). In other words, our central flash observations cannot pin down the absolute vertical scales of the profiles to within the ±3.6 km REX uncertainty.
The fits are displayed in Fig. 10. Their qualities are estimated through the χ 2 value. Depending on the fits, there are M = 1-3 free parameters (the pressure at a prescribed level, off-track displacement of Pluto with respect to its ephemeris, and its contribution φ P to the flux). In all the fits, there are N = 217 adjusted data points. We note that the value of h in Case (4) has been fixed to 1.35 km, i.e., it is not an adjustable parameter. This is discussed further in the points below: 1. The nominal temperature profile T (r) of Sicardy et al. (2016) with surface pressure p surf = 12.7 µbar provides a satisfactory fit with χ 2 = 198 (χ 2 dof = χ 2 /(N − M) = 0.924 per degree of freedom). In this case, the Bootes-3 and Dunedin stations passed 46 km north and 45 km south of the shadow center, respectively.
2. The nominal REX profiles result in flashes that are too high compared to the observations, as noted by a visual inspection of Fig. 10 (and from χ 2 = 326, χ 2 dof = 1.52). This can be fixed by introducing haze absorption. A typical factor of 0.7 must be applied to the Bootes-3 synthetic flash in order to match the data, while a typical factor of 0.76 must be applied to the Dunedin synthetic flash. This corresponds to typical tangential optical depths (along the line of sight) in the range τ T = 0.27−0.35, for rays that went at about 8 km above the REX 1187.4 km radius. Changing the off-track offset of Pluto does not help in this case, as one synthetic flash increases while the other decreases. This could be accommodated by adjusting accordingly the optical depths τ T , but this introduces too many adjustable parameters to be relevant.
3. A satisfactory best fit is obtained (χ 2 = 214, χ 2 dof = 0.999) by uniformly reducing the REX density profiles by a factor of 0.805 and by moving the shadow center cross-track of Pluto by 17 km north with respect to Case (1), the Bootes-3, and Dunedin stations passing 29 km north and 62 km south of the shadow center, respectively. This displacement corresponds to a formal disagreement at 3-σ level for the center position of Pluto between Cases (1) and (3), when accounting for the noise present in the central flashes (Fig. 10). Thus, such a difference remains marginally significant. We also note that a satisfactory fit to the Bootes-3 flash is obtained, while the Dunedin synthetic flash remains slightly too high. As mentioned again below however, a reduction of the density profile by a factor of 0.805 is implausible considering the error bars of the REX profiles.
4. Using again the nominal REX profiles of Case (2), but imposing a topographic feature of height h = 1.35 km on top of the REX exit radius of 1192.4 km, a satisfactory fit to the Bootes-3 flash is obtained (χ 2 = 205, χ 2 dof = 0.959); in fact the best of all fits for that station. Meanwhile, the Dunedin synthetic flash remains slightly too high compared to observations. In this model, the center of the shadow of Pluto has been moved cross-track by 19.5 km north with respect to the first model, meaning that the Bootes-3 and Dunedin stations passed 26.5 km north and 64.5 km south of the shadow center, respectively. Again the discrepancy relative to the center solution of Case (1) is at 3-σ level, and is thus marginally significant. The particular choice of h = 1.35 km stems from the fact that lower values would increase the Dunedin flash even more, while higher values would decrease the Bootes-3 flash too much. We have not explored further values of h by tweaking the density profiles. Therefore, this is again an exercise to show that reasonably high topographic features may explain the observed flash. Figure 4 summarizes our results concerning the evolution of the atmospheric pressure of Pluto with time. It shows that the observed trend can be explained by adjusting the physical parameters of the planet in a rather restrictive way.

Global atmospheric evolution of Pluto
As noted in Sect. 3, this evolution is consistent with the continuous increase of pressure observed since 1988 (a factor of almost three between 1988 and 2016). It results from the heating of the nitrogen ice in Sputnik Planitia and in the northern midlatitudes, when the areas are exposed to the Sun (just after the northern spring equinox in 1989) and when Pluto is near the Sun (Bertrand & Forget 2016). The model also predicts that atmospheric pressure is expected to reach its peak and drop in the next few years, due to (1) the orbitally driven decline of insolation over Sputnik Planitia and the northern mid-latitude deposits, and (2) the fact that nitrogen condenses more intensely in the colder southern part of Sputnik Planitia, thus precipitating and hastening the pressure drop.
In that context, it is important to continue the monitoring of the atmosphere of Pluto using ground-based stellar occultations. Unfortunately, as Pluto moves away from the Galactic plane, such occultations will become increasingly rare.

Lower atmosphere of Pluto
The models presented in the Sect. 4 and illustrated in Fig. 10 are not unique and not mutually exclusive. For instance, one can have at the same time a topographic feature blocking the stellar rays, together with some haze absorption. Also, hazes, if present, will not be uniformly distributed along the limb. Similarly, topographic features will probably not be uniformly distributed along the limb, but will rather have a patchy structure that complicates our analysis. In spite of their limitations, the simple scenarios presented above teach us a few lessons: 1. Although satisfactory in terms of flash fitting, the nominal temperature profile of Sicardy et al. (2016) seems to be ruled out below the planetocentric radius ∼1215 km, since it is clearly at variance with the REX profiles (Fig. 8), while probing essentially the same zones on the surface of the dwarf planet (Fig. 6). As discussed in Sect. 4.2 however, diurnal changes occurring over Sputnik Planitia might explain this discrepancy, with a cooler (sunset) REX temperature profile and a warmer (sunrise) profile more in line with the DO15 solution. However, current GCM models predict that these diurnal changes should occur below the 5 km altitude level, and not as high as the 25 km level observed here. This issue remains an open question that would be worth investigating in future GCM models.
2. The REX profiles taken at face value cannot explain the central flashes observed at Bootes-3 and Dunedin, unless hazes are present around the ∼8 km altitude level, with optical depths along the line of sight in the range τ = 0.27-0.35. This is higher but consistent with the reported value of τ ∼ 0.24 derived from NH image analysis (Gladstone et al. 2016;Cheng et al. 2017). In fact, the two values are obtained by using quite different methods. Cheng et al. (2017) assume tholin-like optical constant, which is not guaranteed. Moreover, their 0.24 value is the scattering optical depth, while we measure the aerosol extinction (absorption plus scattering). Chromatic effects might also be considered to explain those discrepancies, as the Bootes-3, Dunedin, and the NH instruments have different spectral responses. Our data are too fragmentary however to permit such a discussion.
3. An alternative solution is to uniformly reduce the REX density profiles by a factor 0.805. However, this would induce a large disagreement (8-σ level) on the REX density profile at 7 km altitude, and thus appears to be an unrealistic scenario. Moreover, the underdense versions of the REX profiles would then formally disagree (i.e., beyond the internal error bars of the DO15 light-curve-fitting model) when extrapolated to the overlying half-light level around r = 1300 km. A remedy would be to patch up profiles derived from ground-based measurements with the underdense REX profiles, and re-run global fits. This remains out of the scope of the present analysis.
4. The topographic feature hypothesis remains an attractive alternative, as it requires modest elevation (a bit more than 1 km) above the REX exit region, which is known to be higher than the entry region, Sputnik Planitia. A more detailed examination of the elevation maps of Pluto, confronted with the stellar paths shown in Fig. 6, should be undertaken to confirm or reject that hypothesis. That said, such ± 1 km topographic variations are actually observed all over the surface of the planet (Schenk et al. 2018b).
As a final comment, we reiterate that the flashes were generated by assuming a spherical atmosphere near the surface of Pluto. There is no sign of distortion of the Bootes-3 or Dunedin flashes suggesting a departure from sphericity. It would be useful however to assess such departures, or at least establish an upper limit for them in future works. E. Meza et al.: Lower atmosphere and pressure evolution on Pluto from ground-based stellar occultations, 1988-2016