EXOSpy includes physics- and data-based exospheric models described in this section. We refer the reader to the cited references below as well as the EXOSpy’s website https://exospy.readthedocs.io/for further details.

Chamberlain (1963) developed a physics-based exospheric model to determine density distribution of species at altitudes where collisions between them can be neglected. Here, H atoms are assumed to be thermalized through neutral-neutral collisions at the boundary of the collisional and collisionless atmosphere, termed as the exobase ( ∼ 500  km altitude), and then adopt ballistic or escaping trajectories to space. The H density distribution model is reported as an analytical altitude-dependent function that only requires values of H density (n_exo) and neutral temperature (T_exo) at the exobase, both of which can be specified by external models, such as the well-known model NRLMSIS2.0 (Emmert et al., 2021).

Bailey and Gruntman (2011) developed an exospheric model based on parametric fitting to a spherical harmonic function using single-day measurements of scattered Ly-α emission acquired by the Lyman-Alpha Detector (LAD) onboard the Two-Wide angle Imaging Neutral-atom Spectrometers (TWINS) mission. The approach used in this study leverages the linearity between the H density and its emission that occurs when the medium is optically thin, i.e., the density is low enough such that each H atom scatters a solar Ly-α photon only once. At earth, the exospheric optically thin region is located beyond 3 R _E geocentric distance. In addition, the TWINS/LAD apogee of ∼ 7.2 R _E and its nadir viewing geometry allows only partial observation of the exosphere, i.e., a single LAD line-of-sight (LOS) acquires emission starting from a radial position ∼ 7.2 R E towards infinity on the opposite side, but it cannot observe emission behind the instrument where there is still significant exospheric density with approximate values of 10 ∼ 100 cm⁻³ (Cucho-Padin et al., 2022b). Thus, exospheric density reconstructions using TWINS/LAD Ly-α data typically set 8R _E as an upper boundary since the estimation of the region behind the sensor may feature high uncertainty. Hence, retrieved 3-D, global H density distributions for this model are only valid from the 3 to 8 R _E region. Zoennchen et al. (2015) used a similar approach using multi-day data acquired by LAD/TWINS to generate a model for solar minimum conditions (with data acquired in 2008 and 2010) and solar maximum conditions (with data acquired in 2012). The region of validity for these two models is also restricted to a region from 3 to 8 R _E . Furthermore, Zoennchen et al. (2022) estimated an H density radial profile using Ly-α measurements from the Ultraviolet Imaging Spectrograph (UVIS) onboard the CASSINI spacecraft during an earth flyby on 18 August 1999, on its way to Saturn (Esposito et al., 2004; Werner et al., 2004). These radiance data were fitted to a spherical function to yield an exospheric profile valid from 3 to 15 R _E . Recently, Cucho-Padin et al. (2022b) derived a 3-D model of the outer exosphere using a single global image acquired by the Lyman-Alpha Imager CAmera (LAICA) onboard the Proximate Object Close Flyby with Optical Navigation (PROCYON) spacecraft when it was located at ∼2000 R _E from earth. The H density model is reported as spherical harmonics coefficients with a valid region that ranges from 6 to 15 R _E .

Estimation of neutral densities is also possible through remote sensing of soft X-ray emissions, which are generated by the charge exchange interaction between heavy charged solar wind particles (e.g., O⁺⁷ and O⁺⁸) and exospheric neutral atoms. Connor and Carter (2019) estimate H density in the vicinity of the terrestrial magnetosheath, i.e., the subsolar point at ∼10 R _E , using simultaneous measurements of soft X-ray radiance acquired by the XMM-Newton mission and in situ ion density, acquired by Cluster and Geotail missions. The H density for two events during solar maximum conditions (October 2001 and May 2003) are reported as a spherical function. Following a similar methodology, Jung et al. (2022) used XMM-Newton radiance data and THEMIS ion density data to determine an H density profile during solar minimum conditions (November 2008) and it is also reported as a spherical function. Since both studies utilized data from magnetosheath (located within the optically thin exospheric region), we consider the valid region of their models to range from 3 R _E to ∼12 R _E , where the bow shock (the outer boundary of magnetosheath) is typically located (Fairfield, 1971). Table 1 provides a summary of data-based models included in this package.

EXOSpy provides the user with functions to generate and visualize H density distributions. For example, the function draw1DHmodel receives a 1-D model alias (see Table 1) and an altitude range (in units of R _E ) as inputs to produce the corresponding H density profile, as a numpy array, as well as a density vs altitude graph. Similarly, the function draw3DHmodel generates a 3-D H density profile and plots for the meridional and equatorial planes and spherical shell for a given geocentric radius. In Figure 1, panels (A) (B) and (C) show plots for the 3-D H density distribution model derived by Zoennchen et al. (2015) for solar maximum conditions (EXOSpy’s alias Z15MAX). Specifically, panel (A) shows the equatorial plane, panel (B) shows a meridional (pole-to-pole) plane containing the subsolar point, and panel (C) shows a spherical shell at a geocentric radius of 5 R _E . In addition, panel (D) shows a comparison between 1-D H density profiles along the Sun-Earth line generated with the models included in the package. To determine the spatial distribution of densities, EXOSpy uses the Geocentric Solar Ecliptic (GSE) coordinate system, in which the x-axis points sunwards from earth’s center, the z-axis points to the ecliptic north pole, and the y-axis completes the orthogonal triad. Also, we implemented a python notebook (Example1. ipynb located in the project’s github) to show the use of EXOSpy’s functions to create all graphs in Figure 1. The reader can access the file through the link https://github.com/gcucho/EXOSpy/blob/main/Examples/Example1.ipynb.

Estimation of the column-integrated line of sight Ly-α intensity ( I ( n ̂ ) ) observed by a photometric instrument (e.g., a single-pixel photometer) requires information of its spatial position, the direction of its line-of-sight (LOS) vector and the exospheric volumetric emission rate (ɛ) which is generated by the interaction of solar Ly-α photons and atmospheric particles such as atomic H and molecular Oxygen (O₂).

To calculate I ( n ̂ ) within the terrestrial exosphere, EXOSpy solves the following equation (Bhattacharyya et al., 2017): I n ̂ = Ψ n ̂ 1 0 6 ∫ 0 L max ε 0 l T τ H l e − τ O 2 l + ε m l T τ H l e − τ O 2 l d l + I I P n ̂ , (1) where n ̂ represents the 3-D direction of the instrument LOS, the terms ɛ ₀(l) and ɛ _m (l) are the volume emission rates due to single and multiple scattering of Ly-α photons by H atoms given in units of photons/cm³/sec, the terms τ H ( l ) = σ H ∫ 0 L n H ( l ) d l and τ O 2 ( l ) = σ O 2 ∫ 0 L n O 2 ( l ) d l are the optical depth of the medium due to H atoms and O₂ molecules, σ _H represents the Ly-α photon scattering cross-section and σ O 2 represents the Ly-α photon absorption cross-section by O₂. The line integral is evaluated from the beginning of the LOS (l = 0) to an appropriate upper boundary with a significant density of Ly-α emitters (l = L _max). Also, T(τ) is the Holstein function (Holstein, 1947) that estimates the photon transmission probability in a medium with optical depth τ(l) and is given by the relation T τ = 1 π ∫ − ∞ ∞ e − x 2 e τ e − x 2 d x , (2)

In addition, the term Ψ ( n ̂ ) is the scattering function that estimates the probability of photon scattering in the direction n ̂ and is formulated as Ψ n ̂ = 11 12 + 1 4 cos 2 θ n ̂ , (3) where θ ( n ̂ ) is the angle between n ̂ and the +X _GSE axis (Brandt and Chamberlain, 1959). The term I I P ( n ̂ ) is the interplanetary Ly-α background, in units of Rayleigh, emitted by atomic H of interstellar origin.

The use of Eq. 1 requires volume emission rates for single and multiple scattering that can be estimated by a radiative transfer (RT) model. An RT model computes the single and multiple resonance scattering of photons by atmospheric particles and, for this purpose, requires solar photon flux and the cross-section of the particle-photon interaction as inputs to the model. For the analysis of the exosphere, the RT model uses three inputs: an H density profile (1-D for spherically symmetric distribution) to estimate the Ly-α photon scattering rate, an O₂ density profile to estimate the Ly-α photon absorption rate, and the incoming solar Ly-α flux usually expressed as irradiance at 121.56 nm in units of W/m². As a result, the RT model provides a spatial distribution of volume emission rates.

At high altitudes, where molecular oxygen is not found anymore ( τ O 2 = 0 , e − τ O 2 = 1 ), and the medium is optically thin (τ _H ≪ 1, T (τ _H) = 1), i.e., only single scattering of solar Ly-α occurs (ɛ _m = 0 and ɛ ₀ = g* × n _H), Eq. 1 can be simplified as I n ̂ = g * Ψ n ̂ 1 0 6 ∫ 0 L max n H l d l + I I P n ̂ (4) where n _H(l) is the atomic H density in units of 1/cm³ and g* (also referred to as g-factor) is the scattering rate in units of photons/sec given by the relation g * = f 0 σ H δ ν D π (5) where f ₀ represents the incoming line center solar Ly-α flux in photons/cm²/sec/Å, σ _H is the cross-section of the H Ly-α scattering cross section at line center, and δν _D is the Doppler width of the Ly-α emission line. Eq. 4 is only valid for the optically thin region of the exosphere, i.e., altitudes beyond 3 R _E geocentric distance.

Figure 2 depicts a block diagram that explains the steps to calculate line of sight column-integrated Ly-α intensities using EXOSpy. If the user selects a physics-based exospheric model that extends from 500 km altitude to 10 R _E geocentric distance, Eq. 1 is used to calculate intensities. As mentioned before, it requires volume emission rates for single and multiple scattering obtained from an RT model. This current version of EXOSpy does not include an RT model (white box in Figure 2), but it is considered to be part of a future release. To overcome this problem, we have included seven files with distributions of volume emission rates (ɛ ₀ and ɛ _m ) generated with Chamberlain-based H density profiles (see Section 2.1) and exobase parameters obtained from the NRLMSIS2.0 model corresponding to the equinox (March 20) from 2014 to 2020. This period represents solar maximum to solar minimum variations. In addition, O₂ densities, which typically extend up to ∼ 800 km altitude due to their smaller scale height in comparison to H, are also obtained from the NRLMSIS2.0 model. Solar irradiance data is included in each file and are retrieved from the LASP Interactive Solar Irradiance Datacenter (LISIRD) (https://lasp.colorado.edu/lisird/data/composite_lyman_alpha/). EXOSpy’s function Chamberlain implements the physics-based, 1-D Chamberlain model, which was used to generate the volume emission rates included in the package, and the function CalculateLOSfromVER implements Eq. 1 without considering Interplanetary Ly-α background.

On the other hand, if the user selects a data-based exospheric model, Eq. 4 is used to calculate Ly-α intensity since these models are only valid for the optically thin region of the exosphere (see Table 1). In this case, the user should provide solar conditions as irradiance data which can be retrieved from LASP/LISIRD. EXOSpy’s function GenerateIntensityOpticallyThin implements Eq. 4 without considering Interplanetary Ly-α background.

Then, the user should provide a set of spacecraft positions and LOS directions in GSE coordinates to calculate the line of sight column-integrated Ly-α intensities. Also, EXOSpy allows the user to create a synthetic image comprised of individual LOSs when a FOV (in units of degrees) and a pixel resolution (in units of degrees/pixel) are provided.

Hence, EXOSpy calculates only the exospheric Ly-α emission which is the first term in the right hand side part of Eqs. 1, 4. Interplanetary Ly-α background (I _IP ) can be retrieved from physics-based models (e.g. (Pryor et al., 1992)) or data provided by the Solar Wind ANisotropy (SWAN) instrument onboard the SOlar Heliospheric Observatory (SOHO) mission. The SOHO/SWAN instrument is located at Earth-Sun L1 point (Lagragian point), observing the interplanetary medium in Ly-α using a hydrogen cell. Its main product is a 1-day average all-sky map of the interplanetary Ly-α background.

The source code of EXOSpy is located in the following Github repository: https://github.com/gcucho/EXOSpy/. Also, a release of the package has been uploaded to a Zenodo repository and the download link can be found in (Cucho-Padin et al., 2022).

Furthermore, the package is part of the Python Package Index (PyPi) repository and its latest version (v2.3) can be installed to a personal computer using the command: pip install EXOSpy==2.3. For this, the user is required to have python 3.x and the pip application. Additionally, the following packages are needed: NumPy, SciPy, scikit-learn, and they will be installed automatically if they are not present in the user’s computer.

In this work, we have developed four examples to show how to use EXOSpy. These examples are available within the source code and their locations are indicated in the manuscript. The user can also go to our project’s website https://exospy.readthedocs.io/ for more information regarding these examples.

Current and future models of exospheric density distributions can be assessed using EXOSpy and Ly-α radiance data from several space-based instruments such as the Geocorona Imager (GCI) onboard the Imager for Magnetopause-to-Aurora Global Exploration (IMAGE) mission, the Global UltraViolet Imager (GUVI) onboard the Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics (TIMED) mission, and the LADs onboard the TWINS mission, among others.

We propose a validation methodology that is comprised of four steps: (i) select a Ly-α instrument and identify its orbit and instrument LOS(s), (ii) identify the exospheric region (optically thin or optically thick) observed by the instrument, (iii) set the testing exospheric model and calculate column-integrated Ly-α intensities using EXOSpy’s functions, and (iv) calculate the model performance. To demonstrate our approach, we show an example using TWINS/LAD measurements (step i). NASA’s TWINS mission is comprised of two satellites (TWINS1 and TWINS2) with a highly elliptical orbit, an apogee of 7.2 R _E , orbital precession of ≈1 year, and a 90° longitudinal separation between satellites to allow for stereoscopic observations. Each TWINS spacecraft has two LADs (single-pixel photometers) that rotate 180° around a constant nadir pointing vector (McComas et al., 2009). TWINS radiance and ephemeris data are publicly available in NASA’s Space Physics Data Facility (SPDF) repository. Panel (A) in Figure 3 shows TWINS1/LAD1 observations for three rotations acquired on 11 June 2008, 08:20 UT. These observations have been pre-processed to remove stellar contamination and interplanetary Ly-α background as per the procedure detailed in (Cucho-Padin and Waldrop, 2018). The black line indicates the TWINS1 orbit, the green lines are LAD1’s LOSs, and the blue arrow indicates the initial clockwise rotation of the LAD instrument around the nadir vector. Due to the LOSs viewing geometry, these observations serve to evaluate exospheric models beyond 3 R _E geocentric distance, i.e., optically thin region (step ii). EXOSpy enables the user to input a new exospheric model as spherical harmonic coefficients or as discrete 3-D spatial locations in GSE spherical coordinates using the class exosgrid. In this example, we use the model developed by Zoennchen et al. (2015) for solar minimum conditions, which is included in the package (Z15MIN). Then, we estimated the line of sight column-integrated Ly-α intensities using EXOSpy’s function G e n e r a t e I n t e n s i t y O p t i c a l l y T h i n . Panel (B) in Figure 3 shows a comparison between the calculated intensities based on Z15MIN model in blue line (step iii) and the actual data measured by the TWINS1/LAD1 instrument on 19 June 2008, in the dashed red line. Finally, these values can be evaluated using several performance indicators such as minimum square error or the Pearson correlation coefficient, among others (step iv). Panel (B) shows an evident disagreement between TWINS/LAD observations and calculated intensities in the minima regions that are possibly associated with the multi-day approach used to generate this exospheric model. Specifically, Z15MIN was implemented using Ly-α observations during the periods June 19–21, 2008, and June 19–21, 2010. The intensities shown in the blue line were obtained using Eq. 4 and only considering solar flux for 19 June 2008 (used in Eq. 5). Hence, the exospheric model can reflect the average behavior of density distributions for several days in two different years and may exclude the formation of daily structural patterns in the density distributions. Nevertheless, further analysis should be done to draw firm conclusions regarding the validity of this model.

The python notebook Example2. ipynb uploads the pre-processed ephemeris and attitude data from TWINS1/LAD1 and shows the use of EXOSpy’s functions to generate column-integrated Ly-α intensities. The reader can access the file in https://github.com/gcucho/EXOSpy/blob/main/Examples/Example2.ipynb.

Observations of extra-terrestrial bodies in UV might be contaminated by the bright exospheric emission, and its removal or avoidance is crucial for the correct analysis of the measurements.

To demonstrate how EXOSpy may help quantify this background emission, we estimate the Ly-α intensities as seen from a Low earth Orbit (LEO) satellite. Specifically, we simulate the orbit and viewing geometry of the TIMED/GUVI instrument, which are similar to those of the HST mission. Panel (A) in Figure 4 shows in blue line the orbit of the satellite within the XY ecliptic plane at an altitude of 600 km (optically thick region). Red dots indicate positions when the spacecraft acquires data, and green lines are examples of measurements’ LOSs only for two positions for better visualization. TIMED/GUVI is a spectrograph with a rotating mirror that allows LOS observations of the earth’s limb from look angles (LA) that range from 0 to 20°. A LOS has LA = 0° when its direction is perpendicular to the zenith and LA > 0 degrees when its direction is closer to the earth. The satellite location is represented by the angle between + X _GSE axis and the zenith, also referred to as the solar zenith angle (SZA). For each satellite position (red dot), we simulated LOSs covering 20° LA with a 1-degree resolution. Then, we calculate the column-integrated Ly-α intensity using EXOSpy’s function CalculateLOSfromVER as well as the volume emission rates for the optically thick region during solar minimum conditions (year 2020).

Panel (B) in Figure 4 shows the resulting Ly-α intensities. Two typical features of exospheric emission are shown in this plot. First, the Ly-α intensity is reduced as the satellite transits from dayside (SZA = 0°) to dusk (SZA = 90°) due to the thicker atmosphere the solar Ly-α flux encounters on its path, ultimately resulting in fewer photons reaching the dusk/nightside region. Second, the presence of denser molecular oxygen at low altitudes absorbs Ly-α photons such that LOSs with smaller tangential distances may produce lower intensities. For instance, for a satellite position SZA = 0°, a LOS with LA = 0° yields ∼ 10 kR while a LOS with LA = 20° yields ∼ 8 kR since its path crosses lower down in the atmosphere, which has a larger amount of O₂ molecules. The python notebook Example3. ipynb shows the use of EXOSpy’s functions to generate TIMED/GUVI’s LOSs and their column-integrated Ly-α intensities. The reader can access the file in https://github.com/gcucho/EXOSpy/blob/main/Examples/Example3.ipynb.

The EXOSpy package can be used as part of the designing process of space-based UV instruments for several mission objectives such as estimation of exospheric H density distributions or observation of extra-terrestrial bodies in FUV wavelengths other than Ly-α.

Here we describe how the EXOSpy package can be used to design a wide FOV imager to observe the exosphere in Ly-α. This task typically aims to determine several optical parameters such as FOV (in degrees), pixel resolution (in degrees/pixel), sensor responsivity at Ly-α (in counts/sec/R), integration time (in seconds), the optimal ephemeris, and the optimal pointing viewing geometry. In this case, the designing methodology is an iterative process wherein an initial set of the above-listed parameters are used to generate an image that is then evaluated against certain requirements such as signal-to-noise ratio (SNR) or its capability to yield additional scientific products, e.g., H densities with lower uncertainty. If these requirements are not met, the parameters are modified, and a new iteration is performed. For this application, we simulate an image to observe the optically thin region of the exosphere with the following parameters: a square FOV of 18°, a pixel resolution of 18/72°/pixel, a geosynchronous orbit with a radius of 6.6 R _E over the ecliptic plane, and an anti-sunward pointing viewing geometry. Panel (A) in Figure 5 shows a simulation of this instrument where the red dots indicate two acquisition periods, green sectors show the square FOV, and the black line represents the orbit. To generate synthetic intensities, we select H density distributions during solar maximum conditions, i.e., Z15MAX model (see Table 1). The solar irradiance used in this example was 0.0086 W/m² obtained from LISIRD/LASP website for 20 June 2012. Panels (B) and (C) show the resulting 72 × 72 pixel image for each satellite position. The python notebook Example4. ipynb shows the use of EXOSpy’s function to generate a wide field image for the different pixel LOSs by calculating their corresponding column-integrated Ly-α intensities. The reader can access the file in https://github.com/gcucho/EXOSpy/blob/main/Examples/Example4.ipynb. It is worth noting that these images show a dawn/dusk asymmetry in the intensities, which, at this altitude, is directly associated with an asymmetric distribution of H densities. In addition, we adopt a sensor responsivity at Ly-α equal to 5 × 10^–3 counts/sec/R (typical for a CCD imaging sensor) and an integration time of 10 min. As a result, the number of digital counts in the images range from ∼5100 to ∼7200 counts. Assessment of this parameter setting can be done using a threshold SNR. In optical systems, the main source of noise is the Poisson-distributed shot noise which, in turn, simplifies the derivation of the measurement’s SNR to N , where N is the number of digital counts. Hence, the SNR for these images range from ∼71 to ∼84.

Unlike spectrograph-based instruments that can resolve narrow wavelengths bands, imagers rely on optical filters or a combination of them to generate a narrow responsivity function. Nevertheless, imagers are still the best option to obtain a global view of the target. For example, future astrophysics missions might be interested in observing very distant galaxies in Ly-α using imagers. Depending on the galaxy’s distance and its velocity with respect to the observer, they may be red-shifted (z) from 121.56 nm. For instance, observation of a target with a red-shift parameter of z = 0.2 indicates that an emitted 121.56 nm photon would be observed as 145.87 nm at earth. The proximity of these wavelengths requires an effective selection of optical filters to generate a responsivity function able to suppress the exospheric Ly-α emission. For this purpose, an accurate specification of this bright emission under different solar conditions is required. This task can be achieved using the similar approach presented in Section 3.2 to estimate column-integrated Ly-α emissions for a given satellite position and instrument’s LOSs.