Latitudinal Variations in Methane Abundance, Aerosol Opacity and Aerosol Scattering Efficiency in Neptune's Atmosphere Determined From VLT/MUSE

Spectral observations of Neptune made in 2019 with the Multi Unit Spectroscopic Explorer (MUSE) instrument at the Very Large Telescope (VLT) in Chile have been analyzed to determine the spatial variation of aerosol scattering properties and methane abundance in Neptune's atmosphere. The darkening of the South Polar Wave at ∼60°S, and dark spots such as the Voyager 2 Great Dark Spot is concluded to be due to a spectrally dependent darkening (λ < 650 nm) of particles in a deep aerosol layer at ∼5 bar and presumed to be composed of a mixture of photochemically generated haze and H2S ice. We also note a regular latitudinal variation of reflectivity at wavelengths of very low methane absorption longer than ∼650 nm, with bright zones latitudinally separated by ∼25°. This feature, which has similar spectral characteristics to a discrete deep bright spot DBS‐2019 found in our data, is found to be consistent with a brightening of the particles in the same ∼5‐bar aerosol layer at λ > 650 nm. We find the properties of an overlying methane/haze aerosol layer at ∼2 bar are, to first‐order, invariant with latitude, while variations in the opacity of an upper tropospheric haze layer reproduce the observed reflectivity at methane‐absorbing wavelengths, with higher abundances found at the equator and also in a narrow “zone” at 80°S. Finally, we find the mean abundance of methane below its condensation level to be 6%–7% at the equator reducing to ∼3% south of ∼25°S, although the absolute abundances are model dependent.


Introduction
The study of the aerosols in Neptune's atmosphere has been transformed in recent years by the analysis of visible and near-infrared multi-spectral imaging data.The retrieval of vertical profiles of aerosols from such data requires knowledge of the vertical and latitudinal distribution of the main gaseous absorber, methane, which for many years, in the absence of any other information, was assumed to be constant at all latitudes.However, by analysing Hubble Space Telescope (HST) Space Telescope Imaging Spectrograph (STIS) observations of Neptune from 2003, Karkoschka and Tomasko (2011) showed that the abundance of methane at equatorial latitudes (∼ 4%) was approximately twice that detected at polar latitudes (∼ 2%), a result that had significant repercussions on the inferred vertical structure of aerosols, which were subsequently found to vary less significantly with latitude.In 2018, Neptune was observed during commissioning operations with the Multi Unit Spectroscopic Explorer (MUSE) Integral Field Unit (IFU) Spectrometer at the Very Large Telescope (VLT) at the European Southern Observatory in Chile.An initial analysis of these data (Irwin, Toledo, Braude, et al., 2019) found a similar lat-itudinal variation of cloud-top methane mole fraction as the HST/STIS study, and was later revised to include modelling of the limb-darkening effects (Irwin et al., 2021) using the Minnaert approximation (Minnaert, 1941).Irwin et al. (2021) concluded that the 'deep' (i.e., at 2-4 bar) mole fraction of methane varied from 4-6% at the equator to 2-4% at polar latitudes, with a boundary at ∼30 • S. Most recently, a joint analysis of HST/STIS, Gemini/NIFS and IRTF/SpeX observations of Neptune and Uranus from 0.3 to 2.5 µm (Irwin, Teanby, Fletcher, et al., 2022) has developed a 'holistic' aerosol model of both planets comprised of three basic distinct aerosol layers: 1) a deep H 2 S/photochemicalhaze aerosol layer with a base pressure ≥ 5-7 bar (Aerosol-1); 2) a layer of methane/photochemicalhaze just above the methane condensation level at 1-2 bar (Aerosol-2); and 3) an extended layer of small photochemical haze particles extending into the stratosphere (Aerosol-3).For Neptune an additional contribution from upper level (∼ 0.2 bar) methane ice clouds was required to match the IRTF/SpeX observations at λ > 1 µm.
The atmosphere of Neptune occasionally displays dark spots, seen at blue-green visible wavelengths < 650 nm, including the Great Dark Spot (GDS) observed by Voyager 2 in 1989 (Smith et al., 1989), and more recent examples captured in HST Wide Field Camera 3 (WFC3) observations (e.g., Hueso & Sánchez-Lavega, 2019).The most recent dark spot was discovered by HST/WFC3 in 2018 at 23 • N and named NDS-2018 (Simon et al., 2019).NDS-2018 was of a similar size to the GDS and subsequently drifted equatorwards, disappearing in 2022 (Wong et al., 2022).As part of a global effort to observe and analyse NDS-2018, Neptune was observed with VLT/MUSE in late 2019 (Irwin et al., 2023).After spatial deconvolution, Irwin et al. (2023) were able to detect the NDS-2018 spot at ∼15 • N in these observations (the first ground-based detection of such a spot), and were also able to spectrally characterise it at high spectral resolution, the first time that this has ever been achieved for a Neptunian dark spot.Irwin et al. (2023) found that NDS-2018 was caused by a spectrally-dependent (λ < 650 nm) darkening of the particles in the Aerosol-1 layer at ∼5 bar and also found a nearby deep bright spot, DBS-2019 at ∼ 10 • N, which they concluded was caused by a brightening of the same layer at longer wavelengths.
In addition to capturing the reflection spectra of NDS-2018 and DBS-2019, the 2019 VLT/MUSE data also show distinct latitudinal variations over a wide range of wavelengths (473 -933 nm), in particular the South Polar Wave (SPW) at ∼60 • S, first seen in Voyager 2 images in 1989 (Smith et al., 1989), which is dark at blue-green wavelengths, but invisible at longer wavelengths.The spectral features of the SPW are very similar to those of dark spots and Karkoschka (2011a) concluded that it is caused by a darkening of particles at pressures > 3 bar, which was later confirmed by Irwin, Teanby, Fletcher, et al. (2022).The SPW has been visible ever since the Voyager 2 flyby in HST imaging observations and these observations are reviewed in detail by Karkoschka (2011b).The SPW is found to have a southern boundary at ∼ 70 • S, and a clearer northern boundary that varies with longitude as a wavenumber-1 disturbance of amplitude ∼ 5 • , and whose northern extent varies with time from 50 -55 • S. The SPW feature seems to be related to the South Polar Feature (SPF) near 70 • S, which are short-lived, bright clouds that change on a time scale of hours (Hammel et al., 1989).Karkoschka (2011b) notes that the SPF longitude coincides with that of the northernmost extent of the SPW, which suggests there is a dynamical link.The features may also have been dynamically linked with the second Voyager dark spot (DS2) (Sromovsky et al., 1993) and Karkoschka (2011b) suggest that these features are static with respect to a new coordinate system for the deep rotation of Neptune with a rotation rate of 15.9663 hours, compared with 16.108 hours derived from Voyager radio data (Lecacheux et al., 1993).
While the 2019 VLT/MUSE data show the SPW very clearly, they also detected latitudinal variations at several other wavelengths that have not been fully noted before.In this paper, we further analyse the 2019 VLT/MUSE data to revise our understand-ing of the latitudinal variation in aerosol properties and methane abundance in Neptune's atmosphere.

Observations
The observations of Neptune were recorded in October and November 2019 using the Multi Unit Spectroscopic Explorer (MUSE) Integral Field Unit (IFU) spectrometer (Bacon et al., 2010) at the Very Large Telescope (VLT) of the European Southern Observatory at La Paranal, Chile.MUSE records 'cubes' of data, where each point in the 300 × 300 'spaxel' image contains a complete spectrum covering ∼ 3700 wavelengths from 473 to 933 nm at a spectral resolution of 2000 -4000.The observations were recorded in Narrow-Field Mode, which has field of view of 7.5" × 7.5" and each 'spaxel' is of size 0.025" × 0.025" (equivalent to 530 km × 530 km on Neptune's disc).To improve the signal-to-noise ratio, and also to make the observations consistent with our only available source of methane absorption data (Karkoschka & Tomasko, 2010), the spectra were averaged with a triangular instrument lineshape of Full-Width-Half-Maximum (FWHM) 2 nm, sampled every 1 nm (to achieve Nyquist sampling), reducing the number of wavelengths from ∼ 3700 to 459.Neptune was observed on several occasions in this programme (summarised in Table 1), where five exposures on October 18th (Observations 6 -10), included the NDS-2018 dark spot.Although these data were recorded with the GALACSI adaptive optics system (Stuik et al., 2006), using a laser guide star, the achieved spatial resolution near 500 nm was insufficient to resolve the faint NDS-2018 feature, which has a diameter of ∼0.2" and has low contrast.However, with the development of a novel deconvolution technique, modified-clean, the deconvolved 'cubes' had sufficient discrimination to detect and spectrally characterise the NDS-2018 dark spot and also a nearby deep bright spot DBS-2019 (Irwin et al., 2023).The appearance of Neptune at several wavelengths in the best-resolved deconvolved cube in this set, 'Obs-6' is shown in Fig. 1.Irwin et al. (2023) found that the NDS-2018 dark spot (visible at 551 nm in Fig. 1) is formed by a spectrally-dependent (λ < 650 nm) darkening of the aerosols in the deep Aerosol-1 layer at ∼5 bar, presumably coincident with the H 2 S condensation level.Irwin et al. (2023) also found a new deep bright spot near NDS-2018, visible in the 831-nm image, which they named 'Deep Bright Spot -2019' (DBS-2019), and which, from the narrowness of its reflectance peaks, was determined to be caused by a spectrally-dependent brightening of the same ∼5-bar Aerosol-1 layer at wavelengths longer than ∼ 650 nm.

NDS-2018
DBS-2019 SPW While Irwin et al. (2023) describe the characterisation and interpretation of the NDS-2018 and DBS-2019 discrete features, the MUSE data also have sufficient spatial and spec-tral resolution to determine latitudinal variations of aerosol opacity and methane abundance.These latitudinal changes are clearly visible in Fig. 1, including the dark SPW at 60 • S at short wavelengths (here at 551 nm), but also prominent banded features at longer continuum wavelengths (here at 831 nm), and a small bright collar (80 • S) seen about the south pole at 860 nm, which at a wavelength of strong methane absorption must be caused by increased aerosol abundance high in the atmosphere.The banding seen at 831 nm is also just visible in recent 845-nm HST/WFC3 images (e.g., Chavez et al., 2023), which covers the same reflectance peak, but at much lower spectral resolution (∆λ ∼ 84 nm).However, the banding is particularly prominent at this wavelength, which is at the centre of reflectance peak where we can see to deep pressures in the atmosphere, and the narrowness of the spectral signature of these bands, discussed later, indicate that they are likely caused by variations of the deep aerosol layers.
Although the MUSE observations were photometrically corrected by observing a standard star shortly before viewing Neptune, uncertainty remained in the absolute correction.Hence, to ensure the disc-averaged MUSE spectrum was consistent with the discaveraged spectrum of Neptune measured in 2003 with HST/STIS (Karkoschka & Tomasko, 2011), the data were scaled to give the same reflectivity as HST/STIS in the equatorial region of 10 • S -10 • N. Long-term records of the disc-integrated blue and green magnitudes of Neptune (Lockwood, 2019) reveal no notable changes between 2003 and the 2016, which justifies this simple scaling and which also has the merit of allowing direct comparison with retrievals from the HST/STIS observations.

Analysis
The analysis in this paper follows on from a combined analysis of HST/STIS, IRTF/SpeX, and Gemini/NIFS observations of both Uranus and Neptune by Irwin, Teanby, Fletcher, et al. (2022), which we will refer to henceforth as 'IRW22' or the 'holistic' model.Given that we do not know the composition of the aerosols in Neptune's atmosphere, and do not have a definite expectation of the volume mixing ratio profile of the main visible/IR absorber, methane, the simultaneous retrieval of aerosol and methane abundance profiles is a degenerate problem, even when analysing the 800 -860 nm region that can differentiate between aerosols and methane abundance (Karkoschka & Tomasko, 2011).As a result, since previous studies have concentrated on particular narrow wavelength ranges, this has led to a multitude of similar, but not wholly consistent aerosol/methane solutions.One way to reduce the degeneracy is to analyse simultaneously as wide a wavelength range as possible in order that the sizes of the scattering particles can be better constrained (small particles will have cross-sections that fall as 1/λ 4 , while larger particles will have cross-sections that vary more slowly with wavelength).The wavelength range analysed by IRW22 was 0.3 -2.4 µm, which allowed good discrimination against particle size.However, even when analysing a wide range of wavelengths, there are still multiple solutions as we do not know the composition of particles and so we must try to infer the scattering properties (in this case the opacities and single-scattering albedos) directly from the data.IRW22 addressed this degeneracy by attempting to fit simultaneously not only the disc-averaged reflectance spectrum over this range, but also how the spectrum varied with solar and viewing zenith angles, through fitting their 'limbdarkening' or 'centre-to-limb' functions.This was done by analysing all observations on the planets' discs and quantifying how the reflectance at each wavelength, (I/F ), varied towards the limb of the planet.IRW22 found that this variation was well approximated by the Minnaert approximation (Minnaert, 1941): Here, µ and µ 0 are the cosines of the viewing and solar zenith angles, respectively, (I/F ) 0 is the fitted nadir reflectance, and k the fitted limb-darkening parameter.The reflectance, I/F , is the observed radiance at each location, I, divided by the radiance reflected from a perfectly-reflecting Lambertian surface at the same point, and in this study we used the solar spectrum of Chance and Kurucz (2010).IRW22 estimated the scattering properties of the particles in the different aerosol layers using Mie theory and fitted for their complex refractive index spectra.Note that this analysis assumes that the particles in a certain layer, which might be a mixture of ice and haze particles of various size distributions, can be approximated as being composed of a single size distribution of particles, with a single complex refractive index spectrum.The imaginary refractive index spectra were treated as free parameters to be fitted, while the real parts were computed using a Kramers-Kronig analysis, assuming n real = 1.4 at 800 nm; this value of n real is typical of giant planet condensates, such as methane ice, but is not well constrained (other values were tested and gave similar overall results).We note that the Kramers-Kronig approach requires that we know n imag (λ) at all wavelengths to accurately reproduce n real (λ), whereas in this study n imag (λ) could only be estimated from 0.3 to 2.4 µm.Hence, this reconstruction of n real (λ) is only an approximation.Mie theory was then used to compute the extinction cross-section and single-scattering albedo of the aerosols as a function of wavelength.Finally, the Mie-calculated phase functions were approximated with combined Henyey-Greenstein phase functions at each wavelength to average over features particular to spherical particles, namely the 'rainbow' peak and the back-scattered 'glory'.Although containing several approximations, this approach, first presented by Irwin et al. (2015), has the merit of generating self-consistent values of the extinction cross-section, single-scattering albedo, and phase function spectra, and has considerable advantages, we believe, over approaches where the single-scattering albedo (or phase functions, or extinction coefficients) are modified directly and separately.The holistic IRW22 analysis found that the combined HST/STIS, IRTF/SpeX and Gemini/NIFS observations of both Uranus and Neptune were well approximated with an atmospheric model consisting of three main aerosol layers, outlined earlier: 1) 'Aerosol-1', a deep aerosol layer with a base pressure ≥ 5-7 bar, assumed to be composed of a mixture of H 2 S ice and photochemical haze of mean radius 0.05 µm; 2) 'Aerosol-2', a layer of photochemical haze and methane ice of mean radius ∼ 0.5 µm, confined within a layer of high static stability at the methane condensation level at 1-2 bar; and 3) 'Aerosol-3', an extended layer of small photochemical haze particles (r ∼ 0.05 µm), probably of the same composition as the haze in the 1-2-bar layer, extending from this level up through to the stratosphere.An additional thin layer of micron-sized (r ∼ 3 µm) methane ice particles at ∼0.2 bar ('Aerosol-4') was required to explain enhanced reflection at methane-absorbing wavelengths in the Neptune data longer than 1.0 µm, caused probably by discrete clouds at this level being present along the central meridian of the IRTF/SpeX line-averaged data used, and unresolved clouds in the Gemini/NIFS observations.The region of high static stability at the methane condensation level has been noted by several previous authors (e.g., Hueso et al., 2020;Leconte et al., 2017).
We studied the VLT/MUSE data in the same way as IRW22 and concentrated on the best-resolved, deconvolved 'Obs-6' cube, described by Irwin et al. (2023).For this study we used the same simple 'step' model for the methane profile as IRW22, with variable deep abundance, a fixed relative humidity above the condensation level (discussed later) and limiting the mole fraction in the stratosphere to 1.5×10 −3 (Lellouch et al., 2010).Sample profiles using this model are shown in Fig. 2. We first analysed the discaverage of the deconvolved 'Obs-6' cube, and then analysed the data in latitude bands of width 10 • , spaced every 5 • to achieve Nyquist sampling (masking out discrete clouds in both cases).The high spatial resolution of these observations, combined with a remarkably low level of cloud activity at this time (e.g., Chavez et al., 2023), meant that the few high-altitude discrete methane ice clouds that were visible were easily maskedout in our data.This meant that we did not need to include the upper troposphere (∼ 0.2 bar) layer of micron-sized methane ice 'Aerosol-4' particles in our model, which only have a significant contribution at wavelengths greater than 1 µm anyway.Using a Minnaert analysis we extracted the disc-average nadir reflectance spectrum, (I/F ) 0 (λ), and limb-darkening spectrum, k(λ), from the observations, where λ is wavelength.These spectra were used to reconstruct two synthetic observation spectra, calculated for back-scattering conditions, with the zenith angles (both viewing and solar) set to either 0 • (i.e., nadir), or 61.45 • , which were then fitted simultaneously with our NEMESIS radiative transfer and retrieval tool (Irwin et al., 2008).We fitted to these two angles simultaneously in order to fit both the mean reflectance and limb-darkening spectra, and these two particular zenith angles were chosen to coincide with two of the zenith angles in our matrixoperator multiple-scattering scheme (Plass et al., 1973), thus avoiding any interpolation errors.In this analysis (following Irwin et al., 2023) we increased the mean radius of the Aerosol-1 particles to 0.1 µm (standard Gamma distribution of sizes with variance σ = 0.05), consistent with expectation of fog-like particles, and fixed the mean radius of the Aerosol-2 particles to 0.7 µm (σ = 0.3), which was found by IRW22 to be the mean size of the Aerosol-2 particles that was most consistent with their observations.The size of the Aerosol-3 particles was left unchanged, with a mean radius 0.05 µm and variance σ = 0.05.When fitting to the disc-averaged MUSE data we found that the estimated random errors of the reconstructed spectra were smaller than the uncertainties in our radiative transfer model arising from parameterisation choices and the assumed gaseous absorption coefficients of Karkoschka and Tomasko (2010).Hence, following IRW22, the errors on these spectra were set to 1/50 of the maximum disc-averaged nadir reflectance within 100 nm of each point in the spectrum.These errors were reduced by a factor of 2 in the 800 -900 nm range to ensure a good fit to the region best able to differentiate between cloud and methane abundance.In addition, the nadir reflectivity error was limited to not exceed 1% to ensure that NEMESIS fitted well to the wavelengths near the peak of reflection at ∼ 500 nm.With these errors (set to be the same for both the 0 • and 61.45 • spectra) we could fit the spectra to χ 2 /n ∼3, giving retrieved parameters with meaningful error values.For the shorter wavelengths we also had to account for a small amount of Raman-scattering (where photons absorbed at shorter wavelengths can be scattered to longer wavelengths in hydrogen-dominated atmospheres) and polarisation effects, which we did following the procedures of Sromovsky (2005a) and Sromovsky (2005b), respectively, as outlined by IRW22.To account for the Raman scattering, we added 100 points below the MUSE spectrum in each pixel, covering 373 -473 nm, which we reconstructed from a disc-averaged Minnaert analysis of the 2003 HST/STIS observations (Karkoschka & Tomasko, 2011), using the fitted Minnaert (I/F ) 0 and k spectra to simulate the expected spectrum for the observing geometry of each pixel and scaling this to match the MUSE spectrum where they overlap.In our radiative transfer model we simulated the entire spectrum from 373 to 933 nm, including Raman-scattering.Hence, photons Raman-scattered from wavelengths as short as 373 nm were approximated for in the modelled spectra.However, to avoid the retrieval model trying to fit these synthetic 373 -473 nm data, the error on these points was set to 100%.Finally, for consistency with IRW22 and ease of comparison with their HST/STIS retrievals, the MUSE data were analysed at the same set of wavelengths as for the STIS data.The assumed gaseous absorption coefficients and the assumed vertical profiles of temperature and gaseous abundance were also those used by IRW22 and Irwin et al. (2023).
During the initial retrievals from the latitudinally-averaged data, we found that there was considerable degeneracy between the methane abundance and the opacities of the Aerosol-1 and Aerosol-2 layers, which led to these parameters not varying smoothly with latitude, but showing considerable degeneracy, or 'cross-talk'.The reason for this can be understood from Fig. 2, which shows the best-fit cloud model of IRW22 determined from their combined Neptune data set.It can be seen that there is considerable overlap between the Aerosol-1 and Aerosol-2 profiles in the 2-3 bar region.Since the modelled spectra were sensitive to both clouds in this pressure region, the two opacities are able to interfere with each other with, for example, one reducing and the other increasing leading to little change in the overall modelled reflectance.Hence, we revised the parameterisation of the Aerosol-1 layer to be like the Aerosol-2 layer, i.e., a single vertically thin cloud, but centred at a higher, fixed pressure of 5 bar (Fig. 2).This param-Figure 2. Cloud opacity profiles retrieved in this analysis from disc-averaged observations, and assumed example methane abundance profiles.Left: Cloud scheme used by IRW22 (Irwin, Teanby, Fletcher, et al., 2022) in their combined analysis of HST/STIS, IRTF/SpeX and Gemini/NIFS.Middle: Revised cloud scheme used in this analysis of VLT/MUSE observations.The horizontal thickness of the lines used to plot the profiles indicates the formal retrieved error ranges for the cloud profiles.Right: Assumed methane mole fraction profile for different values of the deep abundance, where the relative humidity above the condensation level is fixed to 50% and the stratospheric mole fraction is limited to 1.5 × 10 −3 as described in the main text.
eterisation was found to fit the MUSE data just as well as the original IRW22 parameterisation, but led to much less overlap between the aerosol distributions and thus much less degeneracy between the Aerosol 1 and 2 opacities and methane abundances in the latitudinal retrievals reported later.This revised Aerosol-1 cloud model was also used by Irwin et al. (2023) in their analysis of the NDS-2018 and DBS-2019 features in this data set.We note that our data do not allow us to determine whether there is a clear gap between the Aerosol-1 and Aerosol-2 layers, only that this decoupling is necessary for our retrieval approach to be reliable and stable.

Disc-average retrievals and analysis
Using the deconvolved 'Obs-6' cube, we first determined the disc-averaged limbdarkening properties for each wavelength (masking out discrete cloud features and instrument artefacts) to yield spectra of the fitted parameters (I/F ) 0 (λ) and k(λ).These parameters were then used to compute the synthetic measurement spectra seen in Fig. 3, where the top panel compares two spectra reconstructed from the Minnaert coefficients with both the solar and zenith angles set to either θ = 0 • or θ = 61.45• , using (I/F )(θ) = (I/F ) 0 µ 2k−1 where µ 0 = µ 0 = cos θ.In this plot, the synthetic estimated error lim- its, discussed earlier, are shaded in grey for both spectra.Note that the spectrum reconstructed at θ = 0 • , is simply (I/F ) 0 (λ).Our NEMESIS fits to these spectra, using our modified aerosol model, are overplotted and show excellent agreement.
The bottom panel of Fig. 3 compares the limb-darkening coefficient spectrum of the disc-averaged MUSE data with that extracted from our NEMESIS fits to the reconstructed MUSE spectra at 0 • and 61.45 • , using showing that NEMESIS correctly fits both the mean reflectivity and limb-darkening spectra.
The spectra in both panels are compared with those obtained from an identical analysis by IRW22 of the HST/STIS observations of Neptune from 2003 (Karkoschka & Tomasko, 2011) and shows that the disc-average data are very similar, which is expected since the MUSE data were scaled to match HST/STIS.However, when comparing the limb-darkening spectra with HST/STIS it can be seen that the MUSE spectra show considerably more limb-brightening at methane-absorbing wavelengths (i.e., k(λ) is smaller than for STIS).This may possibly result from the different deconvolution schemes used in the STIS and MUSE datasets, but as we do not have the original undeconvolved STIS observations we are unable to test this.However, this difference in the limb-darkening meant that the The errors on the synthetic measured spectra are indicated in grey.The spectrum calculated at θ = θ0 = 0 • is of course simply (I/F )0.Bottom panel: MUSE disc-average limb-darkening spectrum k(λ) and fitted limb-darkening spectrum, calculated from the NEMESIS fits to the reconstructed spectra at 0 and 61.45 • zenith angle using Eq. 2. In both panels, the corresponding HST/STIS disc-average spectra from 2003, analysed by IRW22, are over-plotted.The wavelengths greyed-out from 585 -602 nm, were reserved for the laser guide star adaptive optics system of MUSE.In the bottom panel, the dashed line at k = 0.5 is added to help differentiate between limb-darkened (k > 0.5) and limb-brightened (k < 0.5) wavelengths.
Mie-scattering properties calculated from the best-fitting n imag spectra of the three aerosol layers derived by IRW22 did not achieve a very good fit to these MUSE data.Hence, the n imag spectra had to be refitted to the MUSE observations, which we did using our modified 'holistic' aerosol model, also fitting the opacities of the three aerosol layers, the pressure of the Aerosol-2 layer, and the deep methane abundance.The methane relative humidity is difficult to constrain unambiguously from the MUSE wavelength range.IRW22 retrieved a methane relative humidity of 35% from their data set, which included longer-wavelength IRTF/SpeX observations, while Karkoschka and Tomasko (2011) used a fixed value of 60% to match their HST/STIS data.Here, we fix the methane relative humidity above the condensation level to 50%.Our final retrieved disc-averaged n imag spectra for all three aerosols are compared with those derived by IRW22 for Neptune in Fig. 4.Here we can see that for VLT/MUSE the Aerosol-1 particles are determined to be more absorbing (i.e., n imag is higher), the Aerosol-2 particles have similar n imag spectra at short wavelengths, but have higher n imag at λ > 600 nm, while the Aerosol-3 particles (most visible at methane-absorbing wavelengths) are required to be considerably more scattering (i.e., n imag is lower) in order to fit the MUSE observations.(Irwin, Teanby, Fletcher, et al., 2022).

Latitudinally-resolved retrievals and analysis
Having re-fitted the n imag spectra (and thus all associated Mie-calculated scattering properties) of the three aerosol layers, we then determined latitudinal changes in opacity and methane abundance by analysing spectra reconstructed from the latitudinallyresolved Minnaert analysis of the MUSE dataset, where the data were averaged in latitudes bins of width 10 • , spaced every 5 • .A summary of the latitudinally-resolved Minnaert nadir reflectivity spectra, (I/F ) 0 (ϕ, λ), and limb-darkening spectra, k(ϕ, λ), where ϕ is the planetographic latitude, is shown in Fig. 5.For this latitudinal analysis there was not a sufficiently large range of zenith angles south of 70 • S and north of 20 • N to extract meaningful values of both (I/F ) 0 (ϕ, λ) and k(ϕ, λ).Hence for latitudes south of 70 • S, the limb-darkening spectrum, k(λ), was fixed to the average for all latitudes south of 70 • S and only (I/F ) 0 (ϕ, λ) fitted.For latitudes north of 20 • N, in the absence if any other information and since we believe the winds to be, to first order, hemispherically symmetric (Sromovsky et al., 1993), we assumed north-south symmetry in the cloud distribution also.Hence, the limb-darkening spectra were fixed to those found at the corresponding southern latitude, i.e., k(ϕ, λ) = k(−ϕ, λ), and again only (I/F ) 0 (ϕ, λ) fitted.In Fig. 5 we can see the signature of the South Polar Wave (SPW) at ∼ 60 • S in the (I/F ) 0 (λ) difference map at λ < 650 nm and also a slight increase in (I/F ) 0 (λ) south of ∼ 25 • S at most wavelengths.Most of the features in the k(λ) difference map are in regions of low k-values where the reflectances are low and are thus hard to interpret.However, we can see in the SPW at ∼ 60 • S that reflectances are more limb darkened, which is indicative of lower single-scattering albedo.With our modified 'holistic' atmospheric model parameterisation we then used NEME-SIS to fit spectra reconstructed from these latitudinally-resolved Minnaert coefficients at the same two zenith angles used in the disc-average analysis (i.e., 0 • and 61.45 • ) and retrieved latitudinal distributions of atmospheric properties using three different model setups.These are summarised in Table 2 and described in detail in the following subsections.
Fixing the n imag spectra of the three aerosol types to those determined from the disc-average analysis, we first fitted the latitudinally-resolved reconstructed spectra to retrieve latitudinal distributions of: 1) the opacity of all three aerosol layers (opacities are quoted at 800 nm); 2) the pressure of the Aerosol-2 layer; and 3) the deep methane abundance (Model 1 of Table 2), i.e., five parameters in total.The upper tropospheric methane relative humidity was again limited to 50%.In these retrievals the errors on the Minnaert-reconstructed spectra were increased by a factor of 2 to reflect the fact that fewer points were averaged to determine the coefficients and to ensure that the fitted χ 2 /n values were ∼ 1.The results of this analysis are shown in the form of projected images in Fig. 6 and as plots against latitude in Fig. 7.These retrievals show several very clear latitudinal dependencies of the fitted parameters, although the agreement between the observed and fitted reflectivities at 511 and 831 nm is not perfect (Fig. 7).The retrieved deep abundance of methane is seen to decrease from 6-7% at equator to ∼ 3% at the pole, a very similar latitudinal dependence to that derived from previous analyses of Neptune's visible/near-IR spectra (Karkoschka & Tomasko, 2011;Luszcz-Cook et al., 2016;Irwin et al., 2021), with a boundary at ∼ 20 − 30 • S.Although the latitudinal dependence of methane abundance is similar to previous estimates, however, the absolute abundances are slightly higher and we note in the Discussion that these depend also on the assumed vertical profile of methane and the modelled aerosol structure, which are here different.In addition to methane, we also find significant latitudinal variation in the opacity of all three aerosol layers, with variations in the opacity of the Aerosol-2 layer appearing to match the latitudinal variation seen in the 831-nm image (Fig. 1), with opacity peaks at 75 • S, 45 • S and 20 • S. Higher in the atmosphere, similar, but offset, peaks in opacity can be seen in the retrieved Aerosol-3 opacity, with notable peaks at 80 • S, 60 • S and at the equator.The Aerosol-3 opacity shows some similarity with features in the 848nm image, which is a wavelength of medium methane absorption, but is better correlated with the features in the 860-nm image, where methane absorption is strong.In particular, the peak in Aerosol-3 opacity at 80 • S corresponds well with a faint ring of brightness near Neptune's south pole at 860 nm and similarly with a brighter zone at the equator.Finally, there is a good correlation between the Aerosol-1 opacity and variations seen at 551 nm, and a sharp reduction in the retrieved opacity near 60 • S, which coincides with the SPW, with the opacity increasing south of this, corresponding to the brighter reflectivities seen south of the SPW at 511 and 831 nm.Otherwise, the Aerosol-1 opacity and deep methane abundance seem moderately correlated between the equator and 60 • S, which might suggest upwelling in the 25 • S -25 • N region, condensing more H 2 S cloud at ∼ 5 bar, and downwelling elsewhere.It should also be noted that the spatial variation in the pressure of the Aerosol-2 layer is very small, as expected from previous cloud and methane retrievals for Uranus and Neptune (Karkoschka & Tomasko, 2009, 2011).Finally the χ 2 /n of the fits is good and generally less than ∼2, although there are clearly deficiencies at some wavelengths such as 511 and 831 nm.The a priori values assumed for the retrieved parameters are overplotted in Fig. 7 for reference.The a priori error for all parameters was set to 100%, except deep methane abundance, for which a more constrained error of 50% was assumed.
Although the spectral properties of the SPW at ∼ 60 • S can, to first order, be explained by a thinning of the Aerosol-1 layer, the agreement with the observed reflectivities at 511 nm in Fig. 7 is not perfect, and the model does not match the 831 nm reflectivities very well.Instead, we wondered if the SPW might be caused by a spectrally-dependent darkening of the Aerosol-1 layer.The values of many properties are varying with latitude towards the pole, but by subtracting from the (I/F ) 0 spectrum at ∼ 60 • S the average of the (I/F ) 0 spectra either side (i.e., those at 50 • S and 70 • S) we can attempt to isolate the signature of the SPW itself, obtaining the difference spectrum shown in Fig. 8.This spectrum is compared with the difference spectrum of the NDS-2018 dark spot (Irwin et al., 2023), which is the difference between the observed I/F spectrum in the dark spot and that expected at the same location from the average limb-darkening at this latitude (15 • N). Figure 8 shows that there is considerable similarity between the SPW and NDS-2018 difference spectra, both showing a darkening at wavelengths < 650 nm, with strong methane absorption features indicating the features lie at considerable depth in Neptune's atmosphere.Irwin et al. (2023) found that the darkness of the NDS-2018 dark spot could only be explained by a spectrally-dependent darkening of the Aerosol-1 particles, reducing their single-scattering albedo at λ < 650 nm, and it is quite possible, and indeed arguably probable, that the SPW darkening is caused by the same mechanism.To test this we repeated our latitudinally-resolved retrievals, but this time we also retrieved the imaginary refractive index spectrum of the Aerosol-1 particles, retrieving this at eight equally-spaced wavelengths from 300 -1000 nm, spaced by 100 nm, and with a correlation length of 100 nm (Model 2 of Table 2).As noted earlier, the real part of the complex refractive index spectrum was approximately reconstructed with a Kramers-Kronig analysis and the particle scattering properties calculated with Mie theory, smoothing the phase functions with combined Henyey-Greenstein functions.as the observed spectra at these locations minus the expected spectra at the same locations calculated from the centre-to-limb reflectivity functions at these latitudes of 15 • N and 10 • N, respectively.The comparability between the two sets of difference spectra points to a likely common origin for all these features, namely spectrally-dependent single-scattering albedo spectrum perturbations of the Aerosol-1 layer at ∼5 bar.
Allowing the Aerosol-1 n imag spectrum to vary leads to greatly improved fits to the Minnaert-reconstructed reflectivities as shown in Fig. 9, with closer fits to the re-constructed nadir reflectances at all wavelengths, including those at 551, 831 and 860 nm, shown here.Although the overall retrievals of Aerosol-3 opacity and methane abundance are similar to those previously determined (Fig. 7) the variation in Aerosol-1 opacity shows a less pronounced minimum at 60 • S and the Aerosol-2 opacity and base pressure both seem less correlated with observable variations, although a reduction in Aerosol-2 opacity and pressure level are retrieved at 80 • S.However, the goodness of fit is greatly improved at all latitudes and especially south of 50 • S, with χ 2 /n at 60 • S reducing from 1.23 to 0.99.Of course, the retrieval including n imag has additional model parameters, so we should use the reduced χ 2 /(n y −n x ) statistic here, where n y is the number of spectral points and n x is the number of model parameters.In these retrievals, the total number of fitted spectral points is 374, while the total number of fitted parameters is 5 for the fixed n imag case, and 13 for the variable Aerosol-1 n imag case.This leads to χ 2 /(n y − n x )=1.24 at 60 • S for the retrieval with fixed n imag and 1.02 for the variable case, which is still clearly a very significant improvement.Fig. 10 shows the latitudinal dependence of the computed single-scattering albedo (SSA) of the Aerosol-1 particles at 500 nm, and we can see a very good correspondence between this and the location of the SPW in the 551-nm image (Fig. 1), with the single-scattering albedo, ϖ, of the Aerosol-1 particles significantly reduced at 60 • S. The reduction of ϖ north of 10 • N probably arises due to contamination with the NDS-2018 dark spot not being wholly deconvolved, combined with the fact that the zenith angles are less well sampled at these latitudes and so Minnaertreconstructed spectra are less reliable.We note here that Karkoschka (2011a) also inferred a very similar latitudinal variation in the IR single scattering albedo (their Fig. 20).
If the dark SPW at ∼ 60 • S can be explained by a darkening of the Aerosol-1 particles, might a similar process help explain the banded structure seen at 831 nm?The right hand panel of Fig. 10 shows the latitudinal variation of the Aerosol-1 particles' singlescattering albedo at 800 nm.Here, apart from a general trend of increasing ϖ towards the south pole, we see local peaks of ϖ at 0 • , 20 • , 45 • and 70 • , which correspond to peaks in the observed reflectivity at 831 nm (Fig. 1).This suggests that these 831-nm features (and similar banding seen at other longer wavelengths of minimal methane absorption near 680, 750 and 930nm) are caused by a brightening/darkening of the Aerosol-1 particles at these wavelengths.If this is the case, then the effect is very similar to that inferred for the spectral properties of the 'Deep Bright Spot', DBS-2019, seen near NDS-2018 and reported by Irwin et al. (2023).To test this hypothesis, Fig. 8 compares the difference spectrum at 20 • S (i.e., the difference between the (I/F ) 0 spectrum at 20 • S and the average of those at 10 • S and 30 • S) with the difference spectrum of DBS-2019.As can be seen there is again a considerable degree of similarity between these difference spectra, with the difference peaks longer than 650 nm being mostly restricted to narrow bands centred on 750 and 830 nm.These spectra are not identical, however, with some differences seen at shorter wavelengths.The DBS-2019 feature is small and it is possible that there remains some spatial mixing of the light from nearby locations left over from incomplete deconvolution.However, the similarity of these spectra at longer wavelengths, especially the narrowness of the reflectivity peaks, leads us to the conclusion that these features likely have a similar cause.Hence, we conclude that the SPW and NDS-2018 features are both caused by a spectrally-dependent darkening of the Aerosol-1 particles at λ < 650 nm, while the bright 'zones' seen at 831 nm and DBS-2019 are both caused by a spectrally-dependent brightening of the same Aerosol-1 particles at λ > 650 nm.Furthermore, the 511 nm and 831 nm features must be caused by spectrally-dependent perturbations of the Aerosol-1 layer, rather than changes in its opacity, since opacity changes would affect both wavelengths while it is clear that the latitude dependence of these features are very different.
Although the formal retrieval errors of the opacity and pressure level of Aerosol-2 shown in Fig. 9 are smaller than the latitudinal variations, we note that these errors do not include cross-correlation with other model parameters and so may be underes-timates.Given that we find that the bright and dark features in our data can mostly be explained by variations in the reflectivity of the Aerosol-1 particles only, we explored whether we could fit the observations by fixing the properties of Aerosol-2 and varying only the other remaining parameters (τ 1 , n imag1 (λ), τ 3 , and deep methane abundance).Fixing τ 2 and p 2 at all latitudes to the mean values found of 1.75 (at 800 nm) and 2.1 bar, respectively (Model 3 of Table 2), we re-ran our retrieval model, whose results are overplotted in red in Fig. 9.The quality of the fits can be seen to be only very slightly reduced, with the same trends seen in all the other parameters.Hence, we conclude that to first order the Aerosol-2 layer is latitudinally invariant and that the vast majority of the latitudinal changes seen in the MUSE data can be attributed to changes in: 1) the methane abundance; 2) the opacity and scattering properties of the Aerosol-1 layer at ∼5 bar; and 3) the opacity of the upper tropospheric Aerosol-3 layer.We will return to discuss retrieval errors in the discussion section.

Discussion
The presence of bright 'zones' and dark 'belts' in the Aerosol-1 layer at wavelengths of minimal methane absorption near 830 nm (and also less clearly at 670, 750, and 930 nm), with the zones separated by ∼ 25 • is not something that has specifically been noted before and is very curious given that the measured wind speeds show a much more slowly varying latitudinal dependence.These finer-scale latitudinal variations in reflectivity are, in fact, just visible in HST/WFC3 F845M images (e.g., Chavez et al., 2023), although much less clear, and at longer wavelengths Sromovsky et al. (2014) note finer-scale structure in H-band (1.5 µm) Keck observations of Uranus.Although just visible in broader wavelength filters, the banding seen here is only very clear at wavelengths of extremely weak methane absorption and is not seen at wavelengths of stronger methane absorption.Only high-resolution IFU spectrometers such as MUSE can discriminate between these deep features and shallower cloud structures that would appear identical in broader filters, and thus isolate their surprising origin.Such very different structure below the visible haze layers is reminiscent of Cassini/VIMS observations of the 5-µm emission from Saturn's atmosphere, which revealed finely detailed bright and dark bands that bore little resemblance with the broader structures seen at lower pressures (e.g., Baines et al., 2006).This suggests that the circulation of Neptune's atmosphere below the 2-3-bar Aerosol-2 ice/haze layer is likely very different from that seen above that layer.
Analysing these VLT/MUSE data we have shown that Neptune's SPW and dark spots are likely to be caused by the same process, i.e., a darkening of the particles in the ∼5-bar-Aerosol-1 layer at wavelengths less than 650 nm.Our analysis of the SPW signature is helped by the zenith-angle coverage allowing us to constrain the limb-darkening properties, but is complicated by the difficulty in extracting the SPW signature from the latitudinal variations in other parameters.In contrast, the NDS-2018 difference spectrum is relative to other locations in the same latitude band, so extracting its signature is not complicated by latitudinal variations.However, during the time that NDS-2018 was visible in our data set, from 'Obs-6' and 'Obs-10', the observed zenith angle only increased from 42.36 • to 45.37 • and thus we do not have sufficient zenith angle data to constrain its limb-darkening characteristics.However, IRW22 analysed the limb-darkening properties of both the SPW and GDS from Voyager-2 imaging observations of Neptune in 1989 at several wavelengths and found they had similar limb-darkening, which again points to a likely common cause of these short-wavelength dark features.At longer wavelengths, the signature of the deep bright spot DBS-2019 is similar to the differences between the bright 'zones' and dark 'belts' seen at longer continuum wavelengths, such as 831 nm, and both seem to be caused by a brightening of the particles in the ∼5 bar-Aerosol-1 layer at wavelengths > 650 nm.It is intriguing that two very distinctive and very different features are caused by modifications of the same Aerosol-1 layer at wavelengths either less than or greater than 650 nm, and we thus looked to determine whether there might be a simple explanation for these very different reflectance signatures.
In Fig. 11 we have again plotted the n imag spectrum of the Aerosol-1 particles retrieved from our disc-average Minnaert limb-darkening analysis.As can be seen we re- trieve low values near 500 nm, making the particles highly scattering (ϖ = 0.999), rising to large values at 800 nm, making the particles poorly scattering (ϖ = 0.777).The increase of n imag with wavelength is commonly seen with photochemically-produced particles as discussed by IRW22, who conclude that the particles in Aerosol-1 layer are likely composed of a mixture of bright H 2 S ice particles and photochemical products mixed down from above.The large single-scattering albedo difference on either side of 650 nm is reminiscent of the difference between deep dark and bright spots and we wondered what would happen to the computed disc-averaged (I/F ) 0 spectrum if we added an additional component of particles to the Aerosol-1 layer with different, fixed n imag spectra.The n imag values we chose were 5 × 10 −5 , 5 × 10 −4 , 5 × 10 −3 , 2 × 10 −3 and 1 × 10 −3 (overplotted in Fig. 11), roughly spanning the range of n imag seen in the retrieved Aerosol-1 spectrum.We also tested whether the particle size might be important and used a size distribution for the additional particles either equal to that assumed for the Aerosol-1 layer with mean radius r 0 = 0.1 µm, or increased to r 0 = 1.0 µm.The particle size variance in both cases was set to σ = 0.05.The additional particles had the same vertical distribution as the existing Aerosol-1 layer and had an opacity (at 800 nm) relative to the existing Aerosol-1 layer of 10% for the 0.1-µm particles or 30% for the less backscattering 1.0-µm particles.
The differences in the computed (I/F ) 0 spectra are shown in Fig. 12, which demonstrates that when the additional particles have larger n imag values, and thus lower singlescattering albedo, we see a darkening at shorter wavelengths, but little effect at longer wavelengths.This is understandable given that the additional darker particles will have a large effect at wavelengths where the surrounding Aerosol-1 particles are highly scattering, but will have a minimal impact at wavelengths where the surrounding particles are already rather absorbing.Similarly, adding highly-scattering particles to the Aerosol-1 layer, with low n imag , will not make much difference at wavelengths where the particles are already highly scattering, but will significantly affect wavelengths where the surrounding particles have lower albedo.For either additional particle size, when n imag is high (e.g., 1×10 −3 ), the difference in the computed (I/F ) 0 spectrum matches closely to the SPW and NDS-2018 difference spectra (Fig. 8).On the other hand, for highly scattering particles (e.g., n imag = 5 × 10 −5 ) we see not only narrow reflectance peaks at longer wavelengths, but also an increase in reflectance at shorter wavelengths.However, Figure 12.Change in calculated disc-averaged (I/F )0 spectrum when an additional opacity of particles with specified nimag is added to the particles in the Aerosol-1 layer.The additional particles either have a mean radius of 0.1 µm (black) or 1.0 µm (red) and the additional opacity (at 800 nm) is 10% or 30% of the existing Aerosol-1 opacity, respectively.In both cases the variance of the additional size distribution in 0.05.
for small particles and n imag in the range of 5×10 −4 to 5×10 −3 the effect on the (I/F ) 0 spectra is very similar to the 20 • S and DBS-2019 difference spectra of Fig. 8, with narrow reflectance peaks at longer wavelengths and little clear signature at shorter wavelengths.Hence, rather than trying to explain why the average scattering spectra of the Aerosol-1 particles in the 831-nm bright regions and 500-nm dark spots are so different from the background, we could instead consider simpler solutions where dark spots are caused by the addition of dark particles, or 'chromophores', that are generally absorbing across the MUSE range and 831-nm bright regions are caused by the addition of other particles that are generally scattering, such as fresh H 2 S ice.
Although this chromophore model is simple and easy to understand, it does then leave the question of what the additional chromophores might be.An alternative explanation for dark regions in Neptune's atmosphere, noted by IRW22, is that warm regions at the Aerosol-1 pressure level might lead to H 2 S ice sublimating off from the particles, revealing their darker photochemically-produced cores.This would have more of an effect at wavelengths where the reference Aerosol-1 particles are bright.Equally, cooler regions might result in more H 2 S ice condensing on to the particles, which would not make them more reflective at wavelengths where they are already bright, but could lead to them becoming brighter at the longer wavelengths where we see belts and zones.Observations of Neptune with VLA and ALMA (Tollefson et al., 2021) note latitudinal variations in brightness temperature at wavelengths sounding 4 -8 bar that are similar to the reflectivity variations seen in our longer wavelength (> 600 nm) methane windows, which suggests there might indeed be a link.However, the direction of such a link is arguable: warmer regions might lead to H 2 S ice sublimating off the haze cores to make the particles darker, but it might also be that darker aerosol particles will absorb more sunlight and thus increase the local atmospheric warming.
While the 'zones' at 831 nm are consistent with a brightening of the Aerosol-1 particles, we can see from Fig. 10 that there is general increase of the 800-nm single-scattering albedo (SSA) towards the south pole.In Fig. 13 selected retrieved parameters (for Model 3, where the Aerosol-1 n imag spectrum was allowed to vary, but the opacity and pressure level of the Aerosol-2 layer fixed) are plotted as disc images.The SSA image at 500 nm shows generally good correspondence with the brightness of the 511-nm image in Fig. 1 and the SSA image at 800 nm shows good correspondence with the 831-nm image south of the equator, except north of the equator, where the Minnaert-reconstructed spectra are less reliable.In particular, the bright 831-nm zone seen at 75 • S (which was noted by Irwin et al. (2023)), just south of the SPW corresponds with a region of high 800-nm SSA, and another such region lies on the SPW's northern edge.Might it be that fresh material is upwelling on either side of the SPW, then moves into the SPW, suffering photochemical alteration along the way and darkening it to form the 500-nm-dark SPW feature?A possible candidate for the fresh material might be H 2 S ice, and gaseous H 2 S was detected in Gemini/NIFS observations of Neptune in 2010 (Irwin, Toledo, Garland, et al., 2019), who found its signature (near 1.5 µm) to be stronger towards Neptune's south pole, indicating higher abundance at 2-3 bar.However, deeper in the atmosphere, observations with the Very Large Array and ALMA (Tollefson et al., 2021) find low thermal emission at equatorial latitudes and high thermal emission at polar latitudes.At the depths sounded at these microwave wavelengths (∼10 bar), low thermal emission regions are interpreted as being regions of high H 2 S abundance and thus Tollefson et al. (2021) concluded that H 2 S is enriched at equatorial latitudes, just as we see in our CH 4 retrievals, indicative of strong and persistent upwelling.It may be that the H 2 S signature seen by Gemini/NIFS was obscured near the equator by the enhanced Aerosol-1 opacity we retrieve in this study, making it easier to detect at more polar latitudes near 1.5 µm.Or is it the case that the latitudinal distribution of H 2 S is different at different pressure levels?
Much higher in the atmosphere the ring of stratospheric haze (Aerosol-3) surrounding the south pole at 80 • S is intriguing.We can see the feature in our MUSE data at all methane-absorbing wavelengths longer than 650 nm and also searched for it in HST/WFC observations made during the same apparition, detecting it in the FQ619N and FQ727N filters (both narrow-band filters, centred on wavelengths of strong methane absorption), but not finding it in the broadband HST/WFC3 filters.The ring of material clearly resides high in the atmosphere, near the tropopause and perhaps marks the edge of the south polar vortex in Neptune's upper atmospheric circulation, although Fletcher et al. (2014) andde Pater et al. (2014) put the edge of the south polar vortex, or the prograde polar jet, at ∼ 70 • S. In addition to the bright ring at 80 • S at methane-absorbing wavelengths we also see a slight increase of reflectivity at the equator.Such a distribution indicates upwelling of air at mid-latitudes followed by photodissociation of methane leading to haze production, with the air then moving the haze north and south and concen-trating the aerosols at the equator and poles before it downwells (de Pater et al., 2014;Tollefson et al., 2021).
The retrieved 'deep' abundance of methane varies from 6-7% at the equator to ∼3% at polar latitudes, with a boundary at 20-30 • S.These abundances are slightly higher than previously retrieved (Karkoschka & Tomasko, 2010;Irwin, Toledo, Braude, et al., 2019;Irwin et al., 2021), but we know that the absolute abundance depends both on the assumed vertical profile of methane and also the assumed vertical profile of clouds/hazes.Our scheme has a methane profile that is fixed up to the condensation level, falling with a prescribed relative humidity above, and is combined with three simply vertically parameterised aerosol layers.Although we show that this model matches the data very well, alternative parameterisations could fit similarly well, but return different absolute methane abundances.In particular, Karkoschka andTomasko (2009, 2011) and Sromovsky et al. (2011) favour a 'descended' methane profile that is depleted in the lower trosposphere (at 1-5 bar), but returns to the same value at all latitudes at great depth.Indeed, if latitudinal variations in methane really extended to great depth, then Sromovsky et al. (2011) notes that there would be a significant gradient of mean molecular weight with latitude and significant gradients of vertical wind shear (e.g., Sun et al., 1991;Tollefson et al., 2018) that are inconsistent with the observed winds.Hence, while we are confident of the shape of the retrieved latitudinal abundance of methane, we can be less sure of the absolute profile and abundances.However, the general increase at longer MUSE wavelengths of the latitudinally-resolved Minnaert nadir reflectivity (I/F ) 0 south of 20−30 • S, relative to the disc-averaged (I/F ) 0 spectrum, noted in Fig. 5, can now simply be interpreted as being caused by the lower methane abundances determined at these latitudes.
The circulation of Neptune's atmosphere has been studied for decades, but is still puzzling.In the horizontal direction, we know that the winds are predominantly zonal (i.e., east-west), with strongly retrograde winds at the equator with speeds of nearly 400 m/s, becoming prograde polewards of ∼45 • N and ∼45 • S, reaching peak prograde speeds of over 200 m/s at ∼75 • N and ∼75 • S, before falling to zero at the poles (Sromovsky et al., 1993;Sánchez-Lavega et al., 2019).In the vertical direction, the presence of thick upper-tropospheric methane ice clouds and cooler retrieved upper-tropospheric temperatures at mid-latitudes (20-40 • N and 20-40 • S) indicates upwelling at these latitudes and cooling as the air adiabatically expands at the tropopause and moves towards the equator and poles (Conrath et al., 1991;Bezard et al., 1991;Fletcher et al., 2014;de Pater et al., 2014;Roman et al., 2022).This circulation is consistent with the distribution of Aerosol-3 particles we detect in our MUSE observations.Deeper down, the significant vertical reduction of the mean molecular weight at the CH 4 condensation level at ∼ 2 bar seems to act as a barrier to vertical motion (e.g., IRW22) and below this the vertical circulation appears to switch with air rising at equatorial latitudes and falling nearer the pole, leading to higher equatorial abundances of CH 4 , seen here, and also high deep H 2 S abundances seen at microwave wavelengths (Tollefson et al., 2021).This simple 'stacked circulation' view (e.g., Tollefson et al., 2019;Fletcher et al., 2020), however, is inconsistent with the banded structure seen here at methane continuum wavelengths longer than 650 nm and may also be inconsistent with retrieved H 2 S abundances (at 1.5 µm) at pressures less than 5 bar (Irwin, Toledo, Garland, et al., 2019), suggesting that the circulation may actually be even more complicated than previously thought.An interesting analogy is the distribution of NH 3 seen in Jupiter's deep atmosphere by the Juno spacecraft (Li et al., 2017), which appears to rise in an equatorial plume, with much greater abundances than at mid-latitudes, although a more zonal NH 3 structure is seen at the 0.5-2-bar level by both Juno (Fletcher et al., 2021), and also in mid-infrared observations (e.g., Fletcher et al., 2016); this is believed to arise and be maintained by a similar double 'stacked cell' system (e.g., Ingersoll et al., 2000;Showman & de Pater, 2005;Duer et al., 2021;Fletcher et al., 2021;Moeckel et al., 2023;de Pater et al., 2023).
The data analysed in this study have been fitted by varying just the following few variables: the opacity of the three aerosol layers, the pressure level of Aerosol-2 layer, the 'deep' methane abundance, and the imaginary refractive index spectra of the three aerosol types.The scattering properties of a distribution of particles in a planetary atmosphere depend on very many factors, including the size distribution, composition, and mix of different particle species.Regrettably, for Neptune we have very little a priori constraint on what the expected particles and their size distributions should be, nor any reliable information on their complex refractive index spectra.As a result, have to try and retrieve these from the data themselves if we are fit the observations satisfactorily.We do this using optimal estimation method (Rodgers, 2000), expanded upon below, which is a deepest descent approach and thus is most stable when the rate of change of radiance with model parameter is a smoothly varying function of the value of that parameter.Unfortunately, we find that the rate of change of radiance does not vary smoothly enough with the mean radius of the size distribution to reliably retrieve this with our model.Hence, instead we conduct separate retrievals for a range of fixed particle size distributions and choose those that match the data best (e.g., IRW22).We do find, however, that the model behaves well when varying the n imag spectra, and we do find that we need for the scattering properties of the particles to vary more rapidly with wavelength than simple changes in r mean or the variance of the size distribution can achieve.We attribute this to the fact that the particles are likely made up of unknown photochemical products that have strong absorption bands.We thus fix the size distribution of the particles in a certain aerosol layer to an acceptable shape (compared with a range of observations) and then retrieve the n imag spectra, reconstructing the n real spectra using the Kramers-Kronig approach.In reality the particles may be a mix of separate ice and photochemical particles, but they could also be combined together by 'riming' or some other combination mechanism.Unfortunately, we just do not have this information and so we instead retrieve the mean scattering properties of the layer by fitting n imag as we describe.Although an approximation, we find that our approach allows us to fit the data well and we can then work on interpreting what we have found.This is precisely the point of the first part of this discussion, where we show that the variation of the n imag spectra with latitude of the 'mean' particles in the Aerosol-1 layer can be interpreted as being due by the presence of a background distribution of particles, presumably rich in photochemical haze, combined with a latitudinally-varying opacity of fresh ice crystals.It may be that these two particle distributions have different mean radii and variances, but these cannot be uniquely separated using these data and so we assume a single mean size distribution.Ideally, we would have the complex refractive index spectra of a set of laboratorymeasured possible condensates to test against the measured spectra, but this data set does not yet exist.In the meantime, we believe that we have developed a method that reliably retrieves the mean scattering properties of the aerosol layers, which can then be interpreted to gain novel insights into the hazes and clouds in Neptune's atmosphere.
Finally, we note that the retrievals presented here were conducted using the optimal estimation (Rodgers, 2000) framework of the NEMESIS model (Irwin et al., 2008).To model these spectral observations requires a full multiple-scattering radiative transfer model, which must be run at many latitudes, many wavelengths, and with many tunable parameters.Hence, this model is computationally very expensive and to combine it with a Bayesian retrieval framework (such as Monte-Carlo Markov Chain or Nested Sampling, e.g., Skilling (2006)), which require tens of thousands of iterations, would be prohibitively computationally expensive.However, although cross-correlation between parameters can to be interpreted in optimal estimation method from the computed covariance matrix, this is much less easy to comprehend and present than the 'corner plots' of Bayesian approaches.The retrieval errors we present in this work are derived from diagonal elements of the retrieved covariance matrix and corrected for the original a priori errors, since if the retrieved error is the same as the a priori error we have not learnt anything new.Hence, the errors shown are σ = 1/ (1/σ 2 ret − 1/σ 2 apr ), where σ 2 ret are diagonal components of the retrieved covariance matrix and σ 2 apr are diagonal compo-nents of the a priori covariance matrix.However, these errors still do not wholly account for cross-correlation effects.Instead, we estimate the magnitude of these effects by fixing some parameters and observing the effect on the other retrieved variables.For example, going from Model 2 to Model 3 we show that by fixing the pressure and opacity of the Aerosol-2 layer we can achieve very similar quality fits and very similar variations with latitude of the other retrieved parameters.How, then, can we be sure that the retrieved latitudinal variations of the remaining parameters in Model 3 are reliable?Figure 8 shows that dark spots, the SPW and the zone/belt signature at longer continuum wavelengths are all most consistent with changes in the Aerosol-1 layer, and that this study and previous analyses (Irwin, Teanby, Fletcher, et al., 2022;Irwin et al., 2023) find that spectrally-dependent changes in the reflectivity of the Aerosol-1 particles provide the best solution, which is what Model 3 achieves (Fig. 10).There are no clear latitudinal variations that have a spectral signature consistent with significant latitude variations in the Aerosol-2 layer, but variations seen at methane-absorbing wavelengths, and which thus must be high in the atmosphere, are matched by latitude variations in the opacity of the Aerosol-3 layer and there is a clear methane variation signature in the observations, noted earlier.Hence, our Model 3 is consistent with all the evidence in the MUSE data set and is, we believe, reliable.We also note that as this model was based on the 'holistic' model of IRW22, which was originally derived from a much wider wavelength range data set, it is likely to be applicable to longer wavelengths also.Observations of Neptune have recently been made with the NIRSpec and MIRI instruments on the James Webb Space Telescope (JWST) and we look forward to seeing if our parameterisation can be extended successfully to the longer wavelengths sampled by JWST also.Latitude circles are overplotted for ease of reference, with a spacing of 30 • .

Conclusions
Our analysis of latitudinally-resolved centre-to-limb spectra of Neptune, observed with VLT/MUSE in October 2019, has revealed new constraints on the mechanism that causes the darkening of the South Polar Wave at ∼ 60 • S, which we conclude is common with the darkening mechanism for discrete dark spots such as Voyager-2's Great Dark Spot and the more recent NDS-2018 feature.Using a modified version of the 'holistic' aerosol model of IRW22 we find that both these features are consistent with being caused by the darkening (at wavelengths less than 650 nm) of particles in an Aerosol-1 layer, based at ∼5 bar, which we suggest is composed of a mixture of photochemically-produced haze, mixed down from their production level above, and H 2 S ice.In addition, we note a latitudinal brightening and darkening (with a scale of ∼ 25 • ) in Neptune's aerosols, which is only visible at regions of very low methane absorption longer than 650 nm.We show that this feature is caused by a brightening (at wavelengths greater than 650 nm) of the particles in the same Aerosol-1 layer at ∼ 5 bar.We conclude that these different features seen at λ < 650 nm and λ > 650 nm must be caused by spectrally-dependent perturbations of the Aerosol-1 scattering properties, rather than opacity changes of this layer, since they have very different spatial structures and changes in the opacity would affect both wavelength ranges similarly.
Above the Aerosol-1 layer at ∼5 bar, we find the properties of the main Aerosol-2 layer, confined to a region of high static stability at the CH 4 -condensation level at ∼2 bar, is to first order invariant with latitude, although the base pressure and opacity may be slightly reduced at 80 • S.Even higher in the atmosphere, variations in the opacity of the Aerosol-3 tropospheric haze are able to reproduce the observed variation in reflectivity at methane-absorbing wavelengths.We find the abundance of this aerosol to be higher at the equator and also in a narrow 'zone' at 80 • S, which supports the hypothesis that at the tropopause air is rising at mid-latitudes and then moving polewards and equatorwards, concentrating the photochemical haze products there.The bright ring at 80 • S would appear to define the edge of a polar cyclone that has been observed on Neptune for the past two decades, and which are often seen in giant planet atmospheres.
The detailed conclusions of this paper are: 7. We find the South Polar Wave at ∼60 • S to be caused by a perturbation of the particles in the Aerosol-1 layer that reduces the single-scattering albedo at shorter wavelengths.The darkening appears identical to that which darkens the dark spot NDS-2018 and could be caused either by the addition of a chromophore in this layer that is darker than the background particles at all MUSE wavelengths, or by local heating at the ∼5-bar level that sublimates H 2 S ice off to reveal the dark haze cores of the Aerosol-1 particles; 8.We note a newly apprehended latitudinal variation in the reflectivity in narrow continuum wavelengths longer than 650 nm.We interpret these changes as being caused by variations in the scattering properties of the Aerosol-1 layer at wavelengths longer than 650 nm, with bright 'zones' at the equator, 20 • S, 45 • S, and 75 • S, interspersed with darker 'belts' at 5 • S, 30 • S, 50 • S, and 80 • S.This coloration could be caused by the addition of brighter, more scattering particles at the zone latitudes, or by local cooling at the ∼5-bar level that condenses more bright H 2 S ice on to the Aerosol-1 particles.
The banding seen in the narrow windows of very low methane absorption longer than 650 nm is suggestive of a finer latitudinal scale variation than has been noted before on Neptune at the level of the H 2 S/haze Aerosol-1 layer at ∼5 bar and indicates that the circulation of Neptune's atmosphere is even more complicated than that suggested by the 'stacked circulation' models of Tollefson et al. (2019) and Fletcher et al. (2020).Further work is needed to explore such changes more fully and recent observations with the James Webb Space Telescope should provide strong new constraints.In addition, the spatial deconvolution techniques developed for this study could be applied to existing ground-based IFU measurements from instruments such as Gemini/NIFS (Irwin et al., 2011) and VLT/SINFONI (Irwin et al., 2016) to extend the high spatial resolution analysis to longer wavelengths, and perhaps better constrain the latitudinal dependence of H 2 S abundance.Further in the future, a space mission to one of the Ice Giants, comprising an orbiter and an entry probe would be invaluable in providing some measure of 'ground truth' for aerosol models.In particular, it would be useful to observe with a probe the penetration depth of the UV flux to establish the pressure levels at which we may expect the photolysis of methane and other components.Meanwhile, an orbiter would provide reflectivity observations at higher phase angle, which will help to further constrain the allowable range of solutions that are consistent with observations.formed using the NEMESIS radiative transfer and retrieval algorithm Irwin et al. (2008) and can be downloaded from Irwin, Teanby, de Kok, et al. (2022a) (or https://github .com/nemesiscode/radtrancode),with supporting website information at Irwin, Teanby, de Kok, et al. (2022b) (or https://github.com/nemesiscode/nemesiscode.github.io).

Figure 1 .
Figure 1.Deconvolved MUSE 'slices' at 551, 831, 848 nm, and 860 nm, respectively, from the 'Obs-6' cube.Also shown is the modelled appearance of Neptune to the naked eye, reconstructed from the MUSE observations using standard colour-matching functions and correctly accounting for gamma corrections.A reference latitude and longitude grid, with spacing of 30 • in each direction, is also shown.In the 511-nm and naked-eye images the South Polar Wave (SPW) is just visible at bottom left, while the dark spot NDS-2018 can just be seen at top right.The DBS-2019 bright spot is only visible at longer wavelengths of extremely low methane absorption, such as 831 nm shown here.The 848 and 860 nm images, at methane-absorbing wavelengths, probe the haze high in the atmosphere, with 860 nm probing slightly higher and revealing a south polar collar at 80 • S. The bright spots on the left edge of Neptune's disc at longer wavelengths are upper tropospheric methane ice clouds at 0.1 -0.6 bar.

Figure 3 .
Figure 3. Top panel: MUSE spectra reconstructed at 0 and 61.45 • zenith angles (viewing and solar) from our disc-average Minnaert analysis, and fits to them using our NEMESIS model.

Figure 5 .
Figure 5. Latitudinally-resolved Minnaert limb-darkening analysis of MUSE data.Panel (a) compares the extracted (I/F )0 spectra at the equator (red) and 60 • S (blue).A contour plot of the extracted (I/F )0 spectra as a function of wavelength and latitude (planetographic) is shown in Panel (b), while Panel (c) shows a contour plot of the difference of (I/F )0 spectra compared with the disc-averaged (I/F )0 spectrum (red regions are darker than disc-average, blue are brighter).In the contour plots the equator and 60 • S are highlighted with coloured, dashed lines.Panels (d-f) show the same plots as Panels (a-c), but for the fitted Minnaert limb-darkening coefficients, k.

Figure 6 .
Figure 6.Zonally-averaged meridional profiles of atmospheric properties retrieved from the VLT/MUSE Neptune observations (Obs-6) using a model where the particle scattering properties were fixed to the disc-average for all aerosol layers (Model 1).The following fitted properties are projected on to Neptune's disc: a) Aerosol-1 opacity; b) Aerosol-2 opacity; c) Aerosol-3 opacity (×100); d) deep methane mole fraction (%); and e) Pressure of base of Aerosol-2 layer (bar).The cloud opacities referred to here and elsewhere are those at 800 nm.Latitude circles are overplotted for ease of reference, with a spacing of 30 • .

Figure 7 .
Figure 7. Zonally-averaged atmospheric properties retrieved from the VLT/MUSE Neptune observations (Obs-6), using a model (Model 1 of Table 2) where the particle scattering properties were fixed to the disc-average for all aerosol layers, showing: a) Aerosol-1 opacity; b) Aerosol-2 opacity; c) Aerosol-3 opacity (×100); d) deep methane mole fraction (%); e) Pressure at the base of Aerosol-2 layer (bar), and f) χ 2 /(ny-nx) of the fits at different latitudes.The dotted lines in these panels show the a priori values assumed.Note that the error bars shown are the formal errors on the retrievals and do not include cross-correlation effects (see Discussion).The bottom row (panels g -i) shows the latitudinal variations of (I/F )0 at 551, 831 and 860 nm, respectively, where the solid lines show the Minnaert-fitted values (and uncertainties) and the dashed lines show the best-fit values from our retrieval.

Figure 9 .
Figure9.As Fig.7, but showing atmospheric properties retrieved from the VLT/MUSE Neptune observations (Obs-6), when the nimag spectrum of the Aerosol-1 particles is also allowed to vary (Model 2).In the middle-right panel, showing the latitudinal (planetographic) variation of χ 2 /(ny − nx), the reduced χ 2 for the case where nimag for Aerosol-1 is fixed is overplotted as a dotted line for comparison.Overplotted in red in all panels are the results when the opacity (at 800 nm) of Aerosol-2 is fixed to 1.75 and its pressure fixed to 2.1 bar (Model 3), showing very similar fits to the data and similar reduced χ 2 .Note that the error bars shown in panels a -e are the formal errors on the retrievals and do not include cross-correlation effects (see Discussion).

Figure 10 .
Figure10.Estimated single-scattering albedoes of the best fit Aerosol-1 particles at 500 nm and 800 nm as a function of latitude (planetographic), derived from our fitted nimag spectra using Mie theory.The formal confidence limits indicated by the grey shading and do not include cross-correlation effects.

Figure 11 .
Figure 11.Imaginary refractive index spectrum of the fitted disc-averaged Aerosol-1 particles (green), compared with the nimag spectra of the additional particles added to the Aerosol-1 layer in the reflectivity calculations shown in Fig. 12.
N.B., n imag tabulated from 0.3 to 1.0 µm in steps of 0.1 µm.∆p for Aerosol-1, Aerosol-2 is FWHM of Gaussian function in ln p. Aerosol-3 is extended up from p 3 with a set fractional scale height (FSH).Note that where n imag (λ) is indicated as 'fixed', it means that the spectrum is fixed to that retrieved by the disc-average analysis.