Optimization of Radial Diffusion Coefficients for the Proton Radiation Belt During the CRRES Era

Proton flux measurements from the Proton Telescope instrument aboard the CRRES satellite are revisited, and used to drive a radial diffusion model of the inner proton belt at 1.1 ≤ L ≤ 1.65. Our model utilizes a physics‐based evaluation of the cosmic ray albedo neutron decay (CRAND) source, and coulomb collisional loss is driven by a drift averaged density model combining results from the International Reference Ionosphere, NRLMSIS‐00 atmosphere and Radio Plasma Imager plasmasphere models, parameterized by solar activity and season. We drive our model using time‐averaged data at L = 1.65 to calculate steady state profiles of equatorial phase space density, and optimize our choice of radial diffusion coefficients based on four defining parameters to minimize the difference between model and data. This is first performed for a quiet period when the belt can be assumed to represent steady state. Additionally, we investigate fitting steady state solutions to time averages taken during active periods where the data exhibits limited deviation from steady state, demonstrated by CRRES measurements following the March 24, 1991 storm. We also discuss a way to make the optimization process more reliable by excluding periods of variability in plasmaspheric density from any time average. Lastly, we compare our resultant diffusion coefficients to those derived via a similar process in previous work, and diffusion coefficients derived for electrons from ground and in situ observations. We find that higher diffusion coefficients are derived compared with previous work, and suggest more work is required to derive proton diffusion coefficients for different geomagnetic activity levels.

In the outer zone, larger and more rapid changes in the geomagnetic field lead to frequent rearrangement, requiring extra modeling to capture. Such changes are not considered here, but include: the direct injection of high energy (≳10 MeV) solar energetic particles (SEPs) through the front-side magnetopause (Kress et al., 2004(Kress et al., , 2005, and loss induced by adiabatic expansion of drift orbits combined with field line curvature scattering (Anderson et al., 1997;Engel et al., 2015Engel et al., , 2016. In addition, ≲1 MeV SEPs provide a low energy source, penetrating the magnetosphere via open field lines in the magnetotail (Blake et al., 2019). Injection of such low energy protons into the outermost trapped population near geostationary orbit has been shown to occur during impulsive reconfiguration of the geomagnetic field (i.e., Baker & Belian, 1985;Wang et al., 2008). It is thought that enhanced radial diffusion is then required during active periods for inward transport to the proton belt (Boscher et al., 1998). Time-dependent models have embraced the challenge of including such dynamical processes, but uncertainty in empirical inputs, including parameters for radial diffusion, is a persistent issue (Selesnick & Albert, 2019).
In this work, we follow a similar method to Albert et al. (1998) and present modeling of equatorial flux within the CRRES era, revisiting measurements taken by CRRES's proton telescope (PROTEL) with a radial diffusion model that includes more modern evaluations of key source/loss processes, as well as improved theoretical modeling of coulomb collisional loss. We investigate the process of steady state optimization to derive new estimates of proton radial diffusion coefficients for the inner zone. After an introduction to the data (Section 2), we introduce our model and discuss its empirical components (Section 3). We then show a method for selecting different average periods of data in order to eliminate variability (Section 4), justifying the use of time-averaged inputs to represent time-varying processes. Modeling results are then presented and discussed, along with our optimized diffusion coefficients which we compare to other works (Section 5). We conclude with a summary of our findings and recommendations for future attempts at steady state modeling (Section 6).

Data Overview
The Proton Telescope (PROTEL) instrument on-board the CRRES satellite made measurements of differential flux on 24 energy channels in the 1-100 MeV range with full pitch angle resolution (see Violet et al., 1993). Measurements were made from elliptical orbit (350 km perigee, 36,000 km apogee) at 18° inclination. Data for the present study is extracted from the 'pad' files made available by the NASA Space Physics Data Facility (https://spdf.gsfc.nasa.gov/pub/data/crres/), which contain flux averaged over 1 minute intervals as a function of equatorial pitch angle, from August 15, 1990until October 11, 1991. These measurements were mapped from local to equatorial pitch angle using the IGRF85 internal and Olson LOZINSKI ET AL. and Pfitzer (1974) external magnetic field models. After excluding orbits which exhibit bad data (Brautigam, 2001), this data set spans 979 orbits.
CRRES observed dynamic changes in trapped flux due to numerous magnetic disturbances, the most significant being the March 24, 1991 storm occurring roughly half way through the mission. In this event, a large interplanetary shock compressed the magnetosphere causing a storm sudden commencement (SSC). This was accompanied by the arrival of an SEP event, in which solar protons were able to penetrate the magnetosphere and become trapped in the outer zone down to L ∼ 2. This trapping was enabled by the combination of two main factors: first, the fast suppression, and subsequent restoration, of geomagnetic cutoff limits over the course of the SSC. These limits define access regions for incoming particles due to attenuation by the geomagnetic field, and their temporary suppression allowed particles of the same magnetic rigidity closer toward Earth. Second, in tandem, the shock's compression of the magnetosphere induced an azimuthal electric field pulse which led to inward acceleration and transport of trapped particles drifting in time with the pulse (Hudson et al., 1997;Li et al., 1993). This led to a large enhancement in trapped flux over the timescale of a drift orbit which persisted until at least the end of the CRRES mission. Both rapid and diffusive inward transport (to lower L) lead to betatron acceleration of trapped particles due to violation of the third adiabatic invariant and simultaneous conservations of the first and second. This betatron acceleration increases a particle's momentum perpendicular to the magnetic field, and can therefore increase pitch angle toward 90° as the ratio p ⊥ to p ∥ increases. One signature of the CRRES enhancement is therefore highly anisotropic pitch angle distributions at certain affected L. The event followed a period of magnetically quiet conditions, and therefore separates the CRRES data into a quiet and active era. Other magnetic disturbances during the CRRES mission are summarized in Table 1 of Hudson et al. (1997). PROTEL flux measurements averaged over two ∼200 day periods encompassing the quiet and active era have been used to build the CRRESPRO Quiet and Active static radiation belt models (Meffert & Gussenhoven, 1994). The method used to prepare flux maps for the CRRESPRO models is described by Gussenhoven et al. (1993), and data for the current study has been processed using a similar method but for different time average periods. These steps are described in Section 2.2.
An overview of PROTEL data is presented in Figure 1 as a time series of weekly averaged differential flux at 90° equatorial pitch angle (left panels), along with a measure of anisotropy (n) of each pitch angle distribution (right panels), at L bins from 1.3 to 2.3 (ΔL = 0.05). Both quantities have been calculated by fitting each weekly averaged pitch angle distribution using a sin n fitting function, described in more detail in Section 2.2.2. The narrowness of each pitch angle distribution center peak is therefore controlled by n, with a higher n indicating a more anisotropic (narrower) distribution. The time axis spans the whole period for which data is available. The quiet and active periods used for the CRRESPRO quiet and active models are indicated in Figure 1 by two vertical dashed white lines, which show the end of the quiet period average and beginning of the active period average. Magnetic disturbances are also indicated in Figure 1 by dotted black vertical lines, marking the arrival of storm sudden commencements as listed in Table 1 of Hudson et al. (1997). Figure 1 demonstrates the large increase in intensity throughout the outer zone (left panels) following the March 1991 storm (and subsequent enhancements), and the effect on particle pitch angle distributions (right), which become more anisotropic during the active period. Some uncertainty, particularly in the n fitting parameter at L ≲ 1.45, arises because one week is not a sufficiently long average period, leading to some fits based on only a few measurements. We have therefore excluded these values from Figure 1 and use much longer flux average periods in later analysis. From the 9.7 MeV channel and below in Figure 1, increases in anisotropy at L ∼ 1.7 are shown during the active period. New particles were injected into the region 2 < L < 3 (Blake et al., 1992), but this change at lower L could be the result of inward transport and betatron acceleration leading to a higher proportion of 90° particles, illustrating the extent of the region affected by enhancements.

Data Processing
The 1 minute averaged pitch angle data files introduced in Section 2.1 make available average differential flux measurements binned by equatorial pitch angle and L. Each average flux measurement in each pitch LOZINSKI ET AL.

10.1029/2020JA028486
3 of 19 angle/L bin corresponds to a number of observations made within a 1 minute interval, and this number is included in the data along with average ephemeris for each interval. The following processing steps are already implemented on the available. pad files: noise spike removal (see Brautigam, 2001); mapping of the data into equatorial pitch angle bins (5° bins from 0 to 90° for 18 bins in total); binning of the data in L (ΔL = 0.05); and a loss cone correction to remove background flux. There is some uncertainty about the effectiveness of the loss cone correction in particular, with implications for data accuracy, and we include some details about this in the next section.
LOZINSKI ET AL.  We follow convention by ignoring the 8.4 MeV channel due to its overlap, and the 15.2 MeV channel due to sensitivity issues. We also ignore the first two channels at 1.5 and 2.1 MeV because these are below the energy range our model can simulate (the CRAND source begins at 2.31 MeV).

Loss Cone Correction
During computer simulations it was found that >60 MeV protons could pass through PROTEL shielding and be erroneously counted, leading to a background contamination . To correct for this, Gussenhoven et al. (1993) describe a loss cone correction whereby the flux at the pitch angle bin just inside the loss cone is treated as the background noise level. Flux at the pitch angle bin of the loss cone and below is set to zero, whilst the background is subtracted from all higher pitch angle bins. Gussenhoven et al. (1993) validated the results of this correction using a Monte Carlo ray tracing code to derive the PROTEL response function. However, the correction has been shown to be inadequate for the 6.8 and 8.5 MeV channels (too large) and for energies >40 MeV (too small). Furthermore, the correction is similar in magnitude to the measurements themselves at L < 1.4.
To deal with these caveats, in addition to the measures listed in Section 2.2, we ignore the 6.8 MeV channel, which appears to cause a particularly strong deviation in the spectrum at low L. We expect that data driving the model at >40 MeV may represent a slight (∼20%) overestimate owing to limitations of the loss cone correction, but as the majority of energy channels influencing our optimization method are below 40 MeV we expect this will not have a large effect on results. We ignore all data below L = 1.35 and draw attention to uncertainty by calculating the standard deviations of PROTEL measurements (in terms of phase space density) and including them in plots of our results. We expect potential uncertainties in the data at L = 1.35 to have a minimal influence on results as our comparison is dominated by measurements at higher L.

Filtering and Fitting
The next step required to process PROTEL. pad files are the filtering out of anomalous records of flux caused by SEP events, to ensure measurements are indicative of trapped flux only. To address this we follow the method of Gussenhoven et al. (1993): for any time average period, 1 minute average flux values are first averaged, weighted by the number of observations; the mean flux and standard deviation in each pitch angle and L bin is then calculated; and finally, for each bin, flux measurements outside two standard deviations are excluded, and the mean/standard deviation is then recomputed.
In addition to occasional periods of data unavailability, there is a shortage of measurements at 85 and 90° equatorial pitch angle during two ∼month long phases of the mission, beginning around February 1991 and September 1991, and in general there are fewer measurements of equatorial flux at L ≲ 1.5 compared to higher L. The process of averaging equatorial flux over time periods, or ranges in L, coinciding with this shortage of measurements, gives averages associated with few observations. The final processing step is taken to mitigate this issue, and involves fitting all time averaged equatorial pitch angle distributions using the function where j is differential, unidirectional flux, α is equatorial pitch angle, and A, n are fitting parameters. Combined with the loss cone correction which has already been applied, this fit has been shown to work well over this range in L (Valot & Engelmann, 1973). Figure 2 shows the application of this fit to equatorial pitch angle distributions at L = 1.7. Time averaged flux distributions are shown at 2.9, 5.7, 10.7, and 30.9 MeV (rows 1 to 4), for three average periods during the CRRES mission (columns 1 to 3), taken before the March 1991 storm (Quiet), and at two intervals after (Active 1 and Active 2). We use these average periods for later analysis and explain their selection in Section 4.
Average flux values in each pitch angle bin (black crosses) are plotted along with the corresponding standard deviation (vertical blue lines). The best fit (red curve) and fitting parameters are also shown for each distribution. Figure 2 shows the advantage of using this fit to get a stable measurement of 90° flux, especially for the Active 2 period where standard deviation is higher. The fit is weighted by flux at lower pitch angle, LOZINSKI ET AL.
10.1029/2020JA028486 5 of 19 for which many more observations exist given that dwell time off-equator is comparatively high. One disadvantage of this fit is that for highly anisotropic distributions (n > 10), the fit can sometimes seem to overpredict 90° flux. Figures 1 and 2 show that during the geomagnetically active period following the March 1991 storm, time averaged pitch angle distributions generally become more anisotropic, with the n fitting parameter increasing through time. Peak flux indicated by the A fitting parameter also undergoes changes at a given L depending on energy channel, with some channels undergoing a flux increase at L = 1.7 (i.e., 2.9 and 5.7 MeV in rows 1 and 2 of Figure 2), and some showing a decrease in flux (i.e., 30.9 MeV in row 4 of Figure 2).

Model Overview
Our model describes relativistically the distribution of equatorially mirroring particles in terms of phase space density. Phase space density is given by j/p 2 with SI units of m −6 s 3 kg −3 , where j is unidirectional differential proton flux and p is momentum. However, we model the quantity f (μ, L) given by  3 2 0 / f m j p , where m 0 is the proton rest mass. This quantity is proportional to phase space density by a constant and is used so that we may later plot our results alongside the results of Albert et al. (1998), who work in terms of a phase space density with units of km −6 s 3 . We solve for f using a proton radial diffusion equation derived by first expanding the Fokker-Planck equation in terms of three adiabatic invariants μ, J, and L and considering only the D LL element of the diffusion tensor (see Beutier et al., 1995, etc.): This equation is then expanded using the result (Farley & Walt, 1971), then J is set to zero to consider equatorial particles only. As extra source and loss mechanisms, we consider CRAND (S n ) and loss from nuclear collisions and charge exchange (Λf) to arrive at our full model equation: The quantity  fric d dt represents the change in μ due to the cumulative effect of many small deflections by free and bound electrons in the atmosphere, ionosphere and plasmasphere (coulomb collisions). It is given by: following the definition of μ, where dT dx is the mean kinetic energy loss averaged along a particle drift orbit, B e is the magnetic field strength at the equator along the particle drift and α eq is equatorial pitch angle.
The quantity dT dx has contributions due to drift averaged electron density (n e ) and gas molecule density (n i ), and is given in SI units by: LOZINSKI ET AL.
10.1029/2020JA028486 6 of 19 where e is electron charge, Z i is the number of bound electrons for each constituent, I i is the mean excitation energy for bound electrons and λ D is the drift average Debye length (Schulz & Lanzerotti, 1974).
Three boundary conditions are applied to the model: where L min = 1.1, L max = 1.65, and f b is the outer boundary spectrum derived from PROTEL data at L max . Inside the energy range of the included PROTEL channels, f b is linearly interpolated. Outside this range, f b is extrapolated by fitting a second order polynomial function p to the outer boundary spectrum such that log j = p(E), where E is energy. The gradient dp/dE at the lowest (highest) energy channel is then used to linearly extrapolate log j to lower (higher) energies, giving a continuous spectrum. The magnetic field is given by a magnetic dipole model with dipole moment 7.83 Am 2 , representative of the era. Whereas Albert et al. (1998) set their model's data-driven outer boundary at L = 1.7, we set ours at L = 1.65 based on an analysis of the variation in fluxes (see Section 4.2) in order to be more selective of the region steady state modeling can be applied to.

Drift Averaged Density
Coulomb collisions with free and bound electrons in the plasmasphere, ionosphere, and atmosphere cause energy loss, which in modeling terms results in the convective flow off to lower μ. The loss rate is controlled by the number density of various constituents in addition to the local spectrum ∂f/∂μ. To calculate Equation 3, we compose a model of drift averaged density throughout the region of interest using a combination of three density models: NRLMSIS-00 (Picone et al., 2002), the International Reference Ionosphere (IRI, Bilitza et al., 2017) and the plasmasphere model of Ozhogin et al. (2012). Using our drift averaged density model, density of 14 constituents can be queried as a function of F10.7a (81 day average solar radio flux) and day of the year for a particular model coordinate μ and L, preserving seasonal and solar cycle effects included in the underlying models. This variability is shown for electron density in Section 3.6.
The first step in calculating a drift averaged density model was to solve for a range of proton drift orbit trajectories covering the modeling domain. Protons were initialized with α eq = 90° at μ and L over the region of interest and their drift trajectories calculated by solving the relativistic equation of motion for a magnetic dipole field. Over each particle's orbit, neutral density of He, O, N 2 , O 2 , Ar, H, and N and anomalous O is given by NRLMSIS-00 below 2,000 km and extrapolated above this height as a straight line of   log10 i n with altitude. Ion density of O + , H + , He + , O 2+ , and NO+ is given by IRI below 2,000 km and is extrapolated the same way. Electron density is also given by IRI below 2,000 km, and above 2,000 km is interpolated between IRI and the Ozhogin et al. (2012) plasmasphere model value at L = 3.25. The latter is based on a fit to Radio Plasma Imager data and is a static value, unlike the IRI and NRLMSIS-00 models which include various dependences. However, historical measurements indicate relatively consistent electron density at L = 3.25 below the plasmapause (see for example Figure 10, Park et al., 1978), and there is a good agreement between empirical models here too (Ozhogin et al., 2012).
To achieve the solar cycle and seasonal dependence of our drift average density model, five different values of F10.7a were considered (60, 100, 140, 180, and 220 sfu) along with four days of the year (DOY; day 70, 170, 270, and 330). 20 past environments matching each possible combination of F10.7a, DOY were found by searching through a time history of indices. Once a date had been found with a F10.7a, DOY coordinate matching each combination (within a margin of ±10 sfu, ±4 days), density of each constituent was drift averaged by querying each model in the fashion described above using the time corresponding to that F10.7a, DOY coordinate as model input. To eliminate day-night variation in density, each drift average was performed six terms (at 4 hour intervals) throughout the day selected. MLT variation is automatically accounted for by considering the particle trajectory around Earth. When querying density at a test point using the NRLMSIS-00 and IRI models, we occasionally observed large time-dependent changes appearing to LOZINSKI ET AL.
10.1029/2020JA028486 7 of 19 correlate with jumps in daily AP index. This variability was excluded from the above drift averages by manually checking that each of the 20 environments corresponded to magnetically quiet conditions. This drift average density model, combined with our model equation, represents a step forward from the treatment of coulomb collisions by Albert et al. (1998), which was limited to the effect of plasmaspheric electrons.

Diffusion Coefficients
Radial diffusion coefficients for equatorial protons are parameterized according to from Claflin and White (1974), where A m and A e are constants proportional to the power spectrum of certain features in the magnetic and electric field perturbations. This formulation allows for the L and μ dependence giving the best fit to data to be determined via optimization of the free parameters α, β, A m , and A e .

Nuclear Collisions and Charge Exchange
Nuclear collisions and charge exchange are taken into account via the loss term Λ in Equation 3 given by Λ = n H σv, where n H is the neutral Hydrogen density from the drift average density model, σ is the proton charge exchange cross section in atomic Hydrogen (taken from Spjeldvik, 1977), and v is proton velocity. Albert et al. (1998) show this term has an inverse timescale several orders of magnitude below radial diffusion and coulomb collisions, and we also find its inclusion to have minimal effect on our results.

CRAND
Following the definition  3 2 0 / f m j p , the rate of change in model phase space density due to the CRAND process is given by considering only the change in phase space density from additional flux of neutrons, where ∂j p /∂t is the rate of change in proton flux at a set of adiabatic coordinates caused by the pickup of newly produced protons. As the set of adiabatic coordinates corresponds to a drift orbit, ∂j p /∂t is given (following Dragt et al., 1966) by where the integral takes place along a proton drift orbit S, ds is an element of length along the trajectory, j n is the flux of neutrons arriving from the atmosphere and decaying into protons along the drift orbit, and γτ n is the relativistic neutron lifetime where τ n = 887 s. At each point along S, j n is determined at the point on top of Earth's atmosphere that is intersected by following the negative tangent to S in the ds direction.
As neutron flux at the top of the atmosphere is nonisotropic and specified with respect to zenith angle, the vector d s is then transformed to local zenith to find the neutron flux pointing toward each drift orbit position. For points along the proton drift orbit where there is no atmospheric intersection, j n = 0. The quantity   n j therefore represents the average atmospheric neutron albedo flux in the direction of a proton along its drift. Equation 8 assumes that the proton produced moves in the same direction, and absorbs all the kinetic energy, of the decaying neutron (Selesnick et al., 2007).
To deal with this difficult computation in modeling, Dragt et al. (1966) introduced an injection coefficient χ(L, E) that allows the approximation of   n j based on some simplifying assumptions. This method is demonstrated by Claflin and White (1974) to calculate an equatorial CRAND source, and relies on measured zenith angle distributions of neutron flux at a point on top of the atmosphere (for example, White et al., 1972), combined with knowledge of the latitudinal variation of total escaping neutron flux (Jenkins LOZINSKI ET AL. 10.1029/2020JA028486 8 of 19 et al., 1970;Lingenfelter, 1963). Albert et al. (1998) use the CRAND source of Claflin and White (1974), but with help from the authors we identified an error causing S n to be overestimated by a factor of pi. Selesnick et al. (2007) demonstrate a more thorough approach, whereby directional neutron albedo flux at the top of the atmosphere has been computed by modeling the interactions of incoming galactic cosmic rays throughout the atmosphere. These calculations give j j E z R W n n n cv    , , , , where E n is kinetic energy, z is zenith angle, cv R is vertical geomagnetic cutoff rigidity, and W is sunspot number. Using this grid of albedo neutron flux, the integral in Equation 8 can be performed numerically via the following steps: (a) initialize a proton at a given value of μ, α eq , and L inside the modeling domain; (b) trace its drift orbit; (c) go back along the drift orbit, calculating the first intersection between the negative tangent and surface representing the top of the atmosphere; (d) if such a point exists, transform its coordinates to be in terms of z, cv R ; (e) interpolate j n from the albedo neutron flux grid; (f) repeat for each step along the proton drift orbit to calculate the drift average   n j at the desired W, then repeat for other drift orbits across the model domain.
Using the neutron albedo data described above,   n j has been calculated according to this method. The simplified model of vertical geomagnetic cutoff rigidity by Smart and Shea (2005) where λ is dipole magnetic latitude, r is the distance from the dipole center in cm, and cv R has units GV. The averaging was performed for five values of F10.7 (60, 100, 140, 180, and 220 sfu), allowing the solar cycle dependence of S n to be taken into account. The top of the atmosphere was approximated using the WGS84 ellipsoid (Decker, 1986) with an extra 100 km radius, tilted with respect to the dipole axis. By drift averaging down to lower energies, we extend the CRAND source of Selesnick et al. (2007) from ∼20 MeV down to 2.31 MeV (the minimum energy for which neutron albedo data is available), with an upper limit of 500 MeV, covering most of the PROTEL energy range.

The Influence of Plasmaspheric Density
The variation in plasmaspheric electron density captured by our drift averaged density computation, along with the impact this variation has on steady state phase space density, is demonstrated in Figure 3. In the top left panel, the density profiles on four days of the year are shown for a fixed F10.7a value over the model L range. In the bottom left panel, density at five values of F10.7a are shown for a fixed day of the year. These plots show the seasonal and solar cycle dependence of electron density respectively. Right-hand panels show the corresponding steady state solutions, calculated with the model using a fixed CRAND source and using the total D LL values for protons of Selesnick and Albert (2019) as an example.
Each profile of electron density in Figure 3 (left panels) is fixed at L = 3.25 to the Ozhogin et al. (2012) model value, but variation across the L range plotted arises from the ionospheric density given by IRI. The dominant component of variation is from solar cycle effects and reflects the IRI dependence on F10.7a as a model input. To produce the steady state phase space density profiles in Figure 3 (right panels), the outer boundary is held constant using PROTEL data averaged over the prestorm quiet era. In reality the outer boundary flux may vary, but the solutions show some differences and highlight the importance of accurately modeling density both for static and dynamic modeling, given the rather short timescale for seasonal changes. The figure suggests the difference in phase space density as a result of density variation is relatively uniform in L but becomes particularly pronounced near L min . In fact, the corresponding change in phase space density at any L also depends on the diffusion coefficients used, as higher D LL tends to limit the region at which coulomb collisions dominate to lower L.

Selecting Average Periods
Using the model described in Section 3.1, we solve Equation 3 for steady state f, assuming ∂/∂t = 0. We aim to find diffusion coefficients that minimize the difference between f and some time-averaged data via LOZINSKI ET AL.
10.1029/2020JA028486 9 of 19 optimization of the four parameters α, β, A m , and A e in Equation 6. To evaluate this optimization, we minimize the mean square deviation in f between the model and PROTEL data over all utilized energy channels from L min to L max , given by There are several assumptions associated with deriving period-correct diffusion coefficients this way, including that (i) all sources and losses are accounted for via accurate modeling; (ii) the time averaged data which we fit the model to represents true steady state; and (iii) the constraints on D LL imposed via its defining parameters in Equation 6 allow for the correct solution.
We must pay particular attention to the first assumption when calculating steady state for a finite average period. A steady state optimization must approximate source and loss processes controlling flux with static inputs over that period, when in fact these processes may be time-varying. For example, even if our calculation of collisional loss is accurate, seasonal changes in density occur. This means averaging density over a longer time period may be averaging over variation, and therefore such an average is less representative in general. Therefore, for best practice, we should ensure that any time averaged inputs used to calculate empirical source or loss terms are truly representative across the whole average period. The ∼200 days long CRRESPRO quiet period may seem like an appropriate time average period over which to optimize model fit, based on the minimal intensity variation shown in Figure 1 before the 24th March storm, suggesting steady state. However, it may be that the proton belt is in a dynamic equilibrium, whereby changes in source and loss rates do occur, but rebalance to maintain the same level of flux. A dynamic equilibrium could also LOZINSKI ET AL.

10.1029/2020JA028486
10 of 19 involve time-dependent radial diffusion, brought about by frequent changes in geomagnetic activity, whereas we can only derive time-independent average coefficients.
Our selection of an average period in which the data represents steady state is therefore influenced by also wanting to eliminate variability in this manner. Our main concern is plasmaspheric density, given seasonal changes can affect the steady state solution (Section 3.6). The variability in the CRAND source is driven by solar cycle and CRAND is therefore relatively constant over the CRRES mission, and inelastic nuclear scattering (driven by neutral density) is known to have only a minor influence on the distribution at our energy range of interest. To gain insight into proton belt variability over the entire CRRES mission, Figure 4 plots electron density through time (first and second panel), along with the change in weekly average flux (plotted as logarithm of the ratio j to j t0 ) throughout the mission according to several PROTEL energy channels at selected L (third and fourth panels). Figure 4 also marks the epoch of all flux measurements near the model outer boundary L over the course of the CRRES mission (bottom panel) to give an idea of data availability. The first step in our analysis method is to use these plots to select suitable steady state optimization periods. The dynamic variation in density shown for electrons in Figure 4 can be produced for all constituents by converting each epoch to its associated F10.7a and day of year index, then linearly interpolating between the drift averaged density values calculated for environments with surrounding F10.7a, day of year coordinates (see Section 3.2). Electron density has been highlighted as it is the main driver of coulomb collisional loss near the model outer boundary.

Figure 4 (top panel) shows a gradual increase in electron density through time until around March 1991
where it peaks before a gradual decrease. The main changes in flux (third and fourth panels) coincide with geomagnetic disturbances during the active period. The third panel, showing flux at L = 1.7, shows two main enhancements: once following the 24th March storm, and once again coinciding with various SEP/ SSC events that occurred near the beginning of June 1991. These enhancements are interesting because they occur below the region of newly injected SEPs and are limited to <10 MeV particles. In fact, the flux in the 13.2 MeV channel and higher appears to decrease following these events.
We choose the ∼five-month period following start of mission as our quiet average period, limiting it to January 28, 1991 to keep electron density variation within 25%. During the active period following the March 1991 storm, it is likely that the proton belt is not in steady state, making data here less suitable to drive our model or optimize against. However, if the deviation from steady state throughout the inner zone during the active period is limited enough (in terms of the range of affected L and energy channels), the minimization of ϒ may still be dominated by comparison with steady state f data , and similar D LL may be derived. We therefore suggest that by filtering out specific energy channels affected by active period enhancements, we may select additional data average periods following the March 1991 storm and repeat our steady state optimization method to get similar results to our quiet time average.
We therefore select three time average periods in total, indicated by the red, blue and amber colors in the bottom panel of Figure 4. The first of these average periods, expected to produce the best fit using a steady state model, is the Quiet period (red) from 15/08/1990 to 28/01/1991. The second average period, Active 1 (blue), spans from 31/03/1991 to 04/06/1991, and the third period, Active 2 (amber), spans from 13/07/1991 to 15/10/1991. Time averaged density for each of these periods is used in Equation 5, and time averaged equatorial flux is used to fix the model outer boundary. To filter out data that is likely to be out of steady state, we limit the Active 1 and Active 2 optimization of D LL to only include data at μ ≥ 105 MeV/G, because this corresponds to ∼6.8 MeV at L max = 1.65, which is the highest energy channel indicating a major disturbance during the active period (red line in Figure 4, fourth panel). At the next energy channel and above, flux is relatively unaffected at L max and below.
We minimize ϒ for each period to derive three sets of diffusion coefficient parameters α, β, A m , and A e . The minimization is performed using the Nelder-Mead method implemented in the SciPy library (Virtanen et al., 2020), with the same initial guess parameters as used by Albert et al. (1998), taken from Schulz (1991): A m0 = 7e−9, A e0 = 1e−4, β 0 = 2, α 0 = 2. We verified convergence on the same result by modifying the initial guess and repeating a set of runs for the Quiet period.
LOZINSKI ET AL.

Selection of Outer Boundary L
Another interesting feature of Figure 4 is the observed decrease in flux across all energy channels at L = 1.65 during the quiet period (fourth panel). Given the region and lack of enhancements during this time, we expect these fluxes to remain near steady state. This change coincides with the gradual increase in electron density given by our drift average model (second panel), which we expect to be the primary driver of loss at this L. Taken together, this appears to show the type of response demonstrated in Figure 3, whereby steady state flux can change due to seasonal changes in density. A decrease in flux would be expected to accompany an increase in density (increasing loss), and the observed fluxes therefore appear to corroborate the results of our F10.7a and day of year dependent drift averaged density calculation.
The gradual decrease in flux is also visible at L = 1.7 (panel three), except in the lowest two energy channels shown (3.6 and 4.3 MeV), where no systematic decrease is observed. The flux of these two energy channels is much less steady and seems to suggest that proton belt flux is subject to dynamic changes near L = 1.7 even during a quiet period. This highlights the approximation in considering an 'inner' and 'outer' zone, given the region of flux affected by outer zone changes is energy dependent. However, given the systematic decrease in flux at L = 1.65 in comparison, we choose L max = 1.65 as our model outer boundary. Although the changes in 3.6 and 4.3 MeV flux at L = 1.7 appear to be small, changes in outer boundary flux can have a large impact on steady state calculations considering this flux is radially diffused inward to form the profile of the inner belt. We therefore consider this choice of L max as strengthening our steady state assumption, compared to using L max = 1.7 or higher.

Optimization Results
The final converged values of α, β, A m , and A e for the Quiet, Active 1 and Active 2 periods are shown in Table 1, along with the corresponding value of ϒ indicating how closely the steady state solution was able to fit the data. ϒ is higher for the Quiet optimization, but this is partly due to the exclusion of μ ≤ 105 MeV/G data during the Active 1 and Active 2 optimizations, meaning the summation in Equation 10 was performed over N = 105 values of f data during the Quiet period, compared with only N = 60 during the Active 1 and Active 2 periods. The decrease in ϒ from Quiet to Active optimizations is despite a similar fit being achieved, and can be explained by the grouping of excluded data: in the Active 1 and 2 case, energy channels across the low energy range of the model are ignored, and the optimization is therefore able to better fit channels at higher energy, lowering the average mean squared difference. In previous investigations before the exclusion of μ ≤ 105 MeV/G data, the Active 1 and Active 2 optimizations were found to produce higher ϒ for the same N = 105, indicating a closer fit for the Quiet period. This variation in ϒ implies potential caveats when comparing the performance of runs based on the value of a minimized parameter calculated with a different set of f data .
The optimization path from initial guess to converged values is shown for the Quiet period in Figure 5. The process was found to complete fastest when optimizing in terms of   log 10 m A and   log 10 e A , and Figure 5 shows a logarithmic axis for A m and A e to reflect this. During convergence, ϒ reduces quickly at first, and after ∼80 iterations the minimization con-LOZINSKI ET AL. tinues but makes negligible changes to the resultant D LL . A similar phenomenon was observed for both the Active 1 and Active 2 optimization runs.
Using the optimized parameters, Equation 6 gives D LL as a combination of an electromagnetic and electrostatic component. According to our optimized values the electrostatic component does not contribute significantly in any case, being a few order of magnitudes lower than the electromagnetic component, and a similar result is shown by the optimized parameters of Albert et al. (1998). However, our optimization process relies on Equation 6 only as a means to derive the μ and L dependence of total D LL to best fit the data. As such there may be a limit on the physical interpretability of this result, because errors in the data, or evaluation of source/loss terms, may slightly alter the D LL that gives best fit, and may therefore also alter the balance between electromagnetic and electrostatic diffusion indicated by the optimized parameters which appears to be quite sensitive. Figure 6 shows the solution in steady state phase space density (blue) for each period as radial profiles, plotted alongside PROTEL data (crosses). Data included in (excluded from) each optimization is marked in blue (red). Excluded data fulfills one of the following criteria: (i) located at L < 1.35; (ii) has μ ≤ 105 MeV/G during the Active 1 or 2 period; or (iii) has E > 45 MeV (this criterion is explained below). The Quiet, Active 1 and Active 2 results occupy columns 1, 2, and 3 respectively, with different values of μ shown in each row. Figure 6 generally shows a good match between the model result and included PROTEL data for each average period, meaning that a set of diffusion coefficients could be found in each case to explain the data assuming steady state. The Quiet optimization (left column) appears to produce the best fit in general, though the difference between average periods is small. In order to give an idea of data stability, the standard deviation is plotted (black bars) for each data point. This represents the standard deviation in 90° flux from each equatorially mapped pitch angle distribution at fixed energy channels, converted to phase space density, and then interpolated at fixed μ. It is therefore an approximate measure, and in fact the data is somewhat more stable than it indicates because the sin n fitting used by the model to read pitch angle distributions is not heavily influenced by a high level of variation in the 90° values alone. Nevertheless, the standard deviations shown indicate that data from the Active periods contain higher levels of variation, particularly at low μ. This is to be expected because the data corresponds to an average following an enhancement, and data availability is also more scattered (bottom panel of Figure 4). In general, the highest standard deviations for data in each period corresponds to L < 1.4, where the data is known to be less reliable due to the loss cone correction (see Section 2.2.1). Comparison between the Quiet and Active 2 period data also shows the enhanced boundary f caused by the series of SEP/SSC events in 1991 at μ = 50 MeV/G. The Active 2 profile at μ = 50 MeV/G appears out of steady state because of the enhanced boundary, and excluding μ ≤ 105 MeV/G data has stopped the optimization process from trying to fit this feature.
Data at μ = 500 MeV/G (row 6 of Figure 6) forms a plateau in the region 1.3 < L < 1.5, and all three solutions show a deviation from data at μ = 500 MeV/G near this feature. Further investigation indicates this deviation is better described in terms of energy, happening at ∼45 MeV and above, and can thus also be seen at μ = 400 MeV/G at low L (row 5 of Figure 6). The same observation was noted by Albert et al. (1998), and we interpret the difficulty that both steady state models have in reproducing this feature as flux at E ≳45 MeV being out of steady state. Therefore, we have excluded the four energy channels in this range from the optimization process to weight our optimization using steady state data as much as possible. The feature could be the result of an injection prior to the modeled period which has diffused inwards, or alternatively could be the result of incorrect data/processing. As Albert et al. (1998) note, the profile at this μ appears to be stable and quite unaffected by the storm on March 24, 1991. Figure 7 shows a direct comparison between the new steady state profiles derived for the Quiet period versus the solution calculated by Albert et al. (1998). As Albert et al. calculated six solutions using different combinations of CRAND and density models, the solution appearing to give the best fit of these six has been chosen for comparison in Figure 7. This solution relies on the CRAND source of Claflin and White (1974) and the Parameterized Ionospheric Model for electron density (Daniell et al., 1995). Our data processing steps and average window are not the same as used by Albert et al. to prepare the data, and therefore we expect some disagreement between the data itself. Figure 7 therefore shows both the data and modeling LOZINSKI ET AL.

10.1029/2020JA028486
14 of 19 results from both works. There is a particularly noticeable discrepancy between the data at high μ, which may arise from the different definitions of phase space density. Figure 7 indicates that a closer fit to the data is achieved by the current model (blue) compared to the Albert et al. model (red), especially at μ ≤ 200 MeV/G, although our modeling region is slightly more limited having an outer boundary at L = 1.65 compared to L = 1.7. At μ = 500 MeV/G, a slightly closer fit is achieved in general using the current model, but neither model is able to recreate the plateau noted previously.
We next compare the diffusion coefficients derived during the steady state optimizations carried out over all three time average periods to diffusion coefficients derived by two other works. Figure 8 shows the diffusion coefficients given by the optimized values in Table 1 at μ = 100 MeV/G and μ = 500 MeV/G (red, blue and amber lines). We compare these D LL to the overall D LL derived by Ozeke et al. (2014) for electrons at L > 2.5 using measurements of ULF wave power, which depend on activity level indicated by Kp index. As Lejosne et al. (2013) show that the energy dependence of radial diffusion is weak, there should be a close correspondence between proton and electron diffusion coefficients. Another difference is that the coefficients given by Ozeke et al. (2014) are derived using the method described by Fei et al. (2006), whereby diffusion rates are split into a magnetic and overall electric component, and this is shown to result in an underestimation by a factor of ∼2 (Lejosne, 2019). In Figure 8, these coefficients are labeled as 'O14' and are shown for three different LOZINSKI ET AL.

10.1029/2020JA028486
16 of 19 Figure 7. Steady state profiles with optimized diffusion coefficients derived using the current model (blue lines), and the model run by Albert et al. (1998), (red lines) which uses the CRAND source of Claflin and White (1974) and electron density from the Parameterized Ionospheric Model of Daniell et al. (1995). Data is plotted as blue and red crosses, corresponding to each model. Data used by the current model represents a time average over the Quiet period of this work (see Figure 4), whereas data used by Albert et al. (1998) represents the time average from the CRRESPRO Quiet model. Figure 8. Diffusion coefficients derived from each optimization over the three average periods, plotted next to those derived in Ozeke et al. (2014) and the model run from Albert et al. (1998) shown in Figure 7 (labeled O14 and A98 respectively). Top and bottom panels correspond to μ = 100 and 500 MeV/G respectively. values of Kp index (light blue, plotted only at L ≥ 2.5). In addition, the result of Albert et al. (1998) from Figure 7 is plotted, labeled as 'A98' (dashed red line). Figure 8 shows that our newly derived D LL is higher than results by Albert et al. (1998) by a factor of ∼2 to ∼5 depending on μ. In Figure 8 we have plotted results up to L = 3, despite the fact our coefficients were only optimized over the inner zone. Analytical expressions for D LL given by Ozeke et al. (2014) were derived using data at L > 2.5. Therefore, it is interesting to note the region of overlap at 1.65 < L < 2.5, into which both sets of diffusion coefficients can be extrapolated. Two key features of the diffusion coefficients given by Ozeke et al. (2014) are energy independence, compared with the fairly weak energy dependence derived in this work, in addition to the dependence on Kp. For our region, this dependence on Kp may be different because the inner zone is for the most part shielded from the geomagnetic activity at larger L that can lead to faster diffusion rates. However, our chosen method is only able to quantify the effects of radial diffusion over a long time scale, as required to form the steady state distribution of protons over the data average period. If the timescale of radial diffusion varies significantly with activity level, as implied by the results of Ozeke et al. (2014) in Figure 8, our resultant radial diffusion coefficients should be considered averages of a more time-varying process corresponding to solar maximum. This could be another reason that we were unable to achieve a perfect fit against the data, as parts of the distribution may have deviated from steady state during each average period due to enhanced radial diffusion over short timescales, given the higher levels of geomagnetic activity around solar maximum. Therefore, an interesting question is how activity dependence should be taken into consideration when extrapolating diffusion coefficients to higher or lower L than the region they were derived in, and how they can be applied to examine shorter timescales. A comparison of our results to diffusion coefficients derived at low L during solar minimum may be a useful starting point to better understand how long term averages can vary.

Conclusions
Steady state phase space density is solved for using a model that includes a physics-based evaluation of CRAND, and a drift averaged density model to drive coulomb collisions capturing solar cycle and seasonal variation, with the outer boundary set by PROTEL data. We optimized the solution following a similar method to Albert et al. (1998) but using our updated model to recalculate time averaged radial diffusion coefficients that govern the transport of relativistic protons in the proton belt during solar maximum at 1.1 ≤ L ≤ 1.65.
We took steps to improve the modeling process conceptually by selecting time averages that exhibited minimal variability, and by carefully excluding data to strengthen our assumption of steady state. In terms of plasmaspheric density, we find that a suitable time average for steady state optimization should be less than six months to avoid potential seasonal variations. However, we find that the diffusion rates required to explain steady state are sensitive to the density. Density is not well constrained by measurements between the topside ionosphere out to L ∼ 1.7 which leads to some uncertainty. A better way to conduct future studies might be to reduce this uncertainty by somehow calibrating density models with spacecraft data during the time period one is required to model if available. In this case, electron density has been inferred using the CRRES Sweep Frequency Receiver (SFR) instrument, but saturation in regions of high density means measurements are not available below L ∼ 3 (Sheeley et al., 2001).
The similarity between diffusion coefficients derived in this study for three time periods indicates that although the 24th March storm caused a large enhancement, steady state optimization could still be performed during the Active 1 and 2 periods with reasonable results. As only a short interval of a few months or less is required to average data, steady state optimization could be performed at regular intervals throughout a long data set to derive updated diffusion coefficients. In general our radial diffusion coefficients are comparable to but higher than the previous work of Albert et al. (1998) by a factor of 2-3 at μ = 100 MeV/G and a factor of 3-5 at 500 MeV/G, and provide a better fit to PROTEL data.
Due to the long time averages, our method for deriving diffusion coefficients is not able to uncover activity-dependent variability as is suggested by diffusion rates derived from in situ and ground based measurements for electrons at L > 2.5, which vary by an order of magnitude or more for a small change of ∼2 in LOZINSKI ET AL.
geomagnetic activity index Kp. We therefore suggest further investigation is required to fully understand the activity dependence of proton diffusion coefficients in our region of interest.