© J. Marques Oliveira et al. 2022

1 Introduction

The large satellite Triton was discovered in 1846, only 17 days after the discovery of its planet, Neptune. An atmosphere was first speculated by Cruikshank & Silvaggio (1979), with the claimed detection of the gaseous CH₄ spectral signature, although in retrospect the features were due to methane ice on the surface. In any case, the presence of this volatile ice did suggest the existence of an atmosphere. This was to be confirmed ten years later, when the NASA Voyager 2 (V2) spacecraft flew by the Neptunian system in August 1989. During this flyby, Triton’s tenuous atmosphere (mainly nitrogen N2) was detected during the Radio Science Subsystem (RSS) occultation, providing its surface density, pressure, and temperature (Tyler et al. 1989). These results were later improved by Gurrola (1995, G95 hereafter), who derived a surface pressure of p_suf = 14 ± 2 µbar. The studies that used RSS data did not, however, provide a thermal profile. Dynamical aspects, such as wind regimes, were studied using V2 images of plumes near Triton’s surface (Yelle et al. 1991), as well as vertical profiles of methane and hazes using V2’s UV images (Strobel et al. 1990; Herbert & Sandel 1991; Krasnopolsky et al. 1993; Krasnopolsky 1993; Krasnopolsky & Cruikshank 1995; Strobel & Summers 1995).

Triton is currently experiencing a rare ‘extreme southern solstice’, a configuration that occurs every ~650 years (Fig. 1). In particular, the sub-solar latitude on the satellite reached about 50° S in 2000. The various measurements of Triton’s atmospheric pressure using occultations bracket that epoch, from the RSS results in 1989 to the ground-based stellar occultation of 2017 discussed here. In this context, it is interesting to look for ongoing seasonal effects (if any) occurring in Triton’s atmosphere, especially large pressure variations in the last three decades. Such seasonal variations (or its absence) can then constrain global climate models (GCMs) and volatile transport models (VTMs) that account for volatile transport induced by insolation changes.

Since 1989, only a handful of Earth-based stellar occultations have been observed, and there have been a few spectroscopic studies that detected CO and CH₄ in the near-IR (Lellouch et al. 2010) and CO and HCN with the Atacama Large Millimeter/submillimeter Array (ALMA) in the millimetre (Gurwell et al. 2019). As discussed later, the deepest layers accessible during Triton Earth-based occultations (when a central flash¹ is observed) typically lie at 8 km altitude (~9 µbar pressure level). Conversely, the best light curves can provide information up to an altitude of about 190 km, corresponding to a few nanobars.

Here, we report on results obtained from the 5 October 2017 ground-based stellar occultation. Attempts to view this rare event were made from more than 100 sites in Europe, northern Africa, and the eastern USA. We extracted 90 occultation light curves from this campaign, among which 42 show a central flash. This is by far the most observed stellar occultation by Triton ever monitored (and among the most observed event of its kind, all Solar System objects combined) both in terms of latitudinal coverage of the satellite and central flash sampling.

Our goals are to (1) provide Triton’s atmospheric profiles (density, pressure, temperature) derived from this event, (2) compare the results with those obtained from previous occultations, including the V2 RSS experiment, (3) constrain the seasonal variations of Triton’s atmosphere, (4) compare the results with current GCMs, and (5) and derive the shape of the central flash layer.

The 5 October 2017 event is discussed in Sect. 2. The methods used are given in Sect. 3, and in Sect. 4 we present the results. We reanalyse previous events, including a new approach to retrieve new information from the V2 experiment in Sect. 5. In Sect. 6 we discuss the atmospheric seasonal variations, using pressure values from our work and from other works. Section 7 focuses on the analysis of the central flashes. We mention some issues that are to be addressed at a future date in Sect. 8. Concluding remarks are provided in Sect. 9.

2 The 5 October 2017 stellar occultation

2.1 Prediction

In the past two decades, Neptune has been crossing regions with low surface-density of stars. In the early 2010s, faint star surveys up to R = 19 were made using the Wide Field Imager attached to the 2.2 m Max-Planck telescope at the European Southern Observatory (ESO; Assafin et al. 2010, 2012; Camargo et al. 2014). They provided many candidates for occultations by Pluto, large Trans-Neptunian Objects and Centaurs, but no suitable Triton events for 2008-2015.

Predicting occultations by Triton is problematic for two reasons: (1) Neptune’s orbit may have systematic errors, causing a systematic shift of Triton’s position with respect to the stars; (2) Triton’s neptunocentric orbit may have systematic errors due to the large brightness and colour differences between Neptune and Triton, and also from changes in Neptune’s magnitude, making their relative colour variable (Schmude et al. 2016). Both points affect differently the ground-based measurements of Triton and Neptune, resulting in a distorted neptunocentric orbit for the satellite. A way to overcome these problems is to use R or I filters to minimise differential refraction during observations and distribute Triton observations evenly along its orbit around Neptune. Triton’s path is then set by the average ephemeris offsets found for right ascension (α) and declination (δ).

In this context, we gathered more than 4700 charge-coupled-device (CCD) images of Triton in R between 1992 and 2016 with the 0.6-m B&C and 1.6-m P&E telescopes at Pico dos Dias Observatory (OPD) in Brazil (IAU code 874). The observations were reduced with the best available astrometric catalogue at that time, the UCAC4 (Zacharias et al. 2013), following the same procedures described in Gomes-Júnior et al. (2015). From the average ephemeris offsets of many nights, and after sigma-clipping, we found an overall offset (Δα cos δ, Δδ) = (+1 ± 45, − 16 ± 45) milliarcsecond (mas), with error bars at 1σ level, with respect to the DE435/NEP081 ephemeris from the Jet Propulsion Laboratory (JPL). Applying this offset and searching for post-2015 events, we uncovered the promising occultation by Triton of a relatively bright star with V = 12.7, G =12.2 (UCAC4 410143659; Gaia DR2 2610107911326516992) that was to occur on 5 October 2017, crossing all Europe and northern Africa and reaching the eastern USA.

Closer in time to the event, a dedicated 8-night run on the OPD 1.6 m telescope was conducted between 15 and 23 September 2017 to further improve the accuracy of the prediction, in particular to pin down the path of the central flash. The field of view of the images was 6.1′ × 6.1′ with a pixel scale of 180 mas/pixel. At that point, digital coronagraphy (Assafin et al. 2009; Camargo et al. 2015), which mitigates Neptune’s scattered light, proved to be unnecessary. Chromatic refraction corrections to Triton’s position were carried out, but proved to be negligible too, as the observations were made in the I band. Observations from two nights were discarded due to bad weather, but we were able to cover a complete orbit of Triton around Neptune (about 6 days) with more than 1000 images.

Prior to Gaia Data Release 2 (DR2), published in April 2018 (Gaia Collaboration 2016, 2018), the Gaia team² released a preliminary Gaia subset DR2 with 431 stars (R = 12−17) surrounding Triton’s path in the sky plane during the eight nights of our run. The main improvement with respect to DR1 is the inclusion of the stars’ proper motions, leading to mas-level accuracy of stellar positions at epoch. Using these reference stars and the PRAIA package (Assafin et al. 2011), we obtained a mean offset (Δα cos δ, Δδ) = (+7.8 ± 5.4, −17.6 ± 2.6) mas with respect to the JPL DE435/NEP081 ephemeris. The corresponding prediction uncertainty was about 60 km cross-track and 8 s in time. Figure 2 displays the mean offset (Δα cos δ, Δδ) for the eight nights, using the Gaia DR2 catalogue. This is a large improvement from a prediction using only the Gaia DR1 catalogue, shifting the shadow path about 370 km to the south in the sky plane, or some 500–700 km when projected on Earth, depending on the station considered, as well as reducing the offsets uncertainties by factors of about 1.5 and 2.2, respectively. Since the run covers a full synchronous revolution, and rotational, period of Triton, the average offset reflects an error in Neptune’s heliocentric position rather than a neptunocentric error in Triton’s ephemeris.

The ± 60 km cross-track uncertainty on the prediction was essential for better planning observations of the central flash, keeping in mind that the width of the region where that central flash is significant (e.g. more than 20% of the unocculted stellar flux; see Fig. 3) typically spans ± 100 km in the sky plane.

Combining the pre-event Triton’s ephemeris offset with the geocentric astrometric Gaia DR2 position of the occulted star³, we obtained the parameters describing the nominal prediction listed in Table 1. This updated prediction was promptly released to the scientific community before the event⁴ and the corresponding map of the shadow track on Earth is displayed in Fig. 4.

Fig. 1 Sub-solar latitude on Triton over time. Upper panel: overview of the sub-solar latitude on Triton versus time over the last millennium. The blue part corresponds to the period from the V2 encounter (August 1989) to the 5 October 2017 stellar occultation. It shows that during this interval, Triton experienced an extreme summer solstice in its southern hemisphere, with a minimum sub-solar latitude of 50° S in 2000. Lower panel: close-up view of the upper panel around the year 2000. The black points correspond to occultations observed, distinguishing from the black triangle (V2 RSS experiment). The larger symbols are from the data that we use in this paper.

Fig. 2 Average offsets in (α, δ) mas of the eight nights of OPD observations using the preliminary Gaia DR2 catalogue for the astro-metric reduction. Positions in lighter colours were not accounted for to compute the average offset, due to lower quality associated with poor weather conditions. Lower panel: number of images taken for each night.

2.2 Observations of the occultation

The event was attempted from over one hundred stations in Europe, northern Africa, and the eastern USA, resulting in a total of 90 occultation light curves. Figure 5 displays the corresponding occultation chords and Table A.1 lists the circumstances of observations for the respective sites.

The first part of this table lists the stations that were eventually used in a simultaneous fit to a Triton’s atmospheric template model. The second part of the table lists other stations that were not used because the light curves had insufficient signal-to-noise ratios (S/N) to provide relevant contribution to that fit. The last part of the table provides information on the sites that were involved in the campaign, but could not gather data due to bad weather or technical problems.

The rule we used for deciding to include or not a light curve in the global fit is as follows. The S/N was first estimated by calculating the standard deviation σ of the flux outside the occultation. The resulting S/N 1/σ per data point was then re-calculated for a fixed time interval of one second. Assuming a gaussian noise, this implies a multiplicative factor of to be applied to the S/N obtained above. Finally, the noise level must be compared to the ‘useful’ information, that is, the actual drop of signal caused by the occultation, not the total signal. This implies a new multiplicative factor of Δϕ applied to the S/N per second obtained above, where Δϕ is the flux drop observed during the occultation. Consequently, some light curves were eliminated due to a low contrast caused by light contamination from Neptune (see for instance the Abington, Caserta, Agerola 50 cm, or Catania 28 cm light curves in Fig. C.2).

The normalisation described above then allows us to compared consistently the various datasets. We used a normalised S/N cut-off of ten prior to the global fit. There are a few exceptions to that rule for light curves containing a strong central flash (see for instance the Le Beausset and Felsina data in Fig. C.1). They have a poor overall S/N, but we have kept them in the global fit (without flash; Sect. 4.2) to test the effect of the central flash inclusion (see Fig. 3 and Sect. 7.2). This approach allowed us to ensure that the central flashes provide a shadow centre witch is consistent with, but more accurate than the one given by the global fit (see Sect. 7.2).

This said, the cut-off of ten remains somehow arbitrary. We have tested other cut-off values, and we did not obtain any significant changes in the results of the fits presented in the rest of the paper.

The analysis of the light curves (as described in the next section) allowed us to reconstruct the geometry of the event by providing Triton’s position with respect to the occulted star (see Table 1). The reconstructed geometry of the occultation implies a shift of (Δα cos δ, Δδ) = (−7.2, +0.6) mas of Triton’s position with respect to the latest prediction described in the previous section. This means that the occultation occurred about 9 s earlier than expected, and that the shadow centre was 12 km north of the predicted path in the sky plane.

This mismatch between observation and prediction is at a ~1.3σ-level, and is thus insignificant at our accuracy level. It shows in particular that the Gaia DR2 astrometry was crucial in getting an accurate prediction that allowed the detection of the central flash at various stations.

Fig. 3 Simultaneous fits of the data (black dots) by synthetic light curves (blue lines), based on the temperature profile displayed in Fig. 10 (black line) and the pressure boundary condition p1400 = 1.18 µbar (Table 4). The green curves are the residuals of the fits. The lower and upper horizontal dotted lines mark the zero-flux level and the total star plus Triton unocculted flux, respectively. We also note that the three bottom light curves are plotted at a different vertical scale from the others to accommodate the presence of a strong central flash. The stations are sorted from left to right and top to bottom from the northernmost track (St Caprais) to the southernmost track (Calar Alto; see next figure). Each panel has a duration of five minutes and is centred around the time of closest approach (or mid-occultation time) of the station to Triton’s shadow centre, as indicated under the lower left panel in each block of six light curves. For reference, the vertical red line marks the time 23:48 UTC. The stations with exposure times of less than 1 s have been smoothed to have a sampling time as close as possible to 1 s, for easier S/N comparison of the various datasets. We note that in this approach, the sampling of the Constância, Le Beausset, and Felsina Observatory light curves (0.64 s) is kept at its original value so that full resolution versions of the corresponding strong flashes at those stations are displayed here. The same kinds of plots showing all the stations, but with the flashes excluded from the fits, are displayed in Fig. C.1.

Fig. 3 continued. NB. ‘JAVA.’ is the abbreviation of Javalambre, used so that the name of the station fits into the plot.

2.3 Occultation light curves

All our occultation light curves are displayed in Figs. C.1-C.2. As was done in Table A.1, the first group of figures corresponds to light curves that had sufficient S/N to be used in Triton’s atmospheric fit, while the second group is for light curves with lower S/N that were not used in the fit. It should be noted, however, that the best synthetic models expected for those light curves (plotted in grey in Fig. C.2) are fully consistent with the observations, in the limit of the noise level.

Among those light curves, three were used for obtaining atmospheric profiles from an Abel inversion procedure (see Sect. 4). Two of them (La Palma and Helmos) are displayed in Fig. 6, and the third one (Calern) is shown in the upper panel of Fig. 7.

Fig. 4 Triton’s shadow path on 5 October 2017. The black dots are spaced by one minute, the arrow indicating the direction of motion of the shadow. The northern and southern limits of the solid body assuming a radius of 1353 km are also represented, with the predicted path as the white lines and the effective path as the blue lines. The grey area represents night on Earth (dark grey for full astronomical night and light grey for twilight). Stations with a successful observation that have been used in our fit are represented by blue dots, while the ones that were not used are shown as red dots. The white dots are the stations that attempted observation but were clouded out or had technical difficulties. Upper panel: overview of all observing stations. The larger black dot along the black line corresponds to the closest approach of the shadow centre to the geocentre, at around 23h52 UTC (see Table 1). Lower panel: closer look at the central flash path across Europe. The grey lines around the centremost line correspond to a spacing of 50 km (corresponding to about 2.5 mas) once projected in the sky plane.

3 Retrieving Triton’s atmospheric structure

3.1 Methodology

We first adopt a bootstrap method to retrieve the molecular density n(r), pressure p(r), and temperature T(r) of Triton’s atmosphere as a function of the distance to Triton’s centre, r. To do so, one approach is the Abel inversion of our refractive occultation light curves that have the highest S/N. Its primary result is the density profiles n(r), from which the p(r) and T(r) profiles are derived by using the hydrostatic and ideal gas equations.

The other approach is a direct one. It is used once the inversion procedures have provided the density profiles n(r). These profiles are smoothed and parameterised according to physical arguments (discussed later). Then, a ray-tracing scheme generates synthetic light curves that are fitted to the occultation light curves, thus describing the global structure of Triton’s atmosphere, as illustrated in Fig. 5. One product of this fit is the location of Triton’s centre, which is used iteratively with the Abel inversion, thus improving the accuracy on the altitude scale. The other product of our approach is the value of the pressure at a prescribed radius for comparison with previous results, aimed at detecting possible long-term seasonal effects.

The final approach is to fit the central flashes observed in some light curves and to measure the departure (if any) from the spherical shape of Triton’s deep atmosphere, which is eventually used to constrain its wind regime. The other goal of this fit is to reveal possible absorbing material (by comparing the height of the central flash in stations that provided observations in different wavelengths) along the line of sight, such as hazes just above Triton’s surface.

Technical details on the inversion technique are given in Vapillon et al. (1973), ray-tracing schemes and central flash fitting are described in Sicardy et al. (1999), and applications to Pluto’s atmosphere are presented in Dias-Oliveira et al. (2015, DO15 hereafter), and Meza et al. (2019). Numerical values used in both our inversion and ray-tracing codes are summarised in Table 2.

Fig. 5 Geometry of the 5 October 2017 stellar occultation by Triton, as seen in the sky plane. The J2000 celestial north (N) and east (E) directions and the scale are indicated in the upper right corner. Triton’s radius is fixed at RT = 1353 km, and the grey arrow near the equator shows the direction of rotation of the satellite. The (Neptune-facing) prime meridian is drawn as a thicker line compared to the other meridians, and the south pole is marked by the label S. The inclined lines are the trajectories of the star relative to Triton (or ‘occultation chords’) as observed from various stations, with the black arrow indicating the direction of motion. We gathered a total of 90 occultation light curves, 52 of which (corresponding to the blue colour, as in Fig. 4) had sufficient S/N to be included in a global atmospheric fit; the remaining 38 (red colour) with lower S/N were not included in the fit.

Fig. 6 Best two light curves obtained during the Triton occultation of 5 October 2017, at La Palma and Helmos stations (see details in Table 1). Both telescopes were equipped with the same E2V CCD 47-20 detector with quantum efficiency peaking at 600 nm and reaching zero near 300 and 1000 nm, respectively. Left panel: full resolution light curve (cycle time 0.635 s) for the La Palma station. Right panel: same for the Helmos station (cycle time 0.674 s). The spectral ranges used for each instrument are indicated in the figures (I+z at La Palma and V+R at Helmos). The blue lines are the best simultaneous fits obtained with our ray-tracing approach (see Sect. 4).

3.2 Assumptions

Our inversion and ray-tracing schemes assume that:

• (1)

  The atmosphere is composed of pure N₂. The next most abundant species (CH₄) has a volume mixing ratio (hereafter referred to as mixing ratio) [CH₄/N₂] of less than 10⁻³ (Strobel & Summers 1995; Lellouch et al. 2010). Our ray-tracing code shows that such an abundance causes a fractional change of the synthetic flux of about 10⁻⁵ near the half-light level, these effects are negligible considering the noise level of the data.

• (2)

  The atmosphere is transparent. The deepest layers reached during Earth-based occultations are those that cause the central flash, at an altitude of about 8 km (see Appendix B). The validity of this assumption will be discussed in Sect. 8.

• (3)

  The upper atmosphere is globally spherical. This hypothesis is supported by the fact that the observed central flashes are consistent with a spherical shape (Sect. 7). Small departures from the spherical model, however, are observed in some central flash shapes, and are discussed in Sect. 7.

The limitations of our approach described above are presented in Appendix B.

4 Results

4.1 Inverted profiles

We used the three datasets with highest S/N to perform ourinversion method, more precisely the light curves from La Palma (2-m Liverpool telescope, Spain), Helmos (2.28-m Aristarchos telescope, Greece) and Calern (1.04-m C2PU telescope, France) (see Table A.1). At half-light times (where the star flux has been reduced by 50%) and for ingress and egress, each of these stations probe different locations in Triton’s atmosphere. The corresponding latitudes, longitudes, and local solar times of the sub-occultation points are provided in Table 3. The paths of the stellar images over Triton’s surface as seen from these stations are plotted in Fig. 8.

The results of the inversions are displayed in Figs. 9–11. We also plot in these figures the profiles retrieved from our analysis of the RSS occultation (see Sect. 5.1 for more details, as well as a discussion on the connection of these profiles with our results). There are five noteworthy features. First, all six n(r) and p(r) inverted profiles are very similar, showing that the stations at La Palma, Helmos and Calern probed essentially identical atmospheric layers at their ingress and egress points. No significant variations versus local time and latitude are observed.

Second, in their common range of probed altitudes, our density profiles and the RSS profile coincide, to within the noise level of the RSS experiment (the noise level in our retrieved profiles being much smaller). Third, the pressure profiles from RSS and from our inversions are also close to each other. However, contrary to the density profile, some small differences appear. This is discussed in Sect. 5.

Fourth, the general positive gradient in the upper parts of our retrieved thermal profiles is a mere result of the choice of our initial conditions. This is intended to match the general temperature profiles obtained independently by Strobel & Zhu (2017), which were constrained by the RSS occultation data taken in 1989 (see Appendix B).

Fifth, all of our six retrieved thermal profiles show, however, a marked turning point in their deepest parts, where the temperature gradient becomes negative. This gradient is always well below (in absolute value) the local dry adiabatic gradient, so that the atmosphere is convectively stable in those parts (Fig. 11).

4.2 Ray-tracing approach

For all the datasets used here, we employed the same procedure as in DO15, which consists of simultaneously fitting the refractive occultation light curves with synthetic profiles generated by the ray-tracing code. For each station, a least-squares fit is performed to adjust the synthetic light curve to the observation. Due to the uncertainties in the determination of ϕ₀ (the fraction of the flux attributed to Triton; see Appendix B) for some stations, and the lack of calibration for most of them, we considered ϕ₀ as a free parameter when performing these fits. We note that this adds one degree of freedom per station to the fit, and thus increases the error bars on the retrieved atmospheric parameters.

Our ray-tracing method is mainly sensitive to the half-light level. It corresponds to a radius of about 1415 km in Triton’s atmosphere (altitude ~60 km) and a pressure level of ~0.55 µbar. For a prescribed temperature profile T(r), this method returns two best fitting parameters. One parameter is the pressure p_ref at a reference radius, here r_ref = 1400 km. This particular choice stems from the fact that this reference radius has been used in previous works (e.g. Olkin et al. 1997; Elliot et al. 2000b), thus allowing consistent comparisons.

To proceed forward, we have defined a template T(r) profile that matches the inverted profiles obtained at the station with the best S/N (La Palma). It has the same functional variation with altitude as in Dias-Oliveira et al. (2015), where it was applied to Pluto’s atmosphere, except for the upper branch, which is not isothermal but rather has a constant thermal gradient that connects the lower atmosphere to an upper thermosphere. The adopted parameters for the template profile T(r) are provided in Appendix B.2.

The lower part of the profile has been adjusted so as to fit the central flashes (see Appendix B.2 for details). That adjustment provides constraints on the thermal profile between the lowest inverted point of La Palma (~20 km altitude) down to the central flash level (~8 km altitude). Finally, below 8 km, the profile has been connected to the surface’s temperature at 38 K. This is a ‘blind part’ of the profile, as it does not contribute significantly to the refracted stellar flux received on Earth.

The resulting synthetic temperature profile is shown as a thin black line in Fig. 10. Starting from the surface, the profile first has a strong positive temperature gradient of 5 K km⁻¹. This gradient decreases rapidly (Fig. 11) and the temperature reaches a maximum value of about 50 K at r = 1363 km (10 km altitude), thus implying an average gradient of 1.2 K km⁻¹ in that lower part. Our data show a hint of a mesosphere with a negative gradient (also seen in Elliot et al. 2003) that reaches −0.2 K km⁻¹ at r = 1375 km (23 km altitude), before connecting with the general positive gradient of the upper branch.

The strong surface temperature gradient at the surface derives from the need to connect our inverted profiles to the surface at 38 K. Since we do not have information in this lower portion of the atmosphere, we employed the simple hyperbolic form of DO15 to connect our profile to the surface, so that our surface gradient does not necessarily reflect the real value at that level.

This said, the general positive gradient can be achieved by considering the heating by CH4 stemming from near IR absorbing bands. For instance, we estimate that a CH₄ mixing ratio of 0.0004 yields T = 52 K at 1363 km, and thus could explain our result. Strobel & Zhu (2017) ran their model for discrete values of the CH₄ mixing ratios not included in their paper and found that a CH₄ surface mixing ratio ~0.00015 would suffice to support a temperature rise ~9 K, and a CH₄ surface mixing ratio ~0.0004 a temperature rise ~ 12 K, in the first 10 km. Because CH₄ is photochemically destroyed in the lower atmosphere, its scale height is roughly half the N₂ scale height and in terms of CH4 column density one needs a higher surface CH₄ mixing ratio to compensate for its smaller scale height. For remote sensing observations it is the column density that is important and not just the surface mixing ratio that is relevant. We note that the 0.0004 value is smaller than, but roughly consistent with the range found by Lellouch et al. (2010) for the CH4 mixing ratio, 0.0005–0.0010. Moreover, some complications may arise, like the existence of a troposphere.

The troposphere on Triton has been shown to be controlled by turbulent mixing above the surface, and to be sensitive to surface thermal contrasts between N2 ice and the volatile free bedrock (due to different surface albedo or thermal inertia, Vangvichith 2013). On Pluto, climate models showed that the sublimation of cold N₂ ice and subsequent transport of the cold N2 air in the impact basin Sputnik Planitia yield a km-thick cold troposphere as observed by New Horizons Forget et al. (2017); Hinson et al. (2017).

The negative gradient in the mesosphere, reminiscent of the more extended mesosphere on Pluto (Lellouch et al. 2017; Young et al. 2018), calls for the existence of a coolant. It must cool the atmosphere above its peak temperature of ~50 K, as well as radiate away the downward thermal heat flux from the upper atmosphere where T ~ 100 K. There are a few candidates for this coolant: haze particles and/or influx of dust particles that may either be pure H₂O ice or with silicate cores and coated with H₂O ice (see Ohno et al. 2021 for more details).

The pressure at any level can then be deduced by using the temperature template described above. In particular, the surface pressure p_surf can be obtained by the relation p_surf = 12.0 × p₁₄₀₀. This ratio will be used to extrapolate p₁₄₀₀ from p_surf or vice versa. The other fitted parameter is Triton’s DE435/NEP081 ephemeris offset perpendicular to its apparent motion projected in the sky, Δρ. We note that the ephemeris offset along Triton’s motion is decoupled from that fit (see DO15 for details).

Error bars are obtained from the classical function , which reflects the noise level σ_i of each of the N data points, where ϕ_i,0bs and ϕ_i,syn are the observed and synthetic fluxes at the i^th data point, respectively.

We simultaneously fitted a selected pool of 52 light curves obtained during the 5 October 2017 occultation. Other light curves were not considered at this stage because they are affected by higher or non-normal noise that would degrade the global fit. In a first step, we excluded from the fit the parts of the light curves where a strong central flash is present. This is to avoid giving too much weight to those parts, while they reflect only the properties of the deepest atmospheric layers. So, the goal of this fit is to get the global properties of the atmosphere, and in particular, to constrain p₁₄₀₀ and the location of the shadow centre with respect to the occultation chords. In a second step, we included the central flashes in the fit to assess the shape of Triton’s atmosphere and to check if the central flash location coincides with the centre found by the global approach.

After exploring a grid of values for Δρ and p₁₄₀₀, we obtained the χ²(Δρ, p₁₄₀₀) map displayed in Fig. 12. A satisfactory fit should provide a minimum value close to N − M, where M is the number of fitted parameters. Besides p₁₄₀₀ and Δρ, we considered Triton’s contribution to the light curve (ϕ₀) as a free parameter for all light curves. This is because no satisfactory values of ϕ₀ have been obtained for any of the light curves (see Appendix B.2). Thus, the fitted parameters are the values of Δρ and p₁₄₀₀, plus the 52 values of the ϕ₀ (i.e. a total of M = 54 fitted parameters). On the other hand, we used N = 68 446 data points. We then obtain a global value of χ² per degree of freedom, , indicating a satisfactory global fit to the data. An examination of values of for individual light curves also show values near unity for all of them. Thus, none of our light curves show significant discrepancies when compared to the synthetic light curves derived from the synthetic density model shown in Fig. 9. This confirms the spherical symmetry of Triton’s atmosphere on a global scale.

Without considering Δρ, the marginal distribution⁵ for 1σ and 3σ error contours on p₁₄₀₀ are estimated by tracing the iso-levels and , respectively, as shown in Fig. 12. The best-fitting value of p₁₄₀₀, its 1σ error bar and the quality of the fit, , are listed in Table 4. The best-fitting value of Δρ (−359.3 ± 1 km, Fig. 12) is used to retrieve the closest geocentric approach distance between Triton and the star (projected in the sky plane) and its corresponding time (see Table 1).

Finally, the best simultaneous fit corresponding to the minimum of χ² is displayed in Fig. C.1–C. For the sake of completeness, Figs. C.2–C show the synthetic light curves superimposed on the light curves that have not been included in the fit. Although they have poorer S/N, they all confirm that our global model satisfactorily fits these data.

Fig. 7 More detailed views taken of the flashes observed at Calern (top panel) and Constância (bottom panel), with the same setup as in the previous figure.

Fig. 8 Tracks of the primary (red dots) and secondary (blue dots) stellar images above Triton’s surface, as seen from Constância, plotted every 0.1 s. The junctions between the red and blue paths correspond to ingress (left) and egress (right) points for the Constância station. The arrows show the direction of the stellar images’ path. The regions probed by the central flash are those where the dots are more spaced. All the other stations probed essentially the same path (or part of it), with the primary and secondary images being swapped (as well as their directions of motion), depending on whether the station probed north or south of the shadow centre. Since the Earth and the Sun are angularly close (-1°) to each other as seen from Triton, the stellar paths essentially mark Triton’s terminator, the night side extending above the terminator in this figure. The two yellow symbols are for La Palma station, with ingress plotted as a star and egress plotted as a diamond. The two green symbols are the same for Helmos station and the two white symbols for Calern station (see Table 3 for the corresponding values of the latitudes and longitudes). The background image is a global colour map of Triton, produced using V2 data and orange, green, and blue filter images in order to obtain an approximation of Triton’s natural colours. Background image credits: image selection, radiometric calibration, geographic registration and photometric correction, and final mosaic assembly were performed by Dr. Paul Schenk at the Lunar and Planetary Institute, Houston, Texas. Image data are from V2 (NASA, JPL).

5 Reanalysis of results from previous events

5.1 The Voyager 2 radio occultation

During its Triton flyby on the 25 August 1989, the V2 spacecraft sent its radio signal (RSS experiment) back to Earth as it passed behind the satellite. Details on the gathering of the V2 RSS data are given in G95.

The main product of this observation was the temperature and pressure at Triton’s surface. However, although it becomes quite noisy above the 20 km altitude level, the RSS phase delay still provides useful constraints on Triton’s lower atmosphere, with some science left to explore. Here, we give a summary of Gurrola’s work, and describe how we use the V2 phase delay to retrieve Triton’s atmospheric structure in the 10-20 km above the surface.

The V2 high gain 3.7-m antenna transmitted to Earth two radio signals at 3.6 and 13 cm (X band and S band, respectively). The phases in these bands are related to one another by (1)

where ϕ_x is the phase in the X band, ϕ_s is the phase in the S band, and Δϕ is the corrected radio phase corresponding to the neutral atmosphere at the X-band wavelength. This calculation is done to remove plasma effects on the phase. Due to problems in fitting the ingress data, as there seemed to be sudden changes in slope, G95 used only the egress data for his analysis.

Gurrola provided us with the corrected phase delay Δϕ(r) versus altitude above Triton’s surface, as well as the results from his models to obtain only the ‘pure atmosphere phase delay’. This corresponds to the phase delay once a general polynomial trend has been subtracted from Δϕ to account for thermal noise and instabilities in the frequency reference on board V2. These polynomials, referred to as baselines in G95, were designated as B₁, B₂, and B₃. B₁ is the linear baseline used by Tyler et al. (1989), determined using 120 km of the data obtained. G95 considered this insufficient to reliably estimate the drift of the instrument over the atmosphere, as it did not extrapolate from high enough altitudes (so that the atmosphere is too thin to affect the signal phase) downwards towards Triton’s surface. On the other hand, baselines B₂ and B₃ used about 700 km of the data, and are, respectively, the second and third-order polynomials of G95’s best fit at egress. The preferred solution of this author is B₂. The resulting Δϕ(r) is displayed in green in Fig. 13.

Using the B₂ solution, we derive the profiles displayed as green curves in Fig. 9. In order to compare this result to ours, we generated for comparison the phase delay at 3.6 cm that would be observed with our best profiles n(r) (the black line in Fig. 9) as if it were obtained by V2: (2)

where λ is the wavelength (3.6 cm), K is the corresponding molecular refractivity of N₂ (see Table 2), and σN₂(r) is now the column density stemming from our best model. The resulting Δϕ(r) profile is shown in black in Fig. 13, together with the phase delays deduced from the inversions of La Palma and Heimos’ light curves (in colours).

Conversely, we used the V2 corrected X-band radio phase to retrieve the refractivity profile from the Abel inversion (3)

using the auxiliary variable to calculate the integral. Finally, the density profile n(r) = v(r)/K can be deduced.

With this, we can directly compare our results to those of V2. We note that the RSS profiles probe altitude interval levels that overlap our ground-based occultation levels. This overlapping region extends from the lowest inverted points of the La Palma station, r = 1373 km (20 km altitude), up to roughly the reference level, r = 1400 km (47 km altitude), at which point the RSS profiles become too noisy to be reliable, reaching a factor of about 2 at that level.

The examination of the upper panel of Fig. 14 shows that no significant difference in density is detected between the 1989 and 2017 profiles, especially at the ‘junction level’ at r = 1373 km. We note that the RSS density profile is rather insensitive to the particular solution B₁, B₂, or B₃ chosen to retrieve n(r).

Integrating the weight of the atmospheric column provides the RSS pressure profile (lower panel of Fig. 9). However, this profile includes a contribution of the weight of the layers above r = 1400 km, where the RSS phase delay is very noisy. We note that the B₁ pressure profile is quite offset in slope with respect to the B₂ and B₃ solutions (lower panel of Fig. 14). So, and contrar-ily to the density profiles, the general slope of log₁₀(p) versus radius r is quite sensitive to the particular choice of the polynomial baseline. In this context, it is difficult to conclude if the break in slope between the RSS profile and our ground-based results of 2017 is real or not.

Using the preferred B₂ model, and using the Abel inversion (Eq. (3)), we obtain p_surf,RSS = 13.6 µbar and p_1373,RSS = 3.77 µbar as of 1989. The main result of the inversion is the density profile; to translate this into pressure, we need an estimate of the surface temperature T_surf. The error bar on p_surF,RSS caused by the uncertainties on T_surf will be discussed next.

In any instance, our results are consistent with the analysis of G95, p_surf,RSS = 14 ± 2 µbar. This is also fully consistent with our estimation of the surface pressure as of 2017, 14.1 ± 0.3 µbar (Table 4). Thus, no significant variations of surface pressure is found when comparing the RSS results of 1989 and the results derived from the ground-based occultation of 2017.

The value p_1373,RSS = 3.77 µbar that we find is 21% smaller than the value we obtained in 2017 at that level, 4.58 µbar. Propagating this 21% difference to the 1400 km radius then yields p1400, RSS = 0.97 µbar. Estimating the error bar on that value is difficult, because the RSS pressure depends on the (noisy) pressure values obtained above, as mentioned earlier. If we adopt the error bar p_surf,RSS = 14 ± 2 µbar of G95 and propagate it upwards, this yields p_1400,RSS = 0.97 ± 0.14 µbar.

Another, more robust way to estimate p_surf,RSS is to use the RSS density profile alone. The counterpart is that we need an independent measurement of the surface temperature T_surf in order to derive the pressure from the ideal gas equation . These temperature measurements (given below) are more accurate than G95’s estimation (T_surf = 42 ± 8 K) and thus reduce the ± 2 µbar uncertainty of G95’s value of p_surf,RSS. However, this approach is valid only if these temperature measurements apply to the N₂ ice surface that the RSS experiment probed, and if the vapour pressure equilibrium between the N₂ ice and the gas is achieved.

Estimations of T_surf are given by various authors (see also Fig. 15): (Conrath et al. 1989), (Tryka et al. 1993), the range 36.5-41 K (Grundy et al. 1993) and 37.5 ± 1 K (Merlin et al. 2018). Adopting a value of n_surf,RSS = 2.4 × 10¹⁵ cm⁻³ derived from our RSS phase delay inversion (Fig. 14), we find surface pressures of , a range 11.5–13.3 µbar and 12.4 ± 0.3 µbar, respectively, for the four choices of surface temperatures. We note that all these values are consistent with the surface being in vapour pressure equilibrium with the atmosphere, as shown in Fig. 15. This supports the hypothesis that the reported temperatures are indeed representative of the N₂ ice surface.

In summary, we estimate from the surface temperatures given above, a safe surface pressure range of 12.5 ± 0.5 µbar can be derived at the V2 epoch. In this case, the error bar essentially stems from the uncertainties on the temperatures. Comparing this value with our estimation p_surf = 14.1 ± 0.4 µbar in 2017 (Table 4), and assuming a constant factor (12.5/14.1) throughout the profile, we formally obtain p₁₄₀₀ = 1.05 ± 0.04 µbar in 1989. This is consistent with our estimate made above, p_1400,RSS = 0.97 ± 0.14 µbar. We thus estimate a conservative range of p_1400,RSS = 1.0 ± 0.2 µbar for the pressure at 1400 km in 1989.

Fig. 9 Density and pressure proflies of the atmosphere of Triton as a function of radius, r (the distance to Triton’s centre). Upper panel: density profiles of Triton’s atmosphere as a function of radius, retrieved by inverting three light curves obtained during the 5 October 2017 occultation and from the V2 radio phase delay at 3.6 cm. The colour codes are indicated in the upper right part of the plot. The same codes are used in Figs. 10–11 and Figs. 14–15. The thin black curve is a smooth synthetic density profile that fits the inverted profiles and is extrapolated down to the surface. It is derived from the smooth temperature profiles shown in Fig. 10. The solid horizontal line marks Triton’s surface (at radius RT = 1353 km), the dashed line indicates the central flash layer (near 1360 km), and the dotted horizontal line marks the reference radius, rref = 1400 km. Lower panel: corresponding pressure profiles.

Fig. 10 Temperature profiles as a function of radius (upper panel) and pressure (lower panel). The oblique dotted line in the lower panel is the wet adiabat, i.e. the vapour pressure equilibrium line for N2, taken from Fray & Schmitt (2009).

Fig. 11 Temperature gradient corresponding to the upper panel of Fig. 10. The dot-dashed line is the dry adiabatic temperature gradient Γ= − g/cp, i.e. the limit of convective instability, where cp is the specific heat at constant pressure for N2 and g = GMT/r2 is the acceleration due to gravity.

Fig. 12 χ2 map for the simultaneous fit of 52 light curves obtained during the occultation of 5 October 2017, using the ϕ0 corresponding to the temperature profile in Fig. 10. The inner green line contour is the 1σ limit of the fit, while the outer green line is the 3σ limit.

5.2 The 18 July 1997 stellar occultation

This campaign involved one station in the USA and three stations in Australia. It was a joint effort between two groups, and, therefore, both have access to the data. The circumstances of observations are listed in Table A.1 and the geometry of the event is displayed in Fig. 16. More details on these observations and their analysis are given in Elliot et al. (2000a).

Here we provide the results of our own approach to constrain p₁₄₀₀. In particular, we adopt the same temperature profile T(r) as for 2017 (see Fig. 10), but varying p₁₄₀₀ to fit the synthetic light curves to the data. The (χ², Δρ) map is displayed in Fig. 17 and the best fit is shown in Fig. 18. This yields µbar and , indicating a satisfactory fit.

Fig. 13 Radio phase delay during the egress of the V2 RSS occultation on 25 August 1989 in the 3.6 cm X band. Left panel: radio phase delay observed (adapted from G95). The crosses are the data, and the three solid lines (labelled B1, B2, and B3) are three polynomial modellings of the phase delay baseline. As discussed by G95, the preferred baseline solution is B2. Right panel: radio phase after subtraction of the B2 baseline polynomial shown in the left panel (green crosses), thus representing the effect of the atmosphere only. The green line marks the smooth version of that radio phase delay, as constructed by G95. The red profile is the phase delay that would be observed at 3.6 cm from the retrieved density profile of La Palma at immersion (see Fig. 9). Other phase delay profiles obtained from La Palma (emersion) and Helmos would be indistinguishable from the red profile and are not plotted here for sake of clarity. The black profile is the phase delay obtained from our best model of Triton’s atmosphere (see text for details).

5.3 The 21 May 2008 stellar occultation

This event was observed from Namibia (two stations) and from La Réunion island (two stations; see Table A.1). Given that each pair of stations are close together, only two effective chords have been obtained (see Fig. 19). Moreover, these chords being grazing, there is a strong correlation between the closest approach distances of the chords to Triton’s shadow centre and the retrieved reference pressure p₁₄₀₀ (see Fig. 20). The best fit is also shown in Fig. 21. As a consequence, the value of p₁₄₀₀ is poorly constrained at the 1σ level, . At the 3σ level, the value is so unconstrained that it does not bring any information on the temporal seasonal variations of the pressure (Fig. 22).

6 Atmospheric seasonal variations

6.1 Occultation results

Table 4 lists our values of p₁₄₀₀ at various epochs, as well as values taken from other works. Extrapolations to the surface have also been included, assuming a constant ratio p_surf/p₁₄₀₀ = 12.0. The corresponding seasonal variations of p₁₄₀₀ with time is displayed in Fig. 22.

The value of Olkin et al. (1997) indicates a 40% increase in pressure between 1989 and 1995, but at a low significance level of 1.8e. From the 18 July 1997 event, Elliot et al. (2000a) obtained p₁₄₀₀ = 2.23 ± 0.28 µbar, whereas with the same dataset, we obtain . The difference between the two results amounts to a factor of 0.85 and stems from the use of a different template model T(r). This said, this difference remains at the 0.6σ level and is statistically insignificant. Using our value of p₁₄₀₀ for 1997 indicates a pressure increase by a factor of 1.9 between 1989 and 1997, but at a marginally significant 2.5σ level only.

The 4 November 1997 value obtained by Elliot et al. (2000a, 2003), p₁₄₀₀ = 1.76 ± 0.02 µbar, has a much lower error bar due to the high S/N of the light curve, obtained with the Hubble Space Telescope. Taken at face value, this implies an increase in pressure by a factor of 1.76 between 1989 and 1997, at a 3.8σ level. However, this observation was a single-chord event, and a model was used to retrieve the astrometry of this event. Consequently, there is an uncertainty that was not accounted for. Since we do not have access to these data, it is impossible for us to verify their result using our own methods, and therefore, confirm this increase.

Finally, the 21 May 2008 event provided only two grazing chords, bringing no new information. Thus, no firm conclusion can be drawn on any change of pressure between 1989 and 2008.

In summary, we estimate that the surge of pressure reported in the 1990s (compared to the V2 epoch) seems to be confirmed by our own analysis, but it remains debatable considering the paucity of data points available, and the lack of a fully consistent analysis of all the observed events. We note that the 2017 data rules out the concept of a monotonic increase in Triton’s pressure over time, but does not rule out the observed increase in 1995–1997. Regardless, the much more accurate value of p₁₄₀₀ that we obtain in 2017 is fully compatible with that derived from the V2 RSS experiment. If we consider the 3σ level, Fig. 22 shows that no increase can be claimed between the two measurements. So, either no surge occurred between 1989 and 2017 or, if it did, the pressure was back to its V2 value in 2017.

From high-resolution spectroscopy in July 2009, Lellouch et al. (2010) obtained the first detection of methane gas in Triton’s atmosphere since V2, and the first CO gas detection. Their analysis yielded a CH₄ gas number density at the surface larger than inferred from V2 (Herbert & Sandel 1991; Strobel & Summers 1995). Assuming that the N₂ pressure would qualitatively follow a similar seasonal variation, they estimated a 40 µbar pressure in 2009. This value (which did not represent a direct measurement of the N2 pressure) is clearly at odds with the picture shown in Fig. 22, in particular with the 21 May 2008 point.

Fig. 14 Close-up view of Fig. 9. Upper panel: close-up view of the density profiles. Various profiles derived from the RSS occultation are shown in green. The dashed line is the profile retrieved by Tyler et al. (1989) using a polynomial extrapolation B1 to correct for the RSS phase instability. The thin solid line is the profile retrieved by G95 using the polynomial extrapolation B2. The dotted line is the same by G95, but using the polynomial extrapolation B3. The thick solid line is the best model of G95, based on the B2 profile. Lower panel: same for the pressure profiles.

Fig. 15 Close-up view of the lower panel of Fig. 10. The width of each coloured box is Triton’s surface temperature (Tsurf), as estimated by the various authors mentioned just above the boxes (see text for details). The heights of the boxes are the range of surface pressures, psurf, using the ideal gas law , where is the surface molecular nitrogen density derived from our inversion of the RSS data (see text). We note that all the boxes intersect the vapour pressure equilibrium, which is plotted as a dotted line. This shows that the RSS surface density and the estimated surface temperatures are mutually consistent with a pressure being controlled by the N2 ice sublimation.

Fig. 16 Geometry of the 18 July 1997 occultation, with the same conventions as in Fig. 5.

Fig. 17 Same as Fig. 12, but for the 18 July 1997 occultation.

6.2 Climatic context from numerical volatile transport modelling

The climatic context of Triton is described and analysed in detail in a recent paper by Bertrand et al. (2022). Bertrand et al. (2022) employed the VTM of Triton, developed at the Laboratoire de Météorologie Dynamique (LMD), to investigate the long-term and seasonal volatile cycles of N₂ and CH₄ on Triton. Their simulations are constrained by the surface pressure derived from the stellar occultations presented in this paper. In this section, we summarise the main results of this paper that are relevant for the interpretation of our observations.

In VTM simulations, the surface pressure peak occurs slightly after the southern summer solstice (2000) between years 2000 and 2010 (see Fig. 23). The surface pressure seasonal variations obtained by Bertrand et al. (2022) is similar to that obtained by Spencer & Moore (1992), when they artificially maintained a permanent large southern cap of bright N₂ (see their Fig. 7). The larger the northern cap, the more it can serve as a condensation area and buffer N₂ sublimation in the southern hemisphere, which results in a lower and earlier surface pressure peak. The amplitude of the surface pressure peak is strongly attenuated if N₂ ice remains between 30° S - 0°, because condensation will dominate over sublimation between the years 1980 and 2020.

According to the model, Triton’s atmospheric surface pressure will remain at 5 µbar during the next solstice season if the north polar cap extends to 60° N and the south polar cap extends to 0°. The amplitude of the pressure peak is attenuated if N₂ ice deposits remain between 30° S - 0° because these latitudes are dominated by condensation rather than sublimation after the year 2000.

These results suggest that a northern cap extending down to at least 45° N - 60° N is needed in 2017 to restore the surface pressure back to the V2 measured value ~14 µbar. Otherwise, the surface pressure will remain higher than 16 µbar in 2017 with no northern cap. A strong increase in surface pressure cannot occur before 2000 if N₂ ice remains between 30° S - 0°. To ensure that the surface pressure remains greater than 5 µbar during the opposite season (southern winter) a permanent northern cap extending down to 45° N is required.

In their simulations, Bertrand et al. (2022) also investigated the CH₄ cycle by taking into account a small amount of pure CH₄ ice at the surface in addition to the N₂-rich mixture (Merlin et al. 2018). In the case where this pure CH₄ ice is placed at the south pole of Triton, covering 2% of the surface of the visible projected disk, they obtain a large increase in the CH₄ gas abundance from 1990 to 2005, without any significant change in N₂ surface pressure. Since CH₄ is not completely mixed with N₂ ice, it implies that the large increase in CH₄ (with relation to V2) reported by Lellouch et al. (2010) could be decoupled from the N₂ seasonal variations and, therefore, does not necessarily represent a measurement of the global pressure of the atmosphere.

For more details on the N₂ and CH₄ seasonal variations as simulated by the VTM, the reader is referred to Bertrand et al. (2022).

Fig. 18 Simultaneous fits to the 18 July 1997 light curves. The panel covers 300 s in time. All the light curves have been shifted so that the mid-occultation times are aligned. The red vertical tick marks indicate 10:10 UTC at the Brownsville station (USA) and 10:18 UTC for the three Australian stations. The blue lines are simultaneous fits to the data (black dots), using the best value found in Fig. 17, p1400 = 1.9 µbar, and the temperature profile shown in Fig. 10. The green dots are the fit residuals. For each light curve, the upper dotted line is the normalised value of the star plus Triton flux, and the lower dotted line is the background flux. We note that the data from the Brownsville station were normalised during the event, as shown in Elliot et al. (2000a).

Fig. 19 Geometry of the 21 May 2008 occultation, with the same conventions as in Fig. 5.

Fig. 20 Same as Fig. 12 but for the 21 May 2008 occultation.

Fig. 21 Fit to the 21 May 2008 light curves. The same conventions as for Fig. 18 are used, except that the panel now covers 720 s in time. The synthetic light curve for Piton Lacroix is plotted in grey because it is not used in the fit, due to the high noise level. The red vertical tick marks indicate 01:51 UTC for Piton Lacroix and Maïdo, and 01:41 UTC for Hakos.

Fig. 22 Triton’s atmospheric pressure seasonal variations with time, using the values from Table 4. Our results are in red, and values taken from other works are plotted in blue. For better viewing, the value derived by Elliot et al. (2000a) from the 18 July 1997 occultation is plotted in a semi-transparent blue colour, using the same dataset that we are using here for that date, with our red diamond-shaped point just below it. For all points, the thick error bars correspond to 1σ confidence levels. We note that for the 4 November 1997 and 5 October 2017 values, the 1σ error bar is smaller than the diamond-shaped symbol and therefore is not visible. For the 18 July 1997 and 21 May 2008 events, we have also plotted for information our 3σ error bars as thinner lines (see text for discussion).

7 Triton’s lower atmosphere: central flash

The detection of a central flash during the 5 October 2017 occultation offered a unique opportunity to study Triton’s lower atmosphere. Our ray-tracing code shows that the flash is caused by a layer having a typical thickness 2 km, lying at about 8 km above Triton’s surface (radius of 1361 km). In that altitude range, the Abel inversion method is no longer valid, due to the coexistence of two stellar images along Triton’s limb (see Fig. B.2). This problem arises at altitude levels of about 20 km, corresponding to the deepest layers probed by the light curve obtained at La Palma. Consequently, the central flash allows us to gain about 12 km downwards (about 0.6 scale height) compared to the inversion method. The explanation for this is shown in Fig. B.1 and further detailed in Appendix B.2.

Fig. 23 Surface pressure cycle on Triton as simulated with the VTM assuming a different fixed N2 ice distribution in both hemispheres. A thermal inertia of 1000 J s−1/2 m−2 K−1 (SI) was assumed. Bertrand et al. (2022) include more simulations in their paper, and their Fig. 9 shows cases with different thermal inertias. It is of note that their simulations show that a lower thermal inertia would delay the peak of the surface pressure, in opposition with the occultation data points. The blue lines refer to a southern cap extending to the equator, while the pink lines are for a southern cap extending to 30° S. Each line, marked with its corresponding value, refers to a different extension of the northern cap: 45° N, 60° N, 75° N, and no cap.

7.1 The central flash: observations

The central flash swept Europe along the lines shown in grey in Fig. 4. Among the 90 light curves shown in Figs. C.1–C.2, 42 show evidence of a stellar flux increase near mid-occultation, and 23 of them have enough S/N to be used in the central flash modelling.

Figure 24 displays the reconstructed intensity map of Triton’s shadow, with in particular the presence of a bright dot (central flash) near the shadow centre.

At Calar Alto, which passed at about 300 km from the shadow centre at closest approach (C/A), the increase in stellar flux is barely noticeable (Fig. 3), while it reaches the full unocculted stellar flux at Calern, which passed at 29 km from the centrality. At Constância (C/A 8.4 km), the maximum of the flash peaks at three times the unocculted stellar flux, and about 3.4 times the unocculted flux at Le Beausset (C/A 6.7 km, the closest of all stations; see Figs. 3 and 7).

The fit of the central flash is described in Sect. 4.2, except that we now allow a departure from sphericity of the layer responsible for the flash (see Sicardy et al. 2006 for details). We note that the ray-tracing code accounts for both the primary and secondary stellar images. Thus, we are not restricted in using this code as would be the case for the Abel inversion scheme (Appendix B.3).

7.2 The central flash: spherical fit

Assuming a spherical flash layer, we obtain the best simultaneous fits (now including the central flashes) displayed in Fig. 3. The quality of the fit is comparable to the quality obtained without the flashes , showing that no departure from sphericity is detected. A more quantitative assessment for the upper limit for such a departure (some 1.5 km along the limb, as projected in the sky plane) is provided in the next subsection. A closer visual examination of the residuals for the strongest flashes with best S/N reveal, however, some minor and localised features, possibly due to atmospheric waves (Fig. 7), but no global departure from the spherical model.

There is another argument supporting the spherical nature of Triton’s atmosphere. The centre of Triton’s shadow, as determined by the simultaneous fit to all the flashes, while excluding the ingress and egress parts of the light curves, coincides to within 0.1 km with the shadow centre determined by a global fit to all 52 light curves, by excluding the central flashes but including the ingress and egress parts. This 0.1 km offset is not significant, considering that the global fit centre has a typical 1σ error of 1 km cross-track (Fig. 12). In other words, the centre of the central flash layer, which is sensitive to the 8 km altitude level, coincides with the global shadow centre, which is sensitive to the 60 km altitude level. It could be that both atmospheric levels are close to spherical, but displaced in the same way with respect to Triton’s centre, but this configuration seems unlikely.

Fig. 24 Stellar flux in Triton’s shadow for the 5 October 2017 event, normalised to unity outside the body (light blue region). The flux reaches a minimum of about 7% of the unocculted flux inside the shadow and then rises sharply at the shadow centre. The direction of Triton’s rotation is indicated, as well as the equatorial rotation velocity, vrot = 17 m s−1 (in an inertial frame), using the parameters of Table 2.

7.3 The central flash: Limit on atmospheric distortions and winds

We now assess a possible departure of Triton’s lower atmosphere from sphericity, restricting ourselves to the simple model of a globally oblate flash layer. Testing more complex shapes will be performed once Triton’s 3D GCMs are available, something that is beyond the scope of this paper. Once projected in the sky plane, an oblate layer appears as an ellipse with apparent semi-major and semi-minor axes a′ and b′, respectively. The centres of curvature of that ellipse form a diamond-shaped caustic curve where abrupt flux variations are observed (see examples in Fig. 25). The equation of the caustic is (Elliot et al. 1977): (4)

where Oxy is a Cartesian reference system whose origin O is fixed at the ellipse centre, and where Ox (resp. Oy) is aligned with a′ (resp. b′). Since the flash layer lies at ~8 km altitude, we have a′ ~ 1360 km. We note that the orientations of the a′ and b′ axes are still to be specified.

We define the apparent oblateness of the flash layer as . We have explored values of ϵ′ from zero (i.e. a spherical flash layer) to some maximum value, and tracked the corresponding variations of χ² stemming from a simultaneous fit to central flashes. In this study, only the central flashes of the Varages, Calern, Constância, Le Beausset, and Felsina Observatory stations have been considered. The other stations are too far away from centrality and/or with lower quality to usefully constrain ϵ′. This is because the four cusps of the diamond-shaped caustic curve extend up to ~2ϵ′a from the shadow centre according to Eq. (4). As we obtain upper limits of ~0.002 for ϵ′ (see below), we have 2ϵ′a <~ 5 km for a ~ 1360 km. Thus, only the immediate vicinity of the shadow centre (typically less than 20 km) is sensitive to departures from sphericity. More distant stations essentially probe flashes that are indistinguishable from a spherical solution.

We first assume that the semi-minor axis b′ is aligned with Triton’s pole. This corresponds to an oblate flash layer maintained by an axisymmetric zonal wind regime that has a constant angular velocity around that axis. By using the (resp. criterion, we find lcr-level (resp. 3σ-level) upper limits of

for the apparent oblateness of the flash layer. The flash intensity map corresponding to this limit is displayed in Fig. 25. This apparent oblateness must be ‘deprojected’ to obtain the actual oblateness, ϵ, through the relation

where B = 40.5° S is the sub-observer latitude (Table 2), and where the approximation holds for . Using ϵ′ < 0.0011, this yields a lσ-level upper limit ϵ < 0.0019 for the deprojected oblateness. This corresponds to a difference between the equatorial and polar radii r_e and r_p of the layer, respectively, of r_e − r_p ~ 3 km, using r_e=a′ = 1360 km.

We assume that the flash layer shape is entirely supported by zonal winds. In particular, we assume the absence of a horizontal temperature gradient, so that the isobar level also corresponds to the isopyenic (constant density) layer. The radius r of the flash layer is given as a function of the latitude φ by the equation (Hubbard et al. 1993; Sicardy et al. 2006)

where ƒ = rv² (φ)/GM cos²(φ). This equation states that the isobar is locally perpendicular to the effective gravity field, where both the gravity field of the (spherical) body and centrifugal forces are accounted for. Introducing in that expression the polar equation of an oblate flash layer,

we obtain to lowest order in ϵ = (r_e − r_p)/r_e the velocity (5)

the value of GM_T is listed in Table 2. Using ϵ < 0.0019, this provides a 1σ -level upper limit for the zonal wind at the equator of .

We note that this motion can be prograde (positive sign) or retrograde (negative sign), and that it is measured in an inertial frame. Thus, noting the zonal wind in a frame rotating with Triton (which thus measures the atmospheric circulation at that level), we have , where v_rot = 17 m s⁻¹ is the equatorial velocity stemming from Triton’s rotation (see Fig. 24). Consequently, a retrograde zonal wind regime implies a 1σ limit , while a prograde regime implies . Those values are respectively 87 m s⁻¹ and 53 m s⁻¹ if the 3σ upper limit is considered.

Information on the atmospheric circulation of Triton was obtained in 1989 by V2: while surface wind streaks suggested eastward surface winds between latitudes of 15°S and 45°S (Hansen et al. 1990), the deflection of plumes showed that in the atmosphere above, at 8 km near 49°S and 57°S, the wind was westward and prograde (Hansen et al. 1990; Yelle et al. 1995). On the basis of theoretical consideration, Ingersoll (1990) proposed that this could result from a temperature contrast between the cold frost-covered pole and the warm un-frosted equator. More realistic GCMs simulations, including the N₂ condensation-sublimation cycle, have been reported in Vangvichith (2013), and additional relevant simulations have been performed to explore the circulation on Pluto, which is similar to Triton in terms of rotation rate and atmospheric composition (see Forget et al. 2021 and reference therein). These models show that if N₂ significantly sublimes in the southern hemisphere and condenses in the northern hemisphere, the circulation should be dominated by a retrograde circulation resulting from the conservation of angular momentum of the flow, with velocities that cannot be higher than the rotation of the planet (17 m s⁻¹). This is significantly less than the upper limits that we derive from our observations. In any case, as mentioned above, global retrograde winds were not observed in 1989: to get prograde rotation in the mid southern latitude as suggested by the V2 plume observations, the inter-hemispheric condensation must be weak. In that case a thermal gradient could create a weak prograde wind as suggested by Ingersoll (1990), reaching a few metres per second in GCM simulations. However, modelling Pluto suggests that in some conditions a regime of super-rotation (like on Venus or Titan) could occur (Forget et al. 2017). This could explain the plume direction on Triton. Such a super-rotation is thought to initially result from the formation of a high-mid-latitude jet (due to thermal balance between a warm equator and a colder pole, or condensation flow from low latitudes to the pole). Then barotropic waves can transport angular momentum to and from the equator and seriously accelerate the entire atmosphere. In their Pluto GCM, Forget et al. (2017) found mean equatorial zonal wind up to 15 m s⁻¹. However, this could be model-dependant. It is not easy to set a theoretical limit to such a super-rotation. Our upper limit on prograde wind near 50 m s⁻¹ provides a constraint for such a hypothetical super-rotation in 2017.

We have considered other orientations for the central flash layer, as projected in the sky plane, by relaxing the condition that b′ should be aligned with Triton’s pole. This might be the case if other causes of distortion than zonal winds are at work, for example tidal forces from Neptune or Triton’s potential anomalies. Those orientations provide more stringent upper limits of the apparent oblateness ϵ′ because the cusps of the caustic can then get closer to the paths of the central-most stations (Fig. 25). For instance, rotating the caustic by 45 degrees imposes the more stringent 1c upper limit, ϵ < 0.00074 (instead of 0.0019). This requires an equatorial wind of ~40 m s⁻¹, which is still quite a bit larger than the values ~10 m s⁻¹ expected from GCMs, and thus not a constraining limit as far as GCMs are concerned.

We now compare our upper limit for the deprojected oblate-ness of the central flash layer (ϵ < 0.0024) with the value obtained by Elliot et al. (1997) from a single cut inside the central flash region during the 14 August 1995 occultation. Two solutions are considered by those authors, an oblate flash layer with ϵ = 0.042 and a prolate one with ϵ = − 0.032. Adopting the ϵ = 0.042 value, we obtain the central map displayed in the right panel of Fig. 25. This would imply the crossing of the caustic by the five stations shown in the map, and result in strong flux variations at those crossings that are not observed in the data (see Fig. 7). A similar conclusion would be drawn by adopting the prolate value ϵ = − 0.032 of Elliot et al. (1997). We note that the closest approach distance to the shadow centre during the central observation flash of 1995 was about 100 km, which is well outside the diamond-shaped caustic displayed in the right panel of Fig. 25. Thus, caustic crossings could not be tested at that epoch.

The main problem of the large oblateness values obtained by Elliot et al. (1997) is that they imply unrealistically large wind velocities to maintain such distortions. For instance, taken at face value, ϵ = 0.042 results in an equatorial wind velocity of , more than twice the sonic velocity near Triton’s surface (~130 m s⁻¹ at ~40 K) and much larger than predicted by GCMs (see above). A possibility considered by Elliot et al. (1997) was that Triton’s atmospheric distortion was restricted to mid-latitude regions (i.e. it was local rather than global). This permits lower values of wind speeds (110−170 m s⁻¹) depending on whether prolate or oblate solutions are considered. This is still essentially supersonic and not expected from circulation models. Moreover, this distortion would also be detected in our dataset, which densely sampled the central flash region.

An alternative explanation proposed by Elliot et al. (1997) is the presence of hazes that absorbed part of the stellar flux at some specific locations along Triton’s limb, thus altering the central flash shape. However, according to those authors, neither the optical depth obtained for those hazes at the corresponding altitude levels, nor its dependence with wavelength, were consistent with the V2 results. The haze problem is discussed in more detail in Sect. 8.

Finally, the shape of Triton’s solid body as observed in V2 images indicates an oblateness smaller than 0.0014 (Thomas 2000), which is too small to explain the claimed distortions. Moreover, non-radial components of Triton’s inner gravitational field might also cause atmospheric distortions, but they would then be permanent and thus, should be observed also in 2017. From all this, we conclude that the large oblatenesses reported by Elliot et al. (1997) are both theoretically unexpected and inconsistent with our observations.

At this point, we think that the mismatches between the flash models and its observation in 1995 (as well as small departures from our spherical model in Fig. 7) could be caused by small local corrugations of the flash layer induced by gravity waves. As the caustic is the locus of the limb centre of curvatures, its shape is very sensitive to local (but small) corrugations of that layer. Examples of such effects have been investigated for explaining flash shapes in the cases of stellar occultations by Neptune (Hubbard et al. 1988) and Titan (Sicardy et al. 2006). As mentioned earlier, this approach remains beyond the scope of this paper, as long as zonal wind and gravity wave regimes are not available for Triton’s lower atmosphere.

Fig. 25 Maps of the central flash intensity. Left panel: map of the central flash intensity, adopting the 1σ- upper limit ϵ′ = 0.0011 for the apparent oblateness of the flash layer (corresponding to a deprojected oblateness ϵ = 0.0019) (see text for details). The black labels along the iso-intensity contours (in curves) indicate the received stellar flux in units of its unocculted value. The grey diamond-shaped feature near the centre is the caustic curve described by Eq. (4) and corresponding to ϵ′ = 0.0011. In the case considered here, the flash layer is assumed to be aligned with the apparent direction of Triton’s poles, indicated by the dash-dotted line. Neptune’s direction is determined from the position angle 286° of Triton with respect to the planet at the moment of the occultation. Right panel: same, but adopting the oblate solution ϵ = 0.042 found by Elliot et al. (1997) from the shape of a central flash observed during the 14 August 1995 occultation.

8 Pending issues

8.1 Hazes

In the present work, we assumed that Triton’s atmosphere is clear (i.e. free of absorbing material). However, V2 observations in 1989 revealed two kind of absorbing features in the lower atmosphere: hazes and clouds. These features were detected in the visible through the imaging system (Rages & Pollack 1992) and during the UV occultation experiment (Krasnopolsky et al. 1992, 1993; Krasnopolsky & Cruikshank 1995).

The general picture that emerges from the V2 observations is reviewed in Yelle et al. (1995) and is described as follows: hazes are detected up to an altitude of about 30 km, they were observed around the entire Triton’s limb, except for a small clear region near east longitude 280° and between latitudes 4 and 18° S. Thus, unless drastic changes in haze formation occurred, they should also be present during our observation of 5 October 2017. Clouds are observed closer to the surface compared to hazes (i.e. below an altitude level of about 8 km). Contrary to hazes, they exhibit a patchier distribution along the limb.

Due to their low altitudes, absorbing materials should be best detected in the central flash structures, caused by a layer at about 8 km altitude. Estimation of the integrated (down to the surface) vertical optical depth of the hazes at 0.47 µm is τ_vis = 0.005 ± 0.001, with a typical scale height of H_h ~ 12 km (Krasnopolsky et al. 1993; Rages & Pollack 1992). Thus, the integrated vertical optical depth down to the 8 km altitude level should be reduced by a factor of exp(−8/H_h) ~ 0.5. Moreover, the slant optical depth (along the line of sight) is amplified by a factor of , yielding a slant optical depth of hazes at 8 km of the order of 0.005 , and a reduction in the flash amplitude by 5–10% if no changes have occurred since 1989.

Clouds are much denser than hazes, with vertical optical depth to the surface of 0.1 or higher. Cloud particles have a vertical distribution with a scale height H_c comparable to or larger than the atmospheric scale height (about 20 km). Thus, the slant optical depth of clouds at 8 km should be in the order of 0.1 or larger, corresponding to a decrease in the flash amplitude by a factor of 4 or more.

To summarise, hazes are expected to have a mild effect on the central flash heights, with an expected reduction of only 5– 10%. We do not see any departure from the model at that level on the best flash profiles. Moreover, we do not detect trends on the observed flash amplitudes versus wavelength. For instance, the strong flash at Constância (Fig. 7) observed at an effective wavelength of ~0.6 µm agrees with the model at the same satisfaction level as the Varages flash, observed at ~1.3 µm (Fig. 3). The same is true with the Calar Alto observation, which was made simultaneously in the visible (0.4–1.0 µm) and the near IR (1.0–1.7 µm) (see Fig. 3 and Table A.1). A flash is observed at that station, but it is too faint to reveal a difference between the two channels. Finally, the dual observation made at Kryoneri (R and I bands) was too far away from the centre line to show a central flash, but does not show significant differences anyway in the fits by the synthetic light curve (Fig. C.1).

Taken at face value, these results indicate that hazes have no detectable effect on the flash shapes. This result is true a fortiori for the clouds, which would reduce the flash amplitudes by a factor of 4 or greater; this is not observed in the data. In other words, our model consistently explains all the flashes using a clear atmosphere.

However, a difficulty arises at this point since the height of the flash actually depends on the assumed template temperature profile. In order to disentangle the haze versus temperature effects, independent information on the thermal profile of the lower atmosphere is required. This requirement could be met with the ALMA results reported by Gurwell et al. (2019). This issue remains, for the moment, beyond the scope of the present paper.

8.2 Troposphere

Yelle et al. (1991) inferred a troposphere from the V2 observations of geysers and clouds. We do not observe this, as our deepest layer probed, at the central flash level, which coincides with the expected altitude of the tropopause. Therefore, our model is consistent with our central flashes, as discussed in Sect. 7, and there is no need to include a troposphere to accommodate for the data. However, this does not mean that we can exclude a troposphere as we do not have information down to this part of the lower atmosphere.

The absorption, caused by hazes and clouds, impacts our data and the model. For a troposphere to be included, we need independent measurements of the lower atmosphere.

8.3 Gravity waves

In Sect. 7 and Fig. 7 we mention that the residuals in these light curves show minor features that are probably attributed to atmospheric waves. They present themselves as small fluctuations, with little effect on the overall light curve shape, and therefore we did not study them in detail. This topic needs further analysis that is beyond the scope of the present paper.

8.4 ϕ₀ of inverted profiles

The determinations of the baseline levels ϕ₀ used for inverting the profiles of La Palma and Helmos remain an open issue, as they provide inconsistent results in the deepest parts of the T(r) profiles (Appendix B.2). This matter must be analysed further and needs independent results, in particular from ALMA, to confirm which ϕ₀ should be used.

9 Concluding remarks

In this paper we present results obtained from the ground-based stellar occultation by Triton observed on 5 October 2017. The main goals were (i) retrieve the general structure of Triton’s atmosphere between altitude levels of ~8 km (~9 µbar) and ~190 km (few-nanobar level) and (ii) compare these results with other ground-based occultations and the V2 radio occultation to assess the pressure seasonal variation over the last three decades.

The 2017 event yielded 90 positive observations, 42 of which showed a central flash. We used Abel inversions to retrieve density, pressure, and temperature profiles from our three best S/N light curves. We find a hint of a mild negative temperature gradient (reaching −0.2 K km⁻¹) at the deepest part of our profiles (i.e. below the altitude of ~30 km). This constitutes a mesosphere just above an expected stratosphere with a positive temperature gradient that connects the atmosphere to the cold surface.

A ray-tracing approach was used for a global fit to the best 52 light curves, providing a pressure p₁₄₀₀ = 1.18 ± 0.03 µbar at radius 1400 km. It also provides a synthetic and smoothed model to extrapolate the density, pressure, and temperature down to the surface.

A new analysis of the V2 radio experiment, with useful information extracted from the surface up to around 1400 km, shows that the pressure retrieved in the 2017 event is consistent with the pressure obtained in 1989. A survey of pressure values obtained between 1989 and 2017 was conducted. The two past occultations (in 1997 and 2008), reanalysed using our methods, indicate that the surface pressure reported in the 1990s is real, but this remains debatable due to the scarcity of high S/N light curves and the lack of a fully consistent analysis of the best datasets used by other teams.

The pressure that we obtain from the 2017 occultation is consistent with that derived from the Voyager radio experiment, meaning that the pressure is back to its 1989 level. Results from a VTM of Triton, described in detail in Bertrand et al. (2022), do not support a strong increase in surface pressure in the last few decades but instead a modest increase, with a surface pressure reaching up to 20 µbar in the intervening 28 yr. The VTM simulations also suggest that (1) a strong increase in surface pressure before 2000 cannot be obtained if N2 is present between latitude 30°S - 0°, and (2) a northern polar cap should have extended down to at least 45° N - 60° N in 2017 to have the surface pressure back at the V2 level from 1989.

Finally, the central flash analysis does not reveal any evidence of an atmospheric distortion. The atmosphere appears as globally spherical, with a 1σ upper limit of 0.0011 for its apparent oblateness near the 8 km altitude. This corresponds to a global difference of less than 1.5 km between the largest and smallest atmospheric radii at that altitude. This is much smaller than values reported in the literature. In particular, this does not support the existence of supersonic winds previously claimed by Elliot et al. (1997).

Open issues, requiring a more specific analysis of our dataset, will be addressed elsewhere. These include the possible presence (or absence) of hazes and of a troposphere just above Triton’s surface, as well as the possible detection of gravity waves.

Acknowledgements.

J.M.O. acknowledges financial support from the Portuguese Foundation for Science and Technology (FCT) and the European Social Fund (ESF) through the PhD grant SFRH/BD/131700/2017. 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 n° 669416 “Lucky Star”. We thank S. Para who supported some travels to observe the 5 October 2017 occultation. T.B. was supported for this research by an appointment to the National Aeronautics and Space Administration (NASA) Post-Doctoral Program at the Ames Research Center administered by Universities Space Research Association (USRA) through a contract with NASA. We acknowledge useful exchanges with Mark Gurwell on the ALMA CO observations. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. J.L.O., P.S.-S., N.M. and R.D. acknowledge financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709), they also acknowledge the financial support by the Spanish grant AYA-2017-84637-R and the Proyecto de Excelencia de la Junta de Andalucía J.A. 2012-FQM1776. The research leading to these results has received funding 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). P.S.-S. acknowledges financial support by the Spanish grant AYA-RTI2018-098657-J-I00 “LEO-SBNAF”. The work was partially based on observations made at the Laboratório Nacional de Astrofísica (LNA), Itajubá-MG, Brazil. The following authors acknowledge the respective CNPq grants: F.B.-R. 309578/2017-5; R.V.-M. 304544/2017-5, 401903/2016-8; J.I.B.C. 308150/2016-3 and 305917/2019-6; M.A. 427700/20183, 310683/2017-3, 473002/2013-2. This study was financed in part by the Coor-denação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001 and the National Institute of Science and Technology of the e-Universe project (INCT do e-Universo, CNPq grant 465376/2014-2). G.B.R. acknowledges CAPES-FAPERJ/PAPDRJ grant E26/203.173/2016 and CAPES-PRINT/UNESP grant 88887.571156/2020-00, M.A. FAPERJ grant E-26/111.488/2013 and A.R.G.Jr. FAPESP grant 2018/11239-8. B.E.M. thanks CNPq 150612/2020-6 and CAPES/Cofecub-394/2016-05 grants. Part of the photometric data used in this study were collected in the frame of the photometric observations with the robotic and remotely controlled telescope at the University of Athens Observatory (UOAO; Gazeas 2016). The 2.3 m Aristarchos telescope is operated on Helmos Observatory by the Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing of the National Observatory of Athens. Observations with the 2.3 m Aristarchos telescope were carried out under OPTI-CON programme. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 730890. This material reflects only the authors views and the Commission is not liable for any use that may be made of the information contained therein. The 1.2 m Kryoneri telescope is operated by the Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing of the National Observatory of Athens. 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 0.81 m Main Telescope at the OAVdA was upgraded thanks to a Shoemaker NEO Grant 2013 from The Planetary Society. D.C. and J.M.C. acknowledge funds from a 2017 ‘Research and Education’ grant from Fondazione CRT-Cassa di Risparmio di Torino. P.M. acknowledges support from the Portuguese Fun-dação para a Ciência e a Tecnologia ref. PTDC/FISAST/29942/2017 through national funds and by FEDER through COMPETE 2020 (ref. POCI010145 FEDER007672). F.J. acknowledges Jean Luc Plouvier for his help. S.J.F. and C.A. would like to thank the UCL student support observers: Helen Dai, Elise Darragh-Ford, Ross Dobson, Max Hipperson, Edward Kerr-Dineen, Isaac Langley, Emese Meder, Roman Gerasimov, Javier Sanjuan, and Manasvee Saraf. We are grateful to the CAHA, OSN and La Hita Observatory staffs. This research is partially based on observations collected at Centro Astronómico Hispano-Alemán (CAHA) at Calar Alto, operated jointly by Junta de Andalucía and Consejo Superior de Investigaciones Científicas (IAA-CSIC). This research was also partially based on observation carried out at the Observatorio de Sierra Nevada (OSN) operated by Instituto de Astrofísica de Andalucía (CSIC). This article is also based on observations made with the Liverpool Telescope operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council. Partially based on observations made with the Tx40 and Excalibur telescopes at the Observatorio Astrofísico de Javalambre in Teruel, a Spanish Infraestructura Cientifico-Técnica Singular (ICTS) owned, managed and operated by the Centro de Estudios de Física del Cosmos de Aragón (CEFCA). Tx40 and Excalibur are funded with the Fondos de Inversiones de Teruel (FITE). A.R.R. would like to thank Gustavo Román for the mechanical adaptation of the camera to the telescope to allow for the observation to be recorded. R.H., J.F.R., S.P.H. and A.S.L. have been supported by the Spanish projects AYA2015-65041-P and PID2019-109467GB-100 (MINECO/FEDER, UE) and Grupos Gobierno Vasco IT1366-19. Our great thanks to Omar Hila and their collaborators in Atlas Golf Marrakech Observatory for providing access to the T60cm telescope. TRAPPIST is a project funded by the Belgian Fonds (National) de la Recherche Scientifique (F.R.S.-FNRS) under grant PDR T.0120.21. TRAPPIST-North is a project funded by the University of Liège, and performed in collaboration with Cadi Ayyad University of Marrakesh. E.J. is a FNRS Senior Research Associate.

Appendix A Circumstances of observations

Appendix B Retrieval of atmospheric structure

Here we discuss some limitations and caveats associated with the approach mentioned in the main text.

Appendix B.1 Upper and lower limits of probed atmosphere

The density n(r) is reliably retrieved up to the level where the flux standard deviation σ_Φ due to photometric noise is comparable to the drop of stellar flux caused by the occultation (DO15). This corresponds to a density of (B.1)

where K is the molecular refractivity of N₂, H = −n/(dn/dr) is the scale-height and D is Triton’s geocentric distance.

Using the values of K and D listed in Table 2, considering that H ~ 30 km around r = 1500 km, and taking σ_Φ ~ 0.011 for our best dataset (La Palma station), we find that reliable density values cannot be obtained above r_upper ~ 1540 km, where n ~ 4 × 10¹¹ cm⁻³, corresponding to pressures of a few nanobars.

The relative error on n(r) reduces as exp[−(r − r_up)/H] when deeper levels are probed. In practice, only the inversion of the best light curves provide useful retrieved profiles. In our case, this concerns the light curves from La Palma and Helmos stations, providing the profiles displayed in Figs. 9-11.

The deepest layers probed when inverting an occultation light curve are those corresponding to the closest approach of the observing station to Triton’s shadow centre. In our case, La Palma’s light curve provides data down to a radius of about r = 1375 km (i.e. just above the 20 km altitude level). However, the central flashes provide constraints that go further down, typically just under the 8 km altitude level (see Sect. 7).

Appendix B.2 Constructing the temperature profile template

A difficulty encountered during the inversion of the light curves is the assessment of Triton’s contribution to the total flux. We define ϕ₀ = F_T/(F_S + F_T), where F_T (resp. F_S) is the flux coming from Triton (resp. the unocculted star). Thus, ϕ₀ corresponds to the zero stellar flux level in the normalised occultation light curves. Changing its value mainly changes the deepest parts of the retrieved profiles (see Fig. B.1).

Measuring ϕ₀ relies on images where Triton and the star are angularly separated. It is generally a difficult task, since a photometric accuracy of better that 1% is necessary to bring useful constraints. To mitigate differential chromatic effects (as the star and Triton have different colours), images were taken at the same elevation as for the occultation, either during the same nights or during nights before or after the event. As a sanity check, it is also desirable to use another reference star with flux F_R and see if the sum of the ratios F_S/F_R and F_T/F_R outside the event matches the ratio (F_S + F_T/R_R during the event. Mismatches may then reveal possible variabilities in any of the objects involved (the target and reference stars, and/or Triton), and serve as an estimator of systematic sources of errors.

Such calibration images were acquired at the La Palma and Helmos sites, which provided the best datasets in terms of S/N. The focal lengths of those two telescopes are large enough to clearly resolve Triton from Neptune and from the occulted star and to avoid, in particular, flux contamination from the planet in the occultation light curve.

The calibration results are self-consistent for both instruments, with small internal error bars on ϕ₀ (i.e. ϕ₀ = 0.3445 ± 0.0003 for La Palma and ϕ₀ = 0.360 ± 0.013 for Helmos). We note that there is no reason why the ϕ₀ should be the same for the two stations, as different filters were used: I+z (>720 nm) at La Palma and V+R at Helmos. Also, different observing conditions were prevailing at the two stations. In particular, the different airmasses (1.3 and 2.4 at La Palma and Helmos, respectively) may also affect the flux ratios star/Triton. However, the different values of ϕ₀ at La Palma and Helmos are mutually inconsistent (to within error bars) in the sense that they provide significantly different bottom parts for the T(r) profiles.

This is particularly problematic for La Palma, where we obtain a temperature profile that peaks at T ~ 58 K just above the flash layer, using the value ϕ₀ = 0.3445 derived from the calibration (Fig. B.1). This imposes a strong inversion layer in order to connect the temperature profile 38 K at the surface (see Fig. 15 and associated discussion). When used in our ray-tracing code, this profile yields inconsistent synthetic central flashes when compared with observations (Fig. 3). In particular, it is not possible to fit simultaneously the flashes lying north and south of the shadow centre.

We have no satisfactory explanations for the inconsistencies induced by the photometric calibration (i.e. by the retrieved values of ϕ₀). It could stem from unaccounted light contamination from Neptune, although the biggest telescopes should be free of this problem, as mentioned above. We tested the effect of digital coronagraphy to remove any such possible contamination. However, no significant changes of ϕ₀ are obtained with and without the use of coronagraphy. Possible variations of Triton’s flux due to rotational light curves of the satellite have also been considered. However, considering the low amplitude of those rotational curves (Buratti et al. 2011), and the fact that calibrations were made about 90 minutes before the occultation at La Palma, Triton’s flux variation should be well below the 1% level, and thus have a negligible effect on the retrieved value of ϕ₀.

In these conditions, we have opted for another approach. We have varied the value of ϕ₀ for the best dataset (i.e. the La Palma light curve). For each value, we have inverted the light curve to obtain T(r), and then derived a template temperature profile using the Dias-Oliveira et al. (2015) modelling, except that the upper branch of T(r) is not isothermal, but has a constant gradient dT/dr ~ 0.1 K km⁻¹ to account for the general temperature increase (thermosphere) described in Strobel & Zhu (2017).

Each particular value of ϕ₀ provides a set of prescribed T(r) boundary conditions (the green dots in Fig. B.1). Then it is possible to interpolate those prescribed values for any ϕ₀ using smooth polynomial functions, and finally get a one-parameter family of T(r) profiles depending only on ϕ₀.

The main goal here is to find the value of ϕ₀ that best matches all the central flashes, thus constraining the deepest part of the T(r) profile. The best T(r) profile (with ϕ₀ = 0.35885) is shown as a thicker line in the upper panel of Fig. B.1, and is adopted throughout this paper every time we use our ray-tracing code. An example of an output of our tracing code is displayed in the lower panel of Fig. B.1.

Using the notations of Dias-Oliveira et al. (2015), the T(r) profile is constructed by adopting the parameters of Table B.1. We note that although both the T and dT/dr profiles are continuous, the d²T/dr² profile is not. This is evident in Fig. 11 at the inflection point labelled ‘3’ in Fig. B.1. The discontinuity of d²T/dr² creates a very small kink at the corresponding point in the synthetic light curve (see the lower panel of Fig. B.1). This kink is well below the noise level of all the observed light curves, and thus, has a negligible effect on the fit to the data.

Fig. B.1 Template temperature profiles and result of the ray-tracing code. Upper panel: Six template temperature profiles, indicated by the thin black lines, obtained taking values ϕ0 = 0.370, 0.365, 0.360, 0.355, 0.350, and 0.3445 (from left to right) of Triton’s contribution to the normalised occultation light curve at La Palma. The rightmost profile (corresponding to ϕ0 = 0.3445) is the profile obtained from the calibration images taken before the event. Exploring ϕ0 with small incremental steps, we find a best fit to the central flash for ϕ0 = 0.35885 (see the thicker black template profile, where the green dots labelled 1 to 4 are prescribed particular points in the Dias-Oliveira et al. 2015 model, as specified in Table B.1). We note that all the profiles go through the boundary condition T = 54.5 K at r = 1453 km (100 km altitude, upper green dot). The inverted profiles for La Palma corresponding to ϕ0 = 0.35885 are shown in red (ingress) and blue (egress). The horizontal dotted line marks the altitude of the central flash layer. Lower panel: Result of the ray-tracing code, using the best T(r) profile of the upper panel and the best fit value of p1400 = 1.18 µbar (Table 4). The normalised stellar flux is plotted against the distance to Triton’s shadow centre. The green curve represents the flux if only one stellar image is present, and the green dots show the correspondence with the points labelled 1 to 4 in the upper panel. The black curve is the sum of the fluxes from two stellar images, i.e. the sum of the green curve and its mirrored version with respect to the z = 0 axis. This black curve is then used to fit the observations.

Admittedly, there is not a unique way to find a template T(r) model that best fits the flashes. However, our solution should capture the main properties of the real profile, with a temperature maximum reached just above the central flash layer, and a mesosphere with a mild negative temperature gradient above that temperature peak. We note that the inversion layer connecting the profile to the surface at 38 K (i.e. below the temperature peak) has a vanishing effect on the synthetic light curve as the surface is approached, thus defining a ‘blind zone’ as far as our data are concerned.

Appendix B.3 The secondary stellar image issue

The Abel inversion assumes that only one stellar image (the near-limb, or primary refracted image) contributes to the recorded flux. In reality, Triton’s atmosphere also produces a far-limb (secondary) stellar image whose flux is added to the light curve⁶. This is a source of error, especially in the central flash region, where the primary and secondary images have comparable fluxes (see for instance the lower panel of Fig. B.1).

We have performed tests to assess this effect on La Palma’s light curve. At closest approach, the station was at about 685 km from the shadow centre, providing information down to the ~20 km altitude level (Fig B.2). A smooth T(r) profile is used to generate synthetic light curves at this station, one with the flux of the primary image only, and one with the sum of the fluxes of the primary and secondary images.

As expected, the inversion of the one-image light curve correctly retrieves the temperature at the 0.2 K accuracy level (Fig B.2), and the density and pressure profiles at the 0.1% accuracy level. Conversely, the inversion of the two-image light curve correctly retrieves the upper parts of the profiles, but it fails in reproducing the lower parts. For instance, at the deepest point reached by La Palma’s light curve, the primary flux is 0.061 (normalising the unocculted stellar flux to unity), while the secondary flux is 0.0072, about 8.5 times fainter than the primary flux. At that point, the temperature is retrieved to within 0.9 K and the pressure at the 1% level, a satisfactory result at our accuracy level.

However, this discrepancy rapidly increases as deeper levels are probed. For instance at r = 1362 km (9 km altitude level, where the temperature locally reaches a maximum; see the upper panel of Fig. B.2), the discrepancy between the original and retrieved temperature is about 3 K. We also note that in this case the retrieved temperature profile has an unrealistic behaviour as it extends below Triton’s surface. Although the pressure profile is satisfactorily retrieved even if the secondary image is present, it suffers nevertheless the same unrealistic behaviour as T(r), as it also extends below Triton’s surface (lower panel of Fig. B.2).

In summary, the inversion procedure cannot provide reliable results below the 20 km altitude level. In particular the central flash region cannot be used in the inversion procedure. We note, however, that the direct (ray-tracing) approach does include the primary and secondary fluxes, and as such, it can be used to constrain the atmospheric profiles in the central flash region.

Fig. B.2 Effect of the secondary stellar image on the inversion results. Upper panel: Effect on the temperature profiles. The smooth black line connects these profiles to the surface at 38 K, where the pressure is set to 14 µbar. The resulting density profile is used to generate synthetic light curves that feed the Abel inversion procedure. The solid green line is the retrieved T(r) profile with only the primary image accounted for. The green dotted line is the retrieved profile where the primary and secondary images have been added. The profiles obtained from the inversion of La Palma’s light curve at ingress and egress (in red and blue, respectively) are shown for comparison (see also Figs. 9– 11). Black horizontal solid and dotted lines show Triton’s surface and reference radius (1400 km). Lower panel: Same, but with the pressure profiles.

Appendix C Fit to the data

Fig. C.1 Data (black dots) fitted simultaneously with synthetic light curves (blue lines), based on the temperature profile displayed in Fig. 10 (black line) and the pressure boundary condition p1400 = 1.18 µbar (Table 4). The green dots are the residuals of the fits. We note that the central flash regions have been excluded from the fit so that we obtain a global fit that is not influenced by the deepest atmospheric layers. The stations with exposure times smaller than 1 s have been smoothed to have a sampling time close to 1 s, allowing a direct visual comparison of S/N of the various datasets. The lower and upper horizontal dotted lines mark the zero-flux level and the total star plus Triton unocculted flux, respectively. We note that the three central-most stations (Constância, Le Beausset, and Felsina observatories in Fig. C) are plotted at a different vertical scale to accommodate the presence of a strong central flash. Each panel has a duration of five minutes and is centred around the time of closest approach (or mid-occultation time) of the station to Triton’s shadow centre. The stations are sorted from left to right and top to bottom from the northernmost track (Newark) to the southernmost track (Athens; see Fig. C), projected on Triton in the sky plane (Fig. 5). For reference, the vertical red lines mark 23:48 UTC for the European and African stations, and 23:55 UTC for the US stations (Newark, Ithaca, and Dark Sky observatories).

Fig. C.1 Continued. NB. ‘JAVA.’ is the abbreviation of Javalambre, used so that the name of the station fits into the plot.

Fig. C.2 The same as Figs. C.1, but for stations that were not used in the simultaneous fit.

References

All Tables

All Figures

Fig. 1 Sub-solar latitude on Triton over time. Upper panel: overview of the sub-solar latitude on Triton versus time over the last millennium. The blue part corresponds to the period from the V2 encounter (August 1989) to the 5 October 2017 stellar occultation. It shows that during this interval, Triton experienced an extreme summer solstice in its southern hemisphere, with a minimum sub-solar latitude of 50° S in 2000. Lower panel: close-up view of the upper panel around the year 2000. The black points correspond to occultations observed, distinguishing from the black triangle (V2 RSS experiment). The larger symbols are from the data that we use in this paper.

In the text

	Fig. 2 Average offsets in (α, δ) mas of the eight nights of OPD observations using the preliminary Gaia DR2 catalogue for the astro-metric reduction. Positions in lighter colours were not accounted for to compute the average offset, due to lower quality associated with poor weather conditions. Lower panel: number of images taken for each night.
In the text

Fig. 3 Simultaneous fits of the data (black dots) by synthetic light curves (blue lines), based on the temperature profile displayed in Fig. 10 (black line) and the pressure boundary condition p1400 = 1.18 µbar (Table 4). The green curves are the residuals of the fits. The lower and upper horizontal dotted lines mark the zero-flux level and the total star plus Triton unocculted flux, respectively. We also note that the three bottom light curves are plotted at a different vertical scale from the others to accommodate the presence of a strong central flash. The stations are sorted from left to right and top to bottom from the northernmost track (St Caprais) to the southernmost track (Calar Alto; see next figure). Each panel has a duration of five minutes and is centred around the time of closest approach (or mid-occultation time) of the station to Triton’s shadow centre, as indicated under the lower left panel in each block of six light curves. For reference, the vertical red line marks the time 23:48 UTC. The stations with exposure times of less than 1 s have been smoothed to have a sampling time as close as possible to 1 s, for easier S/N comparison of the various datasets. We note that in this approach, the sampling of the Constância, Le Beausset, and Felsina Observatory light curves (0.64 s) is kept at its original value so that full resolution versions of the corresponding strong flashes at those stations are displayed here. The same kinds of plots showing all the stations, but with the flashes excluded from the fits, are displayed in Fig. C.1.

In the text

	Fig. 3 continued. NB. ‘JAVA.’ is the abbreviation of Javalambre, used so that the name of the station fits into the plot.
In the text

Fig. 4 Triton’s shadow path on 5 October 2017. The black dots are spaced by one minute, the arrow indicating the direction of motion of the shadow. The northern and southern limits of the solid body assuming a radius of 1353 km are also represented, with the predicted path as the white lines and the effective path as the blue lines. The grey area represents night on Earth (dark grey for full astronomical night and light grey for twilight). Stations with a successful observation that have been used in our fit are represented by blue dots, while the ones that were not used are shown as red dots. The white dots are the stations that attempted observation but were clouded out or had technical difficulties. Upper panel: overview of all observing stations. The larger black dot along the black line corresponds to the closest approach of the shadow centre to the geocentre, at around 23h52 UTC (see Table 1). Lower panel: closer look at the central flash path across Europe. The grey lines around the centremost line correspond to a spacing of 50 km (corresponding to about 2.5 mas) once projected in the sky plane.

In the text

Fig. 5 Geometry of the 5 October 2017 stellar occultation by Triton, as seen in the sky plane. The J2000 celestial north (N) and east (E) directions and the scale are indicated in the upper right corner. Triton’s radius is fixed at RT = 1353 km, and the grey arrow near the equator shows the direction of rotation of the satellite. The (Neptune-facing) prime meridian is drawn as a thicker line compared to the other meridians, and the south pole is marked by the label S. The inclined lines are the trajectories of the star relative to Triton (or ‘occultation chords’) as observed from various stations, with the black arrow indicating the direction of motion. We gathered a total of 90 occultation light curves, 52 of which (corresponding to the blue colour, as in Fig. 4) had sufficient S/N to be included in a global atmospheric fit; the remaining 38 (red colour) with lower S/N were not included in the fit.

In the text

Fig. 6 Best two light curves obtained during the Triton occultation of 5 October 2017, at La Palma and Helmos stations (see details in Table 1). Both telescopes were equipped with the same E2V CCD 47-20 detector with quantum efficiency peaking at 600 nm and reaching zero near 300 and 1000 nm, respectively. Left panel: full resolution light curve (cycle time 0.635 s) for the La Palma station. Right panel: same for the Helmos station (cycle time 0.674 s). The spectral ranges used for each instrument are indicated in the figures (I+z at La Palma and V+R at Helmos). The blue lines are the best simultaneous fits obtained with our ray-tracing approach (see Sect. 4).

In the text

	Fig. 7 More detailed views taken of the flashes observed at Calern (top panel) and Constância (bottom panel), with the same setup as in the previous figure.
In the text

Fig. 8 Tracks of the primary (red dots) and secondary (blue dots) stellar images above Triton’s surface, as seen from Constância, plotted every 0.1 s. The junctions between the red and blue paths correspond to ingress (left) and egress (right) points for the Constância station. The arrows show the direction of the stellar images’ path. The regions probed by the central flash are those where the dots are more spaced. All the other stations probed essentially the same path (or part of it), with the primary and secondary images being swapped (as well as their directions of motion), depending on whether the station probed north or south of the shadow centre. Since the Earth and the Sun are angularly close (-1°) to each other as seen from Triton, the stellar paths essentially mark Triton’s terminator, the night side extending above the terminator in this figure. The two yellow symbols are for La Palma station, with ingress plotted as a star and egress plotted as a diamond. The two green symbols are the same for Helmos station and the two white symbols for Calern station (see Table 3 for the corresponding values of the latitudes and longitudes). The background image is a global colour map of Triton, produced using V2 data and orange, green, and blue filter images in order to obtain an approximation of Triton’s natural colours. Background image credits: image selection, radiometric calibration, geographic registration and photometric correction, and final mosaic assembly were performed by Dr. Paul Schenk at the Lunar and Planetary Institute, Houston, Texas. Image data are from V2 (NASA, JPL).

In the text

Fig. 9 Density and pressure proflies of the atmosphere of Triton as a function of radius, r (the distance to Triton’s centre). Upper panel: density profiles of Triton’s atmosphere as a function of radius, retrieved by inverting three light curves obtained during the 5 October 2017 occultation and from the V2 radio phase delay at 3.6 cm. The colour codes are indicated in the upper right part of the plot. The same codes are used in Figs. 10–11 and Figs. 14–15. The thin black curve is a smooth synthetic density profile that fits the inverted profiles and is extrapolated down to the surface. It is derived from the smooth temperature profiles shown in Fig. 10. The solid horizontal line marks Triton’s surface (at radius RT = 1353 km), the dashed line indicates the central flash layer (near 1360 km), and the dotted horizontal line marks the reference radius, rref = 1400 km. Lower panel: corresponding pressure profiles.

In the text

	Fig. 10 Temperature profiles as a function of radius (upper panel) and pressure (lower panel). The oblique dotted line in the lower panel is the wet adiabat, i.e. the vapour pressure equilibrium line for N2, taken from Fray & Schmitt (2009).
In the text

	Fig. 11 Temperature gradient corresponding to the upper panel of Fig. 10. The dot-dashed line is the dry adiabatic temperature gradient Γ= − g/cp, i.e. the limit of convective instability, where cp is the specific heat at constant pressure for N2 and g = GMT/r2 is the acceleration due to gravity.
In the text

	Fig. 12 χ2 map for the simultaneous fit of 52 light curves obtained during the occultation of 5 October 2017, using the ϕ0 corresponding to the temperature profile in Fig. 10. The inner green line contour is the 1σ limit of the fit, while the outer green line is the 3σ limit.
In the text

Fig. 13 Radio phase delay during the egress of the V2 RSS occultation on 25 August 1989 in the 3.6 cm X band. Left panel: radio phase delay observed (adapted from G95). The crosses are the data, and the three solid lines (labelled B1, B2, and B3) are three polynomial modellings of the phase delay baseline. As discussed by G95, the preferred baseline solution is B2. Right panel: radio phase after subtraction of the B2 baseline polynomial shown in the left panel (green crosses), thus representing the effect of the atmosphere only. The green line marks the smooth version of that radio phase delay, as constructed by G95. The red profile is the phase delay that would be observed at 3.6 cm from the retrieved density profile of La Palma at immersion (see Fig. 9). Other phase delay profiles obtained from La Palma (emersion) and Helmos would be indistinguishable from the red profile and are not plotted here for sake of clarity. The black profile is the phase delay obtained from our best model of Triton’s atmosphere (see text for details).

In the text

Fig. 14 Close-up view of Fig. 9. Upper panel: close-up view of the density profiles. Various profiles derived from the RSS occultation are shown in green. The dashed line is the profile retrieved by Tyler et al. (1989) using a polynomial extrapolation B1 to correct for the RSS phase instability. The thin solid line is the profile retrieved by G95 using the polynomial extrapolation B2. The dotted line is the same by G95, but using the polynomial extrapolation B3. The thick solid line is the best model of G95, based on the B2 profile. Lower panel: same for the pressure profiles.

In the text

Fig. 15 Close-up view of the lower panel of Fig. 10. The width of each coloured box is Triton’s surface temperature (Tsurf), as estimated by the various authors mentioned just above the boxes (see text for details). The heights of the boxes are the range of surface pressures, psurf, using the ideal gas law , where is the surface molecular nitrogen density derived from our inversion of the RSS data (see text). We note that all the boxes intersect the vapour pressure equilibrium, which is plotted as a dotted line. This shows that the RSS surface density and the estimated surface temperatures are mutually consistent with a pressure being controlled by the N2 ice sublimation.

In the text

	Fig. 16 Geometry of the 18 July 1997 occultation, with the same conventions as in Fig. 5.
In the text

	Fig. 17 Same as Fig. 12, but for the 18 July 1997 occultation.
In the text

Fig. 18 Simultaneous fits to the 18 July 1997 light curves. The panel covers 300 s in time. All the light curves have been shifted so that the mid-occultation times are aligned. The red vertical tick marks indicate 10:10 UTC at the Brownsville station (USA) and 10:18 UTC for the three Australian stations. The blue lines are simultaneous fits to the data (black dots), using the best value found in Fig. 17, p1400 = 1.9 µbar, and the temperature profile shown in Fig. 10. The green dots are the fit residuals. For each light curve, the upper dotted line is the normalised value of the star plus Triton flux, and the lower dotted line is the background flux. We note that the data from the Brownsville station were normalised during the event, as shown in Elliot et al. (2000a).

In the text

	Fig. 19 Geometry of the 21 May 2008 occultation, with the same conventions as in Fig. 5.
In the text

	Fig. 20 Same as Fig. 12 but for the 21 May 2008 occultation.
In the text

	Fig. 21 Fit to the 21 May 2008 light curves. The same conventions as for Fig. 18 are used, except that the panel now covers 720 s in time. The synthetic light curve for Piton Lacroix is plotted in grey because it is not used in the fit, due to the high noise level. The red vertical tick marks indicate 01:51 UTC for Piton Lacroix and Maïdo, and 01:41 UTC for Hakos.
In the text

Fig. 22 Triton’s atmospheric pressure seasonal variations with time, using the values from Table 4. Our results are in red, and values taken from other works are plotted in blue. For better viewing, the value derived by Elliot et al. (2000a) from the 18 July 1997 occultation is plotted in a semi-transparent blue colour, using the same dataset that we are using here for that date, with our red diamond-shaped point just below it. For all points, the thick error bars correspond to 1σ confidence levels. We note that for the 4 November 1997 and 5 October 2017 values, the 1σ error bar is smaller than the diamond-shaped symbol and therefore is not visible. For the 18 July 1997 and 21 May 2008 events, we have also plotted for information our 3σ error bars as thinner lines (see text for discussion).

In the text

Fig. 23 Surface pressure cycle on Triton as simulated with the VTM assuming a different fixed N2 ice distribution in both hemispheres. A thermal inertia of 1000 J s−1/2 m−2 K−1 (SI) was assumed. Bertrand et al. (2022) include more simulations in their paper, and their Fig. 9 shows cases with different thermal inertias. It is of note that their simulations show that a lower thermal inertia would delay the peak of the surface pressure, in opposition with the occultation data points. The blue lines refer to a southern cap extending to the equator, while the pink lines are for a southern cap extending to 30° S. Each line, marked with its corresponding value, refers to a different extension of the northern cap: 45° N, 60° N, 75° N, and no cap.

In the text

Fig. 24 Stellar flux in Triton’s shadow for the 5 October 2017 event, normalised to unity outside the body (light blue region). The flux reaches a minimum of about 7% of the unocculted flux inside the shadow and then rises sharply at the shadow centre. The direction of Triton’s rotation is indicated, as well as the equatorial rotation velocity, vrot = 17 m s−1 (in an inertial frame), using the parameters of Table 2.

In the text

Fig. 25 Maps of the central flash intensity. Left panel: map of the central flash intensity, adopting the 1σ- upper limit ϵ′ = 0.0011 for the apparent oblateness of the flash layer (corresponding to a deprojected oblateness ϵ = 0.0019) (see text for details). The black labels along the iso-intensity contours (in curves) indicate the received stellar flux in units of its unocculted value. The grey diamond-shaped feature near the centre is the caustic curve described by Eq. (4) and corresponding to ϵ′ = 0.0011. In the case considered here, the flash layer is assumed to be aligned with the apparent direction of Triton’s poles, indicated by the dash-dotted line. Neptune’s direction is determined from the position angle 286° of Triton with respect to the planet at the moment of the occultation. Right panel: same, but adopting the oblate solution ϵ = 0.042 found by Elliot et al. (1997) from the shape of a central flash observed during the 14 August 1995 occultation.

In the text

Fig. B.1 Template temperature profiles and result of the ray-tracing code. Upper panel: Six template temperature profiles, indicated by the thin black lines, obtained taking values ϕ0 = 0.370, 0.365, 0.360, 0.355, 0.350, and 0.3445 (from left to right) of Triton’s contribution to the normalised occultation light curve at La Palma. The rightmost profile (corresponding to ϕ0 = 0.3445) is the profile obtained from the calibration images taken before the event. Exploring ϕ0 with small incremental steps, we find a best fit to the central flash for ϕ0 = 0.35885 (see the thicker black template profile, where the green dots labelled 1 to 4 are prescribed particular points in the Dias-Oliveira et al. 2015 model, as specified in Table B.1). We note that all the profiles go through the boundary condition T = 54.5 K at r = 1453 km (100 km altitude, upper green dot). The inverted profiles for La Palma corresponding to ϕ0 = 0.35885 are shown in red (ingress) and blue (egress). The horizontal dotted line marks the altitude of the central flash layer. Lower panel: Result of the ray-tracing code, using the best T(r) profile of the upper panel and the best fit value of p1400 = 1.18 µbar (Table 4). The normalised stellar flux is plotted against the distance to Triton’s shadow centre. The green curve represents the flux if only one stellar image is present, and the green dots show the correspondence with the points labelled 1 to 4 in the upper panel. The black curve is the sum of the fluxes from two stellar images, i.e. the sum of the green curve and its mirrored version with respect to the z = 0 axis. This black curve is then used to fit the observations.

In the text

Fig. B.2 Effect of the secondary stellar image on the inversion results. Upper panel: Effect on the temperature profiles. The smooth black line connects these profiles to the surface at 38 K, where the pressure is set to 14 µbar. The resulting density profile is used to generate synthetic light curves that feed the Abel inversion procedure. The solid green line is the retrieved T(r) profile with only the primary image accounted for. The green dotted line is the retrieved profile where the primary and secondary images have been added. The profiles obtained from the inversion of La Palma’s light curve at ingress and egress (in red and blue, respectively) are shown for comparison (see also Figs. 9– 11). Black horizontal solid and dotted lines show Triton’s surface and reference radius (1400 km). Lower panel: Same, but with the pressure profiles.

In the text

Fig. C.1 Data (black dots) fitted simultaneously with synthetic light curves (blue lines), based on the temperature profile displayed in Fig. 10 (black line) and the pressure boundary condition p1400 = 1.18 µbar (Table 4). The green dots are the residuals of the fits. We note that the central flash regions have been excluded from the fit so that we obtain a global fit that is not influenced by the deepest atmospheric layers. The stations with exposure times smaller than 1 s have been smoothed to have a sampling time close to 1 s, allowing a direct visual comparison of S/N of the various datasets. The lower and upper horizontal dotted lines mark the zero-flux level and the total star plus Triton unocculted flux, respectively. We note that the three central-most stations (Constância, Le Beausset, and Felsina observatories in Fig. C) are plotted at a different vertical scale to accommodate the presence of a strong central flash. Each panel has a duration of five minutes and is centred around the time of closest approach (or mid-occultation time) of the station to Triton’s shadow centre. The stations are sorted from left to right and top to bottom from the northernmost track (Newark) to the southernmost track (Athens; see Fig. C), projected on Triton in the sky plane (Fig. 5). For reference, the vertical red lines mark 23:48 UTC for the European and African stations, and 23:55 UTC for the US stations (Newark, Ithaca, and Dark Sky observatories).

In the text

	Fig. C.1 Continued. NB. ‘JAVA.’ is the abbreviation of Javalambre, used so that the name of the station fits into the plot.
In the text

	Fig. C.2 The same as Figs. C.1, but for stations that were not used in the simultaneous fit.
In the text