The origin and fate of O$_2$ in Europa's ice: an atmospheric perspective

The early prediction and subsequent detection of an O$_2$ atmosphere on Europa, coupled with the discovery that Europa has an ocean under its ice mantle, has made this moon a prime astrobiologic target, soon to be visited by the JUICE and Europa Clipper spacecraft. In spite of the considerable number of observational, modeling, and laboratory efforts, understanding the physics leading to the observed morphology of Europa's near surface O$_2$ atmosphere has been problematic. This is the case as the observed emissions depend on the local incident plasma ion flux, the local temperature and composition of the regolith, as well as on the near surface electron temperature and density. Here we rely heavily on earlier reviews briefly summarizing the observational, laboratory and simulation efforts. Although it is agreed that radiolysis of the surface ice by the incident Jovian plasma is the ultimate source of observed O$_2$, a recent, simple model of thermal desorption from a regolith permeated with O$_2$ has changed the usual paradigm. This suggests that the observed orbital dependence of the local source of the near-surface O$_2$ atmosphere is due to thermal release of O$_2$ likely trapped on the ice grains at dangling bonds with a smaller contribution due to direct sputtering. This could also impact our understanding of the suggestion that the radiolytic products in Europa's regolith might be a source of oxidants for its underground ocean.


Introduction
Based on laboratory data, Jupiter's moon Europa was predicted to have a tenuous O 2 atmosphere due to the radiolytic decomposition of its icy surface by the incident Jovian plasma particles Johnson 1990). Although this prediction has been borne out by extensive Hubble Space Telescope (HST) observations (e.g., Hall et al. 1995; affect our ability to model, for instance, how strongly bound, on average, an individual O 2 molecule is in the radiation-damaged ice. However, consistent with laboratory data, O 2 and H 2 O 2 are observed as trapped species in the ice grains on the icy Jovian satellites (Spencer et al. 1995;Carlson et al. 1999a;Spencer & Calvin 2002;Hand and Brown 2013). Such observations are consistent with radiation damage studies in which under irradiation volatiles, such as O 2 , can become trapped in radiation produced 'bubbles' (multiple vacancy sites/voids) in an ice grain (Johnson & Jesser 1997) or at defect sites  with the concomitant H 2 O 2 being bound as more refractory molecules in an irradiated ice (e.g., Carlson et al. 1999a;Carlson et al. 2009). Although a large amount of important work on these topics still needs to be carried out, below we proceed based on our present, but incomplete, understanding.

Observations
By observing atomic oxygen emission lines using HST, Hall et al (1995; confirmed that Europa has an ambient gas consisting of both atomic and molecular oxygen. The column density of O 2 estimated from these observations was very roughly consistent with the amount produced as predicted by combining estimated plasma fluxes with the laboratory data for O 2 yields . These early observations have been extensively expanded and are roughly consistent with the more recent observations (e.g., McGrath et al. 2009;Roth et al. 2016) as well as those of the New Horizon spacecraft (Retherford et al. 2007; see also Plainaki et al 2018 Fig. 4). The concomitant ambient H 2 may have also been observed indirectly (e.g., Mauk et al. 2003) as well as the dissociation of ejected gas phase H 2 O molecules resulting in Lyman alpha emission from the excited atomic hydrogen. The latter was the first evidence that outgassing of H 2 O might be occurring at Europa (Roth 2014) and has been used recently to describe Europa's hydrogen corona (Roth et al. 2017). Here we focus on O 2 observations. The global source and loss rates of oxygen have also been modeled and discussed with the oxygen atmosphere exhibiting a near surface component dominated by O 2 and an extended component eventually dominated by O. Unfortunately the near surface plasma properties are still uncertain, affecting the calculation of both the radiolytic production of oxygen and its emission and loss rates. Therefore, estimates of the global O 2 source rate at Europa have varied due to the uncertainties of the nature of the particle flux as well as the actual physical conditions at Europa's surface as indicated in Fig. 1. Figure 1. Estimates of global source O 2 source rates with the colors indicating the incident radiation type assumed (Cassidy 2016) What is possibly more interesting, and not yet well understood, is the observed spatial morphology of the observed aurora, which is used to roughly extract the spatial distribution of the atmospheric molecular and atomic oxygen. Experiments indicate that the radiolytically produced O 2 leaves the surface with a speed distribution determined in part by its binding in ice and in part by the surface temperature (Johnson et al. 1983). For this reason it is not easy or maybe even possible to observationally distinguish between thermally desorbed and sputtered O 2 . In addition, the distribution of the O 2 ejection speeds of O 2 is found to be much less than its escape speed. This results in an average hop distance of an ejected O 2 that is much smaller than Europa's radius. In this way the observed surface-bounded atmosphere should reflect, to some degree, the local surface and plasma properties. Since Europa's surface composition appears to vary with latitude and longitude, as does the incident plasma temperature and flux (e.g., Nordheim et al. 2018), it is unsurprising that the auroral data suggest that its oxygen atmosphere is non-uniform (Roth et al. 2016).
Variations in the surface reflectance spectra are suggested to be due to trace sulfur species (Carlson et al. 1999b; on the trailing hemisphere, whereas in the Polar Regions and leading hemisphere the ice appears to be cleaner and possibly dominated by trace carbon species (e.g., Carlson et al. 2009;. Indeed, observations suggest that there may be more atmospheric O 2 over the regions of cleaner ice (McGrath et al. 2004;Roth et al. 2016) but the observed north-south differences in the emission rate are suggested to be due primarily to the differences in electron temperature and density (e.g., Roth et al. 2016, Fig. 10). Because of the uncertainties in the effect of the spatial variations in the composition and the incident plasma, both of which affect the O 2 source rate, carrying out simulations of the atmosphere in order to obtain a detailed description of the emission observations has been problematic. This is not only due to the spatial distribution of radiolytic production and loss, but, as stated earlier, is also due to the uncertainty in the binding state of the O 2 formed in the ice, and its eventual fate after it is ejected and then returns to the regolith. Due to these complications, Oza et al (2018a) took a somewhat different approach to be discussed shortly.
The published HST/STIS images have a spatial resolution between ~71 km and ~95 km per pixel depending on Europa's distance, allowing one to separate emission features originating from the atmosphere within ∼500km from the surface and those feature originating from an extended corona as analyzed in great detail in Roth et al. (2016). The excitation rate depends, of course, on electron temperature, as seen nicely in Fig. 12 of Roth et al. (2016). For the expected electron temperatures near Europa, an intensity ratio of ~2 for the 1356A/1304A lines is indicative of electron impact dissociation of molecular oxygen. It was estimated that throughout its orbit 95-99% of the near-surface (r < 1.25r E ) emission is due to electron impact of O 2 , as opposed to electron impact of O (Roth et al. 2016). In addition to north-south asymmetries discussed, longitudinal asymmetries were consistently found at several orbital positions, f orb , from the HST dataset spanning from 1999-2015 (see Fig.1a). The asymmetries shown were extracted by comparing average emission intensity over a range of latitudes and longitudes, but had no clear plasma or atmosphere explanation (e.g., Roth et al. 2016). Such asymmetries were difficult to discern from the patterns of individual images analyzed over the decades, but orbital and spatial averaging mitigates the uncertainties. One such orbital dependence described in Roth et al. (2016) was subsequently examined in Oza et al. (2018a, b): the near-surface dusk-overdawn oxygen emission intensity ratio. It was used to test simulations of the near-surface O 2 column density.
Although there is considerable scatter in the data, it is seen in Fig. 2 that to first order, the dusk-over-dawn oxygen aurorae intensity ratios are not extremely sensitive to orbital position, f orb , with a nearly independent average ~ 1.6, allowing for the considerable uncertainties. In this rough estimate, the 1356A line is used, as the effect of direct solar resonance scattering by the ambient O is negligible. That these observations do exhibit considerable scattering about the average is not surprising considering the spatial and temporal variability discussed above as well as the difficulty of observing near eclipse. In fact, the ratios near eclipse could be ignored. Below we consider the implications of this result. Before proceeding we note that, although the scatter in the data in Fig. 2 is large, especially near eclipse, the near-surface dusk/dawn emission ratio might be slightly larger in the vicinity of f orb ~180 o . In this region of the orbit the dawn half-hemisphere is the darker, trailing hemisphere and the dusk half-hemisphere is the leading, brighter hemisphere. This would appear to be consistent with a somewhat enhanced ratio, possibly consistent with an enhanced source in the region of cleaner ice. In the following we ignore this possibility as more observations would be required to be definitive.

Simulations
Assuming that Europa's atmosphere is marginally collisional or even collisionless, there have been a number of applications of the laboratory data in simulations of its atmosphere (e.g., Smyth & Marconi 2006;Shematovich et al. 2005;Saur et al 1998;Cassidy et al. 2013;Luchetti et al. 2016;Oza et al. 2018b). Laboratory data indicate there is a surface temperature and irradiation time dependence to the formation and ejection of O 2 by UV photons (e.g., Westley et al. 1995), by plasma electrons (e.g., Orlando & Kimmel 1997) and by plasma ions (e.g., Brown et al. 1982;Johnson et al. 1983;Teolis et al. 2017). Therefore, it was pointed out that the thermal and irradiation history must be accounted for in determining the radiolytic source rate (e.g., Teolis et al. 2005;Cassidy et al. 2013).
Unfortunately, it is difficult to carry out laboratory experiments in a sample that accurately represents a porous regolith (Cassidy and Johnson 2005). That is, the vapor deposited laboratory samples have a porosity at the submicron molecular scale, whereas the regoliths exhibit porosities at the 100s of micron scale. In addition, the charged particles used in many of the experiments had penetration depths smaller than the typical regolith grain size, whereas some experiments use incident particles with penetration depths greater than the typical grain size by irradiating relatively thin samples, in which case the substrate might play a role. Therefore, the concept of 'the surface' can differ in these experiments from that understood to be present on Europa or any other icy body. In addition, the thermal plasma electrons and ions at Europa primarily bombard the trailing hemisphere while energetic ions with energies above ~a few hundred keV bombard the surface much more uniformly . Energetic electrons on the other hand, bombard the surface in highly non-uniform patterns, leading to the formation of 'bombardment lenses' centered on the leading and trailing hemispheres at low latitudes (Paranicas et al. 2001;Nordheim et al. 2018) as indicated in Fig. 3. As energy deposition by energetic electrons dominates at cm to m depths, this degree of non-uniformity has important implications for the radiolytic production of O 2 at depth. Convolved with the significant variation in charged particle energy deposition vs. depth and surface location, is the variability of surface materials and trace species at different locations (e.g. as discussed, the leading hemisphere and polar regions appear to contain different trace species than the trailing hemisphere). For these reasons also agreement on the use of the laboratory data in the simulations remains problematic (see also Cassidy 2016;Plainaki et al. 2018). In spite of such uncertainties, there have been a number of efforts applying models of the laboratory data in collisionless Monte Carlo simulations of Europa's atmosphere. In a recent attempt at simulating the observed oxygen aurora, Oza et al. (2018b) tested nominal source distributions, including, for the first time, the effect of Europa's rotation about Jupiter. They also considered the effect of the varied spatial composition of Europa's surface on the sputtered O 2 as well as on its fate when it returned to the surface. In these simulations, the initially desorbed and thermally reemitted O 2 were tracked as they were ejected and ballistically hopped across the surface subject to gravity and the centripetal and Coriolis forces. Although the range of possible source processes and the range of the surface effects could surely be expanded beyond those used, the conclusion was that Europa's rotation about Jupiter was critical in generating the dusk-over-dawn asymmetry. However, the fact that the observed dusk/dawn hemispherical ratios did not vary strongly with orbital position, as seen in Fig. 2, could not be reproduced. That is, the models used gave near surface emission ratios that were strongly dependent on the orbital position, differing at some orbital positions from the HST ratios by an order of magnitude. Because the spatial dependence of the loss rate used in those simulations was a simple global average, it was also pointed out that one could carefully design a spatially variable source and a spatially variable loss rate that could certainly force the simulations to produce dusk/dawn ratios more uniform with orbital longitude as in Fig. 2. Since such a simulation would be artificial, a different approach was tried.
Although the physical processes are complex, the result in Fig. 2a based on the observations suggests some simple aspect was being overlooked. In order to focus on the role of Europa's orbital position with respect to the solar illumination, the surface rotation and the preferential bombardment of the trailing hemisphere, the longitudinal variation in the latitudinally averaged O 2 column density was modeled. This was done in order to directly extract the role of the longitudinal, rotational and thermal dependence of the likely source rate. Therefore, a simple rate equation, balancing source and loss rates, was used to estimate the dependence on the local longitude, f, of the time-variable, latitudinally averaged O 2 column density, N(f, t), vs. Europa's orbital position, f orb (Oza et al. 2018a). As Europa rotates with angular speed, W, phase-locked to Jupiter, f orb is simply related to the location of solar illumination angle, f s , on Europa's surface. Because the average ballistic hop distance for an ejected O 2 is much smaller than Europa's radius and the lifetime of O 2 is shorter than the time for it to diffuse by hopping a significant distance across Europa's surface, a rate equation for the local O 2 column density, N, can be approximated as (1) Here f is the longitudinal position on the surface, F is the local O 2 flux from the surface and n is the local loss rate. This ignores the hopping across the surface described in Oza et al. (2018b), which could be included by a diffusive term. However, in this form Eq. 1 is readily integrated over time as Europa rotates, allowing one to test the steady state spatial and thermal dependence of O 2 that is suggested by the spatial and thermal dependence of the source rate as well as on the loss rate. Although diffusion across the surface, as well as the centripetal and Coriolis forces, are ignored in Eq. 1, the suggested longitudinal and thermally dependent source and loss rates typically used gave near surface O 2 column densities that also depended very strongly on orbital position, as was seen to be the case in the detailed simulations in Oza et al. (2018b). This suggested that Eq. 1 could be used to test a number of assumed source processes. Using a spatially and thermally dependent sputter source in Eq.1 suggested by data from laboratory experiments, it is seen in Fig. 4 that the dusk/dawn ratio versus orbital position disagrees considerably with the ratios in Fig. 2 extracted from the emission observations. This was the case even allowing for the delays in emission suggested by laboratory measurements (Teolis et al. 2005;). However, O 2 trapped in an ice regolith can in principle be released thermally, as discussed frequently (e.g., Loeffler et al. 2006;Bar-Nun et al. 2013;Laufer et al. 2017), and as can also seen experimentally on warming an irradiated ice sample (Teolis et al. 2005). Therefore, a regolith permeated with O 2 that can be released thermally could in principal dominate the O 2 sputter source. Using an average loss rate, n, Oza et al. (2018a) noted that a source rate depending only on the solar illumination angle, f s , could account for the weak dependence of the dawn/dusk emission ratio on the orbital position. That is, writing the surface desorption flux vs. f s in a clearly oversimplified form where F 0 is the maximum sublimation at noon, f s = p and goes to zero in eclipse, the steady-state, latitudinally-averaged column density only depends on f s : N(f) ~ (F 0 /2n)[1-(1+ b 2 ) -1 (cos(f s ) -b sin(f s ) )] with b = W/n. Using this estimate of the column density the calculated half-hemisphere averaged ratio is found to be independent of f orb . That is, on substituting Eq. 2a into Eq. 1, the half-hemisphere average of dusk/ dawn ratio, <R>, of the near surface, O 2 column densities is Remarkably this depends only on b, the ratio of the rotation rate to the average loss rate, indicating that rotation is critical for the dusk/dawn asymmetry at Europa as shown in the simulations in Oza et al. (2018b). That rotation plays a critical role in the surface bounded atmospheres was noted much earlier for the lunar atmosphere (Hodges & Johnson 1968) and Mercury's Na atmosphere (Cassidy et al. 2015;Leblancand Johnson 2010). It is seen that in this approximation <R> is not only constant but is greater than one. That is, for an O 2 source that is determined primarily by the solar flux, the peak column is shifted from near noon toward dusk simply due to Europa's rotation. The scale height of the observed O 2 atmosphere has been discussed extensively and may be affected on the trailing hemisphere by collisions with sputtered water molecules (Oza et al. in prep). Yet it is agreed that the O 2 molecules have a scale height much smaller than Europa's radius. Therefore, it is not unreasonable that the averaged near surface O 2 column densities are roughly proportional to the averaged near surface component of the observed O 2 emission. Although that is a significant assumption, in Fig. 4 we compare the ratio in Eq. 3 to the observational ratios from Fig. 2. Using the loss rates based on photo and electron impact (Turc et al. 2014, Table 2) gives an ionization rate, n ~ 3.1x10 -6 /s for an electron density of 70 cm -3 . Using that value as a globally averaged loss rate and W = 2.1x10 -5 /s at Europa, then b ~ 6.7. Substituting b in Eq. 3 gives a half-hemisphere average dusk/dawn ratio of <R> ~ 1.2. Although <R> in Eq. 3 is independent of f orb , this underestimates the average ratio, ~1.6, obtained from the HST data. However, the sense of this result is consistent with an emission ratio nearly independent of orbital position. Recent modeling of the plasma flow (Dols et al. 2016) suggests that charge transfer with the ionized component of the atmosphere dominates the near surface loss rate of O 2 . If that is the case, including all O 2 loss processes (e.g., Lucchetti et al. 2016), the net, globallyaveraged loss rate can be as large as n ~ 2.0 x10 -5 /s. This would result in b ~ 1 and <R> ~ 1.6. This is, in surprising, and possibly fortuitous, agreement with the average HST result in Fig. 2.   Fig. 5 (a) Orbital modulation of Europa's equatorial sub-observer temperature, due to its albedo and the Jovian shadow, in degrees K vs. orbital longitude (b) temperature map sub-observer longitude 270 o vs. planetographic longitude and latitude at an obliquity roughly consistent with the Roth et al. observational data. The duration of the surface T drop at Jupiter's obliquity is small, but enough to matter. Depending on the season the eclipse duration is ~ 3% of a Europa day. The regolith T in this model has a skin depth ~ 4cm (Spencer 1987) and ignores a solid state greenhouse effect. (From Oza et al. 2018b, Fig. 2) The source rate in Eq. 2 is, of course, overly simplified as there are leading and trailing differences in composition and albedo resulting in much more complicated temperature profiles as seen in Spencer et al. (1999) and in Fig. 5. Therefore, the form in Eq. 2 could be improved by using a sublimation rate that depended on the actual temperature. Assuming that the dominant source is indeed thermal desorption from a regolith permeated with adsorbed O 2 , the source flux above, F(f s ), could be approximated, as is often the case, as proportional to exp(-U s /kT) with U s an average activation/binding energy. Here we crudely approximate an equatorial surface temperature like that in Fig. 5a, averaged over a range of longitudes and latitudes due to thermal inertia, as kT = kT 0 (1 -e cos ϕ s ) with e = DT/T 0 and DT, the temperature excursion. Assuming DT is small relative to T 0 , the average equatorial temperature, we obtain a latitudinally-averaged source rate, F, that is very roughly proportional to [1 -(U s /kT 0 ) e cos ϕ s ]. This is seen to have a dependence like that in expression in Eq. 2a, justifying its use to approximate a thermal desorption source. In instead of Eq. 2a, one could write, where c can be used to fit the difference in the noon/midnight ratio. In Eq. 2b, Df s is included to account for the shift in the peak temperature from noon as seen in Fig. 5b. That is, due to the thermal inertia of the icy surface the peak in the surface temperature of Europa in the near equatorial region is seen to be shifted toward dusk. For a thermal source like that suggested here, the peak temperature is shifted towards dusk slightly increasing the averaged dusk/dawn ratio, <R>, which we estimate to be ~ 15%. In addition, a direct sputter source would add to any thermally-induced source. For instance, adding to Eq. 2a a global sputter source, F s , as suggested in Cassidy et al. (2013), <R> in Eq. 3 becomes where d s = F s /F 0 . Not surprisingly, this expression for the dusk/dawn ratio is also independent of orbital position, but is altered by d s , the ratio of the direct sputter source to the peak-thermal source. Obviously, more realistic expressions for F(f s ) are needed, but only after we have a good description of the physics that is occurring. The above discussion simply shows that even with additional complexity a primarily thermally driven source is consistent with a critical aspect of the near-surface observations. There are, of course, uncertainties in the average loss rate, n, as well as the sputter to thermal source ratio, d s . Therefore, using the available observations to make a definitive statement on their relative importance is not likely at present. However, most previous simulations have assumed that d s is strongly dependent on orbital longitude. Therefore, a source rate that responds primarily to the local temperature as in Oza et al. (2018a) and not, primarily to an enhanced trailing hemisphere sputtering rate, appears to be a much simpler way to account for the lack of significant variation with Europa's orbital position of the averaged dusk/dawn emission ratio observed. Such a source would appear to suggest that Europa's regolith has become permeated with trapped O 2 , which upon thermal desorption populates the near-surface atmosphere, with, of course, a smaller contribution from the ultimate plasma source which replaces and redistributes the O 2 lost by the escape, subduction, dissociation and ionization.
Clearly, additional improvements can be made by accounting for the differences in temperature and composition on the leading and trailing hemispheres as discussed. However, as seen from the data in Fig. 2a the fully illuminated trailing hemisphere (orbital longitude 270 o ) and the fully illuminated leading hemisphere (orbital longitude 90 o ) exhibit surprisingly similar radially-averaged dusk/dawn emission ratios, though with significant uncertainties. This is consistent with the fact that the difference in radial dependence of the emission intensity of the 1356A line at these orbital longitudes is not large, as seen in Fig.  5 (red lines: dashed leading and dotted trailing). Even the off-disk intensity averages in Fig. 6 are surprisingly similar for the 1304A line. The similar radial dependences for the illuminated leading and trailing hemispheres both of which increase with radial distance from the subsolar point is not inconsistent with a predominantly thermally driven source. The fact that the hemispherical differences seen in Fig. 6 for the 1304 A line are larger than those seen in the 1356A suggests that spatial variation in the plasma, affecting the dissociation of O 2 , and the subsequent ionization and excitation of the dissociation product O, does play a role. Therefore, corrections to the average loss rate, n, in Oza et al. (2018a,b), and also used to obtain Eqs. 3 and 4 above should be accounted for. For this reason, a good description of the plasma temperature and particle flux very close to Europa's surface is needed. Although one could, in principle, reproduce the HST results in Figs. 4 and 6 by designing functions for the source and loss rates that depended on both f and f s , below we discuss the likelihood that a sublimation-dominated source is reasonable.

Discussion
As pointed out following the early laboratory measurements of the radiolysis of ice (Brown et al. 1982;Johnson et al. 1983), the preferential loss of the much more volatile H 2 , implies that the irradiated ice becomes oxygen rich. At Europa this effect is enhanced by the fact that H 2 , produced and desorbed during radiolysis, also readily escapes Europa's gravitational field, whereas the ejected O 2 primarily returns to the surface, as does a significant fraction of the O formed following dissociation of O 2 (e.g., Oza et al. 2018b). Therefore, it is clear that the regolith is oxygen rich, at least down to the penetration depth of the incident energetic particle radiation, ~ a meter (e.g., Cooper et al. 2001;Paranicas et al. 2003;Norheim et al. 2018). Because the regolith is porous, diffusion would suggest that over long time periods the ice could be oxygen rich to much larger depths (e.g., Johnson et al. 2003).
Although the state of the oxygen rich ice is not well defined, and likely differs between hemispheres and with depth, O 2 has been observed as trapped in Europa's surface ice grains (Spencer and Calvin 2002). That observation was associated with a band from a solid-like g phase O 2 at 5773Å. Because solid O 2 cannot be sustained in vapor pressure equilibrium at Europa's surface temperatures (e.g., Fray and Schmitt 2009), it was suggested to be trapped in radiation produced voids (Johnson & Jesser 1997). Recently, that band was observed to exhibit an unexplained spatial and temporal variability at Europa as shown (Spencer & Grundy 2017). This is likely related to a combination of plasma irradiation and local temperature affecting trapping (e.g., Shi et al. 2009) and, possibly, water vapor deposition from the plumes.
The size of the grains, which appears to decrease by more than a factor of 10 from the trailing and to the leading hemisphere (see Cassidy et al. 2013, Fig. 10), can influence the O 2 /H 2 O ratio inferred from the 5773Å O 2 observations. Assuming the band depth is due to O 2 -g, based on the expression derived by Hand et al. (2006) and the data in Spencer and Grundy (2016), then the trapped O 2 /H 2 O ratio seems, quite surprisingly, to be within a factor of a few of that emitted from Comet 67P/ Churyumov-Gerasimenko: O 2 /H 2 O (Bieler et al. 2015) and Comet 1P/Halley (Rubin et al. 2015). In addition to the O 2 observed at Europa in the 5773Å band, individual O 2 can be trapped at defect and dangling bond sites produced by the irradiation (e.g., Cassidy & Johnson 2005;Loeffler et al. 2006;Johnson 2011;Johnson et al. 2013) including at grain surfaces (Teolis et al. 2005). Indeed laboratory experiments indicate that the surface properties play a role in controlling the production and subsequent ejection of O 2 (e.g. Teolis et al. 2006;Teolis et al. 2017). The O 2 is seen to trap in the near surface region of the irradiated laboratory ice samples, which on Europa would be the surface of an ice grain exposed to the incident plasma. This O 2 is driven into the gas phase with the amount and rate depending on the irradiation history (Teolis et al. 2005). Indeed on irradiation followed by warming, O 2 is also seen to evolve (e.g., Teolis et al. 2009;2017) and warming reduces dangling bond sites . Therefore, the production and release of O 2 is clearly temperature dependent, although for the ~135K maximum temperatures expected at Europa, outgassing in laboratory samples in which H 2 O and O 2 are co-deposited was estimated to be small (Teolis et al. 2005). Emission of O 2 from non-irradiated samples of co-deposited O 2 and H 2 O has been related to the crystalline phase transitions (e.g., Laufer et al. 2017), which occurs due to the rearrangement of hydrogen bonds.
Although the formation and ejection of O 2 from a laboratory water ice sample due to the incident radiation is initially negligible, it grows with irradiation time (e.g., Brown et al. 1982;Wesley et al. 1995;Teolis et al. 2017). This is unlike what occurs in gas phase H 2 O and is due to the ability of the solid to accumulate trapped precursors (e.g., Johnson 2011). Therefore, a significant amount of solid-state chemistry occurs before the O 2 yield reaches steady state and precursor trapping has been shown to be critical in modeling of the measured thermal and fluence dependences using chemical rate equations (e.g., Johnson et al. 2005;Johnson 2011). Since the nature of the trapping of both the chemical precursors and the O 2 is critical, it depends on the availability of defects/dangling H bonds, and has been shown, not surprisingly, to be enhanced by irradiation .

O 2 Desorption Energy
The laboratory measurements of the steady state yield of O 2 exhibit a roughly exponential temperature dependence with an apparent activation energy ~0.06eV found in a number of experiments (Reiman et al. 1984;Brown et al. 1984;Fama et al. 2008;Teolis et al. 2017). This activation energy, smaller than the typical hydrogen bond energy in ice and even smaller than the sublimation energy for solid O 2 , is not yet explained. Since it is extracted from irradiation measurements, it is likely related to radiation-enhanced activation and/or diffusion.
Comparing the convenient, but overly simplified, model in Eq. 2a to the form for F(f s ) discussed prior to Eq. 2b, would imply that (U s /kT 0 ) e ~ 1 where e = DT/T 0 . Using DT ~ 10K and T 0 ~ 125K estimated from Fig. 5, a very approximate activation/bond energy can be obtained, U s ~ 0.14eV. This is similar to hydrogen bond energies in ice, which, of course, vary somewhat depending on the ice structure. This is roughly 1.5 times as large as the sublimation energy for pure O2. In the temperature range of interest on Europa, the outgassing from the co-deposited H 2 O/O 2 was slow, as discussed, until the temperature was increased to the phase transition ~ 155K. However, over the relevant temperature range for Europa the co-deposited samples were shown to steadily release O 2 (Fig. 4 in Laufer et al 2017). A desorption energy, possibly fortuitously, close to 0.14 eV was extracted from co-deposition experiments in support of the ROSINA measurements (e.g., Laufer et al. 2017; Table I). These crude comparisons suggest that thermal desorption of O 2 attached to dangling H bonds might be the dominant source of the observed near surface O 2 at Europa.
The early measurements of the sputtered O 2 speed distribution (Johnson et al. 1983;Reimann et al. 1984) indicted that most of the ejected O 2 returned to Europa's surface as mentioned earlier. Therefore, in their extensive simulations Oza et al. (2018b) were concerned with the re-condensing O 2 . They showed the heat of adsorption on the recondensing O 2 molecules significantly affected their rate of desorption back into Europa's atmosphere. Because of surface compositional differences, the effective binding can vary, affecting the spatial morphology of local O 2 column density (e.g., Fig. 4 in Oza et al. 2018b). Therefore, new simulations need to be carried out based on the thermal desorption concept suggested in Oza et al. (2018a) in which one treats the nature of the trapping and the temperature vs. depth and time more carefully. That is, based on the above discussion, the local source is likely dominated by thermal desorption of O 2 trapped primarily at dangling hydrogen bonds from a regolith permeated with O 2 and with a smaller contribution from direct sputtering. Figure 7: Schematic diagram of O 2 trapping and thermal desorption: 1) Primary origin of O 2 (and H 2 ) is magnetospheric ion radiolysis. 2) Due to preferential loss of H 2 , the regolith becomes oxygen rich enhancing the production of O 2. Formed and returning O2 can become trapped at incomplete (dangling) H bonds (shown) as well as in voids (as shown and observed by Spencer & Calvin 2002). 3) The accumulated O 2 can then be thermally desorbed from the weak dangling bonds due to solar heating, maintaining a quasivapor pressure equilibrium (Oza et al. 2018a), with a smaller gas-phase contribution from direct sputtering of O 2. A fraction of the trapped O 2. is likely subducted.

Summary
As pointed out, the state of the trapped O 2 in Europa's regolith, not only affects the observation of the gas-phase O 2 , but also affects its ability to diffuse downward and permeate the regolith. This is determined by the temperature variation with depth. That is, the local temperature gradient determines how any produced or returning O 2 diffuses through the ice matrix via defect or grain boundary diffusion (Johnson & Quickenden 1997;Johnson 2011). The presence of oxygen rich material at depth is also affected by the radiation which modifies the, top ~ meter of Europa's ice mantle as discussed above. This radiation not only produces an oxygen rich ice, but also can affect the crystal structure producing dangling bonds and causes radiation-enhanced diffusion. Since the porosity likely decreases with depth (e.g., Pappalardo et al. 2009), it is critical that the modeling needs to be carried out on the bulk processes and not just the near surface processes.
Although direct diffusion to the depth of the ocean is likely problematic, geologic mixing and subduction of oxygen rich ice has been suggested as a possible source of oxidants for putative ocean biology (Johnson et al. 2003;Greenberg 2008;Russell et al. 2017). These are also processes which might bring any ocean organics to the surface where they might be detected (Johnson and Sundqvist 2008). Any oxidants transported to the ocean would not only be O 2 , but other related and stable radiation-induced oxidants such as H 2 O 2 , sulfates, carbonates, carbonic acids, etc. which can yield an ocean with high chemical potential (e.g., Hand et al. 2006;Vance et al. 2016).
The simple paradigm for production of the observed O 2 atmosphere described here needs to be improved upon to assist in the planning of the observations of the JUICE and Europa Clipper missions. However, the concept considered suggests a fact that should not be surprising. That is, due to the considerable radiolytic processing of Europa's surface, and Europa's significant gravity, it is not surprising that Europa's ice mantle is permeated with oxidants and trapped O 2 . This was, of course, the prime motivation for the suggestion that radiolytic processing of the surface ice combined with geologic activity might lead to subduction of oxidants to Europa's ocean (Chyba 2000;Johnson et al. 2003;Greenberg 2008;Russell et al. 2007).. If Europa's regolith is permeated with trapped O 2 , this might at first appear inconsistent with the Cassini observations at Dione and Rhea. Teolis et al. (2016) showed that the effective source of O 2 was much smaller than expected based on combining laboratory data with Cassini spacecraft measurements of the plasma flux (Johnson et al. 2008). However, the plasma fluxes at these objects are orders of magnitude smaller than at Europa, these bodies are on average far colder, and the surfaces appear to be older indicating less geologic activity than observed, for instance, at Enceladus (Porco et al 2006) and Europa (Roth et al. 2014;Sparks et al. 2016;2017). This suggests that the application of the interesting and useful laboratory data on the gas-phase production of species produced by radiolysis must be applied carefully at each body, accounting for the local environment and the surface properties -----a not very surprising conclusion. Therefore, we suggest new simulations need to be carried out based on the premise that Europa's near surface regolith, consisting of ice grains sizes that vary with depth, and with near surface thermal processing that varies with time and depth, is permeated with trapped O 2 that is primarily formed by the incident Jovian plasma. That is, molecules ejected by sputtering primarily return to the surface where they can diffuse to binding sites of varying strength and react with trace species. Subsequent radiation can drive the trapped O 2 back into the gas phase, but, as suggested in Oza et al. (2018a), thermal desorption of weakly trapped O 2 , likely at dangling H bonds, appears to dominate the observed, near surface O 2 atmosphere.