Archimedean Spiral Distribution of Energetic Particles in Earth's Inner Radiation Belt

Archimedean spirals are common in various fields such as biology, engineering, astronomy, and space physics. Here we report the discovery of patterns resembling the general Archimedean spirals in particle distributions in the Earth radiation belt. Our analytic theory and numerical simulations demonstrate that electrons with initially asymmetric spatial distributions form such spirals in the inner magnetosphere, where particles at smaller radial distances move more slowly in angular velocity. These spirals result in time‐varying peaks and valleys in particle fluxes, referred to as “zebra stripes,” which are well consistent with Van Allen Probes measurements. Although the initial asymmetric distribution may be seeded by the electric field in the magnetosphere, the spiral formation does not require them. Furthermore, we show that, due to the same fundamental motion of charged particles in regions dominated by dipole fields, this spiral phenomenon may also appear in the proton distributions, as well as in planetary magnetospheres.


Introduction
The Archimedean spiral, named after the 3rd-century BC Greek mathematician Archimedes, is a spiral that results from the movements of points away from a fixed point at a constant speed along a line that rotates with a constant angular velocity.The general Archimedean form, defined as r = a + bϕ 1/c , has appeared in a variety of phenomena in nature, ranging from sunflower head (Vogel, 1979), Iridogorgia (Thomas, 1979;Tsuji & Müller, 2019) to galaxies arms (Bland-Hawthorn et al., 2019;Issa, 1991) and nebular patterns (Kim et al., 2017).It is also widely applied in antenna design (Kaiser, 1960), supramolecular polymerization (Sasaki et al., 2020) and fiber-based optical systems (Mangini et al., 2021), etc.Throughout the solar system, the average configuration of the background magnetic field also aligns with Archimedean spirals (Parker, 1958), as it is dragged by charged particles that are ejected radially from the surface of the rotating Sun.However, it remains unknown whether an Archimedean spiral exists in the near-Earth space, which is similarly filled with charged particles interacting with magnetic and electric fields.Here, we present compelling evidence for the existence of such a spiral resembling the general Archimedean form in the inner radiation belt through a combination of analytical theory, numerical simulations, and in-situ observations.The radiation belts, discovered by Explorer 1 satellite in 1958, are typically represented by two donut-shaped regions of energetic particles separated by a slot region.The outer belt is highly unstable due to multiple competing mechanisms such as acceleration/deceleration by various magnetospheric waves, magnetotail injections, atmospheric precipitation, and magnetopause loss (Li & Hudson, 2019;Ripoll et al., 2020;Shprits et al., 2006;Ukhorskiy et al., 2006).The inner belt, located from an altitude of 0.2-2 Earth radii (Re), was previously thought to be very stable.However, it has been challenged by several recent study (Hua et al., 2019(Hua et al., , 2020)), revealing its less stability than previously believed.In the dipole field, the motion of electrons can be described as its guiding center drift around Earth in the azimuthal direction and this drift motion is associated with the conservation of an adiabatic invariant (Northrop, 1963;Roederer, 1970).The evolution of the distribution function for those particles can be calculated by solving the advection equation (Ukhorskiy & Sitnov, 2013).In the following, we will present the framework for first deriving the analytic expression of particle distributions conserving the fourth adiabatic invariant (Schulz & Chen, 2008;Wolf, 1983) in Section 2 and then carrying out numerical simulations using the Rice Convection Model (RCM) (Toffoletto et al., 2003) in Section 3, which is based on the same drift physics.In Section 4, we will show comparison of the RCM simulation results with measurements from the Van Allen Probes.Section 5 will have a discussion and conclude the paper.

The Archimedean Spiral-Like Distribution for Particles Conserving the Fourth Adiabatic Invariant
Here we show how the general Archimedean spiral-like distribution forms theoretically.Assuming the plasma in a flux tube has isotropic pitch angle distribution, it follows that the velocity of the magnetic gradient and curvature drift of a particle with charge q and rest mass m 0 is given as (Wolf, 1983): where W k = m 0 v 2 /2 is the nonrelativistic kinetic energy, and B is the magnetic field strength written as B = B eo R E 3 / r 3 for the region inside 3 R E where a dipole filed is a very good approximation.Here B eo is the strength of the Earth's dipole field at zero altitude at the magnetic equator.
When a flux tube is filled with isotropic particles of a certain kinetic energy W k , it is shown that as the particles drift, their energy invariant (also called fourth adiabatic invariant) λ = W k V 2/3 is conserved (Schulz & Chen, 2008;Wolf, 1983), where V = ∫ds/B is the flux tube volume per unit magnetic flux extending over the entire magnetospheric portion of the field line.For the dipole field, the flux tube volume V can be written as (Wolf, 1983): where L = r/R E is the radial distance in R E in the equatorial plane.
The value of ω GC is negative for protons (q > 0) because of its westward drift and positive for electrons (q < 0).drift of particles in the inner region compared to the outer ones.Therefore, we anticipate that local maxima of the particles' distribution function for a given radial distance will arrange themselves into a spiral-like distribution due to the angular velocity difference in radial direction if the initial distribution function is azimuthally asymmetric and increasing outward.As exemplified in Figure 1, with an initial noon-midnight asymmetry in the distribution function (Figure 1a), electrons conserving a given λ will evolve into a spiral-like distribution (Figure 1b) due to the angular velocity difference in radial direction.
Suppose there are particles conserving a certain λ initially distributed radially at t 0 with an initial phase angle ϕ 0 .These electrons will drift azimuthally around the Earth driven by the corotation and convection electric fields and the gradient/curvature of the dipole magnetic field.Considering that the convection electric field is usually small in low latitude region, the azimuthal positions of the particles at time t will be: where ϕ cor = ω cor t is the phase of E × B drift due to the corotation electric field.ω cor is Earth's corotation angular velocity.ϕ 0 is a constant as well as ϕ cor at a fixed time t.Thus, we have: Equation 5 describes two general Archimedean spiral equation (whose standard equation is defined as Two notable features of the spiral pattern in the magnetosphere can be inferred from Equation 5. First, the parameter b increases with time t, indicating an increase in the number of spiral-shaped circles in the magnetosphere.Second, it suggests a decrease in the distance between any two consecutive circles as time processes.As indicated by the black dashed line in Figure 1b, the local maxima of the distribution function resemble a general Archimedean spiral at a given time t.Additionally, we will use a similar idea to provide an analytic theory for particles that drift in the equatorial plane while conserving the first adiabatic invariant (see Text S1 in Supporting Information S1 for the derivation of the Archimedean spiral-like distribution for particles conserving the first adiabatic invariant).

RCM Simulation Results
The RCM is an established first-principles model of the inner magnetosphere (Toffoletto et al., 2003).The model can calculate both E × B and gradient/curvature drift velocities v k = (E λ k ∇V 2/3 /q k ) × B/(B 2 ) for isotropically distributed λ-conserving particles (Wolf, 1983).The distribution function follows the advection equation , 1983), where S( f k ) and L( f k ) represent sources and sinks respectively and f k is the distribution function.The subscript k describes a given charge, mass and energy invariant λ k species.The electric field E and magnetic field B can either be specified or self-consistently computed (Toffoletto et al., 2003;Yang et al., 2019).The potential electric field is computed by solving the conservation of electric current coupling the magnetic-field-line-aligned current between the magnetosphere-ionosphere system.
First, we consider an idealized RCM simulation using a dipolar magnetic field and moving particles only by corotation electric field and gradient/curvature drift.The initial electron distribution function with an energy range from 40 keV to 1 MeV and L-shell from 1.2 to 6.0 is specified as f AE9 g(MLT), where f AE9 is the distribution function based on electron flux provided by the AE9 model (Ginet et al., 2014) with the input parameter of MLT = 12 hr and T = 15:21 UT on 14 December 2015 and g(MLT) = MLT/12 (MLT ≥ 12hr);2 MLT/12 (MLT < 12hr) is to impose the noon to midnight asymmetry to the original distribution function.This manipulation is reflected in the evenly increased number of electrons from noon to midnight, resulting in the final number of electrons at midnight being twice that of noon (Figure 2a), which indicates an initial phase angle ϕ 0 = π in Equation 5.There are 400 energy channels for electrons in this simulation, electrons in each of which conserves a different λ.The obtained spectrogram along MLT = 12 hr is shown in Figure 2e and a smooth background throughout the inner radiation belt.Besides, the electron precipitation is neglected.The most interesting result of the electrons' general Archimedean spiral-like distribution is the formation and evolution of "zebra stripes" (Ukhorskiy et al., 2014) as manifested in the spectra in the inner radiation belts.Taking the corotation electric fields and the satellite location (L sc ,ϕ sc ) into account.We substitute the energy invariant λ with the satellite-measured energy W sc = m 0 c 2 (γ 1) in Equation 3where γ = 1 (v/c) 2 ) 1/ 2 is the relativistic factor.The equation governing the formation of zebra stripes in the inner radiation belt can be presented as: Thus, λ-conserved particles whose energy W sc satisfies Equation 6 will be observed in the form of zebra stripes at (L sc ,ϕ sc ,t).Each integer n corresponds to particles of a certain energy contributing to an individual stripe.If we rewrite Equation 6as ), the energy difference between two adjacent stripes will be It will tend to zero as L sc or t approaches infinity, which indicates when and where the instrument's finite energy channels cannot distinguish the number of stripes observed.This time-decremented energy difference is also consistent with the narrower stripes when observed over consecutive passes (Lejosne & Roederer, 2016;Liu et al., 2016).
The RCM simulated spectra of logarithmic electron differential flux and detrended logarithmic electron differential flux (The detrended logarithmic electron differential flux is defined as log J log J, where log J denotes the nine-energy-channel sliding average of logJ (Liu et al., 2016).It highlights the energies of peak flux and valley flux in the spectra.)are shown in Figures 2e-2l separately with L-shell from 1.2 to 3 in a fixed radial direction of MLT = 12 hr (i.e., ϕ sc = 0).Well-organized zebra stripes gradually emerge, with their number increasing and the energy difference between each two adjacent stripes decreasing over time (see Movie S2).The governing Equation 6 predicts that a satellite at (L sc ,ϕ sc ) in the polar coordinate system at time t would observe a local flux maximum at energy W k , manifesting as zebra stripes.Thus, the peak-flux energies W peak where zebra stripes should emerge at L sc can be confirmed by W peak f (L sc ) = 3qB eo R E 2 π ϕ cor + 2nπ)/t,n = 0,1,2….For a given time t, the constant value on the right side implies that W peak f(L sc ) will remain constant as L sc varies, shown as black solid lines in Figures 2e-2l with their values labeled aside in unit of keV.Each of these black solid lines corresponds to a value of n representing electrons with an additional drift phase of 2π.The lines of constant W peak f (L sc ) clearly pinpoint where zebra stripes should emerge and are perfectly consistent with the simulation results.
As such, the idealized RCM simulation reproduces the formation and evolution of zebra stripes, which in turn validates our theory of general Archimedean spiral-like distribution.In the same vein, numerical simulations such as the Comprehensive Ring Current Model (Fok et al., 2001) may exhibit similar spirals in the electron distribution function and zebra stripes in the electron spectra as predicted by the analytic theory derived under the assumption that particles conserve the first adiabatic invariant while drifting in the equatorial plane (see Figure S1 in Supporting Information S1).

Comparison of the Simulation Results With Satellite Observations
Measurements from VAPs provide better observational support for our theory.Figures 3a-3d show the chronologically observed spectrograms of logarithmic electron differential flux from Radiation Belt Storm Probes Ion Composition Experiment (RBSPICE) onboard VAP-A and VAP-B during two perigee flybys.The 9-hr averaged geomagnetic activity D st index indicate it was a relatively weak magnetic storm during this period.Figures 3e-3h show the corresponding orbits in the GSM XY plane.During the former perigee flyby, the energetic electrons are smoothly distributed without any clear patterns observed in the inner radiation belt (see Figure S2 in Supporting Information S1).They gradually organized into regular patterns and became denser over time as shown in Figures 3a-3d.According to the observed zebra stripes, it can be inferred that the initial time t 0 is approximately 15:55 UT on 14 December 2015, and the initial phase angle ϕ 0 is close to the dawn side, with an estimated value of 1.7π (see Text S2 in Supporting Information S1 for the method to determine the initial drift time t 0 and phase angle ϕ 0 ).The black solid lines in Figures 3a-3d illustrate the estimated energies where zebra stripes are expected to emerge based on Equation 6.Although these calculated energies are slightly lower than the visually observed peak-flux energies in Figures 3a and 3b, they fit well with the color contours in Figures 3c and 3d.The differences between the predicted and the visually observed peak-flux energies are understandable.Nature is more intricate than our simple assumption, due to the non-dipolar nature of the realistic magnetic field and the complexity of the electric field's temporal and spatial variability over the almost 9-hr period from Figures 3a-3d.
This event is modeled with the RCM using a much more realistic setup.The simulation incorporates a selfconsistent electric and magnetic field.The high-resolution M-I coupling scheme in the RCM self-consistently calculates magnetic-field-aligned currents connecting magnetosphere to ionosphere.The RCM is two-way coupled (Yang et al., 2019) with a finite-volume regional MHD code with staggered non-uniform Cartesian grids to satisfy ∇•B = 0 (Silin et al., 2013).The RCM provides plasma pressure to the MHD code, and the MHD code provides the magnetic field to the RCM.We set up the initial conditions using the empirical magnetic field (Tsyganenko, 1989) and plasma model (Tsyganenko & Mukai, 2003) driven by the correpsonding solar wind and geomagnetic conditions.The start time of the simulation is t 0 , estimated using the same theory in previous paragraph.The initial electron distribution function setup is similar to the idealized one and is specified as f VAPB g (MLT).f VAPB is the distribution function based on electron flux measurements by VAP-B from 15:21 UT to 16:15 UT on 14 December 2015 (see Figure S2 in Supporting Information S1), during which the data-inferred initial drift time t 0 lies.g(MLT) is the same function as the idealized one.The 9-hr averaged geomagnetic activity D st index indicates relatively quiet conditions and without apparent stripe patterns observed throughout the inner radiation belt.The potential drop across the high latitude boundary of the RCM is estimated based on the geomagnetic index K P (Maynard & Chen, 1975).The electric precipitation loss rate is parametrized with K P and MLT, based on an empirical model (Orlova et al., 2014).
Figures 3i-3l are depicted following the same footprints and time records of VAP-A and VAP-B as Figures 3e-3h.Though the energy resolution based on the energy channels of VAPs is lower than our idealized simulation, we can reproduce discernible zebra stripes that increase over time in Figures 3i-3l.However, the numbers of zebra stripes in Figures 3i-3l are more than those in Figures 3a-3d, which could be attributed to the discrepancy between our self-consistent electric and magnetic fields and the real ones.In addition, the start time t 0 inferred from the equation based on the dipole magnetic field and the neglected convection electric field, may not accurately reflect the actual start time, which may cause our simulation's start time to be set ahead of the actual start time and further result in an increased number of zebra stripes.Figures 3m-3p show the general Archimedean spiral-like distribution of electrons that conserve the same energy invariant λ as Figures 2a-2d at the time record when VAP-A/VAP-B was located at L = 2 in Figures 3i-3l.Their locations are marked by pink dots.It is obvious that the number of turns increases and the distance between any two consecutive turns decreases over time.

Discussion and Conclusion
Several theories have been proposed to explain the appearance of zebra stripes.A test-particle model suggested that zebra stripes can be produced through resonance between global diurnal oscillations of magnetic and electric fields induced by rotation and electrons with drift period close to 24 hr in the inner radiation belt (Ukhorskiy et al., 2014).However, a single monochromatic uniform electric field is also found to be sufficient for reproducing the key characteristics of the stripes (Liu et al., 2016).Each peak/valley of zebra stripes corresponds to a population of electrons with similar gradient/curvature drift frequencies, and the stripes become narrower when observed over consecutive passes (Lejosne & Roederer, 2016;Liu et al., 2016).These properties are consistent with our theory and simulation results.A statistical analysis was carried out to show that zebra stripes were typically produced during substorm onsets when the penetration electric field was intensified in the magnetosphere (Lejosne & Mozer, 2020).Furthermore, similar patterns were also reported in the radiation belt observations of Jupiter and Saturn (Hao et al., 2020;Mauk et al., 2005;Müller et al., 2010;Sun et al., 2021Sun et al., , 2022)), explained by a global-scale convective electric field with an approximate noon-to-midnight orientation in Saturn Here the subplots follow the same format as Figure 2. The initial proton distribution function with an energy range from 100 keV to 1 MeV and L-shell from 1.2 to 6.0 is specified as f AP9 g(MLT), where f AP9 is the distribution function based on proton flux provided by the AP9 model (Ginet et al., 2014) with the input parameter of MLT = 12 hr and T = 15:21 UT on 14 December 2015 and g(MLT) is the same as the former one to impose a noon to midnight asymmetry to the original distribution function.There are more than 1,000 energy channels for protons in this simulation to ensure the energy resolution from 100 keV to 1 MeV.To optimize the viewing experience, the inner boundaries in (a-d) have been adjusted to commence at L = 2.The protons in (a-d) conserve a much larger λ = 10740 eV(R E /nT) 2/ 3 than the electrons in Figures 2a-2d.Thus, these protons have a larger gradient/curvature drift velocities than the electrons, resulting in an increase in the number of turns at the same time.(Hao et al., 2020;Sun et al., 2021Sun et al., , 2022) ) and dawn-to-dusk in Jupiter (Hao et al., 2020).Although the significance of the electric field is emphasized during the evolution of zebra stripes, our simulation shows the development of the general Archimedean spiral-like distribution and zebra stripes don't depend on the presence of an electric field.
Here the electric field in our simulation doesn't contribute to the acceleration of any particles as in the generation of drift echoes (Hudson et al., 2017;Vampola & Korth, 1992), but plays a role similar to that of the short-lasting pulsed electric field in zebra stripes generation (Liu et al., 2016).The low latitude electric field, such as dawn-todusk electric field or substorm-associated penetration electric field (Zhang et al., 2008(Zhang et al., , 2009) ) will enhance the E × B drift and result in radial displacement of particles, which may seed the initial asymmetric spatial distribution of plasma in the inner magnetosphere.Our study indicates that this preconditioning will facilitate the role of magnetic gradient/curvature drifts, which act as the key element for the generation the general Archimedean spiral-like distribution.
Furthermore, we predict that there are zebra stripes in other species due to the same drift physics.The RCM simulations show the general Archimedean spiral-like distribution and zebra stripes of λ-conserved proton in Figure 4.Although there are no zebra stripes reported in VAPs' ion spectrograms because of the turnoff of the TOF (time of flight) in RBSPICE measurement below 3R E and the ion measurements' contamination by electron measurements (Wang et al., 2018), zebra stripes had been reported in Saturn's proton and water group ion spectrograms within the energy range of 30 keV to ∼1 MeV (Sun et al., 2022).
In conclusion, this study sheds light on the fundamental drift physics underlying the formation of zebra stripes in the inner radiation belt and provides an elegant theory for predicting and understanding the behavior of energetic particles in the inner magnetosphere.It is well known that energetic particles conserving a given adiabatic invariant are subject to gradient/curvature drift.However, we have discovered for the first time that the differential angular velocity of such magnetic drifts with radial distance is crucial for initially asymmetrically distributed particles to rearrange to a general Archimedean spiral-like shape and for the formation of zebra stripes in the spectrogram.The electric field is responsible for creating the initial asymmetry but not directly act on the evolution process.The Rice Convection Model provides a good framework for investigating the particle dynamics in the inner magnetosphere and our simulation results further confirm the validity of our general Archimedean spiral-like distribution theory.This study has significant implications for understanding the complex dynamics of the radiation belts and their potential impact on space weather.

Figure 1 .
Figure 1.Cartoon of the formation of the spiral pattern in electron's distribution due to their magnetic drift.(a) The color contours represent the distribution function of λconserved electrons, in which deep blue is small and deep red is large.There is an obvious noon-midnight asymmetry in the initial electron distribution.(b) Electrons gradually evolve into a spiral-like distribution due to the angular velocity difference (ω 1 > ω 2 > ω 3 ) in the radial direction.

Figures
Figures 2a-2d shows the simulated 12-hr evolution of the distribution function of electrons for a given energy invariant in the magnetic equatorial plane, which illustrates the formation and development of a spiral shape in the local maxima of distribution function.The spiral becomes denser and clearer over time (see Movie S1).The theoretically derived spirals (black dashed lines in Figures2b-2d) resemble the general Archimedean form and are in perfect agreement with the simulation results at different times.

Figure 2 .
Figure 2. The 12-hr idealized RCM simulation results of the general Archimedean spiral-like distribution and zebra stripes for λ-conserved electrons.(a-d) The distribution function of electrons who conserve the energy invariant λ = 1120 eV(R E /nT) 2/ 3 in the RCM equatorial plane.The black dashed lines in (b-d) represent the derived spirals, similar to the general Archimedean form represented by Equation 5 with t 0 = 0, ϕ 0 = π, ϕ cor = ω cor t. (e-h) The energy-L spectrum of logarithmic electron differential flux along the radial direction at MLT = 12 hr (i-l) The detrended logarithmic electron flux.The black solid lines in (f-h) and (j-l) represent W peak f(L sc ) based on Equation 6 with their values labeled aside in unit of keV indicating where the zebra stripes should appear.

Figure 3 .
Figure 3. Chronological measurements of electron spectra from RBSPICE onboard VAPs and corresponding RCM simulation results.(a-d) Spectrograms of energetic electrons from measurements of VAP-A and VAP-B in chronological order during the weak geomagnetic storms.Two peak-flux points at a certain (L sc ,ϕ sc ,t) in (c) are extracted to calculate the initial drift time t 0 and phase angle ϕ 0 .The black solid lines represent the theoretical zebra stripes.(e-h) The VAP orbits in the GSM coordinate system.The thick red lines indicate the part of the orbit throughout the inner radiation belt corresponding to (a-d).(i-l) The RCM simulated spectrograms under selfconsistent electric field and magnetic field following the same footprints and time records as (e-h).(m-p) The general Archimedean spiral-like distribution of electrons for a given λ in the RCM equatorial plane at the time record when VAP-A/VAP-B was located at L = 2 in (i-l).Their locations are marked by pink dots.

Figure 4 .
Figure 4.The 3-hr idealized RCM simulation results of the general Archimedean spiral-like distribution and zebra stripes for λ-conserved protons.Here the subplots follow the same format as Figure2.The initial proton distribution function with an energy range from 100 keV to 1 MeV and L-shell from 1.2 to 6.0 is specified as f AP9 g(MLT), where f AP9 is the distribution function based on proton flux provided by the AP9 model(Ginet et al., 2014) with the input parameter of MLT = 12 hr and T = 15:21 UT on 14 December 2015 and g(MLT) is the same as the former one to impose a noon to midnight asymmetry to the original distribution function.There are more than 1,000 energy channels for protons in this simulation to ensure the energy resolution from 100 keV to 1 MeV.To optimize the viewing experience, the inner boundaries in (a-d) have been adjusted to commence at L = 2.The protons in (a-d) conserve a much larger λ = 10740 eV(R E /nT) 2/ 3 than the electrons in Figures2a-2d.Thus, these protons have a larger gradient/curvature drift velocities than the electrons, resulting in an increase in the number of turns at the same time.