© The Authors 2026

1 Introduction

Jupiter’s moon Europa possesses a tenuous atmosphere, which is constantly replenished from the weathering of the icy surface (McGrath et al. 2009; Johnson et al. 2009). Irradiation by ions from Jupiter’s magnetospheric plasma and radiolysis of the surface water ice is considered to be the main source for the molecular oxygen (O₂) and hydrogen (H₂) in the atmosphere (e.g. Johnson et al. 2009). However, we note that recent studies have suggested electron-induced radiolysis might also be an important source for icy moon atmospheres (Galli et al. 2018; Meier & Loeffler 2020). Overall, O₂ is the globally dominant species near the surface and forms a near-surface layer with an expected scale height of roughly 20-50 km (Shematovich et al. 2005). The lighter H₂ has a much larger scale height of hundreds of kilometers and tends to escape more readily (Smyth & Marconi 2006). Surface ice sublimation might serve as a source for H₂O molecules in the atmosphere (Plainaki et al. 2018), but significant fluxes would be expected only in the subsolar region. Water-group species in the atmosphere such as OH, O, and H can be produced by dissociation of the primary molecular atmosphere, but also directly through radiolysis of the ice. The atmospheric density likely varies with local time, hemisphere, and magnetospheric conditions.

Hubble Space Telescope (HST) observations of far-UV atomic oxygen emissions at 1356 Å and 1304 Å have provided the first detection of Europa’s molecular atmosphere and set new constraints on its characterization (Hall et al. 1995; Hall et al. 1998). The relative intensities of the O I 1304 Å and O I 1356 Å emission provide a means to distinguish between different oxygen-bearing atmospheric species (Roth et al. 2021).

Using HST’s Space Telescope Imaging Spectrograph (STIS), these UV emissions can be spectrally separated and spatially resolved (Roth et al. 2016).

Near Europa’s limb, the observed oxygen emission ratios (i.e., the ratio of the OI1356 Å to the OI1304 Å emissions) are consistent with an O₂-dominated atmosphere. In contrast, in the subsolar region on the trailing hemisphere, the lower OI1356 Å/OI1304 Å ratio suggests that H₂O is the dominant species (Roth 2021). Observations of Europa’s optical aurora are likewise consistent with an O₂-dominated atmosphere and provide weak evidence of the presence of H₂O as well (de Kleer & Brown 2018; de Kleer et al. 2023). However, optical observations are limited to disk-integrated quantities during the cooler eclipse phase. In contrast, spatially resolved STIS UV measurements are obtained in sunlight, when warm ice sublimates more readily.

Further indirect constraints on the atmosphere were recently provided by measurements of H⁺₂, O⁺₂, and O⁺ pickup ions near Europa as well as H⁺₂ detected farther away from the moon by the Juno spacecraft (Szalay et al. 2022, 2024). These measurements provide first evidence of the H₂ component and suggest that the turnover (production and loss) of both H₂ and O₂ in the atmosphere is lower than assumed in prior studies (e.g., Smyth & Marconi 2006).

While molecular hydrogen has not been detected yet with remote sensing observations, a widely extended atomic hydrogen exosphere was measured through attenuation at the Lyα line (1216 Å) in HST/STIS observations of Europa in transit of Jupiter (Roth et al. 2017b). The first estimations from Roth et al. (2017b) suggested that this H exosphere is more tenuous than those that of Ganymede (Barth et al. 1997b; Feldman et al. 2000; Alday et al. 2017; Naesenius et al. 2026) and Callisto (Barth et al. 1997a; Roth et al. 2017a), as the latter were detected through Lyα emissions from resonantly scattered sunlight. However, a comparison of Europa’s transit data to later HST observations of Ganymede in transit (Roth et al. 2023) indicate that Europa’s H exosphere densities were underestimated by Roth et al. (2017b) when using the transit attenuation method. Specifically, Europa’s H abundance should be about five times higher and, thus, also detectable in emission (not only absorption) in HST Lyα data. In the initial study by Roth et al. (2017b), the line-center cross-section was assumed to be the effective cross-section for the attenuation; however, a lower, line-integrated cross-section should have been used for more accurate results.

The spectral UV images of Europa taken by HST/STIS (McGrath et al. 2009; Roth et al. 2016) map emissions at the two oxygen multiplets around 1356 Å and 1304 Å as well as the hydrogen Lyα emission (1216 Å). In one particular observation taken in December 2012, a localized emission surplus was inferred to be consistent with electron-impact dissociative of H₂O as the source (Roth et al. 2014b). Such a localized abundance of H₂O requires a local source because of the short lifetime of water vapor in Europa’s atmosphere, which is the time when an H₂O molecule is aloft before returning to and presumably condensing on the icy surface (roughly some tens of minutes). The detection was interpreted as outgassing into a plume. An apparent surplus of OI1304 Å emissions from the same location, but no enhanced emissions at OI1356 Å in the same observations, which was consistent with a signal from electron-impact dissociative excitation of H₂O, as measured in laboratory studies (Makarov et al. 2004). The OI1304 Å and OI1356 Å emissions elsewhere on Europa during the December 2012 observations could be aptly explained by the global O₂ atmosphere (Roth et al. 2014b; Roth 2021).

A direct follow-up campaign with HST/STIS to test a hypothesized connection of plume activity to periodic tidal stresses did not provide any confirmation (Roth et al. 2014a). Although various studies later claimed further evidence of plumes at different locations on Europa (e.g., Sparks et al. 2016; Jia et al. 2018; Paganini et al. 2020), there has been no confirmed detection or consistent picture of plume activity at Europa (Giono et al. 2020; Roth et al. 2025).

We present an analysis of 23 Lyα (1216 Å) images of Europa extracted from HST/STIS data taken in the spectral G140L setup in 1999 and on various dates between 2012 and 2020. The simultaneously observed oxygen emissions were published in Roth et al. (2016) and Roth (2021) for the datasets taken until 2015. The Lyα emissions from the first five datasets were included in the analyses in Roth et al. (2014b) and Roth et al. (2014a) and are reanalyzed here. The 23 datasets presented were all taken when Europa was sunlit (not eclipsed) and not in transit of Jupiter. We focus on two different aspects in the Lyα signal: (1) resonant scattering emissions from Europa’s H exosphere, which are expected to be globally abundant and persistent; and (2) potential emissions from water vapor plumes, which would be localized and likely transient.

2 Observations and data analysis

We analyzed all observations of Europa in sunlight (and not transiting Jupiter) taken by HST/STIS with the G140L grating and the 52″ × 2″ slit. Table 1 and Figure 1(left) provide an overview of the observing parameters and geometries of the 23 analyzed datasets. During each visit, between two and ten exposures were taken by STIS with varying exposure time, within two to five HST orbits. In the visits after 2014, the capability of HST to observe moving targets over several hours (or several orbits of HST around Earth) was limited. As a result, most visits after 2014 are shorter with fewer exposures. The water vapor signal derived in Roth et al. (2014b) was observed during visit 3 in December 2012 when Europa was near eastern elongation. Orbital positions similar to visit 3 were thus preferred in the follow-up HST campaigns, leading to a clustering of visits near maximum eastern elongation (Figure 1).

We combined all exposures taken during a visit into one observation to maximize the signal. Both the observing geometry and plasma environment change over the course of a visit, as indicated by the ranges in sub-observer W longitude and Jupiter’s System-III longitude at Europa in Table 1. However, the gain in signal by combining the exposures is needed and any possible variability is neglected as done in previous studies (Roth et al. 2014b). Figure 1 (right) shows an example detector image of the superposition of the three exposures taken during visit 22 in 2018.

The Lyα geocorona (i.e., emissions from Earth’s extended H exosphere) fill the entire slit (yellow in Figure 1 right). The geocorona is particularly bright when HST is on the dayside of the Earth. As in previous studies (Roth et al. 2014a), we selected exposures with low geocorona signal or partially cut exposure durations to remove high-geocorona periods. For this purpose, the background brightness along the Lyα slit away from Europa was monitored using the STIS time-tag data and only periods where the measured background is lower than 10 kR (kiloRayleigh) were used. This threshold has been applied to all datasets except for visit 3. To allow for a direct comparison to the results in Roth et al. (2014b) we use the full exposure time as used in that study. The background was however generally low during visit 3, except for the first 8 minutes in three out of the nine exposures in total (24 minutes high geocorona vs. 140 minutes low geocorona). The total exposure times of the data used in the analysis are given for each visit in Table 1.

In addition to the geocoronal background, the slit centered on the Lyα wavelength contains contributions from interplanetary H, sunlight reflected from Europa’s surface, and Europa’s own H exosphere. Localized emissions may also be present near Europa, produced by electron excitation of H-bearing molecules, such as the proposed water vapor aurora associated with a plume (Roth et al. 2014b). Europa’s exact position on the detector is determined from the reflected sunlight signal, following the approach used in previous STIS studies. To disentangle these various Lyα components and accurately determine Europa’s detector position, we constructed a model image (see Section 3) that includes all known Lyα sources except potential localized emissions (e.g., plume aurora). The model image was also used to find the positions of Europa’s disk on the detector.

The model includes several fitting parameters that primarily scale the relative contributions of the different Lyα sources. In addition, various central (x, y) pixel positions for Europa’s disk are tested. Because Europa’s actual position on the STIS detector differs from that specified in the data file header, it must be determined directly from the observed signal. Accurate localization is particularly important when searching for potential localized emissions near the disk’s limb, where a steep gradient exists between on-disk and off-disk signals (Giono et al. 2020). To further constrain the disk center, the surface reflection at longer wavelengths (λ=1430-1530 Å), where no atmospheric emissions are present, is used to independently assess the central y pixel. The x pixel cannot be constrained with the continuum trace (and only at Lyα). Given HST’s pointing stability of ∼0.005 arcsec (or 1/5 of a STIS pixel), the position of the moon on the detector can be assumed to be stable over the course of a visit.

Fig. 1 Overview of the analyzed HST/STIS observations. Left: orbital positions of Europa from the start of the first exposures to end of the last exposures for each of the 23 HST/STIS visits. Note: the gaps between exposures of a visit are not shown. Visit 21 combined exposures before and after transit. Right: complete STIS detector spectral image from visit 22. The Lyα signal from the geocorona and the interplanetary hydrogen fills the complete slit (yellow-green region). Reflected continuum sunlight from Europa’s surface forms the trace around row y=200. Atmospheric oxygen emissions are present at 1304 Å and 1356 Å, exceeding the surface reflection.

3 Model for the Lyα emissions in the STIS images

The constructed model consists for three components considering different sources of Lyα emission. The first component of the model represents the background and foreground emission present throughout the slit. This contribution is assumed to be constant along the x-direction (i.e., from the left to the right edge of the Lyα slit; yellow in Figure 1, right panel). Along the slit (y-direction), the intensity is modeled with a second-order polynomial function to capture variations in brightness. This component includes geocoronal foreground emission as well as emissions from the hydrogen in the interplanetary medium (IPH) both in front of and behind Europa. The interplanetary medium is the interstellar gas cloud passing through Solar System and the hydrogen component scatters the solar Lyα flux (Izmodenov et al. 2013).

Because Europa’s disk blocks the background, a constant value, representing the interplanetary H brightness derived from a model by Pryor et al. (2024), was later subtracted from the disk. The three coefficients of the polynomial along y were fitted (independently) for each dataset.

The second model component is the signal from Europa’s H exosphere. The measurements in attenuation against Jupiter (Roth et al. 2017b) showed that the exospheric H number density, n_H, is consistent with an escaping (comet-like) profile that falls off with $\frac{1}{{\rho }^2}$, where ρ is the radial distance to the moon center. When integrating this number density profile along the line-of-sight from a distant observer (infinity) to the moon surface (for pixels on the disk where r ≤ R_E) or to far beyond the target (to infinity, outside the disk r > R_E), the resulting line-of-sight H column density in the image plane around the moon becomes $$N_{\mathrm{H}}(r)=\{\begin{array}{ll}{n}_{\mathrm{H},0}\frac{R_E^2}{r}\arcsin \left(\frac{r}{R_E}\right)=N_{\mathrm{H},0}\frac{R_E}{r}\arcsin \left(\frac{r}{R_E}\right)& \text{if\hspace{0.17em}}r\le R_E,\\ n_{\mathrm{H},0}\frac{R_E^2}{r}\pi =N_{\mathrm{H},0}\frac{R_E}{r}\pi & \text{if\hspace{0.17em}}r>R_E,\end{array}$$(1)

with the radial distance to Europa’s center in the image plane, r, Europa’s radius, R_E = 1560.8 km, the peak number density at the surface, n_H,0, and the vertical column density, $$N_{\mathrm{H},0}=n_{\mathrm{H},0}{R}_E.$$(2)

This profile (expressed via Eq. (1)) assumes that the H atoms escape to infinity with constant velocity and thus neglects Europa’s (and Jupiter’s) gravity, as well as other processes such as radiation effects. Despite these simplifications, we chose to adopt this profile for its simplicity and reproducibility.

Europa’s H exosphere can further be assumed optically thin, with a line center optical depth of τ ≈ 0.3 at the line center for the vertical H column density and temperature derived here (Section 4.2). Thus, we assume the same radial decrease for the intensity as a function of distance, r, to the disk center in an image given as $$I_{\mathrm{H}}(r)=\{\begin{array}{ll}{I}_{\mathrm{H},0}\frac{R_E}{r}\arcsin \left(\frac{r}{R_E}\right)& \text{\hspace{0.17em}if\hspace{0.17em}}r\le R_E,\\ I_{\mathrm{H},0}\frac{R_E}{r}\pi & \text{\hspace{0.17em}if\hspace{0.17em}}r>R_E,\end{array}$$(3)

with the disk center brightness, I_H,0.

For the third model component, which is the contribution from sunlight reflected off Europa’s surface, we followed the approach of previous studies (e.g., Roth et al. 2014b; Roth 2021). A high-resolution solar spectrum from the SOlar Stellar Irradiance Comparison Experiment (SOLSTICE) on the Solar Radiation and Climate Experiment (SORCE) (McClintock et al. 2005) was mapped onto Europa’s disk, accounting for the spectral dispersion of the STIS configuration. Because the solar flux varies with solar longitude, we selected a SOLSTICE spectrum from a nearby date when the solar longitude visible from Earth corresponded to that illuminating Europa at the time of the observation (Roth et al. 2021).

The surface reflection model incorporates both lateral albedo variations, due to changes in surface composition and texture, and the spherical geometry of Europa. Because the Lyα reflectivity is roughly anti-correlated with the visible albedo (Hendrix et al. 2005; Roth et al. 2014b; Becker et al. 2018), we inverted the USGS visible albedo map of Europa (United States Geological Survey 2002), adjusted the contrast, and projected it onto the disk for the specific viewing geometry. The reflection from the rough, spherical surface was then simulated using the Oren & Nayar (1994) model, as applied to Europa in previous studies (Sparks et al. 2016; Giono et al. 2020; Roth 2021).

By jointly comparing all datasets, we found the best agreement between the model and STIS images for a slightly reduced albedo contrast, where the maximum reflection is about five times the minimum, and a surface roughness parameter of 0.8. This roughness value is somewhat higher than the 0.57 adopted in previous works (Sparks et al. 2016; Roth 2021), resulting in a disk that appears more uniformly bright, namely, it is closer to the uniform disk assumed in earlier studies (e.g., Roth et al. 2014b). The final reflected-sunlight model is scaled by the albedo, which was treated as a free parameter and fitted individually for each dataset. Although significant albedo variations between observations (especially those with similar geometry) were not expected, the albedo was fitted separately to ensure the correct scaling of the reflected sunlight component in each case.

The model images were first generated with a 3× higher resolution than the STIS data resolution (the latter is roughly 40 pixels across Europa), primarily to ensure a smooth transition at the limb from on-disk to off-disk in the resulting model images. When all the contributions were combined, the spatial resolution was adjusted to the STIS data and convolved with the STIS point-spread-function at the Lyα wavelength (Krist et al. 2011). All model components with their parameters (exosphere brightness, surface albedo, and background polynomial) were then fitted simultaneously for the best agreement with the observed 2D STIS image within the slit.

In the observations, Europa is located in the lower part of the slit (except in the 1999 data) to avoid regions on the detector affected by fiducial bars and dark noise (see Figure 1). For the fit, all pixels from the disk center row up to 200-300 pixels away from the disk center in both directions (and at least 5 pixels away from the detector edge) were included, with the chosen range depending on the disk position and noise level.

Figure 2 shows the data and the fitted model in averaged brightnesses in Rayleigh (1 R = $\frac{{10}^6}{4\pi }\frac{\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}\mathrm{s}}{{\mathrm{c}\mathrm{m}}^2\mathrm{s}\mathrm{s}\mathrm{r}}$) along the detector y and x axes for two example visits. There is a certain degeneracy when fitting the H exosphere profile and background simultaneously but we find very similar results for H density for different choices of y ranges for the fit, which lends confidence to the robustness.

We used the model to find a (x, y) position of Europa’s disk center. The best position was assumed to be the one that minimizes the sum over all pixels (within 1.4 R_E of the disk center) of the squared differences between observed and model pixel brightness. Following an independent check, we found that the fitted surface albedo is also highest for the best-fit position in almost all datasets, as expected for the best alignment. For most images, the sum of the squared residuals varies only marginally in a range of 2 pixels around the minimum, suggesting an uncertainty of ±2 pixels in both the x and y positions.

We then used the surface reflection continuum trace at longer wavelengths (λ = 1430-1530 Å; see Figure 1, right), where no atmospheric emissions are expected to be present, to provide an independent constraint on the y position. For this purpose, we fit a linear function plus a Gaussian to the brightness profile along y, integrated over the wavelength range (x), where the center of the Gaussian was adopted as the best-fit y position.

For most datasets, the y position derived from the reflected continuum agrees with the Lyα-based value to within one pixel. In cases where the two estimates differ by two to three pixels, we adjusted the Lyα y position by one pixel toward the continuum-derived value. For example, during visit 22, the Lyα fit suggested a central position at y = 205, while the continuum fit yielded y = 207; therefore, we adopted y = 206 for the subsequent analysis, as shown in Figure 2b.

Fig. 2 Data analysis and model fitting. Top: brightness for the entire Lyα slit along the y-axis, averaged over the x-axis (i.e., over the slit width). Middle: Zoom-in of the top panel. Bottom: brightness along the x-axis, averaged over a y range covering the disk diameter. The STIS data are shown as a black histogram. The fitted model is represented by dashed lines corresponding to different components: green for the background only, blue for the background plus H exosphere contribution, and red for the full model including surface reflection and IPH correction (see text). Note: the profiles are smoothed due to averaging over a detector axis (and, thus, across the disk), as well as by the convolution of the model image with the STIS point spread function. The vertical dotted lines indicate the disk center in all panels. The model was fitted within the central disk and the gray shaded regions (see top panel).

4 Hydrogen exosphere analysis

4.1 H exosphere brightness

The exosphere brightness peaks at the limb, where it is I_H(r = R_E) = I_H,0 π, and decreases to I_H,0 in the disk center (r = 0), as given by Equation (3). For further analysis and interpretation, we used the value for the brightness at disk center I_H(r = 0) = I_H,0 from the fitted exosphere brightness profile as reference. Figure 3a shows the retrieved disk center H exosphere brightness, I_H,0, for all 23 visits, ranging between 18 R (visit 29) and 139 R (visit 6). Notably, a zero brightness was fitted for visit 24.

Fig. 3 Analysis of Europa’s H exosphere signal: (a) H exosphere brightness in the disk center from the fitted exosphere profiles (Eq. (3)) for each visit. (b) Solar Lyα flux (black asterisks) measured at 1 AU near the observation date (see text). The inverse of the squared heliocentric distance of Europa (orange triangles, in unit 0.01 AU−2) indicates the scaling of the flux used for calculating g-factors. (c) g-factors for resonant scattering of the sunlight (blue diamonds), of the IPH illumination (green) and combined (black). (d) H exosphere vertical column densities from conversion of measured brightnesses. (e) Same as in (d) but plotted against the line-of-sight velocity of Europa with respect to the Earth at the time of the observation. (f and g) Same as (e) but only for visits 2-18 taken in 2012-2015 and for visits 21-29 taken in 2018-2020, respectively. Fitted attenuation profiles for different temperatures of the H in Europa’s exosphere, see text for more information. The 1999 visit is shown in grey in (f) and (g) for comparison.

4.2 Conversion from brightness to column density

We assume the H exosphere to be optically thin at the Lyα line and the observed emission to originate from a single resonant scattering process. The resonantly scattered sunlight from the H exosphere measured by STIS is then given by the photon scattering coefficient (or g-factor, gL_yα) times the column density $$I_{\mathrm{H}}(r)=g_{Ly\alpha }{N}_{\mathrm{H}}(r),$$(4)

connecting Equations (1) and (3). The exosphere model part is scaled linearly with one parameter in the individual fit for each dataset, represented by the vertical H column density, N_H,0, in Equation (1).

The strongest illumination source is the incoming solar flux, but we find that the illumination by the IPH can be an additional source and we chose to include it in the analysis. Thus, we can estimate the scatterable spectral flux and the resulting scattering factors for the two sources, g_sun and g_IPH. From these we calculate a combined resonant scattering factor as $$g_{Ly\alpha }=g_{sun}+g_{IPH}=(\pi F_{sun}+\pi F_{IPH})\frac{\pi e^2}{m_{e}c}f\frac{{\lambda }_0^2}{c},$$(5)

where πF_sun and πF_IPH values are the spectral fluxes from the Sun and IPH in the rest frame of Europa, e is the electron charge, m_e is the mass of the electron, c is the speed of light, f is the oscillator strength (0.416), and λ₀ is the Lyα wavelength (1216 Å).

For the solar source, we used the line-integrated composite Lyα photon flux at 1 AU (Machol et al. 2019) from the same date as for the surface reflection (considering the solar longitude variations), as shown in Figure 3b. The line-integrated flux is converted to spectral flux at the line center using the expression by Kretzschmar et al. (2018) (see their Equation (4)). We then considered a small offset from the line center due to the radial velocity component of Europa with respect to the Sun, using a Lorentzian-Gaussian solar line shape as used by Gladstone et al. (2015). Figure 4a shows three example relative velocities and the corresponding changes in the solar spectral flux. The spectral flux at the Doppler shifted wavelengths is slightly higher than the line center spectral flux by up to 6% (Figure 4a). Finally, the flux is adjusted to the Sun-Europa distance (Figure 3b and Table 1).

For the IPH, we simulate the all-sky Lyα brightness using the model of Pryor et al. (2024). The total flux from the dayside sky hemisphere (as seen from the sub-solar point on the surface) is typically near 4 × 10⁸ photons/cm²/s and hence about 2% of the total solar flux. We assume a simple Gaussian line profile and a temperature of 15000 K for the IPH (Izmodenov et al. 2013). The narrower IPH spectral line compared to the solar line leads to a relatively higher line center spectral flux of 3.5 × 10⁹ photons/cm²/s/Å or ∼23% of the solar line center spectral flux (see the maximum of the IPH curve and the value at the dip in the solar line center in Figure 4a). The interplanetary hydrogen moves with a velocity of 20-25 km/s in the direction of ecliptic longitude and latitude of 72.5° and −8.9° (and the IPH profile is thus shown at a slight offset in Figure 4a).

To estimate the Doppler shift of the IPH hydrogen spectral line resulting from Europa’s relative motion, the radial velocity across each point in the sky must be considered (Joshi et al. 2025). For each visit, we generated a sky map of the IPH brightness using the model of Pryor et al. (2024) (see Figure A.1) and computed the relative velocity for each sky direction at the time of observation to estimate the IPH line flux that can be scattered by Europa’s H exosphere. When the brightest region of IPH illumination in the sky (near the Sun) is strongly Doppler-shifted, the scatterable IPH flux is low, as seen during visit 3. Conversely, when the IPH motion is largely perpendicular to the Sun-Europa line, the Doppler shift in the brightest region is minimal, and the scatterable IPH flux is high, as in visit 29.

Figure 3c shows the resulting g-factors for both the Sun and IPH for each visit. The H column densities after conversion with the g-factors (Eq. (4)) are shown in Figure 3d.

Fig. 4 Analysis of the spectral aspects of the Lyα signals: (a) Typical Lyα spectral fluxes from the Sun (solid) and from the IPH (dash-dotted) at the heliocentric distance of Jupiter. The IPH flux is shown as a Gaussian profile with T = 15 000 K, at an offset to the solar line because of the IPH movement relative to the Sun. The scatterable flux by the H in Europa’s exosphere varies with the relative line-of-sight velocity between the moon and the source (i.e., the Sun or IPH in a certain direction, as shown for three example velocities). The vertical dashed lines (purple, grey, orange) show how the contributions from the Sun and the IPH can vary due to the respective Doppler shift. (b) Lyα transmission of the Earth H exosphere in the rest frame of Earth for an assumed column density of 1 × 1013 cm−2 and temperature of 1000 K (solid). Lyα emission profiles from Europa’s H exosphere (T = 1000 K) from resonant scattering before (dotted) and after (dashed) attenuation in the Earth exosphere from the shown transmission. Three example profiles (blue, grey, red) are shown for three different relative line-of-sight velocities of Europa with respect to Earth. The strongly shifted emission (blue) is not affected by the Earth at all.

4.3 Extinction in Earth exosphere and signal measured by HST

Extinction of emissions from planetary H exospheres can occur in Earth’s extended exosphere via scattering by atomic hydrogen (Alday et al. 2017). For simplicity and following Alday et al. (2017), we assumed that both the Europa exosphere emission profile and the profile of the attenuation by H at Earth are only Doppler-broadened due to the thermal motion of the H atoms. We further assume that the velocity profiles are described by Gaussian normal distributions. To convert the velocity profiles to spectral Lyα emission profile as a function of wavelength (or vice versa), the relative velocity, ∆v, is translated to wavelength shift as $\mathrm{\Delta }\lambda =\frac{\mathrm{\Delta }v}{c}{\lambda }_{Ly\alpha }$, with the speed of light, c, and the Lyα vacuum wavelength, λ_Lyα = 1215.667 Å. The velocity profile for the spectral brightness emitted from Europa’s exosphere is then given by $${\tilde{I}}_{\mathrm{H}}(v)=\frac{I_{\mathrm{H},0}}{2\pi k_B{T}_{\mathrm{H},Eur}/m_H}\exp (-\frac{(v-v_0{)}^2}{2k_B{T}_{\mathrm{H},Eur}/m_H})$$(6) $$=\frac{g_{Ly\alpha }{N}_{\mathrm{H},Eur}}{2\pi k_B{T}_{\mathrm{H},Eur}/m_H}\exp (-\frac{(v-v_0{)}^2}{2k_B{T}_{\mathrm{H},Eur}/m_H}),$$(7)

with the reference velocity, v₀, the H atomic mass, m_H, temperature of Europa’s H exosphere, T_H,Eur, the Boltzmann constant, k_B, and the line-integrated brightness, $I_{\mathrm{H},0}=g_{Ly\alpha }{N}_{\mathrm{H},Eur}$. We note that Europa’s atmosphere is a rarified gas environment and thus it is not necessarily distributed thermally, but the emission Doppler width can still inform about an effective temperature (e.g., Lierle et al. 2022). Figure 4b shows three examples of emission profiles (normalized to a maximum of 0.8 for visibility) for T_H,Eur = 1000 K and v₀ values of 0 km/s, +10 km/s and −20 km/s.

The intensity measured by HST after extinction in the Earth exosphere, I_HST, is given by the Beer-Lambert law I_HST(v) = I_H(v)exp(−τ(v)). The optical depth τ is the product of the line-of-sight Earth exosphere column density, N_H,Geo and the velocity dependent attenuation cross-section $${\sigma }_{scat}(v)={\sigma }_0\exp (-\frac{v^2}{2k_B{T}_{\mathrm{H},Geo}/m_H}),$$(8)

with the line-center scattering cross-section, σ₀. Here, we assume the rest-frame of Earth. The transmission for an assumed example column density of N_H,Geo = 10¹³ cm⁻² and temperature T_H,Geo = 1000 K of the Earth exosphere is shown by the solid black in Figure 4b. In the Earth rest frame, the velocity offset of the Europa emission profile (-₀ in Eq. (7)) is given by the line-of-sight component of the velocity of Europa with respect to Earth υ_rad. The attenuated spectral brightness from Europa as measured by HST is then given by $${\tilde{I}}_{HST}(v)=\frac{g_{Ly\alpha }{N}_{\mathrm{H},Eur}}{2\pi k_B{T}_{\mathrm{H},Eur}/m_H}\exp (-\frac{(v-v_{rad}{)}^2}{2k_B{T}_{\mathrm{H},Eur}/m_H})$$(9) $$\times \exp [-N_{\mathrm{H},Geo}{\sigma }_0\exp (-\frac{v^2}{2k_B{T}_{\mathrm{H},Geo}/m_H})].$$(10)

The total measured brightness, I_HST, for a given v_rad is given by integrating the emission profile over the velocity. The total intensity of the example emission profiles in Figure 4b are 31% (gray, no offset), 95% (red, 10 km/s redshift) and >99.9% (blue, 20 km/s blueshift) of the original intensity, still assuming N_H,Geo = 10¹³ cm⁻² and T_H,Geo = 1000 K (solid line, Figure 4b).

For the analysis of the Europa exosphere intensity measured in the HST datasets, we fit the total measured intensity (i.e., Equation (9) integrated over v) to all datasets assuming the corresponding v_rad at the time of the observations and T_H,Geo = 1000 K. The fitted parameters are the vertical column density, N_H,Eur, and temperature, T_H,Eur, of Europa’s exosphere as well as the Earth exosphere column density, N_H,Geo.

4.4 Constraints on Europa’s H exosphere

A signal from Europa’s H exosphere was derived (i.e., I_H > 0) in all visits, with the mentioned exception of visit 24 where the best fit of the model suggests that a measurable H exosphere signal is not present (I_H = 0), see Figure 3a. The highest brightnesses and thus vertical H column densities (Figure 3d; before considering extinction in the Earth’s exosphere) of around 1.5 × 10¹² cm⁻² were derived for visits 6 and 12.

When plotted against the (absolute) radial velocity of Europa with respect to Earth at the time of the observation, it is apparent that the measured signal is systematically lower at low Doppler shifts (Figure 3e). In addition, the values from visits 21-29 (years 2018-2020) appear to be systematically lower than values from earlier visits (2012-2015 and 1999) at similar radial velocities. This difference was even more pronounced in the signal brightnesses (Figure 3a) but has been balanced somewhat by the g factors, which also (like the exosphere brightness value) were systematically lower in the later years (2018-2020, visits 21-29, Figure 3c) due to the low solar flux (Figure 3b). The solar flux at Europa was particularly low in 2018 (visits 21 to 26), when the solar activity was near its minimum and Jupiter was near aphelion (i.e., farthest away from the Sun; see Figure 3b).

In the following analysis, we consider two periods (visits 218 in years 2012-2015, and visits 21-29 in years 2018-2020) separately, fitting the total intensity profiles (Section 4.3) individually for each period. Table 3 shows the best-fit values and the values for the maximum and minimum temperature T_H,Eur considered to provide a resulting curve in reasonable agreement with the observations. The three profiles are shown in green (best-fit T_H,Eur), blue (minimum T_H,Eur), and red (maximum T_H,Eur) in Figures 3f and g.

The fitted H column density for Europa is primarily constrained by the data points that are not affected by the attenuation. Thus, it is hardly correlated with the other fit parameters and well constrained. The resulting N_H,Eur values differ between the two periods beyond the derived uncertainties (Table 3).

Europa’s exospheric H temperature is coupled to some extent to the column density of the Earth H exosphere: For higher Europa exosphere temperature and thus a broader emission line, a higher N_H,Geo is needed to fit the maximum extinction observed at lower Doppler shifts. Furthermore, changes in T_H,Eur most of all define how gradual (or steep) the transition from maximum attenuation to the unattenuated maximum signal is. This transition happens at Doppler shifts between ∼3 km/s and ∼10 km/s, where only few observations were taken (Figure 3e, 2012-2015). We find that even a very low temperature of T_H,Eur = 100 K, as for an exosphere that is in thermal equilibrium with Europa’s surface, yields reasonable fits to the data. Maximum temperatures that yield χ² values near 1.0 are 5100 K and 6500 K, respectively, for periods 1 and 2. The difference in the resulting best-fit T_H,Eur values between the periods (1000 K and 2100 K) is thus well within the uncertainties of the fitted temperature.

5 Search for localized (aurora) emissions

5.1 Global residual images

The model generated as described in Section 3 includes all known Lyα sources from Europa, as well as those from the background and foreground. For the search for localized emissions from, for instance, patchy aurora (from outgassing at plumes or other sources), we subtracted the modeled images from the STIS images and analyze the residual emissions.

The modeled Lyα contributions are all relatively homogeneous, with only the reflected-sunlight component exhibiting structure due to the inverted albedo map used. Small-scale surpluses over regions smaller than Europa’s disk should therefore not appear in the model images (particularly outside the disk) and thus remain as surpluses in the residual-emission images. After subtracting the model, we extract a 72 × 72-pixel image centered on Europa for each visit and rotate the images so that the north pole points upward. All images are shown in Figure 5 using the same brightness color range from −400 R to +400 R. Because the subtracted model was fitted to match the global brightnesses, the residual emissions are distributed around zero.

The differences in appearance of the emissions originate primarily from the different noise levels, mostly determined by the level of the foreground emission (Table 2) and the total used exposure time (Table 1). Higher foreground emission (from Earth exosphere) and shorter exposure times yield higher noise levels. Visit 6 exhibits the highest noise level because the observations were taken long before Jupiter opposition, when the geocorona emissions were high and substantial exposure cuts led to a short usable exposure time.

The signal on Europa’s disk contains a substantial contribution from sunlight reflected off the icy surface, with intensities ranging from ∼0.7 kR on the leading side to roughly ∼1.3 kR on the trailing side (Table 2). Any apparent surplus on the disk may therefore originate from (or at least be influenced by) local variations in the Lyα albedo that are not captured in the smoothed inverted visible maps used in the model. In some images, it appears that our model does not fully capture the spatial variation in surface reflection, leading to systematically brighter or darker residual regions on the disk. For example, in the image from visit 2 (Figure 5, V2), the left (dawn) hemisphere is brighter and the right (dusk) hemisphere is fainter; similar patterns are seen in later trailing-side images, such as visit 11. Consequently, as in our previous studies (Roth et al. 2014b,a), we refrain from interpreting inhomogeneities in the residual Lyα brightness on the disk.

The most notable on-disk surplus emissions appear in visits 11 and 23, near the central longitude slightly below the equator. In visit 11, a larger region in the left part of the disk shows positive signal, while the right side reveals negative residual signal. This indicates a misfit of the reflection model. During visit 23, several high-count pixels are randomly distributed across the detector, indicating that this feature is likely an artifact. Moreover, no co-located surplus is seen in the oxygen 1304 Å emission, which can serve as supporting evidence of Lyα enhancements produced by electron-impact dissociative excitation of H₂O (Roth et al. 2014b; Roth 2021).

5.2 Limb bin analysis for all visits

For a quantitative assessment of positive Lyα emission outliers, we again followed previous studies in focusing on the region above Europa’s limb between 1.0 R_E and 1.25 R_E. As in Roth et al. (2014b), we subdivided this annulus into 18 azimuthal bins, each with an angular width of 20°, and calculate the average residual brightness in each bin.

In Roth et al. (2014b), the detection of a significant emission surplus was defined when the brightness in one bin exceeds the average of the other bins in the same image by more than three times the propagated statistical uncertainty (>3σ detection). Figure 6a shows the bin brightness, I_bin, divided by its statistical uncertainty, σ_bin, for the brightest limb bin in each image. In contrast to the previous studies, we did not subtract the average brightness of the other bins of the same image in our analysis here, because this average is already very close to zero after considering and subtracting the H exosphere signal.

In our analysis, the maximum value does not exceed the threshold of 3 in any of the visits, indicating that no statistically significant local surplus in Lyα emission is detected. This includes the image from visit 3 in December 2012, which was interpreted in our earlier study as showing an emission surplus. We further discuss the results for visit 3 and compare them to those in Roth et al. (2014b) in Section 5.4. We investigated the statistical behavior of the limb-bin brightness first to validate our results.

With 23 visits (or 23 images) and 18 limb bins in each image, we have a total of N=414 limb bin measurements. Figure 6b shows the frequency distribution of the obtained brightnesses significance values I_bin/σ_bin in a histogram with a column width of ∆I_bin/σ_bin =0.2. The distribution is in very good agreement with the theoretical Gaussian distribution of the form $$f(x)=A\exp \left(\frac{(x-\mu {)}^2}{2{\tilde{\sigma }}^2}\right),$$(11)

for the sample size of N (red curve in Figure 6b). We then fit a Gaussian distribution to the histogram of the observations (blue curve in Figure 6b) and find that the fitted parameters for the maximum A, mean μ, and standard deviation ∂^ are almost identical to the theoretical values for this sample, see labels in Figure 6b.

The obtained residual limb bin brightnesses are, hence, fully consistent with statistical random fluctuations. Although the mean of the fitted distribution is slightly positive (μ = 0.08), the only two outliers above 3 are negative and thus not indicative of emissions but would rather indicate absorption. No positive outliers beyond I_bin/σ_bin = 3 are observed. Statistically, a value exceeding 3 would be expected for a sample with size N=414 with a probability of 60%.

Fig. 5 Residual Lyα images after subtraction of the model image. The images were binned (2 × 2 detector pixels) and boxcar smoothed (3 × 3 binned pixels) for display. The rotation from the detector to the Europa frame leads to missing signal in the corners, as seen in, e.g., the noisy image of visit 6 (V6).

5.3 Statistical tests with synthetic data

We tested the statistical robustness of our results by generating synthetic images of Europa’s emissions and applying the full analysis pipeline to them. To do this, we simulated statistical noise in each of the 23 model images by drawing random values from a Poisson distribution, using the model pixel value (in total counts) as the expected mean.

We then applied the same limb-bin analysis to the 23 synthetic images that we applied to the STIS images. The resulting distributions are again very similar to the expected Gaussian behavior. Figure D.1 shows one example set of 23 synthetic data images and can be compared to Figure 6 (bottom). Furthermore, using 10 000 generated realizations of the full set of 23 images (one model per observation), we test the expected frequencies of obtaining outliers at various significance levels, which could otherwise be falsely interpreted as detections.

In the 10 000 samples of image sets, we find 6790 3σ-outliers (i.e., an bin with a brightness that is ≥3σ above the average limb brightness), corresponding to a probability of $$\frac{6790}{10000\times 23\times 18}=1.6\times {10}^{-3},$$

which is close to the theoretical one-sided 3σ detection probability of 1.4 × 10⁻³. This includes several outliers in an image set (or even in one image). In 4970 image sets (i.e., in roughly half of the total 10 000 image sets), there was no 3σ-outlier, as is the case for our STIS image set. In the other 5030 image sets, at least one image contained at least one bin with a 3σ-outlier, but in many cases there were several outliers in one image set thus leading to the 6790 total 3σ-outliers mentioned above.

A 4σ-outlier is found in only 178 image sets (as compared to the 5030 image sets for 3σ), corresponding to a probability of 4.3 × 10⁻⁵, which is again close to but slightly higher than the theoretical value of 3.2 × 10⁻⁵. Thus, a 4σ outlier would indeed be considered significant, but is only found in our STIS observation image in December 2012 if the disk position is assumed as in Roth et al. (2014b).

We also tested the effect of positional offsets by shifting the synthetic noisy image along x and/or y and performing the analysis using the non-shifted model image as the reference. Such an offset introduces a systematic error and effectively broadens the distribution shown in Figure 6 (bottom), increasing the frequency of bin brightnesses at large σ (as expected). For a shift by only one pixel in one direction, the frequency of a 4σ-outlier increases by a factor of ∼3. For a shift of one pixel in both x and y, a 4σ-outlier occurs about six times more often. When shifting two pixels in x and one pixel in y, which is the main difference between this analysis and that of Roth et al. (2014b), the factor is already ∼20 and there is a 35% chance of obtaining one 4σ-outlier in one set of 23 images. This means that the 4.2σ outlier derived when assuming the (possibly false) disk position used in Roth et al. (2014b) is well consistent with statistical variations. In other words, the disk offset is a reasonable explanation for this outlier.

We note that this test is not fully self-consistent, as the model image would slightly change under such a shift as well, because it is adjusted in brightness to best match the observation at the initially assumed position. However, it roughly quantifies how an incorrectly assumed position affects the probabilities of obtaining a false detection.

Fig. 6 Statistical analysis of the limb bins: (Top) Brightness significance Ibin/σbin of the limb bin with the highest value for each visit (blue squares). For visit 3, the orange square shows the value when assuming the disk position from Roth et al. (2014b). (Bottom) Histogram of the frequency of limb bin brightnesses for all 414 bins from all visits. The dotted red line shows the theoretical Gaussian distribution for sample size N=414. The dotted blue lines shows a Gaussian distribution fitted to the histogram. The values for max amplitude, A, mean, μ, and standard deviation, σ, are given for both Gaussians.

5.4 Comparison to previous results for visit 3

While our limb-analysis results are very similar to those of Roth et al. (2014b) and Roth et al. (2014a) for visits 1, 2, 4, and 5, we find no emission surplus in visit 3, in contrast to the previous study. The data processing and analysis presented here differ in several aspects from the approach used in Roth et al. (2014b). We use a slightly updated description of the surface reflection, accounting for the phase-angle dependence following Oren & Nayar (1994), and we apply a reduced contrast to the inverted visible map. However, these differences are minor and have little impact on the off-disk signal in the limb-bin region.

More importantly, the earlier analysis did not account for emissions from the H exosphere, which was detected only later by Roth et al. (2017b). The H exosphere remained undetected in the first five STIS datasets (visits 1-5), partly because the exospheric signal was comparatively low: Europa’s radial velocity happened to be particularly small during these observations (see Table 2 and Figure 3e). From visit 6 on, the exosphere signal was higher and clearly detectable. Nevertheless, even in visit 3, the peak brightness of the H exosphere at the limb was 145 R (value from Table 2 multiplied by π), and the average brightness in the region between 1.0 and 1.25 R_Eu was 130 R.

In the previous analysis, the only considered source of off-disk emissions above the constant background was the spread signal from on-disk solar reflection, which contributes on average up to only about 50 R in the limb region. When the H exosphere signal is taken into account, a higher intrinsic emission is expected in the limb region. This naturally reduces the residual brightness there. More importantly, it influences the determination of Europa’s disk position in the data, which turns out to be the crucial parameter, as it allows for homogeneously enhanced emissions above the limb.

Figure 7 presents different assumed positions of Europa’s disk and the resulting chi-squared values for the comparison between the data and model images. The best-fit position (with a minimum chi-squared of χ² = 1.05) is shown in the center with a blue frame. For the other positions, the relative increase in χ² is shown as a percentage of the minimum value, since the absolute differences are small. The position assumed in Roth et al. (2014b) is indicated by the orange frame. Using the algorithm described in Section 3, the best-fit position differs from that used in Roth et al. (2014b) by 2 pixels in the x-direction and 1 pixel in the y-direction. The crucial difference is the 2-pixel shift in x, as this moves the bright pixel patch in the lower-left quadrant either off or onto the disk (noting that the images in Fig. 7 are shown in the detector frame and are not rotated to north-up).

In Roth et al. (2014b), we assumed that the lower part of the disk (in the detector frame) was affected by the anomaly (i.e., the localized auroral emissions) and therefore used only the upper part of the disk to determine the central x pixel. The details of the position-determination method are provided in the Supplementary Material of Roth et al. (2014b), in particular their Figure S6. In the present analysis, we use the full image to determine the position using the comparison with a 2d model image, without assuming the presence of an anomaly a priori. The contribution from the H exosphere considered here also influences the position search, as it accounts for relatively homogeneous signals above the limb more accurately. Despite providing a worse agreement with our model image, the disk position derived in Roth et al. (2014b) still provides reasonable agreement with the data and cannot be excluded to be correct, given the noise level in the data.

Figure 8 shows the limb-bin analysis for both positions. With the new position (top panel), the residual brightness in the bins is symmetrically distributed around the average limb brightness (black dashed line), ranging from approximately −200 R to +200 R. Bin 13 is the brightest, with I_bin = 244(±98) R, corresponding to a brightness signal-to-noise ratio of 2.5 (see Figure 6a for visit 3). This is the same bin for which the maximum brightness was inferred in Roth et al. (2014b), where a value of I_bin = 604(±140) R was reported.

With the disk position from 2014, the bin brightnesses on one side of the disk (bins 7-15, except 10) are positive, while they are negative on the opposite side (bins 16-18 and 1-6). On the positive side, bin 13 is again the brightest, with I_bin = 430(±102) R. The difference from the value reported for bin 13 in Roth et al. (2014b), 604(±140) R, of 174 R can be partly explained by the subtraction of the H-exosphere brightness in the new analysis, which has an average value of 130 R in the limb region. When using a reference brightness that is lower by 130 R for the limbbin measurements, together with the old disk position (dashed purple line in Figure 8, bottom), the resulting profile becomes very similar to the profile shown in Figure 3 of Roth et al. (2014a).

Even without correcting the reference brightness, the brightness of bin 13, I_bin = 430(±102) R, corresponds to I_bin/σ_bin = 4.2. This would imply a detection with a slightly higher significance than the value of 4.0 reported in Roth et al. (2014b). The absolute brightness of the maximum bin is lower (430 R vs. 604 R), but the brightness significance is higher because the earlier study used a more conservative error estimate, leading to a larger uncertainty of 140 R compared to 102 R here. In the previous error analysis, we included systematic uncertainties when subtracting the background and disk reflected contributions, while we here do not consider systematic uncertainties in the model fit and when subtracting the modeled contributions. The independent statistical test using the full limb-bin sample of all currently available images (Section 5.2 and Figure 6b) supports the robustness of the uncertainties derived in the present analysis. Thus, when adopting the previous position of Europa’s disk on the detector, an emission surplus is detected with a significance similar to (and slightly higher than) that found in the earlier study.

Fig. 7 Overview of the disk positions around the best-fit position for visit 3 shown in the center with blue frame. The (x, y) position of the central pixel is given below in parenthesis and the corresponding disk limb is shown with the dotted yellow circle. Images are identical otherwise and all scaled and smoothed with 3 × 3 boxcar filter to enhance visibility. The disk position assumed in Roth et al. (2014b) is indicted by the orange frame. The reduced chi-squared for all pixels is shown for the best-fit position. For all other positions the relative decrease of chi-squared value is shown in percentage of the best-fit value. The image from the forward model used to find the best-fit position is shown to the right. Note: the images are shown in detector frame and not rotated, unlike in Figure 5.

Fig. 8 Comparison of limb bin analysis for visit 3 with the new disk position (top) and the disk position (bottom) used by Roth et al. (2014b). The left panels show the residual Lyα emission and the top image is identical to Figure 5, panel V3. The 18 limb bin regions are shown by the yellow dotted lines. The right panels show the average residual brightness in each bin with error and the average brightness of all bins (black dashed). The purple dashed line shows a reference average brightness when an H-exosphere brightness of 130 R is subtracted.

6 Discussion

6.1 Updates in processing and analysis

In this work, we present an analysis of the complete set of the Lyα signal in the HST/STIS FUV G140L 52″ × 2″ observations of Europa in sunlight. For this purpose, we built a forward model for the contributions from surface reflection and H exosphere resonant scattering of the solar Lyα flux as well as the emissions from the Earth exosphere and interplanetary medium. The forward model, with its parametrized individual contributions, was then fitted simultaneously to the observations to obtain Lyα albedo, H exosphere brightness, and fore- and background levels. In previous studies of Europa’s atmosphere and surface with the Lyα STIS data (Roth et al. 2014b,a; Becker et al. 2018), we had not considered H exosphere contributions and independently corrected for fore- and background emissions. In addition, here we considered the surface illumination angle using Oren-Nayar reflectance model while previously we assumed a uniform disk. We also used an adjusted contrast map for the surface FUV reflection (adjusted inverted visible maps). Finally, possible inhomogeneities near the limb from potential plumes are constrained by characterizing differences between data and forward model images.

6.2 Europa’s H exosphere: Density, temperature, and production and loss rates

The derived vertical H column densities in Europa’s exosphere of 1.4 × 10¹² cm⁻² (2012-2015) and 1.1 × 10¹² cm⁻² (2018-2020) are 4-5 times higher than the densities reported by Roth et al. (2017b) based on the attenuation of Jupiter’s daylow in the Europa exosphere in transit. However, the Roth et al. (2017b) densities need to be corrected and after the correction our results are actually in agreement: Roth et al. (2017b) underestimated the H density due to an overestimation of the effective cross-section for the attenuation of the Jovian background dayglow, as pointed out in Roth et al. (2023). When correcting the cross-section to an about 5 times lower value (Roth et al. 2023), the inferred H density from the data of Roth et al. (2017b) is higher by a factor of about five. Therefore, the values reported here are consistent with the H exosphere detection in transit (Roth et al. 2017b) when the appropriate attenuation cross-section is used.

We have separately fitted the two periods (2012-2015 vs. 2018-2020) as the resulting densities appeared to converge to systematically different values. The resulting fitted densities and their relatively small uncertainties suggest that Europa’s H exosphere has indeed a different density in the two periods (Table 3), with a ∼25% decrease in density from the first period to the second. This is a result independent of the attenuation of Earth and derived primarily from the observations taken at large Doppler shifts, which almost all happened either in late 2014/early 2015, or in 2018. The measurements taken at radial velocities lower than 10 km/s (from 1999, early 2014, 2012, 2019 and 2020, see Figure 3, right column) are strongly affected by the attenuation at Earth and possible differences in Europa H density might be masked by this. The H exosphere signal is not detected during visit 24 likely because the expected weak exosphere signal (green curve in Figure 3f) could not be differentiated from the background trend in our fitting algorithm.

The best-fit temperatures for Europa’s H exosphere were 1000 K and 2100 K for the two periods, respectively. The temperature is constrained primarily by observations taken at Doppler shifts between 3 km s⁻¹ and 10 km s⁻¹, where the transition from strong attenuation in Earth’s exosphere to full transmission occurs. The temperature derived for the first period is more reliable because several observations were obtained within this radial-velocity transition range only a few months apart, making the simultaneous fit a more reliable measure. In the second period, only the visit 28 observation had an intermediate Doppler shift, which in addition carried a large uncertainty and was taken a year after the high-Doppler-shift observations so that other factors might have changed.

Using the density and temperature, we estimate the hydrogen source rate, Q_H, again assuming a cometary-like escaping exosphere. In this framework, Q_H is given by $$Q_H=4\pi n_{\mathrm{H},0}{R}_E^{2}v=4\pi N_{\mathrm{H}}{R}_{E}v,$$(12)

where v is the outward velocity of H. We adopt the most probable thermal velocity for an ideal gas, $$v=\sqrt{\frac{2k_B{T}_{\mathrm{H}}}{m_H}}.$$(13)

For a temperature of T_H = 1000 K, derived from the 2012-2015 datasets, this corresponds to v = 4.1 km s⁻¹. Using a column density of 1.4 × 10¹² cm⁻², we obtain a nominal production rate of Q_H = 1.1 × 10²⁷ Hs⁻¹.

Assuming a Maxwellian velocity distribution for the upward velocity component, with a characteristic speed of 4.1 km s⁻¹, approximately 93% of the H atoms exceed the Hill escape velocity (v_esc,Hill = 1.96 km s⁻¹) and are therefore able to escape. This implies an H escape rate from Europa of ∼1.0 × 10²⁷ H s⁻¹. The remaining, non-escaping H atoms re-impact the surface and undergo chemical reactions.

Ionization can be neglected for this estimate. It takes roughly 5000 s for an H atom to travel from the surface to Europa’s Hill sphere (a distance of ∼ 12 000 km, or ∼8 Europa radii) at a velocity of 4.1 km/s. Electron-impact ionization rates are well below 10⁻⁵ s⁻¹, so the H atoms escape before being ionized. Rates for charge exchange with ambient ions are again an order of magnitude lower than electron-impact rates, and thus charge exchange is even more negligible (Smith et al. 2019).

With the lower density derived from the 2018 observations (N_H = 1.1 × 10¹² cm⁻²), the source rate is also lower by 25%. A higher H temperature and thus v on the other hand increases the source rate but even more the loss rate. With T_H = 2100 K, v is 5.7 km/s and the escaping fraction of H is 96%. At T_H = 5000 K, 97% of the H escapes directly.

It is important to emphasize that the assumption of a single temperature describing the effective Lyα line width and thus attenuation might not capture the actual velocity distribution of the H atoms. Line-resolved measurements of atomic sodium (Na) emissions at Europa showed variation in line width and suggested an increase in Na temperature with distance from Europa (Lovett et al. 2025). The strongest signal in our data is from within 1 R_E altitude, and our data would not be sensitive to, for example, an increase in temperature at higher altitudes. For example, the model of Smyth & Marconi (2006) yields an H temperature of 1000 K near the surface, which increases to about 3000 K at 1000 km altitude. The authors explain the temperature profile with a mixture of colder surface-sourced H atoms and hotter H atoms produced by molecular dissociation as well as collisions with the thermal O₂-atmosphere. Our results for the H temperature are consistent with this profile from Smyth & Marconi (2006). We discuss the expected densities and temperatures from different sources for H in the following subsection. However, we note already that significantly higher H temperatures are expected from dissociation of atmospheric molecules to produce atomic H due to excess energies, which would suggest even higher escape rates than estimated with the derived H temperatures.

The obtained column densities are very similar to the H column densities measured at Ganymede ((1-2)×10¹² cm⁻², Barth et al. 1997b; Alday et al. 2017; Roth et al. 2023) and Callisto (∼1 × 10¹²cm⁻², Barth et al. 1997a; Roth et al. 2017a). Given these comparable column densities and the fact that the radii of Ganymede and Callisto are 69% and 55% larger than Europa’s radius, the H number density at the surface, n_H,0, is higher at Europa than at the other two moons assuming the same 1/r² density profile in all cases. The larger radii of Ganymede and Callisto, however, lead to larger total content of the global H exosphere and therefore larger source rates, Q_H (Equation (12)). On the other hand, Europa has a lower escape velocity than Ganymede and Callisto and therefore the escaping fraction of the H exosphere will be higher at Europa (in other words, higher gravity prevents escape of H more at Ganymede and Callisto). Since Europa orbits closer to Jupiter than Ganymede and Callisto, the same amount of escaping H will lead to higher densities along Europa’s smaller orbit.

6.3 Europa’s H exosphere: Source mechanisms and implications

We now examine whether the derived constraints on H exosphere density (considering also the difference between the two periods), temperature, and production and loss rates are consistent with potential source mechanisms for the H atoms at Europa. Possible sources for H are:

• Dissociation of H-bearing molecules in the atmosphere, primarily H₂ or H₂O, driven by either UV photons or electron impact.

• Dissociation of the surface ice followed by direct release of H atoms. The can be induced by charged particle sputtering (e.g., Bar-Nun et al. 1985) or photolysis (Johnson & Quickenden 1997).

We first discuss the dissociation of atmospheric molecules and then briefly investigate a direct surface source in the following.

Photodissociation of H₂ or H₂O is unlikely to be the dominant source of atomic hydrogen. First, electron-impact dissociation has been shown to be significantly more efficient than photodissociation at Europa (Shematovich et al. 2005; Smyth & Marconi 2006; Szalay et al. 2024). In addition, H atoms produced by photodissociation of H₂ or H₂O acquire excess energies of ∼1 eV or more (e.g., Table 2 in Marconi 2007), implying temperatures T ≳ 10⁴ K, well above the upper limit derived here. Furthermore, although the solar UV flux decreased from 2015 to 2018 - potentially leading to a reduced photodissociation rate -the magnitude of this change was smaller (Woods et al. 2005) than the observed ∼25% difference in H density. Moreover, the present short-term solar UV variability (∼20%) is not reflected in the inferred H densities at Europa.

Electron-impact dissociation of H₂ or H₂O is expected to be more efficient at Europa than photodissociation. The energies of the H atoms produced by electron-impact dissociation are not well constrained, and the inferred H temperature therefore cannot be used to validate or invalidate this process. Nevertheless, electron impact is likely still insufficient to account for the required H production rate. Because the abundance of H₂O is likely significantly lower than that of H₂, we focus on H₂ for the following estimates. To sustain an H production rate of Q_H = 1.1 × 10²⁷ H s⁻¹ via dissociation of H₂, an H₂ production rate (e.g., from radiolysis) of $\frac{1}{2}{Q}_H=5.5\times {10}^{26}$ H₂ s⁻¹ would be required to compensate for the loss to H. In addition, an independent H₂ source rate of 4.5(±2.4) × 10²⁶ H₂ s⁻¹ is needed to explain the observed H_J pickup ions, assuming escape of neutral H₂ followed by ionization in a torus (Szalay et al. 2024).

Because dissociative recombination of H_J is negligible (Szalay et al. 2024; van der Zande et al. 1996), dissociation into H versus escape and ionization can be considered two independent pathways for H₂. The total required H₂ source rate to sustain the H exosphere and explain magnetospheric H_J pickup ion densities would then be $$5.5\times {10}^{26}{\mathrm{H}}_2{\mathrm{s}}^{-1}+4.5\times {10}^{26}{\mathrm{H}}_2{\mathrm{s}}^{-1}=1.0\times {10}^{27}{\mathrm{H}}_2{\mathrm{s}}^{-1}.$$

Under this scenario, 55% of the H₂ would need to be dissociated to form H, while only 45% would be ionized to form H_J. This partitioning is unlikely, because H₂ escape occurs on a timescale of ∼10⁴ s, whereas the dissociation time scale (inverse rate) is on the order of 10⁶ s. Consequently, only about 1% (not 55%) of the H₂ is expected to dissociate before escaping.

In addition, the electron conditions depend on the local plasma environment, which is modulated primarily on the synodic rotation period of Jupiter of 11.2 h (Roth et al. 2016). The visits were intentionally scheduled at different plasma environments (see System-III longitude in Table 1) and there is no systematic pattern that could explain the difference in H density between the two periods that we inferred here. Europa is deep inside the magnetosphere and seasonal changes or changes in the solar wind do not affect the local plasma environment.

Instead, the difference in density between the two observing periods is more likely explained by changes in the source rates of the parent molecules. The radiolysis yield of surface ice induced by charged-particle impacts has been shown to depend strongly on surface temperature, increasing approximately exponentially with temperature (e.g., Shi et al. 1995; Johnson et al. 2009; Raut & Baragiola 2013). Using the temperature dependence reported by Shi et al. (1995) for sputtering by O⁺ ions (see their Fig. 2), we find that an increase in surface temperature from 100 K to 103 K, representing a 3% change near Europa’s average surface temperature (Spencer et al. 1999), would result in a ∼28% increase in radiolysis yield. This is comparable to the observed difference in exospheric density. Such a temperature-driven change would directly affect the production rate of H₂ as parent species supplying the H exosphere.

However, it is unclear whether the temperature variations of the surface (and the corresponding changes in radiolytic production yield) are indeed sufficiently large to explain the observed density difference. Seasonal surface temperature variations are poorly constrained. During the earlier period, Europa’s distance from the Sun was ∼5.3 AU, compared to ∼5.4 AU in the later period; this change corresponds to only a 3% difference in solar irradiance at Europa. Variations in solar irradiance due to the solar cycle are negligible in this context, because the power in the visible and infrared portions of the solar spectrum remains nearly constant (Woods et al. 2022). Modeling of Europa’s surface temperature Ashkenazy (2019) suggests maximum temperature changes of only ∼5 K over a full Jovian year, despite ∼20% changes in solar flux at Europa between aphelion and perihelion. If radiolytic production yields indeed change significantly due to surface temperature variations, all radiolytically produced atmospheric species - including O₂ - would exhibit strong seasonal variability.

A hint that seasonal variation might indeed happen is provided by the 1999 HST observation. In 1999, Europa was closest to the Sun at 4.96 AU. The H density derived here, although affected by the Earth exosphere, is clearly higher than the best-fit profiles suggest (gray data points vs. solid green lines in Figure 3f,g). Similarly, the oxygen aurora measured in the same 1999 observation and produced primarily by dissociative excitation of O₂ was by far the brightest (Roth et al. 2016, their Figure 6).

A similarly strong temperature dependence is given for ice sublimation as a source for atmospheric H₂O (e.g., Leblanc et al. 2017). However, production of H from purely sublimated H₂O is insufficient due to the low expected H₂O abundance, and this would not produce a global, uniform H abundance.

Finally, direct production of H from erosion of the icy surface might be an alternative explanation. The low inferred H temperatures constitute a strong argument for this source, because production via dissociation of atmospheric molecules likely lead to temperatures at least an order of magnitude higher (for both photodissociation and electron impact dissociation). It is, however, difficult to assess the efficiency of surface H sources.

There is only one published study on direct sputtering of atomic H off the surface from charged particle irradiation by Bar-Nun et al. (1985). Whether the yield would be sufficient to explain the H source rate obtained, and what H temperature would be expected is unclear. Bar-Nun et al. (1985) mention that the yield they find is independent of the ice temperature, which means the density difference between the two periods could hardly be explained by the direct sputtering source (unless the sputtering flux changes). Incident UV photons can also cause direct dissociation of the water ice in the surface via photolysis, leading to production of atomic H (e.g., Johnson & Quickenden 1997). The efficiency and characteristics of such photolysis-induced desorption of H is yet also not known. Similar to photodissociation of atmospheric molecules, changes in solar UV flux between the two periods might partially explain the found H exosphere variability (see discussion above). However, whether it can account for the quantitative change of 25% cannot be estimated due to lack of knowledge about the process.

Regardless of whether H is produced from the surface directly or from atmospheric hydrogen-bearing molecules, the total hydrogen source (H plus H₂) must be roughly a factor of ∼2 higher (in mass or atoms) than recent estimates by Szalay et al. (2022, 2024) based on H₂ only. Assuming that O₂ is produced at the stoichiometric ratio, this also implies an O₂ production rate about twice as high as in Szalay et al. (2024), bringing it closer to the source rates commonly assumed in atmospheric simulation studies (e.g., Shematovich et al. 2005; Smyth & Marconi 2006; Plainaki et al. 2013).

The inferred H density and escape rate are about an order of magnitude higher than assumed in modeling studies of the atmosphere and escape to a neutral torus (Smyth & Marconi 2006; Smith et al. 2019). With higher values for density and escape, atomic H may possibly become the most abundant species in Europa’s neutral torus and thereby relevant for the interaction with the charged particles (Mauk et al. 2003; Lagg et al. 2003).

6.4 Earth exosphere

We also included the density of the Earth’s H exosphere as a fit parameter and obtained column densities (Table 3) within the range expected from exosphere models, in particular from the NRLMSIS 2.0 atmosphere model (Emmert et al. 2020). For our best-fit results, we find Earth-exosphere H column densities of 9.9 × 10¹² cm⁻² and 3.9 × 10¹³ cm⁻² for the two periods 2012-2015 and 2018-2020, respectively. The difference is likely connected to the variation of Earth’s exosphere with the solar cycle.

The approximately four-times higher column density near solar minimum can be explained by the enhanced H abundance near the exobase - in the lowest part of the exosphere - reported in previous studies (Nossal et al. 2012; Waldrop & Paxton 2013). This elevated H abundance may be due to dynamical production and upward transport at low solar activity (Qian et al. 2018).

During solar maximum, the temperature in the upper atmosphere increases, pushing more H to higher altitudes and into the far-extended exosphere. There, the correlation between H abundance and solar activity appears to be inverted compared to the region near the exobase, with more H observed during solar maximum (Zoennchen et al. 2024). However, the attenuation of the source signal occurs in the optically thick H exosphere above the exobase, namely, within the first few hundred kilometers above Hubble (∼540 km). There, the H density is higher during solar minimum, in line with our results.

6.5 Local inhomogeneities

We did not identify any outliers in any of the 23 residual emission images. The distribution of the residual brightness in the limb bins is fully consistent with purely statistical variations.

Notably, we also did not find an emission surplus in the image from visit 3 taken in December 2012, in contrast to our earlier analysis and interpretation (Roth et al. 2014b). The primary reason for the different results is a difference in the assumed position of Europa on the detector. The importance of this a priori unknown position is crucial in the analysis of HST observations, as previously demonstrated by Giono et al. (2020) for STIS transit FUV filter observations of Europa.

The other main difference is the inclusion of an H exosphere signal from Europa when correcting for other sources. In contrast to our earlier analysis (Roth et al. 2014b), in which we assumed an anomaly in the data and used only a part of Europa’s disk signal to determine its position, in this study we applied the same approach and algorithm to visit 3 as to all other visits. Therefore, the analysis presented here is preferred and considered more impartial. This means that the HST/STIS Lyα auroral observations do not provide evidence of a localized abundance of H₂O from outgassing. The conclusions on a stable H₂O atmosphere above the sunlit trailing hemisphere from Roth (2021) are not affected by our results here as they were derived only from OI1304 Å and OI1356 Å images (and not Lyα).

Roth et al. (2014b) discussed several lines of supporting evidence of the detection of localized H₂O aurora. In particular, an excess of OI1304 Å emission was reported at the same location as the Lyman-α excess. In addition, the persistence of the above-limb Lyman-α signal, as well as an apparent correlation between brightness variations and the plasma environment were emphasized. Although these findings were never considered to constitute evidence on their own, we briefly discuss them in light of the new results for completeness.

The OI1304 Å excess was detected at only a 2.4σ significance level in Roth et al. (2014b) and was hence not a statistically significant outlier itself. With the revised disk location (Figure C.1), there is still a some increased OI1304 Å brightness left of the south pole, but the corresponding bin #13 is well within the standard variation of brightness around the limb and other limb bins have higher brightness.

The potential correlation to the plasma environment is still seen but the persistency of an emission surplus near the south pole throughout the visit disappeared: A Lyα excess in the south polar bin is observed during orbits 3 and 4 of the five HST observing orbits, while no excess is detected in orbits 1, 2, and 5 (for comparison see Figure 2, panels F-J, in Roth et al. 2014b). Thus, although the excess in the two orbits still coincides with plasma sheet crossings (and, therefore, with the highest ambient plasma densities, where enhanced H₂O auroral excitation would be expected), the lack of a statistically significant detection means that this does not constitute evidence of an outlier or local H₂O aurora.

Studies of potential plumes at Europa following the initial detection by Roth et al. (2014b) have built on those results to varying degrees. The potential plume signals found in STIS FUV filter images (Sparks et al. 2016) was shown to be possibly explained by statistical fluctuations and a misalignment of the assumed disk position (similar to the results here) by Giono et al. (2020). Jia et al. (2018) presented plasma fluid simulations that reproduce observed plasma wave and magnetic field signatures measured by the Galileo spacecraft during flyby E12, assuming the presence of a localized, jet-like plume of neutral gas. Recently, Paterson & Collinson (2026) showed that data from the Galileo Plasma Instrumentation does not confirm the density enhancement derived from the Galileo wave data by Jia et al. (2018), discounting a plume encounter during flyby E12 according to the authors.

Paganini et al. (2020) reported a single detection of H₂O vapor emission above the leading hemisphere during an infrared observing campaign with the Keck telescope spanning multiple dates; the study interpreted this result as indicative of intermittent activity. Subsequent non-detections in infrared observations by the James Webb Space Telescope (Villanueva et al. 2023) are consistent with the inferred infrequent nature of such events. At present, direct measurements of water-molecule emissions (particularly via infrared observations and ideally with spatial resolution of the source) likely constitute the most sensitive remote sensing observations for placing constraints on water plume activity at Europa. After arrival in 2030, NASA’s Europa Clipper will systematically and globally constrain plume activity at Europa with its diverse instrument suite (Pappalardo et al. 2024; Roth et al. 2025).

7 Summary

We carried out a comprehensive analysis of Lyα observations of Europa obtained in 1999 and between 2012 and 2020. Our analysis includes forward modeling of all known contributions to the Lyα signal. The study focuses on deriving constraints on Europa’s hydrogen exosphere and on performing a systematic search for localized emission outliers. The key findings regarding the H exosphere can be summarized as follows:

• H exosphere emissions are detected in all but one observation, with varying brightness;

• Lyα illumination from interplanetary H can enhance the flux scatterable by Europa’s H exosphere up to 20%;

• A primary factor controlling the observed exosphere brightness is attenuation of the Europa signal in Earth’s atmosphere, which is significant when the line-of-sight velocity of Europa relative to Earth is less than 10 km s⁻¹;

• Europa exospheric H column densities of 1.4 × 10¹² cm⁻² and 1.1 × 10¹² cm⁻² were derived for the 2014/2015 and 2018 observations, respectively;

• An effective temperature of approximately 1000 K was inferred from the 2014/2015 observations, with an upper limit of 5100 K;

• The inferred H source rate of 1.1 × 10²⁷ s⁻¹ is more than twice the H₂ source rate recently derived from H₂⁺ pickup ions. This suggests that a significant fraction of the hydrogen produced by erosion of surface H₂O ice ultimately becomes atomic hydrogen;

• The high H density and source rate, combined with the relatively low H temperature, are difficult to reconcile with the current understanding of possible source processes. The low temperatures seem to preclude a production from dissociation of molecular species in the atmosphere.

In addition to the findings on Europa’s exosphere, our analysis provides evidence that the H column density in Earth’s exosphere above HST’s altitude (550 km) is higher during solar minimum. This result can be tested by the recently launched NASA Carruthers Geocorona Observatory, which will investigate the H exosphere in detail (Waldrop 2024; Joshi et al. 2024). The results of the search for localized emissions are summarized as follows:

• No statistically significant emission surplus is found in the 23 HST images of Europa obtained in sunlight;

• The derived brightness significances in the limb bins are in good agreement with statistical fluctuations;

• The difference from the previous results published in Roth et al. (2014b) arises from differences in the assumed position of Europa on the detector and from the inclusion of the H exosphere in the analysis;

• When adopting the same assumed position for Europa as in Roth et al. (2014b) and neglecting the contribution from the H exosphere, an outlier would again be identified in the December 2012 observations. Using an updated error analysis, the derived significance of this outlier would be 4.2σ, slightly higher than the significance of the outlier reported in the previous study (4.0σ).

Future observational campaigns and planetary missions will be essential to determine whether Europa is presently cryovolcani-cally active.

Acknowledgements

This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program(s) HST-GO-8224, HST-GO-13040, HST-GO/DD-13619, and HST-GO-13679. LR acknowledges support by the Swedish National Space Agency through grants 2021-00153 and 2024-00112. SRCM is supported by the Swedish Research Council through Starting Grant 2025-04693. SJ is supported by Swedish National Space Agency grant 2020-00187.

References

Appendix A IPH Lyα brightness and map

Figure A.1 shows an example map of the Lyα sky brightness produced by sunlight scattered by the IPH. The region of maximum emission is located near the Sun on the upwind side of the IPH flow. Very close to the Sun, the density of the IPH is depleted by charge-exchange processes; this depletion forms a cavity that extends downstream with the IPH flow, reducing the brightness on the downstream side.

The brightness from these modeled maps is used to estimate the scattering of IPH Lyα radiation by Europa’s H exosphere, in addition to the scattering of direct solar illumination, as described in Section 4.2. The g factor is determined using a sky map of the relative radial velocity of the IPH and by calculating the spectral flux in all directions while accounting for Doppler shifts (see Figure 2a). For reference, the average brightness on the dayside and nightside hemispheres for this day were 790 R and 450 R, respectively. Integrating over each hemisphere, one yields total fluxes of 4.1 × 10⁸ photons/cm²/s for the dayside and 2.3 × 10⁸ photons/cm²/s for the nightside.

Fig. A.1 Sky map of the Lyα brightness seen from Europa for visit 12 (Jan 26, 2015) from interplanetary medium H scattering model (Pryor et al. 2024). The direction to the Sun is shown by the small white circle, the anti-Sun direction is shown by the white asterisk.

Appendix B Lyα albedo variation

The derived Lyα albedo values (Table 2) confirm the previous results of Becker et al. (2018). Our values are overall slightly higher because we use the Oren-Nayar reflectance model, which assumes reduced reflectance near the limb compared to a uniformly bright disk. The new 2018-2020 data fill a gap around 160-200° sub-observer west longitude (see Fig. 4A of Becker et al. (2018), yielding the highest overall albedo at 200° (visit 25), see Figure B.1. Fitting a sinusoidal variation to all albedo values, the fitted peak remains near 270°, as in the previous study, but there also appears to be an additional increase towards the anti-Jovian hemisphere (180-200° W), similar to what Becker et al. (2018) found for the albedo at 1335 Å.

Fig. B.1 Lyα albedo values derived from the 23 visits and a fitted sinusoidal periodic variation.

Appendix C Oxygen images from visit 3

Figure C.1 shows the oxygen aurora images with the new position of Europa and the limb bin analysis.

Fig. C.1 Oxygen aurora images centered at 1304 Å and 1356 Å from visit 3 using the new disk position and brightnesses in the 18 limb bins. For comparison see Figure 1 (panels L and O) and Figure 3 (C and E) in Roth et al. (2014b).

Appendix D Example set of synthetic images

Figure D.1 shows an example set of images based on synthetic noise data.

Fig. D.1 Synthetic Lyα images generated by adding Poisson noise to the forward model image for each of the 23 visits. The synthetic images are then processed and displayed in the same way as the STIS observation images shown in Figure 5, including binning (2x2 detector pixels) and boxcar smoothing (3x3 binned pixels) for display. The noise level relates to the total counts in each image reflecting the noise in the observations.

All Tables

All Figures

Fig. 1 Overview of the analyzed HST/STIS observations. Left: orbital positions of Europa from the start of the first exposures to end of the last exposures for each of the 23 HST/STIS visits. Note: the gaps between exposures of a visit are not shown. Visit 21 combined exposures before and after transit. Right: complete STIS detector spectral image from visit 22. The Lyα signal from the geocorona and the interplanetary hydrogen fills the complete slit (yellow-green region). Reflected continuum sunlight from Europa’s surface forms the trace around row y=200. Atmospheric oxygen emissions are present at 1304 Å and 1356 Å, exceeding the surface reflection.

In the text

Fig. 2 Data analysis and model fitting. Top: brightness for the entire Lyα slit along the y-axis, averaged over the x-axis (i.e., over the slit width). Middle: Zoom-in of the top panel. Bottom: brightness along the x-axis, averaged over a y range covering the disk diameter. The STIS data are shown as a black histogram. The fitted model is represented by dashed lines corresponding to different components: green for the background only, blue for the background plus H exosphere contribution, and red for the full model including surface reflection and IPH correction (see text). Note: the profiles are smoothed due to averaging over a detector axis (and, thus, across the disk), as well as by the convolution of the model image with the STIS point spread function. The vertical dotted lines indicate the disk center in all panels. The model was fitted within the central disk and the gray shaded regions (see top panel).

In the text

Fig. 3 Analysis of Europa’s H exosphere signal: (a) H exosphere brightness in the disk center from the fitted exosphere profiles (Eq. (3)) for each visit. (b) Solar Lyα flux (black asterisks) measured at 1 AU near the observation date (see text). The inverse of the squared heliocentric distance of Europa (orange triangles, in unit 0.01 AU−2) indicates the scaling of the flux used for calculating g-factors. (c) g-factors for resonant scattering of the sunlight (blue diamonds), of the IPH illumination (green) and combined (black). (d) H exosphere vertical column densities from conversion of measured brightnesses. (e) Same as in (d) but plotted against the line-of-sight velocity of Europa with respect to the Earth at the time of the observation. (f and g) Same as (e) but only for visits 2-18 taken in 2012-2015 and for visits 21-29 taken in 2018-2020, respectively. Fitted attenuation profiles for different temperatures of the H in Europa’s exosphere, see text for more information. The 1999 visit is shown in grey in (f) and (g) for comparison.

In the text

Fig. 4 Analysis of the spectral aspects of the Lyα signals: (a) Typical Lyα spectral fluxes from the Sun (solid) and from the IPH (dash-dotted) at the heliocentric distance of Jupiter. The IPH flux is shown as a Gaussian profile with T = 15 000 K, at an offset to the solar line because of the IPH movement relative to the Sun. The scatterable flux by the H in Europa’s exosphere varies with the relative line-of-sight velocity between the moon and the source (i.e., the Sun or IPH in a certain direction, as shown for three example velocities). The vertical dashed lines (purple, grey, orange) show how the contributions from the Sun and the IPH can vary due to the respective Doppler shift. (b) Lyα transmission of the Earth H exosphere in the rest frame of Earth for an assumed column density of 1 × 1013 cm−2 and temperature of 1000 K (solid). Lyα emission profiles from Europa’s H exosphere (T = 1000 K) from resonant scattering before (dotted) and after (dashed) attenuation in the Earth exosphere from the shown transmission. Three example profiles (blue, grey, red) are shown for three different relative line-of-sight velocities of Europa with respect to Earth. The strongly shifted emission (blue) is not affected by the Earth at all.

In the text

	Fig. 5 Residual Lyα images after subtraction of the model image. The images were binned (2 × 2 detector pixels) and boxcar smoothed (3 × 3 binned pixels) for display. The rotation from the detector to the Europa frame leads to missing signal in the corners, as seen in, e.g., the noisy image of visit 6 (V6).
In the text

Fig. 6 Statistical analysis of the limb bins: (Top) Brightness significance Ibin/σbin of the limb bin with the highest value for each visit (blue squares). For visit 3, the orange square shows the value when assuming the disk position from Roth et al. (2014b). (Bottom) Histogram of the frequency of limb bin brightnesses for all 414 bins from all visits. The dotted red line shows the theoretical Gaussian distribution for sample size N=414. The dotted blue lines shows a Gaussian distribution fitted to the histogram. The values for max amplitude, A, mean, μ, and standard deviation, σ, are given for both Gaussians.

In the text

Fig. 7 Overview of the disk positions around the best-fit position for visit 3 shown in the center with blue frame. The (x, y) position of the central pixel is given below in parenthesis and the corresponding disk limb is shown with the dotted yellow circle. Images are identical otherwise and all scaled and smoothed with 3 × 3 boxcar filter to enhance visibility. The disk position assumed in Roth et al. (2014b) is indicted by the orange frame. The reduced chi-squared for all pixels is shown for the best-fit position. For all other positions the relative decrease of chi-squared value is shown in percentage of the best-fit value. The image from the forward model used to find the best-fit position is shown to the right. Note: the images are shown in detector frame and not rotated, unlike in Figure 5.

In the text

Fig. 8 Comparison of limb bin analysis for visit 3 with the new disk position (top) and the disk position (bottom) used by Roth et al. (2014b). The left panels show the residual Lyα emission and the top image is identical to Figure 5, panel V3. The 18 limb bin regions are shown by the yellow dotted lines. The right panels show the average residual brightness in each bin with error and the average brightness of all bins (black dashed). The purple dashed line shows a reference average brightness when an H-exosphere brightness of 130 R is subtracted.

In the text

	Fig. A.1 Sky map of the Lyα brightness seen from Europa for visit 12 (Jan 26, 2015) from interplanetary medium H scattering model (Pryor et al. 2024). The direction to the Sun is shown by the small white circle, the anti-Sun direction is shown by the white asterisk.
In the text

	Fig. B.1 Lyα albedo values derived from the 23 visits and a fitted sinusoidal periodic variation.
In the text

	Fig. C.1 Oxygen aurora images centered at 1304 Å and 1356 Å from visit 3 using the new disk position and brightnesses in the 18 limb bins. For comparison see Figure 1 (panels L and O) and Figure 3 (C and E) in Roth et al. (2014b).
In the text

Fig. D.1 Synthetic Lyα images generated by adding Poisson noise to the forward model image for each of the 23 visits. The synthetic images are then processed and displayed in the same way as the STIS observation images shown in Figure 5, including binning (2x2 detector pixels) and boxcar smoothing (3x3 binned pixels) for display. The noise level relates to the total counts in each image reflecting the noise in the observations.

In the text