Hazy Blue Worlds: A Holistic Aerosol Model for Uranus and Neptune, Including Dark Spots

2.1. Uranus Spectral Observations

We analysed HST/STIS observations of Uranus from 2002 (Karkoschka & Tomasko, 2009), 2012 (Sromovsky et al., 2014) and 2015 (Sromovsky et al., 2019). Since our primary intention in this paper was not to revisit the question of latitudinal variability of methane abundance or polar brightening, we concentrated on the 2002 data (Karkoschka & Tomasko, 2009), when the disc of Uranus appeared particularly featureless, allowing us to combine the data from all latitudes together. This observation was made on 19 August 2002, between 01:43 and 10:57 UT. HST/STIS is actually a long slit spectrometer, but the slit was aligned parallel to the central meridian and then stepped from the central meridian to the edge of the planet to form a “cube” of half the planet, where at each location on the disc a complete spectrum covering 300.4–1,020.0 nm was recorded at a resolution of 1 nm, sampled every 0.4 nm.

We also analysed an observation of Uranus made using IRTF/SpeX, another long‐slit spectrometer. In this case, the slit was aligned on the disc centre and the fluxes integrated along the central meridian. This standard reference spectrum is available on the IRTF Spectral Library website (http://irtfweb.ifa.hawaii.edu/∼spex/IRTF_Spectral_Library), reported by Rayner et al. (2009), and was made in SXD mode (0.8–2.5 μm) on 18 May 2000. These data have a spectral resolution of 0.002 μm.

Finally, we analysed observations of Uranus made with Gemini/NIFS in H‐band (i.e., 1.47–1.8 μm) in 2009 and 2010 (Irwin et al., 2011b; Irwin, Teanby, et al., 2012). Gemini/NIFS is an integral field‐unit (IFU) spectrometer that simultaneously records spatial and spectral information, and where each pixel or “spaxel” is composed of a complete spectrum covering the targeted spectral range with a spectral resolving power of R ∼ 5,200. The actual cube used for our retrieval study was recorded on 2 September 2009, which had good spatial resolution and was reasonably clear of discrete features.

To enable easy intercomparison and simulation speed the HST/STIS and Gemini/NIFS data were smoothed to the spectral resolution of IRTF/SpeX of 0.002 μm (i.e., 2 nm), and sampled with a step of 0.001 μm.

During the time period spanned by these observations, Uranus moved through its orbit about the Sun and Table 1 lists the date, sub‐Earth latitude (planetocentric) ϕ _E , and also the apparent target‐centered longitude of the Sun L _s (which gives a convenient measure of season, with 0° being northern spring equinox, 90° being northern summer solstice, etc.) of the observations. In addition, we list the disc‐integrated photometric magnitudes of Uranus as observed and reported by Lockwood (2019) in the y (551 nm) and b (472 nm) filters of the Strömgren photometric system, which gives a measure of the overall disc brightness and blueness. It can be seen that the IRTF/SpeX and HST/STIS observations were both made with L _s  ∼ 335° and ϕ _E  ∼ −25°. However, the Gemini/NIFS observation comes from just after the northern spring equinox in 2007 with ϕ _E  = 7.9°. Hence, while the IRTF/SpeX and HST/STIS observations will be slightly more weighted to conditions in the southern hemisphere, the Gemini/NIFS observations will be more weighted to equatorial regions and in a slightly later season. However, neither set of observations will sample well the hazes in the polar regions. From 2000 to 2009 the disc‐averaged magnitude of Uranus decreased by ∼0.03 and became slightly more blue as the southern polar “hood” or “cap” (e.g., Sromovsky & Fry, 2008) became less visible. However, these changes are relatively small, enabling us to assume that all three sets of observations were observing approximately the same disc of Uranus. However, care should be taken when comparing this work with Voyager‐2 studies of Uranus's aerosol structure, since Uranus was then near southern summer solstice and had a well developed “hood” and markedly higher albedo and thus brightness (Lockwood, 2019).

2.2. Neptune Spectral Observations

We analysed HST/STIS observations of Neptune (Karkoschka & Tomasko, 2011) made on 3 August 2003, between 4:38 and 14:15 UT. As for the Uranus HST/STIS data these form a “cube” of half of Neptune's disc, covering 300.4–1,020.0 nm.

We also analysed a reference observation of Neptune made using IRTF/SpeX and integrated along the central meridian, available on the IRTF Spectral Library (Rayner et al., 2009). The SXD component of this reference spectrum analysed here (0.8–2.5 μm), was observed on 30 June 2000.

Finally, we analysed observations of Neptune made with Gemini/NIFS in 2009 (H‐band) and 2011 (I, J, and H‐band, i.e., 0.94–1.16 μm, 1.14–1.36 μm, and 1.47–1.8 μm, respectively) (Irwin et al., 2011a, 2016). The actual cube used for our retrieval study was recorded in the H‐band on 7 September 2009 and was chosen for its good spatial resolution and reasonably limited distribution of discrete cloud features.

During the time period spanned by these observations, Neptune also moved through its orbit about the Sun and Table 1 again lists the date, sub‐Earth latitude ϕ _E , apparent target‐centered longitude of the Sun L _s , and the disc‐integrated photometric magnitudes (Lockwood, 2019) of the observations. Neptune's southern summer solstice occurred in 2005, but Neptune's slower orbit about the Sun means that L _s and ϕ _E differed little during the total elapsed period, and can be assumed to be ∼270° and ∼−29°, respectively. However, it can be seen that all three data sets will be weighted towards conditions in the southern hemisphere, with the south polar region clearly visible. From 2000 to 2009 the disc‐averaged magnitudes and colour of Neptune were found not to vary significantly (Lockwood, 2019), again enabling us to assume that all three sets of observations were observing approximately the same disc of Neptune. However, care should again be taken when comparing this work with Voyager‐2 studies of Neptune's aerosol structure, and indeed the Voyager‐2/ISS observations described later in this paper, since Neptune was then still approaching southern summer solstice and was noticeably darker (Lockwood, 2019).

Once again, to enable easy intercomparison the HST/STIS and Gemini/NIFS data were smoothed to IRTF/SpeX spectral resolution.

3.1. Atmospheric and Radiative Transfer Modelling

We modelled these observations using the radiative transfer and retrieval tool, NEMESIS (Irwin et al., 2008; Irwin et al., 2022a, 2022b). To account for multiple scattering, NEMESIS employs a plane‐parallel matrix operator model (Plass et al., 1973), where integration over zenith angle is performed with a Gauss‐Lobatto quadrature scheme, while the azimuth integration is done with Fourier decomposition. We have found that five zenith angles are usually sufficient to model the giant planet I/F spectra and the number of azimuth Fourier components is set from the viewing and solar zenith angles as $N_F=⌊\mathrm{\Theta }/3⌋$, where Θ = max[θ,θ ₀] and where θ is the observation zenith angle (in degrees) and θ ₀ is the solar zenith angle (in degrees). This scaling of N _F with Θ is necessary to reconstruct reliably the scattering functions at higher zenith angles, but comes at a considerable computational cost. Calculations at intermediate zenith angles are interpolated between the matrix operator calculations at the nearest tabulated viewing and solar zenith angles. The five‐angle zenith angle quadrature scheme used here is listed in Table 1 of Irwin et al. (2021).

NEMESIS was run in correlated‐k mode, with the methane k‐tables generated from a number of different sources as described in Section 3.3. For H₂–H₂ and H₂–He collision‐induced absorption we used the coefficients of Borysow et al. (1989); Borysow and Frommhold (1989); Borysow et al. (2000), assuming a thermally equilibriated ortho:para hydrogen ratio for both Uranus and Neptune. Rayleigh scattering was included as described in Irwin, Toledo, Braude, et al. (2019), and the effects of polarization and Raman scattering were included as described in Section 3.3. We used the solar spectrum of Chance and Kurucz (2010), extrapolated to longer/shorter wavelengths with the solar spectrum of Kurucz (1993). This combined solar spectrum was smoothed with a triangular line shape of FWHM = 0.002 μm for consistency with the Uranus and Neptune resampled spectra.

3.2. Atmospheric Models

The reference temperature and mole fraction profile for Uranus is the same as that used by Irwin et al. (2018) and is based on the “F1” temperature profile determined from Voyager 2 and HST/STIS observations (Sromovsky et al., 2011), with He:H₂ = 0.131 and including 0.04% mole fraction of neon. We adopted a simple “step” model for the methane profile with variable deep mole fraction (a priori 4%), variable relative humidity above the condensation level and a limiting stratospheric mole fraction of 1 × 10⁻⁴ (Encrenaz et al., 1998; Orton et al., 2014). For the atmospheric static stability study discussed later, we also added saturation‐limited profiles of H₂S and NH₃. When dealing with condensing gases such as CH₄, once the profile for that gas was determined, the abundance of H₂ and He was scaled at each level, keeping He:H₂ = 0.131, to ensure that the mole fractions added up to 1.0.

The reference temperature and mole fraction profile for Neptune is the same as that used by Irwin, Toledo, Braude, et al. (2019) and Irwin et al. (2021) and is based on the “N” profile determined by Voyager‐2 radio‐occultation measurements (Lindal, 1992), with He:H₂ = 0.177 (15:85) and including 0.3% mole fraction of N₂. We again adopted a simple “step” methane model with a variable deep mole fraction, variable relative humidity above the condensation level and limiting the stratospheric mole fraction to 1.5 × 10⁻³ (Lellouch et al., 2010). Again, the abundances of H₂ and He were adjusted at each level to ensure the mole fractions added up to 1.0, keeping He:H₂ = 0.177.

3.3. Spectral Modelling

The HST/STIS and IRTF/SpeX observations of both Uranus and Neptune were first combined to give composite central‐meridian‐averaged spectra of both planets, which are shown in Figure 1. To create these combined spectra, we took the HST/STIS observations of Uranus and Neptune from 2002 to 2003, respectively, averaged them along the central meridian (masking discrete cloud features in the Neptune STIS data) and then scaled the IRTF/SpeX Uranus and Neptune spectra to match at overlapping wavelengths from 0.85 to 1.0 μm. This scaling was necessary as the IRTF/SpeX data are in units of line‐integrated total flux and needed to be converted to I/F. Given uncertainty in the slit size and position it was not possible to do this ab initio with sufficient accuracy and hence re‐scaling was necessary to ensure consistency with HST/STIS data. The Gemini/NIFS data, although not used in this initial analysis, were also averaged along the central meridian and their reflectivities adjusted to match the scaled IRTF/SpeX spectrum at overlapping wavelengths for later use.

Until now, our radiative transfer model NEMESIS has not incorporated Raman scattering, but since the HST/STIS observations extend to 0.3 μm, where Sromovsky (2005a) shows that Raman scattering is significant, it was necessary here to incorporate this effect. To do this, we followed the approach of Sromovsky (2005a) and considered only the S(0), S(1) and combined “Q” transitions, using the cross‐section absorption coefficients of Ford and Browne (1973). Radiation scattered from short to longer wavelengths was introduced at the longer wavelengths as an additional pseudo‐thermal‐emission term, using the assumption of Sromovsky (2005a) that this Raman re‐emitted radiation is effectively isotropic, with wavenumber shifts of 354.69 cm⁻¹, 587.07 cm⁻¹ and 4,160.00 cm⁻¹, respectively for the S(0), S(1) and combined “Q” transitions. In addition to including Raman scattering, we also further revised NEMESIS to incorporate a correction for polarisation effects as described by Sromovsky (2005b).

Using our upgraded NEMESIS model we show in Figure 2 synthetic spectra calculated from our assumed Uranus and Neptune atmospheres for aerosol‐free conditions, including Rayleigh scattering, Raman scattering, H₂‐H₂ and H₂‐He collision‐induced absorption (CIA), and absorption by gaseous methane. As can be seen, the measured reflectivity (I/F) of both Uranus and Neptune decreases rapidly with wavelength and is well matched overall by a 1/λ ⁴ curve, suggesting that any aerosols which are present are small and of low integrated opacity. In Figure 2 we compare calculations made using three different available sources of methane absorption data. The band data of Karkoschka and Tomasko (2010), converted to k‐tables by Irwin et al. (2011a), covers the entire wavelength range and reproduces most of the observed features well. However, at longer wavelengths we can also use k‐tables generated from HITRAN2016 methane line data (Gordon et al., 2017), and a more recent compilation of methane line data from the TheoRETS project (Rey et al., 2018). The HITRAN2016 data can be seen to be in good agreement with the band/k‐data at wavelengths longer than 1.0 μm, while the TheoRETS data extends reasonably well to even shorter wavelengths of 0.75 μm. However, neither set extends to visible wavelengths and since our aim in this study is to find a holistic aerosol model of Uranus and Neptune that matches all wavelengths covered by HST/STIS, IRTF/SpeX and Gemini/NIFS simultaneously we used the band/k‐data of Karkoschka and Tomasko (2010)/Irwin et al. (2011a) in this study. It is also apparent from Figures 1 and 2 that Neptune is more reflective than Uranus at methane‐absorbing wavelengths longer than ∼1 μm, indicating that Neptune has a higher opacity of upper‐atmosphere aerosols, which must be of a size greater than ∼1 μm in order to be visible at these longer wavelengths. Furthermore, it would appear that any aerosols in Uranus's upper atmosphere must be of very low opacity, since a pure Rayleigh/Raman scattering calculation already matches the observed reflectivity at methane‐absorbing wavelengths reasonably well.

In order to analyse the observations, we needed to determine the depths to which sunlight can penetrate at different wavelengths. In Figure 3 we show contour plots of the two‐way transmission from space, for a vertical path, to different levels in the atmospheres of Uranus and Neptune for aerosol‐free conditions, including Rayleigh scattering, Raman scattering, methane absorption and hydrogen‐helium collision‐induced absorption. Here it can be seen that the penetration depth is very similar for Uranus and Neptune, but that a radiative transfer model wishing to simulate the observations needs potentially to cover a very wide range of pressure levels. Hence, the pressure range for both Uranus and Neptune models was set to cover 40–0.001 bar and the atmosphere split into 39 layers (equally spaced in log pressure) for our radiative transfer calculations.

To best analyse the combined HST/STIS, IRTF/SpeX and Gemini/NIFS observations, and place strong constraint on the vertical structure and spectral dependence of particle scattering properties, we needed to develop an efficient way of modelling the observed centre‐to‐limb variations of the HST/STIS and Gemini/NIFS data and combine these with the wider wavelength‐range central‐meridian‐averaged IRTF/SpeX observations, which provide important constraints on particle size and better probe lower pressure levels at the longer wavelengths. To do this we employed the Minnaert limb‐darkening approximation (Minnaert, 1941), where for an observation at a particular wavelength, the reflectivity R is approximated as:

Here, μ and μ ₀ are the cosines of the viewing and solar zenith angles, respectively, R ₀ is the fitted nadir reflectance, and k is the fitted limb‐darkening parameter. Irwin et al. (2021) applied this model to VLT/MUSE observations of Neptune made in 2018 and found it to give a very good approximation to the observations and was also well reproduced by the NEMESIS radiative transfer model. Hence, in this study we used the Minnaert approximation to simplify our calculation of the disc‐averaged limb‐darkening and central‐meridian averaged spectra following the scheme:

1. As previously described, the HST/STIS and Gemini/NIFS observations were first averaged to the IRTF/SpeX resolution of 0.002 μm and sampled on a regular grid of spacing 0.001 μm. The masked HST/STIS data were then averaged along the central meridian and used to scale the IRTF/SpeX observations to make the combined HST/IRTF processed observations self‐consistent. The Gemini/NIFS data were also averaged along the central meridian (masking out discrete cloud features) and themselves scaled to be consistent with the combined HST/STIS and IRTF/SpeX observations.

2. The averaged HST/STIS and Gemini/NIFS “cubes” were then masked to exclude discrete cloud features and the remaining observations at all latitudes Minnaert‐analysed to derive spectra of the disc‐averaged R ₀(λ) and k(λ). This was simplified in this case since Uranus and Neptune are so distant that the solar zenith angle and viewing zenith angle can be approximated to be the same (i.e., μ = μ ₀) and hence R = R ₀ μ ^2k−1.

3. The fitted Minnaert parameters R ₀(λ) and k(λ) were used to generate reconstructed spectra for HST/STIS and Gemini/NIFS at two zenith angles (0°, 61.45°), corresponding to two of the five zenith angles of the Gauss‐Lobatto multiple‐scattering radiative transfer model used in NEMESIS. The higher angle is large enough to ensure we fully capture the observed limb‐darkening/limb‐brightening, and is coincident with one of the Gauss‐Lobatto quadrature angles, thus obviating the need for any interpolation (Irwin et al., 2021). However, it is not so large that the computational cost becomes excessive.

The two resulting reconstructed HST/STIS spectra, two reconstructed Gemini/NIFS spectra and the measured IRTF/SpeX central‐meridian line‐averaged spectrum were then used as a set of “measured” observations to NEMESIS. By fitting an aerosol model to the HST/STIS and Gemini/NIFS processed spectra at 0° and 61.45° zenith angles, and extracting our best fits to the Minnaert parameters R ₀(λ) and k(λ), we could then simply reproduce the observations at any other angle from Equation 1, assuming that the Minnaert approximation holds at all the other zenith angles of our zenith‐angle quadrature scheme, which Irwin et al. (2021) found to be a good approximation for their analysis of VLT/MUSE Neptune observations. However, we still needed an efficient way of simulating the central‐meridian‐averaged IRTF/SpeX spectra. We could have computed spectra at multiple locations along the central meridian and averaged these, but this would have been slow, especially near the disc edges. Instead, for each iteration we calculated spectra at two zenith angles (0° and 42.47°), extracted Minnaert parameters R ₀ and k at each wavelength (using k = (1 + log(R _μ /R ₀)/ log(μ))/2, where R ₀ is the nadir‐calculated radiance and R _μ is the radiance calculated at μ = cos(42.47°)) and used these to compute the central‐meridian line‐average as follows. The reason we chose the middle zenith angle of our quadrature scheme (42.47°) here, instead of the second largest (61.45°) as we did for the HST/STIS observations, is that calculations at lower zenith angles require fewer Fourier azimuth components and are more rapid; since we only needed to approximate the line‐average in this case, rather than limb‐darkening, the lower angle was found to be sufficiently accurate.

In Appendix B we show that the mean radiance integrated along the central meridian of a planet whose limb‐darkening is well represented by the Minnaert approximation can be written (assuming μ = μ ₀) as:

Hence, if we have an estimate of the Minnaert limb‐darkening coefficient, k, and nadir radiance, R ₀, we can easily compute the corresponding central meridian line‐averaged radiance. This function cannot be integrated analytically and so we pre‐computed a table of ${\bar{R}}_{\mathit{line}}/R_0$ versus k and then interpolated this to the value of k derived at each wavelength. Hence, using calculations at just two angles to determine k we were able to accurately simulate the central‐meridian line‐averaged IRTF/SpeX observations using Equation 2. For the HST/STIS and Gemini/NIFS data, we further verified that the central‐meridian‐line‐average calculated from Equation 2 using the disc‐averaged Minnaert parameters was consistent with the measured IRTF/SpeX observations. For the Neptune Gemini/NIFS observations, this test necessitated some iteration on the degree to which discrete clouds were masked in order to generate the best overall set of self‐consistent Minnaert coefficients. All five spectra could then be fitted simultaneously, or individual spectra fitted separately as necessary.

Our starting assumption in this analysis is that the cloud/haze structures of Uranus and Neptune are substantially the same, but that Neptune has thicker upper‐tropospheric/lower‐stratospheric haze. Hence, to begin with we analysed the simpler case of the Uranus observations and extended to Neptune observations later.

3.4. Analysis of HST/STIS Uranus Observations

We first analysed the HST/STIS reconstructed observations of Uranus at zenith angles of 0° and 61.45° from 0.5 to 1.0 μm. We conducted a first‐pass retrieval using a continuous distribution of haze particles, in order to let the data (rather than a priori assumptions) determine the location of aerosol layers. In this first pass, the particle sizes were assumed to follow a Gamma distribution with mean radius 0.1 μm and variance σ = 0.3. We also fitted for the imaginary refractive index spectrum (a priori value set to 0.001) of these particles at six wavelengths from 0.5 to 1.0 μm, step 0.1 μm, and then reconstructed the real refractive index spectrum using the Kramers‐Kronig relation, assuming the real refractive index at 0.8 μm was fixed to a value of 1.4. The resulting complex refractive index spectrum was then used to calculate, using Mie theory, the extinction cross‐section, single‐scattering albedo, and phase function spectra. However, since we expect the particles (ice and haze) to be non‐spherical, the phase function spectra were approximated with combined Henyey‐Greenstein phase function parameters to average over features peculiar to spherical particles such as the “rainbow” and “glory.” For the methane profile we assumed a “step” function with retrievable deep abundance up to the condensation level, but with the relative humidity above the condensation level fixed to 100% and the stratospheric abundance limited to not exceed 1 × 10⁻⁴. The retrieval was attempted with two widely spaced a priori aerosol abundance profiles (assumed to have constant opacity/bar at all pressure levels) to achieve a “bracketed” retrieval, where the retrieved profiles will overlap where they are well constrained by the measurements and tend back to their respective a priori values at pressures where there is little sensitivity. Figure 4 shows the fits to the reconstructed spectra together with the fitted aerosol abundance profiles. As can be seen, the fit struggles slightly at short wavelengths, but otherwise we find that the HST/STIS observations are well fitted by an aerosol structure with a thin, well defined peak at ∼1.5 bar and a second layer at pressures greater than ∼4 bar. Although we can clearly detect the top of this deeper layer, we cannot constrain its base pressure as the bracketed retrievals tend back to their respective a priori values at pressures exceeding approximately 5–7 bar; this makes sense given that the two‐way transmission to space for a vertical path rapidly falls to zero in this pressure region, even for an aerosol‐free atmosphere (Figure 3). This indication of a double‐peaked cloud/haze structure was also noted by Sromovsky et al. (2019), as described earlier.

3.5. Analysis of Combined HST/STIS and IRTF/SpeX Uranus Observations ‐ Snippet Analysis

Having noted that there may be two aerosol layers at depth in Uranus's atmosphere, we attempted to fit the combined HST/STIS and IRTF/SpeX data with three aerosol layers: one based at p > 5–7 bar; one based at 1–2 bar; and an additional extended haze in the upper troposphere/lower stratosphere. We allowed the mean particle sizes of these constituents to vary and also fitted for the imaginary refractive index spectra of each particle type. We found this to be a massively degenerate problem, with multiple combinations of parameters able to match the observations equally well. To attempt to prevent the solutions becoming stuck in local minima, we performed multiple retrievals with randomly varied a priori values of the pressure levels of these aerosol components and also their imaginary and real refractive indices, but were still unable to find a reliable model that was clearly better than any of the others. An additional problem is that since the HST/STIS data sets do not include estimated errors, we needed a way to set appropriate error estimates on these data that were compatible with the published errors on the IRTF/SpeX observations and so provided appropriate weight between the two data sets. Most of all, though, we needed a way to decouple the vertical structure of the best‐fitting aerosol profile from the wavelength‐dependent scattering and absorption properties of the particles necessary to model the spectra over a wide wavelength range.

The method we settled on to do this was to split the HST/STIS and IRTF/SpeX spectra into a series of narrow wavelength‐range “snippets”, each covering a wavelength bin of width 0.1 μm (narrow, but in most cases covering a reasonably large range of methane absorption strengths) and stepped by 0.05 μm to achieve Nyquist sampling. Then for each spectral snippet, we adopted a continuous distribution of aerosols with r _mean  = 0.1 μm and σ = 0.05, assumed that n _real  = 1.4, and fitted for the imaginary refractive index, which we assumed to be constant across the moderately narrow wavelength bins chosen. For the snippet covering 0.8–0.9 μm, which includes CIA bands of H₂–H₂ and H₂–He, we also fitted for the deep CH₄ abundance, which we then kept fixed for all other snippets. Hence, this snippet had to be fitted first.

For the HST/STIS “snippets” we only fitted to the nadir spectrum to keep the retrievals fast, while for the IRTF/SpeX “snippets” we fitted to the central‐meridian line average as described earlier, modelling spectra at two zenith angles (0° and 42.47°) and combining using Equation 2. At this point we considered what errors would be appropriate to assume for the data. Looking at the quality of the fits, we found we were able to fit the HST/STIS reflectivity spectra to a precision of χ ²/n ∼ 1 if we assumed the errors to be equivalent to 1/50 of the peak reflectivity in each bin. This estimated error covers unknown systematic uncertainties from sources such as the Lucy‐Richardson spatial deconvolution applied to the HST/STIS data, the inhomogeneous aerosol structure, and the methane absorption k‐table parameters. While the HST/STIS data do not include error estimates, the IRTF/SpeX data do contain these, and so we added the systematic R/50 errors in quadrature with the measured/scaled errors. Our approach meant that: (a) we achieved a good balance between the STIS and SpeX data; (b) we did not overfit the peak radiance wavelengths; and (c) in regions of low reflectance for the SpeX data we did not attempt to fit noise.

The fitted aerosol profiles for each snippet were then combined together and plotted as a contour function of aerosol opacity/bar against the mean bin wavelength as shown in Figure 5. Also plotted on Figure 5 are the depths for two‐way nadir transmission of 0.2 and 0.8, respectively, for an aerosol‐free atmosphere and including or excluding Rayleigh and Raman scattering to give an idea of how deeply we can see at different wavelengths. We can draw a number of conclusions from Figure 5. First of all these “snippet” retrievals favour a vertically confined aerosol layer at ∼1.5–2.5 bar for all wavelengths that penetrate deeply, except 0.4–0.6 μm, where a much deeper aerosol layer at p > ∼5 bar is indicated. The apparent base of this deeper layer is, as in the low a priori example shown in Figure 4, caused by the retrieved profile tending back to a priori at depth. The highest opacity/bar values of these constituents are generally found between two‐way transmission levels of 0.2 and 0.8, where the retrieval has most sensitivity. It is likely that the deeper p > 5–7 bar aerosol layer is present at all wavelengths, but cannot be detected at longer wavelengths due to lower reflectivity, increased gas absorption, and increased opacity of the overlying 1.5–2.5‐bar aerosol layer. It should also be noted that the aerosol layer at ∼2 bar has a much lower integrated opacity than the deeper layer, since although both have similar peak opacity/bar values, the 2‐bar layer spans a smaller range of pressure. The narrow wavelength range of detectability makes it difficult to say much about the particle size of aerosols in the deepest layer. However, for the 1.5–2.5 bar layer the opacity drops slowly with wavelength, suggesting particle sizes of the order of 1 μm. Finally, it can be seen that at short wavelengths (i.e., 0.3–1.0 μm) there is general indication of an extended upper tropospheric/lower stratospheric haze, where the opacity/bar increases with altitude at pressures less than ∼1 bar, but where the opacity drops with increasing wavelength. The retrieved wavelength dependence of this constituent of ∼1/λ ⁴ suggests a very small particle size, while the increasing opacity/bar with altitude is consistent with microphysical modelling of haze production in Ice Giant atmospheres of Toledo et al. (2019); Toledo et al. (2020), who predict a haze distribution with a fractional scale height of 2–4.

3.6. Simultaneous Retrieval of Combined HST/STIS, IRTF/SpeX, and Gemini/NIFS Uranus Observations

Having established the general vertical aerosol structure of Uranus from the first‐pass retrieval, and the wavelength‐dependent absorption from the snippet analysis, we moved to a parameterized model for the final retrieval. The parameterized model for Uranus included three aerosol layers: (a) a deep cloud/haze layer, “Aerosol‐1” fixed with a base at 10 bar (i.e., well below where we are sensitive to), with variable fractional scale height; (b) a vertically thin cloud/haze at 1–2 bar, “Aerosol‐2”, modelled as a Gaussian distribution with fixed width of ∼0.1 scale heights; and (c) an extended tropospheric/stratospheric haze, “Aerosol‐3”, based at 1–2 bar with a fixed fractional scale height of 2.0. Since the Aerosol‐1 layer is effectively infinitely deep compared with our sensitivity functions, in this study we only quote its integrated opacity from space down to the level at which this sensitivity diminishes, which can be seen from Figure 4 to be ∼5 bar. Also, since the deep Aerosol‐1 layer is only of importance at short wavelengths we assumed it to be composed predominantly of very small particles and assumed a Gamma distribution of sizes with r _mean  = 0.05 μm and variance, σ = 0.05. We assumed the same particle size distribution for the tropospheric/stratospheric haze particles (Aerosol‐3) since the opacity of this component needs to fall as 1/λ ⁴. For the intermediate cloud/haze layer at 1–2 bar (Aerosol‐2) we tried a number of mean particle sizes of size ∼1 μm, with variance σ = 0.3. Given the degenerate nature of this problem, we first fitted to the HST/STIS data alone in the wavelength range 0.8–0.9 μm to constrain the deep methane abundance and pressure of the 1–2 bar Aerosol‐2 component, which we then used as the a priori in subsequent all‐wavelength retrievals. To do this, we generated a set of 30 retrievals with randomly chosen values of the Aerosol‐1 fractional scale heights (in range 0.05–0.15), randomly chosen Aerosol‐2 mean particle radii (in range 0.1–4.0 μm) and then fitted for the opacities of all three constituents, their mean imaginary refractive indices (assumed constant in 0.8–0.9 μm range), together with the deep methane abundance, setting the relative humidity above the condensation level to 100% and limiting the stratospheric mole fraction to not exceed 10⁻⁴. The real refractive index of all components was assumed to be 1.4 (other values for the real refractive index were tried, but were not found to more or less appropriate). Then, for the best‐fitting aerosol/methane solutions to the 0.8–0.9 μm region (achieving χ ²/n < (1), we fitted the entire combined STIS/SpeX/NIFS data set using these retrievals as the a priori, keeping the mean Aerosol‐2 particle size fixed, but fitting again for all other properties and also for the imaginary refractive index spectra from 0.2 to 1.9 μm, with step 0.1 μm and assuming n _real  = 1.4 at 0.8 μm. We again reconstructed n _real at other wavelengths using the Kramers‐Kronig relation and then generated cross‐section spectra, single‐scattering albedo spectra and phase function spectra from Mie theory, approximating the phase functions with double Henyey‐Greenstein functions.

For the HST/STIS spectra at zenith angles of 0° and 61.45° we modelled 500 points from 0.2 to 1.0 μm, with 300 points from 0.2 to 0.5 μm, spaced every 0.001 μm to fully capture the Raman‐scattering component, and the remaining 200 points equally spaced from 0.5 to 1.0 μm. Although the HST/STIS data do not extend below 0.3 μm, we still needed to model the spectra from 0.2 μm upwards in order to compute the effect of sunlight being scattered to the longer wavelengths observed by STIS. We generated synthetic “measured” radiances at λ < 0.3 μm having the same reflectivity as observed at 0.3 μm, but with very large errors in order that NEMESIS did not try to fit to them, but did include their Raman‐scattering effect on the measured wavelengths. For wavelengths longer than 0.3 μm, we set the reflectivity errors to be equivalent to 1/50 of the reflectivity of the nearest peak, but in addition prevented the reflectivity error from exceeding a value of 0.01 in order to force the model to fit the data well near ∼0.5 μm, where the observed reflectances peak. We also halved the errors in the 0.8–0.9 μm range to ensure good fit at the wavelengths best able to differentiate between methane abundance and 1–2‐bar aerosol layer (Aerosol‐2) pressure.

For the Gemini/NIFS data, we selected 100 wavelengths equally spaced between 1.5 and 1.8 μm, again at 0° and 61.45° zenith angles. The errors were set to 1/100 of the maximum reflectivity of the 1.55 μm peak to give these data slightly more weight in the final fit than the IRTF/SpeX data since they contain clearer limb‐darkening information.

Finally, for the IRTF/SpeX data we selected 200 wavelengths, equally spaced between 0.8 and 1.9 μm (longer wavelengths were discounted as they can be seen in Figure 1 to have low SNR), and set the errors to be 1/50 of the nearest peak reflectance, added in quadrature with the published errors.

Figure 6 shows one of our best fits to the combined data set where we can see that we match the observed reflectivity spectra very well. The vertical cloud/haze structure retrieved in this case is shown in Figure 7, both in terms of opacity/bar and opacity/km at 0.8 μm (see Appendix C for an explanation of the different opacity units), which is comprised of a deep aerosol layer (Aerosol‐1) based at 10 bar with fractional scale height 0.11 and integrated opacity to 5 bar (at 0.8 μm) 0.81, a middle vertically thin aerosol layer (Aerosol‐2) at 1.48 bar and opacity 3.1, composed of particles of mean radius 0.6 μm, and an extended haze (Aerosol‐3), based at 1.6 bar with fixed fractional scale height 2.0 and retrieved opacity 0.03. The retrieved methane mole fraction is 3.0% and the relative humidity above the condensation level was fixed to 100%. It should be stressed here that there is no single solution that clearly fits better than all others. As we have performed multiple slightly perturbed retrievals there is a general class of solutions that fit best, of which this is one representative example. The complete set of retrieved parameters is shown in Figure 8 and the range of solutions and errors are listed in Table 2.

The estimated n _imag spectra for the three constituents are shown in Figure 9. Here, we have taken the fitted n _imag spectra for the best‐fitting cases and combined them together to give an idea of the range of acceptable solutions, together with two estimates of the optimal solution and the fitted spectra from our example best‐fitting case, which are listed in Table 3. As can be seen there are certain features common to all the best‐fitting retrieved imaginary refractive index spectra, namely:

1. The lowest aerosol layer (Aerosol‐1) has very low n _imag at 0.4–0.5 μm, which increases either side, but never gets very big. This means this layer is bright and highly scattering at 0.5 μm, but less so either side. Since we chose this layer to be composed of very small particles its contribution diminishes rapidly at longer wavelengths.

2. The 1–2‐bar aerosol layer (Aerosol‐2) has lowest n _imag from 0.4 to 0.6 μm, which then increases sharply on either side. This increase is necessary to give sufficient ultraviolet (UV) absorption and also to reproduce the observed limb‐darkening/limb‐brightening seen at longer wavelengths. The retrieved mean radius of the particles in this layer of ∼1 μm, makes them moderately forward‐scattering at 0.5 μm, lowering their back‐scatter and thus contribution to the reflected spectra, and increasing their “transparency” to allow light to be reflected from the deeper p > 5–7 bar Aerosol‐1 layer.

3. The tropospheric/stratospheric haze (Aerosol‐3) n _imag spectrum is moderately similar to that of the Aerosol‐2 layer, but is less well constrained.

It should be noted that the fitted n _imag spectra can alias all sorts of other deficiencies in the forward model, such as inaccuracies in the methane absorption data, assumed/estimated particle size, or vertical cloud/haze parameterisation, and so not all features in these spectra may necessarily reflect what we would measure were we able to make in situ observations of the particle scattering properties. However, the similarity between the retrieved Aerosol‐2 and Aerosol‐3 spectra and the generally higher values of n _imag leads us to posit that the middle Aerosol‐2 layer is not a pure condensed cloud, but instead is a region of photochemically produced haze transported from the upper atmosphere and perhaps mixed with some methane ice. This is a point we shall return to later. For the lower Aerosol‐1 layer, this has a lower n _imag and we posit that this is actually the H₂S condensation layer (rather than an NH₄SH layer, suggested by Sromovsky et al. (2019)), condensing onto dark haze particles and thus not conservatively scattering. We did attempt to model these data with a conservatively scattering cloud/haze, but found that there was a detectable contribution from this component at longer visible wavelengths, which from the “snippet” analysis we think is unlikely. Instead, this Aerosol‐1 layer is only visible from 0.4 to 0.7 μm due to: (a) this layer being most reflective at these wavelengths; (b) at longer (and shorter) wavelengths it is obscured by the thicker, and increasingly less scattering overlying 1–2‐bar Aerosol‐2 layer; and (c) we have assumed a small particle size to make sure the extinction cross‐section falls as 1/λ ⁴. Of course it is entirely possible that the particles in the Aerosol‐1 layer are actually larger than the value of 0.05 μm we have assumed here, in which case the retrieval model would likely have estimated larger n _imag at longer wavelengths in order to reduce the reflectivity to a sufficiently small values. Typical sizes of droplets in clouds on Earth are many orders of magnitude larger than value of 0.05 μm assumed here, although such clouds are formed of different condensate at very different pressures and temperatures so may not be representative. However, it is worth emphasising that the actual mean radius of the Aerosol‐1 particles is very uncertain and is not unambiguously constrained by these data.

One final point to make here is that the computed opacity/km of the Aerosol‐2 layer at 0.8 μm is rather low, peaking at about 1.0–10 km⁻¹. This underlines our suggestion that, assuming these are layers are homogeneous, they are extended photochemically produced haze layers, rather than opaque cloud decks, since we would be able to see many hundreds of metres in all directions were we to be present in these layers.

3.7. Analysis of IRTF/SpeX Neptune H + K Band Observations

Having found a good fit to our Uranus observations, we then turned our attention to the Neptune observations. Before analysing the combined HST/STIS, IRTF/SpeX and Gemini/NIFS data sets we first of all sought to understand the gross differences between the Uranus and Neptune IRTF/SpeX observations, in particular the high reflectivities seen at methane absorption wavelengths longer than 1 μm for Neptune, but not seen for Uranus.

The 1.65–2.4 μm region is particularly good for determining upper atmospheric aerosol density since the strong methane absorption band at 1.7 μm lies next to strong H₂‐H₂/H₂‐He collision‐induced absorption bands centred at 2.1 μm. The weighting functions of the methane‐absorbing and hydrogen/helium‐CIA‐absorbing regions cover a similar range of pressure levels in the upper troposphere (0.6–0.1 bar) and thus this spectral region can be used to resolve the degeneracy between aerosol opacity and upper tropospheric methane abundance and retrieve reliable estimates of both (Roman et al., 2018).

Figure 10 shows a “bracketed” retrieval to the IRTF/SpeX spectrum in the wavelength range 1.65–2.4 μm, assuming two widely separated continuous distributions of small aerosols of mean radius 0.1 μm, variance σ = 0.3, and fixed n _imag  = 0.1. As can be seen, a sharply peaked aerosol profile is favoured, peaking at ∼0.2 bar, with the local methane relative humidity retrieved to be ∼32%. At first glance this fitted aerosol structure seems contrary the results of microphysical modelling (Toledo et al., 2019, 2020) where, since photochemically produced haze particles are produced in the stratosphere and then progressively mixed to deeper and deeper pressures, we expect a distribution with decreasing opacity/bar with increasing pressure. To understand what was going on here, we needed to examine the reflectivities over a wider wavelength range, which, as we did for Uranus, we achieved with a Neptune “snippet” analysis, which we outline in the next section.

3.8. Analysis of Combined HST/STIS and IRTF/SpeX Neptune Observations ‐ Snippet Analysis

As for our Uranus analysis, we split the HST/STIS and IRTF/SpeX observations into small, manageable “snippets” of wavelength width 0.1 μm, spaced every 0.05 μm and retrieved the vertical aerosol profile for a Gamma distribution of particles with mean radius 0.1 μm and variance σ = 0.05, varying also the imaginary refractive index n _imag and keeping n _real  = 1.4. The resulting retrieved aerosol vertical/wavelength structure is shown in Figure 11. As can be seen, the aerosol structure shares many similarities with the equivalent Uranus results (Figure 5), with a deep Aerosol‐1 layer seen at p > ∼4 bar and a middle Aerosol‐2 layer at ∼2 bar (slightly deeper than for Uranus). At short wavelengths it can also be seen that the opacity/bar increases with height in the upper troposphere/lower stratosphere, and decays rapidly with wavelength, consistent with our expectations with small, photochemically produced haze particles. However, for Neptune there is clearly also a component of vertically confined particles centred at ∼0.2 bar, which must have a larger mean size in order that their opacity does not drop noticeably with wavelength. We believe here that we are detecting in the IRTF/SpeX data the signature of an additional component of larger‐sized particles at ∼0.2 bar, which we surmise to be condensed methane ice. Such upper tropospheric methane condensation clouds are often seen in Neptune observations, especially in the regions at 20–40°N and 20–40°S and are seen to be highly temporally and spatially variable. We have tried to mask out such regions in the HST/STIS and Gemini/NIFS data. However, since we do not have an image of Neptune observed simultaneously with the IRTF/SpeX central‐meridian line average, we do not know how significant this component of upper tropospheric methane clouds is in the recorded IRTF/SpeX spectrum. Hence, we need to be careful when interpreting the IRTF/SpeX Neptune data, especially when used in combination with the other data sets. However, the good correspondence between the scaled IRTF/SpeX and Gemini/NIFS data sets, achieved, as described earlier, by optimising the cut‐off level for bright clouds in the Gemini/NIFS data, reassures us that these data are reasonably self‐consistent.

3.9. Simultaneous Retrieval of Combined HST/STIS, IRTF/SpeX, and Gemini/NIFS Neptune Observations

Having found from the “snippet” analysis that our expected aerosol solution for Neptune was very similar to that for Uranus, we attempted to fit the combined STIS/SpeX/NIFS Neptune data set in the same way, except that we also added a thin layer of methane ice centred at 0.2 bar and with variable mean radius and fixed variance σ = 0.3 (N.B. although we have modelled this layer as being homogeneous, it is also possible that this component is composed of sub‐pixel‐scale clouds). We first fitted the 0.8–0.9 μm region separately with a range of randomly varied a priori Aerosol‐1 fractional scale heights and Aerosol‐2 base pressures and radii and used those achieving χ ²/n < 1 as the a priori for the full retrievals, fixing the Aerosol‐2 radii. Figure 12 shows one of our best fits to the combined Neptune data set, where we can see that we again match the observed reflectivities and limb‐darkening very well. The vertical aerosol structure retrieved for this case is shown in Figure 13, which is comprised of a deep Aerosol‐1 layer at fixed base pressure 10 bar with fractional scale height 0.18 and integrated opacity to 5 bar (at 0.8 μm) 0.81 (coincidentally the same value as for our sampled Uranus retrieval shown earlier), a compact Aerosol‐2 layer at 2.08 bar and opacity 1.8 composed of particles of mean radius 0.68 μm, an extended haze (Aerosol‐3), based at 1.6 bar with fixed fractional scale height 2.0 and opacity 0.05, and a detached methane ice aerosol layer near the tropopause (0.2 bar), of opacity 0.03 comprised of particles of mean radius 2.8 μm. The methane deep mole fraction is 7.7% and relative humidity at the tropopause is 34% in this case. It should be stressed here that once again there is no single solution that fits better than all others; as we have performed multiple slightly perturbed retrievals there is a general class of solutions that fit best, of which this is one representative example. The complete set of retrieved parameters is shown in Figure 14 and the range and errors on our best fit solutions are listed in Table 2. It is worth noting, however, that the retrieved imaginary refractive index spectra of the main aerosol components, shown in Figure 9, and listed in Table 4 are similar to those retrieved for Uranus, suggesting a similar composition. It is also worth noting that the opacity of the middle Aerosol‐2 layer at ∼2 bar is approximately half that determined for Uranus's atmosphere, while the opacity of the deeper Aerosol‐1 layer is also greatly reduced. The lower Aerosol‐2 opacity means that the Aerosol‐1 layer has a greater contribution to the modelled spectra and so we have greater constraint on the Aerosol‐1 imaginary refractive index spectrum as can be seen. For the methane ice layer at 0.2 bar, this is retrieved to be composed of moderately large particles (r = 2–3 μm), which means they are most reflective at the longer wavelengths where the increased scattering at methane‐absorbing wavelengths is seen. It also means that they are effectively conservatively forward‐scattering at visible wavelengths, and thus they have very little effect here (as can be seen later in Figure 18). Finally, since we discriminated against cloudy regions when we compiled our mean data sets, the opacity we retrieve for the Aerosol‐4 layer will of course be significantly less than a true disc‐average that includes such clouds.

One curious feature of these retrievals is the significantly larger deep methane abundance we retrieve for Neptune of 7 ± 1%, versus 3 ± 1% for Uranus. An analysis of recent Neptune/MUSE observations (Irwin et al., 2021) found values varying from 5 ± 1% at the equator to 3 ± 1% at the south pole, which compared reasonably well with previous analyses of these HST/STIS observations and others (e.g., Karkoschka & Tomasko, 2009; Karkoschka & Tomasko, 2011; Sromovsky et al., 2011, 2014) who found values of ∼4% at equatorial latitudes and ∼2% at polar latitudes for both planets. In this study, we have analysed disc‐averaged spectra so we might have expected to retrieve a methane abundance of ∼3%–4% for both Uranus and Neptune, and given that we would see more of the polar regions in the Neptune disc‐averages than in the Uranus observations we might expect to retrieve lower methane abundances for Neptune than we do for Uranus. In practice the determination of “deep” methane abundance (i.e., immediately below the methane condensation level) is almost inextricably tied up with the assumed cloud/haze parameterisation scheme, and also the assumed vertical distribution of methane. Our simple haze model, which is designed to match the observations at all wavelengths simultaneously, leads to fits that do not match the 800–860 nm range as well as other models that are able concentrate on this region alone (e.g., Irwin et al., 2021) and achieve closer fits, lower methane abundances, and slightly different aerosol profiles. Also, in this work we have assumed the same simple “step” model for the methane profile (constant mole fraction up to some fraction of the saturated vapour pressure curve set by the relative humidity) as Irwin et al. (2021), while other authors used a “descended” methane profile (e.g., Sromovsky et al., 2019). Finally, although there is little cross‐correlation seen in our “corner‐plots” (Figures 8 and 14) it is possible that the methane abundance is correlated with some of the cloud parameters, perhaps Aerosol‐1, and hence for Neptune the methane abundance is being slightly overestimated and the Aerosol‐1 parameters under‐ or over‐estimated. Hence, in practice the methane retrieval problem is extremely degenerate and while here we find a difference in the retrieved deep methane abundances for Uranus and Neptune, there is not sufficient evidence to claim that this is a robust result. Perhaps with a revised haze parameterisation scheme, or different Aerosol‐1 particle assumptions, we may determine values that are closer to each other. We hope to return to this question in future work.

3.10. HST/WFC3 Observations of Dark Spots in the Atmospheres of Uranus and Neptune and Their Spectral Characteristics

The analysis of available reflectance spectra of Uranus and Neptune reveals very similar aerosol structures for both planets, with a deep Aerosol‐1 layer based at p > 5–7 bar detectable at visible wavelengths, a middle Aerosol‐2 layer at 1–2 bar seen at other wavelengths, a vertically extended stratospheric haze (Aerosol‐3), detectable mostly at visible wavelengths, and finally, for Neptune, the presence of condensation layer of moderately large methane ice particles near the tropopause. Such a model can then be used simulate the general appearance at all wavelengths of either planet with good accuracy. However, there is one aspect in which Uranus and Neptune seem to differ substantially and that is the presence of “dark spots.”

The first dark spots observed in an Ice Giant atmosphere were the Great Dark Spot (GDS) and Dark Spot 2 (DS2), seen by Voyager 2 in Neptune's atmosphere in 1989 (Smith et al., 1989; Sromovsky et al., 1993). The GDS was most visible at a wavelength of 0.48 μm, but was undetectable longward of 0.7 μm. Since the Voyager flyby several new dark spots have been seen in Neptune's atmosphere from Hubble Space Telescope observations in 1994, 1996, 2015 and 2018 (Hammel et al., 1995; Hsu et al., 2019; Sromovsky et al., 2001; Wong et al., 2018). These have all had similar spectral characteristics to those noted for the GDS, namely that they are visible near 0.5 μm, but not visible longward of 0.7 μm. For Uranus, there has to date only been one dark spot detected, which was observed in 2006 (Hammel et al., 2009). Unlike the Neptune dark spots the Uranus dark spot seems to have been detectable to longer wavelengths.

The dark spots seem to be some sort of vortex structure and the colouration must be due to either a darkening of the aerosol particles, or a change in opacity of the aerosol layer, or perhaps a combination of both effects. Putting aside the Uranus dark spot observation for one moment, the darkening effect needs to be confined to a rather small and particular range of wavelengths. It could very well be caused by a spectrally limited darkening of the aerosol particles, but we would then need to have quite a particular perturbation of the cloud/haze scattering properties and simultaneously confine that change somehow to be within a vortex. On the other hand, if dark spots are caused by cloud/haze opacity changes, we need to explain how its signature is so spectrally confined. We believe that the retrievals described in previous sections provide the answer. The deep Aerosol‐1 layer based at p > 5–7 bar is only visible at wavelengths less than ∼0.7 μm and if that were to have lower opacity, or reduced reflectivity, then we might expect to see a dark spot at just those wavelengths.

To see if this phenomenon was apparent in our modelled Uranus observations we generated synthetic images of Uranus from our best‐fitting cloud/haze model in the seven wavelength channels observed by HST/WFC3 with the OPAL program in 2014 (Simon et al., 2015). The location of the filters used (F467M, F547M, FQ619N, F658N, FQ727N, F845M, and FQ924N) are compared with the central‐meridian line‐averaged HST/STIS Uranus spectrum in Figure 15 and in Figure 16 we compare the observations with synthetic images. The middle and bottom rows of Figure 16 show synthetic images at each filter wavelength, which were constructed using the following procedure:

1. In our baseline simulations we used our optimal fit to the reconstructed spectra at 0° and 61.45° to extract Minnaert parameters I ₀(λ) and k(λ) and, using $I=I_0{\mu }_0^k{\mu }^{k-1}$, constructed synthetic disc images of Uranus at each modelled HST/STIS wavelength, I _reference (λ) for the HST/WFC3 observation zenith angles μ and μ ₀.

2. Then, from our fitted cloud/haze model, we calculated modelled spectra (at 0° and 61.45°) with either: A) the opacity of the Aerosol‐1 layer based at p > 5–7 bar set to zero; or B) the imaginary refractive indices of the particles in this layer set to 0.001 at all wavelengths (which darkens the particles at short wavelengths). We then extracted modified R ₀(λ) and k(λ) Minnaert parameters for both cases, which we used to generate two additional sets of synthetic Uranus images, I _modified (λ).

3. We then constructed a weighting factor, f _w , to simulate a cloud/haze darkening/clearing near the disc centre, with latitude Φ ₀ = 20°N and central meridian longitude Λ ₀ = 0°, setting $f_w=\exp \left({\left(\left(\mathrm{\Phi }-{\mathrm{\Phi }}_0\right)/\mathrm{\Delta }\mathrm{\Phi }\right)}^2+{\left(\left(\mathrm{\Lambda }-{\mathrm{\Lambda }}_0\right)/\mathrm{\Delta }\mathrm{\Lambda }\right)}^2\right)$, where Φ is the latitude on the disc, Λ is the longitude relative to the central meridian, and where ΔΦ = 5° and ΔΛ = 10°.

4. We then made weighted averages of our two sets of images at each HST/STIS wavelength: I _mean (λ) = (1 − f _w )I _reference (λ) + f _w I _modified (λ).

5. Finally, the combined images were convolved with the HST/WFC3 filter functions to create the synthetic images.

From Figure 16 it can be seen that aside from the latitudinal cloud/haze and methane abundance variations visible in the observed images (which we did not attempt to simulate here) we capture reasonably well the observed limb‐darkening/limb‐brightening at all wavelengths. Although the HST/WFC3 OPAL 2014 observations did not contain any dark spots, if there were a hole in the deep Aerosol‐1 layer in Uranus's atmosphere, or a darkening, it can be seen in the middle and bottom rows of Figure 16 that it would produce a dark spot at precisely the same wavelengths as seen by HST/WFC3, that is, at 467 and 547 nm. However, while for the Aerosol‐1‐darkening case the spot is invisible at longer wavelengths, for the clearing case it is still just visible at 658, 845 and 924 nm. In addition, the darkening simulation results in a spot that is darker at 467 nm than in the clearing simulation, which is again more consistent with the observations.

Turning to Neptune, we generated synthetic images from our best‐fitting model in the seven wavelength channels observed by HST/WFC3 with the OPAL program in 2018 (Simon et al., 2019), which reported the detection of a new dark spot, NDS‐2018. The channel filter functions used (F467M, F547M, FQ619N, F657M, FQ727N, F763M and F845M) are shown in Figure 15 and the observed images can be seen in the top row of Figure 17. In the middle and bottom rows, we present our simulations with a dark spot centred at Φ ₀ = 15°N, Λ ₀ = 0°, and also with a clearing/darkening at latitude Φ ₀ = 60°S. Here, we can again see that we predict the appearance of the clearing of the deep Aerosol‐1 layer to be most clearly detectable at the two shortest wavelengths, very similar to the observed spot. However, for the clearing simulation, the modelled spot is still just visible at 657, 763 and 845 nm, which is contrary to the observations and also contrary to the darkening simulation, where the modelled spot is completely invisible at the longer wavelengths. In addition, the spot modelled with the darkening hypothesis is darker at 467 nm than the clearing case, which is again more consistent with the observed properties.

Considering Uranus and Neptune together it can be seen that although a clearing of the Aerosol‐1 layer based at p > 5–7 bar produces almost the right response, it does not predict the spot to be darker at 467 nm than 547 nm, and the simulated spot can still just be seen at longer continuum wavelengths, contrary to the observations. In contrast, our simulated images where we darken the Aerosol‐1 layer produces a feature that much more closely resembles the real NDS‐2018 on Neptune. We find that removing the Aerosol‐1 layer does not darken the spot sufficiently at short wavelengths since there is still significant Rayleigh scattering from the air itself at these depths. However, by darkening this aerosol layer, the contrast of the dark spot at short wavelengths is greatly increased. Add to that the disappearance of the feature at longer wavelengths by darkening the Aerosol‐1 particles, rather than removing them altogether, and on balance we conclude that darkening the particles provides a better match to the observations. Further evidence in support of this hypothesis comes from an analysis of earlier HST/WFC3 Neptune images (Karkoschka, 2011), where the author suggested that the dark belt at 30–60°S was best explained by a darkening of the haze at pressures >3 bar. Karkoschka (2011) went on to speculate that since dark spots have very similar spectral characteristics to the 30–60°S belt, they might be explained by vertical motions bringing dark material from depth up to altitudes where it can be observed. Finally, we also see that our predicted Neptune dark spot is slightly darker than the equivalent spot modelled for Uranus, which can be explained by the fact that the opacity of the overlying ∼2‐bar Aerosol‐2 layer on Neptune is lower than it is on Uranus. The higher opacity of the UV‐absorbing Aerosol‐2 layer on Uranus also explains why the UV reflectivity of Uranus is lower than Neptune.

Since the scattering cross‐section of the approximately micron‐sized Aerosol‐2 particles is found to have a roughly white visible reflectivity spectrum, the higher opacity of the Aerosol‐2 layer on Uranus also mostly explains why Uranus appears to have a paler blue colour to the human eye (to some more greenish) than Neptune, which we show in Figure 18. In this figure we show the simulated colours of the background, masked regions of Uranus and Neptune reconstructed from the Minnaert limb‐darkening approximations to the HST/STIS observations using Equation 1. Here, as a representative case we assumed the solar and emission zenith angles varied across the disc of both Uranus and Neptune as observed for the Uranus 2014 HST/WFC3 observations shown in Figure 16. In these simulations the modelled radiance spectra across the planet were convolved with the CIE‐standard (Commission Internationale de l’Éclairage) red, green, and blue human cone spectral sensitivities of Stockman and Sharpe (2000); Stockman (2019) and are, as near as we can determine the “natural” colours that the background atmospheres of these planets would have if we were to observe them from space. As can be seen Uranus appears paler and less blue than Neptune, although the colour differences are actually rather subtle. In Figure 18 we can see that if we modelled the planets' appearances for aerosol‐free conditions and with no methane absorption then both planets would appear pearly white since even though Rayleigh‐scattering is more effective at blue wavelengths, we eventually reach depths at even red wavelengths where Rayleigh scattering becomes effective. It is only when we add absorption by the red‐absorbing methane that the characteristic blue colour of these planets develop. Then as we add the scattering effects of Aerosol‐1, Aerosol‐3 and Aerosol‐4 components the modelled appearances become paler, but it is not until we add in scattering from Aerosol‐2 that the main colour difference manifests itself. Why the Aerosol‐2 layer on Uranus is thicker than that of Neptune is not clear. Potentially, we suggest that perhaps the more dynamically overturning atmosphere of Neptune is more efficient at clearing this haze layer through methane condensation at its base, as discussed later. Also, we note that during the Voyager‐2 encounters with these planets, Uranus was close to southern summer solstice when the south polar region was covered in a “hood’ of haze, which we presume to be Aerosol‐2. Hence, the planet would likely have looked even paler than we have simulated here from HST/STIS observations in 2002 and 2003.

In summary we find that the aerosol structure that we have retrieved from the combined HST/STIS, IRTF/SpeX and Gemini/NIFS data for Uranus and Neptune are also consistent with HST/WFC3 observations and can also explain the observed colour difference between these planets. We also find that a darkening (or possibly opacity change) of the deepest aerosol layer (Aerosol‐1) leads to a darkening that is consistent with the observed spectral characteristics of dark spots. Hence, we propose that the dark spots sometimes seen in Neptune's atmosphere (and occasionally Uranus's) are caused by a darkening of the deep Aerosol‐1 layer. If dark spots in Uranus's atmosphere really are visible over a wider wavelength range, as reported by Hammel et al. (2009), then this suggests that Uranus's deep Aerosol‐1 layer, presumably containing a significant fraction of H₂S ice, is brighter over a wider range of wavelengths, since retrieval tests conducted with a conservatively scattering lower cloud/haze produced a detectable dark spot signature at all the HST/WFC3 wavelengths when the opacity of this layer was reduced, except in the strong methane bands at 619 and 727 nm.

3.11. Voyager‐2/ISS Observations of Dark Spots in Neptune's Atmosphere

Having found a good correspondence with the spectral properties of dark spots seen in Neptune's atmosphere with HST/WFC3 by darkening the Aerosol‐1 layer, we wondered whether we might also be able to find a good correspondence with the original Voyager‐2 imaging observations, which first detected such phenomena (Smith et al., 1989). During Voyager‐2's approach to the Neptune system in August 1989, the Imaging Science Subsystem (ISS) took a particular sequence of observations, “Full color set for atmospheric dynamics,” consisting of 319 full‐disc observations from 16th to 18th August in each of the main filters of its Narrow‐Angle Camera: UV, Violet, Blue, Green, Orange and also with no spectral filter (Clear). Observations were also made in three filters of its Wide‐Angle Camera: “CH4‐U” (541 nm), “CH4‐JS” (619 nm) and another Orange filter. The transmission functions of all these filters are shown in Figure 15. In these images the Great Dark Spot (GDS) and also Dark Spot 2 (DS2) were observed at several phase angles. We examined the limb‐darkening properties of the GDS and also the 50°S dark belt in these observations to see if they supported or contradicted our hypothesis that dark regions are caused by a clearing, or darkening, of the deep H₂S/haze Aerosol‐1 layer. Representative Voyager‐2 observations from this set are shown in Figure 19. It can be seen that the GDS is particularly prominent in the Blue filter channel, but becomes less so as we move to both shorter and longer wavelengths. The bright companion clouds become more prominent at longer wavelengths, as the Rayleigh scattering becomes less significant, indicating these features to lie at higher altitudes, presumably in the upper troposphere. The contrast between the dark spots and the surrounding regions looks similar to that between the darker belt centred at 50°S and surrounding latitudes and we assume, like Karkoschka (2011), that both are formed a similar way. Looking at the 50°S dark belt it can be seen that this region is most visible near the central meridian and gets less clear towards the limbs.

To interpret these observations more quantitatively, we performed a limb‐darkening analysis of the Voyager‐2 ISS images, which we summarise in Figure 20. In this figure we have analysed the observations for each of the six NAC flters studied, first of all concentrating on images where the GDS and DS2 were not visible and analysing the limb‐darkening of the background atmosphere in two latitude bands: 15–25°S (centred on the GDS) and 45–55°S (centred on the dark belt). Assuming the Minnaert approximation to hold, Equation 1 can be re‐expressed as

where R is the measured reflectance (i.e., I/F), μ and μ ₀ are cosines of the viewing and solar zenith angles respectively, and R ₀ and k are the fitted Minnaert parameters. Hence, plotting μR against μμ ₀ on a log‐log plot should yield a straight line if the Minnaert approximation is good and in Figure 20 we have plotted such curves for the two latitude belts for all filters. We find that the observed limb darkening is well described by the Minnaert model and also that the 45–55°S belt has a lower (I/F)₀ and is noticeably less limb‐darkened than the 15–25°S belt, having a significantly shallower slope (i.e., lower k, but still greater than 0.5, which is the boundary between limb‐darkening and limb‐brightening).

We then analysed individual observations in the data set containing the GDS at different zenith angles, selecting pixels within the GDS and over‐plotting their reflectance in the same way in Figure 20. There was good sampling for the NAC images, but far fewer WAC images, which were also less spatially resolved making the limb‐darkening less clear. In addition, the CH4‐JS channel seems poorly flat‐fielded and photometrically corrected. Hence, the WAC images were discounted in the limb‐darkening analysis. It can be seen in Figure 20 that pixels within the GDS have a similar (I/F)₀ to those in the 50°S belt and the limb‐darkening is again less than the background 20°S region, but less so than for 50°S. However, the fact that the GDS and 50°S belt have both lower (I/F)₀ and lower k than the background at 20°S suggests these dark features really do have a common origin. To determine if this common origin is due to lower reflectance from the deep Aerosol‐1 layer we have also plotted in Figure 20 the simulated limb‐darkening from our best‐fitting cloud/haze retrieval to the combined STIS/SpeX/NIFS Neptune data set, including the contribution from this layer, or modifying it by clearing or darkening. It can clearly be seen that reducing the contribution from the Aerosol‐1 layer lowers both (I/F)₀ and the limb‐darkening parameter k, just as is seen for the Voyager‐2 dark regions, suggesting that this may indeed be the cause of the darkening. On balance, the difference in limb‐darkening is again slightly better simulated by darkening the Aerosol‐1 layer, rather than removing it.

In summary, the Voyager‐2 ISS observations of Neptune support our suggestion that dark regions are caused by a darkening (or possibly clearing) of the deep Aerosol‐1 H₂S/photochemical haze layer. It has sometimes been suggested that the dark spots might be caused by a darkening of the overlying haze. We looked to see if darkening the 1–2‐bar Aerosol‐2 layer could explain the observations, but found that this is not supported by the limb‐darkening evidence. Such a scenario would lead to stronger limb darkening than observed, since as the spot moved towards the limb we would be looking through more and more of this darker haze. Instead, we believe we are looking at the darker deep Aerosol‐1 layer through the fairly uniform 1–2 bar Aerosol‐2 layer, which is reasonably scattering at these wavelengths. Hence, as spots move towards the limb we see more reflectance from the Aerosol‐2 layer at 1–2 bar and less from the deep Aerosol‐1 layer, leading to the reduced limb darkening observed.

4.1. Vertical Aerosol Structure

Excepting the additional upper tropospheric methane aerosol layer on Neptune, we find very similar vertical aerosol structures for both Uranus and Neptune composed of:

1. Aerosol‐1: A deep layer based at p > 5–7 bar of moderately scattering particles (assumed, as previously discussed, to be sub‐micron‐sized) that we assume to be coincident with the main H₂S cloud/haze condensation layer, perhaps composed of a mixture of photochemically produced haze particles and H₂S ice;

2. Aerosol‐2: A vertically thin haze near the methane condensation level at 1–2 bar, composed of approximately micron‐sized particles that are scattering at visible wavelengths, but more absorbing at UV and longer wavelengths; and

3. Aerosol‐3: A vertically extended haze of sub‐micron‐sized particles with a fractional scale height of 2 and moderately similar refractive index spectra to the 1–2 bar Aerosol‐2 layer.

The similarity of the structure (both vertical and spectral scattering) of the two planets should not be surprising given their similar tropospheric gaseous composition and temperatures, but how can we account for the fact that we apparently have a layer of photochemically produced haze at the methane condensation level and not methane clouds? We hypothesise that to account for this structure we need to consider the vertical stability of Uranus's and Neptune's atmospheres.

In Appendix A we outline the theory behind buoyancy and static stability in planetary atmospheres. At pressures where the temperature profile falls more slowly with height than the lapse rate, the air is stable to vertical perturbations and will try, if moved vertically, to return to its original level, giving rise to gravity waves with a frequency defined by the Brünt‐Väisälä frequency, N. For most atmospheres, this frequency just depends on the vertical gradient of the potential temperature $\theta =T{\left(\frac{p_0}{p}\right)}^{\gamma }$, where γ = R/C _p , T and p are the local temperature and pressure, p ₀ is a reference pressure level, R is the gas constant and C _p is the molar heat capacity at constant pressure. The Brünt‐Väisälä frequency can be shown in Appendix A to be:

where θ ₀(z) is the potential temperature of the background air at altitude z.

The latent heat released by the condensation of clouds lowers the rate at which temperature falls with height, making the air more stable for both planets, and this can be seen in our calculations of the Brünt‐Väisälä frequency in Figure 21, with a local region of higher Brünt‐Väisälä frequency seen just above the methane condensation level. However, while Equation 4 is perfectly acceptable for most atmospheres, where the molecular weight does not change greatly with height, this is not the case for Uranus and Neptune near the methane condensation layer. Assuming a deep mole fraction of 4% for methane, the mean molecular weight reduces from ∼2.9 to ∼2.3 g/mol for both planets across the methane condensation region, a decrease of ∼20%. This decrease in molecular weight at the condensation level increases the static stability even further since the heavier air below is naturally less buoyant than the lighter air above. We show in Appendix A that a more appropriate equation for the Brünt‐Väisälä frequency for atmospheres where the molecular weight also changes significantly with height is:

where M ₀(z) is the molecular weight of the background air at altitude z. It can be seen in Figure 21 that including this term greatly increases the static stability of the atmosphere at the methane condensation level.

We hypothesise that this localised vertical region of static stability somehow leads to an enhanced opacity of larger haze particles that, at the lower boundary, act as seed particles for methane condensation. Given the large mole fraction of methane, and ready supply of cloud condensation nuclei (CCN), methane cloud formation will likely be extremely rapid and Carlson et al. (1988) estimate that methane particles at 1–2 bar in Ice Giant atmospheres may grow to a size of ∼5 mm in as little as ∼100s. Hence, these large methane ice particles would almost immediately precipitate, or “snow”, out, which would explain why pure methane ice clouds are almost never seen at these levels. The precipitating methane ice particles will drop to deeper levels before sublimating and releasing their payload of photochemical haze particles, which in turn may serve as CCNs for H₂S ice formation in the 4–10‐bar region. The process is analogous to the “smust” concept that Hunten (2008) proposed to explain a vertical gap in the heavy hydrocarbon gas abundances of Jupiter. The difference is that in the case of Uranus and Neptune, methane ice transports haze particulates, while in the case of Jupiter, Hunten (2008) proposed that haze particulates transported volatile hydrocarbon species. Since the mole fraction of H₂S is expected to be very much lower (∼10⁻⁴) than for methane, particle growth at 4–10‐bar will not be so rapid (e.g., Carlson et al., 1988) and hence a stable layer of H₂S cloud could form. In addition, there are no vertical stability barriers since the latent heat release and change in mean molecular weight at this pressure level are very small, as can be seen in Figure 21.

Returning to the Aerosol‐2 layer, methane condensation acts to create a clearing beneath it in two ways. Condensation onto haze particles at the base of the layer causes them to snow out directly, and the stable layer produced by latent heating and reduced molecular weight reduces K _zz and slows the rate at which haze particles in the Aerosol‐2 layer mix into the deeper levels. Within the layer, the mechanism for the local enhancement of haze opacity is currently unknown and understanding of this newly discovered feature will require a future microphysical modelling study. Toward the top of the compact Aerosol‐2 layer, the decrease in haze density with altitude (in the presence of eddy mixing) implies some sort of particle density source at the base of the layer, rather than above it, because eddy mixing alone can only modify the slope of a gradient from the source region to the sink region, not reverse it. The action of methane, near saturation, could be significant. Considering the likely hydrocarbon composition of haze particles, it is possible that methane condensation, even at subsaturated vapour pressures, could influence the surface chemistry of haze particles (in analogy to hydrophilic aerosols significant in terrestrial hazes). It is also possible that while methane condensation at the base is rapid, leading to the formation of large particles, there will also be some component of smaller particles in the size distribution. The fractional abundance of these smaller haze/ice particles is likely to get larger as we move up to lower pressures and lower temperatures in the Aerosol‐2 layer, where the saturated vapour pressure of methane will be much lower. Hence, what we may be seeing in the Aerosol‐2 layer is a mixture of photochemical haze and methane ice. Microphysics in the unique environment of the Aerosol‐2 layer may thus be complex, and it is clear that the haze particles cannot be considered static tracers of mixing.

While this scenario provides an attractive model for the background aerosol structure of both planets, the atmosphere of Neptune is more energetic and contains frequent cloud features that appear to be methane condensation. For the most part these appear to form in the upwelling regions at 20–40°N and 20–40°S and at pressures of 0.6 to 0.1 bar (e.g., Irwin et al., 2016; Molter et al., 2019). However, Irwin et al. (2011a) report the appearance of several clouds at 60–70°S in 2009 that are deep in the atmosphere at 1–2 bar and may be rare examples of methane clouds forming near the methane condensation level itself.

As to the nature of the photochemical haze, there are several candidates. Khare et al. (1993) report the spectra of “tholins” generated by irradiating H₂O/C₂H₆ mixtures, which have minimal absorption at 0.8 μm and increasing absorption as we move to UV, or longer wavelengths. While the spectrum of these “tholins” is not identical to those derived here it is qualitatively similar. Another candidate is acetylene soot (Dalzell & Sarofim, 1969), which also has increasing absorption at long wavelengths. Indeed, such hazes share some characteristics with the phenomenon of “blue hazes” commonly observed in the Smoky Mountains of eastern Tennessee and the Blue Ridge Mountains of Virginia (Ferman et al., 1981). Examining the scattering cross‐sections of the retrieved aerosol particles, we find that the small particles assumed to be present in the deep Aerosol‐1 component, and found to be necessary in the extended Aerosol‐3 component above ∼1.5 bar, would appear to be blue at visible wavelengths. However, the larger particles in the 1–2 bar Aerosol‐2 layer would appear white in the visible.

The stability of the 1–2‐bar region in the Ice Giant atmospheres has also been concluded by Guillot (1995) and Leconte et al. (2017). Intriguingly, Teanby et al. (2020) noted that a stable layer at this pressure level has important implications for the internal oxygen enrichment of Neptune (Cavalié et al., 2017; Teanby et al., 2020). Limited mixing at 1–2 bar could permit externally sourced CO entering Neptune's atmosphere from comets to be trapped in the upper troposphere, which removes the need to have extreme internal oxygen enrichment in order to explain the wide wings in Neptune's sub‐mm CO lines. Previous work had suggested CO to be present throughout the troposphere, which requires an extremely high internal oxygen enrichment compared with the solar composition by a factor of at least a few hundred (Cavalié et al., 2017; Lellouch et al., 2005; Luszcz‐Cook & de Pater, 2013; Moses et al., 2020) in order to act as a deep source for tropospheric CO. If externally sourced CO can be kept in the upper troposphere by limited mixing this opens up the possibility of a more rock‐rich Uranus and Neptune (Teanby et al., 2020), in keeping with more modest oxygen enrichments of ∼50 inferred from the observed D/H ratio (Feuchtgruber et al., 1997).

Finally, in Appendix D we outline a review of methods for calculating the eddy diffusion coefficient, K _zz , in planetary atmospheres. With some modification of model parameters suggested by Ackerman and Marley (2001) we found we could derive a K _zz profile that was consistent with upper stratospheric determinations of K _zz by Fouchet et al. (2003), as can be seen in Figure 21. With these parameters we find a local minimum of K _zz at the methane condensation level, but also a smaller minimum near the tropopause; this second K _zz minimum may help to explain the detached methane ice layer (Aerosol‐4) seen near the tropopause of Neptune, where the atmosphere may be dynamically vigorous enough to mix methane up to such levels.

4.2. Nature of Dark Regions in the Atmospheres of the Ice Giants

In this paper we have shown that dark regions seen occasionally in Neptune's atmosphere, and rarely in Uranus's atmosphere, can be explained by a darkening, or opacity change of the deep Aerosol‐1 layer, which we have suggested may be composed of a mixture of H₂S ice and photochemically produced haze. Such a condensation level is consistent with the estimate from ground‐based microwave studies that the deep abundance of H₂S in Uranus's and Neptune's atmospheres are $53.8_{-13.4}^{+18.9}\times$ solar (Tollefson et al., 2021) and $37_{-6}^{+13}\times$ solar (Molter et al., 2021), respectively. We believe the reason that such regions are more commonly reported on Neptune than Uranus is that the overlying 1–2–bar Aerosol‐2 layer on Neptune has lower opacity, making the deep Aerosol‐1 layer more visible.

Interpreting the results of our radiative transfer analysis with respect to the structure of dark vortices is complicated, especially if the dark spots are to be explained by a darkening of particles in the deep haze layer. Particle colour (resulting from changes to n _imag ) cannot be unambiguously attributed to composition and/or heterogeneous microphysical structure, and thus linking these quantities to vortex dynamics and structure would be speculative given the limited information currently available. Even for Jupiter, where more detailed information is available, links between vortex structure and aerosol properties remain tentative (e.g., Wong et al., 2011). It is possible that dark spots are generated by secondary‐circulation uplift within vortices dredging up and concentrating a “chromophore” material from warmer depths below, which is dark at visible wavelengths. However, it is not clear what this material might be, and photolysis, which might be thought necessary to generate dark particles, is unlikely to be significant at pressures of ∼5–7 bar since the UV flux will be very low. Rather than interpreting such features as being caused by the addition of a new “chromophore”, an alternative explanation might be that such regions are anomalously warm at the H₂S condensation level and so H₂S ice sublimates into the vapour state, revealing the darker CCN photochemical haze core, similar to the mantling process proposed in (West et al., 1986) to explain belt‐zone color differences in Jupiter's atmosphere. Such a scenario may be consistent with cooling above the anticyclone mid‐plane and warming below, as seen in numerical simulations (Hadland et al., 2020; Lemasquerier et al., 2020; Stratman et al., 2001) and theoretical models (Marcus et al., 2013). However, more observations are needed to fully constrain why the particles in the Aerosol‐1 layer appear to be darker in the centre of dark spots and at dark latitudes.

Our results place interesting constraints on the vertical structure of dark spots. Geostrophic anticyclones in stratified fluids—whether salt lenses in the Earth's oceans or the Great Red Spot in Jupiter's atmosphere—are characterized by high density anomalies in their upper regions, and low density anomalies in their lower regions (e.g., dark contours in Figure 2 of Barceló‐Llull et al. (2017), Figure 3 of Marcus et al. (2013), Figure S8 of Lemasquerier et al. (2020)). Thus, how giant planet atmospheric vortices affect the observable aerosol structure depends on the location of aerosol layers with respect to the vortex mid‐plane. Anticyclones are also more efficiently mixed, compared to the stratified background atmosphere (e.g., Hassanzadeh et al., 2012).

Figure 22 compares differences in model Jupiter and Neptune vortex structures that might be consistent with our deep aerosol results. In Jupiter's case, observable aerosol layers (the NH₃ ice cloud layer plus the tropospheric haze above it) intersect the upper part of the vortex, where the high density anomaly is associated with cooler temperatures (e.g., Cheng et al., 2008). The cooler temperatures enhance haze (e.g., N₂H₄) condensation (West et al., 1986), and efficient mixing delivers cloud and haze precursors to high altitudes, producing the increased aerosols that are seen within the upper portion of Jovian anticyclones (Banfield et al., 1998).

Depleted, or darkened, aerosols in Neptune's dark spots, as possibly suggested by our analysis, may indicate that the vortex midplane is located at higher altitudes than the deep Aerosol‐1 layer. Higher temperatures in the low‐density anomaly would inhibit H₂S condensation. Depending on the vertical distribution of the dark haze CCNs populating the deep layer, mixing within the anticyclone may also play a role by altering the concentration of CCNs at the H₂S condensation layer. This scenario may also constrain the top of the vortex in this model. Efficient mixing within the anticyclone would disrupt the middle Aerosol‐2 layer near 1–2 bar, which is not observed. We have demonstrated that we can model the dark spot appearance at all HST/WFC3 wavelengths (Figure 10) without any changes to the middle Aerosol‐2 layer, which would mean that the high‐density anomaly of the dark vortex may reside entirely below this level. Models where both the Aerosol‐1 and Aerosol‐2 layers were depleted rendered vortices visible at longer HST/WFC3 wavelengths, in contradiction with the observations. If vortices are indeed bounded by approximately the 2‐bar and 10‐bar surfaces, their thickness is less than two atmospheric pressure scale heights. This thickness is much smaller than that for Jupiter, where remote sensing of the visible clouds constrains vortex tops to about the 0.5‐bar level (e.g., Banfield et al., 1998; Cheng et al., 2008; West et al., 2004) and vortex bases to levels near the 10‐bar level, or perhaps as deep as 800 bar in the special case of the Great Red Spot (Bolton et al., 2021; Parisi et al., 2021). Thermal infrared observations, however, (e.g., Fletcher et al., 2010) find that the thermal effects of vortices reach even higher to the tropopause. Jovian vortices thus span ranges of 3–7 scale heights. The potential difference in vortex thicknesses between Jupiter and Neptune may perhaps be related to the deep stratification of these outer planet atmospheres, which is still not known to high precision. Direct measurements at relevant pressure levels are only available for a single time and place: the Galileo Probe entry site at Jupiter (Magalhães et al., 2002; Seiff et al., 1998). Given the dependence of vortex vertical/horizontal aspect ratio on the difference between internal and environmental stratification (Hassanzadeh et al., 2012), and the roughly comparable horizontal dimensions of vortices on Jupiter and Neptune (Li et al., 2004; Wong et al., 2018), the thicker vertical dimensions of Jupiter's vortices may indicate that the environmental stratification is weaker on Jupiter compared to Neptune. However, this tightly constrained thickness scenario runs counter to the observation that such dark spot disturbances are often accompanied by what appear to be orographic clouds near the tropopause, for example, the “scooter” cloud that accompanied Voyager‐2's GDS. Hence, it would seem that the effect of the vortices must actually extend much higher. A possible scenario that might be consistent with all the evidence may be that the mid‐plane of the vortex coincides with the statically stable 1–2‐bar Aerosol‐2 layer, that is, the methane condensation level, which could leave this layer unaffected, but would affect both upper and lower clouds.

4.3. Dependence of Solutions on Methane Absorption Spectra

The main retrievals presented here have been conducted using k‐tables generated from the methane band model coefficients of Karkoschka and Tomasko (2010). We used these coefficients as they cover the entire spectral range of our observations leading, we hope, to self‐consistent results. The line data sets that are available, for example, WKLMC@80K (Campargue et al., 2013) have been shown to provide better fits to the H‐band observations (e.g., Irwin, de Bergh et al., 2012), but are limited in their wavelength coverage and so cannot cover the entire 0.3–2.5 μm region. Later databases, such as HITRAN16 (Gordon et al., 2017), and TheoRETS (Rey et al., 2018) extend the lower wavelength limit to 1.0 and 0.75 μm, respectively, but these still fail to cover the spectral region where dark spots are most visible. We repeated some of our retrievals for 0.76–2.5 μm region alone, using k‐tables generated from the methane band parameters (Karkoschka & Tomasko, 2010) and also from the TheoRETS database and found better fits to the observations with TheoRETS, and broadly similar retrieved aerosol structures and scattering parameters for the two different sources of methane absorption data. However, significant differences were found between the aerosol structures retrieved from the 0.76–2.5 μm using using the methane band/k‐data and those retrieved from the full 0.3–2.5 μm range, with significantly more opacity of the deep Aerosol‐1 layer estimated when using the restricted wavelength range. This may go some way to explaining the difference between the vertical profiles of aerosol retrieved here and those from previous studies, concentrating on the 820 nm or H‐band regions alone, that did not include a separate deep Aerosol‐1 layer and generally put the main Aerosol‐2 layer deeper at 2–3 bar, and for the 820 nm retrievals deduced a lower methane abundance for Neptune (Irwin et al., 2021). It is only the methane band data of Karkoschka and Tomasko (2010) that allows us to analyse the whole range simultaneously and deduce the presence of the deeper layer Aerosol‐1 layer, whose spatial variations are, we suggest, responsible for dark spot features.

We did attempt trying to combine the band and line data together, but there are significant differences in the wavelength regions of overlap and we found, generally, that the fits were worse when using a combination of methane data from different sources. In order to extend this analysis we really need for the line data sources to extend to 0.3 μm so that we are able to analyse the whole spectral range with a self‐consistent set of absorption data. We look forward to the time when this might be possible.

4.4. Spatial Deconvolution

While the HST/STIS observations presented here have been spatially deconvolved with a Lucy‐Richardson deconvolution scheme, this has not yet been possible with the Gemini/NIFS data we have analysed. For Uranus we do not consider this to be too significant a problem since the disc of Uranus seen in 2009 was moderately featureless. However, for the Neptune observations, there were significant levels of high methane ice clouds. Although we masked these regions when extracting the Minnaert coefficients of the background regions, given the likely shape of the Point‐Spread‐Function (PSF) it is possible that the unmasked regions were still contaminated by light from the bright, cloudy regions, which may help to explain why the NIFS Neptune data show very little limb‐brightening at methane‐absorbing wavelengths, in contrast to the expectations from the fits to the combined data and also to the behaviour seen in the Uranus NIFS observations. We hope, when possible, in future work to attempt to spatially deconvolve these Gemini/NIFS observations to explore this.