# PEMEM Percentile: New Plasma Environment Specification Model for Surface Charging Risk Assessment

Stepan Dubyagin, Natalia Ganushkina, Angélica Sicard, Jean-Charles Matéo-vélez, L. Monnin, Daniel Heynderickx, Piers Jiggens, G. Deprez, F. Cipriani

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Stepan Dubyagin, Natalia Ganushkina, Angélica Sicard, Jean-Charles Matéo-vélez, L. Monnin, et al.: PEMEM Percentile: New Plasma Environment Specification Model for Surface Charging Risk Assessment. Journal of Geophysical Research Space Physics, 2024, 129 (2), pp.e2023JA032026. 10.1029/2023ja032026. [hal-04664798](https://hal.science/hal-04664798v1)

HAL Id: hal-04664798

https://hal.science/hal-04664798v1

Submitted on 30 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Space Physics

## METHOD

10.1029/2023JA032026

**Key Points:**

* The new specification model of cis-GEO near-equatorial plasma environment aimed at the surface charging risk assessment is presented

* The model covers the energy ranges 1–100 keV for electrons and 50 eV–50 keV for protons

* For a given orbital scenario, the model outputs flux percentiles that are expected to be encountered during the mission lifetime

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

S. Dubyagin and F. Cipriani,
stepan.dubyagin@fmi.fi;
fabrice.cipriani@esa.int

**Citation:**

Dubyagin, S., Ganushkina, N., Sicard, A., Matéo-Vélez, J.-C., Monnin, L., Heynderickx, D., et al. (2024). PEMEM percentile: New plasma environment specification model for surface charging risk assessment. *Journal of Geophysical Research: Space Physics*, 129, e2023JA032026. https://doi.org/10.1029/2023JA032026

Received 14 SEP 2023
Accepted 17 JAN 2024

©2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_3.png" /></figure>

## PEMEM Percentile: New Plasma Environment Specification Model for Surface Charging Risk Assessment

S. Dubyagin¹ ID, N. Ganushkina¹,² ID, A. Sicard³, J.-C. Matéo-Vélez³ ID, L. Monnin³, D. Heynderickx⁴ ID, P. Jiggens⁵ ID, G. Deprez⁵, and F. Cipriani⁵

¹Finnish Meteorological Institute, Helsinki, Finland, ²Climate and Space Science and Engineering Department, University of Michigan, Ann Arbor, MI, USA, ³ONERA - The French Aerospace Lab, Toulouse, France, ⁴DH Consultancy BV, Leuven, Belgium, ⁵ESTEC, ESA, Noordwijk, The Netherlands

**Abstract** The Plasma Environment Modeling in the Earth's Magnetosphere (PEMEM) is a European Space Agency activity supporting the development of a new specification model for the spacecraft surface charging risk assessment. This paper presents a description of the basic model version: the PEMEM *percentile* model. The model is intended to be used for space missions with near-equatorial orbits. The model is based on the Van Allen Probes particle measurements inside the geostationary orbit. The model's primary input is a planned spacecraft trajectory. It outputs statistical characteristics of the plasma environment which are expected to be encountered during a mission lifetime. These characteristics include differential electron and proton flux percentiles for a set of energies (percentile spectra), and percentiles of the integrated electron flux. The model covers the energy range of 1–100 keV for electrons and 40 eV–51 keV for protons. Since extreme spacecraft charging usually occurs in the eclipse, the same characteristics can be separately output for the periods when the spacecraft is shadowed by the Earth.

### 1. Introduction

Among the environmental factors threatening the success of unmanned space missions, high-energy particle radiation is arguably the most widely recognized by the general public, since the particles with energies greater than ~0.1 MeV can penetrate the spacecraft's outer layers and damage the onboard computer. However, while lower-energy plasma particles may not evoke the same level of fear as high-energy radiation, they should not be underestimated, as under certain conditions they can be deposited on the spacecraft's surface leading to surface charging that also presents a potentially deadly threat to the spacecraft.

The spacecraft surface charging is a phenomenon when the spacecraft charges negatively with respect to the ambient plasma due to a combination of the photoelectric effect, and interaction with incident particle fluxes (Davis et al., 2008; Garrett, 1981; Whipple, 1981). Since, under certain conditions, the spacecraft's negative potential can reach as low as a few thousand volts, this process can distort the spacecraft measurements and also lead to electrostatic discharge between the elements of the spacecraft leading to surface damage and the electronic system malfunctions. Since the first detection of this phenomenon during the early era of space flights (DeForest, 1973), numerous anomalies, failures, and even losses of the missions were attributed to surface charging (Koons et al., 1998; Leach, 1995).

The spacecraft charging potential is determined by the balance of electric currents created by particles emitted from the surface and incident upon it (Garrett & Whittlesey, 2012). In turn, these currents are influenced by the availability of charge carriers and the spacecraft's potential: the negatively charged spacecraft repels low-energy electrons, diminishing the incident electron flux, with the repulsion strength increasing as the potential becomes lower. When the spacecraft is in sunlight, its surface emits low-energy photoelectrons, resulting in a positive spacecraft charge. This positive potential typically does not exceed a couple of tens of volts. However, if the incident electron flux from ambient plasma is very high, the corresponding current can surpass the photoelectron current, leading to negative spacecraft charging. In sunlight, the negative potential rarely drops below a hundred volts. Conversely, if the spacecraft is in shadow, there is no photoelectric current, and the ion flux from the ambient plasma becomes the sole source of charge carriers capable of balancing the electron flux current. However, the ion flux is orders of magnitude lower than the electron flux. As a result, a spacecraft eclipsed by the Earth continues charging until the electric field's repelling effect reduces the electron flux enough to balance the ion flux. Occasionally, this condition is met only after the spacecraft potential becomes several thousand volts
negative. To further complicate the subject, there is one more current created by secondary electrons emitted when an energetic particle strikes the surface. While many of these processes depend on surface materials and spacecraft geometry and cannot be easily quantified, empirical studies have shown that the flux of keV electrons is a primary controlling factor (Davis et al., 2008; Ganushkina et al., 2021; Matéo-Vélez et al., 2016, 2018, 2019; Sarno-Smith et al., 2016; Thomsen et al., 2013).

A few engineering 3D modeling systems numerically simulating the charging of a spacecraft for a given spacecraft design and ambient plasma environment have been developed, including European Space Plasma Interaction System (SPIS) (Sarrailh et al., 2015), https://www.spis.org/, and NASA Charging Analyzer Program NASCAP-2K (Mandell et al., 2006). At the same time, significant efforts were made toward modeling the Earth's magnetosphere plasma environment (Boscher et al., 2003; Ginet et al., 2013; Sawyer & Vette, 1976; Sicard et al., 2018; Sicard-Piet et al., 2006, 2008; Vette, 1991). Among these models, the family of NASA's models AE8, AP8, AE9, AP9 (Ginet et al., 2013; Sawyer & Vette, 1976; Vette, 1991) targeted the high energy particles responsible for the internal charging. As a result, these models lack the MLT dependence which can be ignored for MeV particles. NASA's SPM model (Ginet et al., 2013), designed to address the surface charging and keV energies, still lacks MLT dependence. On the other hand, these models provide reach functionality and can output a variety of statistical characteristics. On the contrary, the recent Global Radiation Earth ENVIronment model (Sicard et al., 2018, 2019) includes MLT dependence and is capable of modeling electron fluxes from 0.19 keV and higher but it only outputs mean flux. There exists another class of empirical models that output average environmental characteristics based on solar wind parameters or geomagnetic activity indices (e.g., Denton et al., 2015; Sillanpää et al., 2017; Stepanov et al., 2021). A subclass within this category employs a machine learning approach for determining plasmaspheric density (Bortnik et al., 2016; Chu et al., 2017; Zhelavskaya et al., 2017; Zhou et al., 2022), and for predicting electron fluxes (Boynton et al., 2016, 2019; Simms et al., 2022, 2023; Swiger et al., 2022). However, these models pose challenges for practical application in risk assessment, as neither the solar wind parameters nor the activity indices are known in advance for future missions.

To fill the aforementioned gaps, a suit of two specification models targeting the plasma environment in the energy range critical for surface charging was developed under ESA contract. The models are based on Van Allen Probes flux measurements but the other data sets can be added later. In this paper, we present the basic model, Plasma Environment Modeling in the Earth's Magnetosphere (PEMEM) *percentile*, providing the percentiles of the electron and proton fluxes for a given orbital scenario. A more advanced PEMEM AEdMC (AE index-driven Monte Carlo) model, which provides the worst-case electron flux for a given confidence level and is geomagnetic activity-dependent, will be presented in a separate paper.

## 2. Instrumentation and Data Preparation

Energetic Particle, Composition, and Thermal Plasma Suite (ECT) onboard Van Allen Probes include three instruments, but in this paper, we use the data of the two covering the lower energies: Helium, Oxygen, Proton, and Electron (HOPE) (Funsten et al., 2013) and Magnetic Electron Ion Spectrometer (MagEIS) (Blake et al., 2013). We use the release 4 data covering the years 2013–2018.

The HOPE detector measures electrons and protons in alternate spin cadence, thus, the data for given species are provided at 23-s resolution while the accumulation cycle lasts 12 s. We use the precomputed omnidirectional electron flux data in the apogee mode which are provided along with the pitch-angle distributions. There are 72 electron energy channels covering the energy range from ~15 eV up to ~52 keV. For some periods, there are data only for a fraction of these channels. There are data quality flags that indicate potential issues with the data. The measurements with nonzero flags are discarded.

The MagEIS L3 data are provided at 12-s resolution. The energy range is from 32 keV to several MeV. We use the data up to 150 keV to ensure the interpolation of the flux at the model's upper limit of ~100 keV. For some periods, the electron flux was corrected for penetrating radiation background and we use these data when available, otherwise, we use uncorrected flux data. There is no precomputed omnidirectional flux in the MagEIS level 3 data files. For each of 25 energy channels, the angular distribution is provided at 11 pitch-angles spanning 0–180° range. To convert from the pitch angle resolved data to the omnidirectional flux, we perform the numerical integration assuming the gyrotropy of the angular distribution.
$$J(E) = \frac{1}{4\pi} \int j(E, \alpha) \sin \alpha d\alpha \quad (1)$$

Here, J is the omnidirectional flux, E is the energy, α is the pitch angle, and j is the unidirectional flux.

A significant portion of the MagEIS data set lacks data for the field-aligned parts of the pitch-angle distributions. When possible, we fill these gaps assuming the symmetry: $j(E, \alpha) = j(E, 180^{\circ} - \alpha)$. However, this procedure still leaves ~20% of distributions with no data for the field-aligned sectors, with the exact percentage varying depending on energy, year, and probe. If the number of problematic sectors exceeds two, the distribution is discarded (constituting 2%–3% of measurements). Otherwise, we extrapolate the flux in these sectors using the flux measurements from three adjacent sectors, assuming a linear dependence of log j on α. We do not allow the extrapolated values to exceed the maximum measured flux (for the given pitch-angle distribution) by more than a factor of 2. If they do, we set their values to this threshold value. It should be noted that a factor of 2 threshold was chosen ad hoc.

Boyd et al. (2019) reported that a disagreement between the MAGEIS and HOPE electron measurements was found for some dates. We used the correction factors provided by the authors for these dates (see Table S1 in Supporting Information of Boyd et al. (2019)). In addition, since the MagEIS measurements can be potentially contaminated in the inner zone (Claudepierre et al., 2019), we use only measurements for L > 3.5.

Following Boyd et al. (2019), we discard the electron measurements of those HOPE channels where/when there are fewer than 10 counts. The only exception was the periods when MAGEIS measurements were not available. It should be noted that since the surface charging is caused by the high electron fluxes, the accuracy of low electron fluxes is not so critical and we prefer to keep the data record with low counts in the data set even if it is not quite accurate.

For the periods when the spacecraft were charged negatively, the combined electron and proton spectra from both HOPE and MAGEIS were corrected for spacecraft potential. The potential has been estimated with the ion peak method (e.g., Davis et al., 2008) from the ion spectra measured by the HOPE instrument following the same approach as in Matéo-Vélez et al. (2018) and Matéo-Vélez et al. (2019). The description of the correction procedure can be found in Text S1 in Supporting Information S1.

To regularize the spectra measured on two probes in different years to a uniform format, we interpolated the spectra to a fixed set of energies using log-log linear interpolation. The energies spanned the ranges 0.9–106 keV for electrons and 0.04–51 keV for protons with a factor of 1.3 increment. This results in 19 and 28 energies for electrons and protons, respectively.

Finally, since the surface charging is controlled by the integrated rather than differential electron flux (Thomsen et al., 2013), we also calculated the integrated electron flux for a set of the lower integration limits represented by the first 13 energies of the regularized spectra covering the range 0.9–22 keV. The upper integration limit was fixed at 100 keV. We use McIlwain L and MLT provided with Van Allen Probes data. These coordinates were calculated using Olson and Pfitzer (1977) (OP77Quiet) model.

### 3. Empirical Statistical Distributions

We start our model development by analyzing the statistical properties of the electron flux measurements. In contrast to MeV electrons, the drift of the lower energy electrons is slow enough for the electric field and loss mechanisms to affect the distribution on the scale of the azimuthal drift period. As a result, the electron fluxes reveal a prominent MLT dependence that is getting stronger during enhanced geomagnetic activity (e.g., Fernandes et al., 2017). Figure 1 shows L, MLT maps of the 50th and 99th percentiles of the electron flux for 1, 10, and 100 keV energies. The percentiles were calculated for the data subsets in 0.25 $R_E \times 1$ hr bins. Note that the color scales are different for 50th and 99th percentiles with the red color corresponding to roughly one order of magnitude higher flux for the 99th percentile in comparison to the 50th percentile for the same energy. The MLT asymmetry can be seen for all plots except for Figure 1e. It is seen that the asymmetry is strong even for 100 keV electrons in Figure 1f. The highest fluxes in Figures 1b and 1d (1, 10 keV) are observed in 21–24–06 MLT sector. It is exactly the sector where the surface charging risk is the highest (Matéo-Vélez et al., 2018, 2019; Sarno-Smith et al., 2016). For 100 keV, the sector of highest fluxes is shifted eastward to 1–14 MLT. These figures illustrate the
<figure><img src="image_4.png" /><figcaption>**Figure 1.** 50th and 99th percentiles of the electron flux calculated for spatial bins for three energies shown as a function of L and MLT. Flux units are s⁻¹ cm⁻² ster⁻¹ keV⁻¹.</figcaption></figure>

importance of the inclusion of MLT dependence in the models for surface charging risk assessment. Although the region of highest fluxes for the 99th percentile is closer to the Earth in comparison to that for the 50th percentile (especially obvious for 10 keV), it is still located outside $L = 4$. As was mentioned in Section 2, the electron data can be affected by the penetrating radiation in the inner zone ($L < 3$), but Figure 1 demonstrates that the surface charging risk is low in this region. It is also confirmed by the direct surface charging observations (Matéo-Vélez et al., 2019; Sarno-Smith et al., 2016). For this reason, we exclude this region from the analysis, assuming that the surface charging risk is low there. Similar flux percentile maps for protons can be found in Supporting Information, Figure S2 in Supporting Information S1.

We further investigate the statistical properties of the electron flux measurements by plotting statistical distributions for 10 keV electron flux for selected spatial bins. The vertical axes in Figures 2a–2d show the exceedance probability: the probability that a value on the horizontal axis will be exceeded. Here, it is calculated as a fraction of the flux measurements exceeding the value given on the horizontal axes. In comparison to Figure 1, we use larger bin size: 0.25 $R_E \times 2$ hr. In Figures 2a and 2b, we demonstrate the MLT dependence (different sectors shown by different colors) for two L parameter intervals at the inner and outer boundaries of the region, respectively. In Figure 2a, the distributions in 10–20 MLT sector (orange—yellow—green curves) reveal the presence of only one statistical population (unimodal distributions; the curve's slope is getting progressively steeper from left to right). It will be later shown that it is well described by a single log-normal distribution. In the other MLT sectors (red—blue), the high-flux tails of the distributions reveal some traces of the additional statistical population (bi-modal distribution): flattening of the curve or even shoulder-like shape.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Exceedance probability functions of the 10 keV electron differential flux for fixed L range and different MLTs (a, b) and fixed MLT sectors and different L values (c, d).</figcaption></figure>

In contrast, for the outer bins (Figure 2b), the curves having shoulder-like shapes (corresponding to bi-modal distributions) can be seen for almost all but 12–18 MLT sectors. This is in agreement with Denton et al. (2015) results who reported bi-modal electron flux distributions at geosynchronous orbit. Similarly to Figure 2a, the additional population is manifested most prominently in 00–04 MLT sector. Comparing the two figures, it can be seen that almost half of the samples can be originated from the additional populations in Figure 2b while their share is less than 1% in Figure 2a.

In Figures 2c and 2d we investigate the variation of the distributions over radial distance for two MLT sectors revealing the most different distribution shapes, 0–2 and 16–18 MLT respectively. It is seen that at the post-midnight MLT sector (Figure 2c), the contribution of the second statistical population gets progressively stronger with distance. On the other hand, for the pre-dusk sector (Figure 2d), there is almost no dependence on L, and the presence of a second population can be barely traced for L values greater than ~5. Note also that there is very subtle dependence on L for the last three L parameter ranges: 5.5–6. This suggests that the model results for a given MLT sector can be extended outward.

Comparing Figures 2a and 2d, it can be seen that the distributions for exceedance probability greater than $10^{-2}$ are almost identical for all MLTs and L. This suggests that there is a background population and there is some process operating sporadically that injects a secondary population. The probability of observing this population is lower in the inner part of the region and in the dusk sector. This secondary distribution is never observed at $L = 3.5$ and 16–18 MLT. All these results imply that this process is a substorm injecting 10 keV electrons near midnight, with injection depth typically not reaching $L = 3.5$, and these electrons drift eastward being partially lost to the atmosphere so that only a small fraction of them reaches the dusk sector. We will not discuss this further since the details of the physical process forming the statistical distributions are not important for the model design.

In order to investigate the variability of the electron flux, we calculated the correlation of the electron flux logarithm with its delayed values. It should be noted that it is impossible to separate the spatial and temporal variation in the time series from a moving spacecraft. For this reason, we calculate autocorrelation only for the outer bins at $L = 5.75–6$, where a probe moves slower near its apogee. We selected all measurements in a given
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Correlation of the logarithmic flux with its lagged values for two spatial bins.</figcaption></figure>

bin and then calculated the correlation for different delays irrespective of whether delayed measurements belonged to that bin or not. Figures 3a and 3b show the autocorrelation for RBSP-B for two spatial bins at 0–2 and 14–16 MLT, respectively. The colors indicate three different energies. It can be seen that for all energies and for both bins the correlation is still higher than 0.6 even for 1 hr delay demonstrating that on average the electron flux variation is slow. The variation is slower for 100 keV energy and for the post-noon bin.

However, the curves in Figures 3a and 3b were obtained using all measurements irrespective of the geomagnetic activity. It is natural to expect that the electron flux is more volatile during active periods. To investigate this, we selected the measurements with an electron flux higher than 95th percentile for a given bin. Then we calculated autocorrelation with delayed measurements irrespective of whether the electron flux was higher than the percentile or not. The results are shown in Figures 3c and 3d. Figure 3c shows that for 0–2 MLT bin for 1, 10 keV energies the correlation drops to ~0.5 for 5-min delay and to nearly zero for 20-min delay. The autocorrelation scale is somewhat longer for 100 keV. Although the autocorrelation scales are comparable or longer in Figure 3d for 14–16 MLT bin, they are still much shorter than those in Figures 3a and 3b. This is an indication that keV electron flux can vary even on the sub-minute scale for extreme events. For this reason, we do not use averaged electron flux data sets like those developed by Boyd et al. (2019, 2021) but use the data at the original resolution. On the other hand, the usage of all data in a given bin also leads to at least two problems. Statistical distributions are supposed to describe the random variables and the randomness implies that each next sample is independent of the previous one. Figure 3 demonstrates that it is not the case even during active periods. The interval between the measurements should be at least 20 min during active periods and more than 1 hr during quiet periods. Another problem is illustrated in Figure 4. The black crosses show the empirical exceedance probability distribution derived from all measurements in a given bin. For measurements sorted in descending order, the exceedance probability for a given flux value $J_i$ can be calculated as $i/N$ where i is the measurement's index and N is the data
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Exceedance probability distribution for all measurements in the bin (black) and for those selected from the individual bin crossings (red). The vertical dashed line marks the flux value above which all measurements correspond to one and the same bin crossing.</figcaption></figure>

set size. Thus, the black symbols in Figure 4 correspond to individual measurements. All symbols to the right from the vertical dashed line correspond to a single event (a single bin crossing by the probe trajectory). Obviously, the flux did not vary rapidly during this event and as a result, the symbols form an almost vertical drop at the distribution's tail. Such a feature is very common and can be seen for many bins and energies. If such a curve is fit by an analytical distribution, this feature can lead to an underestimation of the flux at the distribution tail. In addition, since the duration of the individual bin crossing by a spacecraft trajectory varies from one measurement (if the corner of the bin is crossed) and up to ~1 hr (if a spacecraft moves azimuthally along the bin near the apogee), different substorms can be represented with different statistical weights. To avoid these undesired effects, we resampled our data set randomly picking up a single measurement from each bin crossing. It is equivalent to saying that for each bin one orbital period is represented by a single randomly selected measurement. The ~9 hr Van Allen Probes' orbital period almost definitely provides the independence of the sample values on prehistory. The empirical distribution obtained by this method is shown in Figure 4 by red cross symbols. This approach also makes the data coverage more uniform: Figures 5a and 5b show the number of measurements in a bin and the number of bin crossings, respectively. It is seen that there is no radial dependence between L = 3.5 and L = 6 in Figure 5b.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) The number of measurements in 2h × 0.25 Re bins. (b) The number of bin crossings.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6. Explanatory sketch. Not to scale.</figcaption></figure>

## 4. Model Philosophy

The results of Section 3 demonstrate that in order to adequately predict statistical characteristics of the particle fluxes for a variety of orbits, the model should take into account a number of requirements. First of all, the model should include at least dependencies on L and MLT (we cannot assess the importance of the dependence on $B/B_{eq}$ using the measurements in the near-equatorial region). For this reason, we model the statistical distributions separately for each individual spatial bin (all bins have the 0.25 $R_E \times 2$ hr size). Second, although the RBSP data set is rather large and covers the descending phase of the solar cycle including maximum and minimum, the user might want to get the estimates for the extreme events and it is desirable for the model to be able to extrapolate beyond the limits represented in the parent data set. To address this requirement, we fit the empirical distributions with the analytical distribution functions assuming that they provide a reasonable extrapolation of the distribution tails. Although it was shown that the empirical distributions generally can not be described with an unimodal distribution, using too many distributions results in overfitting and unrealistic extrapolation. For this reason, we use only uni- or bi-modal distributions for the electron flux. Note also that the autocorrelation time scale can be as long as an hour during quiet periods. On the other hand, it can be as short as a minute during disturbed time. Since surface charging is a very fast process

and can develop on a time scale of milliseconds, we use the original 23-s resolution data resampled as described in Section 3.

The model philosophy is schematically illustrated in Figure 6 and can be outlined as follows: The mission designer provides the coordinates of the expected spacecraft trajectory for the modeled period (blue curves in Figure 6). Using these data, the total time ($\Delta t_i$), that the spacecraft will spend in each bin can be evaluated as well as the time ($\Delta T$), that the spacecraft will spend outside the region of the surface charging risk (model validity region). The exceedance probability distribution for the differential or integrated flux encountered during the modeled period can be obtained as a weighted sum of distributions for all bins where weights are proportional to the durations the spacecraft will spend in the corresponding bins:

$$P_k(J) = \sum_i w_i \cdot P_{ki}(J) \quad (2)$$

Here, $P_{ki}(J)$ is the exceedance probability function for the $k$th energy and $i$th spatial bin. The summation is done over all model's spatial bins (index $i$). The index $k$ corresponds to the energy of the differential flux or the lower integration limit for the integrated flux. The weight coefficient for a given bin is calculated as

$$w_i = \frac{\Delta t_i}{\Delta T + \sum_i \Delta t_i}, \quad (3)$$

where the denominator represents a total modeled time.

Recalling the relation between the percentage corresponding to $N$th percentile and the exceedance probability $P = 1 - N/100$, the percentiles of the flux for $k$th energy channel can be calculated by inverting Equation 2.

# 5. Detailed Model Design

## 5.1. Time-In-Bin Calculation

The total time spent by a spacecraft in a given spatial bin is calculated from the L, MLT coordinates of the spacecraft trajectory provided by the user at finite time resolution. The coordinates are converted to 2D cartesian vector using the following equations:
$$ \vec{x} = \begin{pmatrix} L \cdot \cos(\pi MLT/12h) \\ L \cdot \sin(\pi MLT/12h) \end{pmatrix} \qquad (4) $$

The spacecraft motion between two adjacent points of the trajectory is approximated using linear interpolation:

$$ \vec{x}(t) = \vec{x}_1 + (\vec{x}_2 - \vec{x}_1) \frac{(t - t_1)}{(t_2 - t_1)}, \quad t \in [t_1, t_2] \qquad (5) $$

Here, *t* is the time, and indices 1 and 2 correspond to the two consecutive points of the trajectory. Using Equation 5, the times of the trajectory crossings with the bin boundaries can be found and the time in the bin can be calculated.

Since the same environment leads to much higher surface charging if the spacecraft is in the shadow (eclipsed by the Earth), we have added an option to separately output the results for the eclipse periods. To identify the eclipse periods, the GSE coordinates and linear approximation for the motion along the trajectory segment are used:

$$ \vec{r}(t) = \vec{r}_1 + (\vec{r}_2 - \vec{r}_1) \frac{(t - t_1)}{(t_2 - t_1)}, \quad t \in [t_1, t_2] \qquad (6) $$

Here, $\vec{r}$ is the 3D vector in GSE coordinates.

The eclipse condition is defined as:

$$ \sqrt{Y_{\text{GSE}}^2 + Z_{\text{GSE}}^2} < 6371.2\text{km} \qquad (7) \\ X_{\text{GSE}} < 0 $$

Equations 6 and 7 are used to find the times when the spacecraft enters and exits the eclipse. The total time for a given bin is calculated by summing the contributions from all trajectory segments. For those bins where the spacecraft spends part of its time in eclipse, two total times are calculated: in eclipse and in sunlight.

## 5.2. Electron Differential Flux Distribution Fitting

The lognormal distribution has been used by many authors to describe the particle fluxes in near-Earth space (Feynman et al., 1993; Hardy et al., 2008) and we will show that it works well for the electron and proton fluxes in our energy range. The lognormal distribution is defined as

$$ A(J) = \frac{1}{J\sqrt{2\pi}\sigma} \exp\left(-\frac{(\log(J) - \mu)^2}{2\sigma^2}\right), \qquad (8) $$

where *A* is the probability density function of a variable *J* (in this paper, *J* is the flux), and *μ* and *σ* are the free parameters or of the distribution. If a combination of two or more lognormal distributions is used to describe the data, weight coefficients and the normalization condition give 3 · *M* − 1 free parameters, where *M* is the number of lognormal distributions.

$$ A^{\text{mult}}(J) = \sum_{m=1,M} \omega_m \cdot A(J, \mu_m, \sigma_m) \\ \sum_{m=1,M} \omega_m = 1 \qquad (9) $$

Here, $\omega_m$ are the weight coefficients.

For each spatial bin and energy channel, we perform two fittings using a single and a sum of the two lognormal distributions getting values of *μ*, *σ*, and *ω* parameters. There are several ways of fitting the analytical statistical distribution to the data set. However, since the main goal of the modeling is to describe the high-flux tail of the
statistical distribution, which is of interest for the surface charging risk assessment, we would like to give more weight to the high-flux tail of the distribution. It is not possible to do this using the maximum likelihood fitting technique. For this reason, we find the parameters of the analytical distribution minimizing the absolute deviation between the empirical and analytical exceedance probability functions in log-scale:

$$ \Delta = \sum_{i} W_{i} \cdot \left| \log \left( \frac{P_{i}^{\text{empir}}}{P_{i}^{\text{analytic}}} (J_{i}) \right) \right| \quad (10) $$

Here, $P_{i}^{\text{empir}}$ is the empirical probability of exceedance (corresponding to the ith flux measurement: $J_{i}$), $P_{i}^{\text{analytic}}(J_{i})$ is the probability of exceedance given by the analytical distribution at $J_{i}$, and $W_{i}$ is the weight coefficient. The summation could be done over all measurements in the data set (see the details of the empirical exceedance probability computation in Section 3) but for the sake of efficiency, we only do this for the highest 50 flux values (distribution tail), while the rest of the distribution is sampled so that there are 100 points evenly spaced in log-scale. Defining the weight coefficients $W_{i}$ we take into account two aspects: First, we account for the density of points per log J interval to ensure the equal significance of the different parts of the empirical exceedance probability curve. Second, since the empirical estimations of the probability of exceedance became less accurate for the distribution tail, we set the weight coefficient for the last 10 points (highest flux) to be an inverse of the expected empirical distribution error and we fix the error value for the rest of the points as shown below:

$$ W_{i} = \begin{cases} \frac{1}{\rho_{i}\delta_{i}}, & \text{if } i \ge N - 10 \\ \frac{1}{\rho_{i}\delta_{i}=10}, & \text{otherwise} \end{cases} \quad (11) $$

Here $\rho$ is the measurement density (per log J) and $\delta$ is the error of exceedance probability estimation for a given measurement. It is assumed that the flux measurements are sorted in ascending order so $i > N - 10$ corresponds to the 10 highest flux measurements.

Figures 7a and 7b show examples of the empirical exceedance probability distributions and their approximations for two spatial bins. Blue and red solid curves correspond to single- and two-lognormal approximations. It can be seen that the distribution in Figure 7a can be adequately described by a single distribution but two distributions are needed in Figure 7b. Thus, we use a single lognormal distribution if it provides a good approximation and we use two-lognormal distributions otherwise.

### 5.3. Electron Integrated Flux Distribution Fitting

An identical approach was used for fitting the distributions of the integrated electron flux. It should be noted that our data set includes the integrated electron fluxes calculated for 13 lower integration limits in the range between 1 and 20 keV (the upper limit was always set to 100 keV). Thus, we model the integrated flux for different integration limits independently. In other words, the lower integration limit serves as an energy channel in Equation 2 (index k).

### 5.4. Proton Differential Flux Distribution Fitting

The philosophy of the proton percentile flux model is similar to that of the electron model. However, there are a few differences. There are no distributions of the integrated flux for protons. The energy range covers lower energies: 42 eV–51 keV in 28 energy channels. The most important difference: A different distribution fitting technique is used. Similar to the electron flux model, we used multiple log-normal distributions. However, in contrast to electrons, we want both low- and high-flux parts of distribution to be approximated with comparable accuracy. For this reason, we used the Expectation Maximization method (Ghojogh et al., 2020; Hastie et al., 2009) to fit successively 1-, 2-, 3-, and 4-modal distributions. Although the usage of multiple lognormal distributions can potentially lead to overfitting and large errors in the case of extrapolation, we have chosen this approach because the protons are less important for the surface charging risk assessment, and the accuracy of
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The empirical distributions (cross symbols) fitted by log-normal distributions (solid curves) for the (a, b) electron and (c, d) proton differential flux for 10 keV energy. Vertical axes show the (a–c) exceedance probability and (d) cumulative probability. Dashed curves in panels (a, b) show the expected range of the empirical distribution error for 95% confidence.</figcaption></figure>

extrapolation can be traded off for accurate representation within the percentile range covered by the Van Allen Probes' measurements (0.1%–99.9%).

Figures 7c and 7d show the exceedance probability distribution and cumulative distribution of the differential flux for 10 keV protons. Although the usage of four populations does not give an obvious improvement in the fit quality in comparison to the three-population fit, we use the highest number of populations as long as the fitting procedure has converged.

## 5.5. The Model Extension to Beyond Van Allen Probes Apogees

Note that the empirical distributions and their analytical fits are available only within Van Allen Probes apogees distances. We set a limit of 900 data points for a given MLT-L bin above which the fit is considered to be valid. Figure 5b shows that the outer boundary encompassing the region that meets this criterion is located at $L \approx 6.5$ on the nightside and at $L \approx 6$ on the dayside. This region does not fully cover the geostationary orbit and it would be desirable to extend our model to outer distances. Figures 2c and 2d demonstrate that electron flux distributions are almost constant for $L > 5.5$. This can be considered as a justification for using the outermost model distributions to extend the model applicability region farther outward. Note that no real extrapolation of the distribution parameters is used, we just use the distributions for the same MLT sector for $L$ values beyond $L = 6-6.5$. Note that the MLT dependence is still preserved because such extension is done separately for each MLT sector.

# 6. Model Application Example

In this section, we present the example of the model application for the spacecraft with electric orbit raising. The spacecraft is inserted into orbit with low perigee and ~12 $R_E$ apogee. The electric propulsion is then used to circularize the orbit at the geostationary orbit. Orbit coordinates for 250 days were provided by ONERA at 500-s resolution and are shown in Figure 8a.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** (a) The orbit for the mission with EOR in L MLT coordinates shown at 25-min resolution. (b) The amount of time spent in a bin.</figcaption></figure>

The model's outer limit was set to L = 8 (see Section 5.5). Figure 8b shows the total time the spacecraft spends in each bin. It should be noted that of 6,021 hr of the total time only 3,289 hr were spent inside the modeled region 3.5 < L < 8. For demonstration purposes, the time outside the modeled region was set to zero (ΔT = 0 in Equation 2). The times-in-bin from Figure 8b were used to calculate the weights (Equation 2) and then Equation 1 was used to derive the statistical distributions of the electron and proton fluxes for the given orbit. These statistical distributions are shown in Figures 9a–9c. Different colors correspond to different energies. Figure 9a shows the electron differential flux exceedance probability. Note that the presence of two statistical populations can be seen. Figures 9b and 9c show the exceedance probability and cumulative probability distributions for the proton differential flux, respectively. If we take a “horizontal slice” of these distributions at some probability value, we get a spectrum (flux vs. energy) for a given percentile. Such spectra for a few percentiles are shown in Figures 9d–9f. It is seen that the proton spectrum for the 99.999% percentile is more noisy than that for electrons. Since 99.999% percentile is in the extrapolation region (there are ~2,000 measurements in the individual spatial bin, see Figure 5b), we attribute this noise to 3-, 4-population fits for the protons and resulting overfitting of the distributions.

The model also outputs the distribution of the integrated electron flux. This distribution is shown in Figure 10 by a black curve. However, we can also integrate the percentile spectra like those shown in Figures 9d–9f. Since there is a direct correspondence between exceedance probability and percentiles, we can overplot such values on Figure 10. The red curve in Figure 10 shows the integrated electron flux distribution obtained from the percentile spectra. It is seen that this distribution can exceed the modeled one by a factor of ~1.5. This difference results from a percentile spectrum's nature. Since the percentile spectra are essentially the differential flux percentiles at a set of energies computed independently, these are upper envelope spectra. On the other hand, the real spectrum corresponding to the *i*th percentile of the integrated flux is equal to the *i*th percentile spectrum only if the differential fluxes at all energies vary in coherence. It is hardly the case. Since the power-law or even exponential decrease of the typical electron spectrum tail, and given that the realistic electron temperatures are below 10 keV, the ~10 keV electrons are major contributors to the 10–100 keV integrated flux. This means that the real spectrum corresponding to the *i*th percentile of the integrated flux is expected to be close to the *i*th percentile spectrum at 10 keV and is expected to lie below it at higher energies.

Apart from the percentile spectra shown in Figures 9d–9f and the integrated electron flux exceedance probability distribution shown in Figure 10, the model can output the time of exceedance for the user-specified integrated electron flux threshold value. This time is obtained by multiplying the exceedance probability corresponding to the threshold flux value by the total mission time. The model also separately outputs all data products for the part of the orbit when the spacecraft is eclipsed by the Earth.

It is important to emphasize that there is no precise relationship between the percentile spectra or the integrated flux distribution and the resultant spacecraft potential, let alone its spatial distribution. To obtain this information, spacecraft surface charging engineering simulations should be used, and the percentile spectra should be used as an input representing the extreme environment. In its simplest form, this step includes the development of the
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The (a) electron and (b, c) proton flux distributions for different energies. (d–f) Corresponding percentile spectra.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Distribution of the integrated (10–100 keV) electron flux modeled independently (black) and calculated from percentile spectra (red).</figcaption></figure>

CAD model of the spacecraft (specifying the spacecraft geometry and surface materials), and specifying the ambient plasma environment. The examples of the SPIS simulation application to the spacecraft in the GEO environment can be seen in Matéo-Vélez et al. (2014) and Nakamura et al. (2018). Since the preferred format for the plasma environment description in such simulations is the parameters of the standard phase space density distribution functions rather than spectra at discrete energies, the PEMEM percentile spectra are approximated by 1–4 Maxwellian distributions and the corresponding parameters (densities and temperatures) are output along with values of the residual error.

In this section, we showed the example of the PEMEM model application to a single orbit. However, the model can be run for a family of close orbits optimizing the orbital parameters to minimize the surface charging risk. Various scenarios can be envisioned for such optimization. If resources allow, the surface charging simulation can be run for each orbital run for one or a few combinations of the electron and proton percentile spectra, in eclipse and in
sunlight. Then, the orbit with the lowest charging potential can be selected. Alternatively, the orbit with a minimum integrated electron flux (for a certain percentile) can be selected, and then, the surface charging simulation can be run for the percentile spectra for this orbit to make sure that the expected charging potential meets the requirements.

## 7. Discussion and Conclusions

The main limitations of the model are its limited coverage of the radial distances ($L = 3.5–6.5$) and the absence of latitudinal dependence. Technically, the range of the radial distances can be easily extended if a corresponding data set is ingested in the model. However, with the distance increasing, the dipolar configuration of the magnetic field changes to the sheet-like configuration of the plasma sheet. Since the plasma sheet is a dynamic structure its geometry evolves fast in response to the Earth's dipole tilt angle and geomagnetic activity variations, the probability of the spacecraft staying inside/outside the plasma sheet should be modeled separately. Note also that the estimation of McIlwain L becomes progressively less accurate with the distance increasing.

As regards the latitudinal dependence, inside the geosynchronous distance, the latitudinal dimension (e.g., dependence on $B/B_{eq}$) can be added to the model design keeping the same model philosophy (weighted sum of the statistical distributions over the spatial bins). For farther distances, the argumentation mentioned in the previous paragraph is valid and should be addressed.

Except for limitations related to the limited spatial coverage of the parent data set, a few other limitations of the percentile model should be mentioned.

* The percentile spectrum flux values at different energy channels are modeled independently. For this reason, they typically are not physically consistent (it is unlikely that the fluxes at 1 and 100 keV are correlated). It is especially critical for protons because the 50 eV–50 keV energy range covers the interval of proton temperature variation. Figure 9f shows examples of such unrealistic spectrum shapes.

* The model low energy limit for protons is ~50 eV. We do not have reliable measurements of colder ions because the measurements in this energy range are usually distorted by the nonzero spacecraft potential. At the same time, the presence of dense cold ions mitigates the surface charging and is important to the surface charging risk assessment.

* The electron and proton fluxes are modeled independently. At the same time, it is highly likely that these quantities are coupled.

* The PEMEM percentile model design does not include the dependence on the geomagnetic activity. Its output corresponds to the average geomagnetic activity during the years 2013–2018 (Van Allen Probes mission lifespan). Although individual geomagnetic activity intensifications cannot be predicted for future missions, the solar cycle phase can be predicted reasonably well and could be included in the model parameterization.

* Although the model can estimate the total time for which a certain flux value will be exceeded, it cannot answer what is the probability that a certain flux value will be exceeded during the mission for a given confidence level. The analysis of the electron flux autocorrelation (see Section 3) has shown that the flux enhancements tend to be organized as events likely corresponding to the substorms. The following very rough estimation can be used: if the time of exceedance is higher than the typical substorm timescale (~1 hr), it is likely that the threshold value will be exceeded. Extending this reasoning further, the probability can be estimated as exceedance time divided by the substorm timescale. Obviously, confidence levels cannot be taken into account in such a simple estimate.

The last two issues are addressed in the advanced version of the PEMEM model, AE index-driven Monte Carlo (AEdMC) model, which will be published in a separate paper.

Among the PEMEM *percentile* model's advantages, we would like to emphasize its simplicity and robustness and its computational speed. Multiple orbital scenarios can be processed very fast during one run and the output can be used to select the best or worst orbital scenario. After that, more advanced models can be used to model that scenario.

Apart from surface charging risk assessment, the model can be used to statistically validate other particle measurement data sets in the inner magnetosphere. If there is a relatively large data set of the electron or proton measurements whose accuracy is compromised, the user can apply the PEMEM model to the corresponding orbit and compare the statistical distributions provided by the model with that obtained empirically from the data set.
<table><caption>**Table 1**
Summary of the Model's Inputs and Outputs</caption><thead><tr><th>Species</th><th>Input</th><th>Settings</th><th>Output</th></tr></thead><tbody><tr><td>Electrons</td><td>Orbital L<sup>a</sup>, MLT<sup>a</sup>, GSE (for the full modeled period)</td><td><ul><li>L-shell value of the model's outer limit</li><li>Specifications of the required percentiles</li><li>Energy of the lower integration limit (E<sub>low</sub>) for the electron integrated flux</li><li>Integrated electron flux threshold value</li></ul></td><td><ul><li>Percentiles of differential electron flux for 19 energies in 1–100 keV range (percentile spectra)</li><li>Parameters of the multi-Maxwellian approximation for percentile spectra</li><li>Percentiles of the electron flux integrated over the E<sub>low</sub>–100 keV interval</li><li>Total time the threshold value is exceeded</li></ul></td></tr><tr><td>Protons</td><td>Orbital L<sup>a</sup>, MLT<sup>a</sup>, GSE (for the full modeled period)</td><td><ul><li>L-shell value of the model's outer limit</li><li>Specifications of the required percentiles</li></ul></td><td><ul><li>Percentiles of differential proton flux for 28 energies in 50 eV–51 keV range (percentile spectra)</li><li>Parameters of the multi-Maxwellian approximation for percentile spectra</li></ul></td></tr></tbody></table>

## 8. Summary

The PEMEM percentile model is designed for surface charging risk assessment for missions with near-equatorial orbits. The model outputs statistical characteristics of the electron and proton plasma environments that are expected to be encountered on a given orbit. These characteristics can be further used as input to the spacecraft surface charging simulations (e.g., https://www.spis.org/) to get a value of the corresponding charging potential. The model's input and output are summarized in Table 1. Note that the same output parameters can be output for the periods when the spacecraft is eclipsed by the Earth.

## Data Availability Statement

All RBSP-ECT data are publicly available at the Web site https://rbsp-ect.newmexicoconsortium.org/rbsp_ect.php.

## Acknowledgments

This work was supported by ESA Contract 4000128226/19/NL/AS. The work of S.D. and N.G. was partly supported by the Project 339329 of the Academy of Finland. The work of NG at the University of Michigan was partly funded by National Aeronautics and Space Administration Grant 80NSSC20K0353, and National Science Foundation Grant 1663770. Processing and analysis of the (HOPE, MagEIS, REPT, or ECT) data by the Van Allen Probes team was supported by Energetic Particle, Composition, and Thermal Plasma (RBSP-ECT) investigation funded under NASA's Prime Contract NAS5-01072.

## References

Blake, J. B., Carranza, P. A., Claudepierre, S. G., Clemmons, J. H., Crain, W. R., Dotan, Y., et al. (2013). The magnetic electron ion spectrometer (MagEIS) instruments aboard the radiation belt storm probes (RBSP) spacecraft. *Space Science Reviews*, 179(1-4), 383-421. https://doi.org/10.1007/s11214-013-9991-8

Bortnik, J., Li, W., Thorne, R. M., & Angelopoulos, V. (2016). A unified approach to inner magnetospheric state prediction. *Journal of Geophysical Research: Space Physics*, 121(3), 2423-2430. https://doi.org/10.1002/2015JA021733

Boscher, D., Bourdarie, S., Friedel, R., & Belian, R. (2003). Model for the geostationary electron environment: POLE. *IEEE Transactions on Nuclear Science*, 50(6), 2278-2283. https://doi.org/10.1109/TNS.2003.821609

Boyd, A. J., Reeves, G. D., Spence, H. E., Funsten, H. O., Larsen, B. A., Skoug, R. M., et al. (2019). RBSP-ECT combined spin-averaged electron flux data product. *Journal of Geophysical Research: Space Physics*, 124(11), 9124-9136. https://doi.org/10.1029/2019JA026733

Boyd, A. J., Spence, H. E., Reeves, G. D., Funsten, H. O., Skoug, R. M., Larsen, B. A., et al. (2021). RBSP-ECT combined pitch angle resolved electron flux data product. *Journal of Geophysical Research: Space Physics*, 126(3), e2020JA028637. https://doi.org/10.1029/2020JA028637

Boynton, R. J., Amariutei, O. A., Shprits, Y. Y., & Balikhin, M. A. (2019). The system science development of local time-dependent 40-keV electron flux models for geostationary orbit. *Space Weather*, 17(6), 894-906. https://doi.org/10.1029/2018SW002128

Boynton, R. J., Balikhin, M. A., Sibeck, D. G., Walker, S. N., Billings, S. A., & Ganushkina, N. (2016). Electron flux models for different energies at geostationary orbit. *Space Weather*, 14(10), 846-860. https://doi.org/10.1002/2016SW001506

Chu, X., Bortnik, J., Li, W., Ma, Q., Denton, R., Yue, C., et al. (2017). A neural network model of three-dimensional dynamic electron density in the inner magnetosphere. *Journal of Geophysical Research: Space Physics*, 122(9), 9183-9197. https://doi.org/10.1002/2017JA024464

Claudepierre, S. G., O'Brien, T. P., Looper, M. D., Blake, J. B., Fennell, J. F., Roeder, J. L., et al. (2019). A revised look at relativistic electrons in the Earth's inner radiation zone and slot region. *Journal of Geophysical Research: Space Physics*, 124(2), 934-951. https://doi.org/10.1029/2018JA026349

Davis, V. A., Mandell, M. J., & Thomsen, M. F. (2008). Representation of the measured geosynchronous plasma environment in spacecraft charging calculations. *Journal of Geophysical Research*, 113(A10), A10204. https://doi.org/10.1029/2008JA013116

DeForest, S. E. (1973). Electrostatic potentials developed by ATS-5. In R. J. L. Grard (Ed.), *Photon and particle interactions with surfaces in space* (pp. 263-267). Springer Netherlands.

Denton, M. H., Thomsen, M. F., Jordanova, V. K., Henderson, M. G., Borovsky, J. E., Denton, J. S., et al. (2015). An empirical model of electron and ion fluxes derived from observations at geosynchronous orbit. *Space Weather*, 13(4), 233-249. https://doi.org/10.1002/2015SW001168

Fernandes, P. A., Larsen, B. A., Thomsen, M. F., Skoug, R. M., Reeves, G. D., Denton, M. H., et al. (2017). The plasma environment inside geostationary orbit: A Van Allen probes HOPE survey. *Journal of Geophysical Research: Space Physics*, 122(9), 9207-9227. https://doi.org/10.1002/2017JA024160

<footnote><sup>a</sup>Coordinates should be calculated using Olson and Pfitzer (1977) (OP77Quiet) model.</footnote>
Feynman, J., Spitale, G., Wang, J., & Gabriel, S. (1993). Interplanetary proton fluence model: Jpl 1991. *Journal of Geophysical Research*, 98(A8), 13281-13294. https://doi.org/10.1029/92JA02670

Funsten, H. O., Skoug, R. M., Guthrie, A. A., MacDonald, E. A., Baldonado, J. R., Harper, R. W., et al. (2013). Helium, oxygen, proton, and electron (HOPE) mass spectrometer for the radiation belt storm probes mission. *Space Science Reviews*, 179(1-4), 423-484. https://doi.org/10.1007/s11214-013-9968-7

Ganushkina, N. Y., Swiger, B., Dubyagin, S., Matéo-Vélez, J.-C., Liemohn, M. W., Sicard, A., & Payan, D. (2021). Worst-case severe environments for surface charging observed at LANL satellites as dependent on solar wind and geomagnetic conditions. *Space Weather*, 19(9), e2021SW002732. https://doi.org/10.1029/2021SW002732

Garrett, H. B. (1981). The charging of spacecraft surfaces. *Reviews of Geophysics*, 19(4), 577-616. https://doi.org/10.1029/RG019i004p00577

Garrett, H. B., & Whittlesey, A. (2012). Introduction to the physics of charging and discharging. In *Guide to mitigating spacecraft charging effects* (pp. 6-25). John Wiley & Sons, Ltd. https://doi.org/10.1002/9781118241400.ch2

Ghojogh, B., Ghojogh, A., Crowley, M., & Karray, F. (2020). Fitting a mixture distribution to data: Tutorial. https://doi.org/10.48550/arXiv:1901.06708

Ginet, G. P., O'Brien, T. P., Huston, S. L., Johnston, W. R., Guild, T. B., Friedel, R., et al. (2013). AE9, AP9 and SPM: New models for specifying the trapped energetic particle and space plasma environment. *Space Science Reviews*, 179(1-4), 579-615. https://doi.org/10.1007/s11214-013-9964-y

Hardy, D. A., Holeman, E. G., Burke, W. J., Gentile, L. C., & Bounar, K. H. (2008). Probability distributions of electron precipitation at high magnetic latitudes. *Journal of Geophysical Research*, 113(A6), A06305. https://doi.org/10.1029/2007JA012746

Hastie, T., Tibshirani, R., & Friedman, J. (2009). *The elements of statistical learning: Data mining, inference and prediction* (2nd ed.). Springer.

Koons, H., Mazur, J., Selesnick, R., Blake, J., Fennell, J., Roeder, J., & Anderson, P. (1998). The impact of the space environment on space systems. In *Paper presented at 6th spacecraft charging technology* (pp. 7-11). (AFRL-VS-TR-20001578).

Leach, R. D. (1995). Failures and anomalies attributed to spacecraft charging (Tech. Rep.). United States, Alabama, 35812. Retrieved from https://ntrs.nasa.gov/api/citations/19960001539/downloads/19960001539.pdf

Mandell, M., Davis, V., Cooke, D., Wheelock, A., & Roth, C. (2006). Nascap-2k spacecraft charging code overview. *IEEE Transactions on Plasma Science*, 34(5), 2084-2093. https://doi.org/10.1109/TPS.2006.881934

Matéo-Vélez, J.-C., Paulmier, T., Sicard, A., Dirassen, B., & Payan, D. (2019). Experimental investigation of surface potentials of materials under electron spectra representative of GEO and MEO worst case environments. *IEEE Transactions on Plasma Science*, 47(8), 3885-3890. https://doi.org/10.1109/TPS.2019.2925413

Matéo-Vélez, J.-C., Pignal, C., Balcon, N., Payan, D., Sarrailh, P., & Hess, S. (2014). GEO spacecraft worst-case charging estimation by numerical simulation. In *Spacecraft charging technology conference 2014 (13th SCTC)*. Retrieved from https://onera.hal.science/hal-01070320

Matéo-Vélez, J.-C., Sicard, A., Payan, D., Ganushkina, N., Meredith, N. P., & Sillanpää, I. (2018). Spacecraft surface charging induced by severe environments at geosynchronous orbit. *Space Weather*, 16(1), 89-106. https://doi.org/10.1002/2017SW001689

Matéo-Vélez, J.-C., Sicard-Piet, A., Lazaro, D., Inguimbert, V., Sarrailh, P., Hess, S., et al. (2016). Severe geostationary environments: Numerical estimation of spacecraft surface charging from flight data. *Journal of Spacecraft and Rockets*, 53(2), 304-316. https://doi.org/10.2514/1.A33376

Nakamura, M., Nakamura, S., Kawachi, R., & Toyoda, K. (2018). Assessment of worst GEO plasma environmental models for spacecraft surface charging by SPIS. *Transactions of the Japan Society for Aeronautical and Space Sciences, Aerospace Technology Japan*, 16(6), 556-560. https://doi.org/10.2322/tastj.16.556

Olson, W. P., & Pfitzer, K. A. (1977). *Magnetospheric magnetic field modeling* (annual scientific report No. ADA037492). 5301 bolsa avenue. MCDONNELL DOUGLAS ASTRONAUTICS COMPANY.

Sarno-Smith, L. K., Larsen, B. A., Skoug, R. M., Liemohn, M. W., Breneman, A., Wygant, J. R., & Thomsen, M. F. (2016). Spacecraft surface charging within geosynchronous orbit observed by the Van Allen Probes. *Space Weather*, 14(2), 151-164. https://doi.org/10.1002/2015SW001345

Sarrailh, P., Matéo-Vélez, J.-C., Hess, S. L. G., Roussel, J.-F., Thiébault, B., Forest, J., et al. (2015). SPIS 5: New modeling capabilities and methods for scientific missions. *IEEE Transactions on Plasma Science*, 43(9), 2789-2798. https://doi.org/10.1109/TPS.2015.2445384

Sawyer, D. M., & Vette, J. I. (1976). Ap-8 trapped proton environment for solar maximum and solar minimum (Tech. Rep.). United States. Retrieved from http://inis.iaea.org/search/search.aspx?orig_q=RN:09351076

Sicard, A., Boscher, D., Bourdarie, S., Lazaro, D., Standarovski, D., & Ecoffet, R. (2018). GREEN: The new global radiation earth ENvironment model (beta version). *Annales Geophysicae*, 36(4), 953-967. https://doi.org/10.5194/angeo-36-953-2018

Sicard, A., Boscher, D., Lazaro, D., Bourdarie, S., Standarovski, D., & Ecoffet, R. (2019). New model for the plasma electrons fluxes (part of GREEN model). *IEEE Transactions on Nuclear Science*, 66(7), 1738-1745. https://doi.org/10.1109/TNS.2019.2923005

Sicard-Piet, A., Bourdarie, S., Boscher, D., & Friedel, R. (2006). A model for the geostationary electron environment: POLE, from 30 keV to 5.2 MeV. *IEEE Transactions on Nuclear Science*, 53(4), 1844-1850. https://doi.org/10.1109/TNS.2006.877878

Sicard-Piet, A., Bourdarie, S., Boscher, D., Friedel, R. H. W., Thomsen, M., Goka, T., et al. (2008). A new international geostationary electron model: IGE-2006, from 1 keV to 5.2 MeV. *Space Weather*, 6(7), S07003. https://doi.org/10.1029/2007SW000368

Sillanpää, I., Ganushkina, N. Y., Dubyagin, S., & Rodriguez, J. V. (2017). Electron fluxes at geostationary orbit from goes maged data. *Space Weather*, 15(12), 1602-1614. https://doi.org/10.1002/2017SW001698

Simms, L. E., Ganushkina, N. Y., van de Kamp, M., Liemohn, M. W., & Dubyagin, S. (2022). Using ARMAX models to determine the drivers of 40-150 keV GOES electron fluxes. *Journal of Geophysical Research: Space Physics*, 127(9), e2022JA030538. https://doi.org/10.1029/2022JA030538

Simms, L. E., Ganushkina, N. Y., Van der Kamp, M., Balikhin, M., & Liemohn, M. W. (2023). Predicting geostationary 40-150 keV electron flux using ARMAX (an autoregressive moving average transfer function), RNN (a recurrent neural network), and logistic regression: A comparison of models. *Space Weather*, 21(5), e2022SW003263. https://doi.org/10.1029/2022SW003263

Stepanov, N. A., Sergeev, V. A., Sormakov, D. A., Andreeva, V. A., Dubyagin, S. V., Ganushkina, N., et al. (2021). Superthermal proton and electron fluxes in the plasma sheet transition region and their dependence on solar wind parameters. *Journal of Geophysical Research: Space Physics*, 126(4), e2020JA028580. https://doi.org/10.1029/2020JA028580

Swiger, B. M., Liemohn, M. W., Ganushkina, N. Y., & Dubyagin, S. V. (2022). Energetic electron flux predictions in the near-earth plasma sheet from solar wind driving. *Space Weather*, 20(11), e2022SW003150. https://doi.org/10.1029/2022SW003150

Thomsen, M. F., Henderson, M. G., & Jordanova, V. K. (2013). Statistical properties of the surface-charging environment at geosynchronous orbit. *Space Weather*, 11(5), 237-244. https://doi.org/10.1002/swe.20049

Vette, J. I. (1991). The NASA/National space Science data Center trapped radiation environment model program, 1964 - 1991 (Tech. Rep.). United States. Retrieved from http://inis.iaea.org/search/search.aspx?orig_q=RN:24044692
Whipple, E. C. (1981). Potentials of surfaces in space. *Reports on Progress in Physics*, *44*(11), 1197-1250. https://doi.org/10.1088/0034-4885/44/11/002

Zhelavskaya, I. S., Shprits, Y. Y., & Spasojević, M. (2017). Empirical modeling of the plasmasphere dynamics using neural networks. *Journal of Geophysical Research: Space Physics*, *122*(11), 11227-11244. https://doi.org/10.1002/2017JA024406

Zhou, Q., Chen, Y., Xiao, F., Zhang, S., Liu, S., Yang, C., et al. (2022). A machine-learning-based electron density (MLED) model in the inner magnetosphere. *Earth and Planetary Physics*, *6*(4), 350-358. https://doi.org/10.26464/epp2022036