Neptune's Spatial Brightness Temperature Variations from the VLA and ALMA

We present spatially resolved ($0.1'' - 1.0''$) radio maps of Neptune taken from the Very Large Array and Atacama Large Submillimeter/Millimeter Array between $2015-2017$. Combined, these observations probe from just below the main methane cloud deck at $\sim 1$ bar down to the NH$_4$SH cloud at $\sim50$ bar. Prominent latitudinal variations in the brightness temperature are seen across the disk. Depending on wavelength, the south polar region is $5-40$ K brighter than the mid-latitudes and northern equatorial region. We use radiative transfer modeling coupled to Markov Chain Monte Carlo methods to retrieve H$_2$S, NH$_3$, and CH$_4$ abundance profiles across the disk, though only strong constraints can be made for H$_2$S. Below all cloud formation, the data are well fit by $53.8^{+18.9}_{-13.4}\times$ and $3.9^{+2.1}_{-3.1}\times$ protosolar enrichment in the H$_2$S and NH$_3$ abundances, respectively, assuming a dry adiabat. Models in which the radio-cold mid-latitudes and northern equatorial region are supersaturated in H$_2$S are statistically favored over models following strict thermochemical equilibrium. H$_2$S is more abundant at the equatorial region than at the poles, indicative of strong, persistent global circulation. Our results imply that Neptune's sulfur-to-nitrogen ratio exceeds unity as H$_2$S is more abundant than NH$_3$ in every retrieval. The absence of NH$_3$ above 50 bar can be explained either by partial dissolution of NH$_3$ in an ionic ocean at GPa pressures or by a planet formation scenario in which hydrated clathrates preferentially delivered sulfur rather than nitrogen onto planetesimals, or a combination of these hypotheses.


INTRODUCTION
Neptune is the prototypical 'ice giant': a giant planet composed mainly of elements heavier than hydrogen and helium by mass, such as oxygen, nitrogen, carbon, and sulfur. Within the cold environment of Neptune, the products expected to form from these elements in the observable atmosphere include: H 2 O, NH 3 , CH 4 , and H 2 S. These molecules provide clues to the bulk composition of the planet and so constraining their abundances is crucial for understanding Neptune's formation and thermal history.
Bright CH 4 clouds and aerosol hazes pervade throughout Neptune's upper atmosphere. Optical and near-infrared wavelengths are sensitive to these components, limiting views in the atmosphere to pressures less than ∼ 1 bar. H 2 S and NH 3 condense at higher temperatures and pressures than CH 4 , meaning optical and near-infrared observations are blind to their deep abundances. Radio wavelengths probe beyond these shallow features, resolving the structure of Neptune's atmosphere down to ∼ 50 bar. Thus, NH 3 and H 2 S profiles can be constructed by inverting radio spectra. While early radio observations could only obtain disk-averaged measurements of Neptune, the observed high brightness temperatures longward of 10 cm required NH 3 , a prominent microwave absorber, to be significantly depleted (de Pater & Massie 1985). This is possible if an NH 4 SH cloud forms at ∼ 50 bar and if the H 2 S abundance exceeds that of NH 3 , resulting in the complete removal of NH 3 during the cloud's formation. In such an atmosphere, an H 2 S abundance between 30 − 60× solar and an NH 3 abundance of ∼ 1× solar are required to fit the disk-averaged data (de Pater et al. 1991;DeBoer & Steffes 1996). Indeed, the detection of H 2 S spectral features near 1.58 µm in the tropospheres of the ice giants implies that their deep bulk S/N ratio is greater than one (Irwin et al. 2018(Irwin et al. , 2019a.
de  presented centimeter maps of Neptune from 2003, finding that the disk-averaged spectrum agreed with the abundances obtained from earlier radio observations. In addition, they found that the bright south polar cap must be significantly depleted in H 2 S down to ∼ 40 bar in order to match the observed brightness temperature at wavelengths of 0.7 − 6.0 cm. However, this study did not investigate brightness variations at other latitudes, as the sensitivity and resolution were not good enough to detect significant variations apart from those in the south polar cap.
In 2011, an upgrade of the VLA was completed. This expansion improved the continuum sensitivity by 5-to-20-fold and increased the wavelength coverage and bandwidth. This prompted a program to reobserve Neptune at centimeter wavelengths. The resulting maps, presented in this paper, show clear brightness temperature variations across the disk akin to that seen in millimeter maps produced from the Atacama Large Submillimeter/Millimeter Array (ALMA) (Tollefson et al. 2019). The mm and cm probe between ∼ 1 − 50 bar on Neptune, meaning the most complete picture of Neptune's upper atmosphere to date can be reconstructed by synthesizing these data. As NASA and ESA debate the merits of a next-decade ice giant mission, a firm handle on the uncertainties in the composition and dynamics of Neptune's upper atmosphere are critical . Just how well are the quantities of N, S, C, P, and O constrained from ground-based observations, in-situ, and orbital measurements, and how does the amount and distribution of condensibles affect the observed atmospheric dynamics Hueso & Sánchez-Lavega 2019;Fletcher et al. 2020)? Do our uncertainties on these elements necessitate an instrument like Juno's Microwave Radiometer (MWR) (Janssen et al. 2017) on a future spacecraft to Neptune (Rymer et al. 2020)? What bands and configurations would be most useful for Neptune atmospheric science with the next generation VLA (ngVLA) and could it replace an MWR equivalent?
This paper is organized as follows. First, we present longitudinal-smeared maps of Neptune taken with the expanded VLA in 2015 between 0.9 − 9.7 cm (Section 2). Next, we outline the structure of Neptune's upper atmosphere and the free parameters used in our modeling (Section 3). We then combine these new VLA maps with 2003 VLA and ALMA observations of Neptune and use a Markov chain Monte Carlo (MCMC) implementation of the radiative transfer code Radio BErkeley Atmospheric Radiative transfer (Radio-BEAR) to obtain retrievals for the abundance profiles of Neptune's condensibles (Section 4). Finally, we compare our findings to prior results (Section 5) and end with a summary of key takeaways (Section 6).

Data
This work makes use of two primary data sets: 1) ALMA millimeter observations of Neptune taken from 2016−2017, described in Section 2 of Tollefson et al. (2019); 2) centimeter observations of Neptune taken with the upgraded VLA taken in 2015, described below.
We observed Neptune with the expanded VLA, an interferometer located near Socorro, New Mexico, on 1 and 2 September, 2015. Maps were obtained at wavelengths of 0.9 cm (Band Ka), 2.0 cm (Band Ku), 3.0 cm (Band X), 5.1 cm (Band C), and 9.7 cm (Band S). The VLA consists of 27 antennas grouped into three arms of nine antennae to form a 'Y'-shape. Every four months, the configuration is changed by moving the antennae along tracks. The 'A' configuration is the VLA's most extended; the length of each arm is ∼ 21 km, forming a maximum baseline of ∼ 36 km. The maximum baseline is inversely related to the angular resolution, i.e. beam size, meaning variations across Neptune's disk are most distinct in the A configuration. The observation setup strongly impacts the shape of the beam. The symmetry of the VLA three-armed track, Neptune's low declination at the time of the observation, and extended time on source causing the uv-plane sampling to fill out due to Earth's rotation all contribute to the beams' shapes. 1 Neptune was observed on two days, each for 7 hours divided into many 5 minute scans which rotate through all wavelength bands. Thus, the degree of longitudinal smearing is high in the resulting maps as we observe a Resolution at sub-observer location, assuming equatorial and polar radii of 24,766 km and 24,323 km, respectively (Lindal 1992).
Neptune throughout nearly an entire 16.11 hr rotation period (Warwick et al. 1989). Table 1 lists a summary of our observations. We supplement these observations with VLA maps in 2003 from  in order to model Neptune's disk-average temperature. These observations were taken in 5 bands for a total of 8 hours in each wavelength, including calibrators. Three of these bands (2 cm, Ku; 3.6 cm, X; 6 cm, C) were also taken in the 'A' configuration while the other 2 (0.7 cm, Q; 1.3 cm, K) were taken in the 'BnA' configuration, which is a hybrid of the B and A configurations, and was ideal for imaging Neptune at low declination.

Calibration and Imaging
The VLA observations were loaded from the NRAO data archive and converted for use in the MIRIAD software package (Sault et al. 1995). Flagging, calibration, and imaging were performed within MIRIAD. The calibrators used were 3C48 (flux density) and J2246-1206 (phase). In addition, self-calibration was performed to correct for short-term variability in the phases caused by fast atmospheric fluctuations. The self-calibration model used was a limb-darkened disk that best matched the observations. The limb-darkened profile is represented by the peak brightness temperature T b multiplied by cos k θ, where θ is the emission angle and k is a limb-darkening constant. Values for T b and k were found at each wavelength such that the difference between the limb-darkened model disk and observations were minimized.
Our final longitudinally-smeared maps after subtracting the limb-darkened disk of Neptune are shown in Figure  1. While the images are 'longitudinally-smeared', pixel-to-pixel variations are still present on the disk due to the presence of artifacts on the order of the beam size. Figure

Error estimation
The random noise in our VLA maps is calculated by averaging over four regions of the sky with boxes equal to the diameter of Neptune and taking the root-mean-square (RMS). RMS values range from 0.4 − 2.0K and are similar to the RMS at a given latitude within the disk of a residual map. Table 2 lists our estimated errors in each band. The RMS does not include systematic effects, such as errors in the bandpass or flux calibration. Uncertainties in the flux density are estimated at 5% or less in each band so we use this as a conservative estimate for the absolute error in our disk-averaged temperature data.

MODELING
We generate models of Neptune's brightness temperature using the radiative transfer (RT) code Radio-BEAR 2 . Radio-BEAR generates synthetic spectra by solving the equation of radiative transfer through a model atmosphere. For a fuller description of Radio-BEAR, we refer the reader to de Pater et al. (2005, which outline how the model atmosphere is constructed from bottom-up and details the absorption coefficients and line profiles for the species in this work. A description of the temperature profiles, cloud structure, and compositions considered in the model atmosphere is given below. 3.1. Temperature Profile Figure 1. Longitude-smeared maps of Neptune taken with ALMA (top row) and the upgraded VLA (bottom row). The color scale has been chosen to enhance the brightness contrasts across the disk. All maps are residuals, where a uniform limb-darkened model has been subtracted to highlight the temperature contrasts between latitudes. The north pole is indicated by a white line in the ALMA maps; the VLA maps are rotated so that the north pole is pointing upward. Contour lines delineate the latitude transitions between bands, with a grid over a blank disk shown in the upper right for clarity. Neptune's disk is outlined with a white ellipse. The FWHM of the beam is indicated in white in the bottom left of each map.   Fig. 1 and are obtained by subtracting the observed brightness temperatures from best-fitting uniform limbdarkened disks, as described in in Section 2.2 Dashed black lines delineate regions where brightness temperature variations are evident and define the latitude bins used in our modeling. Lines may not extend to 90 • S due to limits in the resolution of the observation. The bump in the residual temperature at 9.7 cm is due to the poor resolution, as the large beam size impacts the limb-darkened fit.
Temperature profiles are calculated from deep in the atmosphere upward, assuming either a dry or wet adiabat such that the profile matches the Voyager 2 temperature of 71.5 K at 1 bar (Lindal 1992). At pressures shallower than 1 bar, the temperature profile follows that derived from mid-infrared inversions by Fletcher et al. (2014). Temperatures, pressures, and altitudes are related through hydrostatic equilibrium. We assume that the atmosphere is in local thermodynamic equilibrium. Figure 3 plots the wet and dry temperature profiles used in this paper. These are derived from a 'nominal' atmosphere with 30× protosolar abundances of H 2 S, CH 4 , and H 2 O, with 1× NH 3 , no PH 3 , 'intermediate' H 2 (see: Sec. 3.5), and 100% relative humidity for each condensible species 3 . Figure 3 also plots the nominal abundances of Neptune's condensibles.

Cloud Structure
The condensibles in Neptune's upper atmosphere are H 2 S, CH 4 , H 2 O, and NH 3 , all of which may condense to form clouds. Clouds expected to form on Neptune, from bottom-up, include: an aqueous ammonia solution (H 2 O-NH 3 -H 2 S), water-ice, ammonium hydrosulfide (NH 4 SH), H 2 S-or NH 3 -ice (whichever is left over after NH 4 SH formation), and CH 4ice (Weidenschilling & Lewis 1973;Atreya & Wong 2005). The cloud density may affect microwave measurements.  . The normalized weighting functions at nadir for the wavelengths used in this study (a), showing the pressure at which each observation is sensitive. Different gas profiles, temperatures, and geometries will alter the peaks of the weighting functions. The clouds expected to form on Neptune as a function of their density and condensation pressures are shown in (b). The cloud density and condensation pressure is determined from thermochemical equilibrium and the nominal gas profiles, pictured in (c), which assume 30× protosolar CH4 (light blue), H2S (green) and H2O (dark blue), and 1× protosolar NH3 (yellow) at the deepest levels (below all cloud formation). The cloud and gas structure is also linked to the assumed temperature profile; both the dry and wet adiabatic temperature profiles for the nominal gas abundances are plotted in (d). However, little is known about the cloud density on Neptune and clouds have been shown to not affect the microwave opacity on Jupiter . Therefore, we ignore the effect of cloud opacity in our modeling and focus instead on the effect of gas opacity.

Condensible Species
The gas opacity of the microwave spectrum for gas giant atmospheres is dominated by H 2 S, NH 3 , H 2 O, and the collision-induced absorption (CIA) of H 2 (we include: H 2 -H 2 , H 2 -He, and H 2 -CH 4 ) (de Pater & Mitchell 1993). To form the NH 4 SH cloud, H 2 S and NH 3 are reduced in equal molar quantities until the product of their partial pressures reaches the equilibrium constant of the reaction forming NH 4 SH. On Uranus and Neptune, the observation that NH 3 is absent while H 2 S is present above the NH 4 SH layer (de Pater et al. 1991;Irwin et al. 2018Irwin et al. , 2019aMolter et al. 2020) implies that this process takes up all of the NH 3 , leaving an excess of H 2 S gas above the layer. Therefore, while NH 3 and H 2 O are strong microwave absorbers, they are only abundant deeper than ∼ 40 bar and only impact Neptune's radio spectrum significantly at wavelengths longer than 10 cm. This is shown in Figure 3, which plots the normalized weighting functions at each observed wavelength assuming a nominal atmosphere. Our ALMA and VLA observations, with wavelengths shorter than 10 cm, are sensitive to pressures between 1 − 50 bar, peaking at altitudes at and above the aqueous NH 3 , H 2 O-ice, and NH 4 SH cloud formations. While NH 3 can not be probed directly with our wavelength coverage, the chemical connection between H 2 S and NH 3 means its abundance can be inferred. We therefore allow the NH 3 profile to vary unlike in Tollefson et al. (2019). In our model, the formation of the NH 4 SH cloud is governed by equilibrium chemistry described in Lewis (1969).
Our forward models allow the 'deep' abundances of gaseous H 2 S, NH 3 , and CH 4 to vary. We define 'deep' as pressures below the NH 4 SH cloud (forming at ∼ 50 bar) but above aqueous solution cloud formation (P > 100 bar). For the nominal abundances, the solution cloud in full thermochemical equilibrium removes about 5% of the interior H 2 S and 25% of the interior NH 3 . From our retrieved values of 'deep' H 2 S and NH 3 , we obtain abundances below all cloud formation. We set no a priori restriction on the 'deep' H 2 S or NH 3 abundance, meaning either may survive above NH 4 SH formation. H 2 S or NH 3 , whichever persists, will then form an ice cloud and we allow its relative humidity to vary (see Fig. 3, which shows H 2 S-ice forming, as the 'deep' H 2 S abundance is larger than that of NH 3 ).

Phosphine (PH3)
PH 3 is an important disequilibrium species on giant planets, tracing both chemistry and convective motion. In the deep atmosphere, PH 3 should oxidize to form P 4 O 6 and dissolve in water (Fegley & Prinn 1986). In the upper atmosphere, PH 3 is photolyzed and subsequent photochemical reactions may form to produce P 4 or complex polymers an compounds. Thus, PH 3 must be rapidly uplifted from the deep atmosphere to exist, making it useful as a passive tracer for both horizontal and vertical motions (Fegley & Prinn 1986;Fletcher et al. 2009).
PH 3 has not been detected on Uranus or Neptune, but if present is an important microwave absorber (DeBoer & Steffes 1996). The effect is most prominent at the PH 3 (1 − 0) rotation line at 266.9 GHz. The radio observations analyzed here only reach frequencies as high as 242 GHz (ALMA Band 6), meaning only the wings of the absorption line are detectable. However, if the abundance of PH 3 is large enough, the pressure-broadened rotation line will be detectable in our highest frequency data. Thus, upper limits on the uplifted PH 3 abundance can be placed.
Our forward model assumes a constant PH 3 abundance until the temperature and pressure become so low that PH 3 condenses (∼ 1 bar). We also only include PH 3 in our retrievals at latitudes where we expect upwelling and enriched condensibles, i.e. latitudes with comparatively cold brightness temperatures.

ortho/para H 2
The ortho-/para-H 2 fraction also influences Neptune's radio brightness temperature by modifying both the adiabatic lapse rate and the gas opacity (Trafton 1967;Wallace 1980;de Pater & Massie 1985;de Pater & Mitchell 1993). The ratio of ortho-to para-hydrogen in equilibrium depends on temperature; however, fast vertical mixing could bring the ratio of ortho and para states of hydrogen away from equilibrium and toward a "normal" ratio of three parts ortho to one part para. At latitudes where fast vertical mixing is unlikely, H 2 is presumed to exist in an "intermediate" state, proposed by Trafton (1967). In this case, the ortho and para states (which define the CIA properties) are set to the equilibrium value at the local temperature while the specific heat is set near that of "normal" hydrogen. "Intermediate" hydrogen is discussed further in Massie & Hunten (1982), who provide physical reasons for the choice of this "intermediate" state. "Normal" hydrogen has been shown to decrease the microwave brightness temperature relative to "intermediate" hydrogen (de Pater & Mitchell 1993;Luszcz-Cook et al. 2013a;Tollefson et al. 2019).
In our retrievals, we parameterize the state of H 2 as between 0.0 and 1.0, where 0.0 represents the fully "normal" state while 1.0 is fully "intermediate." Values in between represent a weighted average of the absorption coefficient between the two states. The specific heat is always set near to that of "normal" hydrogen, regardless of the parameter state.

Retrievals
In order to estimate uncertainties in our model parameters, we couple our Radio-BEAR models to Markov Chain Monte Carlo (MCMC) simulations via a python implementation of the Goodman & Weare (2010) ensemble sampler called emcee (Foreman et al. 2013). emcee has been used in near-infrared analyses of Uranus' and Neptune's hazes (de Kleer et al. 2015;Luszcz-Cook et al. 2016). We use similar log-likelihood Gaussian function and uniform/log-uniform priors as these authors; letting θ represent the free parameters in the model, the likelihood function ln p is: where T b,n is the observed brightness temperature, T b,m is the RT modeled brightness temperature, and σ n is the total uncertainty, all at a given wavelength ν. This approach is in contrast to the methods in Tollefson et al. (2019) who compared the ALMA data to forward models of Neptune's atmosphere. From these forward models, they obtained deep abundances of H 2 S and CH 4 that fit the latitudinally-varying brightness temperatures. A downside to this approach is that one can only rule out models which are improbable from χ 2 -statistics, meaning uncertainties are not retrieved.
We use 30 walkers and let MCMC run for between 4000−5000 steps. This run length ensured that the autocorrelation time was sufficiently shorter than the number of steps. All plots and tables show the distributions after burn-in, where  the range of retrieved values approximately appear like their final probability density. The burn-in phase ends after 50% of the steps have been completed. Observed and modeled brightness temperatures are obtained in each of the identified latitude bins listed in Section 2.2. One issue with RT models of particular regions on Neptune is that the finite size of the point-spread function (PSF) results in blurring of the disk. That is, the temperature within a particular latitude bin is a convolution of that latitude region, nearby latitudes, and sometimes the background sky. This effect is hard to model while conducting MCMC for two reasons: convolving the whole disk with the PSF is computationally expensive; and the model composition of the surrounding latitudes must be simultaneously known in order to obtain the exact temperature distribution across the disk. We circumvent this issue in two steps. First, we generate limb-darkened model disks of the best-fitting disk-averaged models (see Section 4.1) and then determine the brightness temperature in the region of interest with and without convolving the disk with the PSF. The ratio of the PSF convolved disk-averaged model temperature to that without convolution is called the PSF-scale. We multiply each MCMC retrieved model by the PSF-scale to obtain our final model brightness temperature incorporating the effect of convolution. This model temperature is fed into the emcee likelihood function. The effect of this scaling is smallest near the center of the disk and largest near the limb. Second, we add an error term to account for the uncertainty introduced by this approach, where σ T is the total error term used by emcee to calculate the likelihood function, σ CAL is the 5% calibration error and σ P SF is the error introduced by RT-modeling not accounting for the PSF. In our MCMC retrievals, we define an additional free parameter, σ, that is proportional to σ P SF : where T is the observed brightness temperature and S is the product of the beam's semi-major and semi-minor axes. Adding an extra free parameter in the uncertainty is general practice in MCMC, as it encompasses unknown sources of error (for instance, via the introduction of PSF-scale). We find that the maximum retrieved values of σ P SF are never larger than 5% of the product of T and S at disk center and 15% at the limb.

Atmospheric Models and Free Parameters
We consider four types of atmospheric models for these retrievals, depending on the latitude bin and expected dynamics, described below. Table 3 lists each of these models and the allowed free parameters. Figure 4 summarizes the effect of altering these parameters on the disk-averaged microwave spectrum of Neptune relative to the nominal model. Of note, the first two panels highlight the strong interaction between the deep H 2 S and NH 3 abundances as a function of wavelength. NH 3 can exceed H 2 S, resulting in NH 3 surviving NH 4 SH cloud formation, by either sufficiently depleting H 2 S (blue line, 0.5× S, in the H 2 S deep panel) or enriching NH 3 (yellow line, 10× S, in the NH 3 deep panel). The millimeter spectrum becomes brighter, as more NH 3 will reacts with more H 2 S at NH 4 SH formation, thereby reducing the total amount of H 2 S above it. On the other hand, the centimeter spectrum becomes colder, as NH 3 gas is a stronger microwave absorber than H 2 S gas. All other parameters have the strongest impact at λ ≤ 1 cm. Note, however, that our ability to retrieve the deep CH 4 abundance is limited compared to that of H 2 S. This is due to a combination of having less overall impact on the microwave opacity, lower than our calibration errors, and other sources like the H 2 S relative humidity, PH 3 , and the H 2 state having comparable or stronger effects in the mm. The considered models are: 1. Enriched Atmosphere: Used for latitude bands which have relatively cold brightness temperatures, where enriched and upwelling air is expected. The 'deep' 4 abundances of of H 2 S, NH 3 , and CH 4 are allowed to vary. The NH 4 SH, H 2 S-ice and CH 4 -ice clouds are allowed to form and the relative humidity of H 2 S is varied. The PH 3 abundance is allowed to vary; a uniform vertical profile is assumed up to the saturation pressure of PH 3 . The ortho-to para-H 2 fraction is allowed to vary between 0.0 − 1.0, representing the range of fully 'normal' to fully 'intermediate' states. Both dry and wet adiabats are considered.
2. Depleted Atmosphere: Used for latitude bands which have relatively warm brightness temperatures, where depleted and downwelling air is expected. The 'deep' abundance of condensibles is set to that of the enriched atmosphere at altitudes below the 'mixing pressure:' P mix . At shallower altitudes, the abundance of condensibles is varied by some fraction of the deep abundance. The final profiles for the condensibles look like step functions, with the transition at P mix . We allow the formation of H 2 S-and NH 3 -ice. Both dry and wet adiabats are considered.
3. Disk Average: The disk-averaged model for the atmosphere looks like that of the enriched, but we only let the 'deep' abundances of H 2 S, NH 3 , and CH 4 , and the relative humidity of H 2 S to vary. This model is done to compare to prior work and so both PH 3 and the H 2 state are not varied. We caution that disk-averaged retrievals poorly describe the physics within the planet. Both dry and wet adiabats are considered.    (Tollefson et al. 2019) to form our set of data used in MCMC retrievals. A summary of this data set is given in Table 2. Figure 5 plots the best-fitting H 2 S and NH 3 profiles and 500 random retrieved profiles from the posterior distribution for a wet adiabat. Figure 5 also shows the observed disk-averaged brightness temperatures versus wavelength, likewise plotting the best-fitting and 100 random retrieved model spectra assuming a wet adiabat. Table 4 lists the retrieved parameters. H 2 S is more abundant than NH 3 in both the dry and wet thermal profiles. There is also significant positive correlation between the deep NH 3 and H 2 S abundances. This is expected based on the NH 4 SH cloud chemistry; retrievals with larger NH 3 require additional H 2 S to remove the NH 3 during the NH 4 SH reaction to match the observations probing above the NH 4 SH cloud. The dry adiabat permits marginally larger abundances of all condensibles throughout the upper atmosphere due to the need to increase the opacity in order to offset the higher temperatures.
The latitude band between 36 • S−12 • S is dark in both the VLA maps presented here (Fig. 1) and ALMA residual maps. Low brightness temperatures imply increased opacity from Neptune's condensibles. Conversely, Neptune's south polar cap, between 90 • S−66 • S, is warm in both datasets, implying a lower abundance of condensibles. Due to the lower abundance of condensibles, the peaks of the normalized weighting functions are deeper in the atmosphere than for the nominal model (Fig. 3), meaning we are more sensitive to abundances of Neptune's condensibles below the NH 4 SH cloud at the south pole. As a result, we first obtain retrievals for the mixing pressure P mix and abundances of condensibles over the south polar cap. The retrieved probability density functions for the condensible abundances below P mix are then used as priors for the condensible abundance below the NH 4 SH cloud in the cold mid-latitude band between 36 • S−12 • S. The atmospheric models and free parameters used for these bands are described in Table 3 and Section 3.6.
The brightness temperatures fed into MCMC are obtained by averaging over all pixels that are within ±60 degrees of the sub-observer longitude at each latitude band. Likewise, the modeled spectra are obtained using the average emission angle within this zonal average. The uncertainty is dominated by the flux calibration error. σ PSF is also added to this uncertainty, though its effect is small compared to the calibration error.   NH3 below all clouds 4.7 +2.5 −3.7 × 10 −4 3.9 +2.1 −3.1 Table 6. As Table 5, but assuming a dry adiabat in each region.
Data for Neptune's latitudinal variations are not available in the VLA 2003 maps apart from those at the south polar cap. We do not use these data as the temperatures reported in  are the maximum temperatures within the south polar cap at each wavelength instead of the average.
Tables 5 and 6 list the 16th/50th/84th percentiles on the gas profile parameters for the simultaneous 36 • S−12 • S and 90 • S−66 • S models. If a parameter is not well constrained, the 97.5th-percentile (2σ) upper limit is given instead. Corner plots showing the covariance and probability distribution for the free parameters are given in Figure 6, assuming the temperature profile at the south polar cap follows a wet adiabat. As in the disk-average results, there is significant correlation between the deep H 2 S and NH 3 abundances. The presence of additional NH 3 at deep levels will lead to additional H 2 S at the NH 4 SH layer. However, H 2 S is the primary absorber above the NH 4 SH layer and so our data are very sensitive to its abundance. Therefore an increased NH 3 abundance requires a compensating increase in the H 2 S abundance. In addition, our results show that the H 2 S profile is constrained at the cold mid-latitudes but unconstrained at the south pole. The NH 3 abundance below the NH 4 SH cloud is somewhat constrained at each latitude band, though a long-tail forms at low abundances for retrievals with less H 2 S. In addition, P mix is well-constrained at the south polar cap and is inline with the NH 4 SH condensation pressure. The retrievals also indicate a constant amount of NH 3 above P mix at the south polar cap, though the median values are larger than the 12 ppb abundance required by  to best match their own observations. The deep CH 4 abundance is only weakly constrained as its effect on the radio opacity is small compared to the calibration uncertainty (see Fig. 4). Moreover, it competes with effects from PH 3 , the H 2 state, and the relative humidity of H 2 S. At 36 • S−12 • S, we obtain an upper limit of 85.4× protosolar or 4.1% mixing ratio at the 2σ level, assuming a wet adiabat globally. These upper limits are consistent with the disk-average amount of CH 4 reported by Baines et al. (1995) at 2.2%. However, the 4% deep CH 4 abundance favored by Karkoschka & Tomasko (2011) is similar to our 2σ upper limit. A global dry adiaba permits larger abundances of each condensible, with the CH 4 abundance completely unconstrained. Our CH 4 results may be reconciled with differences in the assumed thermal profiles, as the thermal profile assumed by Karkoschka & Tomasko (2011) is slightly warmer than a wet adiabat.
Our results are consistent with no PH 3 with upper limits of ∼ 0.4 − 0.8 ppm, depending on the assumed temperature model. We are completely insensitive to the ortho/para H 2 state.
The relative humidity of H 2 S is also unconstrained, apart from 2σ upper limits. As shown in Figure 4, the difference between relative humidities of a few percent or less on the retrieved spectra are minimal when compared to our uncertainties; depleting more and more H 2 S has an inconsequential effect on the opacity once enough is removed. Figures 7 and 8 plot the best fitting abundance profiles for H 2 S and NH 3 between 36 • S−12 • S and 90 • S−66 • S assuming wet and dry adiabat temperature profiles, respectively. Five hundred random profiles sampled from the burned-in posterior are also plotted, representing the range of allowable profiles. Brightness temperatures are also plotted, comparing one hundered random retrieved model spectra and data, showing an agreeable fit. A comparison  Fig. 7) unless H 2 S is supersaturated. For the dry thermal profile, the H 2 S-ice cloud forms at higher altitudes and larger H 2 S abundances are possible. However, H 2 S needs to be fully saturated or supersaturated to reach detectable upper atmosphere abundances.

Globally-Varying Profiles -Mitigating Systematic Uncertainty
One shortcoming of the analysis in the previous section is that the calibration error masks the evident brightness variations across the disk. In this section, we simultaneously model the atmospheric properties matching both the observed brightness temperature within the 36 • S−12 • S latitudinal band, and the observed temperature differences be- tween 36 • S−12 • S and other latitudes. This approach reveals the cause of spatial trends in the brightness temperature. The new likelihood function to maximize is: The first sum on the right-hand side is the same as equation (1) and represents the fit to the observed brightness temperature at 36 • S−12 • S with errors dominated by systematic calibration. The second sum represents the brightness temperature difference between 36 • S−12 • S and the other latitudinal band: Errors in the temperature variation between latitudes at the same wavelength are mainly due to random fluctuations. The uncertainty σ ∆ is: The first term on the right-hand side is the difference in brightness temperature between the reference and modeled latitude band times the 5% calibration error. σ RM S is the random error on the sky divided by the square root of the number of beams that fit within the reference and modeled latitude bands.
The deep abundances of H 2 S, NH 3 , and CH 4 are allowed to vary for each pair of latitudes. Moreover, H 2 S may either subsaturate or supersaturate. H 2 S may be supersaturated up to a pressure P ss . At altitudes shallower than P ss , the H 2 S profile follows the saturation curve with 100% relative humidity. At the south polar cap, the abundances of H 2 S and NH 3 at pressures shallower than P mix are also allowed to vary, as in the previous section. We assume no PH 3 and fully intermediate H 2 at all latitudes as we do not expect to be sensitive to these parameters given our prior findings (Table 6).
Tables 7 and 8 list the best-fitting values for 16/50/84th uncertainties for constrained parameters assuming a global wet and dry adiabat, respectively. Figures 9 and 10 plot randomly-sampled and the best-fitting retrieved trace gas profiles and residual spectra ∆T . We obtain good agreement between the observed and modeled temperature differences. Moreover, the retrieved parameters for the 36 • S−12 • S band are consistent between model runs to within uncertainties, indicating this approach is self-consistent.
These results show global variability in Neptune's NH 3 and H 2 S abundances. NH 3 , although overall very much depleted due to the formation of NH 4 SH, is enriched at the south polar cap above P mix relative to the rest of the planet. A globally uniform H 2 S abundance below NH 4 SH formation is consistent with our results. However, above NH 4 SH formation, the H 2 S abundance decreases poleward, mirroring brightness temperature increases poleward (Fig.  1). Differences in the H 2 S relative humidity cause the apparent brightness variations in the short-wavelength maps, which probe the H 2 S-ice cloud and are coldest / darkest at 36 • S−12 • S and 4 • N−20 • N. These regions are also where supersaturated H 2 S is needed in order to fit the observed temperature variations. This discussion is visually summarized in Figure 11, which plot the probability density of the retrieved H 2 S mixing ratio in each latitude band relative to that obtained for 36 • S−12 • S in the simultaneous fit. These plots demonstrate the magnitude and significance of these variations at the three pressure regimes discussed: below NH 4 SH formation (100 bar), above H 2 S-ice formation (5 bar), and in-between (20 bar). In the Appendix, we show that latitudinal variations in the thermal profile can only explain the observations if the temperature variations are localized to the H 2 S-ice formation pressures and are on the order of ∼ 8K or greater.
A measure of the improvement in the fit between two models is the Deviance Information Criterion (DIC, Spiegelhalter et al. (2002)). DIC is defined as: θ are the free parameters, D(θ) is the mean deviance of the retrieved parameters, which in this case is the average of the log-likelihood emcee probabilities, defined in Eq. (4). The second term on the right hand side is one-half the variance of these probabilities.
Information criteria like DIC are typically used when analyzing MCMC results in order to select models which provide a good fit to the data (the first term) and minimize the variance in the retrieved probabilities (the second term) resulting from model complexity. The difference between two DIC values can be used to determine the better model assuming the parameters roughly follow a Gaussian distribution. Differences greater than 10 favor the model with the lower DIC score. Table 9 lists the DIC score for the global wet and dry adiabat models at each latitude band. The dry adiabat is statistically favored everywhere, apart from 50 • S−36 • S where neither thermal profile is favored. We additionally show in Section 2 of the Appendix that when the lapse rate is allowed to vary with the trace gas abundances, the resulting retrieved temperature-pressure profiles follow a dry adiabat more closely than a wet adiabat.

DISCUSSION
Retrievals of Neptune's microwave spectrum from radiative transfer modeling require an atmosphere dominant in H 2 S over NH 3 . In every retrieval, all NH 3 is taken up by the formation of the NH 4 SH cloud around 50 bar, leaving H 2 S to persist and condense at higher altitudes. The global thermal profile has a moderate impact on the produced H 2 S and NH 3 profiles. Below all cloud formation, we obtain H 2 S abundances of 8.7 +2.6 −1.7 × 10 −4 (37.3 +15.7 −10.3 × protosolar) Positive values indicate that the modeled latitude band is warmer than 36 • S−12 • S at that wavelength.  Figure 10. As Fig. 9, but assuming a global dry adiabat. . Free parameters marked with a '− were not allowed to vary at that latitude band. 'not supersat.' and 'not subsat.' means that supersaturated and subsaturated H2S were not obtained in any final retrievals, respectively. The listed values are the 50th percentile of the posterior distributions with the error bars corresponding to the 16th and 84th percentiles, i.e. the 1σ uncertainties. Parameters for which the 16th and 84th percentiles vary by more than one order of magnitude have the 2.5th/97.5th percentile value listed instead, representing the 2σ lower/upper limit. If both supersatured and subsaturated H2S are possible in a latitude band, the 2σ limit is listed for both Pss and H2S H rel .
for a wet adiabat and 1.3 +0.4 −0.3 × 10 −3 (53.8 +18.9 −13.4 × protosolar) for a dry adiabat. We are less sensitive to the deep NH 3 abundance as all contribution functions probe at and above the NH 4 SH cloud. However, we do find a clear positive correlation between the H 2 S and NH 3 abundances: more deep NH 3 requires more deep H 2 S. We obtain loose bounds on the NH 3 abundance: 2.5 +2.7 −2.3 ×10 −4 (2.1 +2.3 −1.9 × protosolar) for a wet adiabat and 4.7 +2.5 −3.7 ×10 −4 (3.9 +2.1 −3.1 × protosolar) for a dry adiabat. Models fitting to the brightness temperature variations between latitude bands statistically favor a globally uniform thermal profile following a dry adiabat.
Using the H 2 S/NH 3 ratio as a proxy for the S/N ratio, we completely rule out observable S/N ratios less than unity. This is in agreement with earlier radio work on Neptune, although we obtain better constraints on these abundances (de Pater et al. 1991;DeBoer & Steffes 1996;Luszcz-Cook et al. 2013a;Tollefson et al. 2019). Sub-unity S/N is also ruled out on Uranus, suggesting some commonality in the formation history of the ice giants (Gulkis et al. 1978;de Pater et al. 1991de Pater et al. , 2018bMolter et al. 2020). However, their true bulk S/N ratio may be much smaller than observed or assumed, due to the potential loss of ammonia in an ionic/superionic water ocean at 20 GPa . Indeed, if CO is uplifted into Neptune troposphere, it implies that the water content and O/H ratio is several hundred times protosolar (Luszcz-Cook et al. 2013b). An additional way to explain the absence of N is the preferential delivery of volatiles onto planetesimals via hydrated clathrates (Hersant et al. 2004). In the cold environment of the outer protoplanetary disk, the N 2 to NH 3 fraction is roughly 10, while S is primarily in the form of H 2 S. Both H 2 S and NH 3 are readily trapped in hydrated clathrates while N 2 gas is not, meaning the majority of available N will not be swept into the ice giants. The clathrated hydrates hypothesis also implies that the O/H ratio within Uranus and Neptune is at least 100 times protosolar.
The differences in brightness temperature between latitude regions reveal prevalent global variations in the H 2 S and NH 3 profiles. Brightness temperature variations between latitude bands were modeled by comparing the retrieved spectra to profiles simultaneously fit to the 36 • S−12 • S region. This approach limits the calibration error, placing tight constraints on how the trace gas profiles vary latitudinally.
At the south pole, the data can be fit with a downwelling atmosphere, where air depleted in H 2 S and NH 3 species subsides down to a mixing pressure, P mix , below which the atmosphere becomes well mixed and the trace species equal their deep values. The best-fit retrieved H 2 S and NH 3 mixing ratios at altitudes above P mix are 4.4 × 10 −5 and 4.1 × 10 −7 , respectively, for a dry adiabat. This is in broad agreement with south polar models of 2003 VLA data presented in , whose best fitting depleted model required 3.5 × 10 −5 parts H 2 S and 1.2 × 10 −8 parts NH 3 above the NH 4 SH cloud.
Below NH 4 SH formation, our results are consistent with homogenous NH 3 and H 2 S. Clear latitudinal trends in the H 2 S profile, however, must be present at higher altitudes. Just above NH 4 SH cloud formation, the H 2 S abundance is highest between around the equator and mid-latitudes and diminishes poleward. Irwin et al. (2019a) tentatively detected 1 − 3 ppm H 2 S at Neptune's cloud tops between 2.5 − 3.5 bar. Their detection is more robust near the south pole than the equator and both their retrieved cloud top pressures and abundances increase toward the south pole. Our results assuming a wet adiabat globally are not consistent with their findings (for instance, see the gray rectangles in Fig. 7). The wet thermal profile is about 7K cooler at 3 bar than the Irwin et al. (2019a) profile. Only our dry adiabat models at the south pole allow H 2 S gas abundances close to the values found by Irwin et al. (2019a), as warmer temperatures push the H 2 S cloud base to higher altitudes.
Supersaturated H 2 S between 36 • S−12 • S and possibly between 4 • N−20 • N is required to obtain good fits to the observed brightness temperature variations if the thermal profile is fixed globally. Numerous near-infrared works favor a two-layer cloud/haze structure featuring a shallow tenuous upper haze and an opaque cloud deck between 2 − 4 bar (Irwin et al. 2014;Luszcz-Cook et al. 2016;Molter et al. 2019). This lower cloud is presumably H 2 S-ice, supporting our finding. The ratio of the partial pressure of the condensate in a supersaturated model to the partial pressure of condensate following its saturation curve, φ, is a few hundred percent in our H 2 S supersaturation model. Similar degrees of supersaturation are expected within ammonia plumes on Jupiter . However, on Earth φ does not exceed 10% for water (Young 1993). This is because the timescale of cloud formation is very quick in the presence of cloud condensation nuclei (CCN). Little is known about the amount of CCNs available on the gas giants. But, theoretical calculations by Moses et al. (1992) of homogeneous, heterogeneous, and ion-induced nucleation rates produce large values and ranges for φ, depending on the species. They find φ ≥ 3 − 1000 in order for methane and photolyzed hydrocarbon aerosols to form in Neptune's stratosphere. Neptune may lack substantive amounts of CCN conducive to aerosol formation. Future work should consider both the effect of cloud micro-physics and dynamics on the trace gas distributions within the ice giants.
The above discussion is synthesized in Figure 12, which shows one potential global H 2 S profile. For each latitudinal band, we assume the best-fitting H 2 S abundances and relative humidities from Table 8. Broadly, these meridional trends in gas abundances are related to the global circulation patterns of the planet.
Warm temperatures within Neptune's south polar cap are detected at mid-infrared and radio wavelengths. In the mid-infrared, (Hammel et al. 2007) observe a warm south pole explained by prevalent ethane and methane emission. (Orton et al. 2007) further argue that seasonal warming at the south polar cap may explain the excess of these molecules in the stratosphere; the warm temperatures would overcome the cold-trapping of methane below its condensation point at ∼ 1 bar and allow methane to escape upward, diffuse globally, and form ethane via photolysis. (Orton et al. 2007) also argue that rising air would be expected given the unexpected high abundances of these species. However, warm temperatures relative to the rest of the planet seen in the mid-infrared thermal emission (Fletcher et al. 2014) and a lack of widespread, persistent cloud coverage in the near-infrared (e.g., ) are more indicative of dry subsiding air. Radio observations, including these, observe high brightness temperatures at Neptune's south polar cap (Luszcz-Cook et al. 2013a;Tollefson et al. 2019), which we also interpret as dry, low-opacity, subsiding air. Cloud coverage surrounding the south polar cap and faint, distinct near-infrared clouds located near, but not at, the south pole might be indicative of vigorous convection (Luszcz-Cook et al. 2010), analogous to Saturn's polar activity (Dyudina et al. 2008;Fletcher et al. 2008). This vigorous convection may be a mechanism to explain all observations across the wavelength spectrum.
The mid-latitudes, defined loosely as between 50 • S−15 • S and northward of 15 • N, are where Neptune's brightest and strongest methane cloud activity is observed in the near-infrared (Martin et al. 2012;Fitzpatrick et al. 2014;Tollefson et al. 2018). In contrast, the equatorial region is nearly featureless. In addition, the mid-latitudes are colder than the equator and south pole in the mid-infrared (Conrath et al. 1998;Fletcher et al. 2014). Combining these observations, a global circulation pattern is inferred at altitudes shallower than ∼ 1 bar (where the near-and mid-infrared probe): cold, enriched methane air rises at the mid-latitudes and travels to the equator and poles, where the methane-depleted air subsides and warms via adiabatic compression. However, this picture is complicated by the relative excess of gaseous CH 4 at the equator, which is more consistent with rising air (Karkoschka & Tomasko 2011;Tollefson et al. 2018;Irwin et al. 2019b). We show in this paper that H 2 S is most abundant at the equator and southern mid-latitudes, in line with this meridional trend in CH 4 . At altitudes below methane condensation (∼ 1 bar), the aforementioned circulation scheme is, thus, at odds with the retrieved CH 4 and H 2 S abundances. Therefore, a more complicated picture of global circulation on Neptune is needed to explain each multi-wavelength observation. Fletcher et al. (2020) synthesize decades worth of analysis on the ice giants and argue that vertically-stacked circulation cells are necessary to bridge the observed patterns above and below 1 bar on both Uranus and Neptune.

CONCLUSIONS
We observed Neptune at radio wavelengths with the Very Large Array in five bands between 0.9 − 9.7 cm from 1 − 2 September, 2015. The longitude-smeared maps reveal brightness variations across Neptune's disk. These variations are alternating dark and bright latitudinal bands. Dark (bright) bands are consistent with high (low) opacity sources and cold (warm) brightness temperatures. We model Neptune's brightness temperature distribution using the radiative transfer code Radio-BEAR coupled to MCMC, varying the abundance profiles of Neptune's condensibles to obtain best-fits and retrievals to the observed microwave spectra. Models are fit to data from both the VLA presented in this work and 1 − 3 mm ALMA maps from Tollefson et al. (2019). Combined, these data probe from 1 bar down to >50 bar, where the NH 4 SH cloud forms. • The assumed thermal profile has a moderate impact on the retrieved H 2 S, NH 3 , and CH 4 profiles. A global dry adiabat is preferred over a wet adiabat as a warmer thermal profile is statistically favored in models fitting brightness temperature variations across the disk. All below results are given for the dry adiabat models.
In the mm, the CH 4 contribution to the modeled spectra competes with other opacity sources, namely PH 3 , the H 2 state, and the H 2 S relative humidity, dulling its signal in the retrievals.
• The downwelling south polar cap (90 • S−66 • S) is consistent with depleted H 2 S and NH 3 down to a mixing pressure P mix around the NH 4 SH cloud formation: ∼ 33 bar. Only when H 2 S is fully saturated or supersaturated can retrieved H 2 S abundances exceed 1 ppm at 3.5 bar at the south polar cap, as measured by Irwin et al. (2019a).
Only our dry adiabat model can produce values at pressures close to those detected.
• The observed brightness temperature variations between latitude bins are consistent with decreased H 2 S above the NH 4 SH cloud deck moving away from 36 • S−12 • S.
• The observed brightness temperature variations are more consistent with models supersaturating H 2 S between 36 • S−12 • S, and possibly 4 • N−20 • N, than with models following thermochemical equilibrium with some subsaturation. Models which do not allow H 2 S to supersaturate require latitudinal variations in the kinetic temperature on the order of 8K or greater to fit the data. These variations in the thermal profile are localized to H 2 S-ice formation pressures. Future work is needed to decouple the effects of temperature and gas opacity on Neptune's radio spectrum.
• We find no spectral evidence of PH 3 . Our retrievals are consistent with a PH 3 2σ upper limit of 0.4 ppm (≤ 0.1× protosolar).
• We cannot constrain the ortho/para H 2 fraction anywhere on Neptune.
Advances in radio astronomy are critical for further constraining the abundance of condensibles in the ice giants. A next generation VLA (ngVLA) would improve the resolution ten-fold at 20 cm if the longest 1000 km baselines are utilized (de Pater et al. 2018a;Selina et al. 2018). The resulting high-quality data at long wavelengths would permit a strong constraint on the NH 3 abundance beneath the NH 4 SH layer and perhaps even on H 2 O. In addition, a five-fold improvement in the sensitivity would mean vastly shorter integration times to achieve an equivalent RMS to this work, meaning zonal variations can be detected over a similar observing period and global maps can be obtained. Juno's big advantage over ground-based radio observations of Jupiter is that it flies beneath the synchrotron radiation belts that dominate longward of 6 cm. This is not an issue on the ice giants. Moreover, the resolution and sensitivity between the ngVLA and a MWR-like instrument on a potential future Neptune orbiter are quite similar .
The decision to include an MWR-equivalent on a direct mission to the ice giants will have to weigh whether it competes too heavily with the ground-based capabilities of an optimally running ngVLA. Whatever the choice, it is clear that the ngVLA or a direct mission to the planet are required to improve our maps of Neptune's trace species and elucidate theories regarding the environment in which Neptune formed and evolved.

ACKNOWLEDGMENTS
This paper makes use of the following VLA data from program VLA/15A-118. The National Radio Astronomy Observatory (NRAO) is a facility of NSF operated under cooperative agreement by Associated Universities, Inc.
VLA data used in this report are available from the NRAO Science Data Archive at https://archive.nrao.edu/archive/advquery.jsp.
This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00859.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.
This research was supported by the National Science Foundation, NSF Grant AST-1615004 to the University of California, Berkeley.
E. Molter was supported in part by NRAO Student Observing Support grant #SOSPA6-006.

APPENDIX .1. Globally Uniform Temperature Profiles
In this section, we address the assumption that the temperature profile is globally uniform at a given pressure level. We first show that this assumption leads to reasonable wind speeds according to the thermal wind equation, then demonstrate that not allowing the trace gas abundances to vary meridionally can result in well-modeled spectra only for nonphysical temperature profiles.
We test the validity of the assumption that the meridional temperature gradient is zero for P > 1 bar using an order-of-magnitude comparison. The thermal wind equation relates the vertical wind shear to both the meridional gradient in temperature and composition: Ω is Neptune's rotation rate (1 × 10 −4 rad s −1 ), θ is the latitude (rad), u is the zonal wind velocity (m s −1 ), r is altitude (m), r eq is Neptune's equatorial radius (2.5×10 7 m), g is Neptune's gravitational constant (10 m s −2 ), T is the temperature, q is the molar mass fraction of CH 4 (the most prominent trace gas on Neptune), and C = (1 − )/ where is the molar mass ratio of CH 4 compared to the ambient atmosphere. For the purpose of this exercise, we assume that the ambient atmosphere consists of H 2 , He, and CH 4 ; H 2 S and NH 3 have comparatively smaller abundances than CH 4 and are therefore less important to the compositional term. ≈ 7 and C ≈ −1. At the southern mid-latitudes, sin θ ≈ −0.5 and the above equation becomes:  Figure A1. Residual temperatures as a function of latitude, where the observed temperature at 36 • S-12 • S is subtracted from the temperature at the given latitude band. Colored lines are example model fits to the data assuming a nominal abundance profile and the same colored temperature profile given in Fig. A2.
From Karkoschka & Tomasko (2011), q is about 0.1 at Neptune's south pole (≤ 1 − 2% mixing ratio) and 0.2 toward the southern mid-latitudes and equator (≥ 2 − 4%). So dq/d(θ) ∼ +0.1 at the mid-latitudes for 1 ≤ P ≤ 3 bar. Thus, the right hand side is about 5 × 10 −8 and du/dr ≈ −0.5 m s −1 km −1 . The sign and magnitude of this estimate is similar to that of Tollefson et al. (2018) who observed vertical wind shear at Neptune's equator by tracking nearinfrared clouds. Therefore, the thermal wind equation alone cannot rule out the assumption of constant meridional temperatures.
Since the vertical wind shear is relatively unconstrained, especially below 1 bar, meridional kinetic temperature variations are not precluded from a thermal wind analysis. To test this, we estimate what temperature profile is needed to model the observed residual brightness temperatures. We fix the trace gas abundances globally to equal the nominal model (30×S H 2 S, CH 4 , H 2 O and 1×S NH 3 , with no PH 3 , 'intermediate' H 2 , and 100% relative humidity at all cloud formation). Fig A1 plots the observed and modeled temperature variations versus wavelength, meaning a dry adiabat model for 36 • S-12 • S is subtracted from the model temperature at each other latitude bin. The temperature profiles used to create these models are plotted in A2. To fit the residual temperatures, two adjustments to the dry adiabatic lapse rate are required. If the kinetic temperature were the sole cause of the brightness temperature variations, then substantial local variations on the order of 8K or greater are required. Interestingly, these peak around pressures where the H 2 S-ice cloud forms (∼ 10 bar), suggesting that the natural explanation for meridional brightness variability is due to changes in the H 2 S profile altering the microwave opacity. .

Varying Temperature and Trace Gas Abundances with MCMC
In the following, we consider a combination of kinetic temperature and trace gas abundance variations that may explain the observed brightness temperature latitudinal differences. To vary the kinetic temperature in the forward RT model, we allow the lapse rate to vary: Assuming hydrostatic equilibrium and the ideal gas law:  Figure A2. a) Temperature profiles as a function of latitude (colored lines) producing the residual spectra fits plotted in Fig. A1. b) Residual temperature profile as a function of latitude (colored lines), where the 36 • S-12 • S dry adiabat profile is subtracted from the temperature at the given latitude band.
where R is the ideal gas constant, g is the gravitational constant for Neptune, and for each layer P is its pressure, T its kinetic temperature, and M its molecular mass.
The molecular mass of the layer, M , is set by the mole fraction abundance of H 2 , He, and CH 4 where [He]/[H 2 ] = 18%, in line with infrared results from Burgdorf al. (2003), and [CH 4 ] = 1.44% (30× protosolar) at altitudes deeper than its condensation pressure. Equation (4) is integrated and solved for T (P ). For this work, T (P ) is altered by varying the lapse rate at altitudes shallower than a prescribed boundary pressure. At altitudes deeper than this boundary pressure, T (P ) is set to the dry adiabat profile used throughout this paper. We run three models fitting ∆T b between 66 • S-50 • S and 36 • S-12 • S. Model 1 sets the boundary pressure at the onset of H 2 S-ice formation and does not allow H 2 S to supersaturate. Model 2 sets the boundary pressure at the onset of NH 4 SH formation and does not allow H 2 S to supersaturate. Model 3 sets the boundary pressure at the onset of NH 4 SH formation and does allow H 2 S to supersaturate. Both the lapse rate and condensible profiles are varied independently between latitude bands. Figures A3 and A4 plot the retrieved gas and temperature profiles and their corresponding residual spectra, respectively, for each of the three models. Table A1 lists the retrieved model parameters, showing that similar NH 3 and H 2 S profiles are obtained regardless of the chosen boundary pressure, where the lapse rate begins to vary. Figure  A5 plots the χ 2 goodness-of-fit as a function of wavelength for the best fitting spectra in each of the three models considered. Only when supersaturated is allowed (model 3) can a good fit be obtained at 0.9 cm. Allowing latitudinal temperature variations without H 2 S supersaturation (models 1 and 2) can not explain all the observed data. In fact, the retrieved latitudinal differences in the kinetic temperature are consistent with zero, as shown in Figure A6. While the 0.9 cm datum is just one point of evidence in favor of supersaturated H 2 S, it can not be explained away by an error in calibration, as discussed in Section 4.3. Figure A3 also shows that the retrieved temperature-pressure profiles are consistent with a dry adiabat while few retrievals follow the lapse rate of a wet adiabat, consistent with the findings in Section 4.2. We therefore conclude that latitudinal variations in the opacity due to the condensible species are the cause of Neptune's observed brightness temperature distribution, and the current fitting favors supersaturation of H 2 S at the coldest latitudes. . χ 2 goodness-of-fit versus wavelength corresponding to the best fitting residual spectra ∆T b given in Fig. A4, showing that a good fit is achieved everywhere only when H2S is allowed to supersaturate (Model 3). Varying temperature without allowing supersaturation (models 1 and 2) does not produce a good fit at the 0.9 cm data point. The total χ 2 for each model is given in the plot legend.  Figure A6. Distribution of retrieved kinetic temperature differences between 66 • S-50 • S and 36 • S-12 • S at 5 bar for model 3.
Differences are consistent with zero and no larger than 1 K.