Formation of the mass density peak at the magnetospheric equator triggered by EMIC waves

We report a simultaneous observation of two band electromagnetic ion cyclotron (EMIC) waves and toroidal Alfvén waves by the Van Allen Probe mission. Through wave frequency analyses, the mass density ρ is found to be locally peaked at the magnetic equator. Perpendicular fluxes of ions (< 100 eV) increase simultaneously with the appearances of EMIC waves, indicating a heating of these ions by EMIC waves. In addition, the measured ion distributions also support the equatorial peak formation, which accords with the result of the frequency analyses. The formation of local mass density peaks at the equator should be due to enhancements of equatorial ion concentrations, which are triggered by EMIC waves’ perpendicular heating on low energy ions.


Introduction
Mass density controls the time rate of the change of magnetohydrodynamic (MHD) processes and is important for the space weather. It plays a pivotal role in the dynamic evolution of heavy ion composition in the magnetosphere, which can significantly alter the properties of plasma waves (Rauch and Roux, 1982;Meredith et al., 2003;Shi QQ et al., 2013;Min et al., 2015;Zong Q-G et al., 2017, Huang et al., 2020. The magnetospheric mass density is dependent on the number density and the ion components. However, it is difficult to directly measure ion components because the dominant populations contributing to the mass density often have low temperature, and the spacecraft's potential prevents the low-energy ions from reaching onboard detectors (Denton, 2006;André and Cully, 2012). In order to determine the ion components and the mass density, many researchers resort to characteristic frequencies of plasma waves, e.g. electromagnetic ion cyclotron waves Kim et al., 2015) and toroidal Alfvén waves (Singer et al., 1981;Schulz, 1996;Denton et al., 2004Denton et al., , 2006Denton et al., , 2015Denton, 2006;Takahashi et al., 2006Takahashi et al., , 2011Takahashi and Denton, 2007;Min et al., 2013;Nosé et al., 2015).
Variations of mass density ρ along the magnetic field lines are often modeled by a power law function: where ρ eq denotes the mass density at the magnetic equator, α decides the trend of ρ along the field line, and R E , R, and L denote the Earth radius, the geocentric radius, and the L shell value (in units of R E ), respectively. For a positive (negative) α, ρ will increase (decrease) toward higher latitudes, while for α = 0, ρ will be constant along the field line.
Previous researches have found α = 3 in the plasmasphere for a diffusive equilibrium model (Angerami and Thomas, 1964) and α = 4 in the plasmatrough for a collisionless model (Eviatar et al., 1964). Nosé et al. (2015) chose α = 0.5 to study the formation process of an oxygen torus in the inner magnetosphere at L = 3.0-4.0 and L = 3.7-4.5. Denton et al. (2006) found α = 2 appropriate for L = 4-5 and α = 1 appropriate for L = 5-6 through studying the average values of toroidal Alfvén frequencies. Min et al. (2013) chose α = 1 to derive the global equatorial mass density distribution covering radial distances from 4 to 9R E .
Nevertheless, observations also indicate the existence of α < 0, which means a local peak of ρ at the magnetic equator. Takahashi and Denton (2007) found that there was evidence for α < 0 in the afternoon magnetic local time (MLT) sector through statistically studying the toroidal Alfvén frequencies measured by the Geostationary Operational Environmental Satellites (GOES) at geostationary orbit. Denton et al. (2004) also found a negative α for the power law mass density model through investigating toroidal Alfvén waves observed at MLT ~ 16.6 hr and R ~ 6.4R E .
One possible reason for the peak of ρ eq is that the centrifugal force acts more effectively on the heavier ions and preferentially traps them at the magnetic equator (Lemaire and Gringauz, 1998;Denton et al., 2006). Denton et al. (2015) demonstrated that the more negative α at dusk potentially arises from the greater contribution of trapped ring current particles (with tens of keV temperature), especially O + , that drift westward from the magnetotail. In addition, perpendicular heating of cold ions by ion cyclotron waves may contribute to the formation of the local ρ eq peak (Takahashi and Denton, 2007). However, direct observational evidence has not corroborated the formation of equatorial mass density peaks from perpendicular heating of cold ions by ion cyclotron waves.
In this paper, we present a typical case where toroidal Alfvén waves and two band EMIC waves are simultaneously observed.
Wave frequency analyses illustrate that mass densities are locally enhanced at the equator. Through investigating measured ion fluxes, the formation of the mass density peak is shown to be promoted by the equatorial trapped ions, which are heated by EMIC waves.

Mass Density Derived from Toroidal Alfvén Waves
For the toroidal Alfvén mode waves in an arbitrary field model, the wave equation derived by Singer et al. (1981) is: where ξ is the field displacement in the azimuthal direction, and h is the scale factor for the normal separation between the field lines in the direction of oscillation which equals to Lcos 3 λ in the dipole field, where λ is the magnetic latitude. s is the distance along the field line with two ends fixed at ionosphere altitude (100 km in this work), ω is the frequency of toroidal mode waves, B 0 is the unperturbed magnetic field, and V A is the Alfvén velocity equal to B 0 (μ 0 ρ) -1/2 , where μ 0 is the free space permeability. With a shooting method (Press et al., 1997), Equation (2) is numerically solved for the mass density using a fourth-order Runge-Kutta algorithm.

Mass Density Derived from Cutoff Frequencies
The dispersion relationship for a parallel propagating EMIC wave (L mode) in a multi-ion species cold plasma (H + , He + , O + ) is : where n is the index of refraction, ω is the wave frequency, ω pe is the electron plasma frequency, Ω e is the electron gyrofrequency, ω pj is the ion plasma frequency, Ω j is the ion gyrofrequency, and the suffix j denotes ion species (H + , He + , O + ). After some algebra and approximations, the cutoff frequencies, where n 2 = 0, as functions of ion concentrations (η j ) are as follows : ω He + ,co = Ω H + 32 where η H + , η He + , and η O + are the concentrations of H + , He + , and O + , respectively. ω H + ,co denotes the H band cutoff frequency, while ω He + ,co denotes the He band cutoff frequency. Ω H + is the H + gyrofrequency. Equation (7) maintains the electrical neutrality (electron number density n e equals ion number density n i ). The ion concentrations can easily be obtained by solving the Equations (4), (5), and (7). Then the local mass density ρ 0 is: The equatorial mass density ρ eq is obtained by bringing back ρ 0 to Equation (1).

Observations and Analyses
Toroidal Alfvén waves and EMIC waves were simultaneously observed by the Van Allen Probe-B between 05:30:00 UT and 08:30:00 UT on Aug. 29, 2017. The observations consist of combined data from Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) (Kletzing et al., 2013), Electric Field and Waves Suite (EFW) (Wygant et al., 2013), and Helium, Oxygen, Proton, and Electron (HOPE) mass spectrometer (Funsten et al., 2013) instruments on board the Van Allen Probe-B.

Two Band EMIC Waves
The electron number density (n e ), EMIC wave power spectral density (PSD), ellipticity (ε), and wave normal angle ( Figures 1c and 1d show that the EMIC waves have mixed left-hand and linear polarizations with small normal angles. In order to determine the ion components accurately, we choose data from 07:45:41 UT to 07:50:41 UT corresponding to the time interval (300 s) between two vertical dashed lines in Figures 1b−1d to perform FFT. The time-averaged PSD is used for cutoff frequency determination (results with input spectra from different times are also checked; the resultant average ion mass has a minor change). The corresponding PSD is displayed in Figure 1e. The cutoff frequencies are determined through the following: (1) PSD(f+df)-PSD(f) > 0.5×10 −3 nT 2 /Hz.

Magnetoseismic Analyses
We smooth the 1 Hz resolution magnetic fields with a 90 s sliding window to obtain the background magnetic fields. Then we subtract the background from the magnetic fields to get the magnet-ic fluctuations. The magnetic fluctuations in GSE coordinates are transformed to magnetic field aligned (MFA) coordinate system with azimuthal component B ϕ , radial component B r and compressional component B z . Only the electric field components E y and E z in the MGSE coordinates obtained from the EFW are used here because the component E x along the spin axis is less accurate. We obtain the three dimensional electric field (E ϕ , E r , and E z ) in MFA coordinate system under the assumption of E·B = 0. This method works well because the angle between the magnetic field line and the spacecraft spin plane is larger than 15° in this case. The PSD of ULF wave is obtained by applying a 30 min FFT sliding window and a 5 min moving forward step to fields in MFA coordinates.     (Takahashi et al., 2011), the harmonic at ~10 mHz is an even mode. In addition, for the power law density model, the spacing between the successive harmonics, f n+1 -f n , is approximately identical (Schulz, 1996;Takahashi and Denton, 2007), and one more visible harmonic exists below this even mode (see Figure 2a, 05:00 UT to 06:30 UT, ~6 mHz). Hence, the emission at ~10 mHz is considered to be the second harmonic.   (c)). (f) Wave power spectrum for fields between 07:33:11 UT and 08:03:11 UT (30 min, between two red dashed lines in (a)) estimated by FFT (black curve) and MEM (red curve). We apply FFT to the fields between 07:33:11 UT and 08:03:11 UT (30 min, denoted by two red vertical dashed lines in Figure 2a, the same central time as Figure 1e) when EMIC waves coexisted to determine harmonic frequencies, which is denoted by the black line in Figure 2f. The red solid curve is the estimated PSD using a maximum entropy method (MEM) similar to Takahashi and Denton (2007). We adjust the order of MEM approximation to find the best representative spectral peaks of f 2 and f 3 . The harmonic f 2 ~9.44 mHz and f 3 ~16.67 mHz are selected for magnetoseismic analyses due to their stronger power and distinct harmonic structure.

Earth and Planetary
A dipole magnetic field model is adopted because, during the interval, there is small difference of the field line calculated by the dipole model and the TS05 model (Tsyganenko and Sitnov, 2005), and the magnetic field derived from the dipole model (~170 nT) is close to the observation (~172 nT). Assuming the field line is fixed at an altitude of 100 km, we numerically solve Equation (2) for theoretical harmonic frequencies with ρ eq determined in Section 3.1 through altering α from −6 to 6. The best parameter α best is determined from the best fit between theoretical and observed frequencies (f n,th and f n,obs ) by minimizing the quantity, where N is the number of harmonics used and subscript n is harmonic labels. For a spectrum with resolution of 0.56 mHz (signal length of 30 min), the error in manually harmonic determination is ~0.05. We find α = −5.0 ~ −2.9 (S < 0.05) with α best = −3.8 (minimum S < 0.04). The negative α indicates that the mass density reaches a local peak at the equator and reduces toward larger latitudes. To note, in a real sense the mass density will finally increase with latitude as the magnetic field line approaches the ionosphere. The value of α < 0 exists because the magnetic field line near the equator plays a main role in the Alfvén frequency determination . Based on the WKB formula for the Alfvén frequency , at higher latitudes the magnetic field is stronger and hence V A is larger, placing less effect on the values of the frequencies.f For a dipole magnetic field with a power law distribution, the separation between the successive harmonics does not depend significantly on the harmonic number (Schulz, 1996;Takahashi and Denton, 2007). The relationship between the parameter α and the harmonic frequency separation is introduced by Takahashi and Denton (2007). The normalized frequency separation decreases from ~0.40 for α = −6 to ~0.33 for α = 6, while for α = 0 the separation is ~0.38. Hence, a smaller α corresponds to a broader normalized frequency separation. We apply a 30 min sliding FFT window and a 15 min advancing step to B ϕ to determine the harmonic frequency for frequency separation checking. For every spectral line, peaks with spectral width at 75% of the peak intensity (> 0.1 nT 2 /Hz) less than 3 mHz are selected. Then, we eliminate the outliers and manually choose the most representative frequencies. Figure 3a shows the PSD of B ϕ with harmonic frequencies (black dots) over-plotted. EMIC wave amplitudes are shown in Figure 3b (H band and He band are denoted by red and blue lines, respectively). The normalized frequencyff f separations with error bars are shown in Figure 3c. Obviously, is larger (~0.4-~0.45) when EMIC waves show up (between two red dashed lines), corresponding to smaller α (< 0) than most of those outside of this time interval ( < 0.38 and α > 0). This agrees with the result of the theoretical calculation (Takahashi and Denton, 2007). The presence of EMIC waves and broader harmonic frequencies indicates that the equatorial mass density peaks are possibly triggered by EMIC waves.

Ion Measurements
We divide the case into Event 1 (time before 07:00:00 UT) and Event 2 (time after 07:00:00 UT) according to the EMIC waves' intensities. Figures 4a and 4b show the electron number densities (n e ) derived from upper hybrid frequencies of Event 1 (left) and Event 2 (right), respectively. Figures 4c and 4d show the ion number densities for H + (n H , red lines), He + (n He , magenta lines), O + (n O , cyan lines), and ensemble (n total = n H + n He + n O , black lines) obtained from HOPE flux data with energies < 100 eV.  Figures 4 (i, k, m) and (j, l, n), respectively. All fluxes are summed over the energy interval titled above each panel. The proportion of higher energy (> 100 eV) ion fluxes are very small, hence they are not included. As EMIC waves show up, H + fluxes are significantly enhanced and become more concentrated near pitch angle at 90° (Figures 4c and 4h), which indicates a trapping distribution of these low energy H + . He + fluxes show a weaker relation with the waves in Figures 4k and 4l, and weaker anisotropic distribution is found. No simultaneously enhanced flux of O + is observed, as shown in Figures 4m and 4n. Hence, as in Figures 4c and 4d, n H , n He , and n total increase with enhanced wave amplitudes, while n O remains small and almost constant. In addition, the simultaneous increases of n e and n total are well-correlated with wave emissions and trapping H + fluxes, indicating certain modulation mechanisms among them. The coexistence of enhanced wave amplitudes and trapped H + fluxes indicate that low energy H + ions should be perpendicularly heated by the EMIC waves. Heating of low energy ions by EMIC waves has been observed and simulated by many previous researchers (Gary et al., 1995;Omidi et al., 2010;Bortnik et al., 2010;Kitamura et al., 2018;Yu XD and Yuan ZG, 2019). Cold protons have a non-resonant interaction with enhanced proton-cyclotron-like fluctuations generated near the magnetic equator (Gary et al., 1995). The non-resonant character of this interaction has been confirmed by Anderson and Fuselier (1994), and their observations show that the low-energy protons' perpendicular temperature increase is always observed at 90°p itch angles. The simulation results of Yu XD and Yuan ZG (2019) also show that, under the influence of the EMIC wave, the perpendicular beta (temperature) of the cold protons increases significantly, while the parallel beta (temperature) remains almost unchanged. Thus, the perpendicular heating of the cold protons can be explained by the non-resonance interaction with the EMIC waves. Ions with larger equatorial pitch angles will mirror at lower latitudes and may lead to higher ion abundances near the equator. Therefore, by investigating the measured ion fluxes, we find that the peak formation is most likely due to the enhanced ion abundances at equator, which are triggered by EMIC waves' perpendicular heating on low energy ions. Figures 4e and 4f illustrate that ρ (black lines) increases while M (red lines) decreases when EMC waves exist. Since ρ is the product of M and number density n, the formation of a local mass density peak can be attributed to a dominant increase of n at the equator.

Discussion
According to Liouville's theorem, for a bi-Maxwellian distribution with A > 1 (A < 1) the integral number density will decrease (increase) towards higher latitudes along the field line, while for A = 1 the number density will be constant (Roederer and Zhang H, 2014). A is the temperature anisotropy (perpendicular temperature divided by parallel temperature). For quantitative analyses, we convert the time averaged H + fluxes to phase space densities in two time intervals as denoted by two black vertical lines in Figure 4c (06:08:24 UT to 06:13:24 UT of Event 1) and Figure 4h (07:44:12 UT to 07:47:12 UT of Event 2), respectively. Then, we adopt the bi-Maxwellian distribution to fit the ion phase space densities (Yuan ZG et al., 2019) with energy < ~100 eV and map the number densities along the field line according to Liouville's theorem. Unfortunately, the measured O + fluxes are too deficient to perform fitting (as shown in Figures 4m and 4n). Because no directional preference is observed in low energy O + fluxes, we assume the observed O + distribution has A = 1, n O+ = 0.0649 cm −3 of Event 1 and A = 1, n O+ = 0.0388 cm −3 of Event 2 (densities are derived from HOPE by adding all fluxes with energy level < 100 eV). The fitting parameters for H + , He + , and O + are displayed in Table 1. Liouville's theorem states that the distribution function f(r, v, t) Earth and Planetary Physics doi: 10.26464/epp2021008 37 along a single particle trajectory through any point (r(t), v(t)) in phase space is constant, which means that f(r(t), v(t), t) = f(r 0 , v 0 , 0) with initial condition r(0) = r 0 and v(0) = v 0 . For simplicity, we assume fitted ions are trapped on dipole field lines under adiabatic conditions. We still uphold some basic restrictions, i.e. absence of field-aligned electric fields, no time-variations, and symmetric angular distributions everywhere (Roederer and Zhang H, 2014). We choose the distribution function f(T 0 , μ 0 ) at the magnetic equator and f(T s , μ s ) at point s on the same field line, where T = mv 2 /2, m is the particle rest mass, v is the particle velocity, μ is the cosine value of the pitch angle α. This gives us: where n s and n 0 are the number densities at point s and the equator, respectively. B s and B 0 are the corresponding magnetic fields at two points (dipole model is used here). Particles with μ below the limiting value μ 0 * = (1−B 0 /B s ) 0.5 at the equator will mirror before reaching point s. The integral of f(T 0 , μ 0 ) over 0 < T 0 < ∞ and −1 < μ 0 < 1 is equal to 1. The fitted parameters in Table 1  Evidently, the total mass density in each event has a local peak at equator, which is consistent with the prediction of our wave frequency analyses. Power law distributions with α = −5.0 (black dotted dashed curve), −3.8 (green solid curve), and −2.9 (black dashed curve) obtained from Section 3 are also displayed in Figure 5d for comparison. These illustrate that, though the ions measured by the HOPE occupy only a small part of the ensemble due to the limitation of in situ detection, the distribution of ρ total is approximately in line with the power law distribution derived from frequency analyses. Additionally, the electron number densities derived from upper hybrid frequencies increase together with enhanced H + fluxes, which is denoted by two dashed vertical lines in Figures 4a, 4c, 4f, and 4h, respectively. The electron number densities have increases of 6−10 cm −3 , corresponding to an enhancement of 10%−25% from the ambient electron densities. Under the assumption of quasi-neutrality of plasma, the obvious enhancements of electron densities also support equatorial ion concentrations under the influence of EMIC waves. The observed mass density peaks are located in the postnoon sector with MLT ~13 hr where a negative α is possible to ap- pear (Denton et al., 2004Takahashi and Denton, 2007). If the peaks in Figure 2f are also taken to be the fourth and fifth (or sixth) harmonics, α will be larger than 6, which is unrealistic and, to our knowledge, nobody has reported values of α > 6.

Summary
Through frequency analyses of simultaneously observed toroidal harmonics and EMIC waves, equatorial mass density peaks are determined. The distributions of measured low energy ions support the mass density peak formation at the equator according to Liouville's theorem. An increase of perpendicular fluxes of ions (< 100 eV) heated by EMIC waves can promote mass density peak formation at the equator. This paper thus provides observational evidence that equatorial mass density peaks are formed due to equatorial ion concentrations triggered by EMIC waves. In turn, equatorial mass density peaks influence the excitation and propagation of magnetohydrodynamic waves, which will be investigated in future work.