Retrieval of Thermospheric O and N2 Densities From ICON EUV

As activity in Earth orbit continues to grow, it is important to characterize the environment of near‐Earth space. One means of remotely sensing lower thermospheric neutrals is by measurement of O and N2 density through the observation of far‐ultraviolet (FUV) airglow of atomic oxygen at 135.6 nm and the N2 Lyman‐Birge‐Hopfield (LBH) bands (~130–180 nm), as has been done on the Ionospheric Connection Explorer (ICON), Global‐scale Observations of the Limb and Disk (GOLD), and Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) missions. This technique is not without limitations, however, as the FUV measurements suffer from contamination by ionospheric emissions at low latitudes and auroral emissions excited by precipitating energetic electrons and protons at high latitudes. Previous work has shown the potential for making measurements of O and N2 density in the lower‐middle thermosphere using observations of extreme‐ultraviolet (EUV) airglow. This measurement approach has a potential advantage in that it does not have an inherent ionospheric emission that must be accounted for. Additionally, these emissions are primarily excited directly by solar UV rather than electron impact and thus have the potential to enable expansion of neutral density observations into the auroral zone and polar cap where the FUV measurement cannot be applied. This article demonstrates a new approach and algorithm designed to retrieve thermospheric O and N2 density from 150 to 400 km using measurements from the ICON EUV instrument. The retrieval results throughout 2020 are summarized and compared to measurements from ICON FUV, GOLD, and SWARM.


Introduction
The middle thermosphere (150-400 km altitude) is an important region of near-earth space, containing the lowest orbiting artificial satellites and overlapping with the F-region ionosphere.As such, it is important to characterize this region with density and temperature measurements.In-situ mass spectrometer observations are inherently limited to single-point measurements, and these measurements can be difficult to maintain in the lower thermosphere since the high density causes orbital decay.Observation of far-ultraviolet (FUV) dayglow is a prominent means of remotely measuring thermospheric neutrals.The 135.6 nm O doublet and N 2 Lyman-Birge-Hopfield bands are excited by photoelectron impact on the dayside, and limb measurements of these dayglow profiles from TIMED-GUVI have been used to retrieve profiles of O and N 2 density (Meier et al., 2015).However, it is more common to use disk-viewing measurements of these emissions for the retrieval of ΣO/N 2 , defined as the ratio of the summed O and N 2 column densities from the top of the atmosphere to a depth corresponding to an N 2 column density of 10 17 cm 2 (Equation 8) (Meier, 2021;Strickland et al., 1995).This quantity is important in the study of ionosphere-thermosphere coupling since O + is the dominant ion in the F-region, whose primary source is ionization of O and whose primary sink is loss through chemical reactions with N 2 .ΣO/N 2 has been a data product on several spacecraft missions, most recently the Ionospheric Connection Explorer (ICON) (Stephan et al., 2018) and current Global-scale Observations of the Limb and Disk (Correira et al., 2021) missions.
In support of its mission to measure the connection between Earth's ionosphere and neutral atmosphere, the ICON spacecraft carried an extreme-ultraviolet (EUV) spectrometer to observe O + in the F-region ionosphere (Immel et al., 2018).The ICON EUV spectral coverage spans approximately 54-88 nm and measures limb profiles nominally from 100 to 500 km, at 20 km resolution (Sirk et al., 2017).Profiles of the O + 61.6 nm (sometimes referred to as 61.7 nm) and 83.4 nm features are fitted with a forward model to retrieve the altitude and density of the F2 peak (hmF2 and NmF2, respectively).The former is produced with the photoionization of atomic oxygen and is used to estimate the inital source of O + in the F-region (Stephan et al., 2017).Many additional spectral features other than these two targeted O + features are captured by the EUV instrument.Among these, Tuminello et al. (2022) have demonstrated that the 87.8 nm band is dominated by an N or N + feature at low altitudes and that sub-limb measurements of the 61.6/87.8brightness ratio follows ΣO/N 2 , similar to the 135.6/LBH ratio.One conclusion of this work was that a retrieval of O and N 2 profiles from the ICON EUV measurements is plausibly feasible.
This work presents our results on the development and validation of an algorithm to retrieve O and N 2 altitude profiles using ICON measurements of EUV dayglow.Section 2 introduces the ICON EUV L1 Data.In Section 3, we discuss how we will use Discrete Inverse Theory to perform a forward model retrieval of neutral densities using the EUV data.Section 4 details the forward model itself, including the estimation of g-factors for the spectral features of interest.Section 5 presents an overview of the trends in the retrieval results, and Section 6 follows with a comparison of retrieved densities to other measurements.

ICON EUV Data
The ICON EUV instrument is a 1D imaging spectrometer.The raw output of the instrument is a 2D array of photon counts, with 108 pixels in the spatial (imaged) dimension aligned to local vertical spanning roughly 20-550 km tangent altitude, and 168 pixels in the spectral dimension spanning wavelengths 54-88 nm at 2.4 nm resolution.In order to convert these measured photon counts into emission brightness (in units of Rayleighs), the instrument completes a monthly on-orbit absolute flux calibration via observations of the full moon, and a flatfielding calibration obtained from a nadir-pointing measurement (see Sirk et al. (2017); Korpela et al. (2023)).In order to generate the ICON EUV Level 1 (L1) data product, the calibrated flux is summed across each of 12 the bands as defined to isolate the brightest emissions in the EUV dayglow spectrum.In the ICON EUV data products, these color bands are named based on the brightest or central feature expected within the band, but in some cases the band can contain multiple spectrally adjacent emission features.
Two of these wavelength bands contain the 61.6 and 83.4 nm O + emissions which are used in the ICON daytime ionospheric retrieval.These are both primarily produced by ionization of O into an excited state of O + , followed by naturally occurring decay into a lower-level state that generates an EUV photon.The 83.4 nm photons are resonantly scattered by ambient ionospheric O + , so the shape of the measured altitude profile of the emission brightness contains information about the distribution of the F-region ionosphere.The 61.6 nm emission is optically thin to the ionosphere and can be used to characterize the original source function of 83.4 nm photons.The combination of these emissions thus allows for a retrieval of the peak height and density of the F-region ionosphere, which constitute the EUV Level 2 (L2) product (Stephan et al., 2012(Stephan et al., , 2017)).The ICON ionospheric algorithm scales its forward model of the 61.6 nm emission to include not only corrections for atomic oxygen density, but also for solar flux and potential systematic errors in sensor calibration and thus does not attempt to directly infer O from these emissions.However, the 61.6 nm emission ties directly to O densities, and in this work is used to infer O.
Of the remaining defined bands, five (labeled as 53.7, 55.5, 67.3, 71.8, and 79.1 nm) are by other O + features.The band centered on 87.8 nm has been shown to track with N 2 at low altitudes (Tuminello et al., 2022) and is the other emission of interest to this work.The ratio of 61.6-87.8nm correlates to the FUV ratio 135.6 nm to LBH sub-limb region (50-140 km) where the ICON EUV and FUV fields of view overlap.Tuminello et al. (2022) speculated that these two EUV features could be used to measure thermospheric O and N 2 .In the sections that follow, we will describe our creation of an EUV dayglow model which generates the emission profile of a chosen feature, and demonstrate how we use this model to characterize the different contributing emissions to the EUV 87.8 nm band, yielding a forward model which estimates the total emission profile in the ICON EUV L1 87.8 nm band for a given atmospheric composition.We then use this forward model along with our model of the 61.6 nm emission band as the basis to conduct an iterative fitting to the altitude profiles of these emissions using discrete inverse theory with maximum likelihood to infer the underlying O and N 2 density profiles.

Limb Retrieval Method
In order to retrieve atmospheric densities from ICON EUV airglow profiles, we conduct a forward model inversion using discrete inverse theory (DIT).This is the same method used for the ICON EUV ionosphere retrieval algorithm (Stephan et al., 2017) and has previously been used to retrieve altitude profiles of thermospheric neutral density composition using limb-profile measurements from TIMED-GUVI (Meier et al., 2015).
A forward model (F) predicts the value of observed quantities ( y → ) for given state parameters of a physical system ( x → ) .
In this formalism, (F) is the airglow model which predicts the measured brightness altitude profile of an emission feature ( y → ) for given profiles of O, O 2 , and N 2 density ( x → ) .
Our consideration here is an inverse problem, where we wish to recover the atmospheric state from the airglow measurement.The forward model F is a non-linear operator, and the inverse of F is not well-defined (and is in our case under-determined).The DIT approach finds the parameters whose forward model image is the best fit to the measurement vector within the uncertainties represented for y → .To do so, the parameters of x → are iteratively adjusted to improve the fit to y → until converging to their optimized values.This process involves repeated calls to the forward model to estimate the model Jacobian matrix via finite difference (Menke, 2012).
It is possible to implement the DIT algorithm using atmospheric densities on an altitude grid as the state vector x → ; however this would create a large parameter space with large regions corresponding to physically unrealistic atmospheres.We use MSIS 00 as a driver to reduce the atmospheric state vector to three parameters.The first parameter is a scalar of the solar F10.7 input to MSIS 00, which modifies the shape of the density profiles of each species.The second and third are scalars applied directly and uniformly to the MSIS 00 output O and N 2 density profiles, respectively.While our approach includes the capability to scale the O 2 density profile, our early tests found this to have little impact on the retrieved altitude profiles because it only has relevance to photon absorption at the lowest altitudes below our region of interest, so we do not use this as a parameter of the fit.The rest of the MSIS 00 inputs, such as latitude, longitude, and A p are also left unchanged, defined only by the conditions of each specific measured altitude profile.This effectively yields a catalog of physically reasonable atmospheres which can cover a wide range of realistic atmospheric density profiles by the specification of only three parameters.
Figure 1 shows an example of a single retrieval.The top left plot shows the measured dayglow profiles and the initial forward model output.The three parameters are adjusted until they converge to the values seen in the text at the bottom right.In principle, it is possible for non-convergence to occur; however, very few retrieval attempts resulted in non-convergence, so we simply discard these failed attempts in our analysis.The top right plot shows the forward model updated to the atmosphere given by the scaled parameters, and the O and N 2 densities associated with this atmosphere are shown in the solid lines in the bottom left panel.
For a given airglow feature l with parent A, we express the volume emission rate (VER) where n A is the number density of species A (assumed to vary only with altitude).Note that A is the species prior to excitation.For example, 61.6 nm is created by an excited-state transition of O + ion but has a parent source of O that is ionized to the upper-level excited state.Equation 2 serves as the definition of the g-factor g l , which contains both the photochemistry of A and the availability of the source of excitation (solar UV and photoelectrons).The solar component of the g-factor is given by: where F is the attenuated solar flux of frequency ν at position r → , σ l (ν) is the excitation cross section to the excited state that produces emission l, and ν 0 is the frequency corresponding to the excitation threshold.
The EUV emissions have a small contribution from photoelectron impact excitation, and the photoelectron gfactor follows a similar form: where Φ( r → ,E) is the differential photoelectron flux at energy E at position r → , σ l,e (E) is the photoelectron excitation cross section, and E 0 is the excitation threshold energy.The g-factor is the sum of the solar and photoelectron contributions.
For the emissions under consideration for this work, the production is dominated by solar photoionization, with any photoelectron contribution being less than ̃10% above 200 km.
Because we are evaluating emissions for which routine measurements have not been made, we have adapted g-factors for the well-characterized 83.4 nm (currently used to define the initial ionization of O in the ICON ionosphere algorithm) and N + 108.5 nm (created by the photodissociation of N 2 (Meier, 1991)) emissions as our proxies for O and N 2 parent emission features, respectively.These g-factors are parameterized based on runs of the Atmospheric Ultraviolet Radiance Integrated Code (AURIC) (Strickland et al., 1999), as a function of solar zenith angle, solar F10.7 index, and column density.These g-factors must then be scaled to match the observed emission characteristics of the targeted model emission features (see Sections 4.2 and 4.3).
The g-factors allow us to model the initial production of photons, but these photons will undergo absorption and scattering by the atmosphere as they pass from the location of production to the sensor.At these wavelengths, we assume that only absorption by neutrals is of concern based on the shapes of the profiles (Tuminello et al., 2022).Under this assumption, the total brightness along a single line of sight j is given: where τ l,j (s) is the optical depth of feature l at point s along path j.We consider absorption by neutral O, O 2 , and N 2 .The optical depth is calculated: where σ l,k is the cross-section for absorption of a photon at the wavelength of l by species k and the integral is evaluated between the point s and the spacecraft along line of sight j.Cross-sections for O have been measured at To find the contribution from a single point, the volume emission rate is calculated at that point (Equation 2), and the optical depth is found by integrating neutral density along the red path to the satellite (Equation 7).Another integration is performed to sum up these contributions along the entire line of sight (Equation 6).
high resolution (Soto et al., 2023) and for O 2 and N 2 have been taken from Conway (1988).The calculation is repeated for the line of sight j of each pixel to generate the modeled airglow profile.

Modeling the 87.8 nm Passband
In order for the dayglow model to accurately predict EUV instrument response for given neutral density profiles and observation conditions, the airglow features composing each bin must be well-characterized.The 87.8 nm passband on the ICON EUV instrument has previously been shown to track with N 2 at low altitudes (Tuminello et al., 2022).The nominal 84.5-87.2nm range of this bin contains several N and N + features (Kramida et al., 2021) that are assumed to be generated by the photodissociative excitation of N 2 , and the nearby 87.8 nm O doublet has previously been observed in the limb dayglow (Gentieu et al., 1984).In this section, we show that the 87.8 nm passband is comprised of a blend of these O and N features and describe the process by which we determine the contribution of each.
Recall (see Section 4.1) that our available g-factors are for production of O + 83.4 and N + 108.5 nm dayglow derived from O and N 2 , respectively.We make a simplifying assumption and model the g-factors of the lesserstudied features of interest as scalar multiples of the well-tested 83.4 nm (O parent) or 108.5 nm (N 2 parent) g-factors.This approach has previously been used to model the 61.6 nm O + emission for the ICON EUV ionosphere retrieval (Stephan et al., 2017).To correctly model a given feature or blend of features, we must determine the appropriate scaling of the O + 83.4 and/or N + 108.5 nm g-factors.
To accomplish this, we first fit the ICON EUV L1 data collected over the entirety of year 2020, assuming the MSIS 00 is climatologically correct.We keep the neutral density profiles fixed but adjust the blending and scaling of the contributing emission features within the passband.Previous work showed that MSIS O and N 2 densities should be scaled by ̃85% for low solar conditions (Emmert et al., 2008;Meier et al., 2015), as is the case for all of 2020 -therefore, we have completed this precursory part of the analysis using the MSIS outputs for O and N 2 scaled by this amount.
Initially, we attempted to model the 87.8 nm passband using only spectral features with an N 2 source; we used scaled N + 108.5 nm g-factors and absorption cross-sections at a single wavelength.Figure 3 shows the forward model output of three candidate features with absorption cross-sections at 86.5, 87.0, and 87.5 nm.The forward model is linear in scaling of the g-factor, so each curve has been normalized to peak at the same brightness so that we can compare the shape of the profiles.All three fail to reproduce the modest peak around 200 km and the topside emission scale height of the measured 87.8 nm profiles.Modeling the passband as a blend of these features also fails to reproduce the shape seen in the data.
Figure 4 shows the same ICON EUV profile as in Figure 3, with a linear fit of the forward model output of the N 86.5 nm triplet and the O 87.8 nm doublet (Kramida et al., 2021).The blend of O and N features matches the shape of the profile above the peak and replicates the decreasing brightness below the peak.The O 87.8 nm has previously been observed in the limb dayglow (Gentieu et al., 1984), so we conclude that the EUV 87.8 nm bin is in fact seeing this emission in addition to the previously documented N emissions.As for the N, it is likely that the instrument is seeing contributions from several N and/or N + features, all presumed to have N 2 sources due to the low density of N in the terrestrial atmosphere at these altitudes.The profiles are distinct, as the absorption crosssections differ, but they are too similar for the fitting algorithm to consistently find a plausible fit using contributions at multiple wavelengths.Modeling of the EUV 87.8 nm bin could likely be improved by a detailed study of high resolution spectra, but for the purposes of this analysis, we consider the passband to be comprised of O 87.8 and N 86.5 nm, since this line results in the best fits to the data.
In order to determine the appropriate scaling of the O + 83.4 and N + 108.5 nm g-factors (and thus the blending of respective 87.8 and 86.5 nm features), we have conducted a survey over every fifth day of ICON EUV data throughout 2020, in order to cover a range of geophysical locations, local times, and solar zenith angles.To improve signal-to-noise ratio in these dim emissions, we sum up five consecutive 12-s ICON EUV exposures to create an effective 60 s exposure.Using the observation conditions (local time (LT), latitude, longitude, and solar zenith angle (SZA)) of the middle exposure and central spatial pixel and a scaled MSIS atmosphere, we run the dayglow model for O 87.8 and N 86.5 nm and use DIT (see Section 3) to find only the g-factor scalars which produce the best fit to the EUV data.Figure 4 depicts an example of this.
The median value of the g-factor scalars from the fits are 0.794 and 0.0708 for the 86.5 and 87.8 nm features, respectively.Figure 5 shows the distribution of these scalars normalized by dividing by the median value.The deviation of actual atmospheric profiles from MSIS climatology is one known cause of the spread in these distributions.Another possible source is that scaling the g-factors from another emission only provides an approximation of the g-factors of the desired emission, since the photochemistry differs between the atomic transitions.The most obvious difference between these emission features and the O + 83.4 and N + 108.5 nm features is that the energies of solar EUV photons that create these states will have different cut-offs, as well as different corresponding atmospheric attenuation.While the distribution of the scalars around a central value supports the approach we have used, this scaling remains as a potential source of systematic uncertainty within the density retrieval algorithm.
The g-factor scalars vary with local time and (to a lesser extent) throughout the year.This is likely caused in part because MSIS does not exactly characterize the short-term variation of the atmosphere with local time, season, or solar conditions; or it could be that the g-factor scaling approximation varies in accuracy with atmospheric and solar conditions.It is likely a combination of both.It is also worth noting that some of the variation at early and late LT could be due to more complicated solar geometry at high solar zenith angles.
Local time and coverage cycle are the best predictors of the optimal g-factor scalar, so the data shown in figures are median-averaged by LT bin for each coverage cycle and interpolated to create a lookup table like the one shown in Figure 6.For each ICON exposure, the forward model uses these tables to determine the scaling of the 83.4 and 108.5 nm g-factors.

Modeling the 61.6 nm Passband
The ICON EUV sensor was designed such that the 61.6 nm passband would spectrally isolate the O + 61.6 nm triplet (Sirk et al., 2017), as this feature is a measurement requirement for the ICON daytime ionosphere retrieval (Stephan et al., 2017).We use a similar method to that in Section 4.2 to determine the 83.4 nm g-factor scaling factor to be applied for the O + 61.6 nm emission.As for the 87.8 nm passband (see Section 4.2) we conduct a survey over every fifth day of 2020.For each 60 s exposure, we run the airglow model for 61.6 nm on the (scaled) MSIS densities as before and find the scaling for the 83.4 nm g-factor that best fits the 61.6 nm data.
The median value of the 61.6 nm g-factor scalar over the whole data set is 0.125, consistent with the value used in the ICON ionosphere algorithm (Stephan et al., 2017).Figure 7 shows the distribution of the normalized g-factor scalars.The distribution of the 61.6 nm g-factor scalars is tighter than that of the 87.8 nm scalars (Figure 5).This is largely because the blending of features in the latter is an additional source of systematic uncertainty as an error in the estimation of the scalar for one feature is coupled to the other.
There appears to be a secondary population in the right tail of the histogram in Figure 7. Unlike the 87.8 nm scalars, these 61.6 nm scalars do not show a strong dependence on LT but rather show an increase near the end of 2020.This coincides with an increase in solar F10.7 index; however, the exact cause has not been confirmed at this time.Rather than interpolating over LT as in Figure 6, the 61.6 nm g-factor scalar is estimated by interpolating over DOY.

Results
We ran the EUV neutral density inversion algorithm for all measurements between 6 and 18 LT during every fifth day of 2020, summing five consecutive ICON EUV exposures, now fixing the blending and scaling factors and iterating on our three density-relevant model parameters (see Section 3). Figure 8 shows the distribution of F10.7, O density, and N 2 density scaling factors from all retrievals.
MSIS 00 is used for the driver of the model and was used to determine the g-factor scalars used in the forward model.Thus, direct comparison to MSIS 00 does not provide validation.However, the trends in differences between the retrieval and MSIS 00 do indicate some overarching results.The distribution of F10.7 scalars shown in Figure 8 peaks below 1.0, indicating that the retrieval generally fits a cooler atmosphere to the EUV measurements compared to MSIS climatology.While the O scalar peaks around the training value of 0.85, indicating that the retrieval primarily adjusts O profile shape rather than scale (via F10.7), the N 2 scalar peaks lower, indicating that the shape and scale of N 2 profiles are adjusted.
The left pane of Figure 9 shows how retrieved O and N 2 densities differ from MSIS 00 as a function of altitude.In the case of each species, the scale height is smaller above 200 km, which is why both densities are increasingly lower than MSIS 00 at high altitudes.The low altitude data show that the retrieved O density agrees with the scaled down MSIS 00, while the retrieved N 2 base is generally lower.These observations are in agreement with our interpretation of the scalar distributions in Figure 8.

Validation
The Gaussian-like distribution of the retrieved scalars is evidence that the retrieval is physical, but this alone does not validate the EUV retrieval results.
In this section, we qualify and provide validation for the EUV results by comparing to MSIS 2.0, Swarm accelerometer-derived mass density, and ICON FUV and GOLD O/N 2 .TUMINELLO ET AL.

MSIS 2.0
While MSIS 2.0 is not an entirely independent model from the MSIS 00 model used to drive the retrieval, there is still information to be gained by comparing EUV retrieval measurements to MSIS 2.0.Namely, MSIS 2.0 has incorporated new data and is expected to be a major upgrade to MSIS 00 (Emmert et al., 2021).Figure 10 shows A qualitative interpretation of the scalar distribution in Figure 8 and the left panel of Figure 9 is that the retrieval has identified that N 2 needs to be scaled down by more than O does. Figure 10 shows that this is a change that was made in MSIS 2.0.However, MSIS 2.0 has N 2 decreased by about the assumed 15%, while reducing O by only 10%.
Figure 11 shows how the retrieved O and N 2 densities compare to MSIS 2.0 with local time at 250 km.The O difference has a local minimum around 13 LT and increases when approaching the terminator.The N 2 densities are the lowest relative to MSIS 2.0 in the evening.Of all our major sources of systematic uncertainty, the only one with a local time dependence is the g-factor scaling of the 86.5 and 87.8 nm lines.As the O 87.8 nm scalar increases with local time, it is possible the more of the brightness of the blended feature being attributed as originating from O leads to the decrease in retrieved N 2 density in the afternoon.In the case of O, it is unlikely that this is a result of systematic error in the retrieval; the 61.6 nm g-factor scalar is not LT dependent and 86.5 nm scalar has a LT dependence similar to the O density difference, which would likely increase estimated O density near the terminator.We are very possibly measuring a variation that is not present in MSIS 2.0.

Swarm Mass Density
As discussed in Section 1, recent neutral density measurements in the 150-400 km range of our retrieval are not plentiful.Aside from the FUV derived measurements discussed in Section 6.3, there are virtually no remote sensing measurements of neutrals in this region.Just above this region, however, it is possible to make in-situ measurements of mass density derived from accelerometer or GNSS measured satellite drag.In this section, we compare the EUV retrieved densities to such measurements from Swarm-C.
The Swarm mission utilizes three satellites: Swarm-A, Swarm-B, and Swarm-C, which orbit at altitudes of approximately 460, 530, and 470 km, respectively.The Swarm neutral mass density product is derived from spacecraft accelerometer data and high fidelity gas-surface interaction modeling (March et al., 2021).Of the two The MSIS 00-driven framework provides a convenient and effective means to extend our measurements up to Swarm altitudes, as we need only report the densities on an extended altitude grid using the same retrieved scalars.We define a conjunction as a pair of Swarm and ICON EUV observations that occur within 5°latitude, 5°longitude, and 30 min of each other.We consider the EUV retrieved mass density to be the sum of the contributions from O and N 2 only, which should account for 90%-95% of the mass density for a given exposure (according to MSIS 2.0).
The upper panel of Figure 12 shows the correlation between Swarm and EUV retrieval mass densities with correlation coefficient 0.93.For the majority of the data, the two measurements are very near equality.However, much of this strong correlation can be attributed to the fact that the retrieval has access to the MSIS climatology.
In order to separate the measurement from the driving model, we have plotted the Swarm and EUV measurements less the MSIS 00 mass density (of O and N 2 ) in the lower panel of Figure 12.The correlation is still strong and positive (r = 0.71), indicating that we really are measuring some of the variation that Swarm sees.It is noteworthy that the EUV retrieval is measuring higher densities than Swarm-C, despite the 5%-10% of mass density missing from species other than O and N 2 .

ICON FUV and GOLD ΣO/N 2
The traditional means of satellite remote sensing of the neutral thermosphere is through the measurement of the N 2 Lyman-Birge-Hopfield band and O 135.6 nm doublet in the far-UV (FUV) band of the electromagnetic spectrum.
Most often these dayglow features are observed on the Earth's disk in nadir-viewing geometry and used to calculate the column O/N 2 ratio (ΣO/N 2 ).This quantity is defined as the ratio of the column integrals of O and N 2 from the top of the thermosphere downwards until the N 2 column density reaches 10 17 cm 2 ; usually the altitude of this column base falls in the 130-140 km range.During 2020, this data product is available from GOLD and from ICON FUV.The retrieval algorithms for each of these products utilize a complex connection between the 135.6/LBH ratio and ΣO/N 2 dependent on factors such as solar EUV intensity and viewing geometry (Meier, 2021); we will forego this and calculate an EUV ΣO/N 2 by integrating our retrieved MSIS 00 atm according to Equation 8.This approach has previously been used and validated with data from TIMED-GUVI (Meier et al., 2015).
The ICON FUV instrument is a spectrographic imager with channels for LBH (near 157 nm) and 135.6 nm airglow.The imager measures 1D profiles on Earth's limb and disk (Mende et al., 2017).The ICON 2.4 data product is daytime ΣO/N 2 measured on the disk using a table lookup (Meier, 2021).During the daytime, ICON FUV has similar north-viewing geometry to that of the EUV instrument.The instruments are not synchronized, but do utilize the same 12 s cadence (Stephan et al., 2018).
We have combined the EUV data into 60 s exposures, so for each EUV retrieval we take the closest (in time) FUV ΣO/N 2 measurement as a conjunction if it occurred within 60 s.There is a latitude offset of the EUV and FUV measurements, since the former comes from the limb and the latter the disk, but this remains an acceptable comparison.
The GOLD instrument is a two channel UV spectrograph which observes Earth's disk in the FUV from a geostationary orbit around the mouth of the Amazon river.During the daytime, GOLD takes 12 min scans of the terrestrial disk, alternating between the northern and southern hemisphere.Column brightness of 135.6 nm and LBH 140.5-148.0nm are converted to ΣO/N 2 using a lookup table constructed using MSIS 00, the NRLEUV solar spectrum model, and the AURIC airglow model.The resolution of the GOLD ΣO/N 2 product is about 1.8°at nadir (Correira et al., 2021;Eastes et al., 2017).
We define a conjunction between GOLD and ICON EUV as a pair of observations that occur within 3.6°latitude, 3.6°longitude, and 1 hr of each other.If multiple GOLD pixels satisfy this criterion for a given EUV retrieval, we choose the pixel that minimizes the score defined as follows, where λ indicates longitude and ϕ latitude and Δ indicates the difference between the GOLD and EUV values: Figures 13 and 14 contain ΣO/N 2 scatter density plots of ICON EUV, ICON FUV, and GOLD, respectively.The ICON EUV retrieval correlates well with each of the other measurements, especially given the relatively loose conjunction definition used for ICON EUV and GOLD observations.Interestingly, ICON EUV retrieved ΣO/N 2 typically exceeds the ICON FUV value but is less than the GOLD value.The ICON EUV retrieval mostly differs from GOLD by a positive bias.The EUV ΣO/N 2 seems overly sensitive compared to the FUV, indicated by the slope greater than unity in the fit in Figure 13.
A difference between the EUV and FUV ΣO/N 2 is to be expected since EUV always observes a higher latitude; however, Figure 15 show that the difference does not reverse with season, which would occur if the discrepancy came solely from the latitude offset.Figure 15 provides a useful visualization of how the different ΣO/N 2 measurements vary throughout the year.Each source has its problem days, but for the most part they follow the same trend.Some notable exceptions occur, such as when ICON EUV and FUV trend opposite of GOLD around the beginning of March, or when FUV experiences dips in November and December that are not noticeable in GOLD or EUV.Lastly, we note that for one day, November 1, the EUV median ΣO/N 2 of 1.4 lies beyond the chart limits as presented here.

Discussion
We have discussed the development of an algorithm to measure thermospheric O and N 2 density using limb measurements of airglow at 61.6 and 87.8 nm from the ICON EUV instrument.This began with the development of an EUV airglow forward model, which predicts an EUV airglow profile at a given wavelength.Production is modeled using scaled g-factors from the well-established AURIC model, which form the basis for the existing ICON EUV ionosphere retrieval algorithm.Transport modeling considers pure absorption by neutrals.Neither the g-factors nor the cross sections are published with uncertainties, and each is a significant source of systematic uncertainty in the forward model and retrieval.
The forward model was used to perform a detailed study of the dayglow near 87.8 nm, which has seen little study outside of this work and its predecessor (Tuminello et al., 2022).The emissions are identified to be a blend of the O 87.8 nm line and a nearby N or N + feature.A survey was conducted over the ICON EUV measurements in 2020 to determine the best-fit g-factor scaling as a function of season and local time.The g-factors used represent another major source of systematic uncertainty.The lack of specific identification of the N/N + is a comparatively minor source of uncertainty, since absorption cross-sections vary a non-negligible but small amount over the range of possible wavelengths.
Discrete Inverse Theory is used to find values of three scalars applied to the MSIS 00 atmospheric model (input solar F10.7 flux and output O and N 2 density) which under our forward model best fit the ICON EUV measurements.Here the retrieval inherits systematic uncertainty in the ICON EUV data and from any systematic biases in the MSIS 00 model that would render the inverse method unable to converge to the true atmospheric state.These factors are relatively minor compared to the contributions to systematic uncertainty discussed in the previous paragraphs.
The retrieval output is physically plausible and compares reasonably to the MSIS 00 and MSIS 2.0 models.The most significant trend is that the retrieved atmospheres are cooler with smaller scale heights compared to the models.It is worth noting that thermospheric temperature modeling has not been significantly updated in MSIS since MSIS 00, meaning that MSIS 2.0 lacks temperature input from around the low solar minimum in 2020.
The retrieval compares favorably to coincident measurements from multiple satellite-based instruments.
Comparisons to mass densities derived from accelerometer-determined satellite drag on SWARM-C show that the two measurements show strong correlation before and after correcting for the EUV retrieval's access to the MSIS climatology.EUV derived measurements of ΣO/N 2 correlate strongly with FUV derived measurements on the same spacecraft, although the EUV values respond more sensitively to changes in atmospheric state.ICON EUV measurements of ΣO/N 2 are generally higher than that from ICON FUV; however, they are generally lower than what is seen by the GOLD instrument from GEO.Additionally, GOLD and EUV retrievals of ΣO/N 2 show similar sensitivity.All three measurements of ΣO/N 2 show comparable variation on daily and seasonal timescales.
The initial favorable comparison to other measurements from this retrieval algorithm encourages several areas for continuing development.One key aspect that remains is a rigorous characterization of retrieval uncertainty.New observations and further experimental results could also be used to improve the retrieval.High resolution cross-sections, such as the ones provided to us for atomic O would help to reduce systematic uncertainty, as would higher resolution measurements of EUV airglow that are sensitive in the 87.8 nm region.

Figure 1 .
Figure 1.A single forward model retrieval, showing the initial conditions and forward model output with the final fit and forward model output.The initial atmospheric state from scaled MSIS results in a peak that is too high, and the final fit shows that the peak has been lowered by decreasing the F10.7 input and scaling down O and N 2 densities.

Figure 2 .
Figure 2. Diagram of the forward model geometry.To find the contribution from a single point, the volume emission rate is calculated at that point (Equation2), and the optical depth is found by integrating neutral density along the red path to the satellite (Equation7).Another integration is performed to sum up these contributions along the entire line of sight (Equation6).

Figure 3 .
Figure 3.An example of an ICON EUV 87.8 nm exposure (black points) compared to forward model output of three N features (green lines).

Figure 4 .
Figure 4.An example of an ICON EUV 87.8 nm exposure (black) compared to a fit (purple) blending the N 86.5 nm triplet (red) with the O 87.8 nm doublet (blue).The legend indicates the factors by which the N + 108.5 nm and O + 83.4 g-factors are scaled, respectively.The shading represents the uncertainty in the linear fit due to the uncertainty in the data.The shape of the data is well-modeled, though the data appear to be more noisy at low tangent altitude.

Figure 5 .
Figure 5. Histograms indicating the distribution of the normalized (divided by the median values of 0.794 and 0.0708, respectively) g-factor scaling for the 86.5 and 87.8 nm features over the course of 2020.

Figure 7 .
Figure 7. Histogram indicating the distribution of the normalized (divided by the median value of 0.125) g-factor scaling for the O + 61.6 feature over the course of 2020.

Figure 8 .Figure 9 .
Figure 8. Distribution of F107, O, and N 2 scalars for all retrievals in 2020.The dotted line indicates the initial value, 1.0 for F10.7 and 0.85 for the density scalars.(These values of 1.0 and 0.85 are respectively used in order to determine the g-factors.)

Figure 10 .
Figure 10.Percent difference in O and N 2 density between MSIS 00 and MSIS 2.0.Error bars represent the upper and lower quartiles.The black dashed line represents equality, while the gray dashed line indicates 85% of MSIS 00 density, which was used to determine the g-factor scalars.

Figure 11 .
Figure 11.Scatter density at 250 km altitude of O (left) and N 2 (right) percent difference from MSIS 2.0 versus local time.Note that the y-axis scale differs between the two plots.

Figure 12 .
Figure 12.Swarm-C accelerometer derived density and EUV derived mass density (top) and the same densities minus the mass density of O and N 2 from MSIS 00 (bottom).The dashed blue lines and plotted text indicate the least squares fit and correlation coefficient of the data.The solid blue lines represents equality between the data.The small high density population is clearly driving the least squares fit in each case.

Figure 13 .
Figure 13.Scatter plot showing the correlation between coincident measurements of ICON FUV and EUV O/N 2 .The solid blue line and plotted text indicate the least squares fit and correlation coefficient of the data.The dashed red line represents equality between the data.

Figure 14 .
Figure 14.Scatter plot showing the correlation between coincident measurements of GOLD and ICON EUV O/N 2 .A small number of points fall outside the plot ranges for both GOLD and ICON EUV but are a sufficiently few points to be appropriately ignored.The solid blue line and plotted text indicate the least squares fit and correlation coefficient of the data.The dashed black line represents equality between the data.