Subseasonal Variation in Neptune's Mid-infrared Emission

Combining unpublished observations with archival data, we collected all currently available imaging data that spatially resolved Neptune's disk in the mid-infrared (see Table 1). In total, this amounts to more than 95 images acquired between 2003 and 2020 at wavelengths ranging from 7.7 to 25.5 μm. Data were acquired using infrared instruments on 8–10 m class telescopes at multiple ground-based facilities, listed in Table 1. Images prior to 2011 were previously analyzed and published, in part, in Hammel et al. (2007), Orton et al. (2007, 2012), de Pater et al. (2014), and Sinclair et al. (2020) and further summarized by Fletcher et al. (2014), who investigated temporal variability between the Voyager 2 encounter (1989) and Neptune's southern summer solstice (2005). All subsequent imaging data are analyzed here for the first time, including Subaru-COMICS images from 2011, 2012, and 2020 and Very Large Telescope (VLT)-VISIR images from 2018. Observational details for all individual images are provided in the Appendix.

Collectively, images were acquired in 26 different bandpass filters ¹⁵ transmitting at wavelengths within the N band (7–13 μm) and Q band (17–25 μm), as summarized in Table 2 and Figures 1 and 2. Clustered at wavelengths of enhanced telluric transparency, these filters sense Neptune's emission spectrum in just a few distinct spectral regions, with sensitivity to different molecular transitions and pressures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For clarity, we broadly sort the majority of the images into the following four distinct groups defined by their passband and, correspondingly, the molecules and pressures they sense:

• 1.  
  17–25 μm sensing emission from molecular hydrogen (H₂), primarily at pressures near the tropopause (∼80–200 mbar; 15 images);
• 2.  
  11–13 μm, sensing emission primarily from stratospheric ethane (C₂H₆), with maximum contributions from ∼0.1 to 2 mbar (42 images);
• 3.  
  8–9 μm, sensing emission from stratospheric methane and its isotopologue, deuterated methane (CH₃D), primarily from ∼0.1 to 1 mbar (17 images);
• 4.  
  7–8 μm, sensing emission from stratospheric methane (CH₄) primarily from ∼0.01 to 1 mbar (15 images).

Additionally, there are six images with filtered central wavelengths of 10–11 μm that are weakly sensitive to both stratospheric hydrocarbons (mainly ethane and ethylene (C₂H₄)) and H₂ emission from deeper in the upper troposphere (at ∼1 bar). However, given their broad but weak sensitivity, low signal-to-noise ratio (S/N), and limited temporal coverage, these images provide few constraints and do not factor prominently into our analysis.

The radiances at mid-infrared wavelengths are fundamentally determined by the combination of the kinetic temperature of the atmosphere and the abundance of the emitting molecules at the pressures sensed. The abundance of hydrogen is assumed to be well mixed and known in the atmosphere, and thus emission from H₂ at 17–25 μm provides an unambiguous indication of the atmospheric temperatures in the upper troposphere. However, since the abundances of hydrocarbons are likely spatially and temporally variable given chemical sources, sinks, and transport (Moses et al. 2018, 2020), variation in the observed emission from methane (7–8 μm), deuterated methane (8–9 μm), and ethane (11–13 μm) alone cannot differentiate between variation in stratospheric temperature or composition. To help resolve this ambiguity, we include additional spectroscopic data sensitive to the stratospheric temperatures, as discussed in Section 2.2.

The atmospheric pressures sensed by each filter are estimated using radiative transfer modeling (NEMESIS; Irwin et al. 2008), assuming vertical temperature and chemical profiles based on Greathouse et al. (2011) and Moses et al. (2018). As illustrated in Figure 3, these contribution functions also group by filter passband, with some overlap among the defined groups. The precise shapes and locations of these contribution functions are dependent on the assumed atmospheric model and observational emission angle but broadly show sensitivity to gases at a range of pressures extending from the upper troposphere to the stratosphere. The hydrogen filters sense the deepest but also have additional minor contributions from a wide range of pressures. As noted by Fletcher et al. (2014), the contribution functions for the H₂ filters are somewhat bimodal, with contributions from lower pressures increasing with increasing emission angle. The ethane, deuterated methane, and methane filters sense progressively lower pressures in the stratosphere, with significant contribution from the strongest methane lines continuing at pressures less than 10⁻⁵ bars.

Zoom In Zoom Out Reset image size

The effective spatial resolutions of the images vary owing to differences in the telescope apertures, focal lengths, detector plate scales, filtered wavelengths, and atmospheric seeing. The mean spatial resolution amounts to Neptune's disk subtending ∼28 pixels at 44 of latitude per pixel, although the highest-resolution images—acquired in 2018 with VLT-VISIR at a plate scale of 00453 pixel⁻¹—offer twice this resolution. Subaru-COMICS images have the poorest resolution at only 013 pixel⁻¹, subtending closer to 7° of latitude per pixel at the disk center. The effective seeing disks varied between approximately 03 and 075, with an average of 047 based on the FWHM of the calibration stars (see Figure 4). This limited resolution has the effect of blurring the edges of the disk with the sky, artificially suppressing the emission for the outer 10%–30% of the disk. This must be kept in mind when comparing and interpreting observations.

Zoom In Zoom Out Reset image size

Over the 17 yr covered by ground-based thermal imaging, the subobserver latitude changed by less than 44, from 291S to 247S (see Figure 4). This small change in subobserver angle amounts to a change of 2 pixels or fewer across all images, which implies that most observed changes cannot be attributed to changes in the observing geometry. Changes in emission angle can potentially modify the perceived radiances along the extreme northern and southern edges of the disk, but this is unlikely to be significant over the span of our data (although it will become significant as the south pole appears closer to the limb in the coming years). Likewise, changes in the subsolar latitude are modest, with a minimum value of 292S in 2005 images (near southern summer solstice, L_s = 270°) and maximum of 241 S in 2020 images (L_s ≈ 303°).

As standard for mid-infrared observations, images were acquired using the chopping-and-nodding technique to effectively remove the contribution of the sky and telescope's thermal emission from the image. By this process, hundreds of ∼100 ms exposures—short enough to preclude detector saturation—were combined to build up a sufficient S/N for the target. Total on-target integration times ranged from several minutes to roughly 30 minutes, depending on the filter and observatory. In most but not all cases, accompanying observations of standard stars were acquired for flux calibration.

Radiances for each pixel were calibrated using radiance conversion factors derived from observations of standard stars. These conversion factors are derived by comparing the observed stellar signal in digital units to the expected flux densities for standard stars. Stellar flux densities are often provided by the observatories, but we calculated stellar flux densities from Cohen stellar spectra (e.g., Cohen et al. 1995) for the filter transmissions assuming typical atmospheric conditions for the observatories (Noll et al. 2012). These flux densities were in very close agreement—within 3% or less—of values currently maintained by Gemini ¹⁶ and the European Southern Observatory. ¹⁷ The standard stars used are listed along with each Neptune observation in the Appendix. To account for differences in air mass between Neptune and the calibration star, air-mass extinction corrections were calculated from observations of stars repeated at different air masses, when available, or otherwise estimated from values of Krisciunas et al. (1987). Extinction corrections typically amounted to less than 5%. In the few cases when stellar observations were absent, calibrations were estimated from observations on proximate nights.

Resulting uncertainties in the image radiances are dominated by uncertainties in the radiance conversion factors, which generally vary by star, wavelength, and time. Statistical analysis of mid-infrared conversion factors at VLT by Dobrzycka & Vanzi (2008) has shown typical temporal variability in these conversion factors of up to 10% in filters at 8.6 and 11.2 μm, where atmospheric transparency is greatest, and less than 20% for the Q2 filter, where telluric water vapor decreases transparency. With minutes to hours typically lapsing between Neptune and stellar observations, errors can become significant, particularly at wavelengths at which the emission from the sky is greatest (see Figure 2 and atmospheric transmissions in Table 2).

Similar wavelength dependence in the uncertainties can be inferred from variability in repeated observations of Neptune on relatively short timescales. For example, six VISIR observations acquired at 12.2 μm (NeII_1) on five nights in 2009 August and September varied by less than ±10%, with a standard deviation of only 6%. During the same period, five images at 8.9 μm varied by less than 14% with a standard deviation of 10%, while the methane images at 7.8 μm—where the atmospheric transparency is less than 50%—varied by ±30% of their mean, with a standard deviation of about 25%.

Differences in the spatial resolution and filter passbands can also lead to relative differences in the inferred trends across the disk, even at similar wavelengths. For example, Neptune was observed by both VISIR (12.2 μm, NeII_1) and COMICS (12.4 μm, F09C12.50W1.15) on the same night in 2008 September, but the VISIR image shows significantly greater limb brightening (see Figure 26 in the Appendix). Given the contemporaneity, the difference is observational, likely owing to the relatively finer plate scale resolution, lower air mass, and narrower passband (centered on the ethane emission) of the VISIR observation. As such, even images at very similar wavelengths should be compared with some caution.

Aside from the conversion factors, repeated attempts to calibrate data with plausible alternative choices for parameters regarding the atmospheric model, sky subtraction, area of aperture photometry, and air-mass correction resulted in changes in radiance of less than 5%. As a whole, the data do not exhibit any obvious trends in relative radiances versus air mass or precipitable water vapor. Altogether, we estimate uncertainties of up to 10% in radiance for ethane images (11–13 μm), 20% for the CH₃D images (8–9 μm), 25% for methane images (7–8 μm), and, conservatively, 30% for hydrogen images (17–25 μm).

Correcting for blurring.—Given limited spatial resolution and Neptune's ∼23 angular diameter, a simple direct inversion of the images will provide unreliable temperatures nearer the edges of the disk, where blurring with space surrounding the disk artificially reduces the observed radiances. In theory, a deconvolution of the images with the image point-spread function (PSF; estimated from the stellar seeing disk) in Fourier space can correct for this degradation; however, direct deconvolutions amplify the noise, which we find can become dominant in most images, requiring yet more smoothing to overcome. We instead adopted an alternative approach to effectively deconvolve the images by modeling idealized data and the effects of blurring, following Roman et al. (2020).

We began by extracting radiances from the central meridian of the disk in the unaltered images, but limited to points on the disk with modest emission angles (i.e., μ ≥ 0.5). Selected ethane (11–13 μm), hydrogen (17–25 μm), and methane (7–9 μm) images were used, with multiple images combined to give averages for each filter group representative of each observing epoch. These extracted radiances were then inverted using an optimal estimation retrieval algorithm, NEMESIS (Irwin et al. 2008), to create a rough model of the atmospheric temperature structure versus pressure and latitude. A vertical temperature profile based on Moses et al. (2005) and Greathouse et al. (2011), with chemical abundances of Moses et al. (2018), was used as the prior (as shown in Figure 3). Inferred temperatures at modest emission angles were extrapolated to latitudes viewed at higher emission angles (i.e., μ < 0.5) to complete the disk. Assuming zonal uniformity, we then forward-modeled radiances for all locations on the disk from this initial model of temperatures (as a function of pressure and latitude), using NEMESIS to solve the radiative transfer at the relevant observing geometries for each pixel. Forward models were calculated at the resolution of the spatially normalized data (i.e., for a disk with an equatorial width of 51.8 pixels), with care taken to adequately represent the contributions from the limb. The resulting modeled radiances yielded an initial simulated image of the idealized disk as a first approximation.

We then used this simulated image to estimate the effects of blurring in the data. The simulated image was convolved with its corresponding data's calibration star image (which we assume approximates the effective PSF of the data) to mimic the blurring suffered by the real disk of Neptune in the observations. We then directly compared the idealized and blurred disks to define a flux-conserving multiplicative correction (for each pixel) that converts between the two. These correction factors, determined from the simulated images, therefore represent the transformation of Neptune's disk due to atmospheric blurring. We applied this factor inversely to the real data, with the goal of approximating how the true disk would appear prior to blurring. Then, in order to evaluate whether this corrected disk was truly consistent with the data, the corrected image was artificially blurred (via the same stellar convolution) for direct comparison with the actual observations. Fractional differences between observed and artificially blurred disks were applied to the correction factors, and the process was repeated, iteratively adjusting the factors as needed until the data and model agreed to within 10% or less (see Figure 18). The process is imperfect, as it can introduce subtle spurious structure in the idealized model, but this fine structure can be removed by carefully smoothing the solutions prior to retrieval while preserving the larger latitudinal trend.

Zoom In Zoom Out Reset image size

Retrievals from corrected images.—The resulting idealized images—now effectively our best approximation of the data prior to being blurred by the atmosphere—were then used as the inputs into the retrieval algorithm. We extracted strips of radiances along the central meridians spanning all visible latitudes (averaged over 9 pixels in longitude, amounting to ∼04, or 20° of longitude at the equator), but now extended to emission angles of ∼81° (μ = 0.15). These radiances were inverted using the NEMESIS code (Irwin et al. 2008), again with the same a priori temperature and chemical profiles and vertical correlation length scale of 1.5 pressure scale heights. For clearer interpretation of results, we allowed only temperature to vary in our retrievals, as the dearth of contemporaneous quadrupole measurements makes it impossible to constrain both temperature and chemical abundances independently across all the years. While both temperature and composition likely vary in Neptune's atmosphere, the analysis of Section 3.4 suggests that it is reasonable to attribute the observed changes primarily to variation in temperatures. The seasonal photochemical models of Moses et al. (2018) do predict slight changes in ethane mixing ratios over the observed period (<19% at 0.5 mbar between L_s ∼ 265° and L_s ∼ 304°), but we use the 2009 chemical model (L_s ∼ 278°, roughly corresponding to the average date of our data) for all retrievals and forward models for consistency. We tested the sensitivity of our results to this choice, and we found that retrieved 0.1 mbar temperatures differed by less than 2 K when neglecting time variability in the chemical model over this period, with maximum differences at the south pole. A more thorough study of simultaneous temperature and compositional gradients across the disk using contemporaneous quadrupole and ethane emission (e.g., following Greathouse et al. 2011) is left to future work.

Temperature profiles were modeled from ∼10 to 1 × 10⁻⁸ bars of pressure, over 180 vertical layers. However, the limited sampling and poor vertical resolution of the filter contribution functions mean that only pressures near the tropopause (10⁻¹ bars) and lower-stratosphere (10⁻³–10⁻⁵) bars were strongly constrained; temperatures at other pressures tended toward the a priori profile, and the implied lapse rates are only approximate. Clouds and aerosols were neglected, given their presumed negligible opacity at these infrared wavelengths. Absorption coefficient (k) distributions from line-by-line calculations were used to model the gaseous opacity, calculated from the GEISA 2003 database (Jacquinet-Husson et al. 2005). Collision-induced opacities were taken from Fletcher et al. (2018b) for H₂–H₂ and Orton et al. (2007) and Borysow et al. (1988, 1989) for H₂–He and H₂–CH₄. Images were sorted into six epochs for retrieval—those dating from 2003, 2005, 2006, 2008/2009, 2018, and 2020. All retrievals included selected ethane (11–13 μm) and hydrogen (17–25 μm) images, with multiple images averaged over each epoch to improve S/N. In 2005, since no hydrogen images were available, images from 2006 were substituted. Likewise, stratospheric images from 2008 and 2009 were combined, owing to their limited number and similar radiances, and paired with the hydrogen images from 2008 (no hydrogen images were acquired in 2009). Given the questionable photometric calibrations of the methane images, we perform two sets of retrievals—those with and without representative methane or CH₃D (7–9 μm) images included—to isolate their effect on the inferred temperature structure.

We assume local thermodynamic equilibrium (LTE) in our retrievals, neglecting the effect of non-LTE on the emission. Non-LTE emission becomes increasingly significant as the pressure and wavelength decrease. Appleby (1990) estimated that at 7.7 μm, neglecting non-LTE can result in errors of up to 2 K at 0.1 mbar, increasing to 20 K at 0.1 μbar. While the 12 μm ethane images have a peak contribution around 0.5 mbar, the ∼8 μm methane images do have peak contributions near 0.1 mbar, with lesser contribution from as little as 1 μbar, particularly at low emission angles near the limb. This leads to an additional error of roughly 2 K in temperatures retrieved from ∼8 μm data, growing slightly at extremely high and low latitudes.

Finally, the retrieved meridional temperature structures were used in a forward model, assuming zonal uniformity, to once again produce synthetic images to validate by comparison with the data. The resulting idealized images were convolved with corresponding stellar PSFs and degraded with synthetic noise to allow a direct comparison with the observations. Random noise was modeled as a normal distribution of randomly generated values with a mean of zero and a standard deviation equivalent to that measured in the corresponding images (off of the disk).

As is typically the case with such retrievals, our derived solutions are nonunique, especially considering that potential compositional changes will alter the retrieved temperatures. For example, if methane was relatively depleted in the stratosphere at high latitudes—as it is thought to be in the troposphere (Karkoschka & Tomasko 2011; Irwin et al. 2019)—we would be underestimating the polar temperature in retrievals from the same radiance. Nonetheless, given the limits of our data, we attempted to represent the simplest solutions to the temperature field consistent with the observations. We assessed the accuracy of the solutions by evaluating the reduced χ² of the retrievals and percentage differences between the modeled and true images. For the reduced χ², we assume

$$\begin{eqnarray}&&{\chi }_{\nu }^{2}=\sum _{i=0}^{\nu }\displaystyle \frac{{\left({\mathrm{data}}_{i}-{\mathrm{model}}_{i}\right)}^{2}}{{\sigma }_{i}^{2}}\displaystyle \frac{1}{\nu },\end{eqnarray} \tag{ 1 }$$

where we sum over the number of observed radiances ν, and σ is taken to be the calibration uncertainty.

Comparisons between data and synthetic observations constructed from retrievals are shown in Figure 19. In the figure, we refer to the models constructed from retrievals using only hydrogen and ethane images as model-A and those that additionally include methane images—or, if methane images are absent, CH₃D images—as model-B. The corresponding temperature structures inferred from the retrievals are shown in Figures 20 and 21, along with reduced χ² values of each retrieval versus latitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The goodness of the fits and the temperature structures implied by the images vary significantly by epoch. The agreement between data and model is strongest in the hydrogen-sensing filters, where disk-averaged differences are just a few percent or less in nearly all cases. Models of the stratospheric-sensing images have greater errors, with values depending on whether methane images are included in the retrieval. When methane and CH₃D images are ignored (i.e., model-A), the ethane images are reproduced to within an average error of up to 10%, with the largest errors seen at the northern and southern limbs. The χ² values are generally very small, suggesting that our assumed uncertainties are too large, but χ² values increase significantly at latitudes south of −70° and north of 20° (see Figure 20). These errors reflect the challenges of modeling the radiances at high emission angles, even for a single stratospheric filter, and likely indicate residual errors in our attempts to effectively deconvolve the images, as well as possibly errors due to our plane-parallel approximation of the radiative transport. The inferred temperatures at these poorer-fit locations should be considered somewhat less accurate, but they nonetheless yield a reasonable match between the data and modeled images.

The contour plots of Figure 20 show retrieved temperatures for model-A. Accompanying contribution functions (to the right of each panel) indicate the pressures at which temperatures are actually constrained, and Figure 22 shows the retrieved temperatures versus latitude near the contribution peaks. In the upper troposphere, observations are consistent with the same basic temperature structure inferred from Voyager-IRIS. The tropopause temperatures in our data reach as low as 49–52 K at midlatitudes and rise roughly 4–8 K warmer at the equator and south polar regions (Figure 22(c)). In comparison, the Voyager-IRIS 100 mbar temperatures ranged from ∼51 K at 45°S to ∼57 K at the equator (Conrath et al. 1998; Fletcher et al. 2014). In our measurements, systematic uncertainties in radiances translate to a roughly 3 K systematic uncertainty in temperatures at the ∼100 mbar tropopause, with random noise and retrieval error adding uncertainties of ∼2 K. The precise temperatures retrieved are also dependent on the wavelength of the Q-band filter used, since different filters sense slightly different pressures, as indicated by the contribution functions to the right of each panel in Figure 20. Differences of less than 2 K between epochs are likely not significant.

Zoom In Zoom Out Reset image size

In the stratosphere, the temperature structures are more variable. When the methane/CH₃D emission is ignored, we see a pattern with cooler temperatures at equatorial and southern low latitudes compared to warmer temperatures at the south pole and, in most cases, northern midlatitudes. Accompanying plots in Figure 22(a) show the range in retrieved temperatures versus latitudes at just the 0.5 mbar peak of the ethane contribution function. Latitudinal temperature contrasts between the warmer south pole and cooler equator increased from roughly 8 to 28 K between 2003 and 2020; polar temperatures rose from 152 to 163 K between 2018 and 2020 alone. Uncertainties are roughly ±2 K, increasing to 4 K near the poles and edges of the disk, where fits are poorer. Models in 2008/2009, 2018, and 2020 best match the data, while 2003 and 2005 are poorest but still good at central latitudes.

The retrieved stratospheric temperatures of Figure 20 were constrained by 11–13 μm radiances; they did not consider the radiances from methane or CH₃D images, and their modeled emission poorly reproduced the observed radiances at 7–9 μm. At these shorter wavelengths, the forward-modeled emission is either too dim (2003, 2005, 2006, 2020) or too bright (2008/2009, 2018), with average errors up to 30%. This discrepancy suggests that the simple temperature and chemical structures inferred from the ethane images and the a priori alone are not complex enough in all years.

When the 7–9 μm radiances were additionally included in the retrievals (model-B), the retrievals responded by increasing the 0.01 mbar temperatures in 2003, 2005, 2006, and 2020, particularly at central latitudes, while slightly cooling the same heights in 2008/2009 and 2018 (see Figures 21 and 22(b)). The modeled methane/CH₃D images are improved, with errors reduced to 15% or less, but the ethane fits are in most cases worsened from errors of 10% to 15%, as the retrievals apparently struggle to fit both channels simultaneously. In 2003, in particular, the ethane model appears too bright, while the methane model is still too dim. The corresponding reduced χ² values greatly increase compared to those without methane or CH₃D, with values approaching or exceeding 10 in 2003 and 2006. In contrast, both filter groups can be simultaneously fit to within a few percent in 2018, with reduced χ² values generally remaining ≤1. The reduced χ² values still increase toward the edges of the disk but now also show an increase near the equator in most cases, corresponding to warmer equatorial temperatures. The brightness of the south polar feature also tends to be slightly too dim in our modeled ethane images, while simultaneously too bright in our modeled methane images. This inability to fit both simultaneously may indicate that additional compositional variation may be present, potentially owing to the interplay between photochemistry and dynamical transport. For example, polar downwelling could potentially increase in the mid- and lower-stratospheric mixing ratio of ethane while simultaneously reducing that of methane, given the different vertical gradients in the mixing ratios of each (i.e., ethane increases with height, methane does not; Moses et al. 2005, 2018). Alternatively, it may also simply reveal errors in the calibration or the consequence of neglecting non-LTE emission, which makes it impossible to simultaneously reproduce radiances in multiple filters with similar contribution functions.

The larger calibration uncertainties and lower S/N associated with the 7–9 μm imaging mean that retrieved temperatures are only certain to within ±3 K, but errors may increase to ±8 K near the pole and edges of the disk.

Previously published values of latitudinal temperature and/or chemical abundance gradients in the stratosphere are limited and appear to differ from our results. Voyager mid-IR measurements of C₂H₂ emission in 1989 (L_s ∼ 236°) showed a maximum range of 20 K (or a factor of two in C₂H₂ abundance) at the 0.03–2 mbar level, with minimum at 55°S and maximum between 0° and 20°S (Bézard et al. 1991). We find similarly large contrasts at 0.5 mbar in the 2018 and 2020 data (L_s ∼ 299° and L_s ∼ 303°; Figure 22), but with a minimum at 20°S and a maximum at the pole. A minimum at midlatitudes is found in 2003 (L_s ∼ 266°) and 2008/2009 (L_s ∼ 278°), but with only modest temperature contrast (3.5–5 K) between midlatitudes and the equator. Given that 15 yr or more separate these observations, it is possible that these considerable differences may be explained by genuine temporal variability and/or differences between the ethane and acetylene emission.

In contrast, two subsequent studies of Neptune spectra contemporary to the imaging data show comparatively little variation with latitude. Fletcher et al. (2014) analyzed Keck-LWS spectra from 2003—companion spectroscopy to the imaging analyzed here—and determined that temperatures only varied by 3 K (155–158 ± 2 K) between 20°N and 80°S at 0.5 mbar. Likewise, Greathouse et al. (2011) analyzed multiple 2007 Gemini-TEXES spectra, including methane and ethane emission in addition to the quadrupole emission analyzed here, to determine temperatures at 2.1, 0.12, and 0.007 mbar. To within measurement uncertainties, their retrieved values were consistent with latitudinally uniform temperatures, with variation of no more than 2 K between ∼70°S and 12°N (153.5–155.5 ± 2 K at 0.12 mbar, 122–124 ± 2 K at 2.1 mbar). The variation increased to 5 K if the assumed stratospheric methane volume mixing ratio (VMR) was doubled to 1.5 × 10⁻³. Note that our assumed CH₄ profile from Moses et al. (2018) has a mixing ratio of 1.15 × 10⁻³ at 5 mbar, as derived from the Herschel observations of Lellouch et al. (2015), decreasing to 0.9 × 10⁻³ at 0.5 mbar as a result of vertical diffusion and chemical loss near the lower-pressure homopause. In both studies, the marginally greater temperatures appeared nearer the equator, with a slight minimum at southern midlatitudes (i.e., 55°S in the TEXES analysis). These published values of 2003 and 2007 temperatures do appear relatively uniform in latitude compared to most of our inferred stratospheric temperature gradients at 0.5 mbar, shown in Figure 22(a), but note that our 2003 and 2008/2009 curves have the weakest gradients. The 2003 curve only varies by 3.5 K between the equator and 60°S (144.5–148 ± 2.5 K), while the 2008/2009 curves vary similarly between 20°N and 70°S (141–144.5 ± 2.5 K). Both curves also show a minimum at 40°–50°S and an increase toward the equator over this latitude range. Over this limited range, our results are largely consistent in trends with the previous analyses of separate but nearly contemporaneous data. It is only at latitudes observed nearer the edges of the disk that we detect significantly larger gradients in these same years. This discrepancy is perhaps reasonable to expect, given that the spectra were acquired at lower spatial resolution and with no attempted correction for beam dilution near the edges. As suggested by the comparison between quadrupole spectra and ethane images (Figure 16), beam dilution can suppress the radiances and consequent temperature gradients across the disk.

We also note that our retrieved 2003 temperatures (neglecting methane emission) are also about 10 ± 3 K colder at 0.5 mbar than Fletcher et al. (2014) retrieved from the contemporaneous spectra. Yet the 0.1 mbar disk-integrated temperatures between 2003 and 2008 (as plotted in the inset of Figure 23) are remarkably consistent with contemporaneous 0.1 mbar temperatures reported by Fletcher et al. (2014)—well within the ∼±3 K uncertainties for retrievals from ethane images. The apparent discrepancy at 0.5 mbar may therefore be explained by differences in the vertical resolution of the two data sets. The relatively broad contribution functions of the imaging data create uncertainty in the retrieved temperatures at precise pressures, particularly when the vertical temperature gradient is large.

Zoom In Zoom Out Reset image size

To better assess our vertical temperature gradients, we computed disk-averaged temperature profiles for each epoch, as shown in Figure 23. Compared to the profile of Fletcher et al. (2014), derived from 2007 AKARI-IRC spectra (Fletcher et al. 2010), our temperature profiles are slightly colder and increase more slowly with height at 10–1 mbar, but more rapidly with height at <0.5 mbar, before becoming nearly isothermal at relatively lower pressures. As a result, our retrieved temperatures are comparatively cooler at 0.5 mbar (consistent with the discrepancy noted above) but warmer at 10⁻⁵ bars, such that the integrated thermal emission is comparable. This difference in lower-stratospheric lapse rates appears to follow directly from equivalent differences in the chosen a priori profiles, indicating that the precise lapse rate at these pressures is not strongly constrained by these imaging data. Hence, considering the discrepancy, we may ascribe an uncertainty of roughly 10 K when defining temperatures at precise pressures in this region of strong vertical gradient. This weakness is unsurprising, given the relatively broad contribution functions of the imaging filters. The greater vertical resolution provided by spectroscopic data is better suited for vertical structure analysis, and so the profiles and precise pressures reported here should not be interpreted as contradictory to the previous work. Rather, we emphasize that the primary strength of the images is their ability to clearly reveal variation in radiance with latitude and time, and this variation will significantly affect the disk-averaged temperature retrievals.

We find that our retrieved temperature profiles show variation across the years at pressures less than 1 mbar. If the methane radiances are discounted, the profiles mostly fall between the profiles of Moses et al. (2005) (based on spectral observations from the early 1990s) and Greathouse et al. (2011) (based on 2007 TEXES data), from which studies our prior is derived. Our profiles differ by a maximum of only about 6 ± 4 K at 1 mbar and 8 ± 4 K at 0.1 mbar, from a maximum in 2003 to a minimum in 2018. If the methane radiances are included, the range increases to 20 ± 14 K at 0.1 mbar, owing mostly to much warmer temperatures in 2003–2006 and cooler temperatures in 2008/2009. These results are roughly consistent with differences of less than 10 K at 1 mbar and 6 K at 0.1 mbar among the multiple retrievals from ground-based spectroscopy (2003–2007) presented in Fletcher et al. (2014). However, our 2003 result based on methane radiances appears as a questionable outlier—15^{+ 10} ₋₈ K warmer at 0.04 mbar than was retrieved from the contemporaneous Keck-LWS spectra. This discrepancy possibly indicates a calibration error in the methane radiances, which may be exaggerating measured temperatures in a year that already appears exceptionally warm based on ethane images alone.

Seasonal variation due to changes in the subsolar latitude are expected to produce variability in Neptune's temperatures (e.g., Conrath et al. 1990; Greathouse et al. 2008) and stratospheric photochemistry (Moses et al. 2018), both of which can potentially alter observed mid-IR radiances.

In radiative modeling, seasonal temperature variation is typically characterized by the length of the radiative time constant relative to the orbital period. Radiative time constants near the tropopause are relatively long (approaching a century; Figure 3), and so the observed lack of variation in Q-band images is consistent with expectations. However, given shorter radiative time constants in the stratosphere (Li et al. 2018), the stratospheric temperatures should vary across seasons, reaching a maximum in the summer hemisphere shortly after the peak insolation at summer solstice (Conrath et al. 1990; Greathouse et al. 2008). The imaging data bracket Neptune's 2005 southern summer solstice (L_s ∼ 266–303), but rather than warming, we observe a decline in the disk-integrated stratospheric radiances and implied temperatures. Only between 2018 and 2020 do the polar regions appear to increase in radiance.

Likewise, the photochemical modeling of Moses et al. (2018) predicts a peak in the southern hemisphere ethane abundance following the maximum solar flux at solstice. As methane is photolyzed by high-energy photons, its photolysis leads to the production of ethane and other hydrocarbons. Changes in photochemistry are expected to be greatest at lower pressures, where CH₄ photolysis rates are largest and chemical time constants are shortest, and higher latitudes, where insolation varies most greatly with seasons. As with the temperatures, the variation is muted and delayed with increasing pressure—in this case due to longer advective timescales as the ethane slowly descends. The Moses et al. (2018) model predicts that ethane should vary significantly at pressures less than 1 mbar, varying annually by a factor of nearly ∼2.5 near the poles at 0.1 mbar. Across the southern hemisphere, the 0.1 mbar ethane VMR and column abundance are expected to peak around L_s ∼ 300 (roughly 2019) and fall gradually to a minimum following winter solstice. In contrast, the seasonal variation of methane should appear insignificant in comparison to its overall abundance at the pressures probed by the 7–8 μm observations. Therefore, variation in seasonal photochemistry should cause ethane radiances (11–12 μm) to increase near solstice, while the methane radiances (7–9 μm) remain constant. Yet this appears inconsistent with our observations, which decline in radiance across all stratospheric-sensing wavelengths between 2003 and 2010.

The explanation for the apparent inconsistency between seasonal models and observations may be resolved by the feedback between photochemical and radiative processes. While methane absorbs sunlight and warms the atmosphere, photochemically produced hydrocarbons—primarily ethane and acetylene—are powerful infrared emitters that serve to cool the stratosphere. The balance between this radiative heating and cooling changes as the amount of photochemical hydrocarbons changes. Applying a coupled radiative-chemical model to Saturn's atmosphere, Hue et al. (2016) showed that this interplay can produce a peak in summer temperatures prior to maximum insolation at summer solstice, as the late spring production of infrared emitters overwhelmingly cools the stratosphere, counteracting the growing solar heating. Although the effect was limited to pressures less than 0.1 mbar in the Saturn study, an equivalent process in Neptune's atmosphere may help explain the observed decline in stratospheric radiances beginning in 2003—prior to the solstice. Since the observed radiances are sensitive to both the temperature and abundance, the indirect effect of cooling the atmosphere would need to dominate over the direct effect of increasing the ethane to explain the decreasing trend at 11–13 μm radiances. Indeed, our modeling in Figure 14 shows that expected seasonal variation in ethane would produce only a very slight increase in the disk radiances (≲1 × 10⁻⁸ W cm⁻² sr⁻¹ μm⁻¹) over the observed period. The likely important effect of photochemistry on the seasonal temperatures and radiances suggested by our results should be further investigated with coupled radiative-chemical modeling in future work.

Nonetheless, given Neptune's 165 yr orbital period, any seasonal changes are expected to occur gradually over decades. The rapid changes observed between 2018 and 2020 appear surprisingly swift for seasonal response, particularly considering that the south pole has been constantly illuminated since 1963. Nor can the slowly changing subsolar latitude explain the sudden change in meridional gradients between 2006 and 2008. Additional processes appear to be operating in Neptune's atmosphere on subseasonal timescales, and on both regional and global scales.

Cloud and haze activity.—As evidenced by sporadic clouds and occasional vortices, Neptune clearly exhibits variable weather in the upper troposphere and lower stratosphere on subseasonal timescales; therefore, it is worth considering whether this weather can be related to the stratospheric variability.

In an analysis of Hubble Space Telescope (HST) visible and NIR imaging, Karkoschka (2011) found that trends in discrete cloud activity changed on a timescale of ∼5 yr. These discrete clouds were inferred to form near the tropopause. It was speculated that discrete cloud variation may be related to variation in large-scale dynamics at these pressures. Karkoschka (2011) also reported variation in the albedo of hazes and dark bands on similar timescales, but these were interpreted as changes in the aerosol abundances located much deeper in the atmosphere (>1 bar).

The observed maximum mid-infrared radiances in 2003 roughly coincided with the peak in the discrete cloud coverage in 2002 as determined by Karkoschka (2011). These peaks also roughly corresponded with the apparent seasonal peak in Neptune's disk-integrated visual albedos (472 and 551 nm), which had steadily brightened since the 1950s before leveling off in the early 2000s (Lockwood & Thompson 2002; Sromovsky et al. 2003; Lockwood & Jerzykiewicz 2006; Lockwood 2019). The subsequent decline in the stratospheric emission across multiple filters from 2003 to 2009 then coincided with a period of declining discrete cloud and plateauing visual magnitudes.

These coinciding trends suggest an intriguing connection between the mid-IR, discrete cloud cover, and the visual albedo, but the physical mechanism linking these variables is not obvious. Most easily observed at near-IR wavelengths, the high discrete clouds appear to contribute relatively little to overall visual (∼500 nm) albedo, which is thought to be most sensitive to deeper aerosols (Karkoschka 2011). Correlated variability would therefore suggest a link over many pressure scale heights, linking the upper troposphere and stratosphere. Implied variation near the tropopause is important because it marks the broad, cold boundary across which tropospheric methane must be transported to the stratosphere to produce ethane and other photochemical derivatives. Precisely how methane is transported from the troposphere to the stratosphere is unresolved given the tropopause saturation cold trap (Smith et al. 1989; Baines & Smith 1990; Orton et al. 2007; de Pater et al. 2014), but moist convection has been suggested to play a role (Stoker 1986; Lunine & Hunten 1989; Sinclair et al. 2020). If abundant clouds and hazes are indicative of convection or higher methane humidity, then perhaps some meteorological upwelling mechanism temporarily enriched the tropopause and lower stratosphere in hydrocarbons prior to 2003. As a consequence, the radiative heating rates at these levels could have temporarily increased, until the upwelling ceased, leaving the methane, clouds, and temperatures to decline over the following years. However, this is purely speculative, as no variation in the stratospheric methane abundances has been identified.

Alternatively, inertia–gravity waves originating from intermittent convective plumes, vortices, or other long-lived, cloud-generating disturbances may cause variable heating as waves break in the stratosphere. Indeed, Roques et al. (1994) showed that variable gravity wave heating could explain anomalous, variable structures in the upper-stratospheric temperature profiles inferred from stellar occultations between 1983 and 1990. Periodic, subseasonal variation in this heating may be possible if it is associated with planetary oscillations, such as those seen on Jupiter and Saturn, which perturb stratospheric temperatures at low latitudes (Leovy et al. 1991; Friedson 1999; Orton et al. 2008; Fletcher et al. 2017; Guerlet et al. 2018). If analogous oscillations are operating on Neptune, they may help explain the mid-IR variation at lower latitudes, but many decades of consistent observations would be needed to establish a periodicity with statistical significance. The implied possible correlation between clouds and temperature so far is still based on a limited number of observations and not without exception, as an unusually cloudy outburst at the equator in 2017 (Molter et al. 2019) has shown.

Other recent visible changes in clouds and hazes at greater pressures include a dark vortex in the north tropics, first observed by HST beginning in 2018 (Simon et al. 2019). Dark vortices are examples of occasional disturbances that can potentially generate waves and affect weather over a large vertical range. They appear to be localized regions of reduced scattering or upwelling of low-albedo aerosols from below, often triggering orographic-like companion clouds at higher altitudes (e.g., Smith et al. 1989; Sromovsky et al. 1993; de Pater et al. 2014; Hueso et al. 2017). The precise height and vertical extent of Neptune's dark vortices are unknown, but retrievals suggest that they are centered well below the stratosphere, nearer a pressure of 7 bars or more (Irwin et al.2022).

Dynamical modeling of the dark vortex predicts that these features are associated with a decrease in the methane mixing ratio within the vortex, along with a slight drop in temperature above (Stratman et al. 2001; Hadland et al. 2020), although Voyager detected very little if any temperature anomaly over the Great Dark Spot (Conrath et al. 1989) in 1989. The recent prominent dark vortices were imaged by HST at visible wavelengths multiple times, including in 2020 August, just 3 weeks after our mid-infrared Subaru-COMICS observations (Figure 25). Our imaging data show no obvious localized anomaly in the upper troposphere or lower stratosphere in 2020, but the low S/N and blurring of long integration times render such detection unlikely. However, we do note that the equatorial stratosphere appears anomalously cold over most of the disk in recent years, with the exception of the south pole. Aside from the Voyager-IRIS spectra in 1989, no spatially resolved mid-IR data have been available during any past vortex appearances (see, e.g., Wong et al. 2018), and so it is unclear whether the relatively large stratospheric temperature gradient between the cool equator and warmer, northern midlatitudes is in any way related to the dynamical environment in which these vortices form. We also note that the discrete cloud features—so vividly associated with the dark vortex in Voyager images—appear largely absent from this present spot in 2020 (Figure 25) and 2021, despite earlier companion clouds seen in 2018 and 2019. Although a number of dynamical factors may explain the apparent lack of clouds (e.g., the vertical extent of the vortex), the limited condensation, combined with our temperature findings, may suggest that the overlying atmosphere has grown not only colder but also drier in recent years.

Zoom In Zoom Out Reset image size

Establishing whether any of the meteorological changes discussed above are related to the trends in the stratosphere will require additional observations at mid-IR, near-IR, and visible wavelengths. Potential processes coupling the troposphere and stratosphere should continue to be investigated in future work.

Aside from seasonal modulation, the solar flux varies as part of the roughly 11 yr solar cycle. While the total solar irradiance differs by little more than 0.1% over a typical solar cycle (Kopp 2019), ¹⁹ variation in regions of the UV—as represented by a time series of 121.57 nm Lyα irradiance measured at Earth in Figure 24—exceeds 40%. High-energy photons (of wavelengths less than ∼145 nm) are the main drivers behind the methane photochemistry (Moses et al. 2020), and so modulation in the UV flux can potentially produce observable variation in photochemistry if chemical timescales are sufficiently short (Moses & Greathouse 2005).

A strong solar maximum occurred in ∼2001 (cycle 23), followed by a drop in the Lyα over the subsequent 10 yr. The deep minimum in 2009 (cycle 24) gave way to a weak solar maximum in ∼2014. Currently, solar cycle 25 is rebounding from another deep minimum in 2019. Plotted with the time series of Lyα, Neptune's mid-IR radiances followed a roughly similar trend.

However, as with the seasonal cycle, most of the photochemical response to changing far-UV flux of the solar cycle is expected to occur at relatively lower stratospheric pressures than the mid-IR images sense—nearer to 1 μbar, where H-Lyα and other short-wavelength radiation is absorbed and the radiative, chemical, and transport timescales are shortest. At the 0.05–1 mbar pressures sensed by the ethane filters, the photochemical variation will be dampened (Moses et al. 2018). And although methane filters are sensitive to somewhat lower pressures (0.005–0.5 mbar), the methane mole fraction is not expected to show appreciable change given that its photochemical loss is easily replenished from the abundant source below. Combined with the quadrupole observations, the observed mutual variation in the radiances would therefore again suggest a change in the temperatures, rather than chemistry, with a mechanism possibly linked to solar cycle.

A possible correlation between the solar cycles and upper-stratospheric temperatures at even lower pressures had been reported previously by Roques et al. (1994). Analyzing ground-based stellar occultations dating between 1983 and 1990, Roques et al. (1994) retrieved temperatures at pressures of 0.01–0.03 mbar. They found temporal variations of nearly 60 ± 20 K (∼150–210 K) that appeared roughly correlated in time with the sunspot numbers and Lyα flux of solar cycles 21 and 22, but lagging by roughly 1 yr.

To explain these possibly correlated temperature changes, Roques et al. (1994) proposed a mechanism in which the enhanced UV flux leads to the formation of aerosols, which in turn absorb solar radiation and heat the atmosphere.

Chemical and microphysical models have predicted the presence of hydrocarbon hazes (e.g., Romani & Atreya 1988, 1989; Moses et al. 1992; Romani et al. 1993; Toledo et al. 2020), and the presence of UV-absorbing hazes appears common in stratospheres and tropospheres of all the giant planets (e.g., West et al. 1986; Karkoschka & Tomasko 2009; Karkoschka 2011; Roman et al. 2013; Zhang et al. 2013). On Neptune, stratospheric aerosol layers were initially inferred by Smith et al. (1989) at heights of ∼150 km above the tropopause using Voyager high phase angle limb scans, placing them at the necessary sub-millibar pressure levels. Subsequent analysis of Voyager high phase angle data by Moses et al. (1995) indicated the presence of excess extinction consistent with submicron haze particles at pressures less than 15 mbar, including particles near 0.5 mbar. Although ethane, acetylene, and other relevant hydrocarbons do not condense at sub-millibar pressures, refractory hydrocarbons (such as PAHs) that are theoretically produced from higher-altitude ion chemistry can potentially condense at these levels (Dobrijevic et al. 2016). Stellar occultation data analyzed by Roques et al. (1994) also suggest the possibility of an aerosol layer even higher, extending from 50–100 μbar to 0.001–0.01 μbar. However, it remains to be seen whether photochemically produced aerosols at these pressures actually vary in time, and if so, whether solar heating of such aerosols can even provide heating rates large enough to account for the observed thermal variation. Moses et al. (1995) inferred a cumulative haze extinction optical depth of only ∼3 × 10⁻³ in Voyager's clear filter (280–640 nm) at pressures less than 15 mbar using Voyager imaging. Given these low optical depths, they concluded that, if typical, these hazes may not contribute significantly to the stratospheric heating rates.

Deeper in the atmosphere (≥∼100 mbar), solar-cycle-induced variation in the abundance of absorbing tropospheric aerosols would also be partly consistent with possible correlations with visual (472 and 551 nm) magnitudes (e.g., Hammel & Lockwood 2007b; Lockwood 2019). Detrended to remove the longer-term variability associated with seasonal changes (see Figure 24), Neptune's magnitude appeared dimmer near solar maxima and brighter near solar minima from the 1970s to the mid-1990s (cycles 20–22), with changes lagging about 3 yr behind the Sun (Lockwood & Thompson2002; Lockwood & Jerzykiewicz 2006; Aplin & Harrison 2016; Lockwood 2019). Baines & Smith (1990) proposed that these apparent changes were the result of solid-state chemical "tanning" of stratospheric aerosols by enhanced far-UV flux. The albedo change would naturally result in greater heating of the aerosol layers, and this heating would increase with height, presumably following a similar or longer lag relative to the aerosol production rates.

Interestingly, despite decades of seemingly rough agreement, the purported correlation between Neptune's albedo and the solar cycle seemed to break in solar cycle 23, when Neptune's visible albedo appeared brighter near the 2001 solar maximum—opposite to the previously inferred trend. The cause of this apparent break—if the solar correlation had indeed ever existed—has not yet been explained.

We note that the fluctuations in Neptune's magnitude determined by Aplin & Harrison (2016) (reproduced in Figure 24, top panel) appear very similar to the trend in mid-infrared brightness seen in the imaging. The early 2000s were the warmest in our limited record of the stratospheric temperatures, despite the planet's higher albedo at the time. Additionally, as we noted above, the early 2000s were also brighter at near-IR wavelengths owing to increased coverage of discrete high clouds, seen most clearly in methane band images (Karkoschka 2011; Roman et al. 2013). While the aforementioned albedo studies (Lockwood & Jerzykiewicz 2006; Aplin & Harrison 2017; Lockwood 2019) mostly considered the modulation of clouds and hazes in terms of disk-integrated visual magnitudes, the correlation of discrete cloud coverage was only previously hinted at from a few contemporary near-IR observations (Lockwood & Thompson 2002).

In a study of HST images, Karkoschka (2011) described the temporal variation in the discrete cloud abundance (represented by an albedo enhancing filling factor) as varying on roughly 5 yr timescales; however, as shown in Figure 24 (middle panel), the discrete cloud filling factor appears remarkably well correlated with the solar cycle Lyα, although these published cloud observations are currently limited to a single solar cycle. As noted, HST observations in recent years also show a dearth of discrete cloud cover, roughly at a time when the solar cycle reached its most recent minimum (between Solar Cycles 24 and 25 in 2019 December). The possible correlation between the two is interesting because it appears contrary to theoretical expectations.

While statistical analysis by Aplin & Harrison (2016) claimed to show that the Lyα flux was generally anticorrelated with Neptune's visual brightness, variations in the count of galactic cosmic rays were also claimed to have a statistically significant contribution to the magnitude. Galactic cosmic rays fluxes are largely inversely correlated with the solar cycle, reaching maximum fluxes near solar minimum (see Figure 24, bottom panel). Moses et al. (1989) previously proposed that Neptune's brightness variation could be explained by ion-induced cloud nucleation during solar minima by increasing the cloud coverage—similar to what is seen on Earth (Svensmark & Friis-Christensen 1997). However, the data appear to suggest the opposite in this case. Discrete cloud coverage peaked during the 2001 solar maximum and appears anticorrelated with the cosmic-ray rate ²⁰ for the years in which spatially resolved imaging is available. Regular measurements of cloud coverage for the prior maximum are unavailable, but Lockwood & Thompson (2002) note prominent "outbursts" in near-infrared observations in 1977 and 1986–1989—corresponding to times of solar minima and greater cosmic-ray counts. The most recent solar maximum was weak, but the discrete cloud coverage at the time still appeared greater than what was seen during the following solar minimum of 2019 (with the notable exception of the cloudy outburst in 2017; Molter et al. 2019).

Altogether, this suggests that the stratospheric temperatures and albedo at visible and near-IR wavelengths may potentially be correlated, and the solar cycle may provide the physical link. But observations are too limited to draw conclusions yet, and any solar effects on clouds and hazes have seemingly changed since the late 1990s for reasons unknown. It is possible that the expected UV "tanning" effect on the albedo acts independently of the discrete cloud coverage, but unless we are seeing a purely meteorological coincidence, the potential physical mechanisms possibly linking solar cycle variation, stratospheric temperatures, and cloud activity should be further examined. Given that possible changes occurred in the transition from southern spring to summer (2005 solstice; Figure 24, bottom panel), we might find that the aerosol response is somehow seasonally dependent or altered by sufficiently disruptive weather events. More observations over the next solar cycle will be key to assessing the consistency and physical mechanism behind this intriguing but uncertain potential relationship.