Drift phase resolved diffusive radiation belt model: 1. Theoretical framework

Most physics-based models provide a coarse three-dimensional representation of radiation belt dynamics at low time resolution, of the order of a few drift periods. The description of the effect of trapped particle transport on radiation belt intensity is based on the random phase approximation, and it is in one dimension only: the third adiabatic invariant coordinate, akin to a phase-averaged radial distance. This means that these radiation belt models do not resolve the drift phase or, equivalently, the magnetic local time. Yet, in situ measurements suggest that radiation belt intensity frequently depends on magnetic local time, at least transiently, such as during active times. To include processes generating azimuthal variations in trapped particle fluxes and to quantify their relative importance in radiation belt energization, an improvement in the spatiotemporal resolution of the radiation belt models is required. The objective of this study is to pave the way for a new generation of diffusive radiation belt models capable of retaining drift phase information. Specifically, we highlight a two-dimensional equation for the effects of trapped particle transport on radiation belt intensity. With a theoretical framework that goes beyond the radial diffusion paradigm, the effects of trapped particle bulk motion, as well as diffusion, are quantified in terms of Euler potentials, $\left(\alpha ,\beta \right)$, quantities akin to the radial and azimuthal directions. This work provides the theoretical foundations underlying the drift phase resolved transport equation for radiation belt dynamics. It also brings forward the concept of azimuthal diffusion as a phase-mixing agent.

1 Introduction

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.

(2) A slower and bigger bounce motion between the planet’s hemispheres, along the magnetic field direction.

(3) A slow and large drift motion around the planet in a direction perpendicular to 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:

$\begin{array}{c}\frac{\partial f}{\partial t}={\displaystyle \sum }_{i,j}\frac{\partial }{\partial J_i}\left(D_{i,j}\frac{\partial f}{\partial J_j}\right)+Sources-Losses,\end{array}(1)$

where $f\left(t,J_1,J_2,J_3\right)$ is the phase-averaged phase space density, $J_{i=\mathrm{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 “$Sources$” and “$Losses$” 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 ($MLT$) dimension: the resulting modeled radiation belt intensity, $f\left(t,J_1,J_2,J_3\right)$, 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 $MLT$-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 $MLT$-localized injections, dropout echoes following $MLT$-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 $MLT$-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.

2 Theoretical background

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.

2.1 Derivation of the standard radiation belt radial diffusion equation

In the following section, the third adiabatic invariant, $J_3$, is abbreviated to $J$ out of convenience. The objective is to describe the time evolution of a distribution function, $f$, that quantifies the number of particles per unit of $J$ (assuming conservation of the first two adiabatic invariants). This quantity is proportional to the drift-averaged phase space density by a physical constant (e.g., Roederer and Zhang, 2014, their chapter 4). The usual assumption is that many very small uncorrelated random changes of the variable, $J$, occur between times $t$ and $t+∆t$, with a very small total effect $\left(∆J/J\ll 1;∆t\ll f/\left(\partial f/\partial t\right)\right)$. In this case, the time evolution of the distribution function, $f,$ is provided by a Fokker–Planck equation (e.g., Roederer, 1970; Walt, 1994):

$\begin{array}{c}\frac{\partial f}{\partial t}=-\frac{\partial }{\partial J}\left(〈∆J〉f\right)+\frac{1}{2}\frac{{\partial }^2}{\partial J^2}\left(〈{\left(∆J\right)}^2〉f\right),\end{array}(2)$

where $〈∆J〉=\left[∆J\right]/∆t$ is the rate of change for the expected value of the third invariant variation, $\left[∆J\right]=\left[J\left(t+∆t\right)-J\left(t\right)\right]$, and $〈{\left(∆J\right)}^2〉=\left[{\left(∆J\right)}^2\right]/∆t$ is the rate of change for the expected value of the third invariant squared variation. A rewriting of the right-hand side of Eq. 2 provides a mathematically equivalent formulation:

$\begin{array}{c}\frac{\partial f}{\partial t}=\frac{\partial }{\partial J}\left(-〈∆J〉f+\frac{1}{2}\frac{\partial }{\partial J}\left(〈{\left(∆J\right)}^2〉f\right)\right),\end{array}(3)$

which can also be written as

$\begin{array}{c}\frac{\partial f}{\partial t}=\frac{\partial }{\partial J}\left(\left(-〈∆J〉+\frac{1}{2}\frac{\partial 〈{\left(∆J\right)}^2〉}{\partial J}\right)f+\frac{〈{\left(∆J\right)}^2〉}{2}\frac{\partial f}{\partial J}\right)\end{array}.(4)$

To transform this equation into a radial diffusion equation, we use the fact that the two coefficients $〈∆J〉$ and $〈{\left(∆J\right)}^2〉$ are not independent of each other:

$\begin{array}{c}〈∆J〉=\frac{1}{2}\frac{\partial 〈{\left(∆J\right)}^2〉}{\partial J}\end{array}(5)$

(e.g., Lichtenberg and Lieberman, 1992, their section 5.4a; Lejosne and Kollmann, 2020, their section 2.3.2). This relationship (Eq. 5) relies on the assumption of drift phase homogeneity, also known as random drift phase approximation, meaning that each drift phase location is equiprobable. In this context, the Fokker–Planck Eq. 2 reduces to a diffusion equation:

$\begin{array}{c}\frac{\partial f}{\partial t}=\frac{\partial }{\partial J}\left(D\frac{\partial f}{\partial J}\right),\end{array}(6)$

where $D=〈{\left(∆J\right)}^2〉/2$ is the diffusion coefficient in $J$. The diffusion equation is often rewritten in terms of $L^{*}\propto 1/J$:

$\begin{array}{c}\frac{\partial f}{\partial t}=L^{*2}\frac{\partial }{\partial L^{*}}\left(\frac{D_{LL}}{L^{*2}}\frac{\partial f}{\partial L^{*}}\right),\end{array}(7)$

where $D_{LL}=〈{\left(∆L^{*}\right)}^2〉/2$ is the radial diffusion coefficient.

2.2 Euler potentials

An appropriate coordinate to discuss radial diffusion is the $L^{*}$ coordinate (Roederer, 1967), inversely proportional to the third adiabatic invariant, $J$ (Eq. 7). In the following section, we argue that $L^{*}$ is not suited when the objective is to resolve the drift phase. Instead, we introduce the best-suited coordinate, $\mathbb{L}$ (“double-struck L” or “L-Euler”). We discuss the relationship between $\mathbb{L}$ and $L^{*}$ by detailing the underlying role of the Euler potentials.

2.2.1 Third adiabatic invariant, deconstructed in terms of Euler potentials

The radial diffusion equation, retrieved in Section 2.1 (Eq. 7), describes the time evolution of the number of particles per unit of third adiabatic invariant, $J$, or equivalently, $L^{*}.$ The quantities $J$ and $L^{*}$ are MLT-averaged by design. Indeed, the third invariant of a trapped population, $J$, is proportional to the magnetic flux encompassed by the guiding drift shell:

$\begin{array}{c}J\propto {\iint }_{\mathrm{\Sigma }}\mathit{B}\bullet \mathit{d}\mathit{S}={\oint }_{\mathrm{\Gamma }}\mathit{A}\bullet \mathit{d}\mathit{l},\end{array}(8)$

where $\mathit{A}$ is the magnetic vector potential ($\nabla \times \mathit{A}=\mathit{B}$), and $\mathrm{\Sigma }$ is the surface encompassed by the instantaneous drift contour, $\mathrm{\Gamma }$, of the trapped population. The instantaneous drift contour, $\mathrm{\Gamma }$, can be viewed as the intersection of the guiding drift shell with a surface, such as the minimum B-surface (see also, Roederer, 1970, p. 76–79). In other words, to quantify the third adiabatic invariant, $J$, it is necessary to know the guiding drift shell, that is, the set of guiding center locations at all magnetic local times, treating the electromagnetic fields as stationary.

An important underlying requirement to sort trapped particle fluxes using the third adiabatic invariant is the so-called frozen field condition, where in the presence of magnetic field time variations, the cold (frozen) plasma $\mathit{E}\times \mathit{B}$ drifts to remain on the same magnetic field line (Birmingham and Jones, 1968). This assumption requires the Earth’s surface to be a perfect conductor and no component of the electric field to be parallel to the magnetic field direction. In this context, the footpoints of a magnetic field line are rooted at fixed locations at ionospheric altitudes, while the rest of the field line can “move” (stretch, compress, distort) in the magnetosphere in the presence of magnetic field time variations. Thus, the frozen field condition enables a tempting, yet disputed, concept of field line “flagging” and its corollary, field line “motion” (Fälthammar and Mozer, 2007). It is indeed worth emphasizing that a field line is an imaginary concept that aids to visualize the magnitude and direction of a vector field, so there should be no way of differentiating a field line from the other. We assume nonetheless that we can label field lines based on the locations of their rooted ionospheric footpoints. In this context, to determine a guiding drift shell or an instantaneous drift contour, $\mathrm{\Gamma }$, and to compute the third adiabatic invariant, we now have to know the set of field lines that were scanned by the drifting guiding centers at all magnetic local times. In other words, we need information on the field line label at each magnetic local time. This can be done by leveraging the Euler potentials, as discussed in the following.

The Euler potentials $\left(\alpha ,\beta \right)$ are a convenient tool for labeling field lines. They are analogous to the stream function in an incompressible flow in fluid mechanics. They offer a representation of the magnetic field intrinsically dependent on its topology (e.g., Stern, 1967, 1970). Their characterization relies on the fact that the magnetic field is a solenoidal vector field, i.e., $\nabla \bullet \mathit{B}=0$. The Euler potentials are such that

$\begin{array}{c}\mathit{B}=\nabla \alpha \times \nabla \beta \end{array}.(9)$

Thus, the Euler potentials are constant along the magnetic field lines. Since the vector potential can be viewed as $\mathit{A}=\alpha \nabla \beta$, a reformulation of Eq. 8 in terms of Euler potential yields

$\begin{array}{c}J\propto {\oint }_{\mathrm{\Gamma }}\alpha d\beta .\end{array}(10)$

Although there is no uniformity in the definition of the Euler potentials, a suitable set of Euler potentials in a magnetic dipole field is

$\begin{array}{c}\left\{\begin{array}{c}\alpha =-\frac{B_E{R}_E^3}{r}{\sin }^2\theta \\ \beta =\phi ,\end{array}\right.\end{array}(11)$

where $B_E=\mathrm{30,000}$ nT is the magnetic equatorial field at the surface of the Earth, $R_E=6370km$ is one Earth radius, and ($r,\theta ,\phi$) are the radial distance, magnetic colatitude, and azimuthal (i.e., MLT) location with respect to the center of the dipole magnetic moment, respectively.

In the presence of a distorted magnetic field, the expressions provided in Eq. 11 are not valid anymore. That said, it is possible to leverage the facts that (a) the field line footpoints are rooted at ionospheric altitudes, a region where the ambient magnetic field is mainly dipolar, so the Euler potentials can be described by Eq. 11 at ionospheric altitudes and (b) the Euler potentials are constant along the magnetic field lines. With that in mind, we can define a set of Euler potentials $\left(\alpha ,\beta \right)$ such that at the footpoints ($R_E,{\theta }_E$, ${\phi }_E$), and thus all along the field lines:

$\begin{array}{c}\left\{\begin{array}{c}\alpha =-B_E{R}_E^2{\sin }^2{\theta }_E\\ \beta ={\phi }_E,\end{array}\right.\end{array}(12)$

where (${\theta }_E$, ${\phi }_E$), respectively, indicate the magnetic colatitude and longitude of the field line footpoint at $r=R_E$, the Earth’s surface.

If a distorted magnetic field were to change into a dipole field, each field line would “move” in geospace, adopting a dipolar shape, while its footpoints would stay rooted at fixed ionospheric latitudes. Leveraging Eq. 11 in the newly transformed dipole field, a dipolar field line with footpoints at ($R_E,{\theta }_E$, ${\phi }_E$) would have its equatorial apex $\left(r_o,\theta =\pi /2,{\phi }_o\right)$ such that $\alpha \left(r_o,\theta =\pi /2,{\phi }_o\right)=\alpha \left(R_E,{\theta }_E,{\phi }_E\right)$ and $\beta \left(r_o,\theta =\pi /2,{\phi }_o\right)=\beta \left(R_E,{\theta }_E,{\phi }_E\right)$. Thus, the intersection of the dipolar field line footpoint and the magnetic equator ($\theta =\pi /2$) would be at

$\begin{array}{c}\left\{\begin{array}{c}{r}_o=\frac{R_E}{{\sin }^2{\theta }_E}\\ {\phi }_o={\phi }_E.\end{array}\right.\end{array}(13)$

The physical interpretation of this thought experiment is similar to the physical interpretation of the $L^{*}$ parameter. The $L^{*}$ coordinate corresponds to the normalized equatorial radius of the circular guiding contour on which trapped particles would drift after all non-dipolar contributions to the magnetic field and all electric field components have been turned off adiabatically. Here, we introduce the parameter $\mathbb{L}$ (“double-struck L” or “L-Euler”) such that

$\begin{array}{c}\mathbb{L}=\frac{1}{{\sin }^2{\theta }_E}\end{array},(14)$

where ${\theta }_E$ is the magnetic colatitude of the footpoint at $r=R_E$ for the field line passing through the location considered. It corresponds to the normalized equatorial radius of the field line on which trapped particles would bounce if all non-dipolar contributions to the magnetic field were turned off relatively fast (a few bounce periods). As for the angle variable, $\beta$, one can reasonably assume no significant longitudinal bending of the field lines when the magnetic field is stretched or compressed. Thus, in terms of Euler potentials, we have in general that

$\begin{array}{c}\left\{\begin{array}{c}\alpha =-\frac{B_E{R}_E^2}{\mathbb{L}}\\ \beta ={\phi }_E\cong \phi .\end{array}\right.\end{array}(15)$

Combining Eqs 10, 15, given that ${\oint }_{\mathrm{\Gamma }}\alpha d\beta =-2\pi B_E{R}_E^2/L^{*}$, we obtain

$\begin{array}{c}\frac{1}{L^{*}}=\frac{1}{2\pi }{\oint }_{\mathrm{\Gamma }}\frac{d\phi }{\mathbb{L}}.\end{array}(16)$

The parameter $L^{*}$ is the harmonic mean of the $\mathbb{L}$ coordinate along the guiding contour, $\mathrm{\Gamma }$, a relationship that can be utilized to quantify $L^{*}$ (e.g., Lejosne, 2014). In the presence of quasi-trapped particles, i.e., guiding centers drifting along on open drift contour, the parameter $L^{*}$ cannot be defined. On the other hand, the $\mathbb{L}$ coordinate can still be defined on open drift contour, as long as we are dealing with a closed field line.

An illustration to the concepts discussed here is provided in Figure 1.

FIGURE 1

FIGURE 1. (A) A few magnetic field lines constitutive of a trapped particle drift shell, together with the locations of the field line footpoints, necessary to determine the L* parameter—adapted from Roederer (1970). (B) A stretch magnetic field line (in red) is relaxed into its dipolar shape (in blue). The field line footpoint remains rooted at the same ionospheric colatitude, ${\theta }_E$, a parameter that determines the $\mathbb{L}$ parameter, of use for field line labeling. The parameter $\mathbb{L}$ corresponds to the normalized equatorial radius of the dipolar field line (in blue)—adapted from Lejosne (2014).

2.2.2 Euler potentials as appropriate variables to describe bounce-average drift motion of trapped and quasi-trapped particles

The Euler potentials $\alpha$ and $\beta$ are proportional to canonical variables, that is,

$\begin{array}{c}\left\{\begin{array}{c}\dot{\alpha }=-\frac{\partial H}{\partial \beta }\\ \dot{\beta }=\frac{\partial H}{\partial \alpha },\end{array}\right.\end{array}(17)$

where $H$ is a Hamiltonian proportional to the total energy of the guiding center (Northrop and Teller, 1960; Birmingham et al., 1967):

$\begin{array}{c}H=\frac{T}{q}+V,\end{array}(18)$

where $T$ is the kinetic energy, q is the charge of the population considered, and $V$ is the electric potential.

One consequence of Eq. 17 is that the variations of the Euler potentials are related:

$\begin{array}{c}\frac{\partial }{\partial \alpha }\left(\dot{\alpha }\right)+\frac{\partial }{\partial \beta }\left(\dot{\beta }\right)=0.\end{array}(19)$

This property will be leveraged to transform a two-dimensional Fokker–Planck equation in terms of Euler potentials, $\left(\alpha ,\beta \right)$, in a two-dimensional drift-diffusion equation.

3 New derivation of Birmingham et al.’s transport equation to describe trapped particle transport effects on radiation belt intensity

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\alpha d\beta$. 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 $fdJ={\oint }_{\mathrm{\beta }\in \mathrm{\Gamma }}Fd\alpha \left(\beta \right)d\beta$ (with $dJ={\oint }_{\mathrm{\beta }\in \mathrm{\Gamma }}\left.d\alpha \left(\beta \right)d\beta \right)$. 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

$\begin{array}{ll}\hfill \frac{\partial F}{\partial t}& =-\frac{\partial }{\partial \alpha }\left(〈∆\alpha 〉F\right)-\frac{\partial }{\partial \beta }\left(〈∆\beta 〉F\right)\hfill \\ \hfill & +\frac{1}{2}\frac{{\partial }^2}{\partial {\alpha }^2}\left(〈{\left(∆\alpha \right)}^2〉F\right)+\frac{1}{2}\frac{{\partial }^2}{\partial {\beta }^2}\left(〈{\left(∆\beta \right)}^2〉F\right)\hfill \\ \hfill & +\frac{1}{2}\frac{{\partial }^2}{\partial \alpha \partial \beta }\left(〈∆\alpha ∆\beta 〉F\right)+\frac{1}{2}\frac{{\partial }^2}{\partial \beta \partial \alpha }\left(〈∆\beta ∆\alpha 〉F\right)\hfill \end{array},(20)$

where the angle bracket sign, $〈〉$, indicates the rate of change of the expected value for the bracketed variable and $∆X=X\left(t+∆t\right)-X\left(t\right)$. Just like in Section 2.1 (Eq. 3), we rewrite Eq. 20 as

$\begin{array}{ll}\hfill \frac{\partial F}{\partial t}& =\frac{\partial }{\partial \alpha }\left(-〈∆\alpha 〉F+\frac{1}{2}\frac{\partial }{\partial \alpha }\left(〈{\left(∆\alpha \right)}^2〉F\right)+\frac{1}{2}\frac{\partial }{\partial \beta }\left(〈∆\alpha ∆\beta 〉F\right)\right)\hfill \\ \hfill & +\frac{\partial }{\partial \beta }\left(-\left(〈∆\beta 〉F\right)+\frac{1}{2}\frac{\partial }{\partial \beta }\left(〈{\left(∆\beta \right)}^2〉F\right)+\frac{1}{2}\frac{\partial }{\partial \alpha }\left(〈∆\beta ∆\alpha 〉F\right)\right).\hfill \end{array}(21)$

The terms between the large parentheses in Eq. 21 are

$\begin{array}{ll}\hfill & \hfill -〈∆\alpha 〉F+\frac{1}{2}\frac{\partial }{\partial \alpha }\left(〈{\left(∆\alpha \right)}^2〉F\right)+\frac{1}{2}\frac{\partial }{\partial \beta }\left(〈∆\alpha ∆\beta 〉F\right)\\ \hfill & \hfill =\left(-〈∆\alpha 〉+\frac{1}{2}\frac{\partial 〈{\left(∆\alpha \right)}^2〉}{\partial \alpha }+\frac{1}{2}\frac{\partial 〈∆\alpha ∆\beta 〉}{\partial \beta }\right)F+\frac{〈{\left(∆\alpha \right)}^2〉}{2}\frac{\partial F}{\partial \alpha }+\frac{〈∆\alpha ∆\beta 〉}{2}\frac{\partial F}{\partial \beta }\end{array}(22)$

and

$\begin{array}{ll}\hfill & \hfill -\left(〈∆\beta 〉F\right)+\frac{1}{2}\frac{\partial }{\partial \beta }\left(〈{\left(∆\beta \right)}^2〉F\right)+\frac{1}{2}\frac{\partial }{\partial \alpha }\left(〈∆\beta ∆\alpha 〉F\right)\\ \hfill & \hfill =\left(-〈∆\beta 〉+\frac{1}{2}\frac{\partial 〈{\left(∆\beta \right)}^2〉}{\partial \beta }+\frac{1}{2}\frac{\partial 〈∆\beta ∆\alpha 〉}{\partial \alpha }\right)F+\frac{〈{\left(∆\beta \right)}^2〉}{2}\frac{\partial F}{\partial \beta }+\frac{〈∆\beta ∆\alpha 〉}{2}\frac{\partial F}{\partial \alpha }.\end{array}(23)$

Using the Hamiltonian relationships between the Euler potentials (Eq. 17), we have shown in the Appendix that

$\begin{array}{c}\left\{\begin{array}{c}-〈∆\alpha 〉+\frac{1}{2}\frac{\partial 〈{\left(∆\alpha \right)}^2〉}{\partial \alpha }+\frac{1}{2}\frac{\partial 〈∆\alpha ∆\beta 〉}{\partial \beta }=-\left[\dot{\alpha }\right]\\ -〈∆\beta 〉+\frac{1}{2}\frac{\partial 〈{\left(∆\beta \right)}^2〉}{\partial \beta }+\frac{1}{2}\frac{\partial 〈∆\beta ∆\alpha 〉}{\partial \alpha }=-\left[\dot{\beta }\right],\end{array}\right.\end{array}(24)$

provided that the time interval, $∆t,$ is very small in comparison with the characteristic time for the time variation of the Hamiltonian ($\left.∆t\ll H/\left(\partial H/\partial t\right)\right)$. 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, $\left[\right]$, indicate the expected value of the bracketed variable.

Combining Eqs 21–24, Eq. 20 becomes a drift-diffusion equation:

$\begin{array}{ll}\hfill \frac{\partial F}{\partial t}& =-\frac{\partial }{\partial \alpha }\left(\left[\dot{\alpha }\right]F\right)-\frac{\partial }{\partial \beta }\left(\left[\dot{\beta }\right]F\right)\hfill \\ \hfill & +\frac{\partial }{\partial \alpha }\left(\frac{〈{\left(∆\alpha \right)}^2〉}{2}\frac{\partial F}{\partial \alpha }\right)+\frac{\partial }{\partial \alpha }\left(\frac{〈∆\alpha ∆\beta 〉}{2}\frac{\partial F}{\partial \beta }\right)\hfill \\ \hfill & +\frac{\partial }{\partial \beta }\left(\frac{〈∆\beta ∆\alpha 〉}{2}\frac{\partial F}{\partial \alpha }\right)+\frac{\partial }{\partial \beta }\left(\frac{〈{\left(∆\beta \right)}^2〉}{2}\frac{\partial F}{\partial \beta }\right).\hfill \end{array}(25)$

Given Eq. 19, this simplifies to

$\begin{array}{ll}\hfill \frac{\partial F}{\partial t}& =-\left[\dot{\alpha }\right]\frac{\partial F}{\partial \alpha }-\left[\dot{\beta }\right]\frac{\partial F}{\partial \beta }\hfill \\ \hfill & +\frac{\partial }{\partial \alpha }\left(D_{\alpha \alpha }\frac{\partial F}{\partial \alpha }\right)+\frac{\partial }{\partial \alpha }\left(D_{\alpha \beta }\frac{\partial F}{\partial \beta }\right)\hfill \\ \hfill & +\frac{\partial }{\partial \beta }\left(D_{\beta \alpha }\frac{\partial F}{\partial \alpha }\right)+\frac{\partial }{\partial \beta }\left(D_{\beta \beta }\frac{\partial F}{\partial \beta }\right),\hfill \end{array}(26)$

where $D_{\alpha \alpha }=〈{\left(∆\alpha \right)}^2〉/2$, $D_{\beta \beta }=〈{\left(∆\beta \right)}^2〉/2$, $D_{\alpha \beta }=〈∆\alpha ∆\beta 〉/2$, and $D_{\beta \alpha }=〈∆\beta ∆\alpha 〉/2$ are the diffusion coefficients, and $\left[\dot{\alpha }\right]$ and $\left[\dot{\beta }\right]$ are the mean bounce-averaged time rates of change of $\alpha$ and $\beta$, 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:

$\begin{array}{ll}\hfill \frac{\partial F}{\partial t}& =-\left[\dot{\mathbb{L}}\right]\frac{\partial F}{\partial \mathbb{L}}-\left[\dot{\phi }\right]\frac{\partial F}{\partial \phi }\hfill \\ \hfill & +{\mathbb{L}}^2\frac{\partial }{\partial \mathbb{L}}\left(\frac{D_{\mathbb{L}\mathbb{L}}}{{\mathbb{L}}^2}\frac{\partial F}{\partial \mathbb{L}}\right)+{\mathbb{L}}^2\frac{\partial }{\partial \mathbb{L}}\left(\frac{D_{\mathbb{L}\phi }}{{\mathbb{L}}^2}\frac{\partial F}{\partial \phi }\right)\hfill \\ \hfill & +\frac{\partial }{\partial \phi }\left(D_{\phi \mathbb{L}}\frac{\partial F}{\partial \mathbb{L}}\right)+\frac{\partial }{\partial \phi }\left(D_{\phi \phi }\frac{\partial F}{\partial \phi }\right)\hfill \end{array}.(27)$

The term depending on $D_{\mathbb{L}\mathbb{L}}$ mistakenly resembles the one present in the standard radial diffusion equation (Eq. 7): $D_{\mathbb{L}\mathbb{L}}$ and $D_{LL}$ are different. The distribution function, $f$, the coefficient for the standard radiation diffusion equation (Eq. 7), $D_{LL}$, 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 ($\alpha ,\beta$), or similarly ($\mathbb{L},\phi$), and at each time, $t$.

The transport parameters of Eq. 27 are all statistically averaged quantities. The coefficients $\left[\dot{\mathbb{L}}\right]$ and $\left[\dot{\phi }\right]$ (or equivalently $\left[\dot{\alpha }\right]$ and $\left[\dot{\beta }\right]$) 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 $\left(X,Y\right)$ could be any combination of $\left(\alpha ,\beta \right)$ or $\left(\mathbb{L},\phi \right)$), the diffusion coefficient is

$\begin{array}{c}{D}_{XY}=\frac{\left[\left(\mathrm{X}\left(\mathrm{t}+∆\mathrm{t}\right)-\mathrm{X}\left(\mathrm{t}\right)\right)\left(\mathrm{Y}\left(\mathrm{t}+∆\mathrm{t}\right)-\mathrm{Y}\left(\mathrm{t}\right)\right)\right]}{2∆\mathrm{t}}\end{array}.(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, $∆\mathrm{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 ($\mathbb{L}$) and azimuthal ($\phi$) directions. Local effects acting at smaller scales can be readily reinstated by adding relevant coefficients modeling local diffusion, source, and loss mechanisms:

$\begin{array}{ll}\hfill \frac{\partial F}{\partial t}& =-\left[\dot{\mathbb{L}}\right]\frac{\partial F}{\partial \mathbb{L}}-\left[\dot{\phi }\right]\frac{\partial F}{\partial \phi }\hfill \\ \hfill & +{\mathbb{L}}^2\frac{\partial }{\partial \mathbb{L}}\left(\frac{D_{\mathbb{L}\mathbb{L}}}{{\mathbb{L}}^2}\frac{\partial F}{\partial \mathbb{L}}\right)+{\mathbb{L}}^2\frac{\partial }{\partial \mathbb{L}}\left(\frac{D_{\mathbb{L}\phi }}{{\mathbb{L}}^2}\frac{\partial F}{\partial \phi }\right)\hfill \\ \hfill & +\frac{\partial }{\partial \phi }\left(D_{\phi \mathbb{L}}\frac{\partial F}{\partial \mathbb{L}}\right)+\frac{\partial }{\partial \phi }\left(D_{\phi \phi }\frac{\partial F}{\partial \phi }\right)\hfill \\ \hfill & +{\displaystyle \sum }_{1\le i,j\le 2}\frac{\partial }{\partial J_i}\left(D_{i,j}\frac{\partial F}{\partial J_j}\right)+Sources-Losses\hfill \end{array},(29)$

where all quantities are bounce-averaged quantities that depend on the drift phase.

4 Conclusion

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 $\mathbb{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 ($\mathbb{L},\phi$), is computationally more advantageous than working with the action-angle variables ($J_3,$${\phi }_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., $\left[\dot{\mathbb{L}}\right],\left[\dot{\phi }\right],D_{\mathbb{L}\mathbb{L}},D_{\mathbb{L}\phi },D_{\phi \mathbb{L}},D_{\phi \phi }$) 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.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

We describe contributions to the paper using the CRediT (Contributor Roles Taxonomy) categories (Brand et al., 2015). Conceptualization, writing—original draft, and writing—review and editing: all authors. Visualization: SL. All authors contributed to the article and approved the submitted version.

Acknowledgments

SL’s work was performed under NASA Grants 80NSSC18K1223 and 80NSSC20K1351. The author thanks S.D. Walton for helpful comments on a draft version of this manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Albert, J. M., Meredith, N. P., and Horne, R. B. (2009). Three-dimensional diffusion simulation of outer radiation belt electrons during the 9 October 1990 magnetic storm. J. Geophys. Res.114, A09214. doi:10.1029/2009JA014336

CrossRef Full Text | Google Scholar

Beutier, T., and Boscher, D. (1995). A three-dimensional analysis of the electron radiation belt by the Salammbô code. J. Geophys. Res.100 (8), 14853–14861. doi:10.1029/94JA03066

CrossRef Full Text | Google Scholar

Birmingham, T. J., and Jones, F. C. (1968). Identification of moving magnetic field lines. J. Geophys. Res.73 (17), 5505–5510. doi:10.1029/JA073i017p05505

CrossRef Full Text | Google Scholar

Birmingham, T. J., Northrop, T. G., and Fälthammar, C. G. (1967). Charged particle diffusion by violation of the third adiabatic invariant. Phys. Fluids10 (11), 2389. doi:10.1063/1.1762048

CrossRef Full Text | Google Scholar

Bourdarie, S., Boscher, D., Beutier, T., Sauvaud, J.-A., and Blanc, M. (1997). Electron and proton radiation belt dynamic simulations during storm periods: A new asymmetric convection-diffusion model. J. Geophys. Res.102 (8), 17541–17552. doi:10.1029/97JA01305

CrossRef Full Text | Google Scholar

Brand, A., Allen, L., Altman, M., Hlava, M., and Scott, J. (2015). Beyond authorship: Attribution, contribution, collaboration, and credit. Learn. Pub.28, 151–155. doi:10.1087/20150211

CrossRef Full Text | Google Scholar

Drozdov, A. Y., Usanova, M. E., Hudson, M. K., Allison, H. J., and Shprits, Y. Y. (2020). The role of hiss, chorus, and EMIC waves in the modeling of the dynamics of the multi-MeV radiation belt electrons. J. Geophys. Res. Space Phys.125, e2020JA028282. doi:10.1029/2020JA028282

CrossRef Full Text | Google Scholar

Fälthammar, C.-G., and Mozer, F. S. (2007). [Comment on “Bringing space physics concepts into introductory electromagnetism”] on the concept of moving magnetic fields lines. Eos Trans. AGU88, 169–170. doi:10.1029/2007EO150002

CrossRef Full Text | Google Scholar

Fok, M.-C., Buzulukova, N. Y., Chen, S.-H., Glocer, A., Nagai, T., Valek, P., et al. (2014). The comprehensive inner magnetosphere-ionosphere model. J. Geophys. Res. Space Phys.119, 7522–7540. doi:10.1002/2014JA020239

CrossRef Full Text | Google Scholar

Fok, M. C., and Moore, T. E. (1997). Ring current modeling in a realistic magnetic field configuration. J. Geophys. Res.24 (14), 1775–1778. doi:10.1029/97GL01255

CrossRef Full Text | Google Scholar

Glauert, S. A., Horne, R. B., and Meredith, N. P. (2018). A 30-year simulation of the outer electron radiation belt. Space weather.16, 1498–1522. doi:10.1029/2018SW001981

CrossRef Full Text | Google Scholar

Glauert, S. A., Horne, R. B., and Meredith, N. P. (2014). Three-dimensional electron radiation belt simulations using the BAS Radiation Belt Model with new diffusion models for chorus, plasmaspheric hiss, and lightning-generated whistlers. J. Geophys. Res. Space Phys.119, 268–289. doi:10.1002/2013JA019281

CrossRef Full Text | Google Scholar

Hao, Y. X., Sun, Y. X., Roussos, E., Liu, Y., Kollmann, P., Yuan, C. J., et al. (2020). The formation of saturn’s and jupiter’s electron radiation belts by magnetospheric electric fields. Astrophysical J.905 (1), L10. doi:10.3847/2041-8213/abca3f

CrossRef Full Text | Google Scholar

Hao, Y. X., Zong, Q. G., Zhou, X. Z., Fu, S. Y., Rankin, R., Yuan, C. J., et al. (2016). Electron dropout echoes induced by interplanetary shock: Van Allen Probes observations. Geophys. Res. Lett.43, 5597–5605. doi:10.1002/2016GL069140

CrossRef Full Text | Google Scholar

Horne, R. B., Glauert, S. A., Kirsch, P., Heynderickx, D., Bingham, S., Thorn, P., et al. (2021). The satellite risk prediction and radiation forecast system (SaRIF). Space weather.19, e2021SW002823. doi:10.1029/2021SW002823

CrossRef Full Text | Google Scholar

Hudson, M., Jaynes, A., Kress, B. T., Li, Z., Patel, M., Shen, X.-C., et al. (2017). Simulated prompt acceleration of multi-MeV electrons by the 17March 2015 interplanetary shock. J. Geophys. Research:Space Phys.122, 10,036–10,046. doi:10.1002/2017JA024445

CrossRef Full Text | Google Scholar

Hudson, M. K., Elkington, S. R., Li, Z., and Patel, M. (2020). Drift echoes and flux oscillations: A signature of prompt and diffusive changes in the radiation belts. J. Atmos. Solar-Terrestrial Phys.207, 105332. doi:10.1016/j.jastp.2020.105332

CrossRef Full Text | Google Scholar

Jordanova, V. K., Kozyra, J. U., Nagy, A. F., and Khazanov, G. V. (1997). Kinetic model of the ring current-atmosphere interactions. J. Geophys. Res.102 (7), 14279–14291. doi:10.1029/96JA03699

CrossRef Full Text | Google Scholar

Jordanova, V. K., Morley, S. K., Engel, M. A., Godinez, H. C., Yakymenko, K., Henderson, M. G., et al. (2022). The RAM-SCB model and its applications to advance space weather forecasting. Adv. Space Res.doi:10.1016/j.asr.2022.08.077

CrossRef Full Text | Google Scholar

Lejosne, S., Allison, H. J., Blum, L. W., Drozdov, A. Y., Hartinger, M. D., Hudson, M. K., et al. (2022). Differentiating between the leading processes for electron radiation belt acceleration. Front. Astron. Space Sci.9, 896245. doi:10.3389/fspas.2022.896245

CrossRef Full Text | Google Scholar

Lejosne, S. (2014). An algorithm for approximating the L * invariant coordinate from the real-time tracing of one magnetic field line between mirror points. J. Geophys. Res. Space Phys.119, 6405–6416. doi:10.1002/2014JA020016

CrossRef Full Text | Google Scholar

Lejosne, S., and Kollmann, P. (2020). Radiation belt radial diffusion at Earth and beyond. Space Sci. Rev.216, 19. doi:10.1007/s11214-020-0642-6

CrossRef Full Text | Google Scholar

Lejosne, S., and Mozer, F. S. (2020). Experimental determination of the conditions associated with “zebra stripe” pattern generation in the earth's inner radiation belt and slot region. J. Geophys. Res. Space Phys.125, e2020JA027889. doi:10.1029/2020JA027889

CrossRef Full Text | Google Scholar

Li, X., Roth, I., Temerin, M., Wygant, J. R., Hudson, M. K., and Blake, J. B. (1993). Simulation of the prompt energization and transport of radiation belt particles during the March 24, 1991 SSC. Geophys. Res. Lett.20 (22), 2423–2426. doi:10.1029/93gl02701

CrossRef Full Text | Google Scholar

Lichtenberg, A. J., and Lieberman, M. A. (1992). “Regular and chaotic dynamics,” in Applied mathematical Sciences. 2nd (New York: Springer).

CrossRef Full Text | Google Scholar

Nénon, Q., Sicard, A., and Bourdarie, S. (2017). A new physical model of the electron radiation belts of Jupiter inside Europa's orbit. J. Geophys. Res. Space Phys.122, 5148–5167. doi:10.1002/2017JA023893

CrossRef Full Text | Google Scholar

Northrop, T. G., and Teller, E. (1960). Stability of the adiabatic motion of charged particles in the earth's field. Phys. Rev.117 (1), 215–225. doi:10.1103/PhysRev.117.215

CrossRef Full Text | Google Scholar

Patel, M., Li, Z., Hudson, M., Claudepierre, S., and Wygant, J. (2019). Simulation of prompt acceleration of radiation belt electrons during the 16 July 2017 storm. Geophys. Res. Lett.46, 7222–7229. doi:10.1029/2019GL083257

CrossRef Full Text | Google Scholar

Ripoll, J.-F., Claudepierre, S. G., Ukhorskiy, A. Y., Colpitts, C., Li, X., Fennell, J., et al. (2020). Particle dynamics in the earth's radiation belts: Review of current research and open questions. J. Geophys. Res. Space Phys.125, e2019JA026735. doi:10.1029/2019JA026735

CrossRef Full Text | Google Scholar

Roederer, J. G. (1970). Dynamics of geomagnetically trapped radiation. Heidelberg: Springer Berlin.

Google Scholar

Roederer, J. G. (1967). On the adiabatic motion of energetic particles in a model magnetosphere. J. Geophys. Res.72 (3), 981–992. doi:10.1029/JZ072i003p00981

CrossRef Full Text | Google Scholar

Roederer, J. G., and Zhang, H. (2014). “Dynamics of magnetically trapped particles,” in Foundations of the physics of radiation belts and space plasmas. Astrophysics and space science library (Berlin: Springer).

CrossRef Full Text | Google Scholar

Sandhu, J. K., Rae, I. J., Wygant, J. R., Breneman, A. W., Tian, S., Watt, C. E. J., et al. (2021). ULF wave driven radial diffusion during geomagnetic storms: A statistical analysis of van allen probes observations. J. Geophys. Res. Space Phys.126, e2020JA029024. doi:10.1029/2020JA029024

CrossRef Full Text | Google Scholar

Sauvaud, J.-A., Walt, M., Delcourt, D., Benoist, C., Penou, E., Chen, Y., et al. (2013). Inner radiation belt particle acceleration and energy structuring by drift resonance with ULF waves during geomagnetic storms. J. Geophys. Res. Space Phys.118, 1723–1736. doi:10.1002/jgra.50125

CrossRef Full Text | Google Scholar

Schulz, M., and Lanzerotti, L. J. (1974). Particle diffusion in the radiation belts. Heidelberg: Springer Berlin.

Google Scholar

Shprits, Y. Y., Kellerman, A. C., Drozdov, A. Y., Spence, H. E., Reeves, G. D., and Baker, D. N. (2015), Combined convective and diffusive simulations: VERB-4D comparison with 17 march 2013 van allen probes observations, Geophys. Res. Lett.42, 9600–9608. doi:10.1002/2015GL065230

CrossRef Full Text | Google Scholar

Stern, D. (1967). Geomagnetic euler potentials. J. Geophys. Res.72 (15), 3995–4005. doi:10.1029/jz072i015p03995

CrossRef Full Text | Google Scholar

Stern, D. P. (1970). Euler potentials. Am. J. Phys.38, 494–501. doi:10.1119/1.1976373

CrossRef Full Text | Google Scholar

Su, Z., Xiao, F., Zheng, H., and Wang, S. (2010). Steerb: A three-dimensional code for storm-time evolution of electron radiation belt. J. Geophys. Res.115, A09208. doi:10.1029/2009JA015210

CrossRef Full Text | Google Scholar

Subbotin, D. A., and Shprits, Y. Y. (2009). Three-dimensional modeling of the radiation belts using the Versatile Electron Radiation Belt (VERB) code. Space weather.7, S10001. doi:10.1029/2008SW000452

CrossRef Full Text | Google Scholar

Tu, W., Cunningham, G. S., Chen, Y., Henderson, M. G., Camporeale, E., and Reeves, G. D. (2013). Modeling radiation belt electron dynamics during GEM challenge intervals with the DREAM3D diffusion model. J. Geophys. Res. Space Phys.118, 6197–6211. doi:10.1002/jgra.50560

CrossRef Full Text | Google Scholar

Ukhorskiy, A. Y., and Sitnov, M. I. (2013). Dynamics of radiation belt particles. Space Sci. Rev.179, 545–578. doi:10.1007/s11214-012-9938-5

CrossRef Full Text | Google Scholar

Varotsou, A., Boscher, D., Bourdarie, S., Horne, R. B., Glauert, S. A., and Meredith, N. P. (2005). Simulation of the outer radiation belt electrons near geosynchronous orbit including both radial diffusion and resonant interaction with Whistler-mode chorus waves. Geophys. Res. Lett.32, L19106. doi:10.1029/2005GL023282

CrossRef Full Text | Google Scholar

Walt, M. (1994). Introduction to geomagnetically trapped radiation. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511524981

CrossRef Full Text | Google Scholar

Woodfield, E. E., Horne, R. B., Glauert, S. A., Menietti, J. D., Shprits, Y. Y., and Kurth, W. S. (2018). Formation of electron radiation belts at Saturn by Z-mode wave acceleration. Nat. Commun.9, 5062. doi:10.1038/s41467-018-07549-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Woodfield, E. E., Horne, R. B., Glauert, S. A., Menietti, J. D., and Shprits, Y. Y. (2014). The origin of Jupiter's outer radiation belt. J. Geophys. Res. Space Phys.119, 3490–3502. doi:10.1002/2014JA019891

CrossRef Full Text | Google Scholar

Yu, Y., Jordanova, V. K., Ridley, A. J., Albert, J. M., Horne, R. B., and Jeffery, C. A. (2016). A new ionospheric electron precipitation module coupled with RAM-SCB within the geospace general circulation model. J. Geophys. Res. Space Phys.121, 8554–8575. doi:10.1002/2016JA022585

CrossRef Full Text | Google Scholar

Zhao, H., Sarris, T. E., Li, X., Huckabee, I. G., Baker, D. N., Jaynes, A. J., et al. (2022). Statistics of multi-MeV electron drift-periodic flux oscillations using Van Allen Probes observations. Geophys. Res. Lett.49, e2022GL097995. doi:10.1029/2022GL097995

CrossRef Full Text | Google Scholar

Glossary

Appendix

Here, we detail how to obtain Eq. 24, leveraging the fact that the Euler potentials ($\alpha ,\beta$) are proportional to canonical variables (Eq. 17).

We assume some small variations in $\alpha$ and $\beta$ during $t$ and $t+∆t$. In which case, a Taylor approximation of the time variations of $\alpha$ and $\beta$, to the second order, yields

$\begin{array}{c}\left\{\begin{array}{c}\alpha \left(t+∆t\right)=\alpha \left(t\right)+\dot{\alpha }\left(t\right)∆t+\frac{\ddot{\alpha }\left(t\right)}{2}{∆t}^2\\ \beta \left(t+∆t\right)=\beta \left(t\right)+\dot{\beta }\left(t\right)∆t+\frac{\ddot{\beta }\left(t\right)}{2}{∆t}^2.\end{array}\right.\end{array}(\mathrm{A1})$

Rewriting $\dot{\alpha }$ and $\dot{\beta }$ in terms of Hamiltonian (Eq. 17), the second time derivatives are

$\begin{array}{c}\left\{\begin{array}{c}\ddot{\alpha }=-\frac{d}{dt}\left(\frac{\partial H}{\partial \beta }\right)=\frac{\partial }{\partial \alpha }{\left(\frac{\partial H}{\partial \beta }\right)}^2-\frac{\partial }{\partial \beta }\left(\frac{\partial H}{\partial \alpha }\frac{\partial H}{\partial \beta }+\frac{\partial H}{\partial t}\right)\\ \ddot{\beta }=\frac{d}{dt}\left(\frac{\partial H}{\partial \alpha }\right)=\frac{\partial }{\partial \beta }{\left(\frac{\partial H}{\partial \alpha }\right)}^2-\frac{\partial }{\partial \alpha }\left(\frac{\partial H}{\partial \alpha }\frac{\partial H}{\partial \beta }-\frac{\partial H}{\partial t}\right)\end{array}\right.\end{array}(\mathrm{A2})$

(see also Lichtenberg and Lieberman, 1992; their equation (5.4.10), p. 322).

Combining equations Eqs A1, A2, 17, we have

$\begin{array}{c}\left\{\begin{array}{c}∆\alpha =-\frac{\partial H}{\partial \beta }∆t+\frac{{\left(∆t\right)}^2}{2}\left(\frac{\partial }{\partial \alpha }{\left(\frac{\partial H}{\partial \beta }\right)}^2-\frac{\partial }{\partial \beta }\left(\frac{\partial H}{\partial \alpha }\frac{\partial H}{\partial \beta }+\frac{\partial H}{\partial t}\right)\right)\\ ∆\beta =\frac{\partial H}{\partial \alpha }∆t+\frac{{\left(∆t\right)}^2}{2}\left(\frac{\partial }{\partial \beta }{\left(\frac{\partial H}{\partial \alpha }\right)}^2-\frac{\partial }{\partial \alpha }\left(\frac{\partial H}{\partial \alpha }\frac{\partial H}{\partial \beta }-\frac{\partial H}{\partial t}\right)\right).\end{array}\right.\end{array}(\mathrm{A3})$

To the second order in $∆t$, we also have

$\begin{array}{c}\left\{\begin{array}{c}{\left(∆\alpha \right)}^2={\left(\frac{\partial H}{\partial \beta }\right)}^2{\left(∆t\right)}^2\\ {\left(∆\beta \right)}^2={\left(\frac{\partial H}{\partial \alpha }\right)}^2{\left(∆t\right)}^2\\ ∆\alpha ∆\beta =-\frac{\partial H}{\partial \alpha }\frac{\partial H}{\partial \beta }{\left(∆t\right)}^2.\end{array}\right.\end{array}(\mathrm{A4})$

Thus,

$\begin{array}{c}\left\{\begin{array}{c}\frac{1}{2}\frac{\partial }{\partial \alpha }{\left(∆\alpha \right)}^2=\frac{1}{2}{\frac{\partial }{\partial \alpha }\left(\frac{\partial H}{\partial \beta }\right)}^2{\left(∆t\right)}^2\\ \frac{1}{2}\frac{\partial }{\partial \beta }{\left(∆\beta \right)}^2=\frac{1}{2}\frac{\partial }{\partial \beta }{\left(\frac{\partial H}{\partial \alpha }\right)}^2{\left(∆t\right)}^2\\ \frac{1}{2}\frac{\partial }{\partial \alpha }\left(∆\beta ∆\alpha \right)=-\frac{1}{2}\frac{\partial }{\partial \alpha }\left(\frac{\partial H}{\partial \alpha }\frac{\partial H}{\partial \beta }\right){\left(∆t\right)}^2\\ \frac{1}{2}\frac{\partial }{\partial \beta }\left(∆\alpha ∆\beta \right)=-\frac{1}{2}\frac{\partial }{\partial \beta }\left(\frac{\partial H}{\partial \alpha }\frac{\partial H}{\partial \beta }\right){\left(∆t\right)}^2.\end{array}\right.\end{array}(\mathrm{A5})$

Combining Eqs A3–A5 in terms of expected values for the variations, we have

$\begin{array}{c}-\left[∆\alpha \right]+\frac{1}{2}\frac{\partial \left[{\left(∆\alpha \right)}^2\right]}{\partial \alpha }+\frac{1}{2}\frac{\partial \left[∆\alpha ∆\beta \right]}{\partial \beta }=\left[\frac{\partial }{\partial \beta }\left(H+\frac{∆t}{2}\frac{\partial H}{\partial t}\right)\right]∆t,\end{array}(\mathrm{A6})$

$\begin{array}{c}-\left[∆\beta \right]+\frac{1}{2}\frac{\partial \left[{\left(∆\beta \right)}^2\right]}{\partial \beta }+\frac{1}{2}\frac{\partial \left[∆\beta ∆\alpha \right]}{\partial \alpha }=-\left[\frac{\partial }{\partial \alpha }\left(H+\frac{∆t}{2}\frac{\partial H}{\partial t}\right)\right]∆t.\end{array}(\mathrm{A7})$

Assuming that the time interval, $∆t$, is very small in comparison with the characteristic time for the time variation of the Hamiltonian:

$\begin{array}{c}∆t\ll H/\left(\partial H/\partial t\right),\end{array}(\mathrm{A8})$

with $〈∆\alpha 〉=\left[∆\alpha \right]/∆t$ and $〈∆\beta 〉=\left[∆\beta \right]/∆t$, the rates of change of the expected values for the variations, we obtain

$\begin{array}{c}\left\{\begin{array}{c}-〈∆\alpha 〉+\frac{1}{2}\frac{\partial 〈{\left(∆\alpha \right)}^2〉}{\partial \alpha }+\frac{1}{2}\frac{\partial 〈∆\alpha ∆\beta 〉}{\partial \beta }=\left[\frac{\partial H}{\partial \beta }\right]=-\left[\dot{\alpha }\right]\\ -〈∆\beta 〉+\frac{1}{2}\frac{\partial 〈{\left(∆\beta \right)}^2〉}{\partial \beta }+\frac{1}{2}\frac{\partial 〈∆\beta ∆\alpha 〉}{\partial \alpha }=-\left[\frac{\partial H}{\partial \alpha }\right]=-\left[\dot{\beta }\right].\end{array}\right.\end{array}(\mathrm{A9})$