The Boltzmann equation plus Maxwell’s equations provide a closed set of equations to describe the kinetic behavior of a plasma. Typical solutions to this set of equations involve a perturbation approach where the velocity distribution is in the form of f s ( x → , v → , t ) = f s , 0 ( v → ) + f s , 1 ( x → , v → , t ) , where f s , 0 ( v → ) is a zeroth order, time independent distribution (i.e., a Maxwellian or Kappa distribution), and f s , 1 is a first order plane wave fluctuation. In contrast, quasilinear theory seeks to obtain the time evolution of the zeroth order background distribution as it responds to waves with finite amplitudes (Nicholson, 1983). This is an appropriate solution for wave-particle interaction in the radiation belts as it 1) provides a time evolution of phase space density, and 2) observations show that whistler-mode chorus waves have substantial amplitudes that require nonlinear terms to be retained when solving the Boltzmann equation (Zhang et al., 2019; Malaspina et al., 2021).

The quasilinear solution to the Boltzmann equation for electromagnetic waves and relativistic particle energies is derived in Kennel and Engelmann (1966) and Lyons and Williams (1984). The solution is of the form of a diffusion equation: ∂ f 0 ∂ t = ∇ v → ⋅ ( D ¯ ⋅ ∇ v → f 0 ) . Using velocity coordinates of pitch angle α = atan ( v ⊥ / v ∥ ) and momentum p, and assuming azimuthal symmetry, the resulting diffusion equation is ∂ f 0 ∂ t = 1 sin ⁡ α ∂ ∂ α sin ⁡ α ( D α α ∂ f 0 ∂ α + p D α p ∂ f 0 ∂ p ) + 1 p 2 ∂ ∂ p p 2 ( p D p α ∂ f 0 ∂ α + p 2 D p p ∂ f 0 ∂ p ) . (1) Where the diffusion coefficients are all in units of 1/s, and form the symmetric matrix D ¯ = ( D α α D α p D α p D p p ) . (2)

The calculation of diffusion coefficients is detailed in Glauert and Horne (2005). In this paper we use all 4 components of the diffusion tensor in Eq. 2, however we will only show the equations for D α α for brevity. The pitch angle diffusion coefficient is D α α = ∑ n ∫ d X   X D α α n X , (3) where n is the cyclotron harmonic, and X = tan ⁡ ⁡ ( ϕ ) with ϕ as the wave normal angle. The diffusion coefficient for a given cyclotron harmonic and wave normal angle is D α α n X = q s 2 ω i 2 4 π ( 1 + X 2 ) N ( ω i )   ( n Ω s / ( γ ω i ) − sin ⁡   2 α cos ⁡ α ) ⋅ B 2 ( ω i )   g ( X ) [ | Φ n , k | 2 | v ∥ − ∂ ω ∂ k ∥ | ] k ∥ (4)

There is an implicit summation over all of the wave modes ω i that are obtained from solving the cold plasma dispersion relation (defined in Appendix A). The wave frequency distributions are often assumed to be Gaussian, such as B 2 ( ω ) = A 2 ⁡ exp ( − ( ω − ω m δ ω ) 2 ) , (5) A 2 = B w 2 δ ω 2 π [ erf ( ω m − ω l c δ ω ) + erf ( ω u c − ω m δ ω ) ] − 1 . (6)

The parameters ω m , ω l c , ω u c , and δ ω define the Gaussian distribution, and B w 2 is the peak magnetic field intensity of the chorus waves, which shows that D α α ∝ B w 2 . The parameters N ( ω i ) and Φ n , k are obtained from the dispersion relation (see Glauert and Horne, 2005). Furthermore, the wavenumber is constrained by the Landau/cyclotron resonance condition ω − k ∥ v ∥ = n Ω s γ . (7)

The diffusion coefficient for a particle with a given pitch angle and momentum is calculated by first solving the dispersion relation in Appendix A simultaneously with Eq. 7 to obtain the resonant wave frequency and wavenumber. The unspecified parameters in Eq. 4 are then: the wave amplitude at the resonant frequency, B 2 ( ω i ) ; the wave normal angle distribution g ( X ) ; and the background plasma density n e and magnetic field B 0 which are also needed in solving the dispersion relation. The challenge in calculating event-specific diffusion coefficients is in specifying each of these parameters, which are sparsely measured by satellites.

The basic motion of radiation belt particles is to bounce between mirror points in each hemisphere. Along this trajectory the wave properties and background plasma will change, so the diffusion coefficients are bounce-averaged. For the pitch angle diffusion coefficient, the bounce averaging is done as (Glauert and Horne, 2005) 〈 D α e q α e q 〉 = 1 τ B ∫ 0 τ B D α α   ( ∂ α e q ∂ α ) 2 d t , (8) where α e q is the equatorial pitch angle of the particle, and τ B is the particle bounce period. The diffusion coefficients used in this paper are all bounce-averaged.

The baseline set of diffusion coefficients used in this paper are calculated using the PADIE (Pitch Angle and Energy Diffusion of Ions and Electrons) model (Glauert and Horne, 2005; Horne et al., 2013) with the procedure given in Reidy et al. (2021) and the wave dataset from Meredith et al. (2020). This set of diffusion coefficients is based on statistical aggregation of wave parameters from the DE 1, Double Star TC1, and THEMIS (A, D, E) satellites, and plasma density measurements from the CRRES satellite. The wave and plasma data used in the diffusion coefficients are binned at half integer values of L-shell, 1-hour increments of MLT, and for Kp index values of 0, 1, 2, 3, and 4+. Figure 1 shows the pitch angle, momentum, and mixed diffusion coefficients at different Kp bins for the same MLT and L-shell.

The diffusion coefficients from PADIE will be scaled with MEPED data for wave intensity, B _w , and the plasma-gyrofrequency ratio, f p e / f c e . The main benefit of using these statistical coefficients is that they use chorus frequency distributions, B 2 ( ω ) , that are fit to data (Meredith et al., 2020). The frequency distribution is still proportional to the chorus intensity, B 2 ( ω ) ∝ B w 2 , but the normalization factor is computed numerically for each distribution. With this proportionality and Eqs 3, 4, the diffusion coefficients are related to the chorus intensity as D α α ∝ B w 2 . (9)

This proportionality also holds for D α p and D p p , and is the basis for using MEPED data to scale the diffusion coefficients for event studies.

The method of using MEPED observations of electron fluxes in order to estimate chorus wave intensity was first developed in Li et al. (2013) and Ni et al. (2014). This method compares the observed precipitating/trapped electron flux ratio from MEPED to the ratio obtained by solving a 1-D diffusion equation in Kennel and Petschek (1966). The predicted flux inside the loss cone, corresponding to the 0-degree MEPED telescope, is (Kennel and Petschek, 1966) J 0 = ∫ E 1 E 2 ∫ 0 2 π ∫ 0 β S ( E ) 〈 D α α 〉 L C ⁡ cos ⁡ α e q I 0 ( α e q α L C   z 0 ) z 0 I 1 ( z 0 )   A ⁡ sin ⁡ η   d η   d ψ   d E . (10) where I ₀ and I ₁ are modified Bessel functions. The corresponding flux outside the loss cone, corresponding to the 90-degree MEPED telescope, is (Kennel and Petschek, 1966) J 90 = ∫ E 1 E 2 ∫ 0 2 π ∫ 0 β S ( E ) 〈 D α α 〉 L C ⁡ cos ⁡ α e q [ I 0 (   z 0 ) z 0 I 1 ( z 0 ) + log ( sin ⁡ α e q sin ( α L C ) ) ]   A ⁡ sin ⁡ η   d η   d ψ   d E . (11)

In both equations, α L C is the loss cone angle at the equator, and the equatorial pitch angle α e q of a particle is related to the local pitch angle α by sin ⁡   2 α e q = B e q B L E O sin ⁡   2 α , (12) where B L E O is the magnetic field strength at the satellite’s location in low Earth orbit, nominally 800 km altitude. The parameter z 0 is related to 〈 D a a 〉 L C , the bounce-averaged pitch angle diffusion coefficient evaluated at the loss cone angle: z 0 = 2 α L C τ B   〈 D α α 〉 L C ⁡ cos ⁡ α e (13)

The integrations in Eqs 10, 11 are over a finite energy range, [ E 1 , E 2 ] , which we will choose to correspond to the energy ranges from MEPED. The lowest energy interval from MEPED is 30–100 keV and has the highest count rates. In Section 4.2 we further develop this scaling method to account for f p e / f c e by also examining the energy range of 100–300 keV. The function S ( E ) is the electron energy spectrum, which we choose as the Kappa distribution used in Ni et al. (2014).

The remaining integration variables ψ and η account for the solid angle of the detector’s field of view. The azimuthal angle ψ , polar angle η , and detector tilt angle θ are related to the particle pitch angle α by the equation cos ⁡ α = cos ⁡ θ ⁡ cos ⁡ η + sin ⁡ θ ⁡ sin ⁡ η ⁡ cos ⁡ ψ . (14)

The tilt angle θ is defined such that θ = 0 ° corresponds to the 0-degree detector looking parallel to the Earth’s magnetic field, and θ = 90 ° corresponds to the 90-degree detector looking perpendicular to the magnetic field. The tilt angle is a known quantity from the MEPED data.

Eq. 14 can be used to write the local pitch angle α as a function of the integration variables, η and ψ , for a known detector tilt angle. The detector is axisymmetric, so ψ ranges from 0 to 2 π . As written, Eqs 10, 11 show the integration over the polar angle η is from 0 to β , where β is the half-width of the detector field of view. The MEPED detectors on POES and MetOp are typically reported to have a detector half-width of β = 15 ° . However, as discussed in Section 3, Selesnick et al. (2020) used Monte Carlo simulations to compute accurate angular response functions for each detector that accounted for energetic electrons penetrating the shielding around the instrument. In taking the polar integrals in Eqs 10, 11 we utilize the tabulated angular response functions from Selesnick et al. (2020). This is done through the substitution d η → R S e l ( η ) d η , (15) where we denote R S e l ( η ) as the angular response function from Selesnick et al. (2020), and we now integrate across the entire polar angle range from 0° to 90°.

The Li et al. (2013) and Ni et al. (2014) method for scaling diffusion coefficients for chorus intensity is as follows:

1) Specify the statistical diffusion coefficient to scale, 〈 D α α 〉 , as a function of electron energy and pitch angle.

This method can be applied without ever knowing the chorus wave intensity used in PADIE to calculate ⟨ D α α ⟩ . The inferred chorus intensity B w [ i n f e r r e d ] in the radiation belts is B w [ i n f e r r e d ] = s r   B w [ i n p u t ] , (16) where B w [ i n p u t ] is the input for PADIE, and the square root comes from the fact that 〈 D α α 〉 ∝ B w 2 .

Figure 3 shows an example of this scaling process, with Eq. 16 used to plot J 0 / J 90 as a function of chorus wave intensity instead of the abstract scaling ratio. In Figure 3 the scaling method is applied using observed flux ratios by MEPED, leading to an inferred chorus intensity of B w = 435 pT from the 30–100 keV band, and an inferred chorus intensity of B w = 825 pT from the 100–300 keV band. Since the observed fluxes are caused by the same chorus waves in the radiation belts, the scaling method should produce the same chorus intensity values. In the next section we show that scaling the diffusion coefficients for f p e / f c e addresses this problem.

The method for scaling chorus intensity using MEPED data produces different estimates of the chorus intensity in the 30–100 keV and 100–300 keV energy ranges (Figure 3). This discrepancy was noticed in Ni et al. (2014), where they attributed the difference to 1) lower uncertainties in the measured 30–100 keV fluxes due to higher count rates, and 2) the increased importance of momentum diffusion at higher energies. We offer an alternative explanation: the diffusion coefficients assume an inaccurate plasma-gyrofrequency ratio, f p e / f c e . This section develops a method to scale diffusion coefficients with f p e / f c e , which we use to force the chorus intensity to be the same in both energy ranges.

Eqs 3, 4, 5 show how the pitch angle diffusion coefficient is calculated. The plasma-gyrofrequency ratio does not explicitly enter into those equations, but it does change the dispersion relation for chorus waves. From Appendix A, the chorus dispersion relation for parallel propagating waves ( X = tan ⁡ ψ = 0 ) is ω 2 c 2 k 2 = ( 1 + ω p 2 Ω c 2 − ω 2 ( 1 + | Ω c | ω ) ) − 1 , (17) which is solved simultaneously with the Landau/cyclotron resonance condition in Eq. 5. For parallel propagating waves the resonance condition yields k = 1 v ∥ ( ω − n Ω c γ ) . (18)

Substituting Eq. 18 into Eq. 17, and utilizing γ = ( 1 − v ∥ 2 c 2 ⁡ cos ⁡   2 α ) − 1 / 2 with α as the pitch angle, the resonant chorus frequency ω is ω 3 Ω c 3 ( γ 2 ⁡ sin ⁡   2 α + cos ⁡   2 α ) + ω 2 Ω c 2 ( γ 2 ⁡ sin ⁡   2 α + cos ⁡   2 α − 2 n γ ) + ω Ω c ( n 2 + ω p 2 Ω c 2 ( γ 2 − 1 ) cos ⁡   2 α − 2 n γ ) + n 2 = 0 (19)

This is a cubic equation in ω / Ω c , and the coefficients for the 0^th, 2^nd and 3^rd order terms do not depend on the plasma frequency. Therefore, the same solution for ω is obtained if ω p and γ are changed such that n 2 + ω p 2 Ω c 2 ( γ 2 − 1 ) cos ⁡   2 α − 2 n γ = c o n s t a n t . (20) For parallel propagating waves only the n = − 1 resonance is important, however we keep n in the equation since it is a constant anyways. Defining the plasma-gyrofrequency ratio as R = ω p Ω c . (21)

Eq. 20 shows that two sets of parameters (subscripts 0 and 1) lead to the same solution of ω / Ω c if R 0 2 ( γ 0 2 − 1 ) cos ⁡   2 α 0 − 2 n γ 0 = R 1 2 ( γ 1 2 − 1 ) cos ⁡   2 α 1 − 2 n γ 1 . (22)

This can be solved for γ 1 with the quadratic equation. However, looking at Eq. 19, we see γ is present in the 2^nd and 3^rd order terms, and therefore changing γ in response to ω p / Ω c will change those coefficients and the overall solution for ω . We therefore constrain Eq. 22 to the weakly-relativistic limit, such that 2 n γ 0 ≈ 2 n γ 1 . Then Eq. 22 simplifies to R 0 2 ( γ 0 2 − 1 ) cos ⁡   2 α 0 = R 1 2 ( γ 1 2 − 1 ) cos ⁡   2 α 1 . (23)

Using γ = 1 + E / E r , where E r is the electron rest mass, we express this relation in terms of particle energy, E 1 ( E 1 + 2 E r ) = E 0 ( E 0 + 2 E r ) R 1 2 ⁡ cos ⁡   2 α 1 R 0 2 ⁡ cos ⁡   2 α 0 , (24) with the quadratic solution E 1 = − E r + R 0 2 ⁡ cos ⁡   2 α 0 R 1 2 ⁡ cos ⁡   2 α 1 E 0 ( E 0 + 2 E r ) + E r 2 . (25) Notice that in the nonrelativistic limit, and for α 1 = α 0 , Eq. 24 simplifies to E 1 = E 0 ( ω p 1 Ω c 1 ) 2 ( ω p 0 Ω c 0 ) − 2 . (26)

Eq. 25 provides an approximate relation between two different particles with different energies and plasma frequencies that yields the same solution of the dispersion relation for the chorus frequency. These two sets of particle energy and plasma-gyrofrequency parameters will therefore yield approximately the same diffusion coefficients, providing a fast and simple way of scaling D α α with changing plasma frequency. This is effectively done by setting α 0 = α 1 , and shifting the entire diffusion coefficient D α α up or down in energy. An example of this scaling with f p e / f c e is shown in Figure 4, where the scaling is done by shifting the whole diffusion coefficient up or down in energy based on Eq. 25.

Eq. 25 gives a quick and approximate way to shift the diffusion coefficients up and down in energy when f p e / f c e changes. We use this relation to scale the diffusion coefficients with the MEPED data by iterating over the method in Section 4.1 until the inferred B w from the 30–100 keV band matches the inferred B w from the 100–300 keV band. This is done as follows:

1) Specify the statistical diffusion coefficient to scale, 〈 D α 〉 .