Emission of Energetic Neutral Atoms From the Magnetosphere‐Atmosphere Interactions at Callisto and Europa

We analyze the emission of energetic neutral atom (ENA) flux from charge exchange between Jovian magnetospheric ions and the atmospheres of Callisto and Europa. For this purpose, we combine the draped electromagnetic fields from a hybrid plasma model with a particle tracing tool for the energetic ions. We determine the ENA flux through a concentric sphere located just outside of each moon's atmosphere, thereby capturing the complete physics imprinted in these emission patterns. In order to constrain the modifications to the ENA emissions that arise from the periodic change of the ambient plasma conditions, we calculate the emission morphology at multiple positions during a Jovian synodic rotation. To isolate the influence of field line draping, we compare to the emission patterns in uniform fields. Our major results are: (a) At Europa and Callisto, the majority of detectable ENA emissions are concentrated into a band normal to the Jovian magnetospheric field. (b) The fraction of observable ENA flux that contributes to this band depends on the number of complete gyrations that the parent ions can complete within the moon's atmosphere. (c) Field line draping partially deflects impinging parent ions around both moons, thereby attenuating the ENA flux and driving significant morphological changes to the emission patterns. (d) The band of elevated ENA flux is locally maximized on the opposite (antipodal) side of the moon where the flux is locally minimized. At Europa, detectable ENA emissions are maximized slightly west of the ramside apex. At Callisto, they maximize near the Jupiter‐facing apex.


10.1029/2023JA031931
2 of 45 field becomes increasingly stretched in the radial direction from the growing magnetodisc contribution (e.g., Connerney et al., 2020).Accordingly, near Callisto, the field at times is largely parallel to the moon's orbital plane.However, the ambient magnetospheric field at Europa's orbit is still predominantly southward (Kivelson et al., 1999;Vance et al., 2021).
The periodic migration of Jupiter's magnetic equator sweeping past Callisto and Europa manifests as a time-varying horizontal component and a nearly constant north-south component of Jupiter's magnetospheric field at the orbit of each moon (Connerney et al., 2022;Kivelson et al., 1999).The time variability of the horizontal component generates an induction signal in conducting regions of both moons that possesses multiple excitation periods, the respective moon's synodic orbital period being the dominant cycle by about an order of magnitude (Kivelson et al., 1999;Seufert et al., 2011;Zimmer et al., 2000).Induced magnetic fields at Europa mainly stem from a salty subsurface ocean under its icy crust (Kivelson et al., 2000;Vance et al., 2021).Currents in Europa's ionosphere produce a secondary induction signal a small fraction (<1%) of the amplitude attributed to the oceanic induction signature (Hartkorn & Saur, 2017).Alternatively, the dominant conducting source region at Callisto, whether it be a subsurface ocean (e.g., Khurana et al., 1998;Zimmer et al., 2000), an anisotropic, conductive ionosphere (Hartkorn & Saur, 2017;Vance et al., 2021), or some combination of the two is under debate.The induced field can be approximately described by a dipole centered at each moon, with their magnetic moments antiparallel to the Jovian inducing component (Saur et al., 2010;Vance et al., 2018).
Europa and Callisto are constantly overtaken by a sub-corotating population of thermal magnetospheric plasma (e.g., Bagenal et al., 2015Bagenal et al., , 2016;;Kivelson et al., 2004), confined to a disc-shaped equatorial plasma sheet that is centered near the Jovian magnetic equator (e.g., Bagenal & Delamere, 2011).This plasma is incident upon both moons' trailing hemispheres with a relative bulk velocity on the order of 100 km/s (Bagenal et al., 2015;Kivelson et al., 2004).The ion population within Jupiter's plasma sheet is mostly comprised of protons, oxygen, and sulfur (Bagenal et al., 2015;Kim et al., 2020).The vertical and radial profiles of the plasma density within the sheet has been empirically modeled using Voyager and Galileo observations: the number density of the plasma sheet decreases with a Gaussian profile as a function of north-south distance to its center.
Ionization of Europa's and Callisto's neutral envelopes creates a population of pickup ions.The acceleration of these pickup ions by the magnetospheric fields drains momentum from the incident flow, slowing and diverting the magnetospheric plasma around each moon (e.g., Liuzzo et al., 2016;Rubin et al., 2015).As the impinging flow is slowed, the frozen-in magnetospheric field drapes around the moon and piles up at the ramside of either satellite.At Callisto, the large gyroradii of ionospheric pickup ions impose substantial hemispherical asymmetries on the perturbed electromagnetic field patterns (Liuzzo et al., 2015).At large distances from these moons' interaction regions, a system of Alfvén wings is generated that ultimately connect to the Jovian ionosphere (Neubauer, 1980(Neubauer, , 1998)).The presence of an induced dipole field reduces the cross sections of the Alfvén wing tubes and the strength of the wing-aligned currents, compared to a scenario without induction (Neubauer, 1999).
There also exists a population of magnetically confined energetic (E ≥ 1 keV) electrons and ions that simultaneously drift azimuthally through Jupiter's magnetosphere and bounce along field lines between mirror points (e.g., Kollmann et al., 2017;Paranicas et al., 2000;Paranicas et al., 2009;Shen et al., 2022).The influx pattern of the energetic ion population onto the surfaces of Callisto and Europa differs from that of the thermal population due to the larger gyroradii and significant field-aligned velocity components that grant the energetic ions increased access to both moons' polar regions (e.g., Addison et al., 2021Addison et al., , 2022;;Liuzzo et al., 2019b).Field line draping and the induced dipole significantly reduce and redistribute energetic ion precipitation patterns compared to uniform fields (Addison et al., 2021;Liuzzo et al., 2019bLiuzzo et al., , 2022;;Nordheim et al., 2022).At Callisto, the draping pattern protects the ramside from energetic ion precipitation in the keV range, while MeV ions can still access the moon uniformly.The induced dipole field tends to focus energetic ions into the regions where the axis of the dipole moment intersects Callisto's surface (Liuzzo et al., 2019a).When Europa is located near the center of the Jovian plasma sheet, field line draping acts to shield the ramside from heavy precipitation by magnetospheric ions.In this configuration, the draping is the strongest, and the energetic ion influx around the moon's ramside apex is minimized compared to large distances to the center of the Jovian plasma sheet (Addison et al., 2021).At high polar latitudes, the influx differs between the ramside and wakeside hemispheres by about a factor of two (Addison et al., 2021).The sensitivity of energetic ion trajectories at Europa and Callisto to the local field perturbations necessitates the inclusion of both the plasma interaction and the induced dipole when analyzing energetic ion dynamics in either moon's vicinity.
ENAs associated with the moons of the outer solar system have been observed both from Europa's neutral gas torus (Mauk et al., 2003(Mauk et al., , 2020) ) and as a result of the interaction between Saturn's magnetosphere and Titan's atmosphere (e.g., Brandt et al., 2012;Mitchell et al., 2005;Smith et al., 2009).At the time of this writing, Cassini is the only spacecraft to have visited the outer solar system with an ENA camera aboard (Krimigis et al., 2004;Mitchell et al., 1996).In late 2000, the Ion and Neutral Camera (INCA) imaged the Jovian magnetosphere as the spacecraft reached its closest approach distance of about 140R J .INCA observed ENAs in the (50-80) keV energy channel, primarily generated by charge exchange interactions with magnetospheric H + ions.The observed ENA emission pattern resembled an elliptical cloud possessing three distinct intensity maxima, one of which was produced by charge exchange with the Jovian atmosphere.The two weaker emission peaks, located at a distance of ≈10R J to the planet, were attributed to the limb-brightened ansae of a dense, trans-European neutral gas torus viewed edge-on by INCA (Mauk et al., 2003).
After advancing beyond the Jovian system, INCA observed ENA emissions during more than 100 separate Titan flybys, including the initial pass on 26 October 2004 (denoted TA).These in situ observations provided a plethora of ENA images of Titan's interaction with Saturn's outer magnetosphere (e.g., Dandouras & Amsif, 1999;Garnier et al., 2010;Mitchell et al., 2005).Similar to the Galilean moons, Titan's interaction with the impinging Kronian magnetospheric plasma results in strong magnetic field line draping and pileup (e.g., Neubauer et al., 2006;Simon et al., 2015).Wulms et al. (2010) made a first step in constraining the impact of the draped fields on the morphology of ENA emissions at Titan.These authors combined the perturbed fields from an MHD model of Titan's interaction with a test particle tracer to analyze the motion of energetic protons near Titan and their charge exchange with the moon's atmosphere.Modeling the ENA flux into a point-like detector, akin to INCA, is not feasible due to prohibitively few particle trajectories that would intersect the detector in the simulation.Therefore, Wulms et al. (2010) used a plane detector positioned at infinite distance to Titan that captured all ENAs approaching the detector approximately along its normal, partially emulating the limited field of view of the actual INCA detector.Wulms et al. (2010) modeled Titan's ENA emission morphology both with draped magnetospheric fields and with a uniform background field as a baseline to constrain the contribution of the moon's plasma interaction to the emission signatures.Their modeled ENA flux maps for both uniform and draped fields revealed a crescent-shaped intensity distribution interrupted by a broad gap with nearly zero ENA flux.However, the ENA emissions in draped fields were found to be maximized on the opposite side of Titan compared to the case of uniform fields; that is, the crescent is rotated about Titan's polar axis by ∼180° in longitude when draping is included, in agreement with the emission pattern observed during TA.Thus, the results of Wulms et al. (2010) highlight the strong dependence of the ENA emission morphology on local magnetic field perturbations near the source of the ENAs.Kabanovic et al. (2018) calculated the field perturbations at Titan during TA with the adaptive, kinetic ions, electron fluid (AIKEF) hybrid model (Müller et al., 2011).They then implemented an ENA tracing model similar to that of Wulms et al. (2010), and compared their results using the hybrid model to the ENA emission maps obtained for draped fields from the preceding MHD model from Wulms et al. (2010).In contrast to the MHD picture, the hybrid approach captures the asymmetries in the magnetic draping pattern at Titan that stem from the large gyroradii of ionospheric pick-up ions (Simon et al., 2007).Similar, large-scale asymmetries in the draped fields are also formed at Callisto (Liuzzo et al., 2015).When using the draped fields from the hybrid model and the same detector geometry as Wulms et al. (2010), the modeled ENA emissions were still found to possess a crescent-shaped intensity distribution.However, the broad emission gap that interrupts the crescent in MHD fields (Wulms et al., 2010) is partially filled in.Kabanovic et al. (2018) also demonstrated that their modeled ENA emissions are reduced in magnitude by about a factor of two when the magnetospheric background field near Titan is increased from 3.77 to 6.0 nT.The emission gap captured by their plane detector was found to grow with ambient magnetospheric field strength because the decrease in gyroradii is accompanied with a reduced accessibility of Titan's atmosphere to energetic ions (see also Regoli et al., 2016).Ergo, even rather subtle changes to the draped magnetic fields leave a clearly discernible imprint in the ENA emission pattern at Titan.
Recently, Tippens et al. (2022) studied the impact of varying magnetospheric field orientations on Titan's ENA emission morphology using a combination of the AIKEF hybrid model and a test particle tracer, similar to Kabanovic et al. (2018).These authors employed a spherical ENA detector that encapsulates all of Titan's atmosphere.Such a detector geometry allows a comprehensive analysis of observable ENA generation at Titan without truncating any physical effects: due to its point-like nature (compared to the scale of the plasma interaction) and limited field of view, an actual spacecraft detector would capture only a small fraction of the ENA flux emanated from Titan's atmosphere.Tippens et al. (2022) found that the majority of emitted ENAs are focused into an equatorial band of high flux perpendicular to the magnetospheric background field orientation.The latitudinal extension of this band is determined by the maximum altitude below which the bulk of the detectable ENA production takes place.Tippens et al. (2022) further established that ENA emissions through any segment of a global, spherical detector carry information from the entirety of the interaction region and are not necessarily determined by the regions of the atmosphere in closest proximity.In addition, these authors demonstrated that field line draping around Titan significantly reduces the overall intensity of the ENA emissions.
The Jupiter Energetic Neutrals and Ions (JENI) detector in the Particle Environment Package on the JUpiter ICy moons Explorer (JUICE) spacecraft functions as an ENA camera (Grasset et al., 2013), the energy range of which is from around 1 keV up to several 100 keV (Galli et al., 2022).Charge exchange cross sections for protons with molecular oxygen (the primary constituent of Europa's and Callisto's atmospheres) fall precipitously at incident ion energies above 100 keV (Lindsay & Stebbings, 2005), suggesting that most ENAs emanating from these two moons' atmospheres will be observable between about 1 and 100 keV.The JUICE spacecraft is scheduled to target both Europa and Callisto for multiple flybys (e.g., Grasset et al., 2013;Liuzzo et al., 2018), and may capture ENA images at both moons in analogy to Cassini's persistent ENA imaging of Titan's interaction region.
The Europa neutral torus that was markedly apparent in the INCA observation of the Jovian system is sourced by escaping hydrogen from Europa's atmosphere (Smith et al., 2019).Due to INCA resolution limitations and the large passing distance during Cassini's flyby, there were no distinct observations of ENA emissions from Europa's local interaction region.Europa's neutral torus is approximately 10 to 12 orders of magnitude less in peak number density than the atmosphere at the moon's surface (Plainaki et al., 2018;Roth et al., 2023;Smith et al., 2019), suggesting that Europa's atmosphere may serve as a highly localized source of intense ENA emissions.Callisto's magnetospheric environment shares strong similarities with the conditions at Titan, including pickup ion gyroradii on the order of the respective moon's diameter and an occasionally super-alfvénic upstream flow that gives rise to highly draped fields (Liuzzo et al., 2015;Simon et al., 2015).The atmospheric scale heights at Callisto and Titan are comparable (e.g., Carberry Mogan et al., 2020;Tippens et al., 2022;Vorburger et al., 2015, and references therein).In addition, the density of Titan's upper atmosphere (above altitudes of 1, 100 km) where observable ENAs are generated is comparable to the atmospheric density at Callisto (Liuzzo et al., 2015;Tippens et al., 2022).The upstream energetic proton flux at keV energies is 1-2 orders of magnitude stronger at Callisto than at Titan (e.g., Liuzzo et al., 2022;Mauk et al., 2004;Regoli et al., 2018).This suggests that intense, observable ENA emissions may be emanating from Callisto's atmosphere.
An increase in the ambient magnetospheric field strength by about a factor of two at Titan already produced an extended ENA emission gap in the Saturn-averted hemisphere associated with the reduced gyroradii of the energetic parent ions (Kabanovic et al., 2018).Periodic variation of Callisto and Europa in magnetic latitude relative to the center of the Jovian plasma sheet is accompanied by similar (Europa) or even much stronger (Callisto) changes in the magnitude of the background field (Connerney et al., 2022;Kivelson et al., 1999Kivelson et al., , 2004)).Analogous to Titan, this effect is expected to generate variations in ENA emission intensity and morphology.
To date, modeling of local ENA emissions at Europa and Callisto has not been carried out by any study.In light of the upcoming JUICE mission's ENA detection capabilities, we combine the AIKEF hybrid code (Müller et al., 2011) with a particle tracing tool to model global ENA emissions from Callisto's and Europa's atmospheres.In particular, our study seeks to address the following questions: • What is the effect of magnetic field line draping on the intensity and morphology of ENA emissions originating from the atmospheres of Europa and Callisto?• How does the variability in the moons' plasma environments during a Jovian synodic rotation drive changes in the global ENA emission signatures?• How do ENA emission intensity and morphology at these moons vary as a function of magnetospheric parent ion energy?• Do the anisotropies observed in the pitch angle distribution (PAD) of energetic ions at Europa (Sarkango et al., 2023) map into the ENA emission patterns?• Does Europa's localized water atmosphere, as recently identified by Roth (2021), leave a discernible imprint in the ENA emission morphology?
This work is structured as follows: the employed hybrid plasma and ENA particle tracing models are described in Section 2. In Section 3.1, we present the modeled perturbed field environments at each moon, as obtained from AIKEF.Section 3.2 discusses ENA emission maps generated for different electromagnetic field configurations at Europa.In Sections 3.3 and 3.4, we investigate how anisotropies in the ambient protons' PAD or a local water vapor "bulge" in the atmosphere affect the ENA emission morphology at Europa.Our results for ENA emissions at Callisto in uniform and draped fields are presented in Section 3.5.A summary of analysis and conclusions drawn can be found in Section 4.

Modeling ENA Emissions
In Section 2.1, we introduce the setup of the AIKEF hybrid model (Müller et al., 2011) applied to calculate the three-dimensional structure of the draped fields at each moon.Section 2.2 then describes our ENA generation model.Energetic protons are traced through the perturbed electromagnetic fields from AIKEF in each moon's local environment.As these ions interact with Callisto's or Europa's atmosphere, they undergo charge exchange and produce ENAs.These ENAs propagate along rectilinear ballistic trajectories and may be detectable by a spacecraft.The same combination of AIKEF model and a particle tracing tool to calculate ENA emission patterns has already been applied at Titan (Kabanovic et al., 2018;Tippens et al., 2022).
At both moons, the right-handed Cartesian Satellite Interaction System with coordinates r = (x, y, z) is employed, where the origin lies at the center of the moon.The corresponding unit vectors of the basis are  ( x𝐱, ŷ𝐲, ẑ𝐳) .Here,  ẑ points northward and is aligned with the Jovian spin axis, and

𝐴𝐴
ŷ bears toward Jupiter.The unit vector  x completes the system and points in the direction of corotation.To describe the ENA flux emanating from various regions around the moons, we introduce the spherical West Longitude coordinate system.Longitude is defined as 0° at the moon's sub-Jovian apex, and increases westward.The downstream, anti-Jovian, and upstream apices are therefore located at 90°W, 180°W, and 270°W, respectively.Latitude is given with 90° corresponding to the moon's geographic north pole and −90° to its south pole.

Modeling the Perturbed Electromagnetic Environments of Callisto and Europa
We apply the AIKEF hybrid model (Müller et al., 2011) to calculate the three-dimensional structure of the perturbed electromagnetic fields near each moon.AIKEF has been utilized in numerous preceding studies of both Callisto (Liuzzo et al., 2015(Liuzzo et al., , 2016(Liuzzo et al., , 2017(Liuzzo et al., , 2018(Liuzzo et al., , 2019a(Liuzzo et al., , 2019b(Liuzzo et al., , 2022) ) and Europa (Addison et al., 2021(Addison et al., , 2022 quantitatively replicated fine structures in Galileo magnetometer data from the Europa E26 flyby that were consistent with the presence of a plume (Arnold et al., 2019).Output from the AIKEF model has also been found to agree well with magnetic field and plasma signatures observed by Galileo during Callisto flybys C3, C9, and C10 (Liuzzo et al., 2015(Liuzzo et al., , 2016)).The AIKEF model treats thermal ions as macroparticles and electrons as a massless, charge-neutralizing fluid.Therefore, it is capable of resolving asymmetries associated with the ionospheric Hall effect or flow shear between magnetospheric and ionospheric species, as well as any asymmetries in the plasma interaction due to the large pick-up ion gyroradii at Callisto (e.g., Liuzzo et al., 2015Liuzzo et al., , 2016)).For a detailed description of the computations performed by AIKEF, the reader is referred to Müller et al. (2011) as well as any of our preceding publications on Callisto and Europa.The three-dimensional AIKEF data sets for Callisto that we use in this study have also been applied by Liuzzo et al. (2022) in a companion study of energetic particle precipitation onto the top of that moon's atmosphere.
The extent of the AIKEF simulation domain is chosen such that the electromagnetic fields have returned to their undisturbed background values at its outer boundaries, except for two highly-localized regions within the Alfvén wing tubes.For Callisto's environment, Liuzzo et al. (2022)  We seek to constrain the impact of different sets of magnetospheric upstream conditions on the shape and intensity of ENA emissions, so we consider two model configurations at each of the two moons.We implement ambient plasma conditions that Europa or Callisto encounter at the center of the Jovian plasma sheet and at the maximum distance below the center of the plasma sheet.The background magnetic field vectors B 0 in each case are enumerated in Table 1.At Europa, the background field vectors are determined by a combination of the latest Juno magnetodisc and internal Jovian field models (Connerney et al., 2020(Connerney et al., , 2022).Callisto's orbit falls beyond are computed for thermal magnetospheric ion (s = i) and electron (s = e) species separately using temperatures T s from Kivelson et al. (2004), where subscript s represents species.The Alfvén speed |v as well as Alfvénic (M A ), sonic (M S ), and magnetosonic (M MS ) Mach numbers are also given.

Table 1
Adaptive, Kinetic Ions, Electron Fluid (AIKEF) Model Parameters the realm of applicability of the Juno magnetodisc model.Therefore, the hybrid model setup adopted from Liuzzo et al. (2022) implements a combination of the VIP4 model for Jupiter's internal field (Connerney et al., 1998) superimposed with the Khurana (1997) model of Jupiter's magnetodisc field.
Near the center of the Jovian plasma sheet, the horizontal component of the ambient magnetic field is weak compared to the southward B z,0 component at both moons (Connerney et al., 2020(Connerney et al., , 2022;;Kivelson et al., 1999).At Europa's orbit, the observed |B h | at the center of the plasma sheet is about five times weaker than the southward component (Kivelson et al., 1999).When Callisto is located at the center of the sheet, |B h | is weaker than 1 nT compared to the southward component B z,0 = −4 nT.Thus, at the center of the plasma sheet, the tilt of the background field against the north-south direction at either moon does not exceed arctan(1/4) ≈ 14°.To maintain a diagnostically accessible geometry that allows straightforward analysis of the physics of ENA generation, we do not take into account the weak B h component for either moon at the center of the plasma sheet (see Table 1 and Liuzzo et al., 2022).In other words, we consider the magnetospheric field at the center of the sheet to be purely southward.As shown by Tippens et al. (2022) for Titan, the tilt of the ambient magnetospheric field vector would merely cause a corresponding rotation (by <14°) of the ENA emission pattern against the z axis, but would not affect the physical mechanisms shaping the ENA flux distribution.
While we do consider the plasma interaction of each moon at maximum distance below the center of Jupiter's plasma sheet, we do not present a dedicated analysis of the corresponding scenario at maximum distance above the sheet.Between these two scenarios, only the signs of the B x,0 and B y,0 components flip; the field magnitude remains largely unchanged (Bagenal et al., 2015;Connerney et al., 2020Connerney et al., , 2022;;Kivelson et al., 1999).When moving from above to below the plasma sheet, certain features of the interaction (e.g., the pick-up tail or the magnetic pileup region) are relocated, but the involved physics does not change (e.g., Simon & Motschmann, 2009).Furthermore, the thermal plasma density is symmetric about the center of the plasma sheet, so the upstream density is identical for equivalent positions above and below the center (Bagenal & Delamere, 2011;Roth et al., 2014a).Thus, we expect the morphology of ENA emissions from either moon to possess the same key features at maximum distance above or below the center of the sheet.
When Callisto and Europa are at their maximum distance from the center of the Jovian plasma sheet, the moons' inductive responses are the strongest (assuming no phase lag).We follow numerous preceding studies and assume the inductive response to stem from a spherically symmetric, infinitely conducting region of the same size as the respective moon.Such an approach was found to reproduce key features of Galileo magnetometer observations at both Europa (e.g., Arnold et al., 2019;Harris et al., 2021;Zimmer et al., 2000) and Callisto (e.g., Lindkvist et al., 2015;Liuzzo et al., 2015;Rubin et al., 2015).Under these conditions, the induced magnetic dipole moment M ind at either moon is centered at its core and given by where R is the moon's mean radius (i.e., R = R E at Europa and R = R C at Callisto).A more rigorous treatment of the inductive response would include for example, finite conductivity of the ocean layer, radial gradients in ocean temperature that map into the conductivity profile, as well as an anisotropic geometry of the conducting layer (e.g., Schilling et al., 2007;Seufert et al., 2011;Styczinski & Harnett, 2021;Vance et al., 2021).However, as we will show in Section 3, the induced dipole moment plays only a very subtle role in shaping the ENA emission morphology at Europa and Callisto; that is, we expect our results to be very robust against minor changes of this parameter.
The bulk velocity u 0 of the incident magnetospheric flow at Europa  (0 = 100 km∕s x) is approximately 85% of rigid corotation.At Callisto, a value of  0 = 192 km∕s x has been found to replicate Galileo magnetometer observations (Liuzzo et al., 2016(Liuzzo et al., , 2017)).This is only approximately 60% of the local corotation velocity at Callisto's orbit (Kivelson et al., 2004).At Europa, the distribution of the thermal ions is given by a mixture of protons, atomic oxygen, and doubly-ionized atomic sulfur (Bagenal et al., 2016), with a mean upstream ion mass m 0 ≈ 18.5 amu (Clark et al., 2020).The AIKEF simulations applied by Liuzzo et al. (2022) for Callisto include a thermal ion population that consists exclusively of atomic oxygen that is, m 0 ≈ 16.0 amu (e.g., Kivelson et al., 2004).At both moons, the Gaussian profile, 10.1029/2023JA031931 8 of 45 describes the variation in the number density n 0 of the upstream thermal plasma as a function of vertical distance d to the center of the Jovian plasma sheet, with equatorial density n ps,0 and vertical scale height of the plasma sheet H ps .The parameters n ps,0 and H ps vary as a function of radial distance to Jupiter.Galileo observations from a large number of crossings near Europa's orbit, obtained by the Plasma Subsystem (PLS), suggest n ps,0 = 50 cm −3 with a scale height of H ps ≈ 1.8R J (Bagenal & Delamere, 2011).However, Roth et al. (2014a) found that a base density value of n ps,0 = 200 cm −3 with a vertical scale height of H ps ≈ 0.9R J is most suitable to reproduce Galileo Plasma Wave Subsystem observations of the upstream density from the spacecraft's targeted Europa flybys (Kurth et al., 2001).We therefore apply the values of n ps,0 and H ps from Roth et al. (2014a) to calculate the upstream density for our Europa model.At Callisto's orbit, the thermal plasma density at the center of the sheet has fallen to n ps,0 = 0.15 cm −3 and the scale height H ps has enlarged to 3.7R J , in accordance with Galileo PLS observations (Bagenal & Delamere, 2011;Liuzzo et al., 2022).
We consider the atmospheres of both moons to be dominated by molecular oxygen, with a minor contribution from carbon dioxide at Callisto, consistent with modeling results and Hubble Space Telescope (HST) observations (e.g., Carberry Mogan et al., 2021;Cunningham et al., 2015;Liuzzo et al., 2015;Plainaki et al., 2018;Saur et al., 2011).While additional traces of hydrogen are also present at both moons (e.g., Carberry Mogan et al., 2020;Roth et al., 2017b), these light particles do not appreciably drain momentum from the magnetospheric flow and therefore have negligible effect on the plasma interaction (see also Simon et al., 2007).At both Europa and Callisto, our treatment of the atmosphere includes an asymmetry in the neutral density profile such that the atmospheric density peaks at the ramside apex and is minimized at the wakeside apex, consistent with the models of Rubin et al. (2015), Liuzzo et al. (2015Liuzzo et al. ( , 2016)) (3) The model of Liuzzo et al. (2022) uses a surface density of  O 2 ,0 = 10 16 m −3 and a scale height  O 2 , = 230 km to describe the molecular oxygen component.For the CO 2 constituent, these authors employ a surface density of  CO 2 ,0 = 4 ⋅ 10 13 m −3 and use  CO 2 , = 230 km .Their chosen parameters yield column densities for each species consistent with the observed (Carlson, 1999;Cunningham et al., 2015;Roth et al., 2016) as well as the modeled values (Vorburger et al., 2015).For Europa's molecular oxygen atmosphere, we employ the asymmetric density profile implemented by Arnold et al. (2019), as well as Arnold, Liuzzo, and Simon (2020), and Addison et al. (2021, 2022): Using a surface density of n n,0 = 5 ⋅ 10 13 m −3 and a scale height of h s,E = 100 km, this atmospheric model is concordant with HST observations of the molecular oxygen column density (Hall et al., 1995;Plainaki et al., 2018;Roth et al., 2014b), as well as estimates of the molecular oxygen column density utilizing Keck/HIRES observations of Europa's auroral emissions during Jupiter eclipses (de Kleer et al., 2023).The density of Europa's dilute hydrogen corona is four orders of magnitude below that of the molecular oxygen component (Roth et al., 2017b).Therefore, we include it neither in AIKEF nor in our ENA generation model.
Recently, a persistent H 2 O atmosphere above Europa's subsolar trailing apex was identified in HST observations (Roth, 2021).This water vapor atmosphere is estimated to have a surface density comparable to that of the molecular oxygen constituent.However, it is tightly confined around Europa's subsolar apex (Cervantes & Saur, 2022).The number density profile of this additional component can be expressed as where φ is the angle between a vector from the moon's center to the subsolar apex and the vector r from the center of the moon to a certain point in the atmosphere.We use γ = 6, in analogy to Cervantes and Saur (2022), who found this value to be consistent with HST observations of the H 2 O profile.In addition, we set the scale height to  H 2 O = 138 km , and the surface density at the subsolar apex to  0,H 2 O = 2.14 ⋅ 10 14 m −3 , consistent with the column density derived by Roth (2021).The ratio of the total numbers of H 2 O to O 2 molecules in our model atmosphere is approximately 0.23, consistent with the latest constraints provided by observations of the moon's optical aurorae by de Kleer et al. ( 2023): these authors found a ratio of 0.25.
To constrain a possible contribution of the H 2 O component to ENA emissions, we consider a single test case for Europa (situated at the center of the Jovian plasma sheet) with the inclusion of this profile in addition to the molecular oxygen component.Europa's subsolar apex and ramside apex were coincident during the occultation observations that inferred the presence of the H 2 O exosphere (Roth, 2021), leaving it ambiguous whether the H 2 O bulge follows the subsolar point as Europa orbits Jupiter or remains concentrated around the trailing apex regardless of the moon's orbital position.We thus position the H 2 O exosphere in the same configuration as it was observed, with the subsolar and ramside apices coincident.
At Europa, electron impact ionization dominates photoionization by at least an order of magnitude (Saur et al., 1998), so we include only this process to generate the moon's ionosphere from its neutral envelope.Carberry Mogan et al. (2023) recently estimated a range of 3.3 ⋅ 10 −6 to 13.6 ⋅ 10 −6 s −1 for the rate of electron impact ionization of H 2 O at Europa.This range of rates largely overlaps with that of the rates for electron impact ionization of O 2 at Europa (1.66 ⋅ 10 −6 -8.5 ⋅ 10 −6 s −1 , see Carberry Mogan et al. (2023)).Therefore, we proceed analogous to Arnold et al. (2019) and Cervantes and Saur (2022), who use the same cross sections for O 2 and H 2 O.At Callisto, photoionization plays the primary role in ionospheric generation (Hartkorn & Saur, 2017), mandating its inclusion in the hybrid model.In addition, AIKEF takes into account an isotropic contribution from electron impacts to emulate Callisto's nightside ionosphere (e.g., Liuzzo et al., 2015Liuzzo et al., , 2022)).The dayside and ramside apices are set to be coincident at Callisto for this study (see also Liuzzo et al., 2022).

ENA Emission Model
While the AIKEF model can, in principle, resolve the motion of energetic ion macroparticles, the velocities of particles in the keV regime are more than an order of magnitude greater than the bulk velocity of the thermal ions.Since the timestep in AIKEF is limited by the Courant-Friedrichs-Lewy condition, treating thermal and energetic ion populations simultaneously would require a prohibitively minuscule timestep in the hybrid model.Near the orbits of Europa and Callisto, energetic ions carry only a small portion of the plasma's current density (Kim et al., 2020;Mauk et al., 2004).Thus, their weak contribution to the currents can be neglected when calculating the electromagnetic field perturbations.Instead, the energetic ions can be treated as test particles injected into the steady-state field configuration from AIKEF (analogous to, e.g., Addison et al., 2021Addison et al., , 2022;;Breer et al., 2019;Liuzzo et al., 2019aLiuzzo et al., , 2019bLiuzzo et al., , 2022)).
To determine the impact of field line draping on the morphology of ENA emissions at Europa and Callisto, we model ENA production in both the uniform magnetospheric background field B 0 and in the draped fields at each moon, similar to Wulms et al. (2010), Kabanovic et al. (2018), andTippens et al. (2022) for Titan.In the uniform case, the electric field is equivalent to the convective term, E 0 = −u 0 × B 0 .We consider monoenergetic parent proton populations at discrete energies of E ∈ {1, 10, 40, 70, 100} keV.Above 100 keV, the value of both the charge exchange cross section as well as the ambient energetic proton fluxes at each moon decay rapidly (Cooper et al., 2001;Lindsay & Stebbings, 2005;Paranicas et al., 2002).Thus, this energy range covers the segment of the energetic ion distribution near the moons that makes the strongest contributions to the observable ENA flux.The gyroradii r g of protons in the considered energy regime are provided in Table 2.
The setup of our ENA generation model is displayed in Figure 1, using Callisto as an example.Parent ions are launched from the nodes of a starting grid that is attached to the outer faces of the AIKEF simulation domain, illustrated by the white lattice in Figure 1.This ensures that even for the highest proton energies considered (100 keV) there is still a distance of at least 2.5r g between any plane of the starting grid and the Callisto (see Table 2).Due to the 10-100 times stronger magnetic field at Europa, the distance between the moon and the faces of the starting grid (in units of proton gyroradii) is even larger.The directions from which energetic protons can impinge upon the moon's atmosphere are dependent upon the orientation of B 0 .Velocities  √ 2∕ (where m p is proton mass) of protons within our considered energy regime range from approximately 440 to 4, 400 km/s.As such, the corotation drift velocity (on the order of 100 km/s) provides only a minor contribution to the net velocity, and particles with pitch angles near 0° and 180° will move along magnetic field lines much faster than their drift in the +x direction.
Energetic protons bouncing between mirror points in the Jovian magnetosphere can generate ENAs as they enter Europa's or Callisto's interaction region and interact with the atmosphere.To incorporate such particles in our model scenario, we include faces of the starting grid at the northern (z > 0) and southern (z < 0) as well as the Jupiter-facing (y > 0) and Jupiter-averted (y < 0) boundaries of the AIKEF domain.In general, parent ions launched from the northern/southern and Jupiter-facing/Jupiter-averted faces of the starting grid can potentially reach the moon's atmosphere and produce ENAs whenever the background field B 0 possesses a component orthogonal to the respective face.Energetic parent ions are also launched from the downstream face (x > 0) of the starting grid.Note that protons initiated at the downstream face can reach either moon by traveling along the magnetic field lines only if the background field possesses a component normal to that face.Otherwise, the corotation drift will carry these protons away from the moon and they cannot travel toward upstream.We refer the reader to Section 2.3 of Tippens et al. ( 2022) for further discussion of the starting grid and how it captures bouncing protons capable of interacting with either moon's atmosphere.In Figure 1, the upstream, Jupiter-facing, and southern planes of the starting grid are depicted.
Each face of the starting grid is populated with nodes that are equally spaced along the x, y, and z directions, forming squares of area dA.At each of these nodes we launch energetic protons of a given initial energy, a sample of which is indicated by the intersections of the white lines in Figure 1.If energetic proton gyroradii greatly exceed the moon's atmospheric scale height, these ions can only complete partial gyrations within the moon's atmosphere when initialized at a certain phase of their gyration.At Europa, energetic proton gyroradii are smaller than or comparable to the atmospheric scale height h s,E (see Table 2).For this reason, parent protons can complete multiple gyrations within the atmosphere (and produce ENAs), irrespective of the phase of its gyration when initialized on the starting grid.However, at Callisto, energetic proton gyroradii range from about 0.02R C to 4.7R C within the considered energy range (see Table 2).Thus, depending on their initial gyrophase and the resolution of the starting grid, a large fraction of the newly initiated protons may miss Callisto's atmosphere entirely, even if the trajectories of their guiding centers come close to the moon.We found that this may lead to "blurring" of our modeled ENA emission maps.In order to obtain a sufficiently smooth ENA emission pattern, the starting grid at Callisto therefore requires a higher spatial resolution than at Europa where energetic ions can complete many gyrations within the atmosphere.
The length scales of Callisto's interaction region are comparable to those at Titan, in particular regarding the size of ion gyroradii relative to the diameter of the moon and the spatial extension of the electromagnetic field perturbations (in units of the respective moon's radius), see, for example, Liuzzo et al. (2015) and Simon et al. (2015).Accordingly, at Callisto we define a separation of 0.04R C between adjacent nodes of the starting grid in each direction, similar to Wulms et al. (2010), Kabanovic et al. (2018), andTippens et al. (2022).Since the starting grid of the ENA generation model is attached to the outer faces of the AIKEF domain, this yields  (30 ∕0.04 ) 2 = 562, 500 unique positions of initialized ions on each of the 6 faces at Callisto.At Europa, a coarser starting grid resolution of 0.15R E is sufficient for a smooth distribution of the modeled ENA fluxes.
At each starting grid node, we initialize a population of parent protons whose velocity vectors isotropically cover a sphere in velocity space.The velocity sphere is discretized in zenith and azimuth using polar coordinates ϕ v and θ v , with resolution on the velocity sphere Δθ v = Δϕ v = 5°.We initialize each particle with a velocity magnitude where r is the position of the particle.

Table 2 Energetic Proton Gyroradii in Units of Respective Moon Radius R (R = R E at Europa and R = R C at Callisto)
Each newly initialized parent ion is "weighted" with a flux value derived from angularly resolved proton spectral energy distributions  Ĩ() observed by the Galileo Energetic Particle Detector (EPD) in the uniform plasma environment outside of each moon's interaction region.At Europa, we employ energetic proton distributions  Ĩ() taken during the close Galileo flybys at multiple distances to the center of the Jovian plasma sheet (e.g., Cooper et al., 2001;Paranicas et al., 2002Paranicas et al., , 2009)).For the center of the sheet, we use the fit provided for the observed EPD spectrum from E12.The E12 flyby occurred at a distance of only d ≈ 0.1R J (see Equation 2) to the center of the sheet.For model setups with Europa located at maximum distance below Jupiter's plasma sheet, we implement the fit to the ambient  Ĩ() from the E26 flyby which took place at a distance of d ≈ 1.0R J below the sheet.
When Callisto is at the center of the Jovian plasma sheet, we employ the fit to EPD data obtained from a Galileo sheet crossing near the moon's orbital distance (Mauk et al., 2004, see their Table 1, label G8 PS/A).Callisto was not proximal to the spacecraft during these observations, so this  Ĩ() spectrum provides an adequate representation of the ambient magnetospheric proton flux.Callisto's M shell value (e.g., Liuzzo et al., 2019bLiuzzo et al., , 2022;;Szalay et al., 2017) changes by almost a factor of 3 going from the center of Jupiter's plasma sheet (M = 26.3) to the maximum distance below (M > 70).The intensity  Ĩ() of the ambient energetic proton flux has been suggested to fall by about an order of magnitude moving from the center to maximum distance from the center of the Jovian plasma sheet (Kollmann et al., 2018;Liuzzo et al., 2022;Paranicas et al., 2018).When Callisto is situated at maximum distance below Jupiter's plasma sheet, we proceed analogous to Liuzzo et al. ( 2022) and use the A single parent ion trajectory (light blue) is displayed, initiated at a grid node on the upstream face of the starting grid.This parent ion remains within Callisto's atmosphere (orange shaded region) for nearly an entire gyration and generates ENAs.ENAs that collide with the moon (red) are not detected, but ENAs with a velocity vector bearing away from the moon (green) are observed by the spherical model detector that encompasses the moon's entire atmosphere.The figure is not to scale.energetic proton distribution observed at the center (Mauk et al., 2004), but downscaled by a factor of 10 at all energies.
The energy spectrum  Ĩ() represents the number flux of particles at a certain energy per steradian (e.g., Paranicas et al., 2009).In the classical regime, this quantity is related to the distribution function f of the proton population with speed  = | ̇| in phase space by see, for example, Kollmann et al. (2018).The PAD included in our model was determined from in-situ spacecraft observations of energetic particle fluxes near Europa and Callisto.The observed particle fluxes  Ĩ() and the PAD can be related to the total energetic proton flux I(E), integrated across all viewing directions, through where p(α(θ v , ϕ v )) is the PAD, α denotes the pitch angle, and dΩ v ≡ sin θ v dθ v dϕ v (see also Addison et al., 2021).
The initial weight of each proton macroparticle is given by   = Ĩ()dΩ = Ĩ() sin ΔΔ for initial energy E when launched at the starting grid.The angular factor sin θ v takes into account that the velocity vectors of newly initiated protons are closer to each other near the poles of the velocity sphere than when near its equator.
Galileo observations of the ambient proton fluxes at Callisto's orbit suggest an isotropic PAD in the energy range we consider (Mauk et al., 2004;Paranicas et al., 2002).In addition, Juno observations revealed a nearly isotropic proton PAD near Callisto at energies of 117 keV, that is, slightly above our considered energy range (Shen et al., 2022).A recent survey of Juno data (Sarkango et al., 2023) suggests that the proton PAD inward of Callisto's orbit may have a nonuniform shape at energies below 200 keV; that is, the proton fluxes minimize around pitch angles of α ≈ 90° and maximize at α = 0° and α = 180°.For our study, we treat the PAD near Callisto as isotropic.We do investigate the influence of anisotropies in the PAD on ENA emissions at Europa (see below).However, we expect the effect of such anisotropies to be very similar at both moons.Using an isotropic PAD for all Callisto runs of our ENA generation model, p is unity and Equation 8becomes Galileo EPD spectra obtained farther inward within the Jovian magnetosphere demonstrate approximately uniform fluxes across all pitch angles in Europa's energetic proton environment (Mauk et al., 2004).Notwithstanding this, deviations from isotropy have been detected at Europa for certain proton energies: for instance, Galileo observations in the energy range 80 ≤ E ≤ 220 keV show a depletion in proton flux for pitch angles near 90° (Lagg et al., 2003).Kollmann et al. (2016) identified a persistent minimum in proton flux around α = 70° at energies of 130 keV.Most recently, the Juno spacecraft found evidence of an anisotropic proton PAD near Europa at E ≈ 117 keV: the observed proton flux around α = 90° is approximately a factor of two lower than when α = 0° or α = 180° (Shen et al., 2022).Sarkango et al. (2023) demonstrated the presence of a similar anisotropy in the PAD for the energy range analyzed in our study (E ≤ 100 keV; see their Figure 1a).
For most applications of the ENA model at Europa, we follow Mauk et al. (2004) and treat the PAD of the ambient proton population as isotropic.However, to constrain the influence of this assumption on the ENA emission morphology, we also consider a test case for Europa at the center of the plasma sheet with newly launched parent ions weighted according to a PAD that assigns more flux to the population traveling along the magnetic field than that propagating normal to it.Building upon the findings of Shen et al. (2022) and Sarkango et al. (2023), we assume (across all considered energies) that the energetic proton flux reflects the observations by Juno: particularly, the flux perpendicular to the magnetic field is about 50% smaller than it is for pitch angles of α = 0° or 180°.For B 0 antiparallel to the z axis (as is the case in our setup for Europa at the center of the plasma sheet) we can model this PAD with the form 10.1029/2023JA031931 13 of 45 The factor in the denominator of Equation 8 is given by This yields an initial macroparticle weight upon launch equal to When a newly launched proton macroparticle travels through a cross section dA on the starting grid (depicted by the white grid lines in Figure 1), it represents a particle "passage rate" (in units of 1/s) of where  n is the unit vector normal to the grid element (see also Section 3 of Wulms et al., 2010).The dependence of the value dF on the cell size dA of the starting grid is important to bear in mind when comparing model results obtained with different starting grid resolutions.Consider, for example, two ENA model runs carried out for Europa under identical ambient magnetospheric conditions (i.e., the ambient distribution function f in Equation 13 is the same in both cases).However, the cell size dA 1 of the starting grid in the first scenario is larger than the cell size dA 2 of the second scenario (dA 1 > dA 2 ).Thus, the number of proton macroparticles launched within the larger cell area dA 1 is elevated by a factor of dA 1 /dA 2 in the second scenario compared to the first.This implies that in the second scenario, the amount of differential flux J 2 (in units of  [ cm 2 sr keV s ] −1 ) initiated in the larger cell area dA 1 (covered by multiple cells of area dA 2 ) is also increased, compared to scenario 1, by the same factor of dA 1 /dA 2 : However, the passage rate dF (in units of 1/s) through area dA 1 needs to be equal in both model scenarios to ensure a valid comparison.Merely using Equation 14and   ∝ Ĩ() ∝  (see Equation 7) would yield Thus, the rate dF 2 in scenario 2 would be artificially enhanced compared to the rate dF 1 in scenario 1, solely because of the smaller size of the cells on the starting grid in scenario 2. To account for this artificial increase in rate dF, the modeled rate dF 2 in scenario 2 needs to be rescaled according to Similar to the approach of Tippens et al. ( 2022), we propagate energetic proton macroparticles through uniform fields or the perturbed electromagnetic fields calculated by AIKEF (see the blue trajectory in Figure 1).In the case of perturbed fields, the electromagnetic field vectors at the position of the macroparticle are determined from the eight adjacent nodes of the AIKEF simulation grid through trilinear interpolation.The equation of motion for a proton (charge e) at position r in electric field E and magnetic field B is given by The trajectory of each proton macroparticle is evolved through space with a Runge-Kutta scheme of fourth order.Since the relativistic Lorentz factor for protons even at the highest energy considered (E = 100 keV) does not exceed 0.1% deviation from unity, a relativistic treatment of parent ion dynamics is not necessary.
If the particle's trajectory crosses one of the outer boundaries of the model domain enclosed by the starting grid, the trajectory is allowed to evolve over 2 r g (see Table 2) beyond the boundary before it is removed from the simulation (akin to Regoli et al. (2016) where R is again the radius of the respective moon.If the local magnetic field strength |B| is less than the background field magnitude |B 0 |, then the gyroperiod in the background field is used as a baseline for Δt.In other words, the timestep at either level of resolution is not permitted to grow above what it would be in the uniform magnetospheric background field.
During its evolution after being launched, an energetic proton undergoes acceleration by the electromagnetic fields.To calculate the modified flux that it represents when entering Europa's or Callisto's atmosphere, we need to combine Equation 7 with Liouville's theorem.By equating the distribution function at initialization on the starting grid to its value at the top of the atmosphere, we can derive a relationship between the differential particle fluxes at these locations.Namely, if the macroparticle begins its evolution at the starting grid with velocity  ̇1 and weight J 1 , and reaches the top of the moon's atmosphere with a velocity  ̇2 , the modified weight J 2 becomes A similar rescaling of the flux is done whenever a proton exits and subsequently re-enters Callisto's or Europa's atmosphere, as displayed by the particle trajectory in Figure 1.
Apart from charge exchange interactions inside of a moon's atmosphere that generate ENAs, scattering through wave-particle interactions by, for example, ion cyclotron waves (Nénon et al., 2018) and Whistler-mode waves (Shprits et al., 2018) may cause localized phase space depletions that prevent the conservation of the distribution function.So far, the presence of such waves near Callisto has neither been conclusively established through observations nor modeling results.In contrast to Callisto's environment, the model of Desai et al. (2017) revealed that pick-up of ionospheric ions causes ion cyclotron waves to develop within 3R E of Europa's wakeside hemisphere, consistent with those observed during Galileo flybys E11 and E15 (Volwerk et al., 2001).However, data from the few Galileo flybys of Europa provide only limited information on the locations of regions with strong wave activity.Besides, Shprits et al. (2018) emphasize that these regions are probably highly localized.Therefore, we still consider Liouville's theorem to be valid in good approximation outside of Europa's and Callisto's atmospheres, similar to the particle tracing models of, for example, Addison et al. (2021Addison et al. ( , 2022)), and Liuzzo et al. (2022).
We propagate parent protons from the starting grid to the top of each moon's atmosphere.The atmospheric profiles in the ENA generation model are identical to the atmospheric profiles used to calculate the electromagnetic field configurations with AIKEF at each moon (see Section 2.1).In Figure 1, the proton enters the atmosphere when the blue trajectory intersects the region shaded in orange.At an altitude of one moon radius above the surface, the number density of both Callisto's and Europa's atmosphere in our model has fallen by 5 orders of magnitude relative to the surface density.We therefore do not consider ENA generation beyond an altitude of R above each moon's surface.At altitudes greater than R, the atmospheric profiles in our model contain only 0.002% of the integrated atmospheric number density at Europa, and 0.3% at Callisto.For this reason, we expect the ENA production at altitudes greater than R to generate only a very minor contribution to the emitted ENA flux, concordant with the findings of Tippens et al. (2022) at Titan.
At Callisto, the column density of the CO 2 component is approximately a factor of 25 lower than that of O 2 (see also Liuzzo et al., 2015).Therefore, ENA emissions from the CO 2 constituent are expected to be weak compared to those emitted through charge exchange with the O 2 profile.Further, cross sections for charge exchange between protons and CO 2 molecules have been measured only for energies above 100 keV (Belkić, 2021;Toburen et al., 1968).Ergo, we do not consider ENA generation through charge exchange with CO 2 at Callisto.However, the CO 2 component of Callisto's atmosphere is still present in the AIKEF runs where-due to its large mass-it contributes to flow deceleration and field line draping (Liuzzo et al., 2015).Roth et al. (2017a) ascertained the presence of a dilute atomic hydrogen corona at Callisto.However, these authors estimated a column density of (6- 12) ⋅ 10 11 cm −2 , which is at least three orders of magnitude less than that of the molecular oxygen profile.Therefore, this constituent is included neither in AIKEF nor in our ENA generation model at Callisto.In addition, a recent modeling study suggested the presence of an appreciable molecular hydrogen component in Callisto's atmosphere (Carberry Mogan et al., 2022).This component is not included in the current model setup, and we will elaborate on its potential impact on the ENA emissions in Section 3.5.
Whenever a proton macroparticle enters Callisto's or Europa's atmosphere, it continuously interacts with the neutral gas through charge exchange, analogous to Wulms et al. (2010), Kabanovic et al. (2018), andTippens et al. (2022).During a timestep within the atmosphere, a macroparticle moves along a path element ds with length ds = |ds|.The loss of energetic proton weight dJ within the atmosphere is determined by the attenuation equation (Kabanovic et al., 2018;Tippens et al., 2022;Wulms et al., 2010), where n n (r) is the neutral density evaluated at the midpoint of the spatial step ds.The energy-dependent cross section σ(E) for the ENA-generating interaction  H + + O2 → H + O + 2 (at Europa and Callisto) is given in Table 1 of Lindsay and Stebbings (2005) and displayed in Figure 2 (green).For the H 2 O charge exchange interaction at Europa, that is, H + + H 2 O → H + H 2 O + , the cross section can be found in Table 4 of Wedlund et al. (2019) and is also displayed in Figure 2 (blue).
In addition to protons, an abundance of singly charged atomic sulfur and oxygen (as well as multiply charged ions of the same species) is persistently observed in the ambient plasma at Europa and Callisto (e.g., Clark et al., 2020;Cooper et al., 2001;Mauk et al., 2004;Nénon & André, 2019).These energetic ion species are not included in our ENA generation model for two reasons.First, charge exchange cross sections throughout the energy range 1 keV ≤ E ≤ 100 keV are not available in the literature.For instance, cross sections for charge exchange between singly ionized atomic oxygen and neutral molecular oxygen have been determined experimentally only for energies between 1 and 5 keV (Lindsay & Stebbings, 2005), and at three discrete higher energies (Loand & Tite, 1969).Nénon and André (2019) were able to use observations of Europa's neutral torus to estimate that the ratio of charge exchange cross sections between sulfur ions and atomic oxygen or molecular hydrogen is less than 1.5.Nevertheless, no profiles of the energy dependent cross sections σ(E) were determined, in particular not for the energy range relevant to our study.
Second, the distribution of charge states within the energetic heavy ion populations near Europa and Callisto (i.e., the number ratio between singly and multiply charged sulfur or oxygen ions) is still under debate (e.g., Clark et al., 2016Clark et al., , 2020 and references therein).The measured heavy ion energy spectra (e.g., Paranicas et al., 2009) do not discriminate between singly and multiply charged ions.Since, for example, O 2+ or S 3+ ions are not expected to be neutralized on the small length scales of Europa's or Callisto's local interaction regions (e.g., Nénon & André, 2019), only a certain fraction of the incident oxygen/sulfur ion populations would actually be able to generate ENAs near these moons.Hence, the observed heavy ion energy spectrograms would have to be "downscaled" to consider only the ENA-generating, singly charged portion of the incident ion population (with the ratio of singly and multiply charged ions scarcely constrained).For these two reasons, we defer the analysis of ENA-generation by heavy ions to a future investigation.Including the heavy ion species would mainly enhance the intensity of ENA emissions and their detectability by spacecraft.
While inside the atmosphere, each parent ion produces a single ENA macroparticle during each time step.Equation 20 is used to convert part of the flux J carried by the parent ion macroparticle into ENA flux.During each timestep that an energetic parent ion spends in Europa's or Callisto's atmosphere, the ENA generation model records the position and velocity coordinates of the macroparticle as its weight is attenuated by dJ, as well as the flux Y = |dJ| of the newly generated ENA macroparticle (see also Kabanovic et al., 2018;Tippens et al., 2022;Wulms et al., 2010).Since the considered near-resonant charge exchange interactions are highly forward scattered (Dandouras & Amsif, 1999;Lindsay & Stebbings, 2005), the velocity vector of each newly-emitted ENA macroparticle is identical to the parent proton's velocity the moment the ENA is created.Starting from the position of the parent ion at the moment of charge exchange, the rectilinear trajectory of each ENA is traced until it encounters one of two fates: either it intersects the moon's surface and it is deleted from the simulation (red trajectory in Figure 1), or it leads away from the moon and the ENA would be observable outside the atmosphere (green trajectories in Figure 1).Our study focuses on the latter case; that is, on the detectability of ENA emissions that emanate from Callisto's and Europa's atmospheres.
The detectable ENA flux is processed into a flux map through the surface of a spherical detector (centered at the moon) at altitude h = R enclosing the Callisto's or Europa's entire atmosphere, thereby capturing all ENAs traveling away from the moon.Due to the inherent truncation of observable ENA flux by the finite field of view, an actual (point-like) spacecraft detector captures only a tiny fraction of the emitted ENA population.To preserve the entirety of the physics contained in the ENA pattern emitted from the moon's atmosphere, we employ a detector geometry that encapsulates the entire ENA generation region.The same approach was applied by Tippens et al. (2022) to investigate the global morphology of ENA emissions at Titan.Incorporation of a realistic, pointlike detector geometry is a goal for future investigation, for example, after the specifications of the ENA detector aboard JUICE have been made available.The spherical detector surrounding each moon in our model is binned into longitude and latitude segments 3° × 3° in extension and the total flux associated with all ENAs intersecting each bin is calculated.To evaluate the dependence of the ENA emission morphology on the initial energy of the parent protons, ENA emission maps have been calculated separately for each of the five starting energies studied.
We do not consider subsequent reionization and reneutralization of ENAs, as done, for example, for Titan by Garnier et al. (2008).Including these effects may rather be the subject of a future study.Sputtering of an icy moon's surface by impinging energetic ions can also result in ENA emissions, as was modeled at Ganymede by Pontoni et al. (2022).Our study focuses only on ENA generation from interactions with Callisto's and Europa's atmospheres.Taking into account ENAs emitted from the surface as well is beyond the scope of this study, as it would require a conceptually distinct modeling approach.However, this investigation allows follow-up studies to isolate the contribution of ENAs emitted from their atmosphere, compared to those generated at these moons' surfaces.

Perturbed Electromagnetic Fields at Callisto and Europa
Figures 3 and 4 depict the structure of the perturbed electromagnetic environment in the vicinity of Europa and Callisto for each hybrid model configuration.The left two columns of Figure 3 correspond to Callisto at the center and at maximum distance below the center of the Jovian magnetospheric plasma sheet, respectively.The right two columns follow identical organization for Europa.The first row (Figures 3a-3d) and second row (Figures 3e-3h) show the flow-aligned magnetic field component B x and field strength  || in the y = 0 plane, respectively.The third row (Figures 3i-3l) displays the magnitude of the electric field |E| in the z = 0 plane.All quantities in Figure 3 are normalized to the respective upstream value, that is, |B 0 | for the top two rows and |E 0 | for the bottom row.When either moon is at the center of the Jovian plasma sheet, the y = 0 plane contains the background magnetic field vector B 0 , and the z = 0 plane contains the background convective electric field vector E 0 .When Europa is located at its maximum distance below the Jovian plasma sheet, the background magnetic field is rotated against the y = 0 plane (toward Jupiter) by 28.5°, and the background electric field is correspondingly rotated by the same angle against the z = 0 plane.Analogously, when Callisto is located below the plasma sheet, the background magnetic field is rotated by 72.5° against the y = 0 plane toward Jupiter (see Table 1).For this case, a depiction of plasma and electromagnetic field quantities in planes containing the strongly inclined background fields is given in Figure 1 of Liuzzo et al. (2022).Figure 4 shows results for Europa when an H 2 O atmospheric component is included in addition to molecular oxygen.We include in Figures 4a and 4b  and Addison et al. (2021Addison et al. ( , 2022) ) for Europa, and by Liuzzo et al. (2015), Liuzzo et al. (2016Liuzzo et al. ( , 2017Liuzzo et al. ( , 2018Liuzzo et al. ( , 2019aLiuzzo et al. ( , 2019bLiuzzo et al. ( , 2022) ) for Callisto.Therefore, we provide only a brief discussion of the features relevant to our analysis of ENA emissions.For a more detailed discussion, the reader is referred to any of our preceding publications using the same model.
We begin our analysis with the case of either moon at the center of the Jovian plasma sheet, that is, the two left columns of Figure 3.At both Europa and Callisto, the flow-oriented component B x,0 of the upstream magnetospheric field is set to zero in these model setups.As depicted in Figures 3a and 3c, B x becomes negative (blue) north of each moon, and positive (red) south of it, with perturbations reaching strengths of |B 0 |/2 at Europa and 3|B 0 |/2 at Callisto.Impinging magnetospheric plasma is slowed and deflected around the obstacle due to ionospheric mass-loading (e.g., Liuzzo et al., 2015;Rubin et al., 2015;Saur et al., 1998).The frozen-in magnetic field is draped around the moon, as depicted in the cutting plane containing the background field B 0 (Figures 3a  and 3c).At large distances from either moon, where field contributions from ionospheric currents are negligible, the interaction generates a system of Alfvén wings which propagate along characteristics   ± = 0 ± 0 and connect Jupiter's polar ionosphere to the local interaction region (Neubauer, 1980(Neubauer, , 1998)).At Callisto, where the flow is super-alfvénic when near the center of the plasma sheet, the Alfvén wings are inclined against B 0 by arctan M A ≈ 74°.At Europa, however, the Alfvén wings are only inclined by 34° relative to the background  1).Within each column, the left sub-column corresponds to the respective moon at the center of the Jovian plasma sheet, and the right sub-column corresponds to the moon at its maximum distance below the plasma sheet.To ensure best possible visibility of the field perturbations, the ranges of the color bars are different between panels (a) and (b) as well as panels (e)-(h).magnetic field.An enhancement in magnetic field magnitude around the upstream apex of either moon is visible in Figures 3e and 3g, depicted in yellow.At Callisto, the enhancement in magnetic field magnitude reaches a value larger than 5|B 0 |.However, at Europa, the field strength in the pileup region is intensified only by a factor of 1.75 (see Figure 3g).When Callisto is located at the center of the Jovian plasma sheet, the background magnetospheric field strength is over a factor of 100 weaker than at Europa (see Table 1).The plasma beta is also three orders of magnitude higher than at Europa.Due to the reduced "stiffness" of the magnetospheric field lines at Callisto, the field strength in the pileup region is elevated (relative to the respective background field) much more than at Europa.
The region between the Alfvén wings where B x becomes zero corresponds to a thin neutral sheet in Callisto's wake (Liuzzo et al., 2022).In this region, the magnetic field magnitude falls to approximately |B 0 |/2 (see Figure 3e).In contrast to Callisto where the Alfvén wings are strongly inclined toward the corotational flow direction, the Alfvén characteristics at Europa are only slightly inclined relative to B 0 .Therefore, in the plane containing B 0 , the region of reduced field strength is widened and exhibits a triangular shape.In Figure 3g, this region is depicted in dark purple with a wakeside depletion of about 30% of |B 0 |.The downstream deficit of magnetic pressure is compensated by the enhanced particle pressure of the ionospheric pick-up ions that populate this region at each moon (see, e.g., Liuzzo et al. (2015) or Arnold, Liuzzo, and Simon (2020) for details).
In Figures 3i and 3k, within the z = 0 plane (i.e., parallel to E 0 ) the electric field magnitude falls close to zero (dark purple) in the wake region of either moon (x > 0).At Europa, this depletion region is approximately symmetric between the Jupiter-facing (y > 0) and Jupiter-averted (y < 0) half spaces.Comparatively, the region where E is reduced near Callisto extends about one R C further into the Jupiter-averted half space, and |E| returns to its background value in the y > 0 half space at a much larger x value than it does in the y < 0 half space.At both moons, the wakeside depletion in |E| is a consequence of reduced bulk flow velocity u as the magnetospheric plasma is mass-loaded by slow ionospheric constituents.In Europa's environment, where the Jovian background field strength is |B 0 | ≈ 400 nT, pickup ion gyroradii are less than 0.2R E .Therefore, the pickup tail is approximately symmetric about the y axis in planes parallel to the undisturbed convective electric field (see also Arnold, Liuzzo, & Simon, 2020).At Callisto, the magnetospheric field strength is two orders of magnitude lower than at Europa, so heavy O 2 and CO 2 pickup ions have gyroradii larger than 6R C , yielding a highly asymmetric pickup tail (see also Liuzzo et al., 2015).
In the following, we discuss each moon's interaction at its maximum distance below the Jovian plasma sheet, that is, the two right side sub-columns of Figure 3.When Callisto is located at this position, the background magnetospheric field is rotated by 72.5° against the y = 0 plane, and therefore the field line draping is only weakly visible in the cut shown by Figure 3b.Alternatively, at Europa the magnetospheric background field, and hence, the Alfvén wing characteristics  (  ± = 0 ± 0 ) are rotated out of the y = 0 plane by only 28.5°.In Figure 3d  Immediately north and south of Europa, the B x component (Figure 3d) is negative (blue) and positive (red), respectively.At larger distances, the B x component in this plane reverses its sign (B x > 0, red, in the north, and B x < 0, blue, in the south).In planes (approximately) perpendicular to the wing characteristics, the field lines of an Alfvén wing take the shape of a two-dimensional dipole (Neubauer, 1980).Consequently, outside of the wing tubes (which connect to Europa's ionosphere) the field lines form regions of "anti-draping," that is, the field vectors are antiparallel to those within the wing tubes to ensure field line closure (see, e.g., Figure 8 of Simon et al., 2011).This effect is visible in Figure 3d, since the distance between the center of the wing tubes and the cutting plane grows with increasing |z| value.At Europa, the field strength in the ramside pileup region (Figure 3h) reaches 1.3|B 0 |, accompanied by a depletion in magnitude downstream of the moon.
In the immediate vicinity of Callisto, the induced dipole acts to locally amplify the magnitude of the magnetic field by about 40%, as depicted by the bright yellow ring within the y = 0 plane (see Figure 3f).Since the induced dipole moment mainly points in (−y) direction, the vectors of the induced field possess a dominant component in the (+y) direction (parallel to B 0 ) within the y = 0 plane, thereby causing an approximately radially symmetric enhancement in the magnetic field magnitude within this plane.At Callisto, the electric field magnitude is reduced by up to 60% in a broad region extending downstream as well as into the Jupiter-facing and Jupiter-averted hemispheres, as shown in Figure 3j.This reduction in field strength is caused by the deflection of the upstream plasma flow around the Alfvén wings, whose characteristics   ± are close to the z = 0 plane.This deflection reduces |u| and produces a corresponding reduction in the magnitude of the convective electric field.The region of reduced electric field strength is asymmetric between y > 0 and y < 0 due to the positive B x,0 component of the magnetospheric background field (see also Simon & Motschmann, 2009).In addition, when Callisto is at its maximum distance below the center of the Jovian plasma sheet, the 72.5° tilt of |B 0 | against the y = 0 plane rotates the pickup tail away from the z = 0 plane, and therefore the tail does not contribute appreciably to the large-scale reduction of |E| in Figure 3j.However, at Europa, the tail of escaping ions is still primarily contained to the z = 0 plane, as can be seen from the similarity in the structure of |E| between Figures 3k and 3l.Even below the center of the Jovian plasma sheet, the background magnetic field at Europa is still nearly aligned 10.1029/2023JA031931 20 of 45 with the z axis (see Table 1), and as a result the electric field depletion remains approximately symmetric between the y > 0 and y < 0 half spaces (see Figure 3l).
For the case of an additional water atmosphere included around Europa's subsolar and (coincident) upstream apex (Roth, 2021), Figures 4a and 4b display the number density of the H 2 O + pickup tail in the planes containing the center of Europa, u 0 , and E 0 or B 0 , respectively.In the y = 0 plane (Figure 4b), the H 2 O + tail exhibits a gap in density between Europa's ionosphere (dark red) and the distant tail region (starting around x ≈ 8R E ).The reason for this is revealed in Figure 4a: the H 2 O + ions forming the tail are swept around the Jupiter-averted and Jupiter-oriented faces of Europa as they propagate downstream.Immediately downstream of Europa, the tail takes a two-pronged shape because the H 2 O + ions (which are generated only around the ramside apex) must first be transported around the moon in order to reach the downstream region.These ions are thus carried out of the y = 0 plane around both sides of Europa, only for the separated tails to rejoin about 8R E downstream (Figure 4b, shown as light blue and green).The ionospheric H 2 O + has a factor of 10 greater number density than the ambient magnetospheric plasma (n 0 = 200 cm −3 ) in the region immediately surrounding the moon, and its density falls to values similar to n 0 at approximately 7R E downstream.The structure of the pick-up tail formed by  O + 2 ions from Europa's global atmosphere looks somewhat different.However, the physics involved in shaping the  O + 2 tail have already been discussed in detail by, for example, Arnold, Simon, and Liuzzo (2020), and we refer the reader to that publication for further information.
In Figure 4a, a wave-like pattern is visible in the long, Jovian-averted flank of the H 2 O + pickup tail.In the y < 0 half space (where E 0 points away from Europa), flow shear between ionospheric H 2 O + and magnetospheric thermal plasma results in the excitation of bi-ion waves.The location and shape of the wave signature are consistent with modeled bi-ion wave features in the pickup tails of comets (Bagdonat & Motschmann, 2002) or Pluto (Delamere, 2009;Feyerabend et al., 2017).Specifically, in agreement with our results for Europa, such waves were found to form only in the half space where E 0 points away from the respective object (e.g., Delamere, 2009;Feyerabend et al., 2017).The electric field perturbations in the z = 0 plane (Figure 4e) are similar to the case of Europa at the center of the Jovian plasma sheet without the H 2 O atmosphere (Figure 3k), except for a filamented structure visible along the y < 0 boundary of the depletion region.This structure corresponds to the aforementioned wave signature which maps into the bulk velocity, and hence, leaves a weak imprint on the electric field.
In the y = 0 plane, the magnetic draping pattern remains very similar to that obtained without including the H 2 O atmospheric component (Figures 4c and 4d).This finding is consistent with the MHD modeling results of Cervantes and Saur (2022), who demonstrated that, even when the ram-wake asymmetry of the global O 2 atmosphere is not considered, the water vapor component generates only localized modifications to Europa's magnetic environment.

Variability of ENA Emissions at Europa
Figure 5 displays the global ENA emissions captured by a spherical detector around Europa when the moon is located at the center of the Jovian plasma sheet.The model results for a (hypothetical) case when the magnetospheric field is uniform and southward are given in the left column, and results for motion of the parent protons through draped fields are given in the right column.The ENA flux onto the spherical model detector is projected onto a rectangular surface for visualization, using the West Longitude coordinate system described in Section 2. Each row corresponds to a single initial energy for the parent protons when launched from the starting grid: 1 keV in row (a-b), 10 keV in row (c-d), 40 keV in row (e-f), 70 keV in row (g-h), and 100 keV in row (i-j).Due to acceleration of the parent protons in the draped electromagnetic fields, the energies of the detectable ENAs may differ from the initial proton energies.However, for energetic protons traveling through Europa's environment, the changes in energy are generally small (see Figure 13 of Nordheim et al., 2022).
In the case of the uniform, southward field (Figures 5a, 5c, 5e, 5g, and 5i) a dominant band feature of high ENA flux parallel to Europa's equator is visible across the center of each panel.This band is shaped like an annulus on the detector, wrapping around Europa and oriented normal to the background magnetic field B 0 .For example, 83% of the total ENA flux through the detector in the 40 keV run (Figure 5e) is contained between latitudes of ±33°.Comparing the emission patterns in uniform fields at different energies (Figures 5a, 5c, 5g, and 5i) reveals that the morphology of this equatorial band exhibits little variation as a function of parent ion energy.Within the band, the detected ENA flux is maximized around 300°W longitude at each energy, that is, slightly west of the detector's ramside apex (270°W longitude).However, the magnitude of the peak flux varies non-monotonically as a function of initial proton energy.The PAD of the energetic parent ions is isotropic in this model setup.However, ENA emissions at Europa are concentrated within this equatorial band perpendicular to B 0 because parent ions with pitch angles near 90° make the dominant contribution to the ENA flux.To illustrate this, we consider two test particles initialized on the starting grid at an energy of E = 40 keV: the first of these parent ions is launched from the upstream face of the starting grid and has a pitch angle of 90°.Its guiding center drifts along a straight line parallel to the (+x) axis, confined to the y = 0 plane and displaced northward by z = R E + h s,E , where h s,E is the atmospheric scale height (see Equation 4).Thus, all ENAs emitted by this parent ion hit the detector; none of them are obstructed by Europa's solid body.The "ENA trace" left by this proton on the detector is confined to a latitude of  arcsin ( , that is, it forms a circle.The second parent ion has a pitch angle of 0° and is initialized on the starting grid's northern (+z) face such that its guiding center moves southward within the y = 0 plane and downstream of Europa at x = R E + h s,E .In other words, the trajectory of the guiding center is perpendicular to that of the first particle.This second test particle translates mainly along B 0 because the corotation drift velocity in (+x) direction (|u 0 | = 100 km/s) is slow relative to the speed = 2, 768 km∕s of a parent proton at E = 40 keV.The length of the guiding center's trajectory through the atmosphere is the same for both test particles:  = 2 √ (2) 2 − ( + ℎ) 2 ≈ 3.39 .However, these two protons' capacity to generate ENA flux is dependent on their actual path length within the atmosphere, and not merely the path length l traveled by the guiding center.Gyration causes the actual path within the atmosphere to be extended compared to l only for the first proton.Taking into account its drift in (+x) direction at velocity u 0 , its guiding center traverses Europa's atmosphere in l/|u 0 | ≈ 53 s, equal to the time it takes to complete about 330 gyrations within the atmosphere.Thus, the parent proton's path length inside the atmosphere associated with gyration alone reads 330(2πr g ) = 93.6RE (where r g = 0.045R E is the gyroradius, see Table 2).In contrast, the second parent ion travels along the magnetic field and traverses a path length of l ≈ 3.39R E inside of the atmosphere (in less than 2 s).
Therefore, the first parent ion interacts with the atmosphere for at least 28 times the path length that the second ion does.In other words, parent ions drifting with pitch angles near 90° trace the longest paths through the atmosphere and emit the majority of the ENA flux through the detector, forming the band of high flux perpendicular to B 0 seen in the left column of Figure 5.In addition, the second particle (moving along the magnetic field) emits ENAs into only a small localized patch on the detector, clustered at the point where its trajectory exits the detector sphere.Taking into account parent protons from the entire (isotropic) PAD, our example shows that their path length within the atmosphere maximizes at α = 90° and falls rapidly as the velocity gains a field-aligned component.The intense ENA emissions detectable at Europa stem from long path lengths within the atmosphere, which are a result of the energetic proton's gyroradii being small compared to the extension of the moon's neutral gas envelope.This still holds at the highest energies considered because even a 100 keV proton has a gyroradius of only r g = 0.07R E ≈ 1.1h s,E (see Table 2).
In the purely southward magnetic field used in this model configuration, ions with pitch angles of 90° gyrate in planes of constant z, which is why the outer boundaries of the band are circles with constant z value (i.e., constant latitude) on the detector sphere.In other words, the orientation of the equatorial band of high ENA flux on the detector sphere (left column of Figure 5) is governed by the orientation of the background magnetic field B 0 .The latitudinal extension of the band is determined by the extension of the Europa obstacle along the magnetic field lines (i.e., in ± z direction), which is approximately  = ±( + ℎ) .At these z values, the gyroplanes of energetic protons are approximately tangential to Europa's dense, lower atmosphere (see also Figure 6a in Tippens et al., 2022).ENAs emitted at  = ±( + ℎ) intersect the detector at a latitude θ b given by where R D = 2R E is the radius of the spherical detector (see also Equation 11in Tippens et al., 2022).Substituting h s,E = 100 km yields θ b = ±32.14°,consistent with our modeled results (left column of Figure 5).Europa is located at the center of the Jovian plasma sheet (see Table 1).The West Longitude coordinate system is used in this projection.Each row represents ENA emissions generated by parent ions initialized on the starting grid at a certain kinetic energy, namely 1 keV in row (a-b), 10 keV in row (c-d), 40 keV in row (e-f), 70 keV in row (g-h), and 100 keV in row (i-j).The left column (panels (a), (c), (e), (g), and (i)) illustrates the ENA emission morphology when the energetic parent protons are evolved through the uniform Jovian background fields B 0 and E 0 .The right column (panels (b), (d), (f), (h), and (j)) displays ENA emissions when the parent ions move through the perturbed electromagnetic fields from the adaptive, kinetic ions, electron fluid hybrid model (see panels (c), (g), and (k) of Figure 3).
The band of enhanced ENA flux does not gradually disappear with increasing latitude, but its boundaries are sharply defined around the "cut-off latitude" θ b approximately given by Equation 21.For instance, at an initial parent proton energy of 70 keV, the detected ENA flux transitions from about  1500 [ cm 2 sr keV s ] −1 (turquoise) to  250 [ cm 2 sr keV s ] −1 (dark blue) across less than 5° in latitude at the edges of the band (see Figure 5g).This rapid drop in detected ENA flux with increasing latitude is mainly caused by the exponential decrease of atmospheric density with rising altitude (see Equation 4).Since the gyration planes of parent ions that generate the band are nearly parallel to the z = 0 plane, parent ions with increasing |z| positions interact with a successively less dense region of the atmosphere, which in turn generates reduced ENA emissions.In addition, the spherical shape of the atmosphere contributes to the steep decrease of detected ENA flux with latitude.For ions with pitch near 90° and |z| > 1R E , (i.e., their gyroplanes intersect Europa's atmosphere but not Europa) the path length of the guiding center within the atmosphere decreases with increasing |z| value.This reduces the length of the trajectory segment where such an ion can generate ENAs that would populate the detector at a specific latitude, further amplifying the latitudinal decrease in ENA flux with increasing |z|.These two effects combined result in the stark decrease of detected ENA flux with latitude, as seen in the left column of Figure 5.
In uniform fields, the ENA flux within the equatorial band is maximized at longitudes slightly west of the ramside apex, shifted approximately 20°-30° in longitude toward the Jupiter-facing apex.The observable ENA flux intensity drops with increasing longitudinal distance to the maximum, with the minimum also positioned slightly west of the detector's wakeside apex (see, e.g., Figure 5g).This asymmetry is most clearly visible at an initial parent ion energy of 40 keV in Figure 5e, but it is a persistent trend also apparent in Figures 5c, 5g, and 5i.To investigate the reason for the shift in the longitude of maximal emissions, we conduct a numerical experiment: consider a single, proton trajectory that is initialized with an energy of 40 keV.This proton has a pitch angle of 90° with an initial velocity parallel to  ŷ .It is launched on the starting grid at r(t = 0) = (−8R E , 0, 0), such that its trajectory is given by where Ω denotes the gyrofrequency.Therefore, the proton's trajectory is confined to the z = 0 plane, and the same holds for the trajectories of all ENAs it emits.Most importantly, the proton spends equal path lengths of its trajectory in the y > 0 and the y < 0 half spaces.We remind the reader that in our model, Europa's atmosphere is symmetric between these two half spaces.
The evolution of this parent ion is illustrated in Figure 6.Figures 6a and 6b show projections of its trajectory onto the z = 0 and y = 0 planes, respectively.Figure 6c displays the relative contributions of the detectable ENA flux onto the Jupiter-facing (or Jupiter-averted) hemisphere of the detector, corresponding to cyan (or gold).Figures 6d and 6e then show the atmospheric density n n (r) at the proton's position and the ENA flux generated at each point in time, respectively.The shading in Figures 6d and 6e denotes the sign of   , which reverses each half-gyration.Since each ENA is emitted with the parent proton's momentary velocity, the sign of   determines whether the ENA is emitted toward the y > 0 or the y < 0 half sphere of the detector (see Figures 6d and 6e).The southward orientation of the magnetospheric background field B 0 causes the parent proton to gyrate in the mathematically positive sense.Therefore, as it drifts into the atmosphere and begins to generate ENAs, it is located at lower altitudes (i.e., in a more dense region of the atmosphere) when it emits ENAs toward the Jovian-facing hemisphere of the detector  ( ̇ > 0) , and is located at higher altitudes (i.e., in a less dense region of the atmosphere) when it emits ENAs toward the detector's Jovian-averted hemisphere (i.e.,   < 0 ).As shown in Figures 6d  and 6e, maxima in local atmospheric density at the proton's location coincide with the center of each interval for which   > 0 , and also with maxima of the emitted ENA flux.
This asymmetric emission process causes a larger amount of ENA flux to be detected on the Jovian-facing hemisphere of the detector compared to the Jovian-averted hemisphere (see Figure 6c).In other words, the slight shift of the ENA emission maximum toward Jupiter (see left column of Figure 5) arises from the sense of proton gyration in the southward magnetic field B 0 , combined with the fact that these particles can complete numerous gyrations within Europa's atmosphere (see also Figures 6a and 6b).We note that analogously, the minimum of ENA flux does not occur exactly at the wakeside apex of the detector (90°W longitude), but is shifted toward larger longitudes (see left column of Figure 5).The solid body of Europa absorbs part of the emitted ENA population.Consequently, not all emitted ENAs are ultimately detected in our experiment displayed by Figure 6.However, during each gyroperiod, the drift displacement of an energetic proton toward Europa is smaller than the gyroradius and significantly smaller than R E (see Table 2).Specifically, on scales of a single gyration the proton trajectory is approximately circular.Therefore, absorption of newly generated ENAs at Europa's surface does not introduce any new, substantial asymmetries into the ENA flux through the detector sphere.
When Europa is at the center of the Jovian plasma sheet, the average ENA flux detected in draped electromagnetic fields is diminished by about 50% in intensity across all energies compared to the average detectable ENA flux in uniform fields (see right column of Figure 5).The slight longitudinal shift of the emission maximum toward Jupiter persists when field line draping is considered.The equatorial band of high ENA flux that was visible at all energies in uniform fields (see left column of Figure 5) is no longer present in draped fields.For instance, the region of (slightly) enhanced ENA flux on the wakeside of the detector shrinks in latitudinal extent when moving from uniform to draped fields.In the latter setup, this region is confined between latitudes of ±10° (see, e.g., Figures 5d-5f), compared to ±33° in uniform fields.However, a "bulge-like" enhancement in ENA flux remains slightly west of the detector's ramside apex, reaching, for example,  1750 [ cm 2 sr keV s ] −1 at an initial proton energy of E = 40 keV (see Figure 5f).This corresponds to 44% of the maximum value reached in the uniform case at the same energy (  4000 , see Figure 5f).
An example of an energetic proton trajectory that contributes to the modified ENA flux pattern in fields is depicted in Figure 7. Two energetic protons are propagated with identical initial conditions: one through uniform electromagnetic fields (left column) and another through perturbed fields from AIKEF (right column).In uniform fields, the parent proton can readily penetrate into Europa's wakeside atmosphere, emitting ENA flux onto the detector (green segment of the trajectory).These ENAs contribute to the signature on the wakeside of the detector, visible in the left column of Figure 5.In draped fields, a parent proton launched with the same initial conditions is already deflected while north of Europa, far before it ever reaches the moon's atmosphere (right column of Figure 7).The draped field lines cause the energetic proton's trajectory to descend toward the ramside atmosphere (see blue arrow in Figure 7d).ENAs emitted along this parent proton's trajectory at |z| > R E may populate any longitude of the detector.However, the proton's ability to populate the wakeside of the detector at low latitudes is highly limited, since Europa absorbs most of the ENAs emitted into the x > 0 half space.Thus, in draped fields the proton can populate the wakeside of the detector with ENAs only from high latitudes along its trajectory, that is, from regions where the atmosphere is dilute (and hence, the emitted ENA flux is low).In contrast, when the proton travels through the dense atmosphere near the ramside apex, all detectable ENAs are emitted toward the upstream hemisphere of the detector.This parent proton's trajectory in both field configurations exemplifies how draped fields can redistribute parent protons' ENA emissions from the wakeside of the detector to the ramside, thereby contributing to the thinning of the wakeside band seen in Figure 5.
An example parent proton that illustrates a separate mechanism that acts to attenuate the detectable ENA flux in draped fields is displayed in Figure 8.This 40 keV proton is initialized at that is, upstream of Europa, in both uniform fields (left column) and draped fields (right column).The proton is launched with an initial velocity of x , corresponding to a pitch angle of 90° in the undisturbed upstream field B 0 .Since the field is oriented southward, this ion's gyration in uniform fields remains confined to the z = −R E /2 plane throughout its evolution until impacting Europa.Thus, the detected ENA flux associated with this particle is confined to a single latitude on the detector of −14°.Due to absorption by Europa's solid body, this proton cannot emit any ENAs toward the downstream hemisphere of the detector.Such parent ions provide the dominant contribution to the ramside portion of the band feature seen in the left column of Figure 5.Alternatively, in draped fields, the same parent proton is not capable of generating any ENA flux at all, but instead entirely "misses" Europa's atmosphere.Upon encountering the outer edge of the moons ramside magnetic pileup region near x ≈ −4R E (see Figures 3g and 8d), the parent ion's pitch angle is shifted away from 90° by the mirror force in the enhanced B field.Then, the proton rapidly propagates into the z < 0 half space and passes the moon at large southern Jovian latitudes (i.e., outside of the plotted domain).Similar deflection of many such ions thus contributes to the reduction of the ENA flux intensity in draped fields.
For both uniform and draped fields, the morphology of the detectable ENA emission patterns (see Figure 5) varies only slightly across different energies.However, the ENA flux intensity depends non-monotonically on the parent protons' initial energy.Within the range studied, the intensity of the detected ENA flux is smallest at an initial proton energy of 1 keV (row 5(a-b)), peaks at 40 keV (row 5(e-f)), and subsequently falls off with increasing energy.At initial parent proton energies of 1 keV (see Figures 5a and 5b), the maximum detected ENA flux is less than 5% of the maximum value detected for parent protons initialized at 40 keV.At initial energies of 100 keV, the maximum ENA flux detected has again fallen and is now about 60% of its value for 40 keV parent protons (see Figures 5i and 5j). .We include two separate proton trajectories with the same initial conditions, traced through uniform (left column) and draped fields (right column).The arrow in each panel indicates the direction of guiding center motion during the ion's initial approach of Europa.When outside of the atmosphere, the proton trajectory is plotted in cyan (uniform fields) or magenta (draped fields).Both cases are colored in green when they are inside the atmosphere and generate ENAs.For both uniform and draped fields, the sample ion is initialized with The intensity  Ĩ() of energetic proton flux measured outside of Europa's interaction region first increases with energy, reaches a maximum around 60 keV, and then falls off again at higher energies (see Section 2.2 and Paranicas et al. (2002)).The value of  Ĩ() is more than a factor of 10 greater at the upper edge of our energy regime (100 keV) than at its lower edge (1 keV).The cross section σ(E) for charge exchange between protons and molecular oxygen decreases monotonically as a function of increasing parent ion energy (see Figure 2).According to Equation 20, the intensity of ENA emissions generated by a parent ion of energy E is dependent on the product of  Ĩ() and σ(E).Therefore, the non-monotonic behavior of  Ĩ() -combined with the steady decrease of σ(E)-maps into the ENA emission patterns, causing the flux through the detector to peak in our 40 keV model scenario.
We now discuss the global morphology of ENA emissions at Europa when the moon is located at its maximum distance below the center of the Jovian plasma sheet.In this setup, the ambient magnetospheric field possesses a nonzero horizontal component B h .Therefore, the induced field from Europa's subsurface ocean (see Equation 1) now influences parent ion dynamics in the moon's vicinity.To isolate the impact of the induced field on the ENA emissions, we calculate emission maps (at all five initial energies) for an additional case where the ambient magnetic field is obtained by superimposing the uniform magnetospheric background field B 0 and the induced dipole field, but without including any plasma effects.In Figure 9, we present the ENA fluxes calculated for each of the three model field configurations, that is, the uniform background field alone (left column), the uniform background field superimposed with Europa's induced dipole (center column), and the perturbed electromagnetic fields resulting from the plasma interaction with Europa's ionosphere and induced dipole (right column).
In contrast to ENA emissions at the center of the plasma sheet, in Figure 15 the band of high ENA flux is clearly discernible for all initial proton energies in draped fields.Although the width of the high flux band remains the same as it was at the center of the plasma sheet (left column of Figure 5), the band now shows a sinusoidal variation as a function of longitude.This apparent "waviness" of the band is generated solely by the chosen projection in the plots: the band is still circular and present across all longitudes, centered in latitude around the great circle perpendicular to the magnetospheric background field B 0 .In contrast to the scenario from Figure 5, the background field is no longer aligned with the (−z) axis but possesses a strong component toward Jupiter (see Table 1).The high-flux band is mainly generated by parent ions with pitch angles near 90°.Hence, the tilt of the background field against the north-south direction inherently rotates the gyroplanes of these parent ions.This results in an analogous rotation of the high flux band on the detector sphere (see also Figure 6 of Tippens et al., 2022).The maximum latitudinal deviation of the band's center from 0° is thus approximately equal to the tilt of the background magnetospheric field against the z axis: 28.5° (particularly visible in Figures 9d, 9g, and 9j).
At all five initial proton energies studied, the intensity of ENA emissions at maximum distance below the Jovian plasma sheet is persistently greater than in the corresponding field configuration (uniform or draped) at the center of the sheet.In Figures 9d, 9g, and 9j, ENA flux detected within the band is at or above  3500 [ cm 2 sr keV s ] −1 across the entire ramside hemisphere of the detector.In the analogous case at the center of the plasma sheet, the ENA flux only reaches this value within a tightly confined region of the band.Besides, such high values occur only at an initial proton energy of 40 keV (see left column of Figure 5), that is, near the peak of the ambient proton flux  Ĩ() .For instance, Figure 9g demonstrates that in uniform fields, a flux of over  3500 [ cm 2 sr keV s ] −1 is detected across a longitudinal width of 185° (from 200°W to 25°W) inside the band, compared to a longitudinal width of only 85° at the center of the plasma sheet (see Figure 5e).The intensity  Ĩ() of the energetic proton flux observed outside of Europa's interaction region is elevated by about a factor of 3 at initial energies E = 1 keV, E = 10 keV, and E = 40 keV when moving away from the center of Jupiter's plasma sheet (Paranicas et al., 2002, see their Figure 2a).This drives the enhancement of ENA flux at the lower half of our energy range, compared to the center of the plasma sheet (Figure 5 vs. Figure 9).
As can be seen from the middle column of Figure 9, the induced dipole field (in isolation) plays a minor role in influencing the intensity and shape of ENA emissions at Europa.The rather subtle effect of the dipole is revealed by comparing, for example, Figures 9g and 9h: a portion of the high-flux region in Figure 9g is "cut away" by the dipole (see Figure 9h), and the ENA flux within the wakeside band drops by about 15%.Besides, comparing Figures 9j and 9k illustrates that in uniform fields, the region of maximum flux (about  3500 ) is confined to the center of the band, whereas this value is only reached within a thin segment near the band's edges when the induced dipole field is included.
In the case of draped fields (right column of Figure 9), the intensity of ENA emissions everywhere on the detector is reduced by approximately 40% at initial proton energies from 10 to 70 keV compared to uniform fields, and by about 60% at 100 keV.Below the center of the plasma sheet, the band of high flux still encloses the entire detector sphere when draping is included.When Europa is at the center of the Jovian plasma sheet, the magnitude of the magnetic field reaches 1.7|B 0 | in the ramside pileup region, compared to only 1.3|B 0 | far outside the sheet (see Figures 3g and 3h).This reduced magnetic field pileup upstream of Europa is less effective in diverting impinging parent ions away from the moon, thereby having slightly less effect on the ENA emission morphology than at the center of the sheet.
Regardless of Europa's distance to the center of the Jovian plasma sheet, a spacecraft on a polar flyby is likely to observe a much weaker ENA flux from the atmosphere than one viewing the interaction from slightly west of the ramside pileup region within the moon's equatorial plane.Field line draping-which renders additional "protection" to Europa's polar atmosphere-is substantially weaker below the center of the Jovian plasma sheet compared to within it (see Table 1): perturbations in B x are only about half the magnitude of those in the center case (see Figures 3c and 3d).This increases the accessibility of energetic protons to the moon's atmosphere (see also Addison et al. (2021)), contributing to the overall enhancement in ENA flux through the detector.The empirical model for the thermal magnetospheric density utilized in this study predicts a density value of n 0 = 200 cm −3 at the center of the Jovian plasma sheet, and a density of n 0 = 40 cm −3 at maximum distance below (Roth et al., 2014a).However, Bagenal and Delamere (2011) suggest a much less oscillatory upstream density at Europa, ranging from 50 cm −3 at the center to 37 cm −3 at Europa's maximum distance from the center.Both values provided by Bagenal and Delamere (2011) are similar to the value used in our study for Europa at its maximum distance below the center of the Jovian plasma sheet, where the draping generates only minor modifications to the ENA flux pattern (see Figure 9).In other words, when using the upstream density values from Bagenal and Delamere (2011) the features in the ENA flux morphology at the center of the plasma sheet would be more similar to those from below the sheet than suggested by our results (see Figure 9).

ENA Emissions at Europa With an Anisotropic Pitch Angle Distribution
Figure 10 displays the global ENA emissions through the detector sphere when Europa is located at the center of the Jovian plasma sheet, with two different PADs used to initialize the parent protons on the starting grid.In both cases, the ENA flux maps have been calculated using the draped fields from AIKEF, and the model results for an isotropic PAD (see Figure 5 With the anisotropic PAD, the qualitative morphology of ENA flux through the detector is highly similar to the case of an isotropic PAD.In both instances, an enhancement in ENA flux is formed near the detector's ramside apex, and a narrow band of slightly elevated ENA flux along the equator can be seen at the wakeside (left and center columns of Figure 10).Since B 0 is aligned with the (−z) axis, both features are approximately symmetric about the detector's equator.Within the emission enhancement formed around the ramside apex, the ENA flux detected with the anisotropic PAD is elevated by no more than 7% compared to the isotropic case (see right column of Figure 10).At all initial proton energies above 1 keV, the anisotropic PAD slightly elevates the ENA flux at mid-latitudes, that is, directly north and south of the central band.
The initial velocity and position coordinates for the energetic protons launched on the starting grid are the same for both PADs; only the initial numerical weights J of the parent protons differ between the isotropic and anisotropic PAD population.Therefore, the inclusion of an anisotropic PAD can only alter the intensities of the ENA flux through certain regions of the detector, but it does not change the overall pattern of the detected flux since the atmospheric accessibility is not altered by changes to the PAD.The anisotropic PAD assigns a higher differential flux J to parent protons with pitch angles near 0° and 180° than to protons with α ≈ 90°.Using Equations 8-11 reveals that the parent protons initialized with a pitch angle α between  arcsin  4 ≈ 52 • and  180 • − arcsin  4 ≈ 128 • carry less flux with the anisotropic PAD than the same protons initiated with an isotropic PAD.Thus, for the ENA flux to be lower in a certain region of the detector with the anisotropic PAD (i.e., red regions in the right column of Figure 10), the ENA flux in that region must mainly stem from parent ions with initial pitch angles in the range 52° ≤ α ≤ 128°.Conversely, the blue areas in the right column of Figure 10 imply locations where the detector receives substantial ENA flux from parent ions with initial pitch angles α ≤ 52° or α ≥ 128°.
For an anisotropic PAD, the detected ENA flux near the equator and away from the ramside bulge is reduced in intensity by about  200 [ cm 2 sr keV s ] −1 (or 12%-18%, red region in Figures 10f, 10i, 10l, and 10o) relative to the isotropic case at initial proton energies E ≥ 10 keV.The ENA flux within this region is mainly generated by parent protons with pitch angles near 90° (see Section 3.2).This even holds at the detector's wakeside: protons initialized on the upstream face of the starting grid outside of Europa's disk (i.e.,  √  2 +  2 >  ) with α ≈ 90° emit ENAs toward the wakeside hemisphere wherever the tangent to their helical trajectory intersects the detector sphere near its x = 2R E apex.The anisotropic PAD contains fewer particles with pitch angles around 90°, thereby causing the reduction in detected ENA flux (red band in the right column of Figure 10).
At energies E ≥ 10 keV, the parent proton population initialized with an anisotropic PAD emits a 20%-25% greater ENA flux toward mid-latitude regions on the detector, located between ±33° and ±80° (blue in Figures 10f, 10i, 10l, and 10o).The ENA fluxes in both scenarios again become similar when approaching the detector's poles.For example, at an initial energy of 40 keV, parent protons from an anisotropic PAD emit around  50 [ cm 2 sr keV s ] −1 (or 23%) more ENA flux into the mid-latitude regions of the detector than protons with an isotropic PAD at the same energy (see Figure 10i).This is most visible in Figures 10i and 10l on the detector's upstream hemisphere.ENA flux through the mid-latitude regions on the detector is mainly produced by parent ions that are initialized on the northern and southern faces of the starting grid with pitch angles around 0° (north) or 180° (south).These protons propagate with a field-aligned velocity that is large compared to both the gyration velocity and drift velocity.ENAs generated by these parent protons are mainly emitted along the magnetic field lines and can be detected at all latitudes (see, e.g., Figures 8b and 8d).Since the flux carried by parent protons with pitch angles outside the range 52° ≤ α ≤ 180° is higher for an anisotropic than for an isotropic PAD, the ENA flux through the mid-latitude regions of the detector is slightly elevated in the anisotropic case.Invoking the same mechanism, ENAs that travel toward the detector's polar caps are mainly generated by parent protons moving along magnetic field lines within Europa's disk  √  2 +  2 ≤  .However, such ENAs are intercepted by the moon's surface before they can reach the detector sphere and therefore, cannot contribute to the observable flux.For this reason, the flux through the detector's polar caps is equally low, regardless of the chosen PAD (white in Figure 10).
When switching from an isotropic PAD to an anisotropic PAD, the difference in ENA emission intensity at specific positions on the detector ranges from less than 0.1% to at most 35% at all initial proton energies above E = 1 keV.Globally, for parent protons initialized with an energy of, for example, E = 40 keV, the integrated flux through the detector reads  3.52 ⋅ 10 6 [ cm 2 keV s ] −1 with an isotropic PAD and  3.60 ⋅ 10 6 [ cm 2 keV s ] −1 with the anisotropic PAD.Hence, including the anisotropic PAD only slightly increases the total signal at this energy.The region where the relative change peaks on the detector is located at the mid-latitudes (±66°), above and below the detector's ramside apex and colored in dark blue (Figures 10f,10i,10l,and 10o).Thus, the morphology of ENA emissions at Europa can only be used to constrain the ambient PAD if precise knowledge of the atmospheric structure, the upstream energetic proton distribution  Ĩ() , and the draped magnetospheric field configuration is available.11g-11i).

ENA Emissions at
The presence of the H 2 O exosphere only slightly alters the draped field configuration near Europa (see Section 3.1 and Figure 4).Therefore, energetic parent protons are granted access to Europa's atmosphere in a very similar fashion to the configuration that only includes an O 2 component.Thus, the ENA emission morphology shares a strong resemblance between the two cases: differences are confined to the region of the detector above the confined H 2 O profile.Since energetic proton gyroradii at Europa are 1-2 orders of magnitude smaller than the moon's radius (see Table 2), absorption by Europa makes the detector's wakeside hemisphere largely inaccessible to ENAs generated near the ramside apex.The charge exchange cross section between protons and H 2 O (similar to protons and O 2 ) decreases steadily with increasing energy (see Figure 2).Considering that the ambient differential proton flux peaks around 60 keV, the detected ENA flux at Europa (both with and without the water vapor) is therefore weakest at an initial energy of 1 keV, peaks in the model results at an initial energy of 40 keV, and then falls off monotonically with increasing initial energy (see left and center columns of Figure 11).
The number density in Europa's ramside atmosphere is elevated by up to 43% with the additional H 2 O profile included, but the ENA flux around the detector's ramside apex is kept below a 20% increase at all energies studied.As parent ions propagate through the (now denser) atmosphere, their numerical weight is attenuated by the two neutral species more rapidly than by the O 2 component alone.Therefore, while a proton traveling through the combined O 2 /H 2 O atmosphere may produce more ENAs at higher altitudes than with O 2 alone, its ability to generate ENAs in the more dense lower atmosphere is reduced.Thus, the quantitative changes to the ENA flux between both cases are rather subtle, and the detection of such a confined neutral profile at Europa would rely upon highly constrained upstream parameters in tandem with direct knowledge of the energetic proton flux.
In addition to a persistent H 2 O atmosphere around Europa's ramside apex, several remote and in-situ observations have revealed localized, transient plumes of H 2 O vapor emanating from Europa's surface (e.g., Arnold et al., 2019;Jia et al., 2018;Roth et al., 2014b).Plumes generate a localized number density near the surface that is only about a factor of 4 larger than the atmospheric H 2 O density used in this study (Roth, 2021;Roth et al., 2014b).The scale height of the H 2 O density in such a plume is estimated to be approximately 200 km (e.g., Arnold et al., 2019), that is, only slightly larger than those of the persistent O 2 and H 2 O components.Thus, our results from Figure 11 suggest that such a plume would manifest at Europa only through a subtle quantitative enhancement of the detected ENA flux above its source region.Such a feature would likely be difficult to discern from the "background" ENA emissions generated by protons that undergo charge exchange with the O 2 and H 2 O components of the moon's global atmosphere.

Variability of ENA Emissions at Callisto
In Figure 12, we present maps of the ENA flux through a concentric spherical detector of radius 2R C around Callisto when the moon is at the center of the Jovian plasma sheet.The layout of Figure 12 is analogous to that of Figure 5, with the left and right columns presenting the detectable ENA flux in uniform and draped electromagnetic fields, respectively.Again, each row corresponds to a distinct initial energy of the parent protons.When Callisto is at the center of the Jovian plasma sheet, the magnetospheric field in the vicinity of the moon is oriented approximately southward, and therefore we do not include an induction signal in this treatment (see Table 1).Since the strength of the magnetospheric field at the center of the plasma sheet falls around two orders of magnitude when moving outward from Europa to Callisto, the parent proton gyroradii at Callisto are approximately two orders of magnitude larger at each respective initial energy, becoming comparable to or larger than 1R C .For instance, at E = 10 keV, the gyroradius of a proton with α = 90° is 1.5R C ≈ 3, 615 km at Callisto but only 0.02R E ≈ 31 km at Europa.
As with Europa, we investigate the (hypothetical) case when Callisto is in uniform magnetospheric fields in order to obtain a baseline understanding of the ENA emission patterns.For this scenario, the bulk of the detectable   1).The figure is arranged in the same manner as Figure 5.The left and right columns show results for uniform and perturbed electromagnetic fields, respectively.Note that the colorscale covers a much broader range here than the scales used in Figures 5 and 9.
ENA flux is focused onto the equatorial region of the detector (left column of Figure 12).This band of high flux is most distinct in Figures 12c, 12e and 12g.The equatorial band wraps around Callisto normal to the orientation of the magnetospheric background field, and its latitudinal extension can again be estimated using Equation 21and the scale height  O 2 , = 230 km of the moon's molecular oxygen envelope.We determine an estimate of θ b = 33.2°for the half-width of the band, which is similar to the latitudinal extent of the modeled band on the detector (about θ b ≈ 36°; see, e.g., Figure 12e).
The boundary of the high flux band at Callisto appears blurred when the moon is at the center of the Jovian plasma sheet.Regions receiving significant flux on the detector reach outside of the estimated boundaries of the band up to latitudes of ±55°, particularly visible as the yellow region in Figure 12c.At Europa, the band of high ENA flux has sharply pronounced latitudinal boundaries because parent protons with pitch angles near 90° trace much longer paths through the atmosphere than those that largely travel along the magnetospheric field lines (see Section 3.2).This difference in path lengths stems from the fact that, due to the small gyroradii, parent protons can complete many gyrations within Europa's atmosphere.In contrast, when Callisto is at the center of the plasma sheet, even the gyroradii of 1 keV protons exceed the atmospheric scale height  O 2 , by a factor of 5. Therefore, ENA generation from individual proton macroparticles is limited to small contributions from partial gyrations within the atmosphere.The left column of Figure 13 illustrates the trajectories of two protons that spend only a small portion of their gyrations within the atmosphere.Both protons are initiated with an energy of 10 keV and a pitch angle of α = 90° upstream of the moon.They propagate through uniform electromagnetic fields, remain confined to the z = 0 plane (Figure 13b), and complete two partial gyrations within the atmosphere, each before impacting Callisto.Each of these protons travels a total path length through the atmosphere of no more than one gyration circle, that is, 2πr g ≈ 6.3R C (see Figure 13a).This length is comparable to the longest possible path of a proton with α = 0° through the atmosphere in uniform fields:  2 √ 3  , or the distance a proton travels along a magnetic field line tangent to Callisto's surface while within the detector sphere.This example demonstrates why the band (bounded by latitudes ± 36°) contains only 50%-60% of the total detectable ENA flux at each initial parent ion energy at Callisto (see left column of Figure 12).In comparison, over 80% of the total detectable ENA flux was concentrated inside the band at Europa (see left column of Figure 5).
Across the studied energy range, the detectable ENA flux in uniform fields is maximized slightly downstream of the detector's Jupiter-facing apex.The minimum flux within the band is detected slightly upstream of the Jupiter-averted apex (see, e.g., Figure 12c).For example, at an initial parent proton energy of 40 keV, the ENA flux reaches its maximum of over  37, 000 [ cm 2 sr keV s ] −1 at 24° W longitude near the northern edge of the band.
The minimum flux detected within the band at the same energy is only 45% the intensity of the maximum, and it is clustered around approximately 200° W longitude.In contrast to this, the maxima and minima of the ENA flux in uniform fields at Europa occur close to the detector's ramside and wakeside apices, respectively (see Figure 5).
To understand the locations of the ENA flux extrema at Callisto, we consider again the energetic proton trajectories (solid and dashed lines) depicted in Figure 13a.Of the detectable ENA flux produced by the sample proton initialized in the y < 0 half space (solid line), the ENA flux emitted during the first trajectory segment within the atmosphere (solid green line) is mostly directed toward the Jupiter-facing, upstream quadrant of the detector.The flux generated during the proton's second incursion into the atmosphere is absorbed by Callisto.A complementary proton trajectory initialized in the y > 0 half space also completes two partial gyrations within the atmosphere (dashed green line).
During the first trajectory segment within the atmosphere, this proton only emits detectable ENAs toward the y > 0 half space.During the second segment, all its emitted ENAs are absorbed by Callisto.This example suggests that protons with pitch angles near 90° emit detectable ENA flux preferentially toward the Jupiter-facing hemisphere of the detector irrespective of their initial y coordinate, due to the direction of their gyration and their large gyroradii (see Figures 12c, 12e, and 12g).If the ambient magnetospheric field orientation were inverted, ENAs would preferably be emitted into the y < 0 half space.If the ambient magnetic field magnitude were enhanced, the rotation of the extrema in ENA flux away from the x axis would be weakened (see Section 3.2).For this reason, the locations of the emission maxima/minima at Europa (where |B 0 | is 1-2 orders of magnitude stronger) differ from those at Callisto.
Our modeled ENA emission patterns for Callisto are consistent with the findings of Tippens et al. (2022) at Titan for a southward magnetospheric field, where the energetic proton gyroradii are also comparable to the size of that moon.Using a similar model, these authors found the band of high ENA flux to persistently form normal to the ambient Kronian magnetospheric field, and they determined an enhancement of over 50% in the ENA flux intensity around the Saturn-facing apex of the detector compared to the Saturn-averted apex.In the y = 0 plane (panel (c)), the trajectory of this proton overlaps with that of the particle launched at y < 0. For the sake of clarity, the corresponding trajectory in draped fields is not displayed.The blue arrows again show an approximate bearing of the guiding center motion as the trajectory enters the panel frame.
As shown in the right column of Figure 12, the observable ENA flux at Callisto is reduced in intensity by the field line draping, compared to a uniform southward field.The detected maximum in ENA flux remains near the Jupiter-facing apex when the plasma interaction is considered.Similar to the case of draped fields at Europa (see right column of Figure 5), formation of a band signature of elevated ENA flux along the equator does not occur.Instead, the ENA flux is now diminished in almost the entire Jupiter-averted hemisphere of the detector.When draped fields are included, this hemisphere receives at least 53% less ENA flux at all initial proton energies, compared to reductions by only 15%-40% in the Jupiter-facing hemisphere of the detector.Near the center of the Jovian plasma sheet, Callisto's tail of ionospheric pickup ions is highly asymmetric due to their large gyroradii (1 − 10R C ), extending several R C into the Jupiter-averted half space (see Figure 3i in this study and Figure 2e of Liuzzo et al., 2022).The ramside magnetic pileup region is "stretched" along the upstream flank of this pickup tail, thereby reaching several R C into the Jupiter-averted hemisphere (Liuzzo et al., 2015(Liuzzo et al., , 2022)).Deflection of the impinging parent ions by the field enhancement in the ramside pileup region further diminishes the accessibility of Callisto's Jupiter-averted atmosphere to energetic protons.A sample trajectory that illustrates proton deflection by the stretched pileup region is displayed in the right column of Figure 13. Figure 13b demonstrates how the trajectory is deflected strongly in the (−y) direction when approaching Callisto, entirely avoiding the atmosphere and generating no ENA flux.In addition, this trajectory in uniform fields is confined to the z = 0 plane, but in draped fields the pileup deflects the proton out of the z = 0 plane (see Figures 13c and 13d).
In both uniform and perturbed fields, the detectable ENA flux is minimized at an initial proton energy of 1 keV, peaks around 10-40 keV, and begins to fall off as energy increases beyond 40 keV (see Figure 12).For our studied proton energies, the detectable ENA flux reaches its peak value of  55, 700 [ cm 2 sr keV s ] −1 in uniform fields at 10 keV.The differential proton flux outside of Callisto's interaction region peaks between 10 and 20 keV (see, e.g., Figure 3 of Liuzzo et al., 2022).The product  () ⋅ Ĩ() maximizes at comparable initial proton energies, consistent with the behavior of intensity as a function of initial parent ion energy in our ENA emission maps.
At the center of the Jovian plasma sheet, the ENA emission intensity at Callisto exceeds that at Europa at each initial proton energy (in uniform and draped fields) by a factor of 20-40 (compare Figures 5 and 12).The ambient proton flux at Callisto's orbit differs from that observed at Europa by less than an order of magnitude across the studied energy interval (Liuzzo et al., 2022;Paranicas et al., 2002).Hence, the difference in ENA flux intensity is largely driven by the disparity in atmospheric density: Callisto's atmosphere possesses an O 2 column density 2-3 orders of magnitude greater than the column density of Europa's O 2 envelope (Carberry Mogan et al., 2021;Cunningham et al., 2015;de Kleer et al., 2023;Liuzzo et al., 2015;Roth et al., 2016).Parent protons interacting with Callisto's atmosphere encounter a higher local density and therefore generate more ENA flux than at Europa.In addition, the weaker ambient magnetic field at Callisto allows energetic parent ions to gyrate on much larger scales.As a result, protons whose guiding center trajectories do not even come close to intersecting Callisto's atmosphere may still penetrate into the moon's neutral envelope and produce ENAs.This is not the case at Europa.This is further illustrated in Figure 14, displaying the energy dependence of proton gyroradii at (a) Europa and (b) Callisto in the uniform magnetospheric background fields at the center of and below the Jovian plasma sheet.In contrast to Table 2, the gyroradii are now normalized to the atmospheric scale height at the respective moon.As can be seen, gyroradii at Europa remain smaller than or comparable to the scale height, regardless of the moon's distance to the plasma sheet.However, when Callisto is located at the center of Jupiter's plasma sheet, proton gyroradii in the weak magnetospheric field may exceed the scale height by more than an order of magnitude.
In Figure 15, we present the ENA flux maps when Callisto is situated at its maximum distance below the center of the Jovian plasma sheet.Since the ambient magnetospheric field B 0 is nearly parallel to the z = 0 plane and forms an angle of less than 17° against the y axis, we include an induced dipole moment at Callisto for this setup.Similar to Section 3.2, we compare the detectable ENA emissions for a constant magnetospheric background field (left column), the constant background field superimposed with the induced dipole (center column), and the perturbed electromagnetic fields resultant from Callisto's interaction with the Jovian magnetosphere (right column).When Callisto is at its maximum distance below the sheet, the vast majority of ENA flux is again focused into a band of high flux on the detector, again centered about the great circle perpendicular to the background field, now wrapping around both geographic poles of the detector.The regions depleted of ENA flux are once again found above where the magnetospheric background field is largely normal to Callisto's surface.However, due to the magnetospheric field orientation, the flux is now depleted around the Jupiter-facing and Jupiter-averted apices (see Figure 15).10.1029/2023JA031931 39 of 45 sheet (see Figures 3a, 3b, 3e, and 3f).As a result, the detectable ENA flux patterns shown in the center and right columns of Figure 15 are highly similar.At Callisto, the variability of the plasma interaction over a synodic rotation is large compared to that at Europa (see Figure 3).Hence, this variability plays a greater role in influencing the morphology of ENA emissions at Callisto: only at this moon does the appearance of the high flux band change drastically between the "below" and the "center" case (see Figures 12 and 15).Further, the effect of the plasma interaction may be weakly discernible in detectable ENA emissions emanating from Callisto's atmosphere when the moon is far outside the Jovian plasma sheet.
Though not yet directly confirmed through observations, the model of Carberry Mogan et al. (2022) suggests the presence of an H 2 atmospheric component at Callisto which is not included in our current model setup.Figure 11 of that study implies that this component may have a column density comparable to that of O 2 .Charge exchange cross sections between protons and molecular hydrogen (Tawara et al., 1985) are similar to those for oxygen molecules, which would indicate a non-negligible contribution to ENA production.Due to the lower mass, the scale height of the H 2 component would be elevated compared to that of the O 2 included in our model.The presence of this component may therefore increase ENA production at high altitudes (above the extension of the O 2 atmosphere).However, for a detector that encompasses Callisto's entire atmosphere, the ENA emission  21suggests that inclusion of such an H 2 component would increase the latitudinal width of the band on the detector.Besides, the H 2 envelope may increase the path length spent by the energetic protons within Callisto's atmosphere.Building upon our conclusions for Europa (see Section 3.2), this may rotate the maxima/minima of the detectable ENA emissions closer toward the ram/wake apices, respectively.However, overall we do not expect this component to introduce distinct, novel morphological features into our ENA emission maps.Since in contrast to the H 2 O bulge at Europa-the substantial presence of H 2 at Callisto has not yet been verified through observations, it is not included in the current iteration of our model.

Summary and Concluding Remarks
In this work, we have studied the intensity and morphology of ENA emissions at Callisto and Europa.Charge exchange between energetic magnetospheric ions and cold atmospheric neutrals results in the production of such ENAs.Our modeling framework is based upon a combination of a tracing tool for energetic magnetospheric ions and the AIKEF hybrid model (Müller et al., 2011) which calculates the three-dimensional structure of the draped electromagnetic fields near both moons.We trace billions of energetic parent ions through draped fields near the two moons.We calculate maps of the ENA flux they emit through a spherical detector centered at each moon and encompassing its entire atmosphere.Such a detector allows us to ascertain a complete physical picture of detectable ENA emissions from Europa and Callisto which could not be simultaneously captured by a point-like spacecraft detector with a limited field of view.
At both moons, we have constrained the variability of detectable ENA emissions across a Jovian synodic rotation by comparing maps of the ENA flux when each moon is either at the center or far outside the Jovian plasma sheet.The magnetospheric field line draping at Callisto and Europa gradually weakens with distance to the center of the plasma sheet.For each configuration, we have analyzed the influence of field line draping on parent ion dynamics and the resultant ENA emissions by comparing to a "baseline" scenario with uniform electromagnetic fields.In order to constrain the detectable ENA flux morphology as a function of parent ion energy, we have calculated maps of the ENA flux through the detector sphere at five initial proton energies ranging from 1 to 100 keV.As such, we predict the global distribution of ENA flux from Callisto's and Europa's atmospheres as will be observable during the upcoming JUICE mission (Grasset et al., 2013).
We report our major results as follows: 1.For a spacecraft located outside of Europa's or Callisto's atmosphere, detectable ENA emissions are concentrated into a band of high flux, centered about the great circle on the detector sphere normal to the ambient Jovian magnetospheric field.The latitudinal width of this band is determined by the extension of the obstacle to the flow along the magnetospheric field lines, approximately given by the sum of the moon's diameter and twice the atmospheric scale height.During the course of a synodic rotation, the orientation of this band on the detector sphere oscillates in phase with the time-varying background magnetic field vector.2. At the center of the Jovian plasma sheet where field line draping around Europa and Callisto is strongest, the band of high ENA flux is attenuated in intensity regardless of longitude on the spherical detector.In this case, the bulk of the observable ENA emissions is focused into a cluster of elevated ENA flux near the detector's ramside apex (Europa) or Jupiter-facing apex (Callisto).When either moon is located at its maximum distance below the center of the Jovian plasma sheet, field line draping again attenuates the intensity of the observable ENA emissions, but the emission patterns recorded on the detector are qualitatively similar for uniform and draped fields.3. The fraction of the total observable ENA flux that is confined to the high flux band is determined by the ratio of energetic proton gyroradii and the moon's atmospheric scale height.At Europa, regardless of its distance to the Jovian plasma sheet, parent ions drifting with pitch angles near 90° can trace a path through the atmosphere 10-100 times longer than that of ions traveling along the magnetic field.This results in ≥80% of the ENA flux emitted from Europa's atmosphere toward the detector being concentrated within the high flux band.When Callisto is at the center of the Jovian plasma sheet, energetic proton gyroradii are a factor of 5-50 larger than the atmospheric scale height.In consequence, protons with pitch angles near 90° can only complete partial gyrations within the atmosphere.Therefore, the band contains only 50%-60% of the total ENA flux leaving Callisto's atmosphere.A significant population of ENAs hit the detector above or below the band.Far outside the Jovian plasma sheet, the ratio of parent ion gyroradii and Callisto's atmospheric scale height is similarly small as at Europa.Therefore, the observable ENA emissions are again largely focused into a band perpendicular to the background field.4. At both moons, a local maximum and a local minimum in ENA intensity form within the band of high flux, separated on the spherical detector by approximately 180° in longitude along the great circle that defines the band's equator.At Europa, the maximum ENA flux is detected slightly west of the ramside apex.When gyrating ions drift into the moon's atmosphere from upstream, they emit ENAs toward Jupiter (i.e., west of the ramside apex) if they gyrate deep into the dense, lower atmosphere.They emit ENAs away from Jupiter (i.e., east of the ramside apex) if their gyration carries them back toward upstream, and they undergo charge exchange with the more rarefied neutral gas at higher altitudes.Since the likelihood of charge exchange is proportional to the local neutral density, a larger ENA flux is detectable west of the ramside apex than east of it.When Callisto is at the center of the plasma sheet, the maximum in ENA flux through the detector is further shifted, now located close to the Jupiter-facing apex.In this scenario, the magnetospheric field strength at Callisto is 1-2 orders of magnitude weaker than at Europa.Therefore, ion gyroradii at Callisto are large compared to the moon's atmospheric scale height.Due to their sense of gyration, parent ions mostly emit ENAs along trajectory segments bearing toward Jupiter, thereby causing the longitudinal shift in the location of the emission maximum.If the heavier magnetospheric ion species (oxygen, sulfur) were included in our Europa ENA model, their larger gyroradii would make the emission maps appear more similar to those calculated for protons at Callisto.However, the occurrence of the band-like structure would be largely unaffected.5.The anisotropic PAD observed for protons at Europa by Juno (Sarkango et al., 2023) carries elevated particle fluxes at pitch angles near 0° and 180°.Including such an anisotropic PAD in our model leaves the morphology of the detectable ENA flux largely unchanged, compared to a setup using the approximately isotropic PAD observed by Galileo (Mauk et al., 2004).Including the anisotropic PAD increases the observable ENA flux by up to 25% in the mid-latitude regions of the detector's ramside hemisphere.The detectable ENA flux within the high flux band is diminished from the anisotropic PAD because it contains fewer parent ions with pitch angles near 90° than the isotropic one.6.A localized H 2 O exosphere around Europa's ramside apex (Roth, 2021), in addition to the moon's global O 2 atmosphere, produces only minor quantitative enhancements to the observable ENA flux.This suggests that even more localized and compact populations of water vapor (e.g., a plume) are likely difficult to detect against the "background" of ENA flux from Europa's global atmosphere, unless they are located in regions where the magnetic field is perpendicular to the moons surface (i.e., outside of the band) and their scale height greatly exceeds that of the persistent atmosphere.
The natural next step building upon this study is to incorporate a realistic, point-like spacecraft ENA detector with a limited field of view into our ENA generation model.In this way, we can perform predictive modeling of observable ENA emissions during scheduled JUICE flybys of Callisto and Europa.Our work also forms the basis for assessing the relative contributions of ENAs emitted from the moons' atmospheres and surfaces (Pontoni et al., 2022) in a future study.

Figure 1 .
Figure1.A depiction of the energetic neutral atom generation model at Callisto, for the case with the moon located at the center of the Jovian plasma sheet.The downstream, northern, and Jupiter-averted faces of the starting grid cuboid (white mesh on gray background) are omitted to allow viewing into the model domain, and the radius of the detector (purple sphere around the atmosphere) is enlarged for clarity.At Callisto, a step size of Δx s = Δy s = Δz s = 0.04R C is used in each of the six planes of the starting grid.A single parent ion trajectory (light blue) is displayed, initiated at a grid node on the upstream face of the starting grid.This parent ion remains within Callisto's atmosphere (orange shaded region) for nearly an entire gyration and generates ENAs.ENAs that collide with the moon (red) are not detected, but ENAs with a velocity vector bearing away from the moon (green) are observed by the spherical model detector that encompasses the moon's entire atmosphere.The figure is not to scale.
and Tippens et al. (2022)): such particles can not return to the moon and produce ENAs.Energetic protons exiting the domain are not immediately removed since they may still enter Europa's and Callisto's atmospheres later in time.The uniform magnetospheric background field B 0 and the corresponding convective electric field E 0 = −u 0 × B 0 are used to evolve particle trajectories outside of the AIKEF domain.An adaptive timestep is implemented to solve Equation 17 dependent both on proximity to the moon as well as the gyroperiod   = 2  || in the local magnetic field B. The timestep Δt in the energetic particle tracer is given by

Figure 2 .
Figure 2. Charge exchange cross sections for protons with molecular oxygen and water vapor, given by Lindsay and Stebbings (2005) and Wedlund et al. (2019), respectively.The cross section for  H + + O2 → H + O + 2 is only available for energies 0.05 keV ≤ E ≤ 100 keV; that is, covering the energy range that we analyze in our study.
Figures 3 and 4 depict the structure of the perturbed electromagnetic environment in the vicinity of Europa and Callisto for each hybrid model configuration.The left two columns of Figure3correspond to Callisto at the center and at maximum distance below the center of the Jovian magnetospheric plasma sheet, respectively.The right two columns follow identical organization for Europa.The first row (Figures3a-3d) and second row (Figures3e-3h) show the flow-aligned magnetic field component B x and field strength  || in the y = 0 plane, respectively.The third row (Figures3i-3l) displays the magnitude of the electric field |E| in the z = 0 plane.All quantities in Figure3are normalized to the respective upstream value, that is, |B 0 | for the top two rows and |E 0 | for the bottom row.When either moon is at the center of the Jovian plasma sheet, the y = 0 plane contains the background magnetic field vector B 0 , and the z = 0 plane contains the background convective electric field vector E 0 .When Europa is located at its maximum distance below the Jovian plasma sheet, the background magnetic field is rotated against the y = 0 plane (toward Jupiter) by 28.5°, and the background electric field is correspondingly rotated by the same angle against the z = 0 plane.Analogously, when Callisto is located below the plasma sheet, the background magnetic field is rotated by 72.5° against the y = 0 plane toward Jupiter (see Table1).For this case, a depiction of plasma and electromagnetic field quantities in planes containing the strongly inclined background fields is given in Figure1ofLiuzzo et al. (2022).Figure4shows results for Europa when an H 2 O atmospheric component is included in addition to molecular oxygen.We include in Figures4a and 4bthe H 2 O + number density in the z = 0 and y = 0 planes, respectively.The bottom row of Figure 4 displays the electromagnetic field quantities.Extensive discussions of the perturbed electromagnetic environment within the framework of AIKEF are given by Arnold et al. (2019), Arnold, Liuzzo, and Simon (2020), Arnold, Simon, and Liuzzo (2020), Breer et al. (2019),

Figure 3 .
Figure 3.The perturbed electromagnetic environments around Callisto (left two columns) and Europa (right two columns) in planes containing the center of the moon, computed using the AIKEF hybrid model.Panels (a)-(d) show the flow-aligned component (B x ) of the magnetic field while panels (e)-(h) show the magnetic field magnitude |B| near the respective moon, both in the y = 0 plane.Panels (i)-(l) show the magnitude of the electric field |E| in the z = 0 plane.When either moon is located at the center of the Jovian plasma sheet (when they are exposed to a purely southward magnetospheric background field), the convective electric field E 0 points in (−y) direction and is therefore parallel to the z = 0 plane.Each quantity is plotted on a scale normalized to the respective background field value (i.e., |B 0 | for the top and middle row, |E 0 | = |u 0 × B 0 | for the bottom row, see Table1).Within each column, the left sub-column corresponds to the respective moon at the center of the Jovian plasma sheet, and the right sub-column corresponds to the moon at its maximum distance below the plasma sheet.To ensure best possible visibility of the field perturbations, the ranges of the color bars are different between panels (a) and (b) as well as panels (e)-(h).

Figure 4 .
Figure 4.The perturbed plasma and electromagnetic environment near Europa when the moon is located at the center of the Jovian plasma sheet, including the localized H 2 O exosphere near the moon's ramside/subsolar apex.In panels (a) and (b), the number density of ionospheric H 2 O + is displayed in the z = 0 and y = 0 planes, respectively.The flow-aligned component of the magnetic field B x is plotted in panel (c) in the y = 0 plane.The magnitude of the magnetic field within the same plane is displayed in panel (d).The magnitude of the electric field is plotted in panel (e) for the z = 0 plane, containing the upstream flow velocity u 0 , the background convective electric field vector E 0 = −u 0 × B 0 , and the center of the moon.All quantities are normalized to the background magnetospheric quantities at Europa, that is, n 0 for panels (a) and (b) and |B 0 | for panels (c) and (d).The electric field magnitude (panel (e)) is normalized to |E 0 | = |u 0 × B 0 |.
the northern Alfvén characteristic   − has a component out of the plane of the figure (i.e., it points away from Jupiter), and the southern characteristic   + has a component into the plane (toward Jupiter).At Europa, the induced dipole moment M ind is (nearly) aligned with the (−y) axis.The vector M ind generated by the nonzero B h forms an angle of 14° with the (−y) axis at Callisto, that is, it is tilted toward upstream.At Callisto (Figure3b), an enhancement in B x by |B 0 |/5 is visible above both the north and south pole of the moon, and a reduction of similar magnitude is present near the equator in the y = 0 plane.In this setup, the background field possesses a positive B x,0 component of approximately |B 0 |/4.The induced dipole field takes the form of a "shamrock" in B x viewed in the z = 0 plane (see also Figures1b and 1jofLiuzzo et al., 2022).If B h were completely aligned with  − ŷ , the B x perturbation associated with the induced dipole is positive when x and y have opposite signs, and negative when both x and y have the same sign.A background field B x,0 > 0 rotates the shamrock structure in B x counterclockwise around the z axis.Therefore, the induced field at Callisto is characterized by B x < 0 along the x axis, acting to reduce the ambient B x,0 component.This rotation similarly causes an enhancement in B x near the poles of Callisto in the y = 0 plane.

Figure 5 .
Figure5.Global energetic neutral atom (ENA) emission flux maps from Europa's atmosphere through a concentric spherical detector with radius 2R E .In this scenario, Europa is located at the center of the Jovian plasma sheet (see Table1).The West Longitude coordinate system is used in this projection.Each row represents ENA emissions generated by parent ions initialized on the starting grid at a certain kinetic energy, namely 1 keV in row (a-b), 10 keV in row (c-d), 40 keV in row (e-f), 70 keV in row (g-h), and 100 keV in row (i-j).The left column (panels (a), (c), (e), (g), and (i)) illustrates the ENA emission morphology when the energetic parent protons are evolved through the uniform Jovian background fields B 0 and E 0 .The right column (panels (b), (d), (f), (h), and (j)) displays ENA emissions when the parent ions move through the perturbed electromagnetic fields from the adaptive, kinetic ions, electron fluid hybrid model (see panels (c), (g), and (k) of Figure3).

Figure 6 .
Figure 6.A sample energetic proton trajectory and energetic neutral atom (ENA) emission statistics highlighting the effect of parent ion gyration within the atmosphere on the morphology of ENA emissions at Europa.The parent ion is initialized at  ( =  0) = (−6R E , 0, 0) with energy E = 40 keV, pitch angle α = 90°, and initial velocity  ̇( =  0) =  ( 0, √ 2   , 0 ) , and moves through uniform electromagnetic fields.The projection of the particle trajectory into the z = 0 and y = 0 planes is displayed in panels (a) and (b), respectively, with the blue arrow indicating the approximate direction of guiding center motion.The red rectangle contains the segment of the trajectory from 71 to 73.27 s, as shown in panels (d) and (e).The orange dashed circle represents the outer boundary of the moon's atmosphere in our model, coinciding with the detector sphere.The trajectory is depicted in cyan outside of the atmosphere (where no ENAs are produced) and green inside the atmosphere (where ENA generation takes place).In panel (c), the total ENA flux detected on the Jovian-oriented hemisphere of the detector (i.e., y > 0) is shown in cyan and compared to the total ENA flux detected by the Jovian-averted hemisphere of the detector (i.e., y < 0), shown in gold.Note, only those ENAs which do not bombard Europa itself are actually detected and included in panel (c).The coordinate on the horizontal axes of panels (d) and (e) refers to the time that has elapsed since the ion was initiated at  ( = 0), see Equation 22.The background colors in panels (d) and (e) show  sign ( ̇ ) for the parent proton, with the time window shaded as either cyan when positive or gold when negative.In segments with   > 0 (or   < 0 ), the proton emits ENAs toward the Jupiter-facing (Jupiter-averted) hemispheres of the detector.Panels (d) and (e) display the atmospheric neutral density at the location of the proton and the differential ENA flux Y =  |d | emitted at each timestep as the parent proton travels through the atmosphere, respectively.In panel (e), the emitted ENA flux is normalized to the maximum value  max ≈ 0.02 [ cm 2 sr keV s ] −1 produced around t = 73.27s.Beyond that time (not shown in the figure), the emitted ENA flux during each gyration decreases since the parent proton has already depleted most of its initial numerical weight.

Figure 7 .
Figure 7. Sample energetic parent protons' evolution through Europa's electromagnetic environment.The 3D ion trajectories are projected into both the z = 0 (panels (a) and (b)) and y = 0 planes (panels (c) and (d)).We include two separate proton trajectories with the same initial conditions, traced through uniform (left column) and draped fields (right column).The arrow in each panel indicates the direction of guiding center motion during the ion's initial approach of Europa.When outside of the atmosphere, the proton trajectory is plotted in cyan (uniform fields) or magenta (draped fields).Both cases are colored in green when they are inside the atmosphere and generate ENAs.For both uniform and draped fields, the sample ion is initialized with  ( =  0) = (0.8RE , −0.1R E , 6R E ), that is, outside of the domain shown in the Figure 7. Sample energetic parent protons' evolution through Europa's electromagnetic environment.The 3D ion trajectories are projected into both the z = 0 (panels (a) and (b)) and y = 0 planes (panels (c) and (d)).We include two separate proton trajectories with the same initial conditions, traced through uniform (left column) and draped fields (right column).The arrow in each panel indicates the direction of guiding center motion during the ion's initial approach of Europa.When outside of the atmosphere, the proton trajectory is plotted in cyan (uniform fields) or magenta (draped fields).Both cases are colored in green when they are inside the atmosphere and generate ENAs.For both uniform and draped fields, the sample ion is initialized with  ( =  0) = (0.8RE , −0.1R E , 6R E ), that is, outside of the domain shown in the figure, and  ̇( = 0) =  (√ 2   sin 5 • , 0, − √ 2   cos 5 • ) , with a kinetic energy of E = 40 keV.

Figure 8 .
Figure 8. Sample energetic protons' evolution through Europa's environment, depicting a common fate shared by ions initialized upstream of the moon with pitch angles of 90°.The layout of this figure is identical to that of Figure 7, with analogous coloring.The test particles are initialized with  ( =  0) = (−6R E , −0.1R E , −0.5R E ), that is, outside of the domain plotted here, and  ̇( =  0) =  (√ 2   sin 90 • , 0, − √ 2   cos 90 • ) at an energy of E = 40 keV.

Figure 9 .
Figure 9. Global energetic neutral atom (ENA) flux maps from Europa when at its maximum distance below the center of the Jovian plasma sheet.The spherical detector is again concentric about the moon and located at |r| = 2R E .Analogous to Figure 5, each row represents ENA emissions generated by parent ions initialized at a certain kinetic energy, namely 1 keV in row (a-c), 10 keV in row (d-f), 40 keV in row (g-i), 70 keV in row (j-l), and 100 keV in row (m-o).Moving from left to right, the columns display the ENA flux maps in uniform magnetospheric background fields (panels (a), (d), (g), (j), and (m)), uniform magnetospheric fields superimposed with Europa's induced dipole (panels (b), (e), (h), (k), and (n)), and including the contribution of currents from the plasma interaction (panels (c), (f), (i), (l), and (o)).The color scale is identical to that used in Figure 5 to facilitate comparison.
Figure 10 displays the global ENA emissions through the detector sphere when Europa is located at the center of the Jovian plasma sheet, with two different PADs used to initialize the parent protons on the starting grid.In both cases, the ENA flux maps have been calculated using the draped fields from AIKEF, and the model results for an isotropic PAD (see Figure 5) are shown again in the left column to aid comparison.ENA emissions generated by parent protons with an anisotropic PAD (see Equations 10 and 12) are given in the center column.The maps in the right column have been generated by subtracting the ENA flux maps for an anisotropic PAD (center column) from those calculated with an isotropic PAD (left column) at each respective initial parent ion energy.Positive values (red) then show regions where the anisotropic PAD weakens the detectable ENA flux, and negative values (blue) display where the anisotropic PAD enhances the ENA flux.

Figure 10 .
Figure 10.Influence of an anisotropic pitch angle distribution (PAD) on the global energetic neutral atom (ENA) emission morphology at Europa when situated at the center of the Jovian plasma sheet.The figure layout is analogous to that of Figure 5, where each row represents an initial parent proton energy.The left column displays the ENA flux with an isotropic PAD in draped electromagnetic fields (panels (a), (d), (g), (j), and (m)) to facilitate comparison, also shown in the right column of Figure 5.The center column (panels (b), (e), (h), (k), and (n)) displays detected ENA flux in the same field configuration, but using the anisotropic PAD given by Equation10.The right column (panels (c), (f), (i), (l), and (o)) then shows the difference between the detected ENA flux with an isotropic (left column) and with an anisotropic (center column) PAD, that is, the left column minus the center column.The white regions in maps of the right column correspond to locations on the detector toward which the isotropic and anisotropic PAD emit the same detectable ENA flux.Red (positive) or blue (negative) indicate regions where the isotropic or anisotropic PAD generate more detectable ENA flux, respectively.Note that the range of the right column's colorscale is reduced relative to the scale used in the left and center columns.

Figure 11 .
Figure 11.Influence of a localized atmospheric H 2 O component around Europa's ramside apex on the global energetic neutral atom (ENA) emission maps at the center of the Jovian plasma sheet.The ENA flux is displayed analogously to Figures 5, 9, and 10, where each row represents an initial parent ion energy.The left column again displays the ENA flux in draped electromagnetic fields (adopted from Figure 5) without the neutral H 2 O component (panels (a), (d), (g), (j), and (m)).The center column (panels (b), (e), (h), (k), and (n)) displays detected ENA flux in the slightly modified draped fields (see Figure 4), now including ENA generation through interaction with the additional H 2 O component.The right column (panels (c), (f), (i), (l), and (o)) then shows the difference between the detected ENA flux with a purely O 2 atmosphere (left column) and an atmosphere including both the O 2 and the H 2 O component (center column).The white regions in the right column correspond to locations where the detected ENA flux is identical with and without the H 2 O included.Regions in red show where the ENA flux through the detector is higher without the atmospheric H 2 O, and regions in blue show where the detector receives higher ENA flux from an atmosphere including both O 2 and H 2 O.Note that the scale of the right column covers a much narrower range of values than the scale used in the left and center columns.

Figure 12 .
Figure12.Global energetic neutral atom emission flux maps at Callisto through a spherical detector centered at the origin with radius 2R C .In this setup, Callisto is located at the center of the Jovian plasma sheet (see Table1).The figure is arranged in the same manner as Figure5.The left and right columns show results for uniform and perturbed electromagnetic fields, respectively.Note that the colorscale covers a much broader range here than the scales used in Figures5 and 9.

Figure 13 .
Figure 13.Two-dimensional projections of three representative parent proton trajectories through Callisto's electromagnetic environment.Each parent proton is initialized at an energy of E = 10 keV and travels through either uniform fields (panels (a) and (c)) or perturbed fields (panels (b) and (d)).Projections of the proton trajectories onto the z = 0 and the y = 0 planes are displayed in the top and bottom rows, respectively.The coloring is analogous to Figures 7 and 8, using light blue for the parent proton traced in uniform fields, magenta for the parent proton in draped fields, and green when either trajectory is embedded within the atmosphere (i.e., |r| ≤ 2R C ) and emitting ENAs.The parent protons with a solid trajectory are propagated from outside of the plot window through both uniform and draped fields, with an initial position given by  ( =  0) = (−10R C , −3.3R C , 0R C ) and an initial velocity of  ̇( = 0) =  (√ 2   sin 90 • , 0, − √ 2   cos 90 • ) =  √ 2   x .The parent proton with the dashed trajectory (visible only in panel (a)) is propagated through uniform fields only with an initial position of  ( =  0) = (−10R C , 3.3R C , 0R C ) and an initial velocity equal to  ̇( = 0) =  √

Figure 15 .
Figure15.Maps of the energetic neutral atom flux through a concentric spherical detector around Callisto when the moon is at its maximum distance below the center of the Jovian plasma sheet.See Figure9for details regarding the organizational structure.The range of the colorscale is equivalent to that used in Figure9to facilitate direct comparison between Europa and Callisto.
Liuzzo et al. (2022))ain with side length 30R C , centered at the moon.At Europa, we apply a cuboid-shaped domain with extensions of −8R E ≤ x ≤ 22R E in the x dimension, −10R E ≤ y ≤ 10R E in the y dimension, and −30R E ≤ z ≤ 30R E in the z dimension analogous toAddison et al. (2021).At both moons we apply a hierarchically-gridded mesh to adequately resolve, for example, the atmospheric scale height.Specifically, at Europa we use a grid resolution of 0.02R E (31 km) for |x|, |y|, |z| ≤ 1.5R E , 0.04R E for 1.5R E < |x|, |y|, |z| ≤ 3R E , and 0.08R E for |x|, |y|, |z| > 3R E .From the model ofLiuzzo et al. (2022), we adopt a grid resolution of 0.05R C (120 km), 0.1R C , and 0.2R C in cubes of the same size (but scaled in units of R C instead of R E ) at Callisto.
Liuzzo et al. (2022)021)dison et al. (2021),and Harris et al. (2021).It is thus defined with respect to ψ, the angle between a given position vector   = (  ) and the  − x axis.This angle is given by cos ψ = −x/|r|.At both moons, the neutral envelope is represented by a barometric law for the density as a function of altitude h and a factor to prescribe the strength of the asymmetry as a function of ψ.The hybrid simulation results adapted fromLiuzzo et al. (2022)at Callisto utilize a neutral density profile (for both O 2 and CO 2 ) with scale heights h i,C (for species i), and surface density n i,0 of the form

Europa With a Localized H 2 O Atmosphere Figure
(Roth, 2021)tes the role that a localized H 2 O atmospheric component(Roth, 2021), in addition to the primary molecular oxygen component, plays in influencing ENA emissions at Europa in draped fields when at the center of the Jovian plasma sheet.Both cases use the same, isotropic PAD to initialize parent protons on the starting grid.Displayed are the detectable ENA flux without the atmospheric H 2 O component (left column), with the H 2 O component included (center column), and the difference between the two maps at each energy (right column).Blue regions in the right column indicate where the detectable ENA emissions are elevated by the additional H 2 O profile.The inclusion of the H 2 O profile above the ramside apex results in a quantitative enhancement to the ENA emissions deposited in the low-latitude portion of the detector's ramside hemisphere, amplifying the bulge-like feature in the ENA emissions already present for a pure O 2 atmosphere.For example, at an initial parent proton energy of E = 40 keV, the ENA flux maximum without the H 2 O profile included is