Formation of electron radiation belts at Saturn by Z-mode wave acceleration

At Saturn electrons are trapped in the planet’s magnetic field and accelerated to relativistic energies to form the radiation belts, but how this dramatic increase in electron energy occurs is still unknown. Until now the mechanism of radial diffusion has been assumed but we show here that in-situ acceleration through wave particle interactions, which initial studies dismissed as ineffectual at Saturn, is in fact a vital part of the energetic particle dynamics there. We present evidence from numerical simulations based on Cassini spacecraft data that a particular plasma wave, known as Z-mode, accelerates electrons to MeV energies inside 4 RS (1 RS = 60,330 km) through a Doppler shifted cyclotron resonant interaction. Our results show that the Z-mode waves observed are not oblique as previously assumed and are much better accelerators than O-mode waves, resulting in an electron energy spectrum that closely approaches observed values without any transport effects included.

R adiation belts are formed when charged particles (usually electrons and protons) are trapped by large scale planetary magnetic fields. The particles then undergo significant acceleration up to relativistic energies [1][2][3] . Radial diffusion is usually assumed to be the dominant mechanism accelerating electrons in planetary radiation belts 4 . In recent years, however, acceleration by Doppler shifted cyclotron resonant wave particle interactions has been shown to be a key process in the Earth's radiation belts [5][6][7][8] and at Jupiter 9,10 . However, at Saturn, local wave particle interactions have been dismissed as unlikely to significantly accelerate electrons 11,12 .
Whilst radial diffusion can result in transport of charged particles towards the high magnetic field strength nearer the planet thus increasing their energy (assuming the first two adiabatic invariants are conserved), wave particle interactions are a local process acting in situ. Wave particle interactions transfer energy between circularly polarized plasma waves and charged particles gyrating around the planetary magnetic field lines. Waves can diffuse electrons not only in pitch angle (a change in the particle's velocity vector with respect to the local magnetic field) but also in energy resulting in particle acceleration 5 and an increase in the electron flux at higher energies. The acceleration process depends on the wave type, initial particle energy, pitch angle distribution and the plasma and magnetic field conditions 13 .
Previous attempts to assess electron acceleration by wave particle interactions at Saturn have assumed that whistler mode chorus waves would be most likely to produce acceleration as these waves are very effective at the Earth and Jupiter. However, at Saturn a combination of lower chorus wave power and higher plasma density means that electron acceleration by chorus waves is very weak 11 . We suggest that a different approach is required at Saturn where magnetic and plasma conditions are such that Zmode waves are intense inside 4 R S 14 . Initial studies of Z-mode interactions with electrons at the Earth show they are a very good candidate for electron acceleration 7,15 .
We find that Z-mode waves are very effective at accelerating electrons inside 4 R S and are capable of increasing the electron flux by four orders of magnitude in a year from an essentially empty radiation belt to a value comparable with observations. In comparison to radial diffusion in this region the Z-mode waves are more effective at filling in an empty radiation belt. We also show that Z-mode waves propagate much closer to the magnetic field than previously assumed 16 and that O-mode waves are much less effective at accelerating electrons than the Z-mode.

Results
Z-mode wave observations at Saturn. Saturn's magnetosphere is dominated by the presence of the Enceladus torus, a region of water group particles emanating from the moon Enceladus just inside 4 R S . The plasma density drops off very rapidly with latitude and also decreases notably away from the plasma source of Enceladus. Z-mode waves are frequently observed inside 4 R S where the combination of low plasma density and higher magnetic field strength resulting from proximity to the planet allows an abundance of Z-mode waves to propagate 14 . A case study of Z-mode waves at Saturn 16 assumed that the direction of wave propagation with respect to the magnetic field (known as the wave normal angle (WNA)) was highly oblique. However, Z-mode waves can also propagate along the magnetic field direction 17,18 . Figure 1a shows an example of wave activity from the Cassini radio and plasma wave instrument 19 (RPWS) as the spacecraft crossed from the northern to the southern hemisphere. The strong band of emissions near 5 kHz just before 10:00 occurred below the local plasma frequency f pe (white+) but above the L mode cut-off frequency f lcut (pink x) in a region where f pe < f ce (gyrofrequency). The apparent polarization 20,21 of these emissions was predominantly less than zero (Fig. 1b) indicating lefthand circular polarized waves in the plasma convention. Since only Z-mode waves can propagate with left-hand polarization between f lcut and f pe the emissions have been identified as Z-mode. The emissions at higher frequencies of 10−20 kHz are also left-hand polarized but since in this instance they lie in the range f pe < f < f ce they correspond to left-hand polarized O-mode waves. Z-mode and O-mode waves were also detected after the spacecraft crossed the ring plane in the southern hemisphere, where in this hemisphere an apparent polarization > 0 (red) corresponds to left-hand polarized waves. The 5 kHz wave band remains Z-mode despite changing polarization at~12:30 because f pe falls rapidly at this time such that the 5 kHz band is now in the regime f pe < f < f ce , where Z-mode waves are RH polarized. The intensity of the Z-mode waves in Fig. 1a is higher than the O-mode.
The Cassini RPWS instrument cannot directly observe the propagation angle of the waves with respect to the local magnetic field for waves with frequencies greater than 2.5 kHz but we can infer this from the polarization of the waves. A circularly polarized wave propagating parallel to the magnetic field with a small WNA will appear with a circular polarization, one propagating across the magnetic field or along the field with a highly oblique WNA will appear with a linear polarization. Fig. 1b shows that there is a mixture of strong and weak left-hand circular polarization with components of linear polarization (not shown) indicating that there is a distribution of WNAs which includes parallel propagation. Instability analysis shows that Z-mode waves can be generated in the field-aligned direction by an unstable distribution of electrons 17,18 . An analysis of the polarization data for this case shows that the WNAs are mostly close to the magnetic field with a peak at~22° (Fig. 2).
The effect of the waves on electron energy. To model the effect of Z-mode and O-mode waves on the flux of energetic electrons we calculated the diffusion coefficients for the waves using the pitch angle and energy diffusion of ions and electrons (PADIE) code 15 (described in the Methods section). We have calculated these coefficients for both Z-mode and O-mode waves since these waves often coexist with frequencies near 5 and 20 kHz. Figure 3 shows that electron pitch angle diffusion by Z-mode waves ( Fig. 3a) using the WNA distribution from Fig. 2 is much higher than that for O mode waves (Fig. 3c). This is also true for energy diffusion ( Fig. 3a and d). Pitch angle diffusion extends up to about 60° (Fig. 3a), indicating that Z-mode waves can scatter electrons into the loss cone at small pitch angles and cause electron loss into the atmosphere of Saturn. Note that the diffusion occurs in two energy bands, one extending up to about 1 MeV and another centered on 1-2 MeV. These bands correspond to Z-mode waves in two frequency bands near 20 and 5 kHz, respectively. Figure 3b shows that energy diffusion extends from a few tens of keV up to several MeV. Since electron phase space density falls rapidly with increasing energy, this suggests that electron energy diffusion to higher energies, i.e., acceleration, is very effective.
The local effect of Z-mode electron acceleration. The effect these diffusion coefficients have on the evolution of electron flux at different energies can be simulated using the British Antarctic Survey (BAS) Radiation Belt model 22 adapted for Saturn (this solves the Fokker-Planck equation as described in the Methods section). Figure 4 shows how the electron energy spectrum evolves over a period of 365 days for the Z-mode and O-mode waves. Figure 4a shows that Z-mode waves with a small WNA can accelerate electrons significantly from <0.1 to >5 MeV. The energy spectrum reaches the range of the empirical SATRAD model based on Pioneer 11 and Voyager 1 and 2 data after 365 days. It is remarkable that the spectrum can be reproduced over such a wide range of energies from wave acceleration alone with no transport effects included. In comparison, the effect of O-mode waves is very small (Fig. 4b).
The intensity of the waves is the main consideration for the strength of electron acceleration in this region because the plasma density is so low. The 20 kHz Z-mode waves are particularly important in accelerating electrons in the hundreds of keV range in these simulations where the initial condition is essentially an empty radiation belt with a low-energy seed population. To a good approximation in this case an increase or reduction in the overall Z-mode power will change the rate of acceleration by the same amount, e.g. increasing the wave power by a factor of 10 will produce the same level of flux in the simulation 10 times sooner. Cosmic Ray Albedo Neutron Decay (CRAND) will provide an additional source of electrons in the same energy range as the 20 kHz band resonates as will radial diffusion of electrons.
Z-mode wave acceleration and radial diffusion. The intensity of the Z-mode waves varies with radial distance from Saturn, exponentially increasing towards the planet 14 . Figure 5a, b show how the strength of the acceleration by Z-mode waves varies with radial distance. The changes in plasma density and magnetic field strength parameters counteract the increase in wave strength as we move closer to the planet 14 by shifting upwards the energy range over which electrons resonate with the Z-mode waves. Inside~2.7 R S the acceleration becomes slower, although it is still present all the way to the outer edge of the A-ring at 2.3 R s . The A-ring absorbs all energetic particles which is reflected in the zero increase in 1 MeV electrons over time at L = 2.3.
Radial diffusion is sufficiently slow that it would take some time longer to fill the region inside 2.7 R S by transporting electrons inwards as shown in Fig. 5c, d where the rate of filling inside 2.7 R S is not increased by adding radial diffusion. Since it is not clear how many electrons can be transported past the orbit of Mimas 23 or the other moons in this region we have used a worst case scenario and assumed that all electrons will be absorbed (in reality, at the very least, electrons drifting with the same period as the moon orbit will not be absorbed). Fig. 5e, f show a simulation of the effect of radial diffusion on its own over the same period as Fig. 5a-d. This shows that the Z-mode waves are more effective at filling up an empty radiation belt than just radial diffusion.

Discussion
Our results show that Z-mode waves play a major role in the formation of Saturn's radiation belts inside the orbit of Enceladus. We suggest the localized nature of this acceleration could explain the observed asymmetries in the high-energy electron population in this region 24 as asymmetries in Z-mode wave intensity are also observed 14 . Our results suggest that Z-mode waves may also play an important role in radiation belt dynamics at Jupiter, particularly inside the orbit of Io where the plasma density is low, and also at the other magnetized planets of the solar system.

Methods
Cassini data. The values for the gyro frequency have been determined from magnetic data taken by the Cassini MAG instrument. The plasma frequency is calculated from data taken by the Langmuir probe on the RPWS instrument on Cassini.

Radiation
g α ð Þ ¼ sin2α 1:3802 À 0:3198 sinα þ ffiffiffiffiffiffiffiffi ffi sinα p ð2Þ where f is the phase space density, t is time, α is the equatorial pitch angle, D αα , D EE , and D LL are the drift-averaged and bounce-averaged pitch angle, energy, and radial diffusion coefficients, respectively, E is the energy, L is the L-shell, µ and J are the first and second adiabatic invariants, E 0 is the electron rest energy, and τ is the loss timescale. The diffusion coefficients are calculated using the PADIE code 15 which solves the resonance condition for cold plasma dispersion in a magnetic field. We use a centered dipole magnetic field with an equatorial surface magnetic field strength, B 0 , of 2.1951 × 10 −5 T which is a good approximation this close to Saturn. We also use a model of the plasma density 25 including the density from the rings based on the Saturn Orbit Insertion data from Cassini.
The BAS RBM uses implicit methods to solve Eq. (1) on a two-grid system with one grid in α, E, and L and the other in µ, J, and L (L is the same in both grids). The pitch angle and energy diffusion are calculated as a separate step 26 using a time step of 5 s and then interpolated using cubic splines onto the µ, J, L grid which is solved every 500 s. The grid resolution is 90 × 30 × 30 points (α, E, L) and the energy grid uses equal spacing in the natural log of the energy (from 40 to 8000 keV at the maximum L value for the simulations in Fig. 4, from 40 to 3900 keV in Fig. 5 also at maximum L). We calculate the pitch angle at which all electrons are assumed to be lost to the atmosphere (the loss cone angle) at an altitude of 1000 km above the planet radius of 60,330 km.
The boundary conditions make use of the SATRAD 27 electron distribution model at Saturn. This is code that encapsulates the data from Pioneer 11, and Voyagers 1 and 2 and is freely available to download with documentation at http:// www.openchannelfoundation.org/projects/SATRAD/. We use the SATRAD differential flux value to define the minimum energy boundary. Although the SATRAD model is the best currently available for this region at Saturn there are significant uncertainties due to the small amount of data included. Comparing the SATRAD model to the spacecraft passes that it is derived from using figure A1 in the SATRAD report shows that there is at least an order of magnitude uncertainty in the fluxes from the model at energies of a few hundred to a few MeV. We therefore show a range of fluxes that span from SATRAD divided by 10 to SATRAD multiplied by 10.
The maximum energy boundary is set at a low value of phase space density, f = 10 50 kg −3 s 3 m −6 at the minimum pitch angle. The inner and outer L-shell boundaries required for the runs involving radial diffusion are subject to the condition f = 0 with the inner boundary set just outside each moon orbit and the outer boundary just inside (three separate runs are used to cover the L space from Epimetheus to Enceladus, separated by each moon). The maximum and minimum boundaries in α are set to ∂f/∂α = 0. The initial grid is set to be a straight line in log 10 (f) vs. log 10 (E) with a negative gradient from E min to E max . The precise gradient of the initial grid is unimportant 10 since the initial rise from a step function of f with E at E min is so rapid compared to the overall timescale. For the radial diffusion only run a seed population at maximum L (just inside each moons orbit) is added by including the flux spectrum from the SATRAD model. The initial dependence of f on α is set at (sin α) 2 to match the dependence in SATRAD.
The uncertainty in the SATRAD values also has an important impact on the initial and boundary conditions for the model runs with the flux at minimum energy acting as a source. We have therefore shown the range of possible fluxes that could be reached assuming that the flux and the minimum energy boundary varies within the SATRAD uncertainty range.
We have deliberately included a worst case scenario of moon absorption to demonstrate that the moons do not have a significant impact on the wave acceleration even when combined with radial diffusion. In fact the moons have a much less significant effect on the electrons in this region than they do on the protons 24 . Also some moons are more effective absorbers of energetic electrons than others 23 . To further isolate the effects of the Z-mode and O-mode waves we have not included CRAND as a source, or losses due to rings and other plasma waves.
Diffusion coefficients. The pitch angle diffusion coefficient, D αα , is the sum of the coefficients from the wave-particle interactions and the diffusion introduced by collisions in Saturn's atmosphere and the neutral torus produced by Enceladus. The PADIE 15 code is used to calculate D αα_wave with the information on the wave properties from the Z-mode wave survey 14 . Z-mode waves are assumed to exist in two wave bands, 5 and 20 kHz, with wave power assumed to be a Gaussian shape with frequency. In general the 20 kHz banded emission can be either O-mode, as in Fig. 1, or Z-mode as in the case presented by Gu and co-authors 16  The wave power 14 is also assumed to vary with both latitude, λ, and radial distance, R (Supplementary Figure 1). We have assumed that R is radial distance along the equator and therefore equivalent to L in a dipole field.
The latitude of peak power is λ m = 25.0°and the width is λ w = 14.142°. Note that wave power taken from Figure 7d in the paper by Menietti and co-authors 14 is the average wave power over latitude but the power required in the equation above is the peak power with latitude. We therefore divide the average wave power by the mean of a Gaussian with a peak power of 1 with the peak power latitude and width values given above, over the range of 0-50°latitude (which is the range over which data is available). This results in the factor of 1.75 × 10 −4 in the equation above. The WNA of the waves is also assumed to have a Gaussian shape based on the information in Fig. 2 with waves having a peak power at a WNA of 21.88°, a width of 17.18°and lower and upper cut-offs of 0°and 56.24°.
The O-mode at Saturn tends to appear in the same narrow band structures as the Z-mode. There are no published surveys of the O-mode at Saturn so for the best comparison of the acceleration capabilities of the two wave modes we replicate all of the Z-mode wave parameters for the O-mode.
For the PADIE calculations we use the cold plasma density model for the Enceladus torus 25 combined with the ring plasma density from Cassini Saturn Orbit Insertion data 25 . We use a scale height 28 to extend the equatorial densities to higher latitudes. We use a centered dipole magnetic field model as in the BAS 10 -11 10 -9 10 -7 10 -11 10 -9 10 -7 <D > (s -1 ) 0°F ig. 3 Bounce and drift averaged electron diffusion coefficients. a Pitch angle diffusion from Z-mode waves, 〈D αα 〉. b Energy diffusion from Z-mode waves, 〈D EE 〉/E 2 . c same as a but for O-mode waves. d Same as b but for O-Mode waves. The color scale indicates the strength of the diffusion. The labels above the horizontal axes show the maximum latitude (mirror latitude) an electron with the corresponding pitch angle would reach during its bounce motion RBM. One further consideration when calculating the diffusion coefficients for the Z-mode waves is the presence of a reverse resonance cone which exists where the wave frequency is lower than the upper hybrid frequency but higher than both the electron gyro frequency and plasma frequency.
We have included scattering from particles in the neutral torus created by Enceladus by calculating a value of pitch angle scattering for the presence of water ions. We use a model of water ions from Enceladus 29 to calculate D αα_torus 30 . We have assumed that the neutral torus is azimuthally symmetric around Saturn and consists of H 2 O, OH, and O with a defined scale height away from the equator 31 .
We have included the pitch angle scattering due to the atmosphere of Saturn using the same method 30 to calculate D αα_atmos . This gives a very high value for the total D αα at pitch angles up to the loss cone angle. To avoid cubic spline interpolation overshoot issues near the loss cone we use the value of D αα_atmos at α = 0°and then we use a quadratic curve to join this up to the value of D αα_atmos at the pitch angle grid point closest to the loss cone. The loss timescale within the loss cone is based on the time taken for the electron energy to decay 32 by 1/e. The ratio of D αα_atmos to the loss time is calculated at α = 0°and then the loss time is calculated at other points inside the loss cone using this ratio and the D αα_atmos at each grid point.
The radial diffusion coefficient D LL = 2 × 10 −14 L 7 s −1 is calculated from Cassini electron data 2 . We regard this as an upper limit for the radial diffusion as the value is likely to be contaminated by wave-particle interactions and to some extent the CRAND process.

Data availability
Cassini RPWS and MAG data are archived in calibrated, full resolution at the NASA Planetary Data System website: http://pds.nasa.gov. The Meudon Cassini RPWS/HFR polarization database is located at http://www.lesia.obspm.fr/kronos/ data. Other datasets generated during the current study are available from the corresponding author on reasonable request. showing the smoothing out of the locally accelerated electrons. e, f Radial diffusion only using the same initial condition as a-d but with an added SATRAD energy spectrum just inside the orbit of each moon to form the radial seed population. A worst case absorption has been assumed for the moons, where all electrons are absorbed, in reality the effects of the moons will be smaller