The Statistical Morphology of Saturn’s Equatorial Energetic Neutral Atom Emission

Saturn’s magnetosphere is an efficient emitter of energetic neutral atoms (ENAs), created through charge exchange of energetic ions with the extended neutral cloud originating from the icy moon Enceladus. We present an analysis using the complete image set captured by Cassini’s Ion Neutral Camera to characterize Saturn’s average ENA morphology. Concentric tori are formed around the planet by oxygen and hydrogen ENAs, with intensity peaks between 7 and 10 RS radial distance, with a ∼1–2 RS dayside offset. Nightside intensity is brighter than the dayside, likely the result of enhancements following large‐scale plasma injections from the magnetotail, and influence of the noon‐midnight electric field. Global intensity is clearly modulated with the near‐planetary rotation period. This Cassini‐era profile of Saturn’s ENA emission advances our understanding of how volcanic moons can influence plasma dynamics in giant magnetospheres and is timely ahead of the planned JUICE mission, which carries the first dedicated ENA detector to Jupiter.

Throughout the Cassini mission, INCA imagery has helped to reveal the global plasma dynamics throughout Saturn's magnetosphere. The neutral cloud (primarily hydrogen-and oxygen-based) extends out to at least ∼40 R S (Melin et al., 2009) and is sustained mainly by outgassing from Enceladus, but also the planet's rings and atmosphere. ENA production primarily occurs at radial distances from the outer edge of the E-ring/Rhea's orbit (∼8 R S ) out to Titan's orbit at ∼20 R S (e.g., Carbary et al., 2008a). Note the inherent dependence on the background neutral distribution means that peak ENA intensities are not necessarily coincident with peak ion intensities. A significant advantage of INCA is the global views of plasma circulation patterns it provides, complementing the partial picture obtained through in situ particle surveys of Voyager (Lazarus & McNutt, 1983) and Cassini (e.g., McAndrews et al., 2009;Sergis et al., 2007;Thomsen et al., 2014;Wilson et al., 2017).
In this study, we present a statistical analysis of Saturn's ENA distribution using all available Cassini INCA images. We developed an algorithm to ingest and calibrate the raw data and project it into the equatorial plane, detailed in the accompanying technical report by Bader, Kinrade, et al. (2021). Using multi-dimensional histograms, we are able to filter observations by pixel-specific parameters carrying information about observation geometry and target location in fixed and rotating reference frames, and create long-term mean averages.

Remote Sensing ENA Imagery
INCA was a time of flight detector on the Cassini orbiter which investigated the Saturn system for more than a decade. It was part of the Magnetosphere IMaging Instrument (MIMI) package whose scientific goals were to study the dynamics and configuration of Saturn's magnetosphere and its interaction with the solar wind and Saturn's moons (Krimigis et al., 2004). Toward these goals, INCA was capable of observing either ions or ENAs in keV energy ranges and determining their mass, energy, and direction of motion. Its large field of view covering 120 ° × 90 ° allowed it to measure significant parts of the ion pitch angle distribution when operated in ion mode Mitchell, Kurth, et al., 2009), or to perform global observations of ENAs in Saturn's magnetosphere when operated in neutral mode.
Over its 13+ years in orbit around Saturn, ENA observations were obtained from a variety of different perspectives-the resulting data set of ENA imagery is hence not homogeneous or easy to unify. We countered this difficulty by projecting all observations into Saturn's equatorial plane where most ENAs are created due to increased ion and neutral densities. The processing steps and the newly created data set are described in more detail in a dedicated technical report by Bader, Kinrade, et al. (2021). Characterization of the vertical ENA emission structure away from the magnetic/spin equator requires the use of a different projection routine and remains a possibility for future work.
In this study, we use the new data set of ENA projections to obtain longterm averages of the ENA intensity in Saturn's equatorial plane throughout the Cassini mission. This is done by calculating histograms across the entire data set; we bin the data by spatial location (X/Y KSMAG , radial distance, local time, PPO phase, …) and by the observed ENA intensity. The ENA intensity binning hereby covers differential ENA fluxes between 10 −15 and 10 4 particles/cm 2 /s/sr/keV in 380 logarithmic bins, resulting in a resolution of 20 bins per order of magnitude. Mean intensities are then calculated using the number of observations in each histogram bin. As described in detail in Bader, Kinrade, et al. (2021), many projections are unsuitable for use as the viewing geometry under which they were obtained leads to pixel stretching or significant loss of resolution. Here we chose to only take into account measurements from spacecraft lineof-sight distances between 5 and 30 R S from a given pixel, and from spacecraft elevations between 50° and 90° above a given pixel. Projections showing signs of data contamination (from sunlight or ion leakage in neutral imaging mode) were identified in Bader, Kinrade, et al. (2021) and excluded from this study. Figure 1 shows the resulting total exposure time across the region of interest within the above constraints. The view is in the Saturn-centered Kronocentric Solar Magnetic (KSMAG) frame, looking down on Saturn from above the north hemisphere with sun direction to the right. This distribution is largely a function of the varying Cassini orbit geometries throughout the mission; some pixel areas were observed for over 100 days in accumulation, with most of the projected coverage within 20 R S having at least 10 days exposure. Figure 2 shows the resulting maps of mean ENA intensity following our processing procedure. As with Figure 1, the view is of the KSMAG equatorial plane, as seen from above the north hemisphere, with the Sun to the right of each image. Panels 2a and 2b show the hydrogen 24-55 keV and 55-90 keV images, respectively, and panels 2c and 2d the oxygen 90-170 keV and 170-230 keV images, respectively. Note the logarithmic color scale mapping, and the shift in scale between 2a and 2b and 2c and 2d.

Mean Intensity Distributions
The toroidal morphology is striking in all four cases. This is formed by a region of low-level emission near to the planet within 5 R S (energetic ions are efficiently absorbed within these distances, e.g., Paranicas et al., 2008), and a ring of enhanced emission beyond this point which varies with energy (see Figure 3). The emissions are near-continuous in local time, with the intensity generally appearing brighter on the nightside than the dayside. A slight dayside offset of the tori by several R S is also perceptible, with the emission ring several Saturn radii further away from the planet than on the nightside (most clearly in the 24-55 keV hydrogen). The low-level emission within ∼5 R S may be attributed to a combination of emission from the planet (e.g., Mitchell, Krimigis, et al., 2009)  study. This is based on the sample distribution of projected pixels with LOS distance within 30 R S of the spacecraft and at least 50° inclination angle. Shown is the equatorial plane of Saturn as seen from above the north pole, with the Sun located to the right. While some regions near the edges of the projection have been observed for less than a day, most pixels include tens to hundreds of days of total exposure.

Local Time-Distance Profiles
To quantify the morphology further, Figure 3 unwraps the intensity maps of Figure 2 into a local time-distance frame. We have also applied Gaussian fits to each local time bin to determine the radial distance of the ENA intensity peak (black dotted line in each panel). These Gaussian fit parameters are provided in the Supplementary Information. First we see that the hydrogen ENAs peak at slightly larger radial distances (∼10-13 R S ) than the oxygen ENAs (∼7-10 R S ). These peak distance ranges are consistent with the earlier morphology study by Carbary et al. (2008a) who analyzed INCA imagery from 120 days in 2007 to reveal ENA intensity peaks at ∼11 R S and ∼8 R S for the hydrogen (20-50 keV) and oxygen (64-144 keV), respectively, albeit at slightly different energy ranges to those used here.
Second, we observe more clearly here the variation of peak distance with local time, that is, the dayside offset of the torus apparent in Figure 2. This is clearest in the hydrogen peak fits in images 3a and 3b, being around 10 R S on the nightside, but extending to 12-13 R S on the dayside (post-noon). The offset is less clear in the oxygen distributions. This apparent global shift of the torus toward the dayside could be associated with the fraction of an R S dawn-dusk trajectory shift of trapped energetic particles which has been associated with the global noon-to-midnight electric field, first discovered by looking at changes in moon wake locations (e.g., Andriopoulou et al., 2012;Thomsen et al., 2012). Wilson et al. (2013), analyzing in situ thermal plasma velocity and magnetic field measurements from Cassini, showed that the electric field is actually offset from the noon-midnight meridian by several hours local time, with the dusk-dawn drift velocity imposed on the global plasma circulation in the post-dawn direction out to at least 15 R S (Wilson et al., 2017). This offset has also been observed in other data sets (e.g., Andriopoulou et al., 2014;Roussos et al., 2019;Sun et al., 2019) and may be produced by variation in the radial force balance at different LTs (Jia & Kivelson, 2016). This may explain why the ENA tori are offset outwards toward post-noon, signifying the point at which particles begin to return inwards again in their displaced orbit. The dawn-dusk asymmetry is evident at distances out to ∼30 R S in the lower energy ENAs, suggesting that the E-field is effective KINRADE ET AL. at larger distances than previously considered. The offset appears to be less pronounced on the dayside for the higher energy bands, comparing the hydrogen peak fits in Figures 3a and 3b. This could be because the higher energy ions are less affected in their trajectory by the electric field asymmetry and thus show smaller orbital displacements, as discussed by Thomsen et al. (2012).
Average local time and radial distance profiles may be extracted from Figure 3, which is what we show in Figures 4a and 4b for each imaging energy range. Panel 4a shows the peak value of the average intensity for each local time bin, and panel 4b shows the local time-averaged radial ENA intensity profiles. Note the logarithmic intensity scale in each case. Two general morphological trends are apparent. The ENA intensity is brighter on the nightside than the dayside (by almost half an order of magnitude), and the mean tori distance decreases with species energy (there being a 1-2 R S shift between the respective hydrogen and oxygen energy peaks in panel 4b). The first of these trends is likely a consequence of periodic injections of energetic plasma from Saturn's nightside following magnetotail reconnection events which disperse as they rotate toward the dayside. The intensity asymmetry is also consistent with the effects of the noon-midnight electric field, as the particle energy and differential flux decrease adiabatically when transported outward on the dayside into regions of lower field strength and lower neutral density. Second, the decrease in tori radius with increasing particle energy may be related to the penetration depth of energetic ions and their effective charge exchange cross section which is a function of energy, that is, higher energy ions travel further inwards radially to a region of higher neutral density before producing an ENA (e.g., Paranicas et al., 2008). Shifting energetic ions and electrons to low L-shells is considered difficult because of gradient-curvature drift out of the injection flow channels. This result may indicate that the long-term ENA picture we see is dominated by large-scale plasma flows, wider azimuthally than those typically associated with interchange processes. Oxygen ions may also be penetrating the inner magnetosphere at higher charge states (subject to lower gradient-curvature drift), before breaking down to lower charge states and producing ENAs (e.g.  respectively. Note that different color scales apply to (a and b) and (c and d). For each local time bin, the radial intensity profile has been fitted with a Gaussian distribution to determine the radial distance of the ENA intensity peak (black dots). also be responsible for accelerating electrons to the inner magnetosphere, an effect which also extends to energetic ions.

Planetary Period Modulation
Periodicity in ENAs was observed early in the Cassini mission  and was first explained by Paranicas et al. (2005) as originating from a rotating point source imposed on a constant, global ring-type emission. Carbary, Mitchell, Brandt, Paranicas, and Krimigis (2008) then quantified the modulation periods in a Lomb-Scargle analysis of dawn-dusk ENA intensity variation using imagery taken from low-latitudes between 2004 and 2007. They found that 64-144 keV oxygen ENAs exhibited strong periodicity around 10.8 h, that is, near planetary-period (specifically the period of southern PPO system which was dominant at that time), while the 20-50 keV hydrogen periodicity varied between 8 and 13 h. We know from keogram track analysis of high-latitude imagery that most rotating ENA enhancements move with speeds around 60%-70% corotation (e.g., Carbary & Mitchell, 2014;Carbary et al., 2008b;Kinrade et al., 2020), which may explain the wider range of periodicities observed in the hydrogen emission. So how can we test for the near planetary-period modulation of global ENA emissions?
Any rotational modulation of the mean global ENA intensity should result from the interaction of both north and south PPO current systems as they superpose in the magnetic equator. The main effect of this interaction (in addition to north-south current sheet displacement) is modulation of the plasma sheet thickness, as the perturbation field vector periodically reinforces and weakens the core field depending on phases and relative strengths of the northern and southern systems (e.g., Bradley et al., 2018Bradley et al., , 2020  (c-f) Variation of the mean ENA intensity with planetary period oscillation (PPO) phase. Shown is the mean ENA intensity at radial distances between 6 and 12 R S on Saturn's dayside (6-18 h local time via noon) per PPO bin, relative to the overall mean ENA intensity in this region. Line colors indicate the different particle types and energy ranges as in panels (a and b); solid (dashed) lines indicate the intensity variation with northern (southern) PPO phase. the plasma sheet at PPO longitudes Ψ N = 0° and Ψ S = 180°, when reconnection in the tail (and therefore, possibly, subsequent ENA intensification following injections) is statistically more likely to occur (Bradley et al., 2018;Jackman et al., 2016). Conversely, a thicker plasma sheet would intuitively lead to more ions interacting with the neutrals to produce more ENAs. Such variations in the electron density (i.e., thick and dense vs. thin and tenuous) have been observed by the Cassini Langmuir Probe (Morooka et al., 2009). The thickness modulation is amplified when the north-south PPO systems superpose in antiphase and is reduced when the systems rotate in phase relative to one another.
Here we order all projected INCA pixels by their corresponding magnetic phase to reveal any periodicity present in the global ENA intensity. For a given exposure time, each pixel has a magnetic longitude with respect to both the northern and southern PPO phases, and throughout the Cassini mission the N-S phase systems swept through various periods of beat configuration and relative strengths. In Figures 4c-4f, we show the relative ENA intensities at each energy band within the main toroid region of 6-12 R S on Saturn's dayside (6-18 h local time), as a fraction of the mean intensity, versus the north and south magnetic phase. We chose to only investigate dayside observations to remove noisy intensity fluctuations arising from magnetotail dynamics. Note also that local time dependence is lost when binning the pixels with respect to the rotating phase systems.
The sinusoidal response visible in all panels of Figures 4c-4f indicates that the ENA intensities within the main toroid distances are being modulated within the frame of the rotating phase systems. The modulation is strongest in the 170-230 keV oxygen emission (4f), the intensity varying periodically between almost double (∼140°N, ∼315°S) and half (∼5°N, ∼175°S) the mean level for both phase systems. The 55-90 keV hydrogen (4d) and 90-170 keV oxygen (4e) emission responses also show clear modulation, but with a slight shift in phasing compared to higher energy oxygen. The 90-170 keV oxygen peaks at a northern phase of ∼225°, for example, compared to ∼140° for the higher energy band. This higher energy phase shift is in the right direction for gradient-curvature spread of the ions to be a contributing factor, but the magnitude of the shift (∼90°) seems too large for the energy difference here alone. The radial phase delay across these distances is only minor (∼30° from 9 to 15 R S , Andrews et al., 2010) such that we do not expect a strong smoothing effect from averaging across the 6-12 R S radial range. If we consider that the ENA intensifications may result from a more distant energization and injection process then it would be appropriate to consider the retarded phase at which they originated (e.g., Bradley et al., 2018). Such an investigation is beyond the scope of the current study. The key finding is the antiphase nature of these responses in the north and south rotating frames, indicating that the global ENA production is affected at some level by changes in plasma sheet thickness and density driven by the rotating perturbation fields as described above. The phases where the intensity peaks occur are all around ∼180°N, ∼0°S, indicating maximum ENA emission from the thicker, more dense plasma sheet sector. We are unable to explain here why the strength of the rotational modulation on the ENA intensity increases with energy, but this is an interesting result that warrants further study.

Conclusions
We have analyzed high-latitude Cassini INCA images of Saturn's equatorial ENA emission, using data from the entire mission. The result is the definitive Cassini-era picture of Saturn's equatorial ENA morphology obtained through remote sensing. The emissions are observed across four imaging bands; low and high energy hydrogen-derived ENAs (24-55 keV and 55-90 keV, respectively), and low and high energy oxygen (90-170 keV and 170-230 keV, respectively). We have used a new algorithm to sort, calibrate and project this ENA imagery into the equatorial plane for analysis, documented in Bader, Kinrade, et al. (2021). The final processed images are based on up to hundreds of days of accumulative exposure in some regions of the magnetosphere.
The long time average structure of the ENA emission is a torus around the planet, with intensity peaks at distances between 7 and 10 R S (the hydrogen emission peaking ∼1 − 2 R S further away from Saturn than the oxygen). This is consistent with both the inner boundary of the high particle pressure region between 8 and 14 R S (e.g., Sergis et al., 2007) and where the total neutral density starts to ramp up toward the inner magnetosphere within ∼8 R S (e.g., Cassidy & Johnson, 2010). The neutrals peak around 4 R S near Enceladus' orbit (e.g., Dialynas et al., 2013;Richardson et al., 1998), whereas inward-traveling energetic particles of up to several 100 keV start to drop in number around L-shells of 8-9 , so the ENAs are reflecting the ion population and not just the neutral density profile. ENA intensity is generally higher on the nightside than the dayside by almost half an order of magnitude, likely the result of persistent, largescale plasma injection activity manifesting in these long time average pictures, and the conservation of adiabatic invariants as plasma drifts around to the dayside. Another trend common to both species and energy ranges is the radial offset of the torus toward post-noon local times by up to 3 R S , which may be at least partly attributable to the trajectory shift experienced by trapped energetic particles under the influence of the global E-field asymmetry. Rotational modulation of the ENA emission is clearly present within the main tori distances (6-12 R S ), given the sinusoidal pattern of mean intensity variation with both north and south magnetic phases. The 170-230 keV oxygen ENAs exhibit the strongest modulation, varying periodically between almost double and half the mean emission level. This modulation effect decreases continuously with energy, being weakest in the 24-55 keV hydrogen, for reasons yet to be explained. The phase pattern is broadly consistent across all species/energies, and indicates that periodic modulation of the plasma sheet thickness could be driving these variations in ENA production-a thicker plasma sheet would intuitively lead to higher numbers of energetic ions interacting with the neutrals to produce more ENAs.
A fall-off in projection coverage at distances around 15-20 R S (see Figure 2)-where the main field-aligned currents associated with the PPOs are thought to close in the equatorial plane (Andrews et al., 2019)-prevents us from making a robust test of the ENA intensities further out in the magnetosphere. The main ENA tori peak at distances around the possible plasmapause-like boundary identified by Thomsen et al. (2015), a possible radial limit for injected hot plasma reaching the cooler, denser plasma of the inner magnetosphere. It may be that the modulation we observe here is driven by the combined effects of both plasma sheet thickness variation plus asymmetries in radial plasma flow associated with convection patterns in the inner magnetosphere, such as those identified in electron densities by Gurnett et al. (2007).
Characterizing Saturn's ENA distribution is useful for the development and constraint of chemistry models, and advances our understanding of how tiny volcanic moons can ultimately influence so much of the plasma dynamics in giant magnetospheres. This is particularly timely ahead of the planned JUICE explorer mission to Jupiter, which will fly with the first dedicated ENA detector to investigate the Jovian magnetosphere. Aside from being a proxy for ion loss processes, ENA imagery can be reverse-engineered to simulate the background neutral population. A question yet to be fully answered is how rotating ENA enhancements relate to counterpart auroral emissions, and what is the nature of the transient current system linking them during injection activity? Future work should investigate further the nature of the periodicity apparent in the global ENA intensity, how this relates to large-scale plasma injection dynamics, and if this is linked to magnetotail reconnection signatures observed preferentially at certain magnetic perturbation phases.