Two methods to analyze radial diffusion ensembles: The perils of space- and time-dependent diffusion

Radial diffusion is a phenomenon studied in both space and fusion plasmas. In this work, we investigate how the radial dependence of that diffusion interacts with initial and boundary conditions. We provide here a guide for different readers to navigate this paper. First, in the introduction, we provide a broad introduction of the motivations behind understanding the uncertainty in radial diffusion modeling in near-Earth space. Plasma physicists without a radiation belt background may wish to review the application of radial diffusion in near-Earth space in Sec. II. For all readers, we explicitly list our goals in Sec. III, introducing labels which are used throughout the manuscript as we develop each research question, find relevant results, and then discuss our conclusions. Section IV contains details of the numerical models and develops the properties we require in an analysis metric, thereby motivating our use of time to monotonicity and the mass- and energy-like quantities $N,E$⁠. We find and present the most significant ways in which these quantities vary in Sec. V. Radiation belt modelers and space weather physicists may be particularly interested in the suggestions we make for future based on our findings in Sec. VI, where we also compare our results to current modeling practices. Some of our more significant conclusions relate to the importance of the outer boundary and the particle gradients instead, which are also discussed in Sec. VI. Where possible, each section is self-contained, to enable those with different interest to find the relevant sections useful.

Earth's radiation belts are a region of highly energized particles, magnetically trapped by the Earth's magnetosphere. The trapped particles undergo several types of periodic motion, the slowest of these being the drift around the Earth. Electromagnetic perturbations on timescales of the drift of electrons around the Earth will scatter those electrons onto orbits closer to, or more distant from, the Earth; this is radial diffusion. Drift timescales can vary with particle energy from minutes to days but is typically considered to be on the order of a few hours. Radial diffusion is one of the major drivers of Earth's radiation belts and due to the variation in timescales of radiation belt particles motions, an approximation for radiation belt modeling can be made using solely this mechanism to reproduce the broad dynamics, although not shorter timescale phenomena such as dropouts.^1,2 This simple yet effective model enables us to explore ways of including uncertainty in our radiation belt models using ensembles. A full model of the radiation belts would require us to acknowledge that they are part of a complex, interdependent system, with consequently larger amounts of uncertainty. Radial diffusion modeling using the Fokker–Planck formalism is reviewed in more detail in Sec. II; briefly, radiation belt modeling using Fokker–Planck simulations does not use the motion of individual particles but instead averages over motion on larger scales, where wave–particle interactions across many scales are included using diffusion coefficients. Where the wave–particle interactions are well understood, variability should be properly characterized in order to capture the correct uncertainty. Where these interactions are not well understood, large uncertainties can indicate inadequate modeling. We see uncertainty across orders of magnitude when estimating $DLL$ in units of $days−1$⁠. For example, using different $DLL$ models on diffusion estimated by particle tracing in an MHD simulation spans four orders of magnitude at $L=5$ (i.e., spanning minutes to days),³ while empirical models parameterized by geomagnetic activity still vary by more than one order of magnitude even at the same location and geomagnetic activity level.⁴ Unsurprisingly, varying the spatiotemporal distribution of diffusion coefficients changes the final particle distribution.⁵ Previous work to address the uncertainty in diffusion coefficients shows that probabilistic inference of these coefficients outperforms our current empirical models.⁶ Uncertainty can be both inherent to the system and an indicator of a poorly described system; both these situations can be addressed using appropriate modeling techniques.

Ensembles can be used to quantify uncertainty in radiation belt modeling. Ensemble modeling involves running a simulation multiple times, with a variety of input conditions or parameter settings, to represent the unknowns in a given system. The impact of these unknowns on the final state can be quantified; alternatively, variability across model outputs provide a measure of uncertainty on that final state. Probability distribution of model outputs provide us with a better understanding of the uncertainty associated with our models. There are several sources of uncertainty in radiation belt modeling, including uncertainty due to approximations or physical unknowns, uncertainty in observations (i.e., in the measured value and in the spacecraft location), uncertainty in boundary conditions or model settings, and physical uncertainty inherent to the system. Physical uncertainty arises from the fact that deterministic models may exhibit chaotic behavior if they are particularly sensitive to initial conditions. Given that the magnetosphere is a complex system and that observations are spare compared to spatiotemporal scales of interest, it is extremely likely that we will need some way to include this chaotic deterministic character. Furthermore, the computational requirements of our modeling mean that we have sub-grid physics: physics on smaller scales that must be included, but cannot be fully modeled numerically. Uncertainty in driving parameters, such as upstream solar wind and ultra-low-frequency (ULF) waves driving radial diffusion, can also affect the accuracy of the model. Example sources of uncertainty include the properties of the solar wind and how this drives ULF waves. In turn, the uncertainty in magnitude and spatiotemporal occurrence of ULF waves will result in uncertainty in radiation belt models. In addition to this uncertainty chain, modeling of each physical process adds further sources of uncertainty, for example, to include the impact of ULF waves, it is typically assumed that the ULF azimuthal mode number is $m=1$⁠, which does not hold in in situ observations.^7,8 As radiation belt modeling improves, it becomes increasingly important to understand how all of these sources of uncertainty impact our final output to be able to analyze ensembles.

As uncertainty becomes increasingly important in radiation belt modeling, methods to account for uncertainty such as data assimilation and statistical methods are becoming the state-of-the-art method for modeling radiation belts,^9–11 while complex systems approaches are being applied to understand the underlying physics.¹² Ref. 13 states that the field of space physics needs ensembles and methods to deal with and analyze ensembles, in order to manage uncertainty in parameterizations and in various parts of numerical space weather prediction. Ensembles, and other probabilistic or statistical methods, are already becoming the norm.^6,14–16

Satellite operators and national meteorological organizations are increasingly compelled to include space weather forecasting as part of their services, including the radiation environment. Geomagnetically trapped particles in the radiation belts have been established as a hazard to spacecraft due to processes such as surface charging, deep dielectric charging, and single upset events.^17,18 The loss of services, such as the Internet, would cause severe socio-economic issues such as losses to navigation, finances (including card-based payments), and tracking/allocation of emergency services and commercial aircraft. These economic motivations for adapting modeling this fascinating system have also encouraged the space plasma physics community to adopt techniques mastered by meteorologists. However, the domain dependence of the diffusion coefficient means that using and analyzing ensemble members is not simple, rendering ensembles less meaningful than desired.

To understand ensemble models for radial diffusion (and for radiation belt modeling more generally), we need to understand how simple model changes and initial and boundary conditions change the outputs, before we include variability to represent uncertainty in physical conditions.

Radiation belt modeling is typically done using three adiabatic invariants; three quantities associated with the periodic motions of trapped, highly charged particles in Earth's magnetosphere. At relatively low energies (e.g., electrons of 1-hundreds of keV), particles are better described using the ring current paradigm. At higher energies (e.g., relativistic electrons $≥0.5$ MeV), we use the geomagnetically trapped radiation belt description.^19,20 Periodic motions of particles in a slowly changing conservative field (such as electromagnetic fields) correspond to conserved quantities; adiabatic invariants. The first adiabatic invariant, $μM$⁠, corresponds to the magnetic moment as particles gyrate in an electromagnetic field. The second invariant corresponds to momentum along the bounce path, as the gyrating particle travels up and down the field line, reflected by the stronger magnetic field “magnetic bottle”). These motions are on scales of microseconds and seconds, respectively. Under the adiabatic approximation, if the system is changing slowly, then these quantities are conserved. On the other hand, if the system is perturbed on timescales comparable to these motions, the quantities $μM$ or J would not be conserved. On a much longer timescale, the periodic motion of electrons drifting around the Earth corresponds to the third adiabatic invariant $Φ$⁠, the magnetic flux enclosed by this drift path. This can be written as $Φ=2πB0RE3R0$ for particles remaining at the equator, with $R0$ the equatorial distance, $B0$ the magnetic field strength at the Earth's surface, and $RE$ the radius of the Earth. One can, therefore, model the radiation belts using these conserved quantities, by tracking the distribution of particles in this adiabatic invariant space and assuming that the change in this space can be modeled using diffusion—where diffusion between these occurs when electromagnetic perturbations occur on the temporal and spatial scales of the corresponding periodic motions. Diffusion is described using diffusion coefficients $Dαβ$⁠. Typically, cross-terms corresponding to the third adiabatic invariant are considered to be zero; on the whole, changes to the third adiabatic invariant (“radial diffusion”) is considered to be on a separate timescale and can, therefore, be approximated alone (e.g., Ref. 21) The full diffusion model of the radiation belts has been used in a number of places, e.g., Refs. 21–24.

For the third adiabatic invariant, one can equivalently consider the drift shells instead of the magnetic flux conserved by the drift orbit. Since the drift shell is uniquely defined by its intersection with the equatorial magnetic field, it is often easier and more intuitive to parameterize the third adiabatic invariant using some form of $L*$ parameter, which roughly translates to which nested drift shell a given electron is confined to. (From this point on we will refer to L-parameters, rather than $L*$⁠, for simplicity in notation). The drift shells a given L correspond to will depend on the specific magnetic field model used, and multiple methods of defining L exist depending on both the magnetic field and the approximations used to specify the drift shells (see, e.g., Refs. 19, 20, 25, and 26). If one considers the L parameter as a proxy for the relative radius of each drift shell, it becomes clear why diffusion across the third adiabatic invariant is called radial diffusion: particles diffusing to different L values are diffusing onto drift shells closer to, or further from the Earth. On the whole, radial diffusion is related to the large scale movement of particles toward and away from the Earth, acting to reduce gradients in the phase space density profile. The inward diffusion also corresponds to an increase in energy, if the first adiabatic invariant is conserved. Meanwhile, smaller scale processes breaking the first and second adiabatic invariants can result in local acceleration of particles (e.g., Refs. 27 and 28).

The diffusion coefficient $DLL$ contains the combined effect of electromagnetic perturbations on timescales corresponding to the electron drift. These perturbations consequently break the third adiabatic invariant, causing the diffusion to nearby $Φ$ (or, equivalently, L). A comprehensive review of the current state of knowledge of radial diffusion coefficients can be found in Ref. 29. The first attempts to quantify this into a diffusion coefficient $DLL$ assumed that such perturbations were stochastic; small scale continual ripples in the magnetosphere. The contributions from magnetic and electric potential perturbations were considered separately, using asymmetric and symmetric perturbations from a simple magnetic field model.³⁰ Unfortunately these theoretical diffusion coefficients are problematic to apply; the magnetosphere is significantly more dynamic, rendering these assumptions invalid, while in practice one cannot observe these quantities to estimate them more accurately. There exists a gap between the theory and application of radial diffusion; accurate diffusion coefficients would require knowledge of the entire magnetosphere at each time step to be able to calculate electron drift paths. Theoretical approaches must use magnetic field models and tend to focus on determining the validity of the underlying assumptions in order to derive a more appropriate method of calculating diffusion coefficients.^26,31,32 Any estimations of $DLL$ used in modeling must make a significant number of approximations, most often based on the techniques in Refs. 33 and 34 (particularly operational models).

Reference 34 derives the diffusion coefficients for a particular storm using a formalism based on:³⁰ for a given magnetic field model, find the deviation $dLdt$ from drift contours due to azimuthally symmetric and asymmetric electromagnetic perturbations, and from this find $⟨(ΔL)2⟩$ and hence $DLL.$ Ref. 34 uses a compressed dipole magnetic field, using an asymmetry factor $ΔB$⁠.Ref. 34 also splits the diffusion coefficients into diffusion due to magnetic and electric perturbations, rather than into perturbations from induced electric and electric potential fields.Ref. 34 used the larger component of the two to avoid counting the effect of the same perturbation twice, but the component which was dominant during their case study is not always the strongest one.^35,36 In many subsequent studies, these components are simply added together to estimate $DLL$⁠. This makes the diffusion coefficients much easier to find and apply, but the resultant double-counting of perturbations means the final values could be off by around a factor of 2.³¹ However, uncertainty in $DLL$ spans at least one and sometimes up to four orders of magnitude, e.g., across similar time periods,³⁶ between observations and models (Ref. 35, Fig. 6), and across different models at the same Kp and L (Ref. 3, Fig. 6) (Ref. 4, Fig. 4). Given this uncertainty, and also the uncertainty in the relative size of induced electric and electric potential components, double-counting is considered an acceptable compromise. Other choices made by Ref. 34 based on properties of their storm-time magnetosphere are also replicated in subsequent methods, simply to have any ability to model the radiation belts at all. This may be why so many of these methods perform similarly - well on average, but with little specificity. A list of sources of uncertainty in $DLL$ can be found in Ref. 37.

Here, we pause to explicitly identify research goals. Having a separate section for this allows us to report back on unfruitful research avenues and to directly label our goals to conclusion in the discussion, making the logic easy to navigate across several key results. No research paper has only a single important result, and we hope this will bring out subtleties that may otherwise be missed.

We split these research goals up into:

1. Primary goals (research questions identified when scoping out the project and applying for funding)

2. Secondary goals (research questions that arose when developing our methodology)

3. Additional goals (further questions that arose as part of the analysis which we realized can answer as part of this work)

We choose to present our work in this way as it reflects our cyclic methodology far more accurately, where we constantly revisited concepts and experiments until we reached a coherent structure explaining our results.

[P1] Benchmarking for ensembles.

• How does one analyze ensembles for radiation belt modeling? How can we quantify differences between ensemble runs?

• Do changes in initial condition affect the differences between ensemble members? If so, which aspects of the initial condition have the strongest impacts?

• Do changes in model settings (such as $Louter$⁠) affect the differences between ensemble members?

[P2] What is the timescale of radial diffusion?

• What is a useful definition of “timescale”?

• What is the general radial diffusion timescale?

• Does timescale vary with initial conditions?

[S1] What do we do if we are not using a data-driven outer boundary?

• What boundary conditions would be physical to use? (rather than what is convenient for our observations)

• Where there are multiple potentially physical boundary conditions, how do they affect the timescale and evolution of radial diffusion?

• What boundary conditions can we use in practice to balance physical boundaries (S1a) with observational and operational modeling constraints?

• How does the choice of outer boundary condition and location interact with the current uncertainty on the outer boundary of radiation belt models?

• How does the outer boundary location interact with the increased diffusion coefficient at high L?

[S2] What happens to an enhancement from local acceleration under radial diffusion (rather than just inwards diffusion of the “background” distribution, i.e., a source at high L?)

[S3] How can time to monotonicity/morphology of the particle distribution be used to analyze ensemble members?

[S4] Do we need to include loss from precipitation via pitch-angle scattering to represent radial diffusion?

• [A1] What analytical methods can be adapted to understand radial diffusion?

• [A2] How important are PSD gradients vs the L-dependence of $DLL$?

• [A3] Is diffusion limited by the smallest value of $DLL$ in the domain, i.e., the diffusion coefficient at lower L?

In Sec. VI, we will discuss our findings for each of these questions and outline future questions that are unanswered.

In this Sec., we will outline the scheme used for the numerical experiments in our ensemble and the metrics we use to analyze our experiments. We choose not to compare variability across the final phase space density distributions after a set time period, but to compare how several useful quantities vary across ensemble members. We use a proxy (time to monotonicity, $tm$⁠) that indicates when radial diffusion has finished changing the shape of the distribution, and mass-like and energy-like techniques from analysis of dynamical systems to understand the ongoing evolution of the distribution. These tools allow us to see how radial diffusion is still contributing to radiation belt dynamics.

In Sec. IV A, we discuss our equation for the initial condition, the parameters we will vary in our ensemble and the diffusion model we use. In Sec. IV B, we will motivate time to monotonicity as a metric for timescale. In Sec. IV C, we derive the energy density-like and mass-like quantities we use to verify and interpret our simulations.

Results from our investigation can be found in Sec. V and are brought into context with current understanding in Sec. VI.

Following initial experiments, we investigated whether our simple model required additional terms to validly represent the physics of the outer radiation belt. An important loss mechanism in the outer radiation belt is pitch-angle scattering (e.g., Ref. 44) where the interaction between electromagnetic waves and electrons results in changes to the velocity vector relative to the magnetic field direction. In the relatively dense region of the plasmasphere, observed plasma and whistler-mode wave conditions are such that pitch-angle scattering is enhanced for high-energy electrons (e.g., Ref. 45).

Reference 49 shows that this scheme is unconditionally stable in the case where diffusion coefficients vary in time, but does not consider where the diffusion coefficients vary in space, as the coefficient matrix is then no longer symmetric. However, verification for our model can be found in the supplementary material of Ref. 5 in addition to the initial SpacePy verification of Ref. 48. The time step and spatial resolution of the simulations are 1 s and 0.1L, respectively.

This reflects typical phase space densities using a peak and step, which represents a state where inward radial diffusion has been occurring for some time (the step) and an enhancement of locally energized particles (the Gaussian). These distributions were chose to reflect the phase space density distributions observed in Ref. 50. Individual parameters are

• A amplitude

• B step size (strength of error function)

• $σ$ width of density peak

• $μ$ location of phase space density peak in L-space

Demonstrations of these initial settings can be found in Fig. 1. Note that $B=0$ would result in an initial PSD distribution of only a Gaussian, while increasing B represents a decreasing difference between height of the peak and height of the step. $B=1$ results in a step of comparable size to the enhancement, where the PSD at high L is the same or greater than the pre-peak PSD. Although at $B=1$ the plateau asymptotes to the height of the Gaussian, since the two terms are summed together this is not a flat plateau. Our default settings are $A=9×104,B=0.05,μ=4,σ=0.38,γ=5$ and $Louter=6.$ The inner boundary is always set at $L=2.5$⁠. We vary the location of the outer boundary (⁠ $Louter$⁠) to reflect the fact that radial diffusion models often have different $Louter$ and to investigate the impact of different outer boundary locations when we know that diffusion varies in space as well as time.

The model in Ref. 5 uses a constant value (Dirichlet) inner boundary to characterize the inner edge of the radiation belt, where particles are lost to the atmosphere. Choices of the outer boundary are less clear. Physically, one expects the PSD to be smooth across the boundary; hence, a Neumann (constant gradient/zero flux) boundary may be appropriate. By default, the model uses a constant gradient (Neumann) outer boundary, set at zero, which matches the physics of a slow injection from high-L as represented in our initial condition. However, a Dirichlet boundary allows the modeler to input PSD values, and indeed large-scale models tend to interpolate time-varying outer boundary values from available data. Neither of these methods are designed to use the true edge of the outer radiation belt, which varies considerably when an outer edge is distinct enough to observe at all.⁵¹ Since both methods should be physically appropriate for the underlying plasma, in our results we investigate the impact of both choices of outer boundary condition. We will review this outer boundary in the discussion, following our results [S1a].

Ensemble modeling for predictive purpose use the variability across model runs in the given ensemble, using a given error metric or loss function. We are also interested in examining timescale of radial diffusion, for which we need a quantitative measure.

However, error metrics are not an ideal tool for comparing the evolution of phase space density. They do not tell us about the evolution of the system state, or properties of that state we are interested in. Therefore, one of the secondary goals of this work was to select and investigate potential metrics. Initially, simple error metrics such as mean square error (MSE) were selected. MSE would quantify a scalar difference between distributions, which would be useful for analysis of variation and uncertainty across ensembles. However, the significant variation in scales covered in this problem make it difficult to generalize or compare different cases. We found this regardless of using linear- or log-based scales; thresholds usually would need to be specified for when radial diffusion was “finished enough” or for when two distributions were “similar enough,” and the results became dependent on that choice of threshold. Instead, our experiments were analyzed with a property that captured the physics of the system we were interested in: time to monotonicity, $tm$⁠. This choice of metric is motivated below. Requirements of the metric used for analysis are [A1]:

• Robustness. The metric used must be insensitive to any thresholds used. It must, therefore, be scale independent (i.e., work across multiple orders of magnitude, because of Kp dependence)

• Interpretable The metric must aid in understanding the system.

• Radiation belt system specific. The metric must be related to radiation belt modeling; it should provide insight into either the system state, or specific properties related to physical processes.

• Time-series informative. The metric must enable analysis of the evolution of system.

Initially, potential measures were tested, such as the decreasing MSE between distributions at each time step, the evolution of maximum gradients, the total area, and the proportional change in amplitude. All these required an arbitrary threshold to determine when a given experiment was “finished,” the choice of which strongly impacted the time until distributions became similar, particularly for larger Kp. For example, MSE-based metrics varied by orders of magnitude depending on several model choices. No MSE-based metric could be found that worked robustly. Furthermore, many of these metrics ended up being dominated by the inner boundary condition, whereas we wanted to know about the evolution of the entire phase space density distribution.

To find an appropriate metric of the how the ongoing radial diffusion affected the evolution of the system, we turned to properties of the PSD distribution under radial diffusion. Figure 2(a) shows the phase space density profiles we expect on the timescale of radial diffusion. A quiet time, or background distribution, is shown in solid black lines: a low level, high L source of particles feeds the system from constant substorms. Due to the L dependence of the diffusion, the drop-off at low L is quite sharp. In dashed lines, an enhancement of particles due to local acceleration is shown at an intermediate L. Following a single enhancement, we expect the distribution to gradually return to a monotonic state via radial diffusion.^27,28

Therefore, since diffusion acts to even out gradients in the PSD, we are more interested in changes in the shape of the PSD distribution than the total PSD (i.e., the integrated area under the curve). Hence, we also include the following criteria:

• Insensitive to total particle population. The metric must be insensitive to shifting the distribution up and down the y-axis.

We chose to work with time to monotonicity, $tm$⁠. When the distribution has become monotonic, the PSD distribution no longer has a peak for radial diffusion to smooth out. $tm$⁠, therefore, indicates when the PSD distribution has stopped changing shape, when radial diffusion is no longer changing the properties of the radiation belts. $tm$ is a proxy for whether radial diffusion is still significantly affecting the evolution of the particle distribution. The potential monotonic distributions are shown in Fig. 2(b); the “background” distribution which is nonzero but unchanging in time (i.e., the influx of particles balanced by constant movement of particles inwards), a shrinking version of this when more particles are lost than enter the domain, and finally a zero profile.

In the results section, we use $tm$ to explore how idealized radial diffusion models vary when changing initial settings. $tm$ represents our physical expectations; our intuition that after a localized enhancement, the PSD will eventually relax to monotonic distribution. We can test our expectations and how the changing parameters affect these expectations. See Sec. VI for evaluation of our chosen measure and for alternative approaches.

In this section, we derive the terms that comprise the changing energy (⁠ $Et$⁠) and mass (⁠ $Nt$⁠). The figures comparing these terms, evaluated explicitly for each of our boundary conditions, can be found in Sec. V where they best support our analysis.

The third term in $Et$ is one that will give us insight into how the dynamics of the distribution will evolve to minimize our energy. This term has an integral, which depends on the square of the derivative in L. As this is a non-negative quantity, this ensures that the integral contributes to energy loss (as it is preceded by a negative sign) while gradients in the distribution function exist, up to a point where the contributions from the boundary conditions balance this out. This is an property identical to Cahn–Hilliard diffusive dynamics, where diffusive terms correspond to energetics that penalize the formation of (sharp) gradients.^53,54 This penalty is enhanced for gradients occurring at higher values of L, as their contribution to the integral will be increased by a positive power of L, suggesting that the distribution function can minimize this integral by moving gradients to lower values of L. Thus, $Et$ decreases toward a long-term solution through diffusion to reduce gradients $(fL)2$ and through the population moving toward lower L, and this movement of the population will be increased for higher powers n.

Our methodology was to begin with $tm$⁠, to find out how long it takes a given initial distribution to reach a monotonic state, and how this varies with initial conditions. We first compared $tm$ with Kp for each parameter in Eq. (10). For $tm$⁠, we note that from Eq. (2), we expect timescale to vary significantly with Kp. Kp is a proxy for strength of radial diffusion, even though $Kp>6$ is unlikely. Since we expect time to monotonicity to depend on Kp we primarily use a heatmap for the ensemble used to investigate each parameter, demonstrating how $tm$ varies with each (Kp, parameter) pair. To aid understanding, these results are also presented in an alternate format, where the $tm$ for each Kp are plotted as a line. Each experiment ran for a week; model runs where monotonicity were not reached are left empty.

Following our initial $tm$ analysis, we then investigated any noteworthy results, for example by looking at the evolution of the phase space density of a specific simulation. In general, runs with a Neumann outer boundary condition (zero flux) are shown on the left, while the right hand column corresponds to a Dirichlet outer boundary condition (constant value).

Finally, we incorporated our analysis from Sec. IV C for each parameter to understand and generalize any patterns we saw. Our $N,E$ experiments are also run for a week, using $n=6$ in our diffusion coefficient. Just as for our $tm$ analysis, we calculated $N,E$ for each returned time step (every 6 h).

Each parameter in the initial condition was systematically investigated. Selected results are presented in the main text in an order chosen to best convey our conclusions; the full set of individual $tm,N$⁠, and $E$ experimental figures can be consulted in the supplementary material, labeled as Figure. SX. As a one dimensional simulation, the number of particles in a given phase space bin (PSD) is in units of $(m−2s)1$ [SI units, per (unit speed $×$ unit length)]. This can be scaled to any of the other units systems used for PSD elsewhere.

Overall, it is clear from Fig. S1 that generally, more runs with a Neumann (fixed gradient) outer boundary reach a monotonic state within a week, compared to runs with a Dirichlet (fixed value) outer boundary. We will investigate the two outer boundary conditions before examining the effect of each parameter in the initial condition. To understand the evolution of f in these experiments, in Fig. 3, we show the phase space density for both Neumann and Dirichlet outer boundary conditions, with an outer boundary located at $Louter=6.5$⁠. We use default settings for the initial PSD, a Kp of 4 and run for a week. Panels (a) and (b) show the heatmap, while waterfall plots are in panels (c) and (d).

In both experiments, the peak remains high and moves inwards. However, there is a difference in the outer part of the simulation, as demonstrated in Fig. 4. The plateau becomes significant for the Neumann boundary but not for the Dirichlet boundary. This makes sense as the outer boundary value is able to rise for the Neumann case, reflecting outward radial diffusion. For Dirichlet runs, the outer boundary value cannot rise, but particles can be lost. Note that the amplitude of the Neumann peak is still comparable to the Dirichlet case (4.55 and 4.52  $log10(PSD)$ at $L=3.9$ and 3.7, respectively); it is the plateau that has changed. More Neumann experiments reach monotonicity as the plateau can rise instead of having to wait until the peak diffuses completely inwards. This inability for the high-L (to the right of the peak) PSD to reach (positive) monotonicity independently of the left part is one reason why it takes longer for the Dirichlet runs to reach monotonicity. This corresponds to outward radial diffusion varying depending on the outer boundary condition [S1b].

For a more realistic outer boundary condition, we need to consider these, along with the fact that we do not currently have a clear outer boundary location. We discuss all these factors together in Sec. VI D.

We can also examine how the evolution of our mass- and energy-like densities varies with Neumann and Dirichlet outer boundary conditions. Figure 5 shows $N,E$ across the week, plus the components of $Nt,Et$⁠. Panels (b) and (d) show the absolute value of components on a log scale to make order of magnitude comparisons earlier; in rare cases in later analyses where the change becomes positive (i.e., a gain in mass or energy, rather than a loss), this is specified.

By definition, for Neumann experiments $Nt(1)=0$⁠, mass is only lost through the inner boundary. For the Dirichlet experiments, it is clear that the outer boundary dominates loss, with $Nt(1)$ up to two orders of magnitude larger than $Nt(2)$⁠; hence, the number of particles decreases more rapidly. We find that less than 0.8% of the original mass is lost with a Neumann outer boundary, while around 24% of the mass is lost with a Dirichlet boundary [S1b].

Comparing the individual terms contributing to changes in mass in Fig. 5(b), we see that loss from the inner boundary $N(2)$ is the same across the week regardless of the outer boundary (and therefore independent of the different interior distribution as the system evolves). For inner boundary loss to be big enough to vary, one must run experiments with very large diffusion coefficients [e.g., using $Kp=9$ in Eq. (2)] or for a much longer time.

The results for our $L2$ norm $E$ are somewhat counterintuitive. In total, we know that experiments with a Neumann outer boundary can eventually reach a lower- $E$ state (zero everywhere) than experiments with a fixed outer boundary value, where the minimum-energy state will have the same fixed value at $f|OB$ as the initial condition. However, we find that Neumann simulations appear to be reaching a limiting state, where $Et$ is increasingly smaller and $E$ relatively unchanged.

For a Neumann outer boundary, $Et=−Et(3)$⁠. While the Dirichlet outer boundary can contribute to changing energy (⁠ $Et=−Et(1)−Et(3)$⁠), we can see from Fig. 5(d) that the dominant mechanism for energy loss is mostly from the reconfiguration term $E(3)$⁠. However, this term accounts for more energy loss when the outer boundary is Dirichlet rather than Neumann, because there are more steeper gradients when the outer boundary is fixed. Since $Et$ depends on gradients (⁠ $fL$⁠), $Et$ is larger for Dirichlet runs as there are gradients both sides of the peak, rather than just to a plateau. For Neumann experiments, the gradients rapidly flatten into a plateau at higher L. Remember that our equation for $E$ has a factor of $1L2$ in it: the same PSD at a higher L contributes less to the norm, because it can be moved around more easily. Hence, Neumann runs have a lower $Et(3)$⁠. Once the plateau has been reached, the only way to lose energy is for material to move down the gradient (and then out through the inner boundary). This process is slow and so $E$ is effectively limited. On the other hand, the Dirichlet experiments are more effectively moving material to higher L, and then out of the domain. Having an outer boundary that allows flux in/out means that $L2$ norm is reducing more quickly than when we allow the PSD at the outer boundary to change. Hence, Neumann appears to be reaching a limiting state first, where $Et=0$ and $E$ is unchanging, even though we know that if left forever, it can reach a lower-energy state than Dirichlet [S1b;S1e].

The results of systematically investigating each parameter, using first time to monotonicity $tm$ and then our quantities $N,E,Nt,Et$⁠, are shown here. Those properties we have considered to have a significant impact on time to monotonicity are presented in detail, while remaining properties are covered briefly in Sec. V A 3 [P1b;S3a].

Step size B: The step size B corresponds to a system where the phase space density is higher further out in the radiation belts; i.e., a situation where inwards radial diffusion from a distant source has already occurred. A higher step, therefore, corresponds to radiation belts that have more material before the enhancement occurs.

Figures 6(a) and 6(b) show how time to monotonicity $tm$ varies with both Kp and increasing step size. Blank values indicate monotonicity was not reached within a week, the length of the experiments. When starting with a larger step size, we find that monotonicity is reached sooner for both outer boundary conditions. This behavior is as expected as with an increased B, and f is already closer to a monotonic state. Again, more ensemble runs with Neumann outer boundary reach monotonicity. Looking at this information in the alternative format in panels (c) and (d), we see that for both outer boundary conditions, $tm$ appears to increase exponentially for smaller step sizes.

Our $N,E$ results show that the evolution is not necessarily straightforward, however. Figure 6(e) shows the total mass in f across the week simulated, for the larger step size $B=2$ for both Neumann and Dirichlet outer boundary conditions (the comparison against the default value can be found in Figs. S3 and S4). As expected, the total number of particles changes very little for a zero flux (Neumann) outer boundary. However, the overall mass response changes when mass can flow across the outer boundary (Dirichlet experiments). We see that for a larger step size $B=2$⁠, the experiment starts to gain mass at around 80 h. From the mass change terms $Nt$ in Fig. 6(f), this is clearly from the outer boundary, when the distribution drops below the fixed outer boundary point $f|OB$⁠. At this point, the gradient $fL|OB$ will become positive, and material will flow into the domain. Inner boundary loss varies with B but is again independent of the outer boundary condition.

The $L2$ norm $E$ is much higher for a larger step, and the norm reduces over several orders of magnitude [see Fig. 6(g)]. The Dirichlet run started out losing more energy than Neumann (as we also see using default settings above) but later in the week, energy loss drops off and the Neumann case has lower energy (and is, therefore, closer to a point where the dynamics have stopped changing). This is because there is an increase in the norm (⁠ $E$⁠) with the reversed outer boundary flow; however, the corresponding energy change term $Et(1)$ reaches a comparable magnitude to the dominant reconfiguration term $Et(3)$⁠. Therefore, the total change $Et$ for the Dirichlet case with $B=2$ becomes very small, while the Neumann case is still reducing in norm because the PSD can be reconfigured (diffused).

Physically, this means that with a higher step, we are finding that the Neumann case reaches a lower energy state by the end of the week than the Dirichlet experiment, unlike our default settings Fig. 5. For the Dirichlet case, the constant outer boundary value is higher than the peak, allowing material to come in through the outer boundary. This experiment is in a state where the main dynamic is a constant churn of mass coming in and then being diffused to lower L to reduce the entire distribution to a lower energy state. Physically, this corresponds to an infinite source at the outer boundary if we were to run this indefinitely. See Sec. VI for our overall conclusions on more suitable outer boundary settings.

Enhancement location $μ$: The Gaussian in the first term of Eq. (10) corresponds to an enhancement, and $μ$ corresponds to the location of this enhancement in L.

Considering $tm$ across a variety of enhancement locations in Fig. 7, we find that monotonicity is reached more quickly for higher $μ$ for both outer boundary conditions, although the shape of the dependence is seen to be very different in Figs. 7(c) and 7(d). A sooner $tm$ for higher $μ$ makes sense as $DLL$ will be larger at higher L, so when the peak is located at higher L, diffusion to flatten this peak happens more quickly. Again, far more runs reach monotonicity with a Neumann outer boundary.

Figures 7(a) and 7(b) show us the general relationship between each parameter and $tm$⁠, while (c) and (d) show us the specifics of each relationship. We see that the relationships between $μ$ and $tm$ for each Kp are quite different for the different outer boundary conditions. This disparity could be due to the different mechanism to the right of the enhancement; a zero flux outer boundary condition allows the plateau to rise and reach monotonicity to the right of the peak, while the same region with a constant value outer boundary cannot reach monotonicity until the enhancement has completely diffused, although material can be lost through the boundary.

The bottom two rows of Fig. 7 compare the evolution of $N,E$ across the simulation for the default value of $μ=4$ to a more distant peak at $μ=5$⁠. The total number of particles $N$ is higher when $μ$ is lower, which is as expected since the step extends further inwards. There is little change in mass with $μ$ for Neumann runs, also as expected. There appears to be different amounts of mass lost in Dirichlet simulations, so we consider the outer and inner boundary particles losses $Nt(1)$ and $Nt(2)$ in (f). These are always negative (particles are only lost, not gained) but look quite different. The outer boundary loss dominates for both values of $μ$⁠. With a higher-L enhancement, the inner boundary loss is less. Therefore, while the outer boundary loss is comparable, the higher the $μ$⁠, the more strongly that outer boundary loss dominates over the loss from the inner boundary. Figure 7(g) shows the evolution of $E$⁠. Lower values of $μ$ (enhancements at lower L) actually result in higher $L2$ norms, because you have more PSD total (for the same reason as the mass above) and more of this mass is at lower L. In both cases, the Neumann experiments reach a configuration where $Et$ is very small and $E$ stops changing. The Dirichlet experiment with a higher-L enhancement rapidly loses $E$ but by the end of the week, is no longer changing much. In (h), we can see the terms of $Et$ for each experiment. For readibility, we show only the $Et$ terms for $μ=5$ here. The reconfiguration energy loss (⁠ $Et(3)$⁠) evens out more quickly for Neumann experiment with $μ=5$ than with the standard initial condition, presumably because it is easier for the step to rise up and plateau in a monotonic state. For the Dirichlet experiment, the reconfiguration term $Et(3)$ dominates, but rapidly drops off until it is comparable with energy loss at the outer boundary.

We find that an initial distribution with more material at high L (i.e., step size B) and with an enhancement at high L (i.e., $μ$⁠) diffuse more quickly. B and $μ$ are the most significant initial conditions, yet the specific evolution of the system varies depending on interaction with boundary conditions [P1b].

Enhancement width $σ$: Figures 8(a) and 8(b) show that for both Neumann and Dirichlet boundaries, $tm$ is reached more quickly for narrower peaks, i.e., for smaller values of $σ$⁠. This difference is only slight. There are two aspects at work here; a wider $σ$ will have access to larger diffusion rates at high L, but the gradient $fL$ will be less steep. We can examine these components using $N$ and $E$ to determine which is more significant.

An enhancement across more L has somewhat more mass and appears to lose mass more quickly for both Neumann and Dirichlet outer boundary conditions. Despite the fact that experiments with a larger $σ$ also start with higher $N$⁠, at the end of the week, 98% and 75% of the mass remains from the initial population (for Neumann and Dirichlet experiments, respectively), compared to 99% and 76% remaining from our default experiments. Both experiments with a wider enhancement lose more from the inner boundary [Fig. 8(d)], again showing that the initial condition controls more loss from the inner boundary than changes in the outer boundary condition. Loss from the outer boundary is slightly more with a larger $σ$⁠, but is roughly comparable.

Total $E$ is higher for a larger $σ$⁠, which makes sense as the experiment has more mass and hence larger $∫f2dL$⁠; the distribution f is further from a steady state. All $Et(1,3)$ with $σ=0.5$ are slightly larger, but overall very similar. Unsurprisingly, with similar $Et$ but a larger starting $E$⁠, experiments starting with a larger $σ$ end up retaining proportionally more of the initial energy (37% and 58% for Dirichlet and Neumann, respectively, rather than 30% and 47% of the initial energy when using default $σ$⁠). Overall, with a wider $σ$⁠, the Dirichlet experiment loses more $E$⁠, even though $tm$ is slower. The results from the investigation of $σ$ indicate that the trade-off between gradient in the PSD and the L-dependence of $DLL$ is subtle and nuanced (a wider enhancement has a less steep gradient but samples higher $DLL$⁠).

The $σ$ results suggest that we should compare the role of the spatial (i.e., L) dependence and the gradients in the distribution function on the overall amount of diffusion. We consider their role in reconfiguration term $Et(3)$⁠, since this generally dominates $Et$⁠. With a constant diffusion coefficient, only the gradient term would contribute to the $E(3)$⁠. With an L-dependent $DLL=D0Ln$⁠, both $(fL)2$ and $D0Ln−2$ will contribute to $E(3)$⁠. We simply compare the order of magnitude of these components via the ratio' of $(fL)2$ to $D0Ln−2$ for $n=6$⁠, shown in Fig. 9. Overwhelmingly, it is the gradient component $(fL)2$ that dominates. This is unsurprising once one considers that $D0∼10−10$ s⁻¹.

Throughout this analysis, we have found that one must consider the whole domain; for example, the amount of diffusion is not limited by the smallest $DLL$ but also the shape of the distribution, loss, the choice of domain, etc. This is because $E(3)$ is the dominant component of $Et$⁠, which determines the PSD evolution. $E(3)$ is an integral over the entire simulation domain.

Figure 9 indicates that the gradients have more impact than the L-dependence of the diffusion coefficient. However, with a longer spatial (L) domain, this will begin to change, particularly for a Neumann (zero flux) outer boundary, where there are fewer gradients. Using an idealized $DLL$⁠, gradients almost always dominate over the effect of $DLL$ increasing with L. This will be why wider enhancements (larger $σ$⁠) take longer to reach monotonicity when using our operational (Ozeke) $DLL$⁠: there are consequently shallower gradients [S1e].

We find that the PSD gradient $fL$ contributes more to the evolution of the system than the diffusion coefficient $DLL=D0L6$⁠. Note that this idealized scaling was chosen to be comparable with the Ozeke diffusion coefficient at a given L and Kp; there are many other models, often more sophisticated, yet all have a significant L-dependence. This is discussed further in Sec. VI E [A2].

In this experiment, we varied the domain for the simulation to see what difference it made. A Dirichlet (fixed value) condition is used in the majority of operational radiation belt models, to reflect observations e.g., Refs. 24, 55, and 56. The simulation domain is curtailed to the location of the spacecraft; different $Louter$ values then correspond to using data from different spacecraft missions to set this outer boundary. Both types of outer boundary condition are investigated in this phase of experiments, and the Dirichlet experiments retain the outer boundary value fixed in the initial condition.

For a Neumann outer boundary condition, a more distant outer boundary (larger $Louter$⁠) took longer to reach monotonicity; i.e., a smaller domain reached monotonicity more quickly [Fig. 10(a)]. This suggests that the choice of outer boundary location changes the shape of the PSD distribution, especially the height of the plateau. For Dirichlet conditions, $tm$ was independent of domain size.

To investigate this, we compare two Neumann runs where we vary the outer boundary location to be $Louter=6.5$ and 7.5. (Equivalent plots for Dirichlet can be found in the supplementary material, Fig. S9). Figure 11 shows the PSD distribution at eight equally distant times throughout the week, with Kp = 4 and using Ozeke $DLL$⁠.

The difference in PSD between peak, and plateau edge (both indicated with dotted lines) by the end of the week is much larger for the run with a wider domain, in Fig. 11(b). The effect of this can be seen more clearly by considering the PSD at the final time step, in Fig. 11(c). The run with $Louter=6.5$ is more close to monotonic because the value of the outer boundary has raised higher. Despite the larger $DLL$ at high L, more reconfiguration of the distribution is needed for a longer domain governed by the same equation and as was seen in Sec. V A 4, the gradients are still more important than the $DLL$ up to $L=6.5$⁠. Although the material in the peak being diffused outwards can be spread across more L when there is a more distant $Louter$ and large $DLL$ values at high L encourage this, the material has to travel farther before the plateau rises and monotonicity is reached. The Neumann $tm$ dependence on domain (i.e., on $Louter$ arises because the shape of the distribution varies with changes in the simulation domain).

Analysis of $N$ in Fig. 12(b) indicates that neither the outer boundary location or condition affects loss from the inner boundary. For Dirichlet experiments with varying $Louter$⁠, mass loss from the outer boundary is of comparable order within a few hours, regardless of where that outer boundary is. This corroborates findings in Sec. V A 4 that the gradients in the PSD distribution dominate diffusion over the entire domain, rather than higher diffusion coefficients located in one region. While the extra mass at time $t=0$ was obvious for a longer domain, the change in initial $E$ is negligible, as can be seen in Fig. 12(c). However, the evolution of the $L2$ norm is nuanced; despite having $Et$ terms of similar order [Fig. 12(d)], with a longer domain, the Neumann experiments reaches a lower level of $E$⁠, while the Dirichlet experiment has a greater value of $E$⁠.

The mechanism behind these results is, unsurprisingly, the rising plateau for Neumann and the outer boundary flux for Dirichlet experiments. In order for Neumann experiments to reach a configuration of lower $E$⁠, the high-L plateau rises. With a longer domain, this means that there are more particles at high L; hence, $E$ can be lower than with the same number of particles at a lower L. Additionally, $E(3)$ is slightly larger for $Louter=7.5$ than 6.5, which can be attributed to outward radial diffusion due to the longer domain and higher $DLL$ at $L=7.5$⁠. For the Dirichlet case, the reconfiguration energy change is almost exactly the same. However, more is lost to the outer boundary for $Louter=7.5$ than 6.5. Even though they lose roughly the same each hour after 70 h, this is enough to make $E$ slightly lower for $Louter=6.5$

Using $N,E$⁠, we find that the choice of outer boundary location changes the shape of the PSD distribution for Neumann; in exactly the same manner, $tm$ varies depending on the outer boundary location. For different outer boundary conditions, a different outer boundary location could result in heading faster or slower to a state where the dynamics are minimally changing (i.e., to minimum $E$⁠). In Sec. VI, we discuss outer boundary choices, including the applicability of using a Dirichlet outer boundary that is fixed, but not using observations [P1c].

Loss from pitch angle scattering is significant and should be included to see how it relates to the timescale of radial diffusion. We include this loss by modeling the electron lifetime.

In general, with loss it is no longer always true that more Neumann experiments reach monotonicity; nevertheless, they still have shorter $tm$ than Dirichlet experiments. All the $tm$ plots can be found in supplementary material; again we select the results that inform us about the overall pattern.

Figure 13 show experiments with pitch angle loss for the default initial condition. Neumann and Dirichlet runs look very similar, suggesting that pitch angle loss may control more of the dynamics than the outer boundary condition. The analytic quantities $N,E$ for these runs confirms this; Loss from the pitch angle scattering approximation dominates over loss from the outer or inner boundary Fig. 14(b), and the evolution of $E$ with loss looks similar regardless of outer boundary condition, unlike the runs without loss [Fig. 14(a)]. The experiments including pitch angle loss are heading much more quickly toward a steady state, and after around a hundred hours the experiments are not changing much in $E$⁠, as $Et∼0$⁠.

$tm$ is more complex to analyze. Over the span of the week, enough mass is lost that the Dirichlet run begins to gain mass. Just as without loss, particles entering the domain are contributing to an increase in $E$⁠, which results in a final $E$ that is higher for Dirichlet than Neumann. As a result, the Neumann case is closer to a steady state. Loss from pitch angle scattering effectively mimics the reduction in PSD that would occur over a very long timescale of radial diffusion, as it is strongest near the edge of the plasmapause (which is located around the bulk of the enhancement). With a constant inflow of particles, Dirichlet runs have a higher $E$ but do not come closer to monotonicity. Electron lifetime loss is creating a new local minimum in the PSD, and so the distribution is not monotonic. To demonstrate what is happening in Fig. 14 for both Neumann and Dirichlet runs, the diagram in Fig. 15(a) shows an example phase space density distribution with this additional minimum. This physical profile is corroborated by the fact that $Etloss>Et(3)$ and therefore loss dominates over reconfiguration in the system evolution. We will expand on the consequences of this below. Finally, there is a difference in inner boundary flux; all experiments including lifetime loss have the same, lower flux [S1b;S4].

Although experiments with loss $L$ that reach monotonicity do so quicker than they did without $L$⁠, for all parameters there are several initial values that reached monotonicity without loss but no longer do once loss is included, usually at lower Kp (i.e., weaker radial diffusion relative to the same loss). The reverse is rarely true; only for narrow $σ$ (e.g., 0.2, 0.25) and $Louter=5.0$ is a $tm$ found with $L$ where none was found before. In this section, we explain the physical mechanism behind this general pattern; again, all individual $tm$ plots can be found in the supplementary material, Fig. S5.

The existence of experimental setups that will never reach monotonicity is most clearly demonstrated by following a case where $Lp>μ$ [Fig. 15(b)]. There will be loss everywhere left of $Lp$⁠. In the region $μ<L<Lp$⁠, if loss $L$ is strong enough, it can work against monotonicity by creating a minimum. The second row of Fig. 15 demonstrates how this case can evolve for Neumann and Dirichlet outer boundary conditions, using default initial conditions, the same loss with the default plasmapause at $Lp=5$ and $Kp=4$ in the diffusion coefficients. The Neumann and Dirichlet [Figs. 15(c) and 15(d), respectively] experiments show that even when the distribution is much reduced (over 24 and 72 h, respectively), it is not monotonic [Note that where the loss does not dominate over radial diffusion, then instead the effect of the loss is for the overall distribution to reduce more quickly, but still maintain the characteristic diffusion distribution (e.g., the intermediate distribution in Fig. 2(b))] [S3a;S3b].

In general, initial conditions impact the diffusion in a similar manner to without loss, for example, $tm$ varies significantly with $μ,B$⁠. In this section, variation of each parameter is compared to the case without loss. Again, we select the most significant results here, while all the figures can be found in the supplementary material Figs. S5–S8. In Sec. V A, each $N,E$ was normalized using initial values $N0,E0$ for the phase space density using all default parameter values for the initial condition. To compare the effect of loss, we instead choose normalization values $N0,E0$ from the initial phase space density of the higher parameter value in each experiment pair, e.g., we normalize using the initial mass and energy density with $B=5$ rather than $B=2$⁠, $μ=5$ rather than $μ=4$⁠, etc. This normalization was chosen to ensure that the effect of each parameter is extracted, rather than repeating experiments comparing the default initial condition with and without loss, which was explored in Secs. V B 1 and V B 2 [P1b;S3a;S4].

Step B: Just as without loss, a higher step size B results in a shorter $tm$ as the distribution is already closer to monotonic. However, this effect is much smaller. As can be seen in Figs. 16(a) and 16(b), the $tm$ relationship to changing B is still exponential, but stops changing once $B≤2$⁠. In fact, a step this large quickly results in gains from the outer boundary; in Figs. 16(d) and 16(f), $N(1)$ and $E(1)$ for experiments with loss are always positive. The Dirichlet experiment also finds that at around 72 h, the reconfiguration term and outer boundary terms dominate over the loss term for the $L2$ norm, $E(1),E(3)>E(loss)$⁠. As a result, the high-B Dirichlet experiment reaches a state where the constant churn of material being brought into the domain and diffused inwards dominates over the loss from pitch angle scattering. The Neumann experiment obviously does not experience this. Because the Dirichlet simulation has quickly reached a point where the dynamics are no longer changing (due to the large fixed value of $fOB$⁠), the Neumann and Dirichlet experiments diverge in $E$ even though in general, the outer boundary condition has less effect than loss.

Enhancement location $μ$: As without loss, a higher $μ$ means that $tm$ is reached more quickly. [We also note the minor effect from Fig. S8(f) that for the first few hours, the loss experiments actually find that $Et$ is dominated by reconfiguration $Et(3)$—this will be because the enhancement is at higher L so more diffusion is possible. Then, loss becomes dominant—for Neumann with loss, the reconfiguration and loss contributions to $Et$ become comparable after around 70 h. “Reconfiguration” $Et(3)$ is not as strong as without loss.]

Amplitude A: There is no change with overall amplitude parameter A—the same as without loss. We see the same general patterns in $N,E$ as discussed in Sec. V A 3.

Enhancement width $σ$: Again, a narrower width (i.e., lower $σ$⁠) reaches $tm$ quicker, the same as without loss. From $N,E$ analysis [shown in supplementary material, Figs. S7(g), S7(h), S8(g), and S8(h)] we find the same overall results as discussed in Sec. V B 1. With loss, the Neumann and Dirichlet runs are more similar to each other than to the runs without loss. With loss, less is lost through the inner boundary.

Without loss, we found that $tm$ varied when we changed the location of a Neumann outer boundary condition, but not a Dirichlet condition. With loss, we find that $tm$ dependence on $Louter$ is very small for both a Neumann boundary [Fig. 17(a)] and a Dirichlet boundary, particularly when the domain boundary and the plasmapause are close [Fig. 17(b)][P1c].

When including loss, the extent of the lossy region is a new parameter to consider. The outer limit of this region is the plasmapause, $Lp$⁠. By default our simulations have $Lp=5$⁠, which is higher in L than the default peak location, $μ=4$⁠. Figures 18(a) and 18(b) show time to monotonicity across a variety of plasmapause locations. The effect for a Dirichlet outer boundary condition is small, with an effect for the lowest plasmapause values (for which $tm$ is sooner), which quickly drops off to reach a value of $tm$ that no longer changes with $Lp$⁠. The effect of plasmapause location with Neumann outer boundary condition is more complex; overall, a plasmapause closer to the Earth reaches monotonicity sooner. However, there is a cutoff in L after which monotonicity is not reached at all; our experiments place this around $Lp=6$⁠. This is unlikely to be a “hard” cutoff but instead due to interactions with $μ$ and $Louter$⁠. These relationships are explored in the following paragraphs.

In Sec. V B 2, we noted that $Lp>μ$ could result in the PSD being unable to reach a monotonic state. This relationship is explored by using three plasmapause locations in the $N,E$ analysis: $Lp=3.5,5,6.$ These results are shown in the final three rows of Fig. 18.

As the default enhancement location is $μ=4$⁠, a plasmapause at $L=3.5$ is at a lower L than the enhancement. Figure 18(e) indicates that material lost from pitch angle scattering is quickly similar to loss from the outer boundary, and Fig. 18(h) demonstrates that $Et$ is dominated by the reconfiguration term $E(3)$ (although this becomes comparable to $E(loss)$ for the Neumann experiment by the end of the week). With less overall loss, the Neumann and Dirichlet experiments still have distinct values of $N,E$⁠.

For a plasmapause at $L=5$ or 6, much more material is lost, as expected since the proportion of particles lost increases with L. The norm for $Lp=5,6$ is quickly very low [Fig. 18(d)]. For $Lp=6$⁠, the Neumann experiment reaches a lower energy state, while for $Lp=5$ the Dirichlet experiment has a lower norm. This is due to more material coming in the outer boundary with a higher $Lp$⁠.

The relative location of peak and plasmapause is the determiners of whether a local minimum is at all possible, while details of $μ$⁠, $σ$⁠, and $Lp$ combine to see whether it occurs in a given simulation. For example, despite the entire enhancement lying inside the plasmapause for $Lp=6$⁠, we see that the local minimum in Fig. 19(e) is very small, even though a local minimum very clearly occurs when the peak is just inside the plasmapause (e.g., $μ=4,Lp=5$⁠, Fig. 19 row 2). In this case, the loss from pitch angle scattering results in the peak rapidly moving inwards, and a local PSD minimum arising between $L=4$ and $L=5$ [Figs. 19(c) and 19(d)]. Peak width may change the rate at which a local minimum arises (e.g., a wider enhancement could mean a higher PSD to lose before reaching a local minimum, while a narrower enhancement could mean that diffusion occurs more quickly, offsetting the pitch angle loss), but it is the relative location of $μ$ and $Lp$ that determine whether this is possible.

For a higher plasmapause location, $E(loss)$ becomes increasingly larger. When $E(loss)>E(3)$⁠, diffusion cannot prevent the formation of the extra minimum demonstrated in Sec. V B 2. Indeed, this is the relationship observed when plotting several intermediate time instances of the experiments with $Lp=3.5,5,6$⁠, shown in Fig. 19. For $Lp<μ$⁠, this minimum does not form (first row of Fig. 19) while the minimum is larger as $Lp$ increases (Fig. 19 second and third rows).

Note that since $E(3)t$ depends on the gradient of the PSD, the growth of this minima (i.e., when $Et(loss)>Et(3)$⁠) is dependent on the initial conditions. Loss may or may not dominate the evolution of the system over reconfiguration. Overall, our results emphasize that the plasmapause location relative to the enhancement is important and that a more distant plasmapause has so much more loss that this can totally change the dynamics [P1b; S1b; S2; S4].

Monotonicity was easier to obtain when there was already a significant background PSD, corresponding to a large high-L source (i.e., high B). Therefore, $tm$ corresponds to the timescale of radially diffusing a local enhancement, with respect to the background PSD (the background particle population). If the enhancement location $μ$ occurs at high L, it does not last as long before the distribution becomes monotonic, although including loss from pitch angle scattering means that the relative location of the peak and plasmapause strongly interact to affect $tm$⁠. Furthermore, a narrower enhancement will be reduced more quickly than a wider one and using $Et$⁠, we attributed this to the gradients of the PSD, which will be discussed in Sec. VI E. Less intuitively, we found that the time for our enhancement to “fade” into the background varies significantly with the outer boundary location if we used a Neumann boundary. This is a key result, and we explore the consequences of this and future avenues in Sec. VI D. Finally, we found that loss from pitch angle scattering has a strong effect on the final shape of the distribution (and therefore the timescale for radial diffusion). Indeed, some numerical experiments never reach monotonicity, casting some questions about future appropriateness of $tm$⁠, which we discuss in Sec. VI A [P1b;S1b;S2;S3a;A2].

Using the evolution of $N$⁠, we could work out what processes were going on, and from $E$ we could work out how the system was evolving to reach a maximally diffused state. We found that although theoretically an experiment with a Neumann OBC can reach a state with a lower $L2$ norm (i.e., zero everywhere), the majority of the time Dirichlet experiments were reaching a state of lower $E$⁠, because mass could be lost from both boundaries. The Neumann experiments were reaching a state where $Et∼0$⁠; where the dynamics were changing very little once the plateau had risen up, while the Dirichlet experiments were still diffusing material from the enhancement by the end of the week [S1b;A1].

Other specific results from the use of $N,E$ include the importance of gradients in the distribution; using $Et(3)$ we can estimate the comparative effect of the $L−$ dependence of $DLL$ and the gradients in the distribution on the evolution of f. We found that gradients are more significant, a result worthy of its own discussion section Sec. VI E. We can also see in Fig. 12(c) that although using a different domain (i.e., a different $Louter$⁠) does not significantly impact the starting value of $E$⁠, it does affect the rate at which the system moves toward a steady state (more in Secs. VI D and VI B). Finally, we can quantify that loss from pitch-angle scattering has a stronger effect on the evolution of f than radial diffusion does [S4;A2].

Understanding of the research goals, and the context of our results, has developed throughout the research lifetime of this work. The initial questions shaped the methodology and investigation of the results shaped the narrative of the analysis. Here, we put the results into context with existing literature and with future modeling choices. To aid navigation, relevant paragraphs are labeled with our initial research questions in Sec. III. Each research goal may be addressed in multiple paragraphs.

Error metrics such as the log-accuracy ratio⁵⁷ are often the first tools considered for comparing distributions; this would be a suitable method to compare the deviation between two phase space density (PSD) distributions f, such as between observations and models (although weighting by the Jacobian $1/L2$ may be necessary, as it has been here). When no “truth” is available to compare against, error metrics become an unsuitable tool as it would require a threshold (e.g., when two distributions are “close enough,” or when radial diffusion is “done”), which would be difficult to motivate objectively. In this section, we discuss the performance of tools $tm$ (time for PSD reach a monotonic state), $N$ [mass density, Eq. (11)], and $E$ [Eq. (12), energy density, or “distance from zero state”] when comparing distributions and how the distribution morphology became so integral to our analysis.

$tm$ and $N,E$ are found to be complementary measures of the shape and evolution of the distribution. They are related, as the state of monotonicity and the state of lowest possible $E$ both have PSDs predominantly weighted toward the outer edge of the domain. $tm$ is a state where the only existing gradients are to the left hand side of the enhancement. Indeed, $E$ particularly penalizes high f at low L [e.g., $f=1$ at $L=4$ results in a higher $E$ (1/16) than at $L=5$ (1/25)]. (Note, however, that a monotonic distribution can be gaining particles from the outer boundary and, therefore, be increasing in $E$⁠). Despite the similarity between low- $E$ and monotonic states, evolution is determined by $Et$⁠; the steady state reached may not be the lowest $E$ configuration possible. $Et$ is in turn usually dominated by reconfiguration term $Et(3)$⁠, which is an integral across L that is also dependent on L. This term tells us that gradients in f (and their location in L) are the strongest factors dominating $Et(3)$⁠. As a result, the distribution will change more rapidly toward the low- $E$⁠, monotonic-like state when there are steeper gradients, and when those gradients are situated at low L.

The dominance of $E(3)$ consequently suggests that quantities capturing the shape of the distribution and the location of features in it are going to be the most useful tools when analyzing the output - as we have found. We conclude that measures using the distribution morphology across the entire domain are necessary to capture the system evolution. Indeed, the outer boundary condition affects the monotonicity of the final state. The shape of the distribution and the ongoing evolution change with the location and condition of the outer boundary. The act of “reducing gradients” means that in general, Neumann experiments reach monotonicity quickly, while the steep gradients remaining in the Dirichlet experiments means that they are actually at lower energy states and still changing more rapidly (losing more $L2$ norm and mass) toward a steady state.

We suggest that we finally settled on this pair of tools because they are (a) both domain dependent and (b) include information about the morphology of the PSD distribution. Both tell us about the entire system, i.e., $E,N$ tells us about the system across all L at each time step; $Nt,Et$ inform us about the ongoing evolution of the system at each time step, and $tm$ is a measure of the long-term evolution. Initially, we searched for domain in-dependent measures; but as the diffusion coefficient, $DLL$ is domain dependent, so should be our method of analysis. Both tools relate, predominantly, to gradients. The existence of steep gradients is the largest contributor to a rapidly changing $E$⁠; a distribution rapidly heading toward a steady state. A moderate limitation of $N,E$ is that one can only use idealized diffusion coefficients rather than the empirical ones from Ref. 39 used in calculating $tm$⁠. Unfortunately, a limitation of $tm$ is that some parameter combinations never reach a monotonic state, i.e., once loss is included at higher L values (see Sec. V B 2 and Fig. 19). The extra minimum in these cases suggest that we may need to characterize by convexity (i.e., number and location of extrema) rather than simply monotonicity and $tm$⁠. Nevertheless, we found $tm$ to be an extremely useful measure to compare the different experiments, which could be related back to radiation belt processes, to a “quiet” state and to quantifiable properties of the PSD f. Using both $tm$ and $N,E$ together gives us a nuanced description of timescale (Sec. VI B). These quantities allowed us a clear and efficient means of comparing experiments and understanding the processes behind the evolution of each [S3a;S3b].

A radial diffusion “timescale” is poorly defined when $DLL$ depends on L; one cannot simply take $1DLL$ as the characteristic timescale. This is an open question both for idealized models and the real radiation belts, because of the number of underlying approximations. Potential measures of timescale include the autocorrelation or dimensional analysis of self-similar solutions (e.g., Ref. 32) Unfortunately, self-similar timescales are difficult to find for the 1D radial diffusion equation Eq. (1) with diffusion coefficients Eq. (2) or Eq. (5) because of the high order of L. Typically, one makes “by eye” judgments of storm time radiation belt simulations. One could numerically solve for the equilibrium solution and then find the time taken to reach this (e.g., Ref. 58), but this will also depend on the metric or threshold one chooses to define as close enough to this equilibrium solution. Alternatively, suggestions were to take either the minimum or maximum $1/DLL$ in the system, as this would at least give you a “fastest possible” timescale (various personal communications). We have explored methods of quantifying timescale using our tools.

Using $tm$⁠, the concept of timescale reduces to “how long until an enhancement is diffused away.” The results in Sec. V A, therefore, include how the initial condition (e.g., size and location of enhancement) relates to timescale without loss. There is a large variability with Kp $,B,μ$⁠, etc. From $tm$ at a Kp of 4–5, we conclude that the timescale is on the order of days, or tens of hours (for all the $tm$ plots together, see Fig. S1). Some combinations of initial conditions will keep $tm$ close to one day, others closer to a week. We show some example timescales using each of these measures in Table I. Based on these results, we conclude that time to monotonicity is a more representative measure of timescale, but has practical limitations based on the outer boundary and on the fact that not all distributions can reach a monotonic distribution. This definition is not particularly suited once one considers loss from pitch angle scattering; although the timescale drops to hours or days, monotonicity is no longer guaranteed it is a poor indicator of timescale [P2a,b,c].