Emission from Volcanic SO Gas on Io at High Spectral Resolution

Jupiter's moon Io hosts a dynamic atmosphere that is continually stripped off and replenished through frost sublimation and volcanic outgassing. We observed an emission band at 1.707 $\mu$m thought to be produced by hot SO molecules directly ejected from a volcanic vent; the observations were made the NIRSPEC instrument on the Keck II telescope while Io was in eclipse by Jupiter on three nights in 2012-2016, and included two observations with 10x higher spectral resolution than all prior observations of this band. These high-resolution spectra reveal a contribution to the SO emission from gas reservoirs at both high and low rotational temperatures. The scenario preferred by de Pater et al. (2002) for the source of the SO gas - direct volcanic emission of SO in the excited state - is consistent with these two temperature components if the local gas density is high enough that rotational energy can be lost collisionally before the excited electronic state spontaneously decays. Under this scenario, the required bulk atmospheric gas density and surface pressure are $n\sim 10^{11}$ cm$^{-3}$ and 1-3 nbar, consistent with observations and modeling of Io's dayside atmosphere at altitudes below 10 km. The low-temperature gas component is warmer for observations in the first 20 minutes of eclipse (in Dec 2015) than after Io had been in shadow for 1.5 hours (in May 2016), suggesting cooling of the atmosphere during eclipse. Excess emission is consistently observed at 1.69 $\mu$m, which cannot be matched by two-temperature gas models but can be matched by models that over-populate high rotational states. Finally, a comparison of the total band strengths observed across eight dates from 1999-2016 reveals no significant dependence on thermal hot spot activity (including Loki Patera), on the time since Io has been in shadow, nor on the phase of Io's orbit at the time of observation.


Introduction
Io is the only object in the Solar System with an atmosphere that is dynamically created by active volcanism. Interactions with Jupiter's plasma environment strip material from the atmosphere at a rate of ∼1 ton/second (Schneider & Bagenal, 2007). Meanwhile, a combination of volcanic outgassing and sublimation of surface frost replenish atmospheric species, predominantly SO 2 (Lellouch et al. 2007). Numerous studies have investigated the role of both of these processes in sustaining Io's atmosphere, but have led to inconsistent conclusions regarding which process is dominant (e.g. Jessup and Spencer 2015;Jessup et al. 2004;Feaga et al. 2009;Spencer et al. 2005;Retherford et al. 2007;Roth et al. 2011;Tsang et al. 2012). Recent 19-µm observations have provided strong evidence that Io's SO 2 atmosphere decreases in column density by a factor of ∼5 during eclipse (Tsang et al. 2016). Although the discrepancies between various datasets remain unresolved, these authors suggest that differences could be explained if different atmospheric support mechanisms dominate at different Ionian longitudes.
Although Io's bulk atmosphere can be studied at wavelengths from millimeter through ultraviolet (e.g. Retherford et al. 2007;Tsang et al. 2012;Moullet et al. 2013), the direct volcanic component is challenging to isolate. In 1999, de Pater et al. (2002 detected emission near 1.7 µm while observing Io in Jupiter's shadow, and attributed the emission to the rovibronic band associated with the forbidden a 1 ∆ → X 3 Σ − transition of SO. The shape and the forbidden nature of the emission were indicative of high gas temperatures (≥1000 K), and the emission was attributed to hot, excited SO gas directly released from a volcanic vent. Follow-up observations on four nights between 2000 and 2003 found a tentative correlation between the strength of the emission and activity at Loki Patera ). De  performed the first spatially-resolved spectroscopy of the emission band and detected a hot plume over the volcano Ra Patera, which was undergoing a powerful eruption at the time.
Models of the 1.7-µm emission band have treated the SO source as a single gas in thermodynamic equilibrium, with gas temperatures of 500-1000 K (de Pater et al. 2002;Laver et al. 2007;de Pater et al. 2007). However, such models were unable to fit the shape of the band, and underpredicted emission both at 1.71-1.715 µm and around 1.69 µm. The authors speculated that the source gas may be out of local thermodynamic equilibrium (LTE), but the low signal-to-noise in the wings of the band prevented more detailed modeling.
We observed the 1.7-µm emission band on three occasions between 2012 and 2016, during eclipses of Io by Jupiter, including two observations with a factor of ∼10 higher spectral resolution than all previous data. Leveraging this significant improvement to spectral resolution, and the corresponding improvement in signal-to-noise, we investigate the source of the emitting gas via more complex models than were previously warranted. The observations, data reduction and calibration procedures are described in Section 2. In Sections 3 and 4 we describe our models and results; more details on the statistical methods are given in Appendix A. The implications are discussed in Section 5, and summarized in Section 6.

Observations and data reduction
We observed Io in eclipse by Jupiter with the W. M. Keck II telescope on Mauna Kea, Hawaii, on three dates between 2012: November 5, 2012December 25, 2015;and May 15, 2016 (dates in UT). On November 5 2012, we observed with the NIRSPEC instrument (McLean et al. 1998) in low-resolution mode (R∼1200), and on December 25, 2015 and May 15, 2016 we observed with the same instrument in high-resolution mode (R∼15000), where the structure of the line complex is significantly more resolved (although individual rotational transitions are still not resolved). During all three dates, simultaneous images of Io were obtained to identify the activity level of Io's volcanic hot spots; these data were obtained with the NIRC2 imager on Keck in 2012 and 2015 (Cantrall et al. 2018;de Kleer and de Pater 2016;de Pater et al. 2017a), and with the NIRI imager on Gemini N in 2016. In December 2015, we also obtained simultaneous data with the Keck/OSIRIS integral-field spectrograph, which spatially resolves the gas emission at low spectral resolution (de Pater et al. 2017b). In 2012 and 2015 Io was moving from sunlight into eclipse, while in 2016 Io was emerging from Jupiter occultation, and had been in shadow for an hour prior to the observation. The observing set-up varied somewhat between dates, and is described in detail for each date below. The details on the spectroscopic observations are summarized in Table 1.

November 5, 2012
On UT November 5, 2012 we observed Io in eclipse with the NIRSPEC instrument. NIRSPEC is equipped with a 1024×1024 InSb ALADDIN detector, and was used in low-resolution echelle spectroscopy mode for this observation, utilizing the NIRSPEC-6 spectral setting to cover the 1.62-2.02 µm region simultaneously. The peak spectral resolution in this mode is R∼2500. However, by adopting the 0.76" slit to cover as much of Io as possible, we reduce the effective spectral resolution to R∼1200. Both Io and the standard star were nodded along the slit for background-subtraction.
Thirteen nod pairs were obtained on Io with 30-second integration times per exposure, leading to 13 minutes of on-source time. All target exposures were used in the analysis.
The data reduction, including the spatial and spectral mapping and image rectification, was performed using the REDSPEC pipeline 1 , using telluric OH lines for spectral calibration. However, due to the uneven slit illumination caused by Jupiter, the subtraction of nodded pairs did not effectively remove background contamination, resulting in a poor pipeline extraction. Instead, the signal was extracted from a 30-pixel-wide band in the rectified images, and the 15 rows beside this band on both sides were averaged for background and subtracted from the signal. Figure 1 demonstrates that the Jupiter background is smoothly varying across the slit, and can therefore be characterized well and subtracted. The stellar spectra were not subject to this contamination, and were extracted successfully via the REDSPEC software.
Standard stars HD 31033 (A0) and HD 32333 (G0/1V) were observed for telluric correction and photometric calibration. Model stellar spectra were created from the Castelli and Kurucz (2004) 1 UCLA infrared lab; http://www2.keck.hawaii.edu/inst/nirspec/redspec Figure 1: Demonstration of Jupiter background removal from a pair-subtracted, rectified NIRSPEC image on Nov 5, 2012. (a) Pair-subtracted image; the horizontal dark and bright lines are Io's continuum emission, and the SO line emission is indicated by arrows in both the positive and negative spectra; (b) Interpolated Jupiter background, which is not removed well in the simple pair subtraction because it varies spatially along the slit; (c) Backgroundsubtracted image showing just the Io spectra with some residual Jupiter background. The background model is based on a polynomial fit to each column, and the scaling is the same in each frame. stellar atlas using the python pysynphot package (Lim et al. 2015). The measured stellar intensity occasionally dipped due to imperfect centering of the star within the slit and/or variable seeing conditions. Within each nod pair, if the difference in measured intensity was > 50% of the total measured intensity, the pair was not used in the calibration. Three out of seven stellar nod pairs were rejected based on this criterion.
The photometric calibration uncertainty arises from the variability in the median intensity of the stellar spectra used in the calibration (1-σ variation of 20%), the variability in median intensity of the Io spectra (1-σ variation of 5.8%), and uncertainty in the slit losses (10%; see Section 2.5), resulting in a total 1-σ systematic uncertainty on the absolute flux calibration of 23%. Within the stellar and target spectra, the pixel-to-pixel noise (residuals from a smooth curve) is at or below the 5% level.
The star and Io were both at <1.05 airmasses throughout the observations, and the spectra were not corrected for airmass differences because the effect at 1.7 µm is negligible. The photon flux at Earth, in units of photons/s/cm 2 /µm, was adjusted on all dates to correspond to an Io-Earth distance of 4.08 AU, as was done by Laver et al. (2007), to facilitate comparison with past work.

December 25, 2015
On UT December 25, 2015 we observed Io with both Keck telescopes simultaneously, using the NIRSPEC spectrograph in high-resolution mode (R∼25,000) on Keck II (McLean et al. 1998) and the OSIRIS integral-field spectrograph with adaptive optics on Keck I (Larkin et al. 2006), for which the grating was upgraded in 2013 (Mieda et al. 2014). The latter observations were presented by de Pater et al. (2017b) and described in detail in de Pater et al. (2018). The sky was clear during the observations but seeing was variable, in particular during the standard star observations. The observing date was selected due to the close proximity of Io to Ganymede, which was used for wavefront sensing for the adaptive optics OSIRIS observations. However, toward the end of the eclipse, Ganymede and Io were at their closest and glare from sunlit Ganymede overwhelmed Io's signal on the slit. The usable Io data on this night amounted to four nod pairs totaling twelve minutes of integration time (60-120 seconds per exposure).
The standard A stars HD 93702 and HD 104181 were observed for telluric correction and photometric calibration. Variable seeing during the standard star observations led to variability in the recorded count rate between exposures. Stellar nod pairs where the difference in counts recorded in the two nods was greater than 20% were not used in the calibration. This excluded four of the twelve standard star observations. The seeing was lower and more stable during the observations of Io, and variability between nods was minimal.
The photometric calibration uncertainty arises from the variability in the median intensity of the stellar spectra (1-σ variation of 20%), the variability in median intensity of the Io spectra (1-σ variation of ∼8%), and uncertainty in the slit losses (8%; see Section 2.5), resulting in a total 1-σ uncertainty on the absolute flux calibration of 23%. The noise along the dispersion direction is under 2% on this date, estimated from the stellar spectra. The star and Io were both at <1.25 airmasses throughout the observations, and the spectra were not corrected for airmass differences because the effect at 1.7 µm is at the 1% level, well below the uncertainty from other sources.

May 15, 2016
The observations on UT May 15, 2016 utilized the NIRSPEC spectrograph in high spectral resolution mode, with a set-up identical to that of December 2015, described above. However, during the May observations thick and variable cloud cover made much of the data unusable.
A high background, varying both spatially and temporally, necessitated nod pair subtraction to isolate the target spectra. We define a criterion for usable data in which the average signal-tonoise of the pair of spectra is greater than two, where signal is defined as the median level of Io's continuum above the median background, and noise is defined as the standard deviation of the sky background in the dispersion dimension, averaged over 10 spatial pixels (comparable to the extraction aperture), in the vicinity of 1.7 µm. Using this criterion, three of the seven total nod pairs were usable, amounting to twelve minutes of on-target integration time.
Due to the poor weather conditions, photometric calibration based on standard star data was not possible. Instead, we calibrate the data based on simultaneous Gemini N adaptive optics imaging with NIRI/ALTAIR (see Figure 5). The NIRI field of view included sunlit Callisto (the wavefront sensor guide star), and while Callisto was typically saturated, there were several low-transmission frames where Io was visible and Callisto was not saturated. We use these frames to calibrate the H-band intensity of Io in eclipse to that of Callisto in sunlight. We then use the spectral slope of both targets to extrapolate the 1.7-µm disk-integrated brightness of Io, assuming a 900-K Planck spectrum in the case of Io, consistent with previous spectra of Io obtained at this wavelength ). In the frames in which Callisto was not saturated, Io was also obscured to the extent that only two of the volcanoes were visible above the noise. We thus calibrated the Io images based on these two volcanoes, and used NIRI images from an earlier part of the eclipse when cloud cover was lighter to estimate that these volcanoes contributed ∼2/3 of the emission seen in the less cloudy images. Taking into account these uncertainties, we find a baseline 1.7-µm intensity of Io-in-eclipse of 3-6×10 −12 erg/s/cm 2 /µm at Earth, and scale the NIRSPEC spectra accordingly.
To facilitate comparison with past data, we scale these levels to what they would be if Io were at 4.08 AU, yielding a value of 4.6-9.2×10 −12 erg/s/cm 2 /µm. The range of values represent the 1-σ uncertainty level of ∼33% in absolute flux calibration. Along the spectral axis, the noise level was higher on this date than in the other observations, with a 1-σ noise level of about 10%.

Slit losses
Io's spatial extent exceeds the slit width, leading to significant differences in the slit losses between the calibration star and Io. The calibrated intensity of Io is therefore given by where F is the flux density of the object, N is the number of recorded counts/second, and f is the fraction of the intensity that fell on the slit for both Io and the star ( * ).
The fraction f * f Io is estimated by convolving models for the star and Io with a Gaussian point spread function (PSF) based on the seeing recorded on Mauna Kea during the observations (see   Table 4, which gives the latitude, longitude and identification of each source.    Table 1). The stellar model treats the star as a point source, while a range of models are used for Io to account for the unknown distribution of brightness on the disk, which may differ between the thermal and gas emission contributions. The Io models include a uniform-brightness model, as well as a range of models where 10-20 hot spots are placed randomly on the disk. The resultant distribution of slit loss correction factors are shown in Figure 6 as a function of seeing, for the 2012 dataset (0.76" slit and Io angular diameter of 1.208") and for the 2015 and 2016 datasets (0.72" slit and Io angular diameter of 0.979-1.00"). The shaded region on each plot indicates the uncertainty in the correction factor f * f Io based on the uncertain brightness distribution. These uncertainties are based on the assumption that the slit was aligned along the center of Io's disk; they may be underestimated if Io was not centered on the slit, as could happen if the thermal emission was asymmetric across the disk and biased guiding. Based on the seeing on Mauna Kea during the time of observations, given in Table 1

Doppler correction
For the low-resolution data from 2012, the line-of-sight velocity was around -10 km/s (Io moving toward Earth), and changed by ∼1.5 km/s during the event. The average Doppler correction is therefore around 5×10 −5 µm, much smaller than the wavelength sampling of 4×10 −4 µm. During the high-resolution observations, the line of sight velocity was around 20 km/s (Io moving toward Earth in 2016 and away from Earth in 2015). The Doppler shift of 1×10 −4 µm is several spectral elements (of size 2.35×10 −5 µm). The <1 km/s change in the line-of-sight velocity during each Figure 6: Photometry correction factor due to difference in slit losses between the star and Io. The thick line is for an Io uniform disk model, and the shaded region, centered on the dashed line, is for models that treat Io's brightness as a random distribution of hot spots. The models differ for the two dates due to the different angular size of Io, and different slit width.
observation is still well below the spectral resolution. Nevertheless, on each date, the spectra of Io were each individually tellurically-corrected and Doppler-corrected before being median-combined.

Thermal continuum
In the near-IR, Io's in-eclipse spectrum is dominated by thermal emission from active volcanism.
In the low-resolution dataset, this background is estimated by fitting the 1.6-2.0 µm spectrum excluding the SO emission region (see Figure 7). For the high-resolution datasets, the spectral coverage lacks sufficient baseline regions, and the thermal emission is left as a free parameter in the fits to the gas emission band. The thermal emission is modeled with a Planck function at a single temperature and emitting area (T th , A th ); this simple parameterization provides a good fit to the data, and the resultant temperatures fall within the range seen in previous observations with larger spectral coverage (∼500-1000 K; de Pater et al. 2002;Laver et al. 2007), as well as with the temperature estimated for December 25, 2015 from the OSIRIS dataset (de Pater et al. 2018).

Gas emission models
The SO emission band near 1.7 µm was first identified on Io by de Pater et al. (2002), and is produced by the 0-0 vibrational band of the forbidden electronic transition from the metastable a 1 ∆ excited state to the X 3 Σ − ground state. Setzer et al. (1999) identified 225 individual rovibronic lines over nine rotational branches within this band.
The SO emission line models are calculated for a given gas temperature according to the method described in de Pater et al. (2002) and Laver et al. (2007), and convolved with a Gaussian PSF  One source in thermal equilibrium In the low-resolution data from November 2012, the background thermal emission was well constrained due to the larger spectral coverage, and the parameters A th and T th were determined and surface emission subtracted prior to the gas model fits.
to match the spectral resolution of the instrument. For the high-resolution datasets, the spectral resolving power used for the model convolution is: R 0.43 slit =25,000, leading to a spectral resolution of R∼15,000 for the 0.72" slit. The lowresolution observations have a spectral resolution of R∼1250 after degrading from the nominal 0.38" slit (where R∼2500) to the 0.76" slit.
Three gas emission models are fit to the data: 1T: a single gas temperature, with the entire source in local thermodynamic equilibrium.
2T: two gas sources at different temperatures, each in local thermodynamic equilibrium. A summary of these models and their free parameters is given in Table 5.
One-temperature (1T) models: These models represent the simplest family of gas emission models, and treat the emission as arising from gas in local thermodynamic equilibrium at a single temperature. The free parameters are the gas temperature T and a scale factor c, in addition to the emitting area and temperature of the thermal continuum (A th and T th ). For all models, scale factors are used to match the models to the absolute photon contribution from each temperature component; the quoted results give integrated band strength instead of the value of the scale factors, as this is the more physically meaningful value. On all dates, the 1T models deviated from the data systematically near the band peak, leading to the more complex models described below.
Two-temperature (2T) models: The two-temperature models treat the SO gas as two reservoirs of molecules, each in local thermodynamic equilibrium at a different temperature. Preliminary fits to all datasets using a grid of temperatures demonstrated a consistent preference for a distinct low-temperature ( few hundred K) and high-temperature (>500 K) component. In order to avoid parameter degeneracy, we therefore restrict the potential range of one temperature component to <500 K, and the other to >500 K.
Disequilibrium (nLTE) models: The 2T models provide an excellent fit to the highresolution data, but are not able to match the shape of the emission band farther from the line center, in the spectral region covered by the low-resolution dataset from 2012 but not the highresolution datasets from 2015 and 2016. In order to fit the full band shape in the 2012 data, we postulate a gas out of local thermodynamic equilibrium (non-LTE). In particular, we allow an over-population of molecules at high rotational quantum numbers (J). A thorough exploration of non-LTE models would require more free parameters than can be uniquely fit by the low-resolution dataset, so instead we treat the non-LTE scenario with two simple parameterizations: when Io was also moving from sunlight into eclipse, as well as with previous measurements of Io's atmospheric temperature in sunlight (Lellouch et al. 2015). On top of this, we add a population of molecules at high rotational states only, populated according to a specific temperature that is independent of the background atmosphere temperatures. The inputs to the model are scaling factors for the low-T, high-T, and non-LTE components, the temperature according to which the high-J states are populated, and the J value above which over-population occurs. As in the previous case, we test a grid of temperatures and J-cutoff, fixing the value of each of these parameters at each point in the grid and fitting only for the relative populations. The temperatures are tested at 50 K increments between 50 and 2000 K, and all J-cutoff values from 20 through 59 are tested.
We note that these models are only two of many possible ways to parameterize a non-LTE state, and that the grid method does not allow for a robust statistical assessment of uncertainties on best-fit parameters. The best-fit parameter values should therefore be viewed with some caution, but the resultant models demonstrate that removing the LTE restriction allows for models that are capable of matching the broad-scale shape of the entire emission band.

Model fitting and selection criteria
In all cases where fits are performed, the best-fit parameter values are determined through a likelihood optimization routine, while the uncertainties are estimated using an affine-invariant ensemble Markov Chain Monte Carlo (MCMC) simulation. These methods are described in detail in Appendix A, where the methodology is shown for the specific example of 2T model fits to the high-resolution December 2015 data.
We calculate two metrics for the purposes of assessing the quality of the fit and comparing between models: the χ 2 value, and the autocorrelation of the residuals. The χ 2 value is normalized to N data -N par , where N data is the number of data points and N par is the number of free parameters in a given fit. We present the value of this metric alongside our final models to provide a quantitative evaluation of the goodness-of-fit. However, we note that for non-linear models the number of degrees of freedom does not correspond directly to N data -N par , and the interpretation of this metric for the purposes of model comparison is therefore not straightforward (Andrae et al. 2010).
The second metric we calculate for each model is the autocorrelation of the residuals. Under the assumption that the noise in the data is random and not correlated in wavelength, the residuals should also exhibit no correlation between different wavelengths in the case of a perfect model fit, while an incorrect model will lead to systematic deviations from the data and hence a positive autocorrelation of the residuals. We calculate the autocorrelation between datapoints separated in wavelengths using the Pearson R coefficient, which is the ratio of the covariance between the two datasets and the product of the standard deviations of each dataset: Here, the two datasets are the residuals from a given model fit, and the same residuals offset by n pixels; we calculate this coefficient for offsets from one to five pixels.

Results
The best-fit one-and two-temperature gas models for all three datasets are shown in Figures   8 and 9, and the corresponding best-fit thermal and gas temperatures are given in Table 6. The normalized χ 2 value and R coefficient for autocorrelation of the residuals are reported in Table 7. a For the Nov 2012 observations, T T h was determined from regions of the spectrum free from gas emission, and the thermal surface spectrum was subtracted from the data prior to fitting the gas model. On the other dates, the surface and gas temperatures were fit for simultaneously.
The latter is reported for wavelength intervals of n=1-5 pixels; R values not consistent with zero indicate that the residuals are correlated across n datapoints, an indication of systematic deviations between model and data rather than random noise. For all three dates, the 1T models provide poorer fits, as seen in the higher χ 2 and in the systematic deviations from the data near the band peak, which can be seen qualitatively in the figures and quantitatively through the values of the Pearson R coefficient, which indicates a greater correlation between wavelengths for the 1T models than for the other models tested.
The 2T  The low-resolution data from November 2012 have larger spectral coverage, and show an excess of emission around 1.69 µm also seen in previous datasets (de Pater et al. 2002;Laver et al. 2007) that cannot be captured by the equilibrium models. Two disequilibrium model families are tested, one in which the entire emission band arises from a gas population where the high and The parameters corresponding to the models shown here are given in Table 6, and goodness-of-fit metrics are tabulated in Table 7. On both dates, the 1T models show deviations from the data near the band peaks, which are better matched by the 2T models. The parameters corresponding to the models shown here are given in Table 6, and goodness-of-fit metrics are given in Table 7. Notes: χ 2 N = χ 2 /(N data − Npar), the value of χ 2 normalized by the number of datapoints minus the number of free parameters in the fit. The Pearson R statistic is consistent with zero if there is no significant autocorrelation in the residuals. The low and high confidence intervals are given in brackets, and italics indicate lags with no significant autocorrelation.  low rotational states are populated according to independent temperatures, and one that places an excess population at high rotational states on top of a two-temperature model with fixed high and low temperatures of 1000 and 170 K. Both models provide equally good fits to the data, and improve on the standard 2T models for the November 2012 data (see Table 7). In particular, disequilibrium models that over-populate molecules with rotational quantum numbers above J ∼ 35 do a better job of capturing the overall shape of the band, including the excess of emission around 1.69 µm. In addition, some of these models produce an excess of emission around 1.725-1.73 µm where excess emission at the few-σ level is also seen in the data. These models demonstrate the ability of disequilibrium gas states to resolve the discrepancy between the simpler 2T models and the data.
Future data with higher spectral resolution across this entire band would allow for a more in-depth modeling and fitting procedure.
For the disequilibrium models that populate high and low rotational states according to different temperatures (nLTE Model A), the χ 2 statistic for the best-fit model for each temperature pair is shown in Figure 11, normalized to the number of datapoints minus the number of free parameters in the fit. The figure demonstrates a preference for a lower temperature around 500 K, although the temperature of the high-T population is poorly constrained. Figure 12 shows the model with the lowest χ 2 value from this family, as well as the lowest χ 2 model from the fixed background disequilibrium model family.

Integrated band strength
Integrating over wavelength across the SO band yields a total photon rate on each date under the assumption of isotropic emission; these values are tabulated in Table 8  respectively. Note that the numbers from the datasets taken at different spectral resolutions are not directly comparable because the high-resolution datasets cover a smaller portion of the band. By using the November 2012 data, with full spectral coverage of the band, we estimate that the values derived from the high-resolution observations should be increased by ∼25% to yield the integrated emission across the entire band. Table 8 also gives the breakdown in emission between the two temperature components in the 2T models, and shows that the emission from both components varies considerably between dates. A direct comparison of the high-resolution spectra between 2015 and 2016 is shown in Figure 13 with thermal continuum removed (as determined from the 2T model fits).

Discussion
Model fits to the high-spectral-resolution datasets indicate that emission arises from SO gas reservoirs at both high and low temperatures, with high temperatures in the 1000-1500 K range and low temperatures around 75-200 K. For the a 1 ∆ excited state of SO, de-excitation due to spontaneous emission is rapid, and the consistent presence of the 1.7-µm band requires an ongoing production mechanism. De Pater et al. (2002) investigated several excitation mechanisms, and Figure 13: Comparison of high-resolution spectra between the two dates of observation, after subtraction of the thermal continuum according to the best-fit parameters for the 2T models (see Table 6). Despite the poor weather conditions during the 2016 observations, and the resultant noisy data, the band structure even outside of the main peaks is similar on both dates. a Integrated emission for the observed region of the spectrum. The high-resolution data (Dec 2015 and May 2016) do not cover the entire band, and underestimate the total emission by ∼25%. The thermal continuum for these two datasets, and the relative contributions from the two temperature components in all cases, are as derived from the 2T models. The total SO emission is as measured directly from the observed spectra.
concluded that the most likely scenario involved SO molecules being directly ejected from a volcanic vent in the excited state. This scenario predicts molecules in the excited a 1 ∆ state and with high rotational temperatures, matching the high-T component.
We now address potential origins for the population of molecules at low rotational temperatures yet in the excited electronic state. The temperatures found for the low-T component are consistent with temperatures previously measured for Io's bulk atmosphere. In Dec 2015 when Io had just entered eclipse, the 172-217 K temperature range is consistent with the 170±20 K measured for Io's dayside atmosphere at 4 µm (Lellouch et al. 2015), though somewhat higher than the <150 K temperatures inferred from mid-infrared data (Spencer et al. 2005;Tsang et al. 2012). In May 2016, Io had been in eclipse for more than an hour prior to the start of observation, and the gas temperature of 76-119 K is consistent with Io's in-eclipse surface temperature of ∼105 K derived from Io's 19-µm brightness (Tsang et al. 2016).
At temperatures of 100-200 K, the thermal energy of the gas (∼0.02 eV) is much smaller than the energy difference between the ground and excited SO gas (0.73 eV), and the cool SO reservoir will not populate the excited state in thermodynamic equilibrium. Molecules with rotational temperatures matching the bulk atmosphere could end up in the a 1 ∆ state in two ways: either the electronic state is excited without raising the rotational temperature, or else a fraction of the gas population cools rotationally before de-exciting to the ground electronic state. The former case requires a non-volcanic excitation mechanism, but de Pater et al. (2002) demonstrate that other processes (solar photoexcitation, electron impact excitation of SO or dissociative excitation of SO 2 , and ionospheric recombination of SO + 2 ) fail to explain the data. The latter scenario is plausible if the timescale to equilibrate the rotational temperature to the surrounding temperature is on the same order as the timescale of decay from the excited electronic state.
We test the consistency of this scenario as follows: we assume that all excited SO gas was initially emitted from a volcanic vent at high rotational temperature, and that the low-temperature emission is sourced from molecules that have equilibrated rotationally to the bulk atmosphere before de-exciting. We then make estimates of the local gas density and pressure from the relative emission from high-and low-T states, and determine if the resultant gas pressures are consistent with past observations.
In the scenario described above, the cool excited SO molecules are produced solely through cooling of the hot SO, and are lost solely through spontaneous de-excitation, so that where n H and n C are the number densities of hot and cold excited SO molecules, A is the Einstein coefficient for spontaneous emission from the a 1 ∆ state (∼2.2/sec; Klotz et al. 1984), and τ is the timescale for collisional cooling. Assuming that the hot and cold SO gas densities remain constant over timescales comparable to the duration of the observations (an assumption supported by the lack of measurable changes in individual spectra over a single observation), the production rate and loss rates of the cool excited SO must be equal: The photon rate is equal to the loss rate from spontaneous de-excitation, and the ratio of total photon rate in the 1.7-µm band to the photon rate from cool SO molecules only is therefore (1+n H /n C ), or (1+Aτ ) using Equation 6. Using the measured photon rates given in Table 8, we derive the timescale for collisional cooling τ . This timescale is the product of the time between collisions τ coll and the number of collisions to equilibrate (∼tens for exchange of translational and rotational energy). The former is defined as the ratio of the mean free path λ to the mean velocity of molecules in the gas v th , while the mean free path is determined by the collisional cross-section σ and the density of the gas n bulk . In this case, n bulk is the density of the gas in Io's bulk atmosphere, which is the reservoir responsible for the collisional cooling of the SO: For the 2015 and 2016 observations respectively, the ratio of the photon rate from all molecules in the a 1 ∆ state to just the low-temperature molecules is 2.4 +1.4 −0.9 and 4.1 +4.0 −2.0 , indicating a time to equilibration of around 0.6 and 1.4 seconds, or 0.06 and 0.14 seconds between collisions (assuming ∼10 collisions are required to reach rotational equilibrium with the surroundings). These numbers can be translated into number densities using the kinetic velocity of the gas: where k B is the Boltzmann constant and m is the molecular mass of SO 2 , the dominant atmospheric constituent. This yields number densities of 1.3 +2.5 −0.7 ×10 11 and 8.6 +17 −5.3 ×10 10 cm −3 , or surface pressures around 3.3 and 1.1 nbar, for the 2015 and 2016 observations respectively. Despite the approximate nature of these calculations, these estimates are consistent with past measurements of Io's dayside atmospheric pressure, which have typically been in the few nbar range (Lellouch et al. 1990;, as well as with the density of Io's near-surface dayside atmosphere predicted by numerical models based on observational constraints on column density (Walker et al. 2010).
However, these densities are inconsistent with the orders-of-magnitude decrease in both gas density and surface pressure expected for a purely sublimation-supported SO 2 atmosphere (Walker et al. 2010). Indeed, recent mid-infrared observations of Io's atmosphere have shown that while the atmospheric SO 2 density does drop at night, it only does so by a factor of 5±2 in density, similarly less than the several orders-of-magnitude decrease predicted by pure vapor pressure equilibrium.
During the daytime, Io's scale height of H s ∼10 km is significantly larger than the mean free path (tens of meters). At night, whether or not the atmosphere is collisional, and hence whether the equations above are applicable, depends on the factor by which the gas density decreases. If the density drop is on the order of a factor of five (Tsang et al. 2016), the 5× longer mean free path is still much smaller than the atmospheric scale height, even accounting for the ∼50% decrease in scale height due to the temperature decrease.
The factor of three difference in atmospheric pressure for the post-ingress and pre-egress atmosphere inferred from the 2015 and 2016 observations, respectively, is due largely to the decrease in temperature; the inferred gas density decreases by only 30%. Note that the lower emission from the low-T component in the May 2016 data (see Table 8) can likewise be explained by a lower gas temperature leading to longer time between collisions, and does not require a significant difference in gas density. Although the best-fit densities therefore imply an atmosphere with minimal change (n ingress /n egress =1.5), the ratio of densities between ingress and egress is fully consistent with the ratio of 5±2 measured by Tsang et al. (2016) within uncertainties. Moreover, while Tsang et al. (2016) isolated the density of the SO 2 gas, our observations are sensitive to the total bulk atmosphere, which includes other species. The higher n ingress /n egress ratios consistent with our observations correspond to the low end of the acceptable post-ingress temperature range (172 K) and high end of the pre-egress range (119 K), which are also the values most consistent with past work.
While the temperature of the low-T gas component in the 2015 and 2016 datasets is therefore consistent with SO gas in rotational equilibrium with the surrounding atmosphere for an atmosphere that cools from ∼170 to ∼100 K during eclipse, this interpretation is complicated by the fact that fits to individual spectra taken during the eclipse do not show a systematic trend in gas temperature (although signal-to-noise is low in individual spectra). Moreover, even the data taken shortly after Io enters eclipse span a long enough time period that the atmosphere should have cooled if the atmosphere follows the same cooling timescale as the surface. Specifically, the November 2012 data were taken over the first 32 minutes of eclipse, and the December 2015 data during the first 17 minutes of eclipse, while Tsang et al. (2016) find that the change in atmospheric density occurs within the first 20-30 minutes, and that the surface cools from around 130 K to below 110 K within the first 10 minutes of eclipse. If the difference between the low-T gas temperatures derived from the Dec 2015 and May 2016 observations does reflect a cooling of Io's bulk atmosphere, the majority of the cooling must therefore be occurring after Io has been in shadow for >20 minutes, implying a delayed cooling of the atmosphere relative to the surface.
The 2T models discussed here for the high-resolution datasets are still unable to match the shape of the emission band farther from the line center, as seen in the low-resolution data presented here and in past work (de Pater et al. 2002;Laver et al. 2007;de Pater et al. 2007). In particular, an excess of emission at 1.69 µm is observed that cannot be fit with equilibrium models. Models that postulate an over-population of high rotational states do a better job at resolving this mismatch.
At low densities, such an over-population could arise from either preferential radiative excitation of the high-J levels or a difference in Hönl-London factors for high-and low-J lines. The latter would lead to different critical densities for each of the transitions, such that it would be possible for some levels to thermalize while others remain out of equilibrium. We note that at 1.69 µm emission arises almost exclusively from the S R branch of the band, and it is also possible that the over-population is occurring in just this branch, or that the laboratory measurements of oscillator strength are incorrect. As the high-resolution data do not cover these portions of the spectrum and the low-resolution data are unable to discriminate between potential disequilibrium models, a more in-depth characterization of the gas temperature and state in the case of non-LTE is not possible from these datasets. Future observations covering the entire emission band from ∼1.67-1.74 µm at high spectral resolution are required to discriminate between more realistic gas models.

Investigating the source of time-variability in total SO emission
The three new datasets presented here, combined with past work (de Pater et al. 2002;Laver et al. 2007), bring the total number of observations of this emission band to eight, spanning the period from 1999 to 2016. Using all eight observations, we explore possible origins of the emission by determining whether the total amount of gas emission is correlated with (1) Activity at Loki Patera, as was suggested by early observations (de Pater et al. 2002;Laver et al. 2007 For each property, we calculate the Pearson R coefficient to test for correlation, and give both the coefficient and the corresponding p-value. A p-value below 0.05 indicates that the null hypothesis that the SO emission is not correlated with the given property is rejected at the 95% level. Positive coefficients indicate positive correlations, and negative coefficients indicate anti-correlations. 1. Activity at Loki Patera: Initial observations of the 1.7-µm SO band were suggestive of a correlation between band strength and the level of thermal activity at Loki Patera (de Pater et al. 2002;Laver et al. 2007). Observations that spatially resolved the SO emission did not find a significant contribution from the vicinity of Loki Patera, but the authors note that the observations were made during a quiescent period at this volcano . Figure 14a shows the integrated band intensity as a function of Loki Patera's brightness near a wavelength of 3.5 µm, for all eight observations. While seven of the eight observations follow a suggestive trend, the strong emission seen in November 2012 when Loki Patera was quiescent (and no other volcanoes showed bright thermal activity) deviates significantly from the trend, and there is no statistically-significant correlation (R=0.47; p=0.28). Note that observations at 3.5 and 3.8 µm are plotted together in Figure 14a, although the flux density at 3.5 µm may be lower than that at 3.8 µm by up to 50% for continuum thermal emission at the cool 300-500 K temperatures typically seen at Loki Patera.

Bright eruptions:
Thermally-bright eruptions can be accompanied by massive plumes of gas and dust (Geissler and Goldstein, 2007). However, observations of the 1.7-µm SO band made while a large eruption at Ra Patera was in progress in November 2002 found that only 10-15% of the total emission was localized to the area around this volcano , and the integrated band strength on this date was the weakest of all observations to date. On December 25, 2015 a large, transient eruption was in progress at Amaterasu Patera, and the integrated band strength on this date was also not anomalously large. On the dates with the highest measured SO emission (Sep 24, 1999 andNov 5, 2012), there was no especially bright activity, although Loki Patera was in the early stages of a bright event during the 1999 observations (de Pater et al. 2017a).
We therefore conclude that the thermal emission is not a good indicator of the amount of excited SO emission, and if the SO is of volcanic origin it is not dominated by events with large thermal signatures.
3. Io's mean anomaly: Since the amount of 1.7-µm SO emission is not dominated by bright volcanic eruptions, perhaps it is driven by volcanic gas emission from fissures that don't have active lava flow components. In such a case, the amount of emission could conceivably be modulated by tidal stresses opening and closing cracks such as in the case of Enceladus (Hedman et al. 2013;Hurford et al. 2007). We test this hypothesis by looking at integrated band strength as a function of Io's mean anomaly; Figure 14b illustrates that such a correlation is not present at a detectable level (R=0.47; p=0.24).

Io's solar distance:
The density of Io's bulk SO 2 atmosphere is correlated with solar distance, which suggests that there is a significant component that is controlled by sublimation (Tsang et al. 2012). If the 1.7-µm SO emission arises from SO molecules in Io's bulk atmosphere, which is primarily sourced from SO 2 , we might expect to see a similar correlation. Figure 14c shows the integrated SO band intensity as a function of Jupiter's solar distance; although there is a weak trend in the expected direction, a best-fit line to the points is consistent with zero slope to within 1σ and there is no statistically-significant correlation (R=-0.46; p=0.25).

Time since in sunlight:
Of the eight observations, all but two were taken shortly after Io had moved from sunlight into eclipse. For the two datasets taken with Io moving from Jupiter occultation into eclipse (i.e. having been in shadow for ≥1 hour prior to observation), the integrated line strengths were 1.6±0.5 and 1.45±0.38 ×10 27 photons/sec (datasets from 2003 and 2016 respectively), which are both average line strengths relative to the other datasets. It therefore appears that the time that Io has been in shadow prior to observation is not the dominant control on the intensity of emission and hence that, if the bulk SO atmosphere is sourced photolytically and is lost to surface interactions over the timescale of an eclipse, then the observed SO emission must be primarily of volcanic origin.
6. Jupiter's System III longitude: If SO molecules were excited to the a 1 ∆ state via electron impact processes (a possibility considered and rejected by de Pater et al. 2002), we might expect to see a correlation in emission with Jupiter's System III longitude. The SO band strength is plotted against Io's position in Jupiter's System III Longitude in Figure 14d, illustrating that such a correlation is not observed (R=0.29; p=0.49). Figure 14: Correlation between the total measured SO emission and various properties of Io: (a) Infrared brightness of Loki Patera. Note that the infrared observations were made at slightly different wavelengths, and the brightness at 3.8 µm is 1-2× the emission at 3.5 µm for typical temperatures of 300-500 K; (b) Io's mean anomaly; (c) Io's solar distance; and (d) Io's position in Jupiter's System III longitude. No significant correlation is seen in any plot, indicating that none of these factors individually is the dominant driver behind the total SO emission at 1.7 µm. Filled symbols correspond to dates when Io was moving into eclipse from sunlight, and open symbols to dates when Io was moving into eclipse from occultation.
We observed Io in Jupiter eclipse on three occasions in 2012-2016 and detected the 1.7-µm rovibronic band associated with the forbidden a 1 ∆ → X 3 Σ − transition of SO. On two of these dates, observations were made at a spectral resolution of R∼15,000, an order of magnitude improvement over all prior datasets. The high spectral resolution of these data permits a more detailed modeling treatment than has been possible in the past. Analysis of the spectra indicates a contribution from gas components at both high and low rotational temperatures. The high-temperature component is consistent with volcanic gas emission (1000-1500 K), while the low-temperature component is consistent with the temperature of Io's bulk atmosphere: 172-217 K during observations of Io just entering eclipse (in December 2015), and 76-119 K after Io had been in shadow for an hour (in May 2016). This is suggestive of a hot volcanic gas that has partially equilibrated rotationally to the bulk atmospheric temperature prior to de-exciting to the ground electronic state, a scenario which is consistent with the radiative properties of SO for a gas density of ∼ 10 11 cm −3 and a pressure of ∼1-3 nbar in the emitting region, if the atmosphere cools through the eclipse. These values are consistent with past observational and modeling work for Io's dayside temperature (e.g. Walker et al. 2010;Lellouch et al. 2007). This interpretation is challenged by the fact that Io's surface cools after entering eclipse over a timescale of ∼10 minutes (Tsang et al. 2016), while individual spectra taken throughout the first 30 minutes of eclipse do not show any systematic trends in gas temperature (though we note that if such a change were subtle, it would not be distinguished in the individual spectra due to low signal-to-noise). If the low-temperature SO emission is indeed probing the bulk atmospheric temperature, the atmospheric cooling must be delayed relative to that of the surface. In addition, the overall band strengths and best-fit gas densities are very similar between the two observing dates. This would seem to contradict recent evidence for the collapse of Io's bulk SO 2 atmosphere in eclipse (Tsang et al. 2016). However, we note that uncertainties on the inferred densities are high, and rely on several assumptions; a collapse by a factor of 5±2, as found by Tsang et al. (2016), is fully consistent with our results.
These three datasets bring the total number of detections of this emission band to eight. Analysis of the total band strength across all eight dates does not find any significant correlation between the observed SO emission and incident sunlight, Io's orbital phase, time since Io was last in sunlight, Jupiter's System III longitude, nor thermal hot spot activity. If the SO producing the emission band is indeed of volcanic origin, this indicates that thermal hot spot activity is not a good indicator of the gas emission.
The two-temperature models that provide a good fit to the high-resolution data are unable to reproduce an excess of emission near 1.69 µm that is seen in the low-resolution data, which cover a wider spectral range. Simple non-LTE models are able to match this emission by over-populating high rotational states, but a detailed analysis is limited by the coarse spectral resolution. Future observations with high spectral resolution across the entire 1.67-1.74 µm region will allow for a more in-depth characterization of the thermodynamic equilibrium state of the gas, which reflects on the gas origins and potentially the vent conditions. Continued observations with high spectral resolution after Io has been in shadow for differing time intervals would confirm or refute the tentative correlation between the temperature of the low-T gas component and the amount of time Io has been in shadow; if the correlation is confirmed, such studies would then provide a new way of studying how Io's atmospheric temperature and density drop at night. parameter. Because these simulations recover the probability distribution for each parameter under the modeling assumptions, the expectation of each parameter and its 1-σ uncertainty can be read from the distribution directly. However, in the case of asymmetric probability distributions, the 50th quantile may deviate from the maximum likelihood value. The 50th quantile is therefore not quoted, but the 16th and 84th quantile are given in Table 6 as the credible interval for each parameter.
Appendix A.1. Modeling the thermal continuum On each date, the continuum level in the vicinity of 1.7 µm is set by thermal emission from hightemperature volcanism, and its intensity and spectral shape vary between observations according to Io's volcanic activity. The low-resolution data from 2012 contain sufficient spectral coverage outside of the SO band to fit for the thermal background and subtract it from the data prior to the gas emission retrievals (see Figure 7). This is not the case for the high-resolution datasets from 2015 and 2016; for the data from these dates, the thermal continuum is fit simultaneously with the gas properties.
In all cases, the thermal background is modeled as a single emitting area at a single temperature; these models fit the data well, justifying the use of such a simplified formulation. The initial fits to the high-resolution data treat the emitting temperature and area (T th and A th ) as independent free parameters in the fits. However, on both dates the simulations find these parameters to be completely degenerate, indicating that for the best-fitting temperature regime, the minor changes in spectral shape between neighboring temperatures do not significantly affect the spectrum. The emitting area as a function of temperature is well-described by a fourth degree polynomial; after determining the polynomial coefficients for each dataset and model, the emitting area is directly calculated from the temperature in the fits, reducing the number of free parameters by one.   The joint probability distribution of the thermal continuum temperature and area from MCMC simulations using the 2T model and the Dec 2015 data. The parameters are found to be degenerate; the modeled A th (T th ) is shown in red; this curve is used to calculate the emitting area directly from the temperature, reducing the number of free parameters in the fits.