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

TL;DR

This study analyzes over a decade of stellar occultations to track Pluto's atmospheric pressure changes and constrains its lower atmosphere structure, revealing a monotonic pressure increase and insights into surface and atmospheric properties.

Contribution

It provides the first long-term calibrated record of Pluto's surface pressure evolution and constrains the lower atmosphere structure using occultation data and a central flash analysis.

Findings

Pluto's atmospheric pressure increased monotonically from 1988 to 2016.

The central flash data suggest possible large diurnal temperature variations or hazes at 4-7 km altitude.

Constraints on N2 ice albedo and emissivity are derived from pressure evolution modeling.

Abstract

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.

Tables

Table 1. Table 1: Adopted physical parameter

Pluto’s mass111111$1$Stern et al. (2015), where $G$ is the constant of gravitation.	$G{M}_P=8.696\times {10}^{11}$ m3 sec-2
Pluto’s radius111111$1$Stern et al. (2015), where $G$ is the constant of gravitation.	$R_P=1187$ km
N2 molecular mass	$\mu =4.652\times {10}^{-26}$ kg
N2 molecular	$K=1.091\times {10}^{-23}$
refractivity222222$2$Washburn (1930).	$+\left(6.282\times {10}^{-26}/{\lambda }_{\mu m}^2\right)$ cm3 molecule-1
Boltzmann constant	$k=1.380626\times {10}^{-23}$ J K-1
Pluto pole position333333$3$Tholen et al. (2008).	${\alpha }_p$= 08h 52m 12.94s
(J2000)	${\delta }_p$= -06d 10’ 04.8”

Table 2. Table 2: Pluto’s atmospheric pressure

	Surface	Pressure at	Fit quality
Date	pressure $p_{\mathrm{surf}}$	1215 km $p_{1215}$	${\chi }_{\mathrm{dof}}^2$
	($\mu$bar)	($\mu$bar)
1988 Jun 09	4.28 $\pm$ 0.44	2.33 $\pm$ 0.24111111$1$The value $p_{1215}$ is taken from Yelle & Elliot (1997). The ratio $p_{\mathrm{surf}}/p_{1215}=1.84$ of DO15’s fitting model was applied to derive $p_{\mathrm{surf}}$. Thus, the surface pressures (and their error bars) are mere scalings of the values at 1215 km. They do not account for systematic uncertainties caused by using an assumed profile (DO15 model), see discussion in subsection 3.2. The qualities of the fits (values of ${\chi }_{\mathrm{dof}}^2$) are commented on in subsection 2.3.	NA
2002 Aug 21	8.08 $\pm$ 0.18	4.42 $\pm$ 0.093	1.52
2007 Jun 14	10.29 $\pm$ 0.44	5.6 $\pm$ 0.24	1.56
2008 Jun 22	11.11 $\pm$ 0.59	6.05 $\pm$ 0.32	0.93
2008 Jun 24	10.52 $\pm$ 0.51	5.73 $\pm$ 0.21	1.15
2010 Feb 14	10.36 $\pm$ 0.4	5.64 $\pm$ 0.22	0.98
2010 Jun 04	11.24 $\pm$ 0.96	6.12 $\pm$ 0.52	1.02
2011 Jun 04	9.39 $\pm$ 0.70	5.11 $\pm$ 0.38	1.04
2012 Jul 18	11.05 $\pm$ 0.08	6.07 $\pm$ 0.044	0.61
2013 May 04	12.0 $\pm$ 0.09	6.53 $\pm$ 0.049	1.20
2015 Jun 29	12.71 $\pm$ 0.14	6.92 $\pm$ 0.076	0.84
2016 Jul 19	12.04 $\pm$ 0.41	6.61 $\pm$ 0.22	0.86

Table 3. Table 3: Regions probed by the central flash (June 29, 2015) and REX experiment (July 14, 2015)

	Time (UT)111111$1$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).	Location on surface	Local solar time222222$2$One “hour” corresponds to a rotation of Pluto of 15∘. A local time before (resp. after) 12.0 h means morning (resp. evening) limb.
June 29, 2015
Bootes-3, primary image	16:52:54.8	186.8∘E, 18.5∘S	7.67 (sunrise)
Bootes-3, secondary image	16:52:54.8	6.8∘E, 18.5∘N	19.67 (sunset)
Dunedin, primary image	16:52:56.0	8.6∘E, 19.7∘N	19.79 (sunset)
Dunedin, secondary image	16:52:56.0	188.6∘E, 19.7∘S	7.79 (sunrise)
NH radio experiment (REX), July 14, 2015
entry	12:45:15.4	193.5∘E, 17.0∘S	16.52 (sunset)
exit	12:56:29.0	15.7∘E, 15.1∘N	4.70 (sunrise)

Equations

$${ϵ}_{N2}\sigma T^4=\left(1-A_{N2}\right)\frac{F}{4},$$

Full text

*Aims. *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.

*Methods. *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 $\sim$10 $\mu$bar to 10 nbar).

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 $\sim$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 $\sim$20% and/or (d) higher terrains block part of the flash in the Charon facing hemisphere.

1 Introduction

Pluto’s tenuous atmosphere was glimpsed during a ground-based stellar occultation observed on 1985 August 19 (Brosch, 1995), and fully confirmed on 1988 June 09 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, density, see the review by Yelle & Elliot (1997).

Since then, Earth-based stellar occultations have been quite an efficient method to study Pluto’s atmosphere. It yields, in the best cases, information from a few kilometers above the surface (pressure $\sim$10 $\mu$bar) up to 380 km altitude ($\sim$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 Pluto’s atmospheric structure and evolution.

Ground-based occultations also provided a decadal monitoring of the atmosphere. Pluto has a large obliquity ($\sim 120\degr$, 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), Pluto’s heliocentric distance 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 zero degree at equinox to $54\degr$ 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 09 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 Pluto’s lower atmosphere structure derived from the flash with profiles obtained by the Radio Science EXperiment instrument on board of NH (REX hereafter) below an altitude of about 115 km

Observations, data analysis and primary results are presented in Section 2. Implications for volatile transport models are discussed in Section 3. The analysis of the 2015 June 29 central flash is detailed in Section 4, together with its consequences for Pluto’s lower atmosphere structure. Concluding remarks are provided in Section 5.

2 Observations and data analysis

2.1 Occultation campaigns

Table LABEL:tab_sites lists the circumstances of all the Pluto stellar occultation campaigns that our group have organized between 2002 and 2016. The first part of this table lists the eleven events that were used in 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-ratio 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. 2010, 2012; Benedetti-Rossi et al. 2014. Some of those campaigns are already documented and analyzed in previous publications, namely the 2002 July 20, 2002 August 21, 2007 June 14, 2008 June 22, 2012 July 18, 2013 May 04 and 2015 June 29 events. They were used to constrain Pluto’s global atmospheric structure and evolution (Sicardy et al., 2003; Dias-Oliveira et al., 2015; French et al., 2015; Olkin et al., 2015; Sicardy et al., 2016), the structure and composition (CH4, CO and HCN abundances) of the lower atmosphere by combination with spectroscopic IR and sub-mm data (Lellouch et al., 2009, 2015, 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 LABEL:tab_sites (2006 April 10). It did not provide any chord on Pluto, but was used to put an upper limit of Pluto’s rings (Boissel et al., 2014).

Note that we include here five more (yet unpublished) data sets obtained on the following dates: 2008 June 24, 2010 February 14, 2010 June 04, 2011 June 04 and 2016 July 19.

2.2 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). It 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.

Note in particular that our adopted Pluto’s radius 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\pm 4$ km. Nimmo et al. (2017) improve that value to $R_{P}=1188.3\pm 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 Section 4. Consequently, it is physically more relevant here when discussing Pluto’s lower atmospheric structure.

We assume a pure N2 atmosphere, which is justified by the fact that the next most important species (CH4) has an abundance of about 0.5% (Lellouch et al., 2009, 2015; Gladstone 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 Section 4, the tangential (line-of-sight) optical depth of hazes found by NH for the rays that graze the surface is $\tau_{T}\sim 0.24$, with a scale height of $\sim$ 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 Section 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 $\sim$35 km, see Hinson et al. (2017) and Fig. 7. This is in line with Global Climate Models (GCMs), which predict that wind velocities in the lower atmosphere should not exceed $v\sim$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 $\Delta r\sim(R_{P}v)^{2}/4GM_{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 constant. Here, the radius $r$ is counted from Pluto’s center, while Pluto’s radius found by NH is 1187 km (Table 1). This will be 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)$ constant with time is justified by the fact that the pressure is far more sensitive to Pluto’s surface temperature – through 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 N2 ice at the surface is enough to multiply the equilibrium pressure by a factor of 1.7 (Fray & Schmitt, 2009). 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 $\sim$5 km (pressure level $\sim$10 $\mu$bar) to $\sim$380 km ($\sim$10 nbar level), see DO15. Rays coming from below $\sim$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 Section 4.

Conversely, rays coming from above 380 km cause too small stellar drops ($<\sim$1%) to be of any use under usual ground-based observing conditions. This 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 $\sim$110 km and pressure $\sim$1.6 $\mu$bar).

2.3 Primary results

The ray tracing code returns the best fitting parameters, in particular the pressure at a prescribed radius (e.g. the pressure $p_{\rm surf}$ at the surface, at radius $R_{P}=1187$ km) and Pluto’s ephemeris offset perpendicular to its apparent motion, $\Delta\rho$. The ephemeris offset along the motion is treated separately, see DO15 for details. Error bars are obtained from the classical function $\chi^{2}=\sum_{1}^{N}[(\phi_{\rm i,obs}-\phi_{\rm i,syn})/\sigma_{i}]^{2}$ that reflects the noise level $\sigma_{i}$ of each of the $N$ data points, where $\phi_{\rm i,obs}$ and $\phi_{\rm i,syn}$ are the observed and synthetic fluxes, respectively. An example of $\chi^{2}(\Delta\rho,p_{\rm 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, $\chi^{2}_{\rm dof}\sim$1.06. Table 2 lists the values of $\chi^{2}_{\rm dof}$ for the other occultations, also showing satisfactory fits. Note the slightly higher values obtained for the 2002 August 21 and 2007 June 14 events (1.52 and 1.56, respectively). The presence of spikes in the light curve for the 2002 August 21 event (on top of the regular photometric noise) explains this higher value, see Fig. 2. From the same figure, we see that the 2007 June 14 light curves at Paranal were contaminated by clouds, also resulting in a slightly higher value of $\chi^{2}_{\rm dof}$. All together, those values validate a posteriori the assumptions of pure N2, transparent, spherical atmosphere with temperature profile constant in time.

In total, we collected and analyzed in a consistent manner 45 occultation light-curves obtained from eleven separate ground-based stellar occultations in the interval 2002-2016 (Table LABEL:tab_sites). The synthetic fits to the light curves are displayed in Figs 2 and 3. Fig. 11 shows the occulting chords and Pluto’s aspect for each event as seen from Earth.

Two main consequences of those results are now discussed in turn: (1) the temporal evolution of Pluto’s atmospheric pressure; (2) the structure of Pluto’s lower atmosphere using the central flash of June 29, 2015. A third product of these results is the update of Pluto’s ephemeris using the occultation geometries between 2002 and 2016. It will be presented in a separate paper (Desmars et al., in preparation).

3 Pluto’s atmospheric evolution

3.1 Constraints from occultations

In 2002, a ground-based stellar occultation revealed that Pluto’s atmospheric pressure 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. In fact, models using global volatile transport did predict this seasonal effect, among different possible scenarios (Binzel, 1990; Hansen & Paige, 1996).

Those models explored nitrogen cycles, and have been improved subsequently (Young, 2012, 2013; Hansen et al., 2015). Meanwhile, new models were developed to simulate possible scenarios for Pluto’s 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 vs. 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 Pluto’s current seasonal cycle. 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 Section 4 (Fig. 7), so that those biases remain small. 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.

3.2 Pressure evolution vs. a volatile transport model

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 subsection, even if the use of the DO15 model induces biases on $p_{\rm surf}$, it should be a good proxy for the global evolution of the atmosphere, and as such, provides relevant constrain for Pluto’s seasonal models.

We interpret our occultation results in the frame of the Pluto volatile transport model developed at the Laboratoire de Météorologie Dynamique (LMD). It is designed to simulate the volatile cycles over seasonal and astronomical times scales 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).

Fig. 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:

(1) the orbitally-driven decline of insolation over Sputnik Planitia and the northern mid-latitude deposits;

(2) 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 N2 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 little mid-latitudinal nitrogen deposits (or no deposit at all) in the southern hemisphere.

In our simulation, nitrogen does not condense much 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 so that nitrogen does not condense much in the southern polar night. Using a thermal inertia between 700-900 J s*-1/2* m*-2* K*-1* permits us to obtain a surface pressure ratio $(p_{\rm surf,2015}/p_{\rm surf,1988})$ of around 2.5-3, as observed. Higher (resp. lower) thermal inertia tend 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 $\epsilon_{N2}$, and only a small range for these parameters allows for a satisfactory match to the observations. Fig. 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:

[TABLE]

where $F$ is the solar constant at Pluto and $\sigma=5.67\times 10^{-8}$ W m*-2* K4 is the Stefan-Boltzmann constant. The surface pressure $p_{\rm surf}$ is then estimated from the surface temperature $T_{\rm surf}$ assuming N2 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})/\epsilon_{N2}$. In practice, in the model, we assume large grains for N2 ice and we fix the emissivity at a relatively high value $\epsilon_{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_{\rm surf}=37.3$ K, and a corresponding vapor pressure $p_{\rm surf}=$ 14.8 $\mu$bar for the N2 ice at the surface. With $A_{N2}=0.73$, we obtain $T_{\rm surf}=37.0$ K and $p_{\rm surf}=12.0$ $\mu$bar. Thus, the simple equation above provides pressure values that are consistent with the volatile transport model displayed in Fig. 4. It then can be used to show that decreasing the nitrogen ice albedo by only 0.01 leads to an increase of surface pressure in 2015 by a large amount of 25%.

4 Pluto’s lower atmosphere

4.1 The June 29, 2015 occultation

The June 29, 2015 event provided seven chords across Pluto’s atmosphere, see Table LABEL:tab_sites and Fig. 11. A first analysis of this event is presented in Sicardy et al. (2016). The two southernmost stations (Bootes-3 and Dunedin) probed the central flash region (Fig. 5). This was a unique opportunity to study Pluto’s lower atmosphere a mere fortnight before the NH flyby (July 14, 2015). During this short time lapse, we may assume that the atmosphere did not suffer significant global changes.

For a spherical atmosphere, there are at any moment two stellar images, a primary (near limb) image and a secondary (far limb) image that are aligned with Pluto’s center 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 Pluto’s limb can be determined at any time (Fig. 5), and then projected onto Pluto’s surface (Fig. 6).

4.2 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 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. 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.

Note also the (serendipitous) proximity of the regions scanned by the June 29, 2015 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, that 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 km and 115 km (Figs. 7 and 8), thus validating our approach. However, we see discrepancies at altitudes below $\sim$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 to the fact that a diurnal sublimation/condensation cycle of N2 occurs over Sputnik Planitia. Then, 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. 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 $\sim$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.

Note that the entry REX profile goes deeper than the exit profile. This reflects the fact that the nominal Pluto’s radii are at $1187.4\pm 3.6$ km at entry and $1192.4\pm 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 will 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\pm 3.6$ km) is fully consistent with the radius derived from NH stereo images at the same location, $1186.5\pm 1.6$ km (Hinson et al., 2017). This said, note that our data do not have enough sensitivity to constrain the absolute vertical scale of the density profiles at a better level than the REX solution ($\pm 3.6$ km), see next subsection.

4.3 The June 29, 2015 central flash

The REX profiles extend from the surface (with pressures of $12.8\pm 0.7$ and $10.2\pm 0.7$ $\mu$bar at entry and exit, respectively) up to about 115 km, where the pressure drops to $\sim$1.2 $\mu$bar. Meanwhile, Sicardy et al. (2016) derive a consistent surface pressure of 12.7 $\mu$bar, with error domains that are discussed later.

This said, 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 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 subsection, however, the N2 diurnal cycle might induce a warmer temperature profile (after nighttime condensation) at a few km altitude 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 that 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 Pluto’s limb probe different density profiles. To simplify as much as possible the problem, 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 rather 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 meters. 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 scale properly 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 Pluto’s contribution $\phi_{P}$ to the observed flux, a linear least-square fit of the synthetic flux to the data has been performed before calculating the residuals. This introduces a supplementary adjustable parameter, $\phi_{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 Pluto’s ephemeris offset as determined in Case (1). (3) We apply an adjustable, uniform scaling factor $f$ to the two REX density profiles (which thus also applies to the pressure profile since the temperature is fixed), and we adjust Pluto’s 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, i.e. 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 $\pm 3.6$ km uncertainty discussed in the previous subsection. 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 $\pm 3.6$ km REX uncertainty.

The fits are displayed in Fig. 10. Their qualities are estimated through the $\chi^{2}$ value. Depending on the fits, there are $M=1$ to 3 free parameters (the pressure at a prescribed level, off-track displacement of Pluto with respect to its ephemeris and Pluto’s contribution $\phi_{P}$ to the flux). In all the fits, there are $N=217$ data points adjusted. Note that the value of $h$ in Case (4) has been fixed to 1.35 km, i.e. is not an adjustable parameter. This is discussed further in the points below:

The nominal temperature profile $T(r)$ of Sicardy et al. (2016) with surface pressure $p_{\rm surf}=12.7$ $\mu$bar provides a satisfactory fit with $\chi^{2}=198$ ($\chi^{2}_{\rm dof}=\chi^{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. 2.

The nominal REX profiles result in flashes that are too high compared to the observations, as noted by a visual inspection of the figure (and from $\chi^{2}=326$, $\chi^{2}_{\rm 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 $\tau_{T}=0.27-0.35$, for rays that went at about 8 km above the REX 1187.4 km radius. Changing Pluto’s off-track offset 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 $\tau_{T}$, but this introduces too many adjustable parameters to be relevant. 3. 3.

A satisfactory best fit is obtained ($\chi^{2}=214$, $\chi^{2}_{\rm dof}=0.999$) by reducing uniformly the REX density profiles by a factor of 0.805 and by moving Pluto’s shadow center cross-track 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-$\sigma$ level for Pluto’s center position between Case (1) and (3), when accounting for the noise present in the central flashes (Fig. 10). Thus, such difference remains marginally significant. Note also that a satisfactory fit to the Bootes-3 flash is obtained, while the Dunedin synthetic flash remains a bit too high. As commented in the concluding Section, however, a reduction of the density profile by a factor of 0.805 is implausible considering the error bars of the REX profiles. 4. 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 ($\chi^{2}=205$, $\chi^{2}_{\rm dof}=0.959$), in fact the best of all fits for that station. Meanwhile, the Dunedin synthetic flash remains a bit too high compared to observations. In this model, Pluto’s shadow center has been moved cross-track by 19.5 km north with respect to the first model, so 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 Pluto’s center solution of Case (1) is at 3-$\sigma$ level, and thus marginally significant. The particular choice of $h=1.35$ km stems from the fact that lower values would increase even more the Dunedin flash, while higher values would decrease too much the Bootes-3 flash. We have not explored further other values of $h$ by tweaking the density profiles. So, this is again an exercise to show that reasonably high topographic features may explain the observed flash.

5 Concluding remarks

5.1 Pluto’s global atmospheric evolution

Fig. 4 summarizes our results concerning the evolution of Pluto’s atmospheric pressure with time. It shows that the observed trend can be explained by adjusting Pluto’s physical parameters in a rather restrictive way.

As noted in Section 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 mid-latitudes, 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 Pluto’s atmosphere using ground-based stellar occultations. Unfortunately, as Pluto moves away from the Galactic plane, such occultations will become rarer and rarer.

5.2 Pluto’s lower atmosphere

The models presented in the Section 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 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 $\sim 1215$ km, since it is clearly at variance with the REX profiles (Fig. 8), while probing essentially the same zones on Pluto’s surface (Fig. 6). As discussed in Section 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 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 $\sim 8$ km altitude level, with optical depths along the line of sight in the range $\tau=$ 0.27-0.35. This is higher but consistent with the reported value of $\tau\sim 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, though, to permit such a discussion.

(3) An alternative solution is to reduce uniformly the REX density profiles by a factor 0.805. However, this would induce a large disagreement (8-$\sigma$ 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 disagree formally (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 ground-based-derived profiles 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, that is known to be higher than the entry region, Sputnik Planitia. A more detailed examination of Pluto’s elevation maps, confronted with the stellar paths shown in Fig. 6, should be undertaken to confirm or reject that hypothesis. This said, such $\pm$ 1 km topographic variations are actually observed all over Pluto’s surface (Schenk et al., 2018b).

As a final comment, we recall that the flashes have been generated by assuming a spherical atmosphere near Pluto’s surface. There is no sign of distortion of the Bootes-3 and Dunedin flashes that suggests 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.

Acknowledgements.

This article is dedicated to the memory of H.-J. Bode, J. G. Greenhill and O. Faragó for their long-standing support and participation to occultation campaigns. The work leading to these results has received funding from the European Research Council under the European Community’s H2020 2014-2020 ERC Grant Agreement n∘ 669416 “Lucky Star”. EM thanks support from Concytec-Fondecyt-PE and GA, FC-UNI for providing support during the 2012 July 18 occultation. BS thanks S. Para for partly supporting this research though a donation, J. P. Beaulieu for helping us accessing to the Hobart Observatory facilities and B. Warner, B. L. Gary, C. Erickson, H. Reitsema, L. Albert, P. J. Merritt, T. Hall, W. J. Romanishin, Y. J. Choi for providing data during the 2007 March 18 occultation. MA thanks CNPq (Grants 427700/2018-3, 310683/2017-3 and 473002/2013-2) and FAPERJ (Grant E-26/111.488/2013). JLO thanks support from grant AYA2017-89637-R. PSS acknowledges financial support from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement no 687378, as part of the project “Small Bodies Near and Far” (SBNAF). JLO, RD, PSS and NM acknowledge financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). FBR acknowledges CNPq support process 309578/2017-5. GBR thanks support from the grant CAPES-FAPERJ/PAPDRJ (E26/203.173/2016). JIBC acknowledges CNPq grant 308150/2016-3. RVM thanks the grants: CNPq-304544/2017-5, 401903/2016-8, and Faperj: PAPDRJ-45/2013 and E-26/203.026/2015. BM thanks the CAPES/Cofecub-394/2016-05 grant and CAPES/Brazil - Finance Code 001. BM and ARGJ were financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. TRAPPIST-North is a project funded by the University of Liège, in collaboration with Cadi Ayyad University of Marrakech (Morocco). TRAPPIST-South is a project funded by the Belgian Fonds (National) de la Recherche Scientifique (F.R.S.-FNRS) under grant FRFC 2.5.594.09.F, with the participation of the Swiss National Science Foundation (FNS/SNSF). VSD, SPL, TRM and ULTRACAM are all supported by the STFC. KG acknowledges help from the team of Archenhold-Observatory, Berlin, and AR thanks G. Román (Granada) for help during the observation of the 2016 July 19 occultation. AJCT acknowledges support from the Spanish Ministry Project AYA2015-71718-R (including EU funds). We thank Caisey Harlingten for the repeated use of his 50 cm telescopes in San Pedro de Atacama, Chile. We thank the Italian Telescopio Nazionale Galileo (TNG), operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. LM acknowledges support from the Italian Minister of Instruction, University and Research (MIUR) through FFABR 2017 fund and support from the University of Rome Tor Vergata through “Mission: Sustainability 2016” fund. The Astronomical Observatory of the Autonomous Region of the Aosta Valley (OAVdA) is managed by the Fondazione Clément Fillietroz-ONLUS, which is supported by the Regional Government of the Aosta Valley, the Town Municipality of Nus and the “Unité des Communes valdôtaines Mont-Émilius”. The research was partially funded by a 2016 “Research and Education” grant from Fondazione CRT. We thank D.P. Hinson for his constructive and detailed comments that helped to improve this article.

6 Circumstances of Observations

Appendix A Reconstructed geometries of the occultations