Simulating the Earth's Outer Radiation Belt Electron Fluxes and Their Upper Limit: A Unified Physics‐Based Model Driven by the AL Index

Using particle and wave measurements from the Van Allen Probes, a 2‐D Fokker‐Planck simulation model driven by the time‐integrated auroral index (AL) value is developed. Simulations for a large sample of 186 storm‐time events are conducted, demonstrating that the AL‐driven model can reproduce flux enhancement of the MeV electrons. More importantly, the relativistic electron flux enhancement is determined by the sustained strong substorm activity. Enhanced substorm activity results in increased chorus wave intensity and reduced background electron density, which creates the required condition for local electron acceleration by chorus waves to MeV energies. The appearance of higher energy electrons in radiation belts requires a higher level of cumulative AL activity after the storm commencement, which acts as a type of switch, turning on progressively higher energies for longer and more intense substorms, at critical thresholds.


Introduction
The behavior of high-energy electrons, typically ranging from 100s keV to 10 MeV, in the Earth's dynamic radiation belts has been studied extensively, particularly following the comprehensive measurements of these electron fluxes from the Van Allen Probes (Baker et al., 2012;Mauk et al., 2014).The key question to answer is how these particles are accelerated to such high energies as they are known to be particularly hazardous to spacecraft.Among various physical mechanisms, the most crucial may be chorus waves, which can provide efficient energy diffusion to accelerate electrons from 100s keV to MeVs (e.g., Horne & Thorne, 1998;Hua, Bortnik, & Ma, 2022;Li & Hudson, 2019;Li et al., 2014;Summers et al., 1998;Tu et al., 2014).Intense whistlermode chorus waves are excited in the magnetosphere during geomagnetically active periods (e.g., Li et al., 2009), and the flux level of energetic electrons varies significantly in response to changes in the properties of the solar wind or geomagnetic conditions.Thus, linking the solar wind or the geomagnetic perturbations caused by the solar wind with the dynamics of electrons has become an essential method for understanding and predicting the radiation belt dynamics.
Statistical analyses examined the correlation between high energy electrons and solar wind conditions and geomagnetic activities, for instance, Meredith et al. (2003) found that the most significant electron flux enhancements outside the plasmapause were associated with periods of prolonged substorm activity.Zhao et al. (2017) demonstrated there is a good relation between ultra-relativistic electron enhancement (>3 MeV) and solar wind speed, as well as the AL index.Zhao, Baker, Li, Jaynes, and Kanekal (2019) further showed that even small and moderate storms can accelerate electrons to multi-MeV energies.However, statistical studies alone have difficulty explaining the exact causative physical mechanisms involved, as different driving factors typically occur simultaneously (e.g., Tsurutani et al., 2006).Event studies have also been performed to investigate the underlying physical mechanisms causing the relativistic electron acceleration.The mechanisms of particle acceleration are generally divided into two categories or a combination of the two effects: local acceleration by chorus waves (e.g., Reeves et al., 2013;Thorne et al., 2013) and inward radial diffusion from higher L-shells (e.g., Ozeke et al., 2019Ozeke et al., , 2020)).Drozdov et al. (2015) used the Kp index to parameterize the diffusion coefficients to drive a long-term Fokker-Planck simulation, but the results show a significant discrepancy for the ultrarelativistic electrons.Our previous research using machine learning (ML) methods successfully simulated electrons at relativistic energy levels and achieved better results than general physical models (Chu et al., 2021;D. Ma et al., 2022), but it was less effective in simulating electrons in the ultra-relativistic range, likely due to a much smaller set of events at these high energies.Thus, the present study focuses on understanding electron flux behavior, informed by the results of our previous machine-learning-based models, especially for the ultrarelativistic range.Agapitov et al. (2019) and Allison et al. (2021) suggested that during geomagnetically active times, the electron density drops to very low values and creates favorable conditions for ultra-relativistic electron acceleration by chorus waves, consistent with earlier simulations by Thorne et al. (2013).Hua, Bortnik, Chu, et al. (2022), andHua et al. (2023) show the maximum electron flux at different energy levels is strongly correlated to the cumulative effects of substorms instead of storms.Our previous ML models were driven only by geomagnetic indices and solar wind parameters, and with interpretable ML methods, we automatically discovered the integrated AL is a key parameter in the acceleration mechanism for relativistic and sub-relativistic electrons (D.Ma et al., 2023Ma et al., , 2024)).These studies indicate that we can use the AL index to model the time-varying physical parameters needed for traditional Fokker-Planck simulations, thereby obtaining a general model driven by the auroral index.We use this idea as a basis for further investigation into the acceleration mechanisms of ultrarelativistic electrons, and perform a series of simulations that cover all the 2012-2018 storm time events.We describe our methods in Section 2 and demonstrate the result in Section 3. We discuss our conclusions in Section 4.

Events Selection and Simulation Methodology
We attempt to simulate all the storm-time enhancement events occurring during the Van Allen Probe era and investigate the conditions that led to the acceleration of relativistic and ultra-relativistic electrons.The criteria used to select simulation events are that the minimum SYM-H < 25 nT, and there is at least a 2-fold increase at 597 keV electron flux during the storm event.We selected 186 events from September 2012 to September 2018 based on the above criteria, and then we performed simulations of each event.
Figures 1a-1o presents one of our selected storms occurring on 2012-10-09 which has been well studied in the past (Thorne et al., 2013).The time t 0 is defined as SYM-H minimum (same for all the events) on 10-09, the event start date is defined as t 0 1 day as the left vertical dashed line, and the simulation period is from t 0 to t 0 + 2 days (taken to be the same for all the events) as shown in the shaded area.The most significant challenge for us is how to perform 186 simulations since the traditional Fokker-Planck simulation approach requires a large amount of computational sources to calculate the diffusion coefficients (Ni et al., 2008).We use a diffusion matrix precalculated from statistical results of wave properties with different density values and a wave amplitude B w = 100pT which are then scaled by the modeled wave amplitude and interpolated as a function of density, as done in the lookup table method to calculate diffusion coefficients (Hua, Bortnik, Kellerman, et al., 2022).The statistical wave properties used in this study are based on the Van Allen Probes survey (Li et al., 2016) and listed in Table 1.
We use the integral hourly AL index value obtained from the OMNI data set summed over the previous 3 hr as the driver for each hourly step to model the wave amplitude and density profile.Figure 2a shows the correlation between the cumulative AL and the power of the magnetic field fluctuations.The data is obtained from Van Allen Probes observations near the heart of the outer belt acceleration region, at L = 5, collected over L from 4.8 to 5.2 and observed within 15°from the magnetic equator when the probes are outside of plasmapause.The statistical results are divided into 6 levels of activity and four different MLT sectors which we will use to calculate the average diffusion coefficients.The acceleration over the dayside MLT sector 12-20 hr is ignored due to the low chorus wave activity that occurs near the equator in that region.It is clear that as geomagnetic activity increases, the associated wave magnetic intensity increases, especially in the night-dawn sectors, consistent with previous studies (e.g., Teng et al., 2018).It is worth mentioning that the vertical lines extending from each data point represent only 10% of the standard deviation in the measurements for the sake of clarity, which means that there are huge uncertainties for large wave amplitude.Figure 2b shows the electron density measurement outside the plasmapause near L = 5, similar to Figure 2a.Importantly, a global decrease in density is observed during high geomagnetic activity.Such low-density observations (below 10 cm 3 ) can contribute to the ultra-relativistic electron acceleration as the wave phase velocity and electron resonance energy increase with a density reduction (Agapitov et al., 2019;Allison et al., 2021).The combination of wave power enhancement and density decrease during the intense substorm period can result in a factor of 10 enhancement in the average pitch angle ) with cutoffs at θ min and θ max .and momentum diffusion coefficients, as shown in Figures 2c and 2d.The detail of the database of the wave and density can be found in Q. Ma et al. (2023).
The chorus waves typically cause both pitch angle scattering loss and acceleration of electrons in the radiation belts.The scattering and acceleration rates are represented by the bounce-averaged pitch angle 〈D α eq α eq 〉 , momentum 〈D pp 〉 , and mixed pitch angle-momentum 〈D α eq p 〉 diffusion coefficients.To evaluate the electron flux variations due to chorus, we perform a local acceleration simulation driven by chorus waves by numerically solving the 2D Fokker-Planck equation (Q.Ma et al., 2012), expressed as: (S(α eq ) sinα eq cos α eq 〈D α eq α eq 〉 ∂f ∂α eq ) + 1 S(α eq ) sinα eq cosα eq ∂ ∂α eq (S(α eq ) sinα eq cos α eq p〈D α eq p 〉 ∂f ∂p ) where f is the phase space density (PSD), t is the time, and S(α eq ) represents the bounce period approximated by S(α eq ) = 1.38 0.32 sin(α eq ) 0.32 ̅̅̅̅̅̅̅ ̅ sin √ (α eq ) (Lenchek et al., 1961).The lower energy boundary is set as the electron PSD at 235 keV provided by our ML model (D.Ma et al., 2022), which is primarily driven by the AL index during the post-storm period (D.Ma et al., 2023Ma et al., , 2024)), as shown in the solid line in Figure 1e.The lower boundary is set to be higher than in previous simulations as we focus on the relativistic electrons.The initial PSD distribution is adopted from a spline fitting from 235 keV to 7.7 MeV for each individual event, and the upper energy boundary is set at 10 MeV with a constant value 10 13 cMeV 1 cm 1 ) 3 .The initial pitch angle distribution is assumed to be f (α eq ,p) = f (α eq = 90°,p) sin α eq .With the time-varying wave amplitude and total electron density from statistical results, we now update the diffusion coefficients every hour with the changing integrated AL driver.Figures 1p and 1q demonstrate the observation and simulation results.The model captures the rapid enhancement at high-energy levels clearly and the simulation result is consistent with a few previous studies on the same event (Hua, Bortnik, & Ma, 2022;Thorne et al., 2013).

Results
Our focus in this study is not to reproduce a single event but to capture the collective behavior of the model on radiation belt acceleration of relativistic and ultra-relativistic electron fluxes occurring consistently over many events.Figure 3 shows the results of all 186 simulations at different energy channels and comparison with observation at the end of simulation time.Notably, the model underestimates the electron PSD below 2 MeV shown in Figures 3a-3c.Such a modeling bias is expected, as our simulation only considers local acceleration and does not take radial diffusion into account.Specifically, the radial diffusion process will transport energetic electrons from higher L shells inward, resulting in the additional enhancement at our simulation position L = 5 (Zhao, Baker, Li, Malaspina, et al., 2019).The radial diffusion coefficients may be higher at lower energies (Liu et al., 2016), which leads to a more significant bias in the lower energy range.In addition, the bias may also come from the fast injection process as the injected electron energy during storm time can be higher than the lower boundary (235 keV) in the simulation.Another possible issue is the difference between the fluxes measured by MagEIS and the REPT instruments at 1-2 MeV energies.For example, the REPT flux measurement at 1.8 MeV is about 2-10 times higher than the MagEIS at a similar energy channel 1.728 MeV (Boyd et al., 2019).Figures 3d-3g show a very good agreement between the simulation and observation at the relativistic energy range.Importantly, multiple years of observations show that the ultrarelativistic electron flux is only sometimes increased after storm times, especially at higher energy levels (Baker et al., 2019).Figure 3i shows the model performance at 6.3 MeV with red lines separating the results into different regions.It is clearly seen that the model captures not only the enhancement but also most non-enhancement events where the post-storm flux level remains at the noise level.These results demonstrate that our method parameterized by the AL index can successfully capture the threshold between enhancement and non-enhancement for electron fluxes.
We now aim to understand what leads to the critical enhancement threshold levels that control high energy flux levels.Given that our model is driven by the AL index, we investigate the relation between the result electron flux and integral AL (integrated with 5-min resolution AL from t 0 1 day to t 0 + 2 days) and compare the results from observations and simulations.Figure 4a shows the flux observations at the end of each event.The flux is normalized by the minimum and maximum flux for each energy channel among all the events.Each dot represents one flux observation color-coded by the normalized flux.At lower energies (<1 MeV), the fluxes have a strong linear relationship with the integral of the AL index, which is seen more clearly in Figure 4b.At higher energies (>1 MeV), the flux remains at a very low level when the integral AL is low, but becomes enhanced only when the integral AL surpasses a threshold x 0 .We thus define a criterion to find such a critical threshold x 0 : where the x is the integral AL of each event, and y(x) is the corresponding flux value.The "switch" function I is 1 when y(x) > 0.05, otherwise is 0. So, such criterion can find the critical threshold x 0 where 70% of the events can end with a clear enhancement when the integral AL is larger than x 0 .The observation threshold is shown as dashed lines in Figures 4b-4d and as a dash-cross line for each energy channel in Figure 4a, and the threshold increases with increasing energy as expected.Figures 4e-4g shows the model results.We calculate the threshold with the same criterion for the modeled flux in Figures 4e-4g, and the results compare very well with the observation, which confirms our notion that electron flux enhancements only become significant when substorms of a certain intensity persist for a sufficient duration, acting as a type of "turn on switch."Moreover, our physics-based model demonstrates that this is due to the controlling effect of prolonged substorm activity on increasing the chorus wave intensity and lowering the electron density in the radiation belt acceleration region (Figure 2).

Conclusions and Discussion
In this article, we investigated the local acceleration process of radiation belt electron fluxes driven by chorus waves and electron density provided by our models under various geomagnetic conditions.Using a lookup table method, we performed a series of event-specific radiation belt diffusion simulations to investigate electron acceleration during storm time during the Van Allen Probe era using an ensemble of 186 storm events.Our findings reveal that: 1.The chorus wave intensity and plasma density are both controlled by accumulated substorm activity: the stronger the substorm, the greater the wave intensity and the lower the electron density.This phenomenon is most pronounced in the night to the dawn MLT region which is critical for electron acceleration due to the intense, equatorial chorus waves that occur there.(Boyd et al., 2019) and 60 events that are above it, thus identifying a nature boundary between enhancement and non-enhancement events.

Geophysical Research Letters
10.1029/2024GL109169 2. We develop a general substorm-driven, local acceleration model using a lookup table method and simulate a large number of strom-time acceleration events at a range of energies, with good performance, illustrating that the accumulated substorm activities is the key parameter controlling relativistic and ultra-relativistic electron acceleration.3. The electron fluxes following geomagnetic disturbances show a clear threshold integral AL, which controls whether a particular electron flux energy channel will be enhanced or not.These results indicate that the continuously elevated substorm activity is the determining factor for ultra-relativistic electron acceleration, and we demonstrate that this is due to the influence of substorms on both the electron density and wave intensity.
This study sheds new light on modeling and understanding the acceleration mechanism of ultra-relativistic electron fluxes.As far as we are aware, this is the first time that the critical geomagnetic activity threshold levels, physical models, and the underlying conditions under which electron flux can increase at high energies have been directly linked.The clear relationship between relativistic electron flux enhancement and sustained substorm activity sheds new light on the drivers of acceleration leading to potentially more effective and accurate model development in the future.For the sake of simplicity, we only consider local acceleration by lower band chorus waves in the heart of the outer radiation belt acceleration region at L = 5, which is close to the peak of the chorus wave intensity (Aryan et al., 2021) and the maximum electron fluxes in the radial profile in the outer belt (Hua, Bortnik, Chu, et al., 2022).However, the model's good performance suggests that this processes may be the predominant mechanism and we can gradually add other mechanisms based on this model in the future studies.
Although in our model, substorms only affect the wave amplitude and the total electron density, it is important to note that enhanced injections are implicitly included in the boundary and initial conditions, and their relationship with substorms has been generally known and appreciated (e.g., Reeves et al., 1990).When selecting events, we did not impose strict constraints on different Dst index, but we cannot conclude that Dst has no significant impact on relativistic electron enhancement.We followed the results obtained from ML to start with the AL index and studied the electron acceleration outside the plasmapause.The relationship between relativistic electron acceleration and geomagnetic storm strength still merits careful investigation.

Figure 1 .
Figure 1.An example of a selected case and its associated FP simulation which occurred on 2012-10-09.(a-d) Geomagnetic indices such as AL and, SYM-H, and solar wind dynamic pressure P SW , and solar-wind speed V SW , (e) black solid line shows the machine learning result of 235 keV electron flux at L = 5 and MLT = 6 (black line) together with observations (blue dots).(f-o) The observation phase space density (PSD) of different energy channels, where black horizontal dashed lines represent L = 5. (p) Observed PSD color-coded with time.(q) FP simulation result of the event.The vertical dashed lines in (p) and (q) represent 0.6 MeV (red), 3.4 MeV (yellow), and 6.3 MeV (blue) energies.

Figure 2 .
Figure 2. Statsical results of AL driven wave and AL driven density models at L = 5 and the corresponding average diffusion coefficients.(a) The statistical result of wave magnetic intensity at four different MLT sectors.(b) The statistical result of electron density.The dots represent the mean value and the error bar represent std/10.(c) The average pitch angle diffusion coefficients 〈D αα 〉 and (d) the average momentum diffusion coefficients as a function of pitch angle at different energies.

Figure 3 .
Figure3.Model results compared with the observations at different energy channels.The model performance for the occurrence of 6.3 MeV flux enhancement is marked by the red lines and event numbers.The model captures 94 events that are below the noise level at 10 11(Boyd et al., 2019) and 60 events that are above it, thus identifying a nature boundary between enhancement and non-enhancement events.

Figure 4 .
Figure 4.The modeled electron flux threshold comparison with observation.(a) The observed flux, normalized and colorcoded by each channel's maximum and minimum flux.The dashed line and cross show the observed threshold (black) and modeled threshold (red).(b)-(d) The observation results and the threshold at different energies.(e)-(g) The simulation results and the threshold at different energies.

Table 1
The Wave Properties Used in the Diffusion Coefficient Matrix Calculation w 30°20-24 MLT 0°-20°N ote.The wave frequency is determined by a Gaussian distribution B 2