The motion of energetic particles trapped in planetary radiation belts is a superposition of three quasi-periodic motions, each evolving on a very distinct spatiotemporal scale, with an amplitude quantified by an adiabatic invariant (e.g., Northrop and Teller, 1960; Schulz and Lanzerotti, 1974):

(1) A very fast and small motion of gyration around the magnetic field direction.

The scale separation between these three quasi-periodic motions spans several orders of magnitude in time and space.

Combining adiabatic invariant theory with Fokker–Planck formalism yields the theoretical framework for a probabilistic model of radiation belt dynamics (e.g., Roederer and Zhang, 2014). The Fokker–Planck formalism accounts for uncertainties in electromagnetic field characterization. The adiabatic theory allows for a three-dimensional phase-averaged representation of radiation belt dynamics rather than a full six-dimensional description in phase space.

The description of radiation belt dynamics as a three-dimensional Fokker–Planck equation reduced to a diffusion equation requires minimal computational resources. This quality has enabled the development of many radiation belt computer codes over the years: Salammbô (e.g., Beutier and Boscher, 1995; Nénon et al., 2017), Diffusion in (I,L,B) Energetic Radiation Tracker (DILBERT) (Albert et al., 2009), Versatile Electron Radiation Belt (VERB) (Subbotin and Shprits, 2009), Storm-Time Evolution of Electron Radiation Belt (STEERB) (Su et al., 2010), DREAM3D, as part of the Dynamic Radiation Environment Assimilation Model (DREAM) project (Tu et al., 2013), and British Antarctic Survey Radiation Belt Model (BAS RBM) (Glauert et al., 2014; Woodfield et al., 2014) are all examples of radiation belt codes relying on the same theoretical basis. While first implemented in the case of terrestrial radiation belts, the three-dimensional Fokker–Planck equation has also been transposed to the radiation belts of Jupiter and Saturn. The resulting codes are widely used for scientific research (e.g., Varotsou et al., 2005; Woodfield et al., 2018; Drozdov et al., 2020) and for space weather purposes (e.g., Glauert et al., 2018; Horne et al., 2021).

On the technical side, these computer codes consist of solving a diffusion equation that provides an approximate description for the time evolution of the radiation belts: ∂ f ∂ t = ∑ i , j ∂ ∂ J i D i , j ∂ f ∂ J j + S o u r c e s − L o s s e s , (1) where f t , J 1 , J 2 , J 3 is the phase-averaged phase space density, J i = 1,2,3 are the action variables, which are proportional to the adiabatic invariants by physical constants, and D i , j are the phase-averaged diffusion coefficients. According to Eq. 1, radiation belts are primarily driven by very small, uncorrelated perturbations to the particle trajectories, at all spatiotemporal scales, from the gyro-scale up to the drift scale. The “ S o u r c e s ” and “ L o s s e s ” terms account for other non-diffusive processes affecting the distribution function (e.g., Schulz and Lanzerotti, 1974). It is worth emphasizing that all quantities in Eq. 1 are drift-averaged, i.e., they are phase-averaged over all three phases. It means that this theoretical formulation cannot resolve the drift phase of trapped particles, or equivalently, the magnetic local time ( M L T ) dimension: the resulting modeled radiation belt intensity, f t , J 1 , J 2 , J 3 , is independent of magnetic local time.

From a theoretical standpoint, it is a reasonable first approximation to consider that radiation belt intensity is independent of magnetic local time: any M L T -dependent structure is expected to dissipate rapidly, on a timescale of a few drift periods, because of the mechanism of phase mixing (e.g., Schulz and Lanzerotti, 1974; Ukhorskiy and Sitnov, 2013). Yet, in practice, in situ measurements of trapped particle fluxes suggest that radiation belt intensity frequently depends on the magnetic local time, at least transiently. Both inner and outer terrestrial radiation belt fluxes typically display drift-periodic oscillations. Depending on the situation, these drift-periodic signatures can be interpreted as drift echoes following M L T -localized injections, dropout echoes following M L T -localized losses, or evidence of trapped particles’ drift resonance with ULF waves (e.g., Sauvaud et al., 2013; Hao et al., 2016; Patel et al., 2019; Lejosne and Mozer, 2020; Zhao et al., 2022). Drift echoes have also been reported in Saturnian radiation belt fluxes (e.g., Hao et al., 2020).

In all cases, processes generating drift-periodic signatures are important due to their connection to radiation belt energization (e.g., Hudson et al., 2020). Yet, three-dimensional radiation belt models cannot account for the generation of drift-periodic signatures. Instead, drift-periodic signatures are usually modeled independently of other processes, by tracking the drift motion of test particles (guiding centers) in prescribed electric and magnetic fields, omitting local processes occurring along the gyration and bounce motions (such as local acceleration by chorus waves for instance) (e.g., Li et al., 1993; Hudson et al., 2017).

In that context, it is necessary to introduce a general equation for radiation belt dynamics that includes M L T -localized effects, and that can account for both local processes, at the gyro-scale, and large-scale effects associated with the radial transport. An equation that meets these requirements is detailed in the following section. It relies on the work by Birmingham et al. (1967), in which a two-dimensional drift-diffusion equation was derived assuming conservation of the first two adiabatic invariants. It is straightforward to generalize the proposed equation to include diffusion in the first two adiabatic invariants. We present a compact way to retrieve the equation proposed by Birmingham et al. (1967), combining Fokker–Planck formalism with relationships derived from the Hamiltonian theory. While adjustments to the three-dimensional diffusion Eq. 1 have already been proposed to resolve the drift phase in radiation belt models (e.g., Bourdarie et al., 1997; Shprits et al., 2015) and ring current models can resolve local time (e.g., Jordanova et al., 1997; 2022; Fok et al., 2014), we propose an alternative from the first principles and describe its underlying theoretical assumptions. Similar to the theoretical framework for ring current models (e.g., Fok and Moore, 1997; Yu et al., 2016), the work discussed thereafter relies on the representation of the inner magnetosphere in terms of Euler potentials (e.g., Stern, 1967). That is why the outline of the remainder is as follows: in Section 2, we provide the theoretical background necessary to derive the equation proposed by Birmingham et al. (1967). In particular, we recall how to derive the standard radial diffusion equation before deconstructing it. We introduce the Euler potential coordinates and relate the Euler coordinates to the third adiabatic invariant. In Section 3, we show how the Fokker–Planck equation in terms of Euler potential coordinates yields a two-dimensional drift-diffusion equation when Hamiltonian relationships between the Euler coordinates are taken into account.

Since this work focuses on improving the modeling of drift effects on radiation belt intensity, we first assume conservation of the first two adiabatic invariants. Thus, all considered quantities are bounce-averaged. We also omit any significant source or loss mechanism. A generalization of the resulting transport equation to include diffusion of the first two adiabatic invariants is straightforward. It is provided at the end of Section 3.

We briefly recall how to derive the standard radial diffusion equation. This informs how to derive the same equation as the one proposed by Birmingham et al. (1967) (Section 3). We also detail the concept of Euler potentials and highlight their connection to the third adiabatic invariant.

Here, we present a compact way to retrieve the equation proposed by Birmingham et al. (1967). This equation represents the time evolution of radiation belt intensity due to transport processes. We describe the time evolution of a distribution function, F , that quantifies the number of particles per unit of Euler potential surface d α d β . This function, F , is proportional to the phase space density averaged over both gyration and bounce phases by a physical constant. It relates to the drift-averaged distribution function, f , introduced in Section 2.1, since the number of particles per unit of third invariant, J , is f d J = ∮ β ∈ Γ F d α β d β (with d J = ∮ β ∈ Γ d α β d β . We assume that many very small random changes of the Euler coordinates occur between times t and t + ∆ t , with a very small total effect. The resulting two-dimensional Fokker–Planck equation is ∂ F ∂ t = − ∂ ∂ α ∆ α F − ∂ ∂ β ∆ β F + 1 2 ∂ 2 ∂ α 2 ∆ α 2 F + 1 2 ∂ 2 ∂ β 2 ∆ β 2 F + 1 2 ∂ 2 ∂ α ∂ β ∆ α ∆ β F + 1 2 ∂ 2 ∂ β ∂ α ∆ β ∆ α F , (20) where the angle bracket sign,   , indicates the rate of change of the expected value for the bracketed variable and ∆ X = X t + ∆ t − X t . Just like in Section 2.1 (Eq. 3), we rewrite Eq. 20 as ∂ F ∂ t = ∂ ∂ α − ∆ α F + 1 2 ∂ ∂ α ∆ α 2 F + 1 2 ∂ ∂ β ∆ α ∆ β F + ∂ ∂ β − ∆ β F + 1 2 ∂ ∂ β ∆ β 2 F + 1 2 ∂ ∂ α ∆ β ∆ α F . (21)

The terms between the large parentheses in Eq. 21 are − ∆ α F + 1 2 ∂ ∂ α ∆ α 2 F + 1 2 ∂ ∂ β ∆ α ∆ β F = − ∆ α + 1 2 ∂ ∆ α 2 ∂ α + 1 2 ∂ ∆ α ∆ β ∂ β F + ∆ α 2 2 ∂ F ∂ α + ∆ α ∆ β 2 ∂ F ∂ β (22) and − ∆ β F + 1 2 ∂ ∂ β ∆ β 2 F + 1 2 ∂ ∂ α ∆ β ∆ α F = − ∆ β + 1 2 ∂ ∆ β 2 ∂ β + 1 2 ∂ ∆ β ∆ α ∂ α F + ∆ β 2 2 ∂ F ∂ β + ∆ β ∆ α 2 ∂ F ∂ α . (23)

Using the Hamiltonian relationships between the Euler potentials (Eq. 17), we have shown in the Appendix that − ∆ α + 1 2 ∂ ∆ α 2 ∂ α + 1 2 ∂ ∆ α ∆ β ∂ β = − α ˙ − ∆ β + 1 2 ∂ ∆ β 2 ∂ β + 1 2 ∂ ∆ β ∆ α ∂ α = − β ˙ , (24) provided that the time interval, ∆ t , is very small in comparison with the characteristic time for the time variation of the Hamiltonian ( ∆ t ≪ H / ∂ H / ∂ t . In practice, the time interval, ∆ t , is of the order of a few bounce periods, that is, very small in comparison with the drift period. The squared brackets,   , indicate the expected value of the bracketed variable.

Combining Eqs 21–24, Eq. 20 becomes a drift-diffusion equation: ∂ F ∂ t = − ∂ ∂ α α ˙ F − ∂ ∂ β β ˙ F + ∂ ∂ α ∆ α 2 2 ∂ F ∂ α + ∂ ∂ α ∆ α ∆ β 2 ∂ F ∂ β + ∂ ∂ β ∆ β ∆ α 2 ∂ F ∂ α + ∂ ∂ β ∆ β 2 2 ∂ F ∂ β . (25)

Given Eq. 19, this simplifies to ∂ F ∂ t = − α ˙ ∂ F ∂ α − β ˙ ∂ F ∂ β + ∂ ∂ α D α α ∂ F ∂ α + ∂ ∂ α D α β ∂ F ∂ β + ∂ ∂ β D β α ∂ F ∂ α + ∂ ∂ β D β β ∂ F ∂ β , (26) where D α α = ∆ α 2 / 2 , D β β = ∆ β 2 / 2 , D α β = ∆ α ∆ β / 2 , and D β α = ∆ β ∆ α / 2 are the diffusion coefficients, and α ˙ and β ˙ are the mean bounce-averaged time rates of change of α and β , respectively. This transport equation coincides with the one provided by Birmingham et al. (1967), their equation (4.11). A change of variables (using Eq. 15) yields: ∂ F ∂ t = − L ˙ ∂ F ∂ L − φ ˙ ∂ F ∂ φ + L 2 ∂ ∂ L D L L L 2 ∂ F ∂ L + L 2 ∂ ∂ L D L φ L 2 ∂ F ∂ φ + ∂ ∂ φ D φ L ∂ F ∂ L + ∂ ∂ φ D φ φ ∂ F ∂ φ . (27)

The term depending on D L L mistakenly resembles the one present in the standard radial diffusion equation (Eq. 7): D L L and D L L are different. The distribution function, f , the coefficient for the standard radiation diffusion equation (Eq. 7), D L L , and more generally, the quantities used for the three-dimensional equation for radiation belt dynamics (Eq. 1) are drift-averaged, i.e., they are independent of the drift phase. Here, the drift phase is resolved: the distribution function and coefficients are bounce-averaged quantities that depend on the drift phase. Thus, they must be evaluated at each location ( α , β ), or similarly ( L , φ ), and at each time, t .

The transport parameters of Eq. 27 are all statistically averaged quantities. The coefficients L ˙ and φ ˙ (or equivalently α ˙ and β ˙ ) indicate ensemble averages of time derivatives for the quantities considered. The ensemble averages are computed at each location and at each time, t , over an ensemble of field fluctuations. The diffusion coefficients are proportional to the time rates of change of the covariances for the quantities considered. Specifically, when considering two variables X and Y (where X , Y could be any combination of α , β or L , φ ), the diffusion coefficient is D X Y = X t + ∆ t − X t Y t + ∆ t − Y t 2 ∆ t . (28)

That is, it is half the time rate of change of the ensemble average for the product of the time variations of X and Y during a time interval, ∆ t . A worked example will be provided in the second part of this work. It will detail how to compute all transport parameters of Eq. 27 in a particular model of field fluctuations.

According to Eq. 26, variations in the distribution function are due to the bulk motion of the plasma in the presence of density gradients and to diffusive effects in both the localized radial ( L ) and azimuthal ( φ ) directions. Local effects acting at smaller scales can be readily reinstated by adding relevant coefficients modeling local diffusion, source, and loss mechanisms: ∂ F ∂ t = − L ˙ ∂ F ∂ L − φ ˙ ∂ F ∂ φ + L 2 ∂ ∂ L D L L L 2 ∂ F ∂ L + L 2 ∂ ∂ L D L φ L 2 ∂ F ∂ φ + ∂ ∂ φ D φ L ∂ F ∂ L + ∂ ∂ φ D φ φ ∂ F ∂ φ + ∑ 1 ≤ i , j ≤ 2 ∂ ∂ J i D i , j ∂ F ∂ J j + S o u r c e s − L o s s e s , (29) where all quantities are bounce-averaged quantities that depend on the drift phase.

The objective of this work is to contribute toward improving the spatiotemporal resolution of physics-based diffusive radiation belt models. The resulting transport Eq. 27 can resolve the drift phase, and the outputs are bounce-averaged rather than drift-averaged. This is of use when the objective is to model fast radiation belt dynamics, such as times of fast radiation belt acceleration or losses occurring during the main phase of geomagnetic storms (e.g., Ripoll et al., 2020; Lejosne et al., 2022). It can also be used to increase the energy range modeled, by including ring current energies.

Although Eq. 27 contains some localized (in L , MLT) diffusion coefficients, its scope is beyond the long-established radial diffusion paradigm used to summarize transport effects on radiation belt intensity. The inclusion of the effects of bulk motion and the diffusion in the azimuthal coordinate enable the modeling of MLT-localized structures, drift-periodic flux oscillations, and their subsequent attenuation due to phase-mixing processes.

Current works leveraging in situ measurements to quantify radial diffusion coefficients require information on average over all magnetic local times of a drift shell. Yet, a spacecraft can only scan the electromagnetic environment along its orbit, limiting the accuracy with which the outputs can be determined (e.g., Sandhu et al., 2021). Because the coefficients introduced in this work depend on magnetic local time, these may be easier to quantify experimentally. Furthermore, describing the effect of drift motion on radiation belt intensity in terms of Euler potentials, or similarly with ( L , φ ), is computationally more advantageous than working with the action-angle variables ( J 3 , φ 3 ): the latter requires tracing the instantaneous drift contour at every time step, while the former only requires local field line tracing. In addition, the definition of the Euler potentials only requires closed field lines, while the definition of the action-angle variables is more restrictive, requiring a closed instantaneous drift contour. Thus, working in terms of Euler potentials allows for the inclusion of quasi-trapped particles from the drift loss cone.

The second part of this work will deal with characterizing the coefficients introduced in Eq. 27 (i.e., L ˙ , φ ˙ , D L L , D L φ , D φ L , D φ φ ) in the special case of electric potential fluctuations in a magnetic dipole field. It will show how to implement the theoretical framework presented in this work.