2.1 Chorus Wave Field Model

To get meaningful results from test particle simulations, we need to construct a wavefield model that resembles chorus elements as observed by spacecraft. In Figure 1, we show an example of a rising lower band chorus element, which we used as a reference for our model. The spectrogram of the total magnetic field power spectral density in Figure 1a was obtained from the multicomponent continuous burst-mode measurements of the EMFISIS Waves instrument on the Van Allen Probes B spacecraft. In this mode, waveforms of all six electromagnetic field components were captured at a sampling rate of 35,000 samples per second. The waveforms were analyzed using 512-point Fast Fourier Transform with the Hann window, 50% overlap, and averaging over three neighboring spectra. Calibrated magnetic field waveforms were transformed into the coordinates linked to the background magnetic field (Figure 1b). They were subsequently band-pass filtered between 0.1 and 0.49 of the local electron cyclotron frequency. The instantaneous frequency (Figure 1c) was then obtained as the time derivative of the instantaneous phase of an analytic signal reconstructed using the Hilbert transform. Frequency data from the time intervals where the amplitude dropped below were cut off (Santolík, Kletzing, et al., 2014). The plasma frequency was obtained from measurements of the upper hybrid frequency (Kurth et al., 2015).

The example element is about long and rises from frequency of about to . After the initial long and weak subpacket, two well separated, high-amplitude packets appear, reaching about of the background field magnitude. With increasing wave frequency, the subpacket structure becomes less regular. The frequency follows a linear trend with irregularities appearing not only in between, but also inside subpackets, suggesting that this behavior cannot be caused purely by Hilbert transform applied on low-amplitude sections. The wave remains quasiparallel during the measurement, with wave normal angle at peak amplitudes starting at about and steadily decreasing.

To capture the observed behavior of the depicted chorus element, we use the parallel whistler-mode chorus model of Hanzelka et al. (2020), which is based on the nonlinear growth theory of chorus emissions (Omura, 2021; Omura et al., 2009). This model solves the so-called chorus equations to compute the evolution of plane wave frequency and amplitude in the source of a subpacket, and after the subpacket reaches saturation amplitude and decays, a new subpacket starts at a point in time and space determined by the propagation of the resonant current to upstream regions (Helliwell, 1967; Trakhtengerts et al., 2003). To obtain the electromagnetic field in the downstream, we used transport equations for wave magnetic field and angular frequency to propagate the wave in one spatial dimension along the magnetic field line. These two equations are (Nunn, 1974; Omura et al., 2008)

(1)

(2)

The evolution of the frequency is a pure advection with group velocity , while the amplitude equation is modified by a nonlinear term on the right-hand side which includes the resonant current component parallel to the wave electric field. These equations are solved by the upwind method with a time step and a spatial step . On this grid, the Courant-Friedrichs-Lewy condition is satisfied over the whole domain while the numerical diffusion is kept reasonably slow (LeVeque, 2007).

The initial boundary values for each subpacket are given by the coupled set of modified chorus equations (Omura et al., 2009; Hanzelka et al., 2020). For detailed description of these equations, the notation and the complete set of input parameters, see Appendix A, where the results of Hanzelka et al. (2020) are briefly summarized.

To obtain the phase of the wave magnetic field vector, frequencies and wave numbers are integrated for each subpacket (we first integrate over time in the source, then over space away from the source). During particle simulations, these integrated phase angles are bilinearly interpolated at the position of each particle to obtain the angle between and . The phase difference between adjacent subpackets is not specified, as it changes with the movement of the source region.

The final model is presented in Figure 2 as plots of magnetic field amplitudes and wave frequencies. In the following particle simulations, we assume waves propagating in both directions from the magnetic equator, as shown in the amplitude plot in Figure 2b. Figures 2c and 2d show the frequency growth within the subpackets, ranging from to , where is the equatorial angular gyrofrequency. The simulation was stopped when the starting frequency of the next subpacket would become higher than ; this upper frequency limit corresponds to observations of lower band chorus with a spectral gap, but it is arbitrary in our simulation procedure as there is no natural cutoff for strictly parallel propagating whistler-mode waves. The simulated chorus element features a realistic subpacket structure with an irregular growth of the wave frequency (compare with Figure 1c) and with an upstream shift of the source region—these properties have been observed in numerical simulations (Hikishima et al., 2009) and also in spacecraft measurements (Demekhov et al., 2020; Foster et al., 2017; Santolík, Kletzing, et al., 2014).

2.2 Backward-in-Time Test Particle Simulation

To obtain a high-resolution velocity distribution at a given point in space and time, we use the backward-in-time simulation method (Nunn & Omura, 2015). In this approach, the equation of motion for electrons in an electromagnetic field is solved from a fixed final point (after the interaction with chorus) to an arbitrary initial point at which the wavefield disappears. The phase space densities corresponding to the velocity of electrons at the end of the simulation are mapped to the fixed final point with the use of Liouville's theorem (Goldstein et al., 2001): the velocity distribution function is constant along the trajectories of a Hamiltonian system. Therefore, if we know the initial, unperturbed distribution at all points along the magnetic field line and we choose the velocities of electrons so that we cover uniformly the velocity space at the final point, we can obtain the phase space density (PSD) of the electron distribution after the interaction with a chorus wave (see Figure 3 for a schematic explanation of the ideas behind the method). This method is obviously not self-consistent, as it requires a precalculated wavefield which does not respond to the changes in the electron velocity distribution. However, unlike in self-consistent kinetic simulations, we can choose to analyze an arbitrary section of the velocity space, thus increasing the resolution. And, maybe even more importantly, the initial distribution can be chosen a posteriori, meaning that we can obtain the perturbations for all possible distributions with just one simulation.

The test particles are advanced backward-in-time by a type of the Boris algorithm (Higuera & Cary, 2017) which does not use the common particle phase approximation (Zenitani & Umeda, 2018). The algorithm preserves phase space volume, which is an important property that allows us to safely apply Liouville's theorem in the numerical simulation. The effect of the mirror force on electrons is included through additional perpendicular components of the background magnetic field (Katoh & Omura, 2006). We also tested that a finer grid in the wave model has only negligible influence on the structure and magnitude of observed perturbations in the electron velocity distribution (peak value of relative PSD change in the simulation with one subpacket changed by less than ). However, some grid-dependent changes can be observed, especially at the edges of electromagnetic holes, as the chaotic trajectory of the resonant electrons strongly depends on their gyrophase at the onset of the resonant interaction.

Since, the phase space density is supposed to be preserved in our physical system, the initial velocity distribution function must evolve adiabatically along field lines. We chose a distribution of the form (Summers et al., 2012)

(3)

(4)

which is a bi-Maxwellian distribution in momenta , . In velocities , , the distribution takes form

(5)

(6)

The initial distribution in the particle phases is set to be uniform. stands here for the number density of hot electrons at the equator, (or ) and (or ) are the equatorial thermal velocities (or relativistic momenta) and is the magnetic dipole field strength along the field line. As noted by Kuzichev et al. (2019), this distribution has a different asymptotic behavior at relativistic velocities than an anisotropic Jüttner distribution. However, we are interested in particles at resonant velocities, which are only weakly relativistic. The loss cone is also not included in the model, because particle behavior at low pitch angles plays no significant role in the nonlinear growth theory.

The input parameters common to all simulation runs are as follows. The ratio of the plasma frequency to the equatorial gyrofrequency is (based on measurements shown in Figure 1) and the ratio of the density of hot electrons to the density of cold electrons is , which is a reasonable value for populations with thermal velocities of a few units to tens of keV (Juhász et al., 2019; Walsh et al., 2020). The parallel and perpendicular thermal relativistic velocities at the equator are and , resulting in an equatorial temperature anisotropy . The background magnetic field is assumed to be dipolar with equatorial strength at the Earth's surface of and the particles are propagating on the -shell with . The corresponding equatorial electron gyrofrequency evaluates to . These parameters are the same as in the wave model (see Appendix A).

All test particle simulations were done on a uniform grid in () space, where is the angle between the wave magnetic field vector and the perpendicular velocity of electrons. The parallel velocities were sampled in the range with 256 points, perpendicular velocities in the range with 256 points, and phases with 64 points (particles with initial velocities larger than are excluded from the sample). To confirm that this resolution is sufficient (Voitcu et al., 2012), we limited the range in velocities to focus on the resonance region of the first subpacket, and doubled the number of gyrophase grid points. Apart from some finer features associated with the phase mixing of lower and higher density regions, there was no noticeable difference. In Movies S1–S3, which depict the evolution of phase space hole, the parallel velocities are sampled on a much smaller range and with 1,024 points, the phases are sampled with 512 points and the perpendicular velocities are fixed to a single value. Each frame in these animations requires a new simulation run.

2.3 Particle Count Analysis

Particle count measured by a particle detector with a geometric factor over a time interval can be expressed as

(7)

where is the 3D velocity distribution. In our analysis, we work with a velocity distribution averaged over gyrophases, , since it can be safely expected that during the time interval , the detector will uniformly sample particles across all gyrophases. We define the number flux as

(8)

The width of the azimuthal angle and polar angle of a single pixel of an electrostatic particle detector is included in the geometric factor , as well as the relative width of a velocity (energy) level . All detectors listed in Table 1 in Section 3.3 have an approximately constant ratio of the width of the energy bin to the energy at the center of the bin across all levels. The energy geometric factor given in the specification of the detectors is related to the velocity geometric factor as .

Table 1. The Geometric Factor and Resolution in Time, Energy and Polar Angle of Particle Detectors on Selected Spacecraft

	FPI	HOPE	PEACE HEEA	MEP-e
MMS	RBSP	Cluster	Arase
Angular resolution (deg)	11.25	18	15 (3.75)	22.5 (3.5)
Energy range (keV)	0.01–30	0.001–50	0.01–26.4	7–87
Energy levels per sweep	32	36	72	16
Time per energy level (ms)	0.23	21	4	15.6
Geometric factor per bin ()	2.1	0.2	6.0 (1.5)	0.7

• Note. The listed parameters are based on Pollock et al. (2016) for FPI (Funsten et al., 2013), for HOPE (Johnstone et al., 1997), for PEACE, and on (Kasahara et al., 2018) for medium-energy particle (MEP). The cylindrical MEP instrument has seemingly the best resolution, but it provides only 16 azimuthal channels with large dead zones, thus providing data from of each channel. For the rest of the instruments, the angular resolution refers to the polar angle bins.

The number of particles required to measure a significant decrease in phase space density can be defined in a variety of ways. Here we use the one-sigma (68%) confidence interval for a Poisson distribution approximated as , where is the number of particles measured on the perturbed distribution in a given energy-angle bin, and we require that the number of particles measured in the corresponding bin of the initial (unperturbed) distribution is at least by larger. Since relative change in number flux or PSD is directly proportional to the difference in number of measured particles, we get , where the subscript “0” denotes the initial state. The required number of particles measured in a single bin of the perturbed distribution can be then defined as . If we were to measure an increase in phase space density, we would relate the difference to the higher value of PSD, which would be the perturbed value, so we can generalize the requirement to .

We must also be mindful of the absolute number of measured particles. If we calculate the so-called exact confidence interval (Meeker et al., 2017), it will have different values for different , for example: for , for , and for . We can see that our approximated upper bound is noticeably underestimated for low , so we should always require in our assessment of measurability to make the approximations valid.

If we decided to work with a two-sigma (95%) confidence interval, the required number of particles would be about four times larger, with the upper bound of the exact confidence interval showing underestimation with nearly the same relative magnitude as in the one-sigma case.

3.1 Response of the Electron Distribution to a Single Chorus Subpacket

We start by analyzing the evolution of the equatorial hot electron distribution during interaction with the first subpacket of our wave model from Section 2.1. To numerically simulate the wave-particle interactions, we use backward-in-time test particle simulations as described in Section 2.2. The final point of particle propagation (initial in the backward sense) is , , lying on the boundary between the first and the second subpacket. The time step is chosen to be a fixed value of . The particles remain near the equator where the magnetic field is almost constant, meaning that each electron gyroperiod is sampled with about 120 points.

During the resonant interaction, an electromagnetic hole forms in the space (see Movie S1 where the evolution of the electromagnetic hole in the first two subpackets is captured). Particles trapped in the hole are transported to lower parallel velocities, while resonant particles that flow around the hole increase their parallel velocity. For a highly anisotropic distribution, this translates to an increase of the pitch angle and the kinetic energy for trapped particles, and a decrease of and for untrapped particles (see Figures 4a and 4b for a sketch of the particle scattering and transport, and see Movies S2 and S3 for changes in the kinetic energy and pitch angle during the propagation of electrons through the first two subpackets). After integrating over particle phases, we obtain a 2D velocity distribution where the hole is not apparent anymore, and the perturbed part of the distribution appears as a stripe along the relativistic resonance velocity curve (see Figure 4c). The resonance velocity is defined as , where is the Lorentz factor and is the wave number. As expected from the inspection of the electromagnetic hole, the stripe consists of an increased phase space density part which is located close to resonance velocity corresponding to the initial frequency of the subpacket, and a decreased density part around resonance velocity corresponding to the wave frequency at the end of the subpacket (see Figure 4d). After integration over perpendicular velocities, we observe a smooth step-like feature in the reduced 1D distribution located around the resonance velocity (see Figure 4e).

3.2 Response of the Electron Distribution to a Full Chorus Element

In the analysis of the resonant interaction of electrons with a chorus element, we proceed in the same fashion as in the previous case of interaction with a single subpacket. Certain particles may now reach regions (or, in the forward time flow, come from regions) with significantly stronger magnetic field, so the time step must be appropriately smaller. A fixed value of was chosen, which translates to about 80 steps per electron gyroperiod at the boundary of the simulation domain where the magnetic field is about three times stronger than at the equator. Furthermore, some particles may reflect at their adiabatic mirror point and interact with the chorus element that propagates toward lower values.

We set the final point to , , which is right after the end of the last subpacket at the equator. The stripes created by the interaction with the subpackets are now distorted by adiabatic motion of particles that resonated with these subpackets farther from the equator (see Figure 5a; also see Movie S4 for subpacket-by-subpacket evolution of the perturbed distribution at the equator). In the 1D reduced distribution, the stripe-like perturbation related to single subpackets remains visible in those cases where the adiabatic motion of resonant electrons straightened the stripes and aligned them thus with the lines of constant parallel velocity (see Figure 5b). Subtraction of the initial distribution reveals an increase in phase space density at higher parallel velocities and a more prominent density decrease at lower velocities, especially around the resonance velocity corresponding to the last subpacket (see Figure 5d). Further analysis of the perturbations in the 2D distribution reveals that the stripes are overlapping, which makes them hard to distinguish, especially at resonance velocities corresponding to higher frequencies (see Figure 5c; also see the second half of Movie S1 to observe the mixing of phase space density between the first and second subpacket).

3.3 Analysis of the Required Resolution of Particle Measurements

The first step in this measurability analysis is to determine the time interval and the range of energies and pitch angles in which we observe large phase space density perturbations and, at the same time, large particle fluxes. If we assume the number of particles measured in each energy-angle bin to be Poisson distributed, then the statistical uncertainty of the measurement can be represented by the relative standard deviation . So, in order to make the measurement of a PSD perturbation significant at the one-sigma level, we need at least particles in each bin, where is the number flux (see Section 2 for a more detailed explanation of confidence intervals and the definition of number flux). This means that unlike in Figure 4d, we need to look at the relative perturbations of the phase space density and associated fluxes. The relative change in PSD after interaction with the first subpacket is shown in Figure 6a. The values of the number flux shown in Figure 6b, which are proportional to the particle counts, reach maximum at , . The quantity we use to determine the most suitable phase space region for measurements is . A plot of this quantity is presented in Figure 6c for the case of the first subpacket. By looking at the state of the perturbed distribution after interaction with each subpacket, we determined that the highest values of are found between subpackets 4 and 7. The case of the fifth subpacket is shown in Figure 6d. It needs to be stressed here that in our nonself-consistent simulation, we cannot determine how fast will the stripe on the left (associated with the first subpacket) decay, and so we should focus only on the newly formed generated right-hand stripe where the nonlinear interaction is taking place.

From Figure 6d we conclude that the optimal energy range for measurement of the phase space density decrease associated with growth of chorus subpackets starts at around and ends at about . In this energy range, the pitch angle width of the stripe is about 6°. A detector that always resolves the stripe should have a pitch angle resolution of half this width, which is 3°. Because the stripe does not lie on a constant pitch angle, the width of energy bins is also limited, but does not pose a major constraint.

In the second step of the measurability analysis, we compare the size and magnitude of the simulated perturbations with the resolution of an actual spacecraft instrument. Table 1 lists electrostatic particle detectors on several currently or recently operating spacecraft with science objectives related to wave-particle interactions in the magnetosphere. We can see that only the medium-energy particle-electrons (MEP-e) particle analyzer on the Arase spacecraft (Kasahara et al., 2018) has a suitable energy range, so we used this detector as a reference in further analysis. The MEP-e is a disklike electron spectrometer which filters incoming electrons by energy with an electrostatic analyzer and then uses avalanche photodiodes for particle detection. The time resolution of is close to the time duration of one chorus subpacket, meaning that we can use data from only one energy level to avoid smearing of the density stripe we wish to observe. With 16 logarithmically spaced energy levels from to , only the highest energy level is suitable to observe the structure depicted in Figure 6d (marked by a black arrow). The FPI detector on the MMS spacecraft and PEACE on the Cluster spacecraft provide better time resolution, but it compromises the particle counts they can obtain. For the purpose of comparison with our simulation, we assumed that the detector is oriented in such a way that it covers full in the pitch angle. Thus, the angular resolution in Table 1 is taken to be equivalent to pitch angle bins, although better resolution could be achieved with a different orientation and a more limited range of pitch angles. Changes of effectivity of the detector as a function of energy are not significant and are not considered here.

Since the measuring capabilities of electron analyzers on spacecraft are restricted to recording the particle count as a function of their kinetic energy and pitch angles, we will from now on inspect the perturbations of the electron distribution using these variables. Furthermore, all aforementioned simulation results represent the distribution at a single instant in time, while spacecraft particle analyzers typically require tens to hundreds of milliseconds to sweep through the relevant range of electron energies. To mimic the particle instrument behavior, we performed 16 backward particle simulations with the final time point being at the end of the fifth subpacket, and in each simulation, the final time point moved back by , adding up to a total time interval of . The chosen energy bin from to was covered with 16 logarithmically spaced grid points and we used 256 linearly spaced points for the pitch angle quadrant (64 points per single MEP-e angular sector) and 128 linearly spaced points for particle velocity phases.

Unfortunately, due to the large dead zones of the azimuthal channels of the MEP-e analyzer, it is immediately clear that we cannot resolve the perturbation in flux, because we have only one angular bin per a channel (see Table 1). Therefore, we had to base our simulation on a hypothetical detector which preserves the parameters of MEP-e, but has more azimuthal channels and negligible dead zones. We proceeded to refine the pitch angle bins to find out a compromise between high particle counts and good angular resolution. The ratio fluctuated as we were refining the angular grid, depending on whether the energy-angle bin fell into the center of the stripe or on the edge. However, with bin size or smaller, the ratio was growing steadily, signifying that by further resizing the bins, we are only losing particle counts and not improving the resolution of the stripe anymore. For bins of size between and , the ratio fluctuated between 1 and 2, with expected number of particles falling down to 3 or 4 per bin in the region of the PSD decrease. To achieve reasonably large particle counts, and also to prevent errors in estimation of confidence intervals, we further had to increase the geometric factor 10 times with respect to the original pitch angle resolution (i.e., the effective aperture of the detector has been increased 10 times), making it comparable to the PEACE detector.

In Figures 7a–7c, we can see an example with angular bin size of . A decrease in phase space density is observed between , followed by fluctuations at larger values of pitch angle. Thus, with the increased geometric factor and angular resolution, the phase space density stripe can be resolved, and the particle counts shown in Figure 7b are above 50. As shown in Figure 7c, the required particle counts for one-sigma significance are met at the lowest point of the dip. We can reach even , corresponding to a two-sigma significance. The points where the nonlinear wave growth is taking place are marked by green diamonds in the plot.

To get a more complete picture of the stripe as it could be seen by particle detectors, we further assumed four detectors operating simultaneously at four logarithmically spaced energy levels from to . We preserved the accumulation time step and energy level spacing as given in the MEP-e specifications, and we used the improved pitch angle resolution of and the increased geometric factor of as in Figures 7a-7c above. As expected from the above analysis, the PSD decrease is well resolved in both energies and pitch angles and the ratio reaches values much smaller than one at each energy level.

We conclude that with the particle detectors mounted on currently or recently operating spacecraft represented by Arase, Van Allen Probes, Cluster and MMS, a direct measurement of perturbations in electron velocity distribution due to interaction with a chorus subpackets is not possible, mainly due to insufficient angular resolution and low upper energy limits. With a dedicated experiment, however, an in situ observation of these perturbations should be within the current technical possibilities. Such an experiment would consist of particle analyzers measuring electrons with energies up to at least with a polar angle resolution better than about , and with a geometric factor of at least (i.e., 10-times better than MEP-e). Unlike with FPI on MMS, obtaining a fast 3D distribution would not be of essence: during high chorus wave activity, the analyzer would operate in a fixed high energy mode with accumulation times comparable to the theoretical time duration of a single subpacket. We believe that with this type of experiment, it should be possible to obtain direct confirmation of the processes described by the nonlinear growth theory of chorus rising tone emissions.