Pluto's Atmosphere in Plateau Phase Since 2015 from a Stellar Occultation at Devasthal

Owing to its high obliquity (120°) and high orbital eccentricity (0.25), Pluto suffers intense seasonal episodes. Its poles remain, for decades, in permanent sunlight or darkness over its 248 yr heliocentric revolution. This leads to strong effects on its N₂ atmosphere which is mainly controlled by vapor pressure equilibrium with the surface N₂ ice. The NASA New Horizons flyby in 2015 July revealed a large depression, Sputnik Planitia, filled by N₂ ice (Stern et al. 2015), which appears to be the main engine that controls the seasonal variation of atmospheric pressure during one seasonal cycle (Bertrand & Forget 2016; Bertrand et al. 2018; Johnson et al. 2021). Apart from these crucial results, a comprehensive review of the composition, photochemistry, atmospheric dynamics, circulation, and escape processes derived from the New Horizons data are presented in Gladstone & Young (2019).

In parallel, ground-based observations of stellar occultations allowed various teams to accurately monitor Pluto's atmospheric pressure since 1988. A compilation of twelve occultations observed between 1988 and 2016 shows a three-fold monotonic increase of pressure during this period that can be explained by the progression of summer over the northern hemisphere, exposing Sputnik Planitia to solar radiation (Meza et al. 2019). This increase can be explained consistently by a Pluto volatile transport model, which predicts that the pressure should peak around 2020. A gradual decline should then last for two centuries under the combined effects of Pluto's recession from the Sun and the prevalence of the winter season over Sputnik Planitia.

Here we present the results of a stellar occultation by Pluto that occurred on 2020 June 6. It was observed in the near-infrared by two large telescopes at the Devasthal station, Nainital, India. The high signal-to-noise ratio light curves obtained with these instruments allow us to derive an accurate value of Pluto's atmospheric pressure using the same approach as in Dias-Oliveira et al. (2015) and Meza et al. (2019). This occultation was particularly timely as it can test the validity of the current models of Pluto's atmosphere evolution. Moreover, as Pluto is now moving away from the Galactic plane as seen from Earth, stellar occultations by the dwarf planet are becoming increasingly rare, making this event a decisive one.

2.1. Occultation

The 2020 June 6 occultation campaign was organized within the Lucky Star project. ¹⁵ The prediction used the Gaia DR2 position at the epoch of occultation (Table 1) and the NIMAv8/PLU055 Pluto's ephemeris derived from previous occultations observed since 1988 (Desmars et al. 2019). More information on the event (shadow path, charts, photometry, etc.) is available from a dedicated web page. ¹⁶

Table 1. Circumstances of Observations, Adopted Parameters, and Result of the Atmospheric Fit

Observation Log
3.6 m (DOT)
Coordinates, altitude	79° $41^{\prime}$ 36 E, 29° $21^{\prime}$ 394 N, 2450 m
Camera	TIRCAM2 Raytheon InSb array
Filter (λ/Δλ, μm)	H (1.60/0.30)
Exposure time/Cycle time (s)	5./5.336
1.3 m (DFOT)
Coordinates, altitude	79° $41^{\prime}$ 61 E, 29° $21^{\prime}$ 415 N, 2450 m
Camera	ANDOR DZ436
Filter (λ/Δλ, μm)	I (0.85/0.15)
Exposure time/Cycle time (s)	1.7/2.507
Occulted star
Identification (Gaia DR2)	6864932072159710592
J2000 position at epoch (ICRF)	α = 19h45m339079, $\delta =-22^\circ 10^{\prime} 19\buildrel{\prime\prime}\over{.} 128$
Pluto's body
Mass a	GMP = 8.696 × 1011 m3 s−2
Radius a	RP =1187 km
Geocentric distance	4.97407 × 109 km
Pluto's atmosphere
N2molecular mass	μ = 4.652 × 10−26 kg
N2 molecular	K = 1.091 × 10−23
Refractivity b	$+6.282\times {10}^{-26}/{\lambda }_{\mu {\rm{m}}}^{2})$ cm3 molecule−1
Boltzmann constant	k = 1.380626 × 10−23 J K−1
Results of atmospheric fit (with 1σ error bars)
Pressure at radius 1215 km	${p}_{1215}={6.655}_{-0.21}^{+0.35}$ μbar
Surface pressure c	${p}_{\mathrm{surf}}={12.23}_{-0.38}^{+0.65}$ μbar
Closest approach distance of Devasthal to shadow center d	${\rho }_{{\rm{C}}/{\rm{A}},\ {\rm{D}}}=-{735}_{-15}^{+7}$ km
Closest approach time of Devasthal to shadow center	tC/A, D = 19:02:43.0 ± 0.14 s UT
Geocentric closest approach distance to shadow center d	${\rho }_{{\rm{C}}/{\rm{A}},\ {\rm{G}}}=+{6044}_{-7}^{+15}$ km
Geocentric closest approach time to shadow center e	tC/A, G = 19:01:01.7 ± 0.14 s UT

Notes.

^a Stern et al. (2015), where G is the constant of gravitation. ^b Washburn (1930), where λ_μ m is the wavelength expressed in microns. ^c Using a ratio p_surf/p₁₂₁₅ = 1.837 given by the template model Meza et al. (2019). ^d Negative (resp. positive) values mean that the point considered went south (resp. north) of the shadow center. ^e Although the quoted error bar is small, a systematic error of about 1 s may be present in t_{C/A, G}; see text.

Download table as:  ASCIITypeset image

The event was successfully recorded in the I and H bands using the 1.3 m Devasthal Fast Optical Telescope (DFOT) and the 3.6 m Devasthal Optical Telescope (DOT), respectively. The I- and H-band magnitudes of the occulted star are ∼12.3 and 11.6, respectively, while those of Pluto during the epoch of occultation were ∼13.8 and 13.3, respectively. Thus, the occulted star was more than 1.5 mag brighter than the combined Pluto + Charon system, ensuring a high-contrast event, i.e., a significant drop of the total flux, which combines the fluxes from the star, Pluto, and Charon.

Both telescopes are operated by the Aryabhatta Research Institute of Observational Sciences (ARIES) located at Nainital, India. Observations were also planned with the 2 m Himalayan Chandra Telescope (HCT), Hanle, operated by the Indian Institute of Astrophysics, Bangalore, India. However, the event was clouded out at that site.

The event was observed in the I band with the DFOT-Andor DZ436 camera (2048 × 2048-pixel; plate scale ∼05 pixel⁻¹). The central 401 × 401 pixel region was used in 2 × 2 binning mode. With a readout rate of 1 MHz, shift speed of 16 μs, and exposure time of 1.7 s, a total cycle time of 2.507 s was achieved. The final acquired image was a FITS data cube of 600 frames. A simultaneous observation was carried out with DOT in the H band with the TIRCAM2 instrument (Naik et al. 2012; Baug et al. 2018), an imaging camera housing a cryo-cooled Raytheon InSb Aladdin III Quadrant focal plane infrared array (512 × 512 pixel; plate scale=0167 pixel⁻¹). The full frame was used with a readout rate of 1 MHz, exposure time of 5 s, and total cycle time of 5.336 s. The final acquired image was a FITS data cube of 255 frames.

From the light-curve fitting described below, we reconstructed Pluto's shadow path on Earth and the geometry of the occultation (Figure 1). Note in this figure that two stellar images (primary and secondary, see Sicardy et al. 2016) actually scanned Pluto's limb. However, the flux of the secondary image was always fainter than that of the primary by a factor larger than 25, making it negligible in our case. Consequently, this event essentially scanned the northern summer hemisphere of Pluto.

Zoom In Zoom Out Reset image size

2.2. Light-curve Fitting

The observed data sets were fitted with synthetic light curves using the method described in Dias-Oliveira et al. (2015), in particular with the same template temperature profile, T(r), where r is the distance to Pluto's center. The approach involves the simultaneous fitting of the observed refractive occultation light curves by synthetic profiles that are generated by a ray-tracing code that uses the Snell–Descartes law. The various parameters used in our fitting procedure are listed in Table 1.

There are M = 4 adjusted parameters in our model: (1) p_surf, the pressure at Pluto's surface, (2) Δρ, the cross track offset to Pluto's ephemeris, and (3–4) the two Pluto+Charon contributions ϕ₀ to each of the two ARIES light curves. Owing to less than optimal sky conditions prevailing, observations to separately measure the occulted star and Pluto's system the nights prior to or after the event were not possible. Due to this unavailability of calibration data, the ϕ₀s in our analysis are not known.

There is a fifth adjusted parameter that is completely uncorrelated with the other four, the time shift, Δt, which needs to be applied to Pluto's ephemeris to best fit the data. This parameter accounts for the ephemeris offset to apply along Pluto's apparent motion and for errors in the star position. It finally provides, t_{C/A, G}, the time of closest approach of Pluto to the star in the sky plane, as seen from the Geocenter; see Table 1.

When fitting the two Devasthal light curves, a discrepancy of 2.4 s appeared in the best-fitting Δt derived from each telescope, thus revealing a problem in the recording of the absolute times at one (or both) telescopes. Since it is difficult to decipher the origin of this discrepancy, any attempt to correct for the same would be futile. Hence, we have chosen to apply independent Δt to each telescope, and calculate the final t_{C/A, G} as an average of the two obtained values, weighted by the quality of the two light curves (measured by the noise in the data). With this approach, although the internal error bar on t_{C/A, G} is small (±0.14 s, Table 1), a systematic uncertainty of the order of 1 s still remains in the quoted value of t_{C/A, G}.

The best fit to the data are shown in Figure 2. The function ${\chi }^{2}\,=\,{\sum }_{1}^{N}\,{[({\phi }_{{\rm{i}},\mathrm{obs}}-{\phi }_{{\rm{i}},\mathrm{syn}})/{\sigma }_{i}]}^{2}$ was used to assess the quality of the fit, where σ_i reflects the noise level of each of the N = 362 data points, and ϕ_i,obs and ϕ_i,syn are the observed and synthetic fluxes at the i^th data point, respectively. Satisfactory fits are obtained for a χ² value per degree of freedom ${\chi }_{\mathrm{dof}}^{2}={\chi }_{\min }^{2}/(N-M)\sim 1$, where ${\chi }_{\min }^{2}$ is the minimum value of χ² obtained in the fitting procedure. This is the case here, with individual values ${\chi }_{\mathrm{dof}}^{2}=0.984$ and ${\chi }_{\mathrm{dof}}^{2}=1.07$ at DOT and DFOT, respectively, and a global value ${\chi }_{\mathrm{dof}}^{2}=1.03$ corresponding to the simultaneous fit to the two light curves.

Zoom In Zoom Out Reset image size

Our fit is mainly sensitive to regions above 30 km altitude, so our primary result is the pressure near the radius 1215 km, ${p}_{1215}={6.665}_{-0.21}^{+0.35}$ μbar (Table 1). A factor of 1.837 is then applied to convert this into p_surf; see Meza et al. (2019). Figure 3 shows the χ² map plotted as a function of the two parameters p_surf and Δρ. Note that because there is only one occultation chord at hand (Figure 1), a correlation between p_surf and Δρ is observed. Using the ${\chi }_{\min }^{2}+1$ criterion, we obtain the best-fitting value of ${p}_{\mathrm{surf}}={12.23}_{-0.38}^{+0.65}$ μbar. The marginal error bar quoted here is estimated by ignoring the value of Δρ.

Zoom In Zoom Out Reset image size

Besides the observations presented here, the 2020 June 6 event was observed by an independent team; see Poro et al. 2021. The event was observed at low elevation (6°) near the city of Karaj in Iran with a 60 cm telescope. The authors mention that they use the ray-tracing method of Dias-Oliveira et al. (2015), i.e., a procedure that should be fully consistent with our own approach.

However, their derived value p_surf = 12.36 ± 0.38 μbar, is questionable. First, the time axis for the Karaj light curve in Poro et al. (2021) is wrong by a large factor of about two; see Figure 2. This makes it impossible to obtain any realistic value of p₁₂₁₅ from this light curve. Second, assuming that the Karaj time axis has some (undocumented) problems, considering the occultation geometry at that station and adopting the p₁₂₁₅ = 6.72 μbar value of Poro et al. (2021), we generated a synthetic model using our own ray-tracing code. By shrinking the timescale in an attempt to superimpose our model (green curve in Figure 2) onto the results of Poro et al. (2021; black curve), we see a clear discrepancy between our models in the deepest part of the occultation. This shows that the ray-tracing code of Poro et al. (2021) is inconsistent with ours. Consequently, the results of Poro et al. (2021) are impossible to obtain considering their published light curve, and probably stems from improper use of their ray-tracing code. Finally, we note that the error bar for p_surf is inconsistent with the error bar that the authors obtain for the pressure at radius 1215 km, p₁₂₁₅ = 6.72 ± 0.48 μbar, see their Figure 3. As the error bar scales like the value of the pressure, the error bar for the surface pressure should be ±0.88 μbar, not ±0.38 μbar.

In Figure 3, we plot our measurement of p_surf in 2020 (red point) along with other published values (Hinson et al. 2017; Meza et al. 2019; Arimatsu et al. 2020; Young et al. 2021). The 2020 June 6 occultation shows that the pressure increase prevailing between 1988 and 2013 stopped and has reached a stationary regime since 2015. This is in line with the Pluto volatile transport model described in Meza et al. (2019).

Our results do not support the rapid pressure decrease claimed by Arimatsu et al. (2020) (who also used the ray-tracing method of Dias-Oliveira et al. 2015) from an occultation observed on 2019 July 17; see the point "A20" in Figure 3. With the closest approach to Pluto's shadow center of 1008 km, the work of Arimatsu et al. (2020) is based on an occultation that was more grazing than the event reported here (with closest approach of 735 km). This induces a larger correlation between the parameters p_surf and Δρ. In particular, Arimatsu et al. (2020) mention that the pressure drop they find between 2016 and 2019 is actually detected at the 2.4σ level, and thus remains marginally significant. We thus estimate that the 2019 data point lacks accuracy to claim that a large decrease occurred in 2019, followed by a return to a pressure close to that of 2015 during the year 2020 (this work).

We do not confirm either the decrease of pressure reported by Young et al. (2021), based on the 2018 August 15 occultation observed from about two dozen stations in USA and Mexico, from which Young et al. (2021) derive a value p_surf = 11.13 ± 0.4 μbar (the point "Y21" in Figure 3), using a method that is not described by these authors.

It is important to note that all the points derived by us between 2002 and 2020 (Meza et al. 2019 and present work) are obtained using a unique template temperature profile T(r). This assumption is backed up by the fact that although the pressure increased by a factor of about three between 1988 and 2016 (Meza et al. 2019), the retrieved temperature profiles in 1988, 2002, 2012, and 2015 (Yelle & Elliot 1997; Hinson et al. 2017; Sicardy et al. 2003; Young et al. 2008; Dias-Oliveira et al. 2015, respectively) are all similar, with a strong positive thermal gradient in the lower part of the atmosphere that peaks at T ∼ 110 K near r = 1215 km, followed by a roughly isothermal upper branch with a mild negative thermal gradient. This globally fits the methane-thermostat model of Yelle & Lunine (1989), where the upper-atmosphere temperature is robustly maintained near 100 K through the radiative properties of atmospheric CH4, almost independently of its abundance, while the lower part is forced by heat conduction with the cold surface near 38 K. So even if small systematic errors are introduced at this stage, a consistent comparison using a constant T(r) between events is possible, so that general trends on pressure evolution can be monitored. Thus, at this stage, the methodologies used by both Arimatsu et al. (2020) and Young et al. (2021) should be compared with ours in some well-defined test cases to see if our approaches are consistent, and thus, fully comparable.

Note that our results could be compared with those of New Horizons, derived from the radio occultation experiment (REX) in 2015 July. Values of p_surf = 12.8 ± 0.7 and 10.2 ± 0.7 μbar at entry and exit, respectively, were obtained (Hinson et al. 2017). The difference between the two values can be attributed to a 5 km difference in radius at the two locations probed by REX, the lower value at exit corresponding to higher terrains on Pluto. As the entry value probed a point over Sputnik Planitia, where sublimation of N₂ takes place, it should be more representative of p_surf than the exit value. We see that our value of p_surf = 12.71 ± 0.14 μbar derived from the 2015 June 29 occultation (Meza et al. 2019) is in excellent agreement with the REX value of 2015 July (point "H17" in Figure 3), giving confidence that our approach not only provides general trends but also good estimates of p_surf.

The 2020 June 6 stellar occultation allowed us to constrain Pluto's atmospheric evolution. The surface pressure that we obtain, ${p}_{\mathrm{surf}}={12.23}_{-0.38}^{+0.65}$ μbar, shows that Pluto's atmosphere has reached a plateau since mid-2015, a result which is in line with the Pluto volatile transport model discussed in Meza et al. (2019). Our result does not support the drops of pressure reported by Young et al. (2021) and Arimatsu et al. (2020) in 2018 and 2019, respectively. These inconsistencies call for careful comparisons between methodologies before any conclusions based on independent teams can be drawn.

We note that if the model presented in Meza et al. (2019) is correct, and considering the typical error bars derived from occultations, it will be difficult to firmly confirm a pressure drop before 2025. Meanwhile, observations should be organized whenever possible, as unaccounted processes may cause pressure changes not predicted by models.

The authors are grateful to the Directors of ARIES and IIA for granting time under Directors' Discretionary Time allotment for observing this event. The authors also thank the staff of the Devasthal Observatory and the IR Group at TIFR for their help during observations. The work leading to these results has received funding from the European Research Council under the European Community's H2020 2014-2021 ERC Grant Agreement no. 669416 "Lucky Star." Research at PRL and IIST is supported by the Department of Space, Government of India. M.A., G.B.-R., F.B.-R., and R.V.-M. thank CNPq and CAPES for grants 313144/2020-6, 314772/2020-0, 465376/2014-2, and Process 88887.310463/2018-00—Project 88887.571156/2020-00. R.V.-M. also thanks grant CNPq 304544/2017-6. P.S-S. acknowledges financial support by the Spanish grant AYA-RTI2018-098657-J-I00 "LEO-SBNAF" (MCIU/AEI/FEDER, UE).

Facilities: Devasthal: 3.6 m - , 1.3 m - .