Observations of drift phase structures in the radiation belts have multiplied over the last decade, facilitated by the use of instruments with high energy resolution channels (Krimigis et al., 2004; Sauvaud et al., 2006; 2013; Mitchell et al., 2013; Ukhorskiy et al., 2014; Hartinger et al., 2018). They have been the object of statistical analyses in the Earth’s inner belt (e.g., Lejosne and Mozer, 2020), in the outer belt (e.g., Zhao et al., 2022), and even at Saturn (e.g., Sun et al., 2021). Drift phase structures in radiation belt measurements are indicative of a transient magnetic local time (MLT) dependence in phase space density (PSD). As such, they are a direct challenge to the purely diffusive framework commonly utilized in radiation belt modeling and data analysis.

Indeed, the current picture for radiation belt acceleration (e.g., Jaynes et al., 2015) relies on the assumption that the effect of radial transport on radiation belt intensity is well captured by a diffusive equation (e.g., Lejosne et al., 2022a). That said, the validity of the radial diffusion equation relies on the assumption that radiation belts are fully phase mixed at all scales, including at the drift scale (e.g., Lejosne and Albert, 2023). As such, it excludes the possibility of any MLT dependence along a drift shell. One could expect any MLT-dependent PSD fluctuation to dissipate rapidly thanks to some efficient phase-mixing process (e.g., Schulz and Lanzerotti, 1974). Yet, ubiquitous observations of drift phase structures in the radiation belts suggest that: a) processes generating significant MLT-dependent structures in radiation belt fluxes (i.e., non-diffusive radial transport events) occur frequently and/or that b) the characteristic time for phase-mixing at the drift scale can be significant. Quantifying the latter is the object of this paper. Determining the amount of time it takes for a drift phase structure to dissipate and fully phase mix at the drift scale is important, because it provides information on the amount of time during which the radial diffusion equation cannot fully represent the effect of radial transport on radiation belt intensity once a MLT-dependent perturbation has occurred.

In this paper, we propose a thorough analysis of the characteristic time for drift phase mixing in the Earth’s radiation belts. While drift phase mixing is usually viewed as an effect related to the finite energy resolution of particle detectors at the drift scale (e.g., Schulz and Lanzerotti, 1974), we show that field perturbations also lead to natural phase mixing. The results of this analysis have theoretical and practical implications. From the theoretical standpoint, they contribute to clarifying the limits of the radial diffusion framework. From the practical standpoint, they provide analytical grounds to improve the analysis of measured drift phase structures and to specify resolution requirements for future particle detectors designed for radiation belt measurements.

Drift phase mixing at the drift scale corresponds to a process by which trapped particles with similar characteristics (adiabatic invariants, charge) initially located in a limited MLT sector along a drift shell end up covering all MLT sectors uniformly. We are interested in determining how long it takes for this phase homogenization process to occur. This defines the characteristic time for drift phase mixing, and the calculation of this time is the focus of this work. Within this section, we describe the framework by which fluctuations or uncertainties within the system drive phase mixing and the procedure by which the resulting phase mixing time can be extracted.

In the theoretical case of a perfectly resolved instrument (measuring one exact set of kinetic energy and pitch angle values), the characteristic times for phase mixing computed in Section 3, t E and t A , would both be infinitely long ( t E = t A = + ∞ ). Yet, we still expect the drift echoes to dissipate. This is because field fluctuations naturally present in space generate perturbations in drift frequency, and ultimately, phase mixing. We present a first quantification of these characteristic times for natural drift phase mixing in the following.

Although we have explored and characterized each mixing process individually, it is clear that in practice none of them operate in isolation. The “true” mixing time, therefore, is the cumulative effect of all operating in tandem. The focus of this section will be to use the analytic expressions for each individual process to construct an overall mixing time for the particles in MLT.

Physics-based radiation belt models usually consist of solving a three-dimensional Fokker-Planck equation reduced to a three-dimensional fully diffusive equation (e.g., Beutier and Boscher, 1995; Subbotin and Shprits, 2009; Su et al., 2010; Tu et al., 2013; Glauert et al., 2014).

While first implemented in the case of the Earth’s radiation belts, this theoretical framework has also been applied to the radiation belts of Jupiter and Saturn (e.g., Woodfield et al., 2014; 2018; Nénon et al., 2017). All these models rely on the assumption that radiation belts are fully phase mixed at all scales, including at the drift scale.

When the theoretical framework underlying these models was first put forward (e.g., Schulz and Lanzerotti, 1974), drift echoes had already been measured by the ATS-1 satellite at the geostationary orbit (e.g., Brown, 1968; Lanzerotti et al., 1971). Yet, it was assumed that rendering MLT-dependences in radiation belts could be omitted to first order. The radial diffusion framework was indeed designed to render radiation belt dynamics on long timescales (longer than the drift period). A new generation of particle instruments, equipped with very high-energy resolution channels, highlighted the omnipresence of drift phase structures in the radiation belts. This further challenged the validity of the radial diffusion approximation to model the effect of radial transport on radiation belt dynamics, a long-standing issue in radiation belt science which is currently the object of renewed scientific interest (e.g., Riley and Wolf, 1992; Ukhorskiy et al., 2006; Ukhorskiy and Sitnov, 2008; Degeling et al., 2008; Obrien et al., 2022; Osmane et al., 2023).

The objective of this work was to quantify the time it takes for a MLT-dependent structure to phase-mix at the drift scale, that is, the time during which the (purely diffusive) radial diffusion equation does not provide an accurate description of the effect of radial transport on radiation belt intensity once a MLT-dependent perturbation has occurred. Using relatively simple first assumptions to describe instrumental response and field perturbations, we found that the time it takes for particles initially localized in local time to phase-mix is measured in hours in the Earth’s radiation belts, determining timescales during which drift remains an important driver of the dynamics. Future studies could consist of describing the effects of radial transport on radiation belt intensity using a drift-diffusion model, following an approach similar to the one presented by Lejosne et al. (2023) for instance, to determine the role played by PSD’s radial gradients (e.g., Hartinger et al., 2020) in drift echoes’ lifetimes. Several analyses of electrostatic drift echoes present in the Earth’s inner radiation belt (also known as “zebra stripes”) showed that these structures can be observed up to 14 h after generation (e.g., Lejosne and Roederer, 2016; Liu et al., 2016; Lejosne and Mozer, 2020) by the Radiation Belt Storm Probes Ion Composition Experiment (RBSPICE) instruments (Gkioulidou et al., 2023) onboard the Van Allen Probes ( r ≈ 6 % ). This is consistent with the order of magnitude obtained in Section 4.3.1 (considering t F ≪ t E and t T ≈ t F .

Regarding theoretical studies, the parameter of interest should be independent of the instrument. Thus, it is the characteristic time for natural phase mixing associated with radial diffusion and energy diffusion, t o . In general, it takes hours to phase-mix in the Earth’s radiation belts, and this time is longer than the characteristic time for phase mixing due to the finite resolution of the particle instrument. This implies that non-diffusive radial transport events could have a greater and longer lasting effect on radiation belt dynamics than what was previously thought. Because the observed magnitude of drift echoes is dampened by the finite resolution of the instrument, this also means that care needs to be taken when interpreting drift echoes’ magnitude in terms of radial transport (e.g., Lejosne et al., 2022b; Obrien et al., 2022). Analysis of measured drift echoes’ magnitude should consider correcting for observational bias before interpreting drift periodic fluctuations in terms of radial transport.

The validity of the radial diffusion equation relies on the assumption that radiation belts are fully phase mixed at all scales, including at the drift scale, even though in reality, this is not necessarily the case. A more realistic implementation of the effects of radial transport on radiation belt intensity using a drift-diffusion model (e.g., Birmingham et al. (1967); Shprits et al., 2015; Lejosne and Albert, 2023) will contribute to determining how different from the radial diffusion paradigm the picture for large-scale radiation belt dynamics is when non-diffusive radial transport events are accounted for.

We assume that the drift phase locations of the trapped particles along the drift shell are distributed randomly, and described by a normal (Gaussian) distribution, N μ , σ 2 on the real axis. That is, the probability distribution function is e − 1 2 x − μ σ 2 / σ 2 π , with x ∈ R a random drift phase location. Because the MLT locations are 2π-periodic, the resulting probability distribution is actually a wrapped normal distribution. As a result, there exists a variance limit ( σ lim 2 ) for which the number of particles becomes uniform over all MLTs. To determine the value of the minimum variance, σ lim 2 , that defines a phase-mixed state, we compute the normalized number of particles at each hour MLT, h , (+/- 0.5 hr) assuming a constant average location, μ , and different values of the variance, σ 2 .

For illustrative purposes, we set μ = π (in radians), or, equivalently, noon in hour MLT. The quantity to compute is: n h ; σ 2 = 1 σ 2 π ∑ k = − ∞ + ∞ ∫ y = h − 0.5 π 12 h + 0.5 π 12 e − 1 2 y + 2 k π − π σ 2 d y (26)

In terms of error function, with: erf z = 2 π ∫ 0 z e − t 2 d t (27)

The distribution, n , is also: n h ; σ 2 = 1 2 ∑ k = − ∞ + ∞ erf h + 0.5 π 12 + 2 k π − π 2 σ − ⁡ erf h − 0.5 π 12 + 2 k π − π 2 σ (28)

Introducing the notations: ν k = h 12 + 2 k − 1 π 2 σ Δ ν = 1 24 π 2 σ (29) yields n h ; σ 2 = 1 2 ∑ k = − ∞ + ∞ erf ν k + Δ ν − ⁡ erf ν k − Δ ν (30)

When the standard deviation, σ , increases enough so that: Δ ν ≪ 1 , a Taylor expansion of Eq. 30 to first order in Δ ν provides the following approximation: n h ; σ 2 = 1 12 σ π 2 ∑ k = − ∞ + ∞ e − ν k 2 (31)

The contribution of the e − ν k 2 terms to the sum decreases rapidly as k increases, indicating that the distribution converges towards a limit. It is possible to approximate the expression of the distribution function even further by conserving only the contribution of k = 0 , k = 1 , and k = − 1 , in that case: n h ; σ 2 ≅ 1 12 σ π 2 e − ν − 1 2 + e − ν 0 2 + e − ν 1 2 (32)

For a phase-mixed distribution, we expect the resulting distribution, n , to be homogenous. In other words, we expect n = c o n s t a n t = 1 / 24 at all 24 MLT bins.

Figure 9A represents the normalized number of particles per MLT bin (with a size of one hour MLT), for three different values of the variance, σ 2 . It shows that the distribution, n , tends towards a phase-mixed state as the variance increases. Figure 9B represents the maximum relative distance between the normalized number, n , and the phase-mixed value over all 24 MLT bins, i.e, d = max n ∙ , σ 2 − 1 / 24 × 24 , as a function of the value of the variance, σ 2 . It shows that the convergence towards the phase mixed state is exponential. For the sake of simplicity, we choose π 2 radians as the minimum variance for which the resulting particle distribution is homogenous over all MLTs (d < 2 % for σ 2 < π 2 ). This choice is independent of the average location, μ .

In this section, we consider a population with kinetic energies randomly distributed around E , and described by a Gaussian distribution of standard deviation, σ E . We show that the variance of the drift phase locations is: σ 2 t = σ E 2 E 2 γ 2 + 1 γ γ + 1 2 Ω 2 t 2 (33)

The average phase location of the population at time, t, φ t , is: φ t = 1 σ E 2 π ∫ u = 0 t ∫ η Ω E + η e − η 2 2 σ E 2 d η d u (34) where the integral over η corresponds to the computation of the population’s average angular drift velocity. It accounts for the fact that the angular drift velocity depends on kinetic energy, and the population’s kinetic energies are randomly distributed around E . We describe the energy variable as E + η , where η indicates the random variation from the average energy, E . The probability distribution function for the random energy variable, η , is a Gaussian distribution, and it is included in the averaging. The integral over u from start time ( u = 0 ) to time t converts velocity in location.

A first order Taylor expansion of Ω E + η , leveraging Eq. 1, is: Ω E + η = Ω E + Ω E η E γ 2 + 1 γ γ + 1 (35) and Eq. 34 can be reformulated as: φ t = Ω E t + Ω E E σ E 2 π γ 2 + 1 γ γ + 1 ∫ t ∫ η η e − η 2 2 σ E 2 d η d u (36)

Because the integral in Eq. 36 is 0, the square of the average phase location is: φ t 2 = Ω 2 t 2 (37)

To compute the average of the square of the phase locations, we follow a similar approach. Given that: Ω 2 E + x = Ω 2 E + 2 Ω E η E γ 2 + 1 γ γ + 1 + Ω 2 E η 2 E 2 γ 2 + 1 γ γ + 1 2 (38) and 1 σ E 2 π ∫ η 2 e − η 2 2 σ E 2 d η = σ E 2 (39) φ 2 t = Ω 2 t 2 + σ E 2 E 2 γ 2 + 1 γ γ + 1 2 Ω 2 t 2 (40)

The result Eq. 33 is obtained by subtracting Eq. 37 from Eq. 40.

In this section, we show that the variance of the drift phase locations associated with radial diffusion, is, to first order: σ 2 t = 2 3 1 4 + 3 4 γ 2 2 D L L Ω 2 L 2 t 3 (41) for equatorial particles. For the sake of simplicity, let us first consider particles at high enough energy that Ω T ∼   Ω ≫ Ω E . In the case of equatorial particles drifting in a background magnetic field, Ω = − 3 M γ q R E 2 L 2 (42) where M is the first adiabatic invariant. Any small and slow enough (i.e., M conserving) field fluctuation generating a perturbation in radial location, δ L , leads to a perturbation in the angular drift velocity, δ Ω , due to the direct dependence of the angular drift velocity with radial location, L, and due to the dependence of the Lorentz factor with L , γ L . Given that: δ γ γ = M d B d L γ 2 E o δ L = − 3 M B E δ L γ 2 E o L 4 (43) δ Ω = − Ω 2 L 1 + 3 γ 2 δ L (44)

The variation in phase, Δ φ , for a cluster of particles starting from the same location, but that are scattered (including in the radial direction) with time by random field fluctuations, is, after a time interval, t : Δ φ t = ∫ 0 t δ Ω u du (45)

Given Eq. 44, this is also: Δ φ t = − Ω 2 L 1 + 3 γ 2 ∫ 0 t δ L u d u (46)

In a first approximation, δ L is a random process, which we view as a random walk. Thus, on average over many scenarios, the variation in phase is zero. For the variance, we focus on computing the average of the square of the phase variation, i.e., the statistical average of: Δ φ 2 t = Ω 2 4 L 2 1 + 3 γ 2 2 ∫ 0 t ∫ 0 t δ L u δ L v d u d v (47)

The statistical properties of random walks (e.g., Brockwell and Davis, 2002, p.14) are such that the autocovariance of δ L is: c o v δ L u , δ L v = 2 D L L ⁡ min u , v (48) where min u , v is the minimum of the variables u and v .

Thus: V a r Δ φ = Ω 2 D L L 2 L 2 1 + 3 γ 2 2 ∫ 0 t ∫ 0 t ⁡ min u , v d u d v (49)

Because ∫ 0 t ∫ 0 t ⁡ min u , v d u d v = t 3 / 3 , V a r Δ φ = Ω 2 D L L 6 L 2 1 + 3 γ 2 2 t 3 = 8 3 1 4 + 3 4 γ 2 2 Ω 2 D L L L 2 t 3 (50)

At this point in the derivation, it is important to realize that the population of particles for which we are quantifying the variance of the drift phase locations is scattering in radial distance as it is scattering in azimuthal location. To compute the characteristic time for phase mixing, it is preferable to stay on the same drift shell instead, meaning, it is preferable to focus on L constant. In that case, it is important to account for the fact that there is a correlation between variation in phase, Δ φ , and variation in radial location, Δ L , as quantified below.

Leveraging Eqs 46, 48, the covariance between Δ L and Δ φ is: c o v Δ L , Δ φ = − 2 D L L Ω L 1 4 + 3 4 γ 2 t 2 (51)

Given that the variance in radial location, V a r Δ L , is typically defined as: V a r Δ L = 2 D L L t (52) the correlation coefficient, ρ , between Δ L and Δ φ is (combining Eqs 50–52): ρ = c o v Δ L , Δ φ V a r Δ L V a r Δ φ = − 3 2 q q (53)

And the variance of the drift phase locations at L = c o n s t a n t , σ 2 , is related to the total variance, V a r Δ φ , through the correlation coefficient: σ 2 = 1 − ρ 2 V a r Δ φ = V a r Δ φ 4 (54)

This relationship can be obtained by considering the joint probability density function of a bivariate normal distribution (in Δ L and Δ φ ) (e.g., Hogg and Craig, 1970, p. 111 and sq) and noting that, for Δ L set to its average value ( Δ L = 0 , the distribution is still normal, with a variance equal to 1 − ρ 2 V a r Δ φ .The combination of Eqs 50, 54 yields Eq. 41.

Small perturbations of the fields of electrostatic origin do not yield first-order perturbations in Ω E . Thus, Eq. 41 remains valid even when considering Ω T = Ω + Ω E , in the electrostatic case.

Finally, let us mention that this derivation relies on the standard assumption that the variance in radial location grows linearly in time (e.g., eq. (C-12)). While this is a typical assumption in radial diffusion research (e.g., Ukhorskiy et al., 2005, their Figure 5; Lejosne and Kollmann, 2020, their equation (5.19)), super-diffusive and/or sub-diffusive regimes could also exist (e.g., Desai et al., 2021; Osmane and Lejosne, 2021). Assuming different statistical models for radial transport would of course lead to different results in this analysis.

In this section, we show that the variance of the drift phase locations associated with energy diffusion, is, to first order: σ 2 t = 1 6 γ 2 + 1 γ γ + 1 2 D E E E 2 Ω 2 t 3 (55) for equatorial particles. Any perturbation in energy, d E , drives a perturbation in drift frequency, δ Ω / 2 π , given by Eq. 35: δ Ω = γ 2 + 1 γ γ + 1 Ω E   d E (56)

The variation in phase, Δ φ , for a cluster of particles starting from the same location, but that are scattered in energy with time by random field fluctuations, is, after a time interval, t : Δ φ t = γ 2 + 1 γ γ + 1 Ω E   ∫ 0 t δ E u d u (57)

In a first approximation, δ E is a random process, which we view as a random walk. Thus, on average over many scenarios, the variation in phase is zero. For the variance, we focus on computing the average of the square of the phase variation, i.e., the statistical average of: Δ φ 2 t = Ω 2 E 2 γ 2 + 1 γ γ + 1 2 ∫ 0 t ∫ 0 t δ E u δ E v d u d v (58)

By definition of energy diffusion and the random walk, the autocovariance of δ E is c o v δ E u , δ E v = 2 D E E ⁡ min u , v . Because ∫ 0 t ∫ 0 t ⁡ min u , v d u d v = t 3 / 3 : V a r Δ φ = 2 3 D E E E 2 γ 2 + 1 γ γ + 1 2 Ω 2 t 3 (59)

At this point in the derivation, it is important to realize that the population of particles for which we are quantifying the variance of the drift phase locations is scattering in energy as it is scattering in azimuthal location. To compute the characteristic time for phase mixing, it is preferable to focus on a given energy, meaning, it is preferable to focus E constant. In that case, it is important to account for the fact that there is a correlation between variation in phase, Δ φ , and variation in energy, Δ E , as quantified below.

Leveraging Eqs 56, 57, the covariance between Δ E and Δ φ is: c o v Δ E , Δ φ = D E E E γ 2 + 1 γ γ + 1 Ω t 2 (60)

Given that the variance in energy, V a r Δ E , is typically defined as: V a r Δ E = 2 D E E t (61)

The correlation coefficient, ρ , between Δ E and Δ φ is, combining Eqs 59–61: ρ = c o v Δ E , Δ φ V a r Δ E V a r Δ φ = − 3 2 q q (62)

And the variance of the drift phase locations at E = c o n s t a n t , σ 2 , is related to the total variance, V a r Δ φ , through the correlation coefficient: σ 2 = 1 − ρ 2 V a r Δ φ = V a r Δ φ 4 (63)

The combination of Eqs 59, 63 yields Eq. 55.