A Wave Model and Diffusion Coefficients for Plasmaspheric Hiss Parameterized by Plasmapause Location

The scattering of electrons via plasmaspheric hiss whistler‐mode plasma waves has profound consequences for the dynamics of electrons in the inner terrestrial magnetosphere, including the radiation belts. Consequently, simulations of inner magnetospheric electron dynamics incorporate hiss wave models, though these models are often parameterized by quantities convenient to describe particle populations (e.g., L shell). However, recent studies have revealed that the spatial distribution of plasmaspheric hiss wave power is only weakly dependent on L shell. Instead, it is dictated by the density structure of the plasmasphere (including radial extent and azimuthal structure). In this work, we create a plasmaspheric hiss wave model, and corresponding particle diffusion coefficients, parameterized by plasmapause location instead of L shell, in order to quantify the importance of including plasmapause organization of hiss waves for inner magnetosphere models. Noticeable differences in electron scattering lifetimes are found when comparing L shell parameterized hiss and plasmapause‐parameterized hiss wave models on the timescales of days for 1 MeV electrons. Comparison of the lifetimes calculated from diffusion around the loss cone also show differences between L shell‐ and plasmapause‐parameterized hiss models. This implies that plasmapause parameterization of hiss waves may be important for modeling specific geomagnetic events.


Introduction
Plasmaspheric hiss is a broadband superposition of whistler-mode plasma waves located within and nearby the plasmasphere, a torus of cold plasma surrounding Earth. Hiss scatters electrons in pitch angle, facilitating their loss to the atmosphere and thereby playing a significant role in shaping inner magnetospheric electron populations, including relativistic radiation belt electrons (e.g., see Millan & Thorne, 2007). For this reason, electron loss by hiss wave scattering is a critical component of simulations of the inner magnetosphere (Albert et al., 2009;Fok et al., 2014;Jordanova & Miyoshi, 2005;Miyoshi et al., 2006;Shprits et al., 2008;Subbotin & Shprits, 2009). Simulations often include the physics of hiss-induced particle scattering using statistical maps of hiss wave characteristics (e.g., intensity, spectral shape) parameterized by L shell, magnetic local time (MLT), geomagnetic activity level (Glauert et al., 2014;Meredith et al., 2007;Orlova et al., 2014Orlova et al., , 2016Tsurutani et al., 2015Tsurutani et al., , 2014, or solar wind parameters (Aryan et al., 2017).

Journal of Geophysical Research: Space Physics
10.1029/2019JA027415 determined by a balance of solar wind-driven convection, corotation with Earth's magnetic field, refilling from ionospheric outflow, and the time history of all these processes (Carpenter & Lemaire, 2004).
In this work, we begin quantifying the difference between hiss parameterized by plasmapause location (L pp -sorted) and hiss parameterized by L shell (L-sorted) for models of hiss wave scattering by (i) using Van Allen Probes plasma wave observations to create a database of hiss wave power parameterized by wave frequency, magnetic local time, plasmapause location (using both L PP , the distance of the plasmapause from Earth, and ΔL PP , the distance from the plasmapause), and the Kp geomagnetic index; (ii) using this observational database to produce a L pp -sorted hiss wave model amenable to calculation of diffusion coefficients; (iii) producing diffusion coefficients based on the wave model; and (iv) applying the L-sorted and L pp -sorted diffusion coefficients to model an idealized geomagnetic storm using the Versatile Electron Radiation Belt (VERB) code (Subbotin & Shprits, 2009), performing one-dimensional pitch angle diffusion simulations, and comparing code outputs.
This work focuses on plasmaspheric hiss, defined here as hiss found within the plasmasphere at frequencies between 150 Hz and ∼2 kHz. Other hiss types will be parameterized in future work, including plume hiss (Li et al., 2019), low-frequency hiss (Li et al., 2013), lightning hiss (Meredith et al., 2007), and exohiss (Zhu et al., 2015(Zhu et al., , 2019. By building a wave model for each hiss type separately, they can be included or excluded from VERB code runs to quantify their relative importance to inner magnetospheric dynamics.
Section 2 describes the wave data and its processing. Section 3 describes the parameterized wave model built from the observations. Section 4 describes the calculation of diffusion coefficients, and section 5 treats the VERB modeling. Conclusions are presented in section 6.

Data Set Creation
This study uses data from the Van Allen Probes mission (Mauk et al., 2013). These two identical spacecraft have elliptical orbits about the Earth with perigee near 600 km and apogee near 6 R E (Earth radii). Their orbits are within 20 • of the geomagnetic equator, and each has an orbital period close to 9 hr. Over ∼2 years, their orbital line of apsides precesses through all MLT. The spacecraft are spin-stabilized with ∼11 s spin period.
Instrument data used in this study are from the Electric Fields and Waves (EFW) instrument (Wygant et al., 2013) and the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) suite . These instruments record and process measurements made by six electric field probes, a three-axis search coil magnetometer (SCM), and a three-axis fluxgate magnetometer. The data products used by this study include the spacecraft potential (32 samples/s), the DC-coupled magnetic field (64 samples/s), density determined from a combination of the upper hybrid frequency (1 sample/6.5 s) and spacecraft potential, and onboard calculated wave power spectra of each SCM axis (65 pseudologarithmically spaced frequency bins from ∼ 2 Hz to ∼12 kHz, 1 spectra/6 s). Wave planarity and ellipticity data, derived from onboard calculated cross-spectral data (1 spectra/6 s), are also used.
Data from both Van Allen Probes are used. For Van Allen Probe B, data from 1 November 2012 to 31 January 2018 are used. For Van Allen Probe A, data from 1 November 2012 to 31 May 2016 are used. Data after May 2016 on Van Allen Probe A are not used, as accumulated radiation damage to that spacecraft's electric field sensor preamplifiers compromised their ability to accurately measure spacecraft potential soon after that date.
The Olson-Pfitzer quiet time magnetic field model (Olson & Pfitzer, 1974) is used to determine L shell values at any given time and spacecraft location. In general, this model is appropriate at high L shells (L > 4.5) during geomagnetically quiet times and at low L shells (L < 4.5) during both active and quiet times. Plasmaspheric hiss, by definition, remains within the plasmasphere and so is present at high L shells during quiet times (extended plasmasphere) and at low L shells during active times (eroded plasmasphere). Therefore, the Olson-Pfitzer quiet time model is appropriate for plasmaspheric hiss studies. L shell is used instead of L* because the relevant quantity for plasma waves is their radial distance from Earth (or the plasmapause) at the geomagnetic equator rather than a particle drift invariant (e.g., Koller et al., 2009, and references therein).
The plasmapause is identified using the method described in Malaspina et al. (2016). This method uses plasma density derived from spacecraft potential measurements calibrated each orbit against the measured upper hybrid resonance frequency, combined with the Moldwin et al. (2002) criteria that density change by 5 times or more over 0.5 L shell. When multiple density gradients satisfying this criteria were found, the one closest to Earth is designated as the plasmapause.
Isolating plasmaspheric hiss wave power from all other phenomenon detected by the Van Allen Probes SCM requires excluding some data from the analysis. Data recorded at L < 1.6 were not considered. Data from half orbits where no plasmapause was detected are excluded because ΔL pp is undefined for those data. Times during spacecraft maneuvers, significant spacecraft surface charging (|V sc | > 20 V), and times when the spacecraft were in Earth eclipse were removed. Wave power > 2 kHz was excluded, as those higher frequencies contain significant contributions from lightning-generated whistler-mode waves (Meredith et al., 2007), which we do not wish to include in the current model. Wave data outside the identified plasmapause L shell or when the corresponding plasma density measurement was < 50 cm −3 were excluded from consideration, to aggressively filter out chorus wave power.
Several filters were applied to the remaining spectral wave data to exclude wave modes not under consideration and to separate signal from noise. Spectral bins dominated by magnetosonic wave power were excluded by removing from consideration spectral data (by time and frequency bin) with high compressability B wave || /B wave total > 0.6. Spectral data had to meet the following criteria: planarity > 0.2, ellipticity > 0.7 (Li et al., 2015), and signal to noise ≥ 5 (using the empirical SCM noise function derived in Malaspina et al., 2017). The low planarity threshold (0.2) is justified here because planarity is being used only to exclude nonhiss waves (such as magnetosonic waves). A higher planarity threshold (0.5 or greater) would be required if the onboard cross-spectral data were being used to estimate wave vectors. After these exclusions, the remaining database includes 1.9 × 10 8 spectral data samples.

Plasmaspheric Hiss Data
The plasmaspheric hiss data examined here are wave power spectral density (units of nT 2 /Hz) binned by five quantities: frequency, distance of the plasmapause from Earth (L pp ), distance from the plasmapause (ΔL pp ), MLT, and the Kp geomagnetic activity index. The 65 pseudologarithmic spaced frequency bins defined by the EMFISIS onboard spectra are used for frequency binning. Four L pp bins are used, covering the range of possible plasmapause locations observable by the Van Allen Probes: 2 ≤ L pp < 3, 3 ≤ L pp < 4, 4 ≤ L pp < 5, and 5 ≤ L pp 6. The 25 bins for ΔL pp span a range of −5 to 0, with a bin width of 0.2 L. Six MLT bins with 4 hr width are used: 0 < MLT ≤ 4, 4 < MLT ≤ 8, and so on. Finally, six Kp bins are used, where the first five have a span of 1 Kp (e.g., 0 ≤ Kp < 1, 1 ≤ Kp < 2). The final Kp bin includes Kp≥ 5. This final bin width was selected to ensure sufficient data for meaningful statistics in this bin. This was necessary because geomagnetic conditions resulting in Kp≥ 5 are rare during the Van Allen Probes era.  Figure 1d, 2 < L PP ≤ 3. The data shown are for 8 < MLT ≤ 12 and 1 ≤ Kp < 2, but the data look similar for other MLT and Kp bins. Contours of amplitude are plotted over the amplitude range shown. In each case, the amplitude contours trace an elliptical shape centered at a few hundred Hz in frequency. The data show properties consistent with prior studies (e.g., Malaspina et al., 2017): the wave power peaks near 400 Hz, at a radial distance approximately between the Earth and the plasmapause. When the plasmapause is eroded, the wave power is compressed into a smaller radial extent, and the amplitude increases.
The thick black horizontal lines indicate 150 Hz. While wave power below 150 Hz is plotted in Figure 1, those data are not used for the plasmaspheric hiss wave model derived in this work. Hiss wave power below 150 Hz is considered low-frequency hiss (Li et al., 2013;Malaspina et al., 2017;Ni et al., 2014) and will be considered separately in future work.

Plasmapause-Parameterized Wave Model
Using the plasma wave data collected by the Van Allen Probes as described above, we create a plasmaspheric hiss wave model, parameterized by frequency, L pp , ΔL pp , and Kp.
The shape of hiss wave power spectral density (PSD) profiles with respect to wave frequency are found to depend strongly on L pp and ΔL pp and weakly on MLT and Kp. PSD profile amplitude is found to vary with L pp , ΔL pp , MLT, and Kp. Therefore, a parameterization is chosen such that the wave model PSD frequency profile shape is determined by L pp and ΔL pp , while the appropriate amplitude scaling for that shape is determined by L pp , ΔL pp , MLT, and Kp.
To create the wave model, PSD data are first averaged over MLT and Kp. Figure 2 shows these averaged PSD versus frequency profiles (solid lines) as a function of L pp and ΔL pp . When PSD profiles have similar peak frequencies for different ranges of ΔL pp , they are combined into a single profile. The resulting PSD frequency profiles were fit with an analytic function (dotted lines) (a piecewise seventh-order polynomial function) to facilitate diffusion coefficient calculations. Fits are carried out separately for frequency ranges before and after the maxima of the PSD profiles peak as follows: We then normalize over the obtained PSD profiles such that the wave amplitude is unity when integrated over the frequency range from 150 to 2,000 Hz.
All values of peak and fitted polynomial coefficients can be found in Table 1 and on the ftp server listed in the Acknowledgments). Note that the fitted PSD profiles may be negative at > 1, 000 Hz where the observed     Figures 2 and 3 ΔL pp L pp range for PSD on Figure 3 L pp range for PSD on Figure 2 < 3 3-4 4-5 PSD value are small. When using the normalized PSD profiles, a minimum value of 1×10 −10 is imposed. The wave amplitude scaling appropriate to each normalized PSD profiles is obtained by comparison with the observational hiss database. This combination of PSD profile fitting and wave amplitude scaling allows us to fully parameterize the hiss wave distributions. Table 2 contains values of the integrated wave amplitudes required to scale the PSD profiles in order to reproduce Figures 2 and 3. Figure 3 shows the comparison between the observed and parameterized hiss wave power distributions and their normalized difference as functions of ΔL pp and L pp . The normalized difference is defined as It is shown that the observed PSD distributions of hiss wave power are well modeled by the fitting results, and the normalized differences between the fits and observations are close to 0.
In addition to the L pp -sorted wave model, a traditional hiss wave model was constructed based on the same wave data (parameterized by L shell, instead of ΔL pp and L pp ). PSD frequency profiles and amplitude scalings were obtained using methodology analogous to that described for the L pp -sorted model.

Diffusion Coefficient Calculations
The Full Diffusion Code (Ni et al., 2008;Shprits & Ni, 2009) is used to calculate plasmaspheric hiss diffusion coefficients. The polynomial fits of the observed L pp -sorted hiss PSD profiles described in section 3 are used as frequency inputs into the code. The Denton et al. (2006) model is used to define plasma number density. The wave normal distribution is assumed to be a Gaussian distribution (Glauert & Horne, 2005) with a peak at tan (  3 < L pp ≤ 4. As expected, the energy diffusion coefficient of the hiss waves is relatively small in comparison to the pitch angle diffusion coefficient. Figure 4c shows the pitch angle diffusion coefficient at the same L shell and Kp index but for a different plasmapause location 4 < L pp ≤ 5. The difference between pitch angle diffusion coefficients (Figure 4g, difference between Figures 4a and 4c) indicates that stronger scattering by hiss waves when the plasmapause is closer to Earth primarily affects electrons with energies from 100 keV to 1 MeV. Figure 4d shows pitch angle diffusion coefficients for the hiss model constructed using traditional L-sorted methodology. The difference between pitch angle diffusion coefficients from L pp -sorted and L-sorted models (Figure 4h, difference between Figures 4a and 4d) is clear, for the same range of electron energies. The differences reach ∼ 2 day −1 . Thus, one can expect that electron distribution dynamics in diffusion simulations will depend on the wave data parameterization (sorting) approach of the hiss model (L-sorted vs. L pp -sorted).
In addition, the variation can be significant on the timescales of geomagnetic storms (a few days).
The calculation of diffusion coefficients allows us to directly analyze electron lifetimes. The electron lifetime is calculated as = D ( LC ) −1 , where LC is the loss cone (Shprits et al., 2006). Electron lifetimes, as functions of kinetic energy, for L pp -sorted plasmaspheric hiss waves at L = 2.4, L = 3.4, and L = 4.4 under different Kp conditions, are shown in Figure 5. Electron lifetimes calculated using the L-sorted hiss wave model are plotted as dotted lines for comparison. There are a few features of interest in the calculated electron lifetimes. First, the shape of the lifetime curves for different L pp locations and for the L-sorted model are similar. Their shapes significantly depend on L rather than Kp, consistent with previous studies (Orlova et al., 2014). Second, electron lifetime curves for different L pp locations have quantitative differences, which suggest the importance of using L pp -sorted hiss wave models.

Pitch Angle Diffusion Simulations
To quantify the potential impact of the wave data sorting approach on modeling results, we perform VERB code simulations of an idealized geomagnetic storm in 1-D mode using pitch angle diffusion only. This approach allows us to focus on the impact that different hiss wave models may have on the simulation, as we neglect other processes such as radial diffusion or acceleration and scattering by chorus waves. Since the electron acceleration by hiss waves is ineffective, we ignore energy diffusion. Neglecting radial and energy diffusion, the Fokker-Planck equation (Schulz & Lanzerotti, 1974) that describes the evolution phase space density ( ) can be written as where T(sin( )) ≈ 1.38 − 0.32(sin( ) + √ sin( )) is a function that corresponds to the bounce frequency, approximated following Lenchek et al. (1961); is the equatorial pitch angle; D is the pitch angle diffusion coefficient; and defines the lifetime of the particle inside the loss cone where it is equal to quarter of the bounce period.
We perform simulations on a grid of ∈ [0.1 • , 89.5 • ] linearly distributed among 101 points. The Dirichlet boundary condition is equal to 0 at = 0.1 • , and the Neumann boundary condition is the derivative equal to 0 at = 89.5 • . The initial condition is an isotropic phase space density distribution ( ( ) = 1). The energy of the electrons is 1 MeV. A similar simulation setup was used in previous studies (e.g., ). Figure 6 shows the results of two simulations at fixed L shell equal to 3.5. These simulations are distinguished by the use of two models for hiss waves described above (L-sorted and L pp -sorted). Different diffusion coefficients are applied for each wave model. Both simulations were performed for 7 days, using a prescribed Both simulations are given an initial 4 day period of constant Kp = 2, which corresponds to a constant location of the plasmapause at L pp = 4.68. During this period, the diffusion coefficients are held constant in both simulations as the simulations reach steady state conditions (no change in phase space density gradient; see Figures 6c and 6d). Hence, the electron dynamics of the remainder of the modeled time period are defined by variation of the diffusion coefficients, which are determined by the wave data sorting approach (L pp -sorted vs. L-sorted). The idealized geomagnetic storm starts on Day 4 with increasing Kp index and corresponding compression of the plasmapause spanning 1 day. On Day 5, the Kp index is allowed to decrease, followed by an expansion of the plasmapause on Day 6. Such conditions are typical for geomagnetic storms.
Comparison of the simulation results shows that the electron decay rate evolution during the storm time is noticeably different (see Figures 6c and 6f). The difference is also visible in evolution of the phase space density profiles for 1 MeV electrons (see Figures 6d and 6g). The step-like increases of decay rate in Figure 6f are the consequence of the small number (4) of discrete L pp bins that are used to calculate the diffusion coefficients. The small number of bins is due to statistical limitations of the wave data. The presence of these steps does not alter the clear difference in the simulated electron dynamics found by using different hiss models.
Based on these results, it is expected that a simulation of a realistic storm using the L pp -sorted hiss model will result in different electron dynamics on the time scale of the duration of the storm (days) compared to simulations using traditional L-sorted hiss wave models. Simulations using the L pp -sorted hiss model may reveal otherwise hidden variation of the electron distribution during storms and may also lead to a different simulated balance between acceleration and loss processes due to changes in the phase space density gradient. However, pitch angle scattering by hiss waves is only one of many processes that define the dynamics of the electrons. Future simulations will include contributions from other very low frequency (VLF) waves via pitch angle, energy, mixed diffusion, and radial diffusion driven by ultralow-frequency (ULF) waves.

Conclusion
In this study, Van Allen Probes observations of plasmaspheric hiss (organized by frequency, ΔL pp , L pp , MLT, and Kp) were compiled for the time period 2012-2018. From these data, an empirical hiss wave model was constructed for hiss parameterized by ΔL pp , L PP , and Kp. Corresponding pitch angle and energy diffusion coefficients (including mixed terms) were calculated, as were corresponding electron lifetimes.
The pitch angle diffusion coefficients for the L pp -sorted empirical hiss model showed noticeable differences when compared with diffusion coefficients calculated for the same L shell but using traditional L-sorted hiss parameterization.
A 1-D mode of the VERB code with idealized geomagnetic storm conditions was used to quantify differences in electron lifetimes as determined using diffusion coefficients calculated using the L pp -sorted hiss wave model and the L-sorted hiss wave model. Clear differences were found over time timescales of the geomagnetic storm (few days). It is expected that implementation of the developed L pp -sorted hiss model into simulations of realistic storms or periods of quiet geomagnetic activity will result in different electron dynamics in comparison with traditional hiss models. However, realistic simulations must include additional processes (e.g., radial diffusion, chorus wave acceleration, and scattering). Based on the results of this work, the L pp -sorted hiss model will contribute to the balance between acceleration and loss in a qualitatively different manner than traditional L-sorted hiss models.
Future studies will expand upon the current model by (i) simulating more realistic geomagnetic variation time histories, (ii) utilizing the 3-D mode of the VERB code, and (iii) developing L pp -sorted parameterizations for other hiss types such as exohiss and low-frequency hiss.