# Latitudinal variation in the abundance of methane (CH4) above the clouds in Neptune's atmosphere from VLT/MUSE Narrow Field Mode Observations

Patrick G. J. Irwin, Daniel Toledo, Ashwin Braude, Roland Bacon, Peter M. Weilbacher, Nicholas A. Teanby, Leigh N. Fletcher, Glenn S. Orton

<figure><img src="image_1.png" /></figure>

## To cite this version:

Patrick G. J. Irwin, Daniel Toledo, Ashwin Braude, Roland Bacon, Peter M. Weilbacher, et al. Latitudinal variation in the abundance of methane (CH4) above the clouds in Neptune's atmosphere from VLT/MUSE Narrow Field Mode Observations. Icarus, 2019, 331 (October), pp.69-82. 10.1016/j.icarus.2019.05.011. insu-03711413

HAL Id: insu-03711413

https://insu.hal.science/insu-03711413

Submitted on 2 Jul 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Latitudinal variation in the abundance of methane (CH₄) above the clouds in Neptune's atmosphere from VLT/MUSE Narrow Field Mode Observations

<figure><img src="image_4.png" /></figure>

Check for
updates

Patrick G.J. Irwin<sup>a,*</sup>, Daniel Toledo<sup>a</sup>, Ashwin S. Braude<sup>a</sup>, Roland Bacon<sup>b</sup>, Peter M. Weilbacher<sup>c</sup>, Nicholas A. Teanby<sup>d</sup>, Leigh N. Fletcher<sup>e</sup>, Glenn S. Orton<sup>f</sup>

## ARTICLE INFO

### Keywords:

Neptune

Neptune

Atmosphere

Atmospheres

Composition

## ABSTRACT

Observations of Neptune, made in 2018 using the new Narrow Field Adaptive Optics mode of the Multi Unit Spectroscopic Explorer (MUSE) instrument at the Very Large Telescope (VLT) from 0.48 to 0.93 µm, are analysed here to determine the latitudinal and vertical distribution of cloud opacity and methane abundance in Neptune's observable troposphere (0.1–~ 3bar). Previous observations at these wavelengths in 2003 by HST/STIS (Karkoschka and Tomasko 2011, Icarus 205, 674–694) found that the mole fraction of methane above the cloud tops (at ~ 2 bar) varied from ~ 4% at equatorial latitudes to ~ 2% at southern polar latitudes, by comparing the observed reflectivity at wavelengths near 825 nm controlled primarily by either methane absorption or H₂–H₂/H₂–He collision-induced absorption. We find a similar variation in cloud-top methane abundance in 2018, which suggests that this depletion of methane towards Neptune's pole is potentially a long-lived feature, indicative of long-term upwelling at mid-equatorial latitudes and subsidence near the poles. By analysing these MUSE observations along the central meridian with a retrieval model, we demonstrate that a broad boundary between the nominal and depleted methane abundances occurs at between 20 and 40°S. We also find a small depletion of methane near the equator, perhaps indicating subsidence there, and a local enhancement near 60–70°S, which we suggest may be associated with South Polar Features (SPFs) seen in Neptune's atmosphere at these latitudes. Finally, by the use of both reflectivity analysis and a principal component analysis, we demonstrate that this depletion of methane towards the pole is apparent at all locations on Neptune's disc, and not just along its central meridian.

## 1. Introduction

The visible and near-infrared spectra of both Uranus and Neptune are formed by the reflection of sunlight off its cloud layers, modulated mostly by the absorption of gaseous methane and also Rayleigh scattering at blue wavelengths. At wavelengths of low gaseous absorption, sunlight can penetrate to, and be reflected from, the deepest cloud layers, while at wavelengths of high absorption only sunlight reflected from the upper layers can be observed. Hence, such spectra can be inverted to determine the cloud structure as a function of location and altitude, but only if we know the vertical and latitudinal distribution of methane. For many years, planetary astronomers assumed that the vertical profiles of methane determined from Voyager 2 radio-

occultation observations could be used at all locations on these planets. However, HST/STIS observations of Uranus recorded in 2002 (Karkoschka and Tomasko, 2009) and similar observations of Neptune recorded in 2003 (Karkoschka and Tomasko, 2011) both revealed that the abundance of methane actually varies significantly with latitude on both planets. HST/STIS observes the 300–1000 nm range, which includes collision-induced absorption (CIA) bands of H₂–H₂ and H₂–He near 825 nm. Such spectra allow the cloud structure to be probed at wavelengths of either mostly CH₄ absorption or mostly H₂–H₂/H₂–He collision-induced absorption, allowing variations of CH₄ to be differentiated from cloud-top pressure variations (Fig. 1). Karkoschka and Tomasko (2009, 2011) found that the methane abundance near the main observable H₂S cloud tops at 2–3 bar varies from ~ 4% at

<footnote><sup>a</sup>Department of Physics, University of Oxford, Parks Rd, Oxford OX1 3PU, UK</footnote>

<footnote><sup>b</sup>CRAL - Observatoire de Lyon, 9 Avenue Charles André, 69230 Saint-Genis-Laval, France</footnote>

<footnote><sup>c</sup>Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany</footnote>

<footnote><sup>d</sup>School of Earth Sciences, University of Bristol, Wills Memorial Building, Queens Road, Bristol BS8 1RJ, UK</footnote>

<footnote><sup>e</sup>Department of Physics & Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK</footnote>

<footnote><sup>f</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: patrick.irwin@physics.ox.ac.uk (P.G.J. Irwin).</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Reproduction of part of Fig. 5 of Karkoschka and Tomasko (2011) showing HST/STIS observations of Neptune made in 2003 at two wavelengths near the 825 nm H<sub>2</sub>-H<sub>2</sub>/H<sub>2</sub>-He CIA band, one where H<sub>2</sub>-H<sub>2</sub>/H<sub>2</sub>-He opacity accounts for 87% of the opacity (left) and one where the opacity is entirely from gaseous methane (right). The images here have been processed to remove discrete features (Karkoschka and Tomasko, 2011). At wavelengths of low H<sub>2</sub>-H<sub>2</sub> and H<sub>2</sub>-He absorption (i.e. right) we see that the polar latitudes appear to be brighter than equatorial ones. If methane abundance were assumed to be constant with latitude, the brighter appearance of polar latitudes in the methane-image could be interpreted as being due to higher clouds there. However, the more uniform brightness seen in the H<sub>2</sub>-H<sub>2</sub>/H<sub>2</sub>-He-absorbing image (left) shows that the appearance of Neptune in the methane-absorbing images is actually due to lower methane abundance at polar latitudes.</figcaption></figure>

equatorial latitudes to ~ 2% at polar latitudes for both planets. These observations point to a very different overturning circulation of air at these pressure levels than that seen in the upper troposphere, where mid-infrared observations find cooler temperatures near the tropopause at mid-latitudes than at the equator and poles, indicative of air rising at mid-latitudes in the upper troposphere and subsiding elsewhere (e.g. de Pater et al., 2014; Fletcher et al., 2014). These observed latitudinal dependencies of methane also mean that cloud retrievals that do not include these variations will likely retrieve erroneously low cloud-top pressures near the poles of these planets.

Since the HST/STIS observations of Uranus in 2002 (Karkoschka and Tomasko, 2009), latitudinal variation in the abundance of methane above the clouds in Uranus's atmosphere has been found in: 1) a re-analysis of HST/STIS 2002 observations (Sromovsky et al., 2011); 2) IRTF/SpEx observations (800–850 nm) made in 2009 (Tice et al., 2013); 3) Gemini/NIFS observations made in 2010 (Irwin et al., 2012) at 1.45–1.85 µm, where the degeneracy between cloud-top pressure and methane was partially broken by the use of high spectral resolution; 4) HST/STIS observations made in 2012 combined with Keck/NIRC2 observations (2007–2011) (Sromovsky et al., 2014), where methane depletion at high latitudes was modelled as regions of downwelling air, using “proportionally descended gas” profiles (found by Sromovsky et al., 2011 to provide a good match to the HST/STIS 2002 observations of Karkoschka and Tomasko, 2009); and 5) further HST/STIS observations made in 2015, again combined with Keck/NIRC2 observations (Sromovsky et al., 2019), which again found depleted methane towards Uranus's north pole, but found that the absolute methane abundance was dependent on the assumed aerosol properties.

While the abundance of methane in Uranus's atmosphere has been extensively studied, the latitudinal variation of methane in Neptune's atmosphere has not been determined since HST/STIS observations in

2003 (Karkoschka and Tomasko, 2011). Hence, it is of great interest to see if the latitudinal variation of methane observed by HST/STIS in Neptune's atmosphere over 15 years ago remains the same today. In this paper, we present new visible/near-infrared ground-based spectral observations of Neptune made in 2018 with VLT/MUSE and use them to make an initial estimation, for the first time with ground-based observations, of the spatial distribution of methane above the clouds in Neptune's atmosphere. A more exhaustive analysis of clouds and cloud-top methane abundance, combining the VLT/MUSE observations with longer wavelength near-infrared observations from the VLT/SINFONI and Gemini/NIFS will be the topic of a future, follow-on paper.

## 2. VLT/MUSE observations

The Multi Unit Spectroscopic Explorer (MUSE) instrument (Bacon et al., 2010) at ESO's Very Large Telescope (VLT) in Chile is an integral-field spectrograph, which records 300 × 300 pixel 'cubes', where each pixel, or 'spaxel', contains a complete visible/near-infrared spectrum (0.48–0.93 µm) with a spectral resolving power of 2000–4000. MUSE has two viewing modes: a Wide-Field Mode (WFM) with a 60" × 60" field of view (FOV) and a Narrow-Field Mode (NFM), with a 7.5" × 7.5" FOV (giving a pixel size of 0.025"), which utilises Adaptive Optics (Arsenault et al., 2012; Ströbele et al., 2012) to improve the spatial resolution to < 0.1". At such a spatial resolution it is important to correct for differential atmospheric refraction between different wavelengths, for which MUSE utilises an atmospheric dispersion compensator that reduces residual shifts between wavelengths to less than a pixel. Commissioning observations of Neptune were made in this new NFM mode on June 19th 2018, which we present and analyse here. The NFM mode needs to reserve some wavelengths for its laser guide star and thus MUSE/NFM observations do not include the wavelength range 0.578–0.605 µm. Six observations of Neptune were made, two with exposures of 10 s each and four with exposures of 120 s, as summarised in Table 1. The appearance of Neptune at several wavelengths as seen in one of these exposures is shown in Fig. 2. The performance of the adaptive optics system is estimated for these observations to achieve a point-spread-function (PSF) with a full-width-half-maximum (FWHM) of 0.06" at the wavelengths considered here. No attempt was made to quantify the shape of the PSF from these observations (which would require the analysis of point object observations, which we did not have) and hence it is not known how significant a fraction of the total observed flux comes from outside the 0.06" FWHM. Hence, we avoided analysing observations made too close to the planet's limb, which could contain an unquantified contribution from space, and observed no biases at the extremes of the latitude range considered.

In contrast to the MUSE spectral resolution of 2.45 Å in NFM, the best available source of methane absorption in this spectral range are the band-model coefficients of Karkoschka and Tomasko (2010), to which we have fitted k-distributions using exponential-sum fitting. These data have a spectral resolution of 25 cm<sup>-1</sup> between 19,300 and 25,000 cm<sup>-1</sup> (0.518–0.4 µm), and 10 cm<sup>-1</sup> at wavenumbers less than 19,300 cm<sup>-1</sup> (0.518 µm). These resolutions equate to roughly 5 Å (i.e. 0.0005 µm) across the visible range and so the MUSE data cannot be

<table><caption>**Table 1**
MUSE/NFM Neptune observations.</caption><thead><tr><th>Obs. ID.</th><th>Date</th><th>Time (UT)</th><th>Exposure time</th><th>Airmass</th></tr></thead><tbody><tr><td>1</td><td>June 19th, 2018</td><td>09:37:41</td><td>10 s</td><td>1.0590</td></tr><tr><td>2</td><td>June 19th, 2018</td><td>09:40:10</td><td>10 s</td><td>1.0579</td></tr><tr><td>3</td><td>June 19th, 2018</td><td>09:43:21</td><td>120 s</td><td>1.0563</td></tr><tr><td>4</td><td>June 19th, 2018</td><td>09:47:44</td><td>120 s</td><td>1.0550</td></tr><tr><td>5</td><td>June 19th, 2018</td><td>09:52:06</td><td>120 s</td><td>1.0541</td></tr><tr><td>6</td><td>June 19th, 2018</td><td>09:56:26</td><td>120 s</td><td>1.0535</td></tr></tbody></table>

The full-width-half-maximum (FWHM) of the point spread function (PSF) of these adaptive optics observations is estimated to be 0.06" at 0.8 µm.
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** VLT/MUSE observation ‘3’ analysed in this study averaged over all wavelengths (i.e. white) and also over three different photometric spectral channels: a) Johnson V (507–595 nm); b) Cousins R (578–716 nm); and c) Cousins I (712–861 nm). In these plots, Neptune’s south pole is at bottom left. It can be seen that deeper cloud structure becomes increasingly clear at longer wavelengths, where we are less sensitive to Rayleigh scattering from overlying gas molecules.</figcaption></figure>

analysed using these coefficients at their native resolution. Hence, the data were first smoothed to make them more compatible with the available methane data and also improve the signal-to-noise (SNR) ratio. We chose to smooth the MUSE data to the resolution of the IRTF/SpEx instrument, which has a triangular instrument function with FWHM = 20 Å (i.e 0.002 µm), sampled at 0.001 µm. We find that this spectral resolution captures the essential shape of the observed spectra very well and has significantly better SNR than MUSE’s native resolution, giving SNRs ranging from ~ 80 at the reflectance peak to ~ 40 at methane-absorbing wavelengths.

The MUSE data were reduced with the VLT/MUSE pipeline v2.5.1 (Weilbacher et al., 2014) and flux-calibrated using the spectro-photometric standard star EG274, observed on the same night. For this study, spectra were extracted from the MUSE data from 64 individual spectral pixels, or ‘spaxels’, spread along the central meridian of Neptune, but not too close to the limb to experience significant mixing with space, as shown in Fig. 3. Twelve spaxels, roughly equally spaced in latitude, were selected for detailed description, indicated by the white dots in Fig. 3 and listed in Table 2. Cube ‘3’ was chosen for this analysis as it was representative of all the 120 s-integrated observations (‘3’–‘6’), and the SNR was sufficient that we did not need to co-add multiple cubes. We originally attempted to fit these data using the noise values provided by the pipeline, but found we could not fit to this level of precision, achieving $\chi^2/n_y$ values of ~ 4 at best. Hence, we doubled these errors in the fits presented here, where we achieve a minimum $\chi^2/n_y$ of ~ 1. We attribute this inability to fit to within the original noise estimate to ‘forward-modelling’ deficiencies in our model, arising from factors ranging from errors including: 1) errors in our gaseous absorption data; 2) systematic deficiencies in our model set up; 3) the smoothing process itself introducing unaccounted-for covariances into the smoothed spectra; and 4) any residual calibration errors.

To check the flux calibration, Fig. 4 compares the measured disc-averaged reflectance spectrum for Cube ‘3’ with the disc-averaged albedo spectrum of Karkoschka (1994), showing the generally good agreement (~ 10%). Neptune has a highly dynamic atmosphere and it is not surprising that the disc-averaged albedo spectrum we have measured in 2018 differs slightly from the spectrum of Karkoschka (1994), which was observed in 1993. We can see in Fig. 3 that the reflectance spectrum changes significantly with both position on the disc and local cloud conditions, as can be seen when we compare these spectra with the observed centre-of-disc reflectance spectrum, which is noticeably brighter. Fig. 4 also compares the MUSE spectra with recently acquired (September 2018) HST/WFC3 observations near the disc centre (HST/OPAL (e.g. Simon et al., 2015) programme; Simon — private communication). Here, we see reasonably good agreement between the VLT/MUSE and HST/WFC3 observations, which together with our reasonable correspondence with the albedo spectrum of Karkoschka (1994) leads us to conclude that the photometric calibration of our data is reliable and consistent with previous observations.

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** VLT/MUSE observation ‘3’ averaged over the Cousins I photometric channel (712–861 nm), plotted as apparent size and angle on the sky, and showing the position of the sampled ‘spaxels’ along the central meridian (black squares). The spaxels selected for detailed description in subsequent plots and tables are indicated by the white squares. Also plotted is a reference latitude-longitude grid, with lines of latitude (planetocentric) separated by 20° and lines of longitude separated by 45°. The equator is indicated by the dashed line, rather than dotted line, and the south pole is at bottom left.</figcaption></figure>

### 3. Radiative-transfer analysis

The vertical cloud structure and methane abundance were retrieved from these observations using the NEMESIS radiative transfer and retrieval code (Irwin et al., 2008). NEMESIS models planetary spectra either by using a line-by-line model, or by using the correlated-k approximation (e.g. Lacis and Oinas, 1991). Since our methane absorption data is from a band model (Karkoschka and Tomasko, 2010) we used the method of correlated-k here, constructing k-tables from the band data using exponential-sum fitting. As with our previous analyses (e.g. Irwin et al., 2016, 2011), we used the matrix-operator multiple-scattering model of Plass et al. (1973) to simulate these reflected sunlight spectra, modelling the atmosphere with 39 levels spaced between ~ 10 and ~ 0.01 bar and using either 5 or 9 angles (both upward and downward) in a Gauss-Lobatto zenith-angle quadrature scheme, with the number of required azimuth components in the Fourier decomposition determined from the maximum of the reflected or incident-solar zenith angles. The collision-induced absorption of H$_2$–H$_2$ and H$_2$–He near 825 nm was modelled with the coefficients of Borysow et al. (1989a,b, 2000)<fnref n="1" />. In retrieval tests, we found that the H$_2$–H$_2$/H$_2$–He CIA coefficients for hydrogen with thermally-equilibrated ortho:para hydrogen ratio were much more consistent with the observed spectra than the ‘normal’ 3:1 ortho:para hydrogen ratio (more appropriate for planets such as Jupiter, where air is upwelled rapidly from deeper, hotter levels) and so we used these coefficients in our retrievals. A similar conclusion was reached by Baines and Smith (1990) and Gautier et al. (1995). Rayleigh scattering was also included, where for each gas the contribution was calculated as described in the Appendix. At these wavelengths, however, effects of polarization and Raman scattering were justifiably ignored (e.g. Sromovsky, 2005). To analyse the measured radiance spectra with our radiative-transfer model we used the solar spectrum of Chance and Kurucz (2010), which was first smoothed with a triangular line shape of FWHM = 0.002 µm to make it compatible with our smoothed spectra. To calculate the reflectivity, we also needed Neptune’s distance from the Sun on the date of observation,

<footnote n="1">¹ Fortran programs for calculating the H$_2$–H$_2$ and H$_2$–He coefficients were downloaded from https://www.astro.ku.dk/~aborysow/programs/index.html.</footnote>
<table><caption>**Table 2**

Retrieval results at selected areas (12 of 64) considered on Neptune's disc.</caption><thead><tr><th>Area</th><th>Latitude</th><th>p<sub>1</sub><br>(bar)</th><th>f<sub>CH<sub>4</sub></sub><br>%</th><th>p<sub>C</sub><sup>*</sup><br>(bar)</th><th>τ<sub>C</sub></th><th>FSH<sub>C</sub></th><th>p<sub>H</sub><sup>*</sup><br>(bar)</th><th>τ<sub>H</sub></th><th>FSH<sub>H</sub><sup>*</sup></th><th>χ<sup>2</sup>/n<sub>y</sub></th></tr></thead><tbody><tr><td>1</td><td>85.5'S</td><td>2.06</td><td>3.06 ± 0.57</td><td>4.3</td><td>14.13</td><td>0.312</td><td>0.038</td><td>0.20</td><td>0.1</td><td>1.99</td></tr><tr><td>2</td><td>74.4'S</td><td>1.88</td><td>4.06 ± 0.49</td><td>4.3</td><td>24.81</td><td>0.289</td><td>0.038</td><td>0.27</td><td>0.1</td><td>2.62</td></tr><tr><td>3</td><td>65.4'S</td><td>2.26</td><td>2.62 ± 0.57</td><td>4.3</td><td>6.71</td><td>0.395</td><td>0.038</td><td>0.23</td><td>0.1</td><td>2.84</td></tr><tr><td>4</td><td>54.2'S</td><td>2.14</td><td>3.03 ± 0.46</td><td>4.3</td><td>7.40</td><td>0.413</td><td>0.038</td><td>0.29</td><td>0.1</td><td>3.00</td></tr><tr><td>5</td><td>45.2'S</td><td>2.10</td><td>3.54 ± 0.57</td><td>4.3</td><td>7.30</td><td>0.424</td><td>0.038</td><td>0.30</td><td>0.1</td><td>3.17</td></tr><tr><td>6</td><td>35.8'S</td><td>1.70</td><td>> 4.74</td><td>4.3</td><td>33.75</td><td>0.297</td><td>0.038</td><td>0.47</td><td>0.1</td><td>2.82</td></tr><tr><td>7</td><td>25.9'S</td><td>1.80</td><td>5.90 ± 1.07</td><td>4.3</td><td>21.56</td><td>0.325</td><td>0.038</td><td>0.44</td><td>0.1</td><td>3.18</td></tr><tr><td>8</td><td>15.1'S</td><td>2.20</td><td>5.89 ± 1.00</td><td>4.3</td><td>6.91</td><td>0.399</td><td>0.038</td><td>0.31</td><td>0.1</td><td>1.63</td></tr><tr><td>9</td><td>5.1'S</td><td>2.20</td><td>5.55 ± 1.13</td><td>4.3</td><td>6.65</td><td>0.411</td><td>0.038</td><td>0.31</td><td>0.1</td><td>1.85</td></tr><tr><td>10</td><td>4.8'N</td><td>2.49</td><td>4.46 ± 1.16</td><td>4.3</td><td>4.95</td><td>0.392</td><td>0.038</td><td>0.23</td><td>0.1</td><td>1.68</td></tr><tr><td>11</td><td>14.6'N</td><td>2.58</td><td>4.88 ± 1.04</td><td>4.3</td><td>4.51</td><td>0.387</td><td>0.038</td><td>0.20</td><td>0.1</td><td>1.21</td></tr><tr><td>12</td><td>26.1'N</td><td>2.07</td><td>> 9.67</td><td>4.3</td><td>14.51</td><td>0.305</td><td>0.038</td><td>0.24</td><td>0.1</td><td>1.17</td></tr></tbody></table>

Notes: variables marked with an asterisk * were all fixed to values found through our limb-darkening analysis; $p_1$ is the pressure (bar) where the cloud opacity to space is unity; $f_{CH_4}$ is the retrieved cloud-top CH$_4$ mole fraction (%), i.e. the mole fraction at pressure $p_1$; $p_C$ is the base pressure (bar) of the lower cloud; $\tau_C$ is the opacity of the cloud (at 0.8 µm); $FSH_C$ is the fractional scale height of the cloud; $p_H$ is the base pressure (bar) of the haze; $\tau_H$ is the opacity of the haze (at 0.8 µm); $FSH_H$ is the fractional scale height of the haze; $\chi^2/n_y$ is the chi-squared statistic of the fit. Where the methane mole fraction at the cloud-top pressure $p_1$ is already limited to not exceed the saturated vapour mole fraction, we instead quote a lower limit for the methane abundance just below the cloud. The assumed refractive index spectra (deduced from the limb-darkening analysis) of the cloud and haze particles are shown in Fig. 6.

<figure><img src="image_8.png" /></figure>

which was 29.94 AU.

The reference temperature and abundance profile used in this study is the same as that used by Irwin et al. (2014) and is based on the ‘N’ profile determined by Voyager-2 radio-occultation measurements (Lindal, 1992), with $He:H_2 = 0.177$ (15:85) including 0.3% mole fraction of $N_2$. This profile was determined at equatorial latitudes and thus its use as a ‘global’ *a priori* is moot. At the level of the main clouds

**Fig. 4.** Reflectivity spectra of Neptune observed with VLT/MUSE compared with the albedo spectrum of Karkoschka (1994) and recent observations with HST/WFC3 (September 2018). The top plot shows the whole VLT/MUSE spectral range, while the bottom plot shows the 0.77–0.93 µm region only. Three different MUSE spectra are shown: 1) a disc-averaged spectrum at MUSE’s native spectral resolution; 2) a disc-averaged spectrum smoothed to a resolution of 0.002 µm; and 3) a measured, smoothed MUSE spectrum at the centre of Neptune’s disc. Recent photometric observations by HST/WFC3 (September 2018) are shown in cyan for comparison, with the approximate width of the filter channels indicated by the horizontal lines.

($p \sim 2$ bar), it has for many years been thought that Neptune's atmosphere should become barotropic, in which case we would not expect large latitudinal variations in temperature at these levels. However, Tollefson et al. (2018) have suggested that a vertical shear seen in the zonal wind speeds inferred from tracking cloud features in H-band (1.4–1.8 µm) and K-band (1.8–2.2 µm) observations, may be caused by a significant variation in atmospheric temperatures from equator to
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Fit to the average observed I/F spectra in Observation ‘3’ in the 5–10°S latitude band sampled at six viewing zenith angles corresponding with angles in the nine-zenith-angle Gauss-Lobatto quadrature scheme used in our analysis. Spectra are plotted with same vertical scale but successively offset by a reflectivity difference of 0.2 for clarity. The black lines are the observed spectra (with errors indicated in grey) while our fits are shown in red. At each angle, the difference is shown below in blue, together with the estimated measurement errors, again in grey. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

mid-latitudes to poles. Depending on their assumptions for the latitudinal variation of methane, Tollefson et al. (2018) find that latitudinal variations in temperature of up to 2–15 K may be present at pressures greater than 1 bar. We will return to the potential effect of such temperature variations on our methane retrievals later in this section. Finally, our *a priori* methane profile was set to have a deep mole fraction of 4% (Karkoschka and Tomasko, 2011), limited to 100% relative humidity in the upper troposphere and constrained from exceeding a mole fraction of $1.5 \times 10^{-3}$ (Lellouch, 2010) in the stratosphere.

Rather than fit the entire spectrum we focussed only on the 0.77–0.933 µm region. This was done for a number of reasons: 1) the number of wavelengths was reduced from 459 to 161, which increased retrieval speeds; 2) by limiting the wavelength range to the single reflectance peak containing the H₂–H₂ CIA feature we could more reasonably assume single particle size distributions; and 3) by analysing longer wavelengths only we are much less sensitive to the effects of Raman scattering, Rayleigh scattering (from air and stratospheric aerosols), and the effects of polarization (Sromovsky, 2005).

We initially attempted to model the spectrum using a continuous profile of identical particles, and found good fits to individual spectra

with either small (0.1 µm) or larger (1.0 µm) particles. However, we found this model was not able to match the observed limb-darkening/ limb-brightening seen at different wavelengths very well and so we instead adopted the parameterised model used by Irwin et al. (2016) to model VLT/SINFONI and Gemini/NIFS H-band observations of Neptune, which was found to provide good limb-darkening/limb-brightening behaviour. In this model, particles in the troposphere are modelled with a cloud with a variable base pressure (set to an *a priori* pressure of ~ 3–4.5 bar) and a scale height retrieved as a fraction of the pressure scale height, while scattering from haze particles is modelled with a second layer, with variable base pressure at ~ 0.1 bar and fixed fractional scale height of 0.1. The particles in the tropospheric cloud were modelled here with a standard Gamma size distribution with mean radius 1.0 µm and variance 0.05, while the stratospheric haze particles were modelled with a Gamma size distribution with mean radius 0.1 µm and the same variance of 0.05. Following Irwin et al. (2015), the real refractive index of both particle types was set to 1.4 (broadly typical refractive index of H₂S and NH₃ liquid/ice: Havriliak et al., 1955, Martonchik and Orton, 1994, Wong and Anderson, 1974) at a reference wavelength of 0.8 µm and NEMESIS used to retrieve the imaginary refractive index spectrum. The *a priori* imaginary refractive index spectrum was set to 0.01 ± 0.005 and sampled at every 0.05 µm between 0.75 and 0.95 µm, with a ‘correlation length’ of 0.1 µm set in the covariance matrix to ensure that the retrieved spectrum varied reasonably smoothly with wavelength. At each iteration of the model, the real part of the particles’ refractive index spectrum at wavelengths other than the reference wavelength of 0.8 µm was computed using the Kramers-Kronig relation (Sheik-Bahae, 2005). Self-consistent scattering properties were then calculated using Mie theory, but the Mie-calculated phase functions were approximated with combined Henyey-Greenstein functions at each wavelength to smooth over features peculiar to perfectly spherical scatterers such as the ‘rainbow’ and ‘glory’, as justified by Irwin et al. (2018) (although such features are not actually present for the particle sizes we have assumed at these wavelengths in our Mie-modelled phase functions). We found that this two-cloud-component model was able to adequately represent the limb-darkening/limb-brightening of Neptune’s atmosphere at different wavelengths and so would not confuse such effects with the latitudinal changes in our retrieved composition and cloud estimates, since we observe higher latitudes at higher emission angles.

To set up a baseline model able to account correctly for limb-darkening/limb-brightening, we first of all selected observations from Cube ‘3’ in the latitude band 10°S–5°S, plotted these as a function of observing zenith angle and extracted the general limb-darkening/limb-brightening curves at all wavelengths. At each wavelength, the observations were fitted with smooth reflectance versus emission angle curves, and sampled at six emission zenith angles between 12.5 and 65° (Fig. 5), corresponding to six of the zenith angles in our nine-zenith-angle quadrature scheme, and not coming too near the disc edge to minimise PSF uncertainty errors. In our retrieval we fitted for: the base pressures and opacities of both cloud/haze layers, the fractional scale height of the cloud, the imaginary refractive index spectra of both particle types and a scaling factor for the deep methane abundance. The methane profile was limited to not exceed 100% relative humidity and to not exceed a mole fraction of $1.5 \times 10^{-3}$ in the stratosphere. This gave 16 freely variable parameters in all. We should stress here that although in this scheme we retrieve a ‘deep’ methane mole fraction, what we are actually sensitive to (and what we actually report) is the abundance of methane just above the main cloud deck, i.e. at 2–3 bar, rather than the ‘bulk’ deep abundance. It is not plausible for the deep methane mixing ratio on Neptune to have a large latitude dependence as this would result in a large gradient in density if not compensated for by a horizontal thermal gradient. The thermal wind equation is actually a density wind equation and if molecular weight gradients existed to great depths then there would also be vertical wind shears extending to great depths, which if integrated to the cloud tops would lead to a huge
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Retrieved cloud structure and imaginary refractive index spectra from our limb-brightening/limb-darkening analysis. The left hand panel shows our fitted cloud/aerosol vertical profiles (in opacity/bar at 0.8 µm), with the vertical distribution of 1.0-µm cloud particles shown as the solid curve, and the distribution of the 0.1-µm haze particles shown as the dotted-line curve. The middle and right hand panels show the retrieved imaginary refractive index spectra of the cloud and haze particles, respectively. Here, the retrieved spectra and errors are shown by the solid lines and light grey regions. The dark grey regions show the assumed *a priori* imaginary refractive index range. Here, we see that the cloud particles are retrieved to have low imaginary refractive index (and are thus quite scattering), while the haze particles are retrieved to have high imaginary refractive index (and are thus quite absorbing).</figcaption></figure>

latitude dependence in the cloud-top wind speed, which is not ob-
served. It is thus likely that any retrieved latitudinal variations in the
‘deep’ methane abundance are restricted to the upper troposphere, as
suggested by Karkoschka and Tomasko (2011) for both Neptune and
Uranus in their Fig. 10. The case for a similar depletion profile on Ur-
anus is made by Sromovsky et al. (2019, 2014), based both on theo-
retical grounds and on the improved fit quality obtained for upper
tropospheric depletion rather than deep depletion.

Our fits to the measured spectra at the six sampled zenith angles are shown in Fig. 5 and it can be seen that very good fits were achieved at all zenith angles ($\chi^2/n_y = 1.62$) giving confidence that our model was able to account for limb-darkening/limb-brightening. Our retrieved cloud/haze structure is shown in Fig. 6, together with our retrieved imaginary refractive index spectra. Here we see, as was found in the H-band from VLT/SINFONI observations (Irwin et al., 2016), that the limb-darkening/limb-brightening is best matched with tropospheric particles with low imaginary refractive index (giving reasonably high single-scattering albedos of $\bar{\omega} \sim 0.9$), while in the stratosphere, the particles are retrieved to have higher imaginary refractive index, giving lower single-scattering albedos of $\bar{\omega} \sim 0.74$; this higher imaginary refractive index for the haze particles is necessary to lower the single-scattering albedo sufficiently to match the observed limb-darkening at methane-absorbing wavelengths, as previously noted for H-band observations by Irwin et al. (2011) and Irwin et al. (2016). While our model was able to constrain the pressure level of the haze reasonably well, it was less able to constrain the base pressure of the cloud. We thus attempted retrievals with a number of different base pressures and finally settled on a *priori* pressure of 4.2 bar, which gave the reasonably good fit (i.e. $\chi^2/n_y = 1.62$) to the observations shown in Fig. 5.

Using the assumption that the same cloud scattering properties can
be used at all latitudes, we then fixed the complex refractive indices of
both particle types, fixed the cloud/haze base pressures and fitted the
central meridian spectra with just four parameters: the opacities of the
cloud and haze layers, the fractional scale height of the tropospheric
cloud, and the deep methane mole fraction scaling parameter. In this
analysis we used our 5-zenith-angle quadrature scheme (for reasons of
speed) and Fig. 7 shows the fits we were able to achieve to the observed
reflectance spectra with our model at our twelve representative loca-
tions, spread along the central meridian and summarised in Table 2.

Here we can see that we achieve a reasonably good fit to the observed
spectra with $\chi^2/n_y \sim 1-3$ at most locations, increasing to larger values
between 20 and 40°S, where the mid-latitude cloud belt is seen. Also
shown in Fig. 7 are lines joining the observed reflectivities at 827 and
833 nm, extended slightly to enable the gradient to be more clearly
seen. These lines will be explained later.

Fig. 8 shows our fitted methane profiles. In most cases, we have plotted the retrieved deep methane abundance and errors directly, but at some locations we retrieved very large methane abundances and very large errors. In these locations, we found that the pressure level where the retrieved deep methane abundance intersects the saturated methane mole fraction profile lay at pressures greater than the pressure where the cloud optical depth (i.e. opacity to space) is unity. Hence, at these locations, since the methane profile is constrained to not exceed 100% relative humidity, the retrieval becomes insensitive to the deep methane abundance and we can only derive a lower limit on the deep methane abundance, which we have set here to be the saturated mole fraction at the pressure where the cloud optical depth is unity. Apart from these few locations we can see that the deep methane abundances are reasonably well constrained in this analysis and decrease significantly towards the south pole, with the retrieved latitudinal differences mostly exceeding the retrieval errors. As we noted earlier, we used a single temperature profile at all latitudes in this analysis. If the real temperatures were slightly warmer in the ~ 0.5–1.5 bar level, the saturated vapour pressure would be higher and thus our retrieved methane profiles would be constrained by different saturated mole fraction profiles. As a result, we would retrieve different deep mole fractions to achieve the same column amount of methane above the cloud top. To test this we computed the saturated methane mole fraction profiles using the latitude-dependent temperature profiles fitted from mid-infrared observations by Fletcher et al. (2014), which we also show in Fig. 8. Here, we see that since the temperature profiles retrieved by Fletcher et al. (2014) are generally cooler than our a priori assumptions the saturated methane mole fractions are even smaller. Hence, to achieve the same column abundance of methane above the clouds, we would have to either lower the cloud top altitude, or increase further the deep methane abundance. However, mid-IR observations are not very sensitive to the temperatures near the 1-bar level and so, even though Fletcher et al. (2014) retrieve cooler
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Fits to the observed spectra at our twelve sample locations (summarised in Table 2). Here the observed spectra, including error bars, are shown in grey, while our fitted spectra are shown in red. At the top of each plot, we also show the difference between the measured and fitted spectra (in blue), multiplied by 5 and offset for clarity, together with the error range in grey (also multiplied by 5). It can be seen that we achieve excellent fits with $\chi^2/n$ of the order of 1.0–3.0. The near- horizontal lines near the reflectance peak link the observed reflectivities at our chosen test wavelengths of 827 and 833 nm, as explained in the text. Positive slopes indicate less methane, negative slopes indicate more methane. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

tropopause/upper-troposphere temperatures, we believe that the cooler temperatures inferred near 1 bar are most likely to be caused more by the vertical smoothing imposed on the profiles, rather than true lowering of the temperatures here although, as we noted earlier, Tollefson et al. (2018) find that the temperatures in the 1–3 bar region may vary with latitude by several K, depending on the assumed latitudinal variation of methane. If this is the case, then this would have an impact on our retrieved deep mole fractions of methane. Hence, although we

detect a clear decrease in methane abundance above the cloud as we go from the equator to the south pole, the absolute abundances are more difficult to constrain.

Fig. 9 shows the latitudinal dependence of the retrieved cloud and methane profiles at all 64 ‘spaxels’ analysed along the central meridian as contour plots, where we have plotted the results for all sampled latitudes, not just the twelve shown in Figs. 7 and 8. Fig. 9 also shows the fitted $\chi^2/n_y$ at all latitudes together with the latitudinal variation in the
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Fitted methane mole fraction profiles at our sample locations, summarised in Table 2. Our methane profile is parameterised to have a constant mole fraction at deep pressures limited to not exceed 100% relative humidity at lower pressures, with the saturated mole fraction calculated from our assumed reference temperature profile and shown here as the blue line. To aid comparison, the methane profile retrieved for the reference pixel area ‘9’ at 5°S is shown in red in all panels. In these plots the uncertainty of the profiles (modified as described in the text) are indicated in grey. A darker grey has been used to indicate the profile error for the reference pixel area ‘9’ at 5.1°S. We can see that polewards of 35°S there is significantly less methane than at equatorial latitudes. The dashed blue lines show the saturated methane mole fraction corresponding to the temperature profiles retrieved at these latitudes by Fletcher et al. (2014). Since these temperature profiles are generally cooler than the reference temperature profile, the saturated methane mole fractions are generally smaller. The horizontal black lines show the pressure where the cloud optical depth at 0.8 µm (i.e. opacity to space) is unity, while the horizontal, dashed red line shows this pressure level for the reference pixel area ‘9’. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

retrieved deep abundance of methane. We can see that we achieve good fits to these data at all locations and that we detect a clear and significant decrease in cloud-top methane abundance from 4–5 % at near-equatorial latitudes to 3–4 % polewards of ~ 40°S. In addition, Fig. 9 also shows the methane latitudinal variation determined from HST/STIS observations on 2003 by Karkoschka and Tomasko (2011),

showing reasonable correspondence. This indicates that our results are broadly consistent with this earlier determination.

4. Reflectivity spectral analysis

To understand how the abundance of methane can be extracted
<figure><img src="image_13.png" /></figure>

**Fig. 9.** Latitudinal dependence of retrieved cloud opacity/bar (at 0.8 µm), deep methane mole fraction and $\chi^2/n_y$. The top panel shows a contour plot of the fitted cloud opacity/bar profiles at all 'spaxels' sampled along the central meridian (brighter regions indicate greater cloud density, with black indicating an opacity/bar of zero and white indicated an opacity/bar of 50). Additional contour lines are drawn for opacity/bar values of 0.01, 0.1, 10 and 50. The bottom panel shows the latitudinal variation of the retrieved deep methane mole fraction, with error bars, and the $\chi^2/n_y$ of the fits. Where the cloud top intersects with the vapour pressure curve of methane, we can only estimate a lower limit of the deeper mole fraction, and indicate this with an arrow symbol instead of an error bar. We see a clear decrease in methane abundance from 5–6% at near-equatorial latitudes to 3–4% polewards of 20–40°S. The red line in the top plot indicates the cloud top pressure (i.e. level where overlying cloud opacity at 0.8 µm is unity). The cloud contour map indicates the main cloud top to lie at similar pressure levels at all latitudes and has a cloud-top pressure of ~ 2 bar. We can also see an increase in stratospheric haze at 20–40°S, associated with the cloudy zone and then clearing towards the north and south. The bottom plot also shows the methane abundance estimated by Karkoschka and Tomasko (2011) (red), together with the scaled central meridian variation of the contribution of the 3rd Empirical Orthogonal Function we derive with our PCA analysis (blue) in Section 5. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

<figure><img src="image_14.png" /></figure>

**Fig. 10.** Top: Observed MUSE spectrum at disc centre (~ 20°S, cross symbols) near 825 nm and fitted spectrum (black line), where the cloud opacity profiles and cloud imaginary refractive indices have been allowed to vary, but the methane profile has been fixed to the *a priori* methane abundance profile, which has a deep mole fraction of 4%. The modelled spectrum recalculated with half or double this deep methane abundance is shown by the red and blue lines, respectively. Also indicated by the vertical lines are the wavelengths chosen for reflectance differencing, which probe regions of similar reflectivity (and thus similar pressure levels) for the reference methane abundance of 4%, but whose opacities are predominantly determined by H$_2$–H$_2$/H$_2$–He CIA and methane absorption, respectively. Bottom: change in modelled reflectivity when the methane abundance is increased (blue) or decreased (red) by 10% from the nominal value of 4%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

from these data, Fig. 10 shows our fit to the spectrum near the disc centre assuming 4% deep methane mole fraction and then spectra re-calculated with the methane abundance either halved or doubled. It can

be seen that as we do this, the modelled reflectivity near the centre of the H$_2$–H$_2$/H$_2$–He CIA bands does not change much, but the reflectivity either side of the 0.820–0.835 µm reflectance peak is significantly
<figure><img src="image_15.png" /><figcaption>Fig. 11. Reflectivity differencing method applied to all six individual MUSE observations, indicated by the labels ‘1’–‘6’, showing the reflectivity at 827 nm minus the reflectivity at 833 nm, i.e. $R_{827} - R_{833}$. Here, as in Fig. 3 the south pole is at bottom left. The noisiness of the first two observations (which were only integrated for 10 s, rather than 120 s) is clear. For the remaining observations ‘3’–‘6’, these were measured with Neptune centred at different positions within the MUSE FOV and with different position angles (corrected for here). It is clear that residual photometric calibration ‘striping’ artefacts are apparent when comparing just two wavelengths as we do here (i.e. the horizontal and vertical stripes). However, for all these observations, it can be seen that the images are slightly darker at the bottom left than top right, indicated a lower methane spectral signal near the south pole. Approaching the south pole a faint arc (or perhaps even part of a ring) of enhanced $R_{827} - R_{833}$ can just be discerned in cubes ‘3’ and ‘6’, indicating possibly enhanced methane at 60–70° S.</figcaption></figure>

altered. In essence, as we increase the abundance of CH₄ the peak becomes slightly narrower, especially on the long-wave side. We wondered if it might be possible to elucidate this phenomenon directly from the observations, without having to run a full and time-consuming retrieval at every location on Neptune's disc. We can see in Fig. 10 that at wavelengths 827 and 833 nm there is roughly equal reflectance in our reference spectrum, indicating that we probe reflectance from roughly the same pressure level, but with the opacity of the former wavelength dominated by H₂–H₂/H₂–He CIA absorption and the latter wavelength dominated by CH₄ absorption. Hence, as the methane abundance rises, $R_{827} - R_{833}$ increases, and as the methane abundance falls, $R_{827} - R_{833}$ decreases. The observed $R_{827} - R_{833}$ differences are indicated by the short, straight lines shown in Fig. 7, where we can see that $R_{827} - R_{833}$ is roughly zero near the equator, but becomes noticeably negative as we approach the south pole, consistent with the lower retrieved abundance of CH₄ found here. Fig. 11 shows mapped images of the observed reflectivity $R_{827} - R_{833}$ difference applied to all six sets of observations. In this scaling we expect larger reflectivity differences to be associated with enhanced abundances of methane and we can indeed see larger signals near mid-latitudes, falling towards the south pole. We can also see a faint enhancement circling the south pole. As can be seen ‘striping’ artefacts are visible in some of these images, due to calibration errors and there was initially a concern that this might have affected our retrievals. However, it can be seen that the striping is not too severe in observation ‘3’, chosen for our retrieval analysis, and we could find no significant impact in our retrievals along the central meridian when crossing from one ‘stripe’ to the next. Hence, we conclude that while these small calibration errors are quite noticeable in these radiance-differencing images, which compare two wavelengths, they do not significantly affect our retrievals that simultaneously fit to multiple wavelengths.

In Fig. 12, we show the averaged reflectivity image at 827 nm, averaged over observations ‘3’, ‘4’, and ‘6’, compared with the averaged image at 833 nm and then the averaged difference image. We again see larger $R_{827} - R_{833}$ signals near mid latitudes, falling towards the south

pole, and also a slight brightening of the difference at 60–70° S, indicating possible local enhancement of methane at this latitude. Looking at the reflectivity images (which we have here enhanced to try to bring out the weaker features) we see a slightly brighter zone at 60–70° S and then lower reflection polewards of 75°S in the 827 nm image, whereas the methane image at 833 nm is relatively featureless. Hence, we appear to detect either a slight thickening or vertical extension of the main cloud from 60 to 70° S which is masked in the methane image by an accompanying increase in methane abundance and thus absorption. We note that an enhancement in methane abundance at 60–70° S can be discerned in our retrieved deep CH₄ abundances in Fig. 9.

## 5. Principal Component Analysis

The weak and rather noisy nature of the CH₄ reflectivity difference map (Fig. 11) led us to explore alternative methods of mapping the spectral reflectivity signal of CH₄ in Neptune's atmosphere. We turned to the technique of Principal Component Analysis (PCA) (e.g. Murtagh and Heck, 1987), used with great success in modelling visible/near-IR Jovian spectra by Dyudina et al. (2001) and Irwin and Dyudina (2002). The basic principle of Principal Component Analysis is that the variance of a set of observed spectra, in this case the varying spectra observed over Neptune's disc, can be decomposed into a set of Empirical Orthogonal Functions (EOFs) that form a linear basis from which any spectrum in the set, $y(\lambda)$, can be reconstructed from a linear combination of the EOFs, $E_i(\lambda)$, as $y(\lambda) = \sum_i \alpha_i E_i(\lambda)$, where the coefficients, $\alpha_i$, describe the relative proportions of the different EOFs in the combined spectrum. The derived EOFs have with them an eigenvalue, $e_i$, and the EOFs are usually ranked in order of decreasing $e_i$. With this ordering it is found that most of the variance can be accounted for by the first EOF (i.e. the one with the largest derived eigenvalue), with decreasingly significant contributions from higher EOFs. The derived EOFs do not necessarily correspond to anything physically significant, but under certain circumstances they can sometimes correspond to real variables.

In this case, since we are interested in searching for the spectral signatures of H₂–H₂/H₂–He and methane absorption, which have similar absorption strengths near 0.83 µm, we analysed the observed Neptune spectra covering the wavelength range 0.80–0.86 µm. The results are shown in Fig. 13. As can be seen the eigenvalues of the fitted EOFs fall rapidly with increasing EOF number and we can also see that the spatial distribution of the fitted weighting coefficients, $\alpha_i$, become increasingly noisy. In fact, we found that the first three EOFs effectively encapsulate all the significant information. We can see that EOF 1 is mostly flat, and that its spatial map corresponds almost exactly with the I/F appearance of Neptune over these wavelengths. Hence, this EOF appears to encapsulate the overall observed reflectivity variation. EOF 2 has a similar spectrum (albeit inverted) and its spatial distribution is almost the inverse of the spatial distribution for EOF 1, although its spatial structure exhibits finer structure, including a dark ‘halo’ about the south pole at ~ 80° S. EOF 3 also contains significant spectral variation, but its spatial distribution is very different from that of EOFs 1 and 2, with significant contribution from Neptune's equatorial latitudes, but low contributions towards the south pole. The spectral shape of EOF 3 looks similar to the expected spectral signature of changing the CH₄ scaling factor, shown earlier in Fig. 10 and reproduced in Fig. 13 for reference. Hence, this EOF seems to correspond rather well with the retrieved cloud-top methane abundance and we can see low values polewards of 40°S, just as we see in our retrieved methane profiles (Fig. 9) and faintly in our reflectivity difference maps (Figs. 11 and 12). In Fig. 9, we compare our retrieved methane abundances with the variation of $\alpha_3$ along the central meridian, scaled to match the overall variations, and find that it broadly matches the retrieved reduction in cloud-top methane abundance towards the south pole. We
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Averaged and contrast-enhanced reflectivity observations in the H<sub>2</sub>-H<sub>2</sub>/H<sub>2</sub>-He CIA band at 827 nm (left-hand panel) and methane-band at 833 nm (middle panel) and the difference between these averages (i.e. R<sub>827</sub> − R<sub>833</sub>) (right-hand panel), all averaging over observations ‘3’, ‘4’, and ‘6’ (observation ‘5’ was omitted as the striping artefacts were too noticeable). The south pole of Neptune is at bottom left, and latitude circles with spacing of 20° are indicated. It can be seen that the reflectivity difference image (right hand panel) is slightly darker at the bottom left than top right, indicated a lower methane spectral signal near the south pole. We can also see a slightly brighter region circling the south pole from 60 to 70°S, indicating slightly higher methane abundance there, which is also just seen in our retrieved latitudinal dependence (Fig. 9). This feature is just apparent in the H<sub>2</sub>-H<sub>2</sub>/H<sub>2</sub>-He image (left-hand panel), but not the methane-band image (middle panel), indicating that at 60–70°S we seem to have a slight thickening or vertical extension of the main cloud, which is masked in the methane image by an accompanying increase in methane abundance and thus absorption.</figcaption></figure>

can also see from both Figs. 9 and 13 that the contribution from EOF 3
has slightly lower values near the equator, similar to the very slight
reduction in retrieved methane abundance near the equator seen in
Fig. 9. Hence, applying the PCA technique in this spectral range seems
to recover one component that roughly matches the retrieved latitu-
dinal variation of methane cloud-top abundance and which adds
credibility to our conclusion that methane at Neptune's cloud tops is
diminished at latitudes polewards of 40°S at all longitudes and not just
along the central meridian as we have formally retrieved.

**6. Conclusion**

We have made an initial estimation of the latitudinal variation of
cloud-top methane abundance in Neptune's atmosphere, for the first
time from ground-based observations, using the new Narrow Field
Mode of the MUSE instrument at VLT. We find that methane varies with
latitude broadly similarly to the variation seen in HST/STIS observa-
tions in 2003 (Karkoschka and Tomasko, 2011). Hence, this distribution
with latitude would appear to be a long-lived feature (or is at least
slowly-varying) and our retrieved abundances appear to be consistent
with this earlier determination, with a retrieved mole fraction of 4–5 %
at equatorial latitudes, reducing to 3–4 % at polar latitudes. Our maps
of CH₄ abundance variations show considerable fine detail, which could
be useful for constraining the detailed structure of Neptune's atmo-
spheric dynamics. The lower values of cloud-top CH₄ over Neptune's
south pole is consistent both with the HST/STIS determinations of
Karkoschka and Tomasko (2011) and also with a multi-wavelength
analysis of de Pater et al. (2014), which concludes that the increased
brightness of Neptune's south pole at radio wavelengths is caused by a
lower humidity of hydrogen sulphide (H₂S). This low humidity may be
part of a wider overturning circulation in Neptune's atmosphere, which
has air rising at mid-latitudes (detected by cooler tropopause tem-
peratures observed in the mid-IR, e.g. Fletcher et al., 2014) and leads to
subsiding motion of ‘dry’ air at polar latitudes, extending all the way
from the stratosphere down to the cloud tops. In our observations, this
upwelling branch is coincident with the cloudy regions seen at 20–40°,
but what happens at equatorial latitudes is less clear. However, a

slightly lower methane abundance apparent at the equator in our re-
trieved parameters (Fig. 9) and also our EOF 3 contribution map
(Fig. 13), might indicate local downwelling at the equator also, which
may explain why the main cloud here is retrieved to have lower opacity
and is less vertically extended. Since we retrieve similar cloud-top
pressures at all latitudes, but very different methane abundances, we
surmise that this cloud is not composed of methane ice, but is most
likely composed of hydrogen sulphide ice, as is probably also the case in
Uranus's atmosphere (Irwin et al., 2018, 2019). Interestingly, de Pater
et al. (2014) note that South Polar Features (SPFs) seen in Neptune's
atmosphere at 60–70°, are convective storms, produced by baroclinic
instabilities at the edge of the south polar prograde jet. The fact that we
detect very slightly enhanced methane at these latitudes appears to
support this conclusion. Further analysis, combining these VLT/MUSE
observations with existing longer-wavelength observations from in-
struments such as VLT/SINFONI and Gemini/MUSE, is planned to more
tightly constrain the cloud properties (using a wider wavelength range)
and also explore alternative methods of parameterising the methane
vertical profile such as the “proportionally descended gas” profiles,
found by Sromovsky et al. (2011) to provide a good match to the HST/
STIS 2002 observations of Uranus. However, this extended analysis is
beyond the scope of this ‘initial-results’ paper, reporting the new
MUSE/NFM Neptune observations.

Facilities: VLT (MUSE).

**Acknowledgments**

We are grateful to the United Kingdom Science and Technology Facilities Council for funding this research. We thank Larry Sromovsky for providing the code used to generate our Rayleigh-scattering opacities. Glenn Orton was supported by NASA funding to the Jet Propulsion Laboratory, California Institute of Technology. Leigh Fletcher was supported by a Royal Society Research Fellowship at the University of Leicester. PMW received support from BMBF Verbundforschung (project MUSE-NFM, grant 05A17BAA). The observations reported in this paper have the ESO ID: 60.A-9100(K).
<figure><img src="image_17.png" /><figcaption>Fig. 13. Principal component analysis applied to MUSE observation ‘3’, showing the contribution of the first five Empirical Orthogonal Functions (EOF). The images in the left-hand column show the spatial variation of contribution of each EOF to the reconstructed spectra, while the right hand column shows the associated EOF itself (i.e. spectrum). In the images, Neptune’s south pole is again at bottom left. It can be seen that the eigenvalues of the EOFs fall rapidly with EOF number, and that most meaningful spatial variation in the observation is encapsulated in the first three EOFs. The shape of EOF 1 is approximately flat and this eigenfunction seems to encapsulate the overall reflectivity as can be seen in the associated image. EOF 2 appears to be almost the inverse of EOF 1, both in terms of its spatial distribution and also spectral shape, although its spatial structure exhibits finer structure, including a dark ‘halo’ about the south pole. The spatial distribution of EOF 3 has low values near Neptune’s south pole, but high values at lower latitudes, and its spectral shape is rather distinct. The spectral shape of EOF 3 is compared here with the computed change in modelled spectrum when the abundance of CH₄ is increased (green), taken from Fig. 10, where we can see a rough correspondence. Hence, we believe EOF 3 can be taken, to a first approximation, as a proxy for the abundance of CH₄ immediately above Neptune’s main cloud top at ~ 2 bar. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
## Appendix A. Calculation of Rayleigh-scattering optical depth

The Rayleigh-scattering cross-section of Neptune's air (i.e. $cm^2$) in our radiative transfer model is calculated as (e.g. Allen, 1976):

$$ \sigma(\lambda) = \frac{8\pi^3 \sum_i x_i (n_i(\lambda)^2 - 1)^2 \delta_i}{3\lambda^4 N^2 \sum_i x_i} $$

where the summation is over each gas present, with mole fraction $x_i$, refractive index spectrum $n_i(\lambda)$ (at standard temperature and pressure) and depolarization factor $\delta_i$. N is the total number density of air molecules (i.e. $cm^{-3}$). The depolarization factor for each gas is calculated as:

$$ \delta_i = \frac{6 + 3\Delta_i}{6 - 7\Delta_i} $$

where the factors $\Delta_i$ are listed in Table 3. Similarly, the refractive index of each gaseous component (at STP) is taken to vary with wavelength as

$$ n_i(\lambda) = 1 + A_i(1 + B_i/\lambda^2) $$

The approach of summing over mole fraction like this was suggested by Larry Sromovsky [private communication], who also shared the parameters used. The coefficients $A_i$ and $B_i$ are also listed in Table 3. Here, the refractive index parameters are taken from Allen (1976), although for CH$_4$, which is not listed, they have been set to be the same as NH$_3$. The depolarization coefficients, $\Delta_i$, are from Penndorf (1957), after Parthasarathy (1951); CO$_2$ values have been assumed for CH$_4$ and NH$_3$.

<table><caption>Table 3
Assumed Rayleigh-scattering parameters.</caption><thead><tr><td>Gas</td><td>A<sub>i</sub></td><td>B<sub>i</sub></td><td>Δ<sub>i</sub></td></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>13.58 × 10<sup>-5</sup></td><td>7.52 × 10<sup>-3</sup></td><td>0.0221</td></tr><tr><td>He</td><td>3.48 × 10<sup>-5</sup></td><td>2.3 × 10<sup>-3</sup></td><td>0.025</td></tr><tr><td>CH<sub>4</sub></td><td>37 × 10<sup>-5</sup></td><td>12 × 10<sup>-3</sup></td><td>0.0922</td></tr><tr><td>NH<sub>3</sub></td><td>37 × 10<sup>-5</sup></td><td>12 × 10<sup>-3</sup></td><td>0.0922</td></tr></tbody></table>

Air refractive index calculated as $n_i = 1 + A_i(1 + B_i/\lambda^2)$, where $\lambda$ is the wavelength in µm. $\Delta_i$ terms are the assumed depolarization factors. The literature sources for the coefficients presented here are described in the Appendix.

## References

* Allen, C.W., 1976. Astrophysical Quantities. Athlone Press, University of London.

* Arsenaull, R., et al., 2012. ESO adaptive optics facility progress report. Proc. SPIE. Adaptive Optics Systems III. 8447 https://doi.org/10.1117/12.926074.

* Bacon, R., et al., 2010. The MUSE second-generation VLT instrument. Proc. SPIE 7735, 773508. https://doi.org/10.1117/12.856027. Ground-based and Airborne Instrumentation for Astronomy III.

* Baines, K.H., Smith, W.H., 1990. The atmospheric structure and dynamical properties of Neptune derived from ground-based and IUE spectrophotometry. Icarus 85, 65–108.

* Borysow, A., Frommhold, L., Moraldi, M., 1989a. Collision-induced infrared spectra of H$_2$–He pairs involving 0→1 vibrational transitions and temperatures from 18 to 7000 K. ApJ 336, 495–503.

* Borysow, A., Frommhold, L., 1989b. Collision-induced infrared spectra of H$_2$–He pairs at temperatures from 18 to 7000 K. Overtone and hot bands. ApJ 341, 549–555.

* Borysow, A., Borysow, J., Fu, Y., 2000. Semi-empirical model of collision-induced absorption spectra of H$_2$–He complexes in the second overtone band of Hydrogen at temperatures from 50 to 500 K. Icarus 145, 601–608.

* Chance, K., Kurucz, R.L., 2010. An improved high-resolution solar reference spectrum for Earth's atmosphere measurements in the ultraviolet, visible, and near infrared. J. Quant. Spec. Radiat. Transf. 111, 1289–1295.

* de Pater, I., Fletcher, L.N., Luszcz-Cook, S., DeBoer, D., Butler, B., Hammel, H.B., Sitko, M.L., Orton, G., Marcus, P.S., 2014. Neptune's global circulation deduced from multi-wavelength observations. Icarus 237, 211–238.

* Dyudina, U.A., Ingersoll, A.P., Danielson, G.E., Baines, K.H., Carlson, R.W., The Galileo NIMS, SSI teams, 2001. Interpretation of NIMS and SSI images on the Jovian cloud structure. Icarus 150, 219–233.

* Fletcher, L.N., de Pater, I., Orton, G.S., Hammel, H.B., Sitko, M.L., Irwin, P.G.J., 2014. Neptune at summer solstice: zonal mean temperatures from ground-based observations, 2003–2007. Icarus 231, 146–167.

* Gautier, D., Conrath, B.J., Owen, T., de Pater, I., Atreya, S.K., 1995. The troposphere of Neptune. In: Cruikshank, D.P. (Ed.), In Neptune and Triton. University of Arizona Press, pp. 547–613.

* Havriliak, S., Swenson, R.W., Cole, R.H., 1955. Dielectric constants of liquid and solid hydrogen sulfide. J. Chem. Phys. 23, 134–135.

* Irwin, P.G.J., de Bergh, C., Courtin, R., Bézard, B., Teanby, N.A., Davis, G.R., Fletcher, L.N., Orton, G.S., Calcutt, S.B., Tice, D., Hurley, J., 2012. The application of new methane line absorption data to Gemini-N/NIFS and KPNO/FTS observations of Uranus' near-infrared spectrum. Icarus 220, 369–382.

* Irwin, P.G.J., Dyudina, U., 2002. The retrieval of cloud structure maps in the equatorial region of Jupiter using a principal component analysis of Galileo/NIMS. Icarus 156, 52–63.

* Irwin, P.G.J., Fletcher, L.N., Tice, D., Owen, S.J., Orton, G.S., Teanby, N.A., Davis, G.R., 2016. Time variability of Neptune's horizontal and vertical cloud structure revealed by VLT/SINFONI and Gemini/NIFS from 2009 to 2013. Icarus 271, 418–437.

* Irwin, P.G.J., Lellouch, E., de Bergh, C., Courtin, R., Bézard, B., Fletcher, L.N., Orton, G.S.,

* Teanby, N.A., Calcutt, S.B., Tice, D., Hurley, J., Davis, G.R., 2014. Line-by-line analysis of Neptune's near-IR spectrum observed with Gemini/NIFS and VLT/CRIRES. Icarus 227, 37–48.

* Irwin, P.G.J., Teanby, N.A., Davis, G.R., Fletcher, L.N., Orton, G.S., Tice, D., Hurley, J., Calcutt, S.B., 2011. Multispectral imaging observations of Neptune's cloud structure with Gemini-North. Icarus 216, 141–158.

* Irwin, P.G.J., Teanby, N.A., de Kok, R., Fletcher, L.N., Howett, C.J.A., Tsang, C.C.C., Wilson, C.F., Calcutt, S.B., Nixon, C.A., Parrish, P.D., 2008. The NEMESIS planetary atmosphere radiative transfer and retrieval tool. J. Quant. Spec. Radiat. Transf. 109, 1136–1150.

* Irwin, P.G.J., Tice, D.S., Fletcher, L.N., Barstow, J.K., Teanby, N.A., Orton, G.S., Davis, G.R., 2015. Reanalysis of Uranus' cloud scattering properties from IRTF/spex observations using a self-consistent scattering cloud retrieval scheme. Icarus 250, 462–476.

* Irwin, P.G.J., Toledo, D., Garland, R., Teanby, N.A., Fletcher, L.N., Orton, G.S., Bézard, B., 2018. Detection of hydrogen sulfide above the clouds in Uranus's atmosphere. Nat. Astron. 2, 420–427.

* Irwin, P.G.J., Toledo, D., Garland, R., Teanby, N.A., Fletcher, L.N., Orton, G.S., Bézard, B., 2019. Probable detection of hydrogen sulphide (H$_2$S) in Neptune's atmosphere. Icarus 321, 550–563.

* Karkoschka, E., 1994. Spectrophotometry of the Jovian planets and Titan at 300 to 1000 nm wavelength: the methane spectrum. Icarus 111, 174–192.

* Karkoschka, E., Tomasko, M., 2009. The haze and methane distributions on Uranus from HST-STIS spectroscopy. Icarus 202, 287–309.

* Karkoschka, E., Tomasko, M., 2010. Methane absorption coefficients for the Jovian planets from laboratory, Huygens, and HST data. Icarus 205, 674–694.

* Karkoschka, E., Tomasko, M.G., 2011. The haze and methane distributions on Neptune from HST-STIS spectroscopy. Icarus 211, 780–797.

* Lacis, A.A., Oinas, V., 1991. A description of the correlated-k distribution method for modelling nongray gaseous absorption, thermal emission, and multiple scattering in vertically inhomogeneous atmospheres. J. Geophys. Res. 96, 9027–9063.

* Lindal, G.F., 1992. The atmosphere of Neptune: an analysis of radio occultation data acquired by Voyager 2. Astron. J. 103, 967–982.

* Lellouch, et al., 2010. First results of Herschel-PACS observations of Neptune. A&A 518, L152.

* Martonchik, J.V., Orton, G.S., 1994. Optical constants of liquid and solid methane. Appl. Opt. 33, 8306–8313.

* Murtagh, F., Heck, A., 1987. Multivariate Data Analysis. Reidel, Dordrecht/Boston/Lancaster/Tokyo.

* Parthasarathy, S., 1951. Indian J. Phys. 25, 21.

* Penndorf, R., 1957. Tables of the refractive index for standard air and the Rayleigh scattering coefficient for the spectral region between 0.2 and 20.0 µm and their application to atmospheric optics. J. Opt. Soc. Am. 47, 176–182.

* Plass, G.N., Kattawar, G.W., Catchings, F.E., 1973. Matrix operator method of radiative transfer. 1: Rayleigh scattering. Appl. Opt. 12, 314–329.

* Sheik-Bahae, M., 2005. Nonlinear optics basics. Kramers-Kronig relations in nonlinear optics. In: Guenther, R.D. (Ed.), Encyclopedia of Modern Optics. Academic Press,
Amsterdam.

Simon, A.A., Wong, M.H., Orton, G.S., 2015. First results from the Hubble OPAL program: Jupiter in 2015. ApJ 812, 55.

Sromovsky, L.A., 2005. Accurate and approximate calculations of Raman scattering in the atmosphere of Neptune. Icarus 173, 254–283.

Sromovsky, L.A., Fry, P.M., Kim, J.H., 2011. Methane on Uranus: the case for a compact CH₄ cloud layer at low latitudes and a severe CH₄ depletion at high latitudes based on re-analysis of Voyager occultation measurements and STIS spectroscopy. Icarus 215, 292–312.

Sromovsky, L.A., Karkoschka, E., Fry, P.M., de Pater, I., Hammel, H.B., 2019. The methane distribution and polar brightening on Uranus based on HST/TSTIS, Keck/NIRC2, and IRTF/SPEX observations through 2015. Icarus 317, 266–306.

Sromovsky, L.A., Karkoschka, E., Fry, P.M., Hammel, H.B., de Pater, I., Rages, K., 2014. Methane depletion in both polar regions of Uranus inferred from HST/STIS and KECK/NIRC2 observations. Icarus 238, 137–155.

Ströbele, S., La Penna, P., Arsenault, R., Conzelmann, R.D., Delabre, B., Duchateau, M.,

Dorn, R., Fedrigo, E., Hubin, N., Quentin, J., et al., 2012. GALACSI system design and analysis. Proc. SPIE. Adaptive Optics Systems III. 8447https://doi.org/10.1117/12.926110.

Tice, D.S., Irwin, P.G.J., Fletcher, L.N., Teanby, N.A., Hurley, J., Orton, G.S., Davis, G.R., 2013. Uranus' cloud particle properties and latitudinal methane variation from IRTF SpeX observations. Icarus 223, 684–698.

Tollefson, J., de Pater, I., Marcus, P.S., Luszcz-Cook, S., Sromovsky, L.A., Fry, P.M., Fletcher, L.N., Wong, M.H., 2018. Vertical wind shear in Neptune's upper atmosphere explained with a modified thermal wind equation. Icarus 311, 317–339.

Weilbacher, P.M., Streicher, O., Urrutia, T., Pécontal-Rousset, A., Jarno, A., Bacon, R., 2014. The MUSE data reduction pipeline: status after preliminary acceptance europe. ASP conference series In: Manset, N., Forshay, P. (Eds.), Astronomical Data Analysis Software and Systems XXIII. Proceedings of a meeting held 29 September–3 October 2013 at Waikoloa Beach Marriott, Hawaii, USA. 485. pp. p.451.

Wong, L.Y., Anderson, A., 1974. Refractive-index measurements of some low-temperature liquids. J. Opt. Soc. Am. 62, 219–222.