Wave-driven butterfly distribution of Van Allen belt relativistic electrons

Van Allen radiation belts consist of relativistic electrons trapped by Earth's magnetic field. Trapped electrons often drift azimuthally around Earth and display a butterfly pitch angle distribution of a minimum at 90° further out than geostationary orbit. This is usually attributed to drift shell splitting resulting from day–night asymmetry in Earth's magnetic field. However, direct observation of a butterfly distribution well inside of geostationary orbit and the origin of this phenomenon have not been provided so far. Here we report high-resolution observation that a unusual butterfly pitch angle distribution of relativistic electrons occurred within 5 Earth radii during the 28 June 2013 geomagnetic storm. Simulation results show that combined acceleration by chorus and magnetosonic waves can successfully explain the electron flux evolution both in the energy and butterfly pitch angle distribution. The current provides a great support for the mechanism of wave-driven butterfly distribution of relativistic electrons.

T he Earth's outer Van Allen radiation belt is composed of trapped electrons with relativistic energy (E k Z1 MeV). Fluxes and pitch angle distributions of relativistic electrons often exhibit dramatic and highly dynamic changes during geomagnetic storms or substorms, which are associated with different physical mechanisms. Those mechanisms include transport, energization and loss processes. Wave-particle interaction plays an important role in energy exchange between various modes of plasma waves and Van Allen radiation belt relativistic electrons. Two types of plasma waves, chorus and magnetosonic (MS) wave (also known as 'equatorial noise'), can efficiently accelerate electrons up to relativistic energies. Relativistic electrons can pose serious damage to satellites and astronauts in space, it is therefore important to understand and ultimately predict Van Allen radiation belt electron dynamic variations.
Radiation belt trapped electrons experience azimuthal driftmotion around Earth when they bounce between northern and southern magnetic mirror points. Due to the day-night azimuthal asymmetry in Earth's magnetic field, their drift shells which are traced out by their guiding centres can separate radially for different pitch angles, which is known as drift shell splitting 1 . Significant drift shell splitting occurs only in the outer radiation belt L46 (with L being the radial distance in Earth radii R E ), where the asymmetry becomes substantial. Previous theoretical 1,2 and observational [3][4][5] works have long confirmed drift shell splitting.
Drift shell splitting separates the high and low pitch angle particles in nightside injections as they move to the dayside magnetosphere. The higher pitch angle particles drift to larger radial distance beyond the magnetopause on the dayside and may consequently be lost from the distribution 6 . Therefore, drift shell splitting can be interpreted in terms of butterfly pitch angle distributions, namely, a sharp dropout in the flux of 90°electrons. Another drift shell splitting is in combination with a negative radial flux gradient as equatorially mirroring particles drift around Earth at locations L ¼ 6-12 (refs 7,8). Drift shell splitting will most easily affect trapped particle populations further out than geosynchronous orbit. Occurrence of drift shell splitting inside of geosynchronous orbit would need a very large magnetopause compression during very large geomagnetic storms.
A fundamental problem, both theoretically and observationally, is that there is any butterfly pitch angle distribution of relativistic electrons below L ¼ 5 associated with a new mechanism instead of drift shell splitting. Previous studies have analysed the characteristics and evolution of pitch angle distributions of the outer radiation belt electrons 9,10 . Chorus and MS waves were proposed to produce butterfly distributions by preferentially accelerating off-equator electrons 11,12 . Adiabatic processes could potentially yield the formation of butterfly distributions at locations LZ6 (refs 13,14). Direct confirmation of wave-driven butterfly distributions below L ¼ 5 requires simultaneous highresolution data but this is generally unavailable before the launch of NASA's Van Allen Radiation Belt Storm Probes in 2012 (ref. 15). Fortunately, the unique events on the 28 June 2013 geomagnetic storm observed from Van Allen Probes provide an excellent opportunity to identify such mechanisms.
Here we report the formation of a unusual butterfly pitch angle distribution of relativistic electrons around L ¼ 4.8, corresponding to the occurrence of strong chorus and MS waves at the same time. Using two-dimensional simulation, we show that flux enhancements of relativistic electrons are most pronounced between the medium pitch angles 30°and 60°due to the dominant acceleration process for combined chorus and MS waves. Meanwhile, the pitch angle distributions close to 90°are increased due to relatively larger acceleration process for chorus alone. The combined acceleration by chorus and MS waves substantially modifies the whole population of relativistic electrons, producing the butterfly distribution. Our detailed modelling, together with the correlated Van Allen Probes data, presents a further understanding on how chorus and MS waves play different roles in Earth's outer radiation belts.

Results
Correlated Van Allen Probe data. On 28 June 2013, a moderate storm with DstE À 100 nT (Fig. 1a) was triggered by an interplanetary coronal mass ejection 16 . This is a relatively minor magnetopause compression, thus drift shell splitting is not expected to occur inside of geosynchronous orbit. A large negative B z (the z component of the interplanetary magnetic field) occurs during the period 1100 hours on June 28 to 1200 UT on June 29 (Fig. 1b), leading to efficient coupling with Earth's magnetosphere and prolonged geomagnetic activity 17 . Distinct whistler-mode chorus and MS emissions (Fig. 1d) were observed by the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) Waves instrument 18,19 onboard Van Allen Probe A. Chorus and MS waves are right hand polarized electromagnetic waves. Chorus waves are excited by the injection of energetic (tens of keV) plasma sheet electrons into the inner magnetosphere during the period of enhanced plasma convection associated with the intervals of negative B z (refs 20-22). MS waves are generated by a ring distribution of low-energy (B10 keV) protons at frequencies close to the harmonics of the proton gyrofrequency 23 .
The fluxes of relativistic (2-3.6 MeV) electrons ( Fig. 1e-h) from the Relativistic Electron-Proton Telescope (REPT) instrument 24 onboard Van Allen Probes decreased during the main phase of the storm when the magnetopause was compressed down to L ¼ 7.6, and dramatically increased during the recovery phase after 29 June when the solar wind dynamic pressure decreased and the magnetopause moved beyond L ¼ 10 (Fig. 1c).
In Fig. 2, we show data during the period 000 to 1600 UT on 29 June 2013 (corresponding to the shaded area in Fig. 1) from the EMFISIS instrument for magnetic and electric spectral intensity of chorus and MS waves (Fig. 2a,b), wave normal angle y (Fig. 2c) and wave ellipticity (Fig. 2d). The observed chorus wave has a low normal angle (yE0°) and circular polarization (ellipticityE0) 25 . Meanwhile, the observed MS wave is highly oblique (yE90°) and linearly polarized (ellipticityE0) 26 . This indicates that the MS or chorus wave k vector is approximately perpendicular or parallel to the ambient magnetic field direction 12 .
A very interesting feature here is that distinct butterfly distributions of fluxes of relativistic (2-3.6 MeV) electrons from the REPT instrument occurred around L ¼ 4.8 over the interval 1228-1312 UT (Fig. 2e-h), corresponding to the occurrence of enhanced chorus and MS waves. Electron fluxes have minima around pitch angle 90°, and peaks around the pitch angle range of 30°-60°or 120°-150°. Observations of butterfly distributions at locations Lo5 are unlikely due to drift shell splitting as Earth's magnetic field is more dipolar and only distorted/stretched under extreme geomagnetic conditions. As this was a relatively minor compression and only a moderate geomagnetic storm, it is unlikely that drift shell splitting is the cause of these observations. Here we perform a two-dimensional simulation, together with the high-resolution observations from the Van Allen Probes, to examine whether chorus and MS waves can be responsible for evolution of both the energy and butterfly pitch angle distribution of the observed relativistic electron flux increase.
Numerical modelling. To model the temporal evolution of the electron distribution function, we need to solve the Fokker-Planck diffusion equation, which is associated with pitchangle and momentum diffusion coefficients driven by waveparticle interaction. Then evaluation of diffusion coefficients needs specific information of wave amplitudes and spectral properties. In this event, due to the brief UT intervals of Van Allen Probes, the in situ measurements of chorus or MS waves were confined to a limited range of magnetic local time (MLT) and L-shells: L ¼ 4-6 and MLT ¼ 22-24 for chorus; L ¼ 2-3.5 and MLT ¼ 10-18 for MS.
Previous works have demonstrated that MS waves can propagate over a broad region of MLT and L-shells 23,27,28 . Here using the previously developed programmes 29,30 , we present raytracing results of MS waves based on the wave data. MS waves are launched at L ¼ 5.6 for different parameters and MLT regions at the geomagnetic equator (Supplementary Table 1). The modelled results ( Supplementary Fig. 1) confirm that MS waves can propagate either into or out of the plasmasphere through the plasmapause, covering a broad region of L ¼ 2-5.6, particularly the observed butterfly distribution location L ¼ 4.8 and MLT ¼ 19.
In this event, Van Allen Probes stayed relatively deep inside the plasmasphere in the day-evening sector. Consequently, those MS waves which potentially occurred in the day-evening sector can not be directly observed by Van Allen Probes in their orbits. Since relativistic electrons drift eastward around Earth approximately in a circular orbit, the observed butterfly distribution should come from the whole contribution of resonances with those MS waves in different MLT regions.
We have also checked the Van Allen Probes data on 28 July 2013 and found that MS waves were indeed observed at L ¼ 4.8 and MLT ¼ 19. However, Van Allen Probes only passed the observed location in a brief UT interval every day, they may not observe the MS waves each time. In addition, using the Gaussian fitting method (not shown for brevity), we find that those MS waves display a similar Gaussian distribution to those at L ¼ 3 and MLT ¼ 17.3 on 29 July. It is therefore possibly reasonable to expect that the MS waves potentially existing at L ¼ 4.8 and MLT ¼ 19 on 29 July to follow the similar Gaussian distribution to that at L ¼ 3 and MLT ¼ 17.3. ARTICLE Meanwhile, for calculating the diffusion coefficients of MS waves, we assume that MS waves propagate very highly oblique (y ¼ 86°-89°), distributed in a standard Gaussian form (X ¼ tany) with a peak value X m ¼ tan89°. We noted that a concise modelling method by using maxima wave intensities was adopted in the previous work 26 .
Moreover, chorus occurs over a broad range of MLT from the nightside through dawn to the dayside 31,32 , potentially producing efficient scattering precipitation of energetic (tens of keV) electrons into the atmosphere. Considering that the ratio between the precipitated and trapped electron fluxes measured by Polar Orbiting Environmental Satellites (POES) is approximately proportional to the chorus power spectral intensity 33 , previous studies have developed a novel technique to obtain a dynamic global model of chorus wave amplitudes as a function of L, MLT and time by converting the measured POES flux ratios at different MLT 34,35 . Here we use the same approach (described in Methods), together with the previous parametric study 36 , to obtain the global distribution of chorus wave amplitude as a function of L, MLT in this event. Specifically, we remove the proton contamination by the same method in the previous work 37 . The model parameters listed in Table 1  In general, chorus waves which consist of substructures are coherent at the equator but become less coherent off the equator, resonating with relativistic electrons 38 . The consequences of coherent wave-particle interactions involving relativistic electrons and chorus waves off the equator have been presented in the previous work 39 . Here the quasi-linear theory to treat wave-particle interaction is adopted and the corresponding conditions will be discussed below. We assume chorus waves to obey the same Gaussian distributions in wave frequency and wave normal angle at and off the equator. Considering that relativistic electrons move along the geomagnetic field line and bounce back between mirror points, we consider wave-particle interaction at each location and calculate bounce-averaged diffusion coefficients.
The obtained bounce-averaged diffusion coefficients (Fig. 3) cover an entire region of pitch angles 0°-90°for chorus but a limited region 30°-60°for MS waves. This can allow for significant increases in fluxes of relativistic (2-3.6 MeV) electrons at medium pitch angles 30°-60°by MS waves 12 and at high pitch angles up to 90°by chorus waves 40 , leading to the resultant butterfly distribution on a timescale comparable to 10 h.
Using aforementioned bounce-averaged diffusion coefficients, we calculate phase space density (PSD) f t evolution of electrons in solving a two-dimensional Fokker-Planck diffusion equation 41 . The evolution of differential flux j is then simulated in the pitch angle region 0°-90°by the subsequent conversion j ¼ p 2 f t and the results are extended to the range 90°-180°due to the mirror symmetry. We show in Fig. 4a remarkable agreement between the simulated pitch angle distribution (Fig. 4c,d) and the REPT observation (Fig. 4a,b) during the acceleration interval. The most pronounced flux enhancement occurs over the medium pitch angles 30°-60°, where momentum diffusion rates for combined chorus and MS waves dominate (Fig. 3). Furthermore, the pitch angle distributions at higher pitch angles close to 90°are enhanced due to larger momentum diffusion rates of chorus waves. Finally, the combined acceleration by chorus and MS waves significantly alters the whole population of relativistic electrons, yielding the butterfly distribution in about 9 h.

Discussion
The present modelling provides a further confirmation of wavedriven butterfly pitch angle distribution of relativistic electrons observed by the REPT instrument at lower locations (Lo5) of the outer radiation belt. This is in contrast to the formation of butterfly distributions at higher locations (L46) induced by drift shell splitting due to local magnetic field asymmetry. The associated physical processes are schematically presented in Fig. 5. Relativistic electrons at Lo5 drift azimuthally around Earth inside the magnetosphere without loss to the magnetopause, continuously resonating with chorus and MS waves. The combined acceleration by chorus and MS waves occurs preferably between the medium pitch angles 30°and 60°, leading to the unusual butterfly distribution. The combined acceleration process described here is a universal physical process, which should also be effective in the magnetospheres of Jupiter, Saturn and other magnetized plasma environments in the cosmos.
It should be mentioned that, in a departure from the coherent chorus-electron interaction approach 38 , we use the quasi-linear theory of wave-particle interaction, which has been frequently adopted by radiation belt research community to treat waveparticle interaction 12,40,42 . Previous work 43 has shown that energetic electron pitch angle scattering by coherent chorus waves is 3 orders more rapid than by incoherent chorus waves. Coherent chorus waves can produce particle diffusion, phase trapping and/or phase bunching, which is determined primarily by the competing effects of wave amplitude and ambient magnetic field inhomogeneity at resonance. However, application of coherent wave-particle interaction to the stromtime global evolution of energetic particle distributions has not  Higher pitch angle electrons at L46 drift to larger radial distance beyond the dayside magnetopause and are consequently lost due to the day-night magnetic field asymmetry (2), yielding the regular butterfly distribution (3). (4) Electrons at Lo5 drift inside the magnetosphere without loss to the dayside magnetopause because of a very small day-night magnetic field asymmetry, continuously resonating with chorus and MS waves. (5) Wave acceleration can enhance electron PSD primarily within the medium pitch angles, leading to the formation of the unusual butterfly distribution (6). yet been established so far. We therefore expect the adopted quasi-linear theory here to be a basis for comparison with future developments in nonlinear modelling.
The basic conditions for the quasi-linear theory are that each individual particle randomly moves in velocity space, resonates with a succession of uncorrelated and small amplitude waves, and is scattered in a small amount in pitch angle and energy each time. Those conditions are basically satisfied in Earth's radiation belts for naturally generated MS and chorus waves, where the MS wave bandwidth is generally above the proton gyrofrequency up to the lower hybrid frequency and the chorus wave roughly lies in the frequency range 0.1-0.8 f ce (f ce being the electron gyrofrequency).
In general, there could be a distorted/stretched geomagnetic field on the nightside during a geomagnetic storm. We have examined the geomagnetic field data from Van Allen Probes and found that the observed data are close to the dipole field model during this storm period. The largest distortion of the dipole field is about 20% in the simulation period. Moreover, we perform a test-particle simulation of trajectories of trapped relativistic (2 MeV) electrons by the TS04 magnetic model 44 for different high pitch angles. Simulation results (Fig. 6) show that, starting at the location L ¼ 4.8 and MLT ¼ 24, relativistic electrons drift eastward around Earth and approach a farthest location LE5.1 on the dayside. Then they pass the observed location without loss to the magnetopause because the magnetopause locates above L ¼ 7.6 in the simulation period. Simulations for different energies or pitch angles (not shown for brevity) indeed show similar results. This indicates that drift shell splitting is unlikely to play a role in the observed butterfly electron distribution in this event.
A basic assumption is adopted here that the cold plasma density remains constant and unchanging during the 9-h simulation period. Using the upper hybrid resonance frequency observation from Van Allen Probes, we infer the ambient electron density at locations along Van Allen Probes' orbits and find that the electron density is comparable to the adopted plasma density and does not change much in 9 h. Analogous to the previous work 40 , this assumption is probably reasonable in the absence of simultaneous and continuous local time observations of cold plasma density in the 9-h period.
To check whether the plasmaspheric plume exists, we analyse the potential data from the Time History of Events and Macroscale Interactions During Substorms (THEMIS) spacecraft at L ¼ 4.8, and find no distinct presence of plasmaspheric plume except a slightly higher electron density in the region around/after MLT ¼ 16 in the simulation period. The plasmaspheric plume resulting from time-varying convection generally forms with a more structured distribution of plasma in the dusk region in the storm main phase. Meanwhile, the plume moves eastward (towards later MLTs) in Earth's rotation eastward, and gradually fades/erodes as the strength of convection decreases. As shown in Fig. 1, the simulation period is in the storm recovery phase on 29 June, 424 h after the onset of the storm on 28 June. Hence, the plume probably already faded and moved to the later MLT region. However, further research is still required in the future because there are numerous uncertainties in accurate determination of either the upper hybrid resonance frequency from Van Allen Probes or the THEMIS spacecraft potential.

Methods
Ray tracing of MS waves. We perform the ray tracing of MS waves with wave vector k and frequency o by the following standard Hamiltonian equations 45 : where R is the position vector of a point on the ray path, t is the group time, D represents the standard wave dispersion relation obeying D(R, o, k) ¼ 0 at every point along the ray path. The spatial variation in D can be written: where B 0 is the ambient magnetic field and N c is the background plasma density. We adopt a dipole magnetic field model and the MLT-dependent plasmatrough density model 46 . The Earth-centred Cartesian and local Cartesian coordinate systems for the ray-tracing calculation are described in Supplementary Note 1.
Calculation of the global chorus wave amplitude. Recent works 34,35 have constructed a global chorus wave model based on the data of low-altitude electron population collected by multiple POES satellites. The Medium Energy Proton and Electron Detector (MEPED) onboard POES has two electron solid-state detector (0°and 90°) telescopes to measure electron fluxes in three energy channels (4 30 keV,4 100 keV and4300 keV) 47 . The 90°telescope largely measures the trapped flux over the invariant latitude range between 55°and 70°, and the 0°t elescope measures precipitating flux inside the bounce loss cone at L41.4 (ref. 48).
Using the electron distribution function near the loss cone 33 , the chorus wave amplitude can be calculated from the ratio between the measured precipitated and trapped electron fluxes (30-100 keV and 100-300 keV). The electron energy spectrum is assumed to follow a kappa-type function 49,50 . The basic equation for linking the ratio of electron count rates measured by the 0°and 90°telescopes to chorus wave amplitude is shown in Supplementary Note 2. The ratios of precipitated and trapped electron fluxes measured by multiple POES satellites are used to construct the global chorus wave model over a broad range in L and MLT, and the ratios obtained in four distinct MLT sectors are shown in Supplementary  Fig. 2.
Calculation of diffusion coefficients due to chorus and MS waves. We assume that the wave spectral density B f 2 follows a typical Gaussian frequency distribution with a center f m , a half-width df, a band between f 1 and f 2 (ref. 51).
here B 2 t is the wave amplitude in units of Tesla and erf is the error function. To allow the data modelling to be more reliable, we average the observed wave magnetic field intensity over the indicated time period in this event and then apply the corresponding Gaussian fit for MS waves as shown in Supplementary Fig. 3. ARTICLE We choose the wave normal angle distribution to also satisfy a standard Gaussian form: where X ¼ tany (y 1 ryry 2 , X 1,2 ¼ tany 1,2 ), with a half-width X o and a peak X m . Based on the observation, we choose X m ¼ tan 89°, X o ¼ tan 86°, X 1 ¼ X m À X o , X 2 ¼ X m þ X o ; and the maximum latitude for the presence of MS wave l m ¼ 10°.
Based on the ray-tracing results ( Supplementary Fig. 1), we assume the MS wave spectral intensity at L ¼ 4.8 and MLT ¼ 19 hours to follow the similar Gaussian distribution to that as shown in Supplementary Fig. 3, which are then used to calculate the MS-driven bounce-averaged diffusion coefficients. For the chorus waves, we adopt the wave parameters as shown in Table 1 to calculate the chorusdriven bounce-averaged diffusion coefficients at the location L ¼ 4.8. We consider contribution from harmonic resonances up to n ¼ ± 5 for both chorus and MS waves. The ambient plasma density is obtained by the MLTdependent plasmatrough density model 46 and further assumed to remain latitudinally constant and unchanged during the 9-h simulation period. Note that there is no realistic data available either in 08-12 or in 16-20 MLT. We assume the chorus wave power in 08-12 MLT to follow the latitude-dependent model 36 . The corresponding adopted wave amplitude, wave spectrum and latitudinal occurrence of chorus waves, the ambient electron density and the equatorial ratio of plasma frequency to gyrofrequency (f pe /f ce ) are shown in Table 1.
Evaluation of the relativistic electron PSD evolution. The evolution of the electron PSD f t is calculated by solving the bounce-averaged pitch angle and momentum diffusion equation here p is the electron momentum, G ¼ p 2 T(a e )sina e cosa e with a e being the equatorial pitch angle, the normalized bounce time T(a e )E1.30 À 0.56sina e ; hD aa i, hD pp i, and hD ap i ¼ hD pa i denote bounce-averaged diffusion coefficients in pitch angle, momentum and cross pitch angle-momentum. The explicit expressions of those bounce-averaged diffusion coefficients can be found in the previous work. 41 The initial condition is modelled by a bi-modal kappa-type distribution function of energetic electrons 49 here n h is the number density of energetic electrons, l indicates the loss-cone index, G is the gamma function, k and y 2 k are the spectral index and effective thermal energy scaled by the electron rest mass energy m e c 2 (B0.5 MeV).
For the pitch angle boundary condition, f t ¼ 0 at the loss-cone a e ¼ a L (sina L ¼ L À 3/2 (4 À 3/L) À 1/4 ) to simulate a rapid precipitation of electrons inside the loss cone ( Fig. 2e-h), and qf t /qa e ¼ 0 at a e ¼ 90°. For the energy diffusion boundary conditions, f t ¼ f k 0 ð0:2 MeVÞ ¼ const at the lower boundary 0.2 MeV, and f t ¼ f k 0 ð10 MeVÞ ¼ const at the upper boundary 10 MeV. Based on the observation, we choose the following values of parameters: y 2 k ¼ 0:15 (B75 keV), k ¼ 4, n h ¼ 0.16 cm À 3 ; a 1 ¼ 0.6, l 1 ¼ 0.6; a 2 ¼ 0.4, l 2 ¼ 1.5. We solve the diffusion equation (6) using the recently developed hybrid difference method 41 , which is efficient, stable and easily parallel programmed. The numerical grid is set to be 91 Â 101 and uniform in pitch angle and natural logarithm of momentum. The drifting average is taken 25% for MS wave and 1/6 for chorus wave in each aforementioned MLT sector (Table 1).